Partitions.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_PARTITIONS_H
#define LARCOREALG_GEOMETRY_PARTITIONS_H

// LArSoft libraries
#include "larcorealg/CoreUtils/DebugUtils.h" // lar::debug::demangle(), ...
#include "larcorealg/Geometry/SimpleGeo.h"   // lar::util::simple_geo::Rectangle

// C/C++ standard libraries
#include <algorithm> // std::upper_bound()
#include <cassert>
#include <cstdlib>  // std::size_t
#include <iterator> // std::next(), ...
#include <memory>   // std::unique_ptr<>
#include <sstream>
#include <string>
#include <type_traits> // std::declval(), std::true_type...
#include <utility>     // std::move(), std::forward()
#include <vector>

namespace geo::part {

  //-------------------------------------------------------------------------
  class AreaOwner {
  public:
    using Area_t = lar::util::simple_geo::Rectangle<double>;

    using AreaRangeMember_t = Area_t::Range_t(Area_t::*);

    AreaOwner(Area_t const& area) : myArea(area) {}

    bool contains(double w, double d) const { return area().contains(w, d); }

    Area_t const& area() const { return myArea; }

    template <typename Stream>
    void dumpArea(Stream&& out) const
    {
      std::forward<Stream>(out) << area();
    }

  private:
    Area_t myArea; 

  }; // class AreaOwner

  namespace details {

    template <AreaOwner::AreaRangeMember_t Range>
    struct PartitionSorterByAreaRangeLower;

  } // namespace details

  //*************************************************************************
  //***  Fundamental classes
  //-------------------------------------------------------------------------

  template <typename Data>
  struct PartitionDataDescriber {
    template <typename Stream>
    PartitionDataDescriber(Stream&& out,
                           Data const* data,
                           std::string indent = "",
                           std::string firstIndent = "");
  }; // struct PartitionDataDescriber<>

  template <typename Stream, typename Data>
  void describePartitionData(Stream&& out,
                             Data const* data,
                             std::string indent = "",
                             std::string firstIndent = "");

  //--------------------------------------------------------------------------
  class PartitionBase : public AreaOwner {
  public:
    // imported types
    using Area_t = AreaOwner::Area_t;
    using AreaRangeMember_t = AreaOwner::AreaRangeMember_t;

    PartitionBase(Area_t const& area) : AreaOwner(area) {}

  protected:
    std::string describeArea(std::string indent, std::string firstIndent) const;

  }; // class PartitionBase

  //--------------------------------------------------------------------------
  template <typename Data>
  class Partition : public PartitionBase {
  public:
    using Data_t = Data;                  
    using Partition_t = Partition<Data>;  
    using Area_t = PartitionBase::Area_t; 

    using Subpartitions_t = std::vector<std::unique_ptr<Partition_t const>>;

    Partition(Area_t const& area) : PartitionBase(area) {}

    virtual ~Partition() = default;

    virtual Data_t* data() const { return nullptr; }

    virtual Data_t* atPoint(double w, double d) const = 0;

    std::string describe(std::string indent, std::string firstIndent) const
    {
      return doDescribe(indent, firstIndent);
    }

    std::string describe(std::string indent = "") const { return describe(indent, indent); }

    template <typename Pred>
    void walk(Pred&& pred) const
    {
      walk(this, pred);
    }

    std::size_t nParts() const { return parts().size(); }

  protected:
    static Subpartitions_t const NoSubparts; 

    virtual Subpartitions_t const& parts() const { return NoSubparts; }

    virtual std::string doDescribe(std::string indent, std::string firstIndent) const
    {
      return PartitionBase::describeArea(indent, firstIndent);
    }

    template <typename Pred>
    static void walk(Partition_t const* start, Pred&& pred);

  }; // class Partition<>

  //**************************************************************************
  //***  Partition class hierarchy
  //--------------------------------------------------------------------------
  template <typename Data>
  class PartitionWithData : public Partition<Data> {

  public:
    using Base_t = Partition<Data>;              
    using Partition_t = Partition<Data>;         
    using Data_t = typename Partition_t::Data_t; 
    using Area_t = typename Partition_t::Area_t; 

    PartitionWithData(Area_t const& area, Data_t* myData) : Base_t(area), myData(myData) {}

    virtual Data_t* data() const override { return myData; }

    virtual Data_t* atPoint(double w, double d) const override
    {
      return Base_t::contains(w, d) ? myData : nullptr;
    }

  private:
    Data_t* myData; 

    virtual std::string doDescribe(std::string indent, std::string firstIndent) const override;

  }; // class PartitionWithData

  //--------------------------------------------------------------------------
  template <typename Data>
  class PartitionElement : public PartitionWithData<Data> {
  public:
    using Base_t = PartitionWithData<Data>; 

    // Import constructors
    using Base_t::Base_t;
  };

  //--------------------------------------------------------------------------
  template <typename Data>
  class PartitionContainer : public PartitionWithData<Data> {
  public:
    using Base_t = PartitionWithData<Data>;
    using Partition_t = Partition<Data>; 

    // inherited types
    using Data_t = typename Partition_t::Data_t;
    using Area_t = typename Partition_t::Area_t;
    using Subpartitions_t = typename Partition_t::Subpartitions_t;

    virtual Data_t* atPoint(double w, double d) const override;

  protected:
    Subpartitions_t myParts; 

    std::size_t size() const { return parts().size(); }

    virtual Subpartitions_t const& parts() const override { return myParts; }

    PartitionContainer(Area_t const& area,
                       Subpartitions_t&& subpartitions,
                       Data_t* defData = nullptr)
      : Base_t(area, defData), myParts(std::move(subpartitions))
    {}

    virtual Partition_t const* findPart(double w, double d) const = 0;

    virtual std::string doDescribe(std::string indent, std::string firstIndent) const override;

    virtual std::string describeIntro() const;

  }; // class PartitionContainer<>

  //--------------------------------------------------------------------------
  template <typename Data, typename Sorter>
  class SortedPartition : public PartitionContainer<Data> {
  public:
    using Base_t = PartitionContainer<Data>;
    using Partition_t = Partition<Data>; 
    using Sorter_t = Sorter;             

    // inherited data types
    using Data_t = typename Partition_t::Data_t;
    using Area_t = typename Partition_t::Area_t;
    using Subpartitions_t = typename Base_t::Subpartitions_t;

    SortedPartition(Area_t const& area,
                    Subpartitions_t&& subpartitions,
                    Data_t* defData = nullptr,
                    Sorter_t sorter = {})
      : Base_t(area, std::move(subpartitions), defData), sorter(sorter)
    {
      initParts();
    }

  protected:
    Sorter_t sorter; 

    Partition_t const* findPartWithKey(double key) const;

    void initParts();

  }; // class SortedPartition<>

  //--------------------------------------------------------------------------
  template <typename Data, PartitionBase::AreaRangeMember_t Range>
  class PartitionSortedByRange
    : public SortedPartition<Data, details::PartitionSorterByAreaRangeLower<Range>> {
  public:
    using Base_t = SortedPartition<Data, details::PartitionSorterByAreaRangeLower<Range>>;
    using Base_t::Base_t; // import inherited constructors
  };

  //--------------------------------------------------------------------------
  template <typename Data>
  class DepthPartition : public PartitionSortedByRange<Data, &PartitionBase::Area_t::depth> {
  public:
    using Base_t = PartitionSortedByRange<Data, &PartitionBase::Area_t::depth>;

    using Base_t::Base_t; // import inherited constructors

    virtual typename Base_t::Partition_t const* findPart(double /* w */, double d) const override
    {
      return Base_t::findPartWithKey(d);
    }

  private:
    virtual std::string describeIntro() const override;

  }; // class DepthPartition

  //--------------------------------------------------------------------------
  template <typename Data>
  class WidthPartition : public PartitionSortedByRange<Data, &PartitionBase::Area_t::width> {
  public:
    using Base_t = PartitionSortedByRange<Data, &PartitionBase::Area_t::width>;

    using Base_t::Base_t; // inherited constructors

    virtual typename Base_t::Partition_t const* findPart(double w, double /* d */) const override
    {
      return Base_t::findPartWithKey(w);
    }

  private:
    virtual std::string describeIntro() const override;

  }; // class WidthPartition

  //--------------------------------------------------------------------------
  template <typename Data>
  class GridPartition : public PartitionContainer<Data> {
  public:
    using Partition_t = Partition<Data>;     
    using Base_t = PartitionContainer<Data>; 

    // Inherited types
    using Data_t = typename Partition_t::Data_t; 
    using Area_t = typename Partition_t::Area_t; 
    using Subpartitions_t = typename Base_t::Subpartitions_t;

    GridPartition(Area_t const& area,
                  Subpartitions_t&& subpartitions,
                  unsigned int nDepthPartitions,
                  unsigned int nWidthPartitions,
                  Data_t* defData = nullptr);

    GridPartition(Area_t const& area,
                  Subpartitions_t&& subpartitions,
                  unsigned int nDepthPartitions,
                  Data_t* defData = nullptr);

  private:
    std::vector<double> widthSeps; 
    std::vector<double> depthSeps; 

    std::size_t nWidthParts() const { return widthSeps.size(); }
    std::size_t nDepthParts() const { return depthSeps.size(); }

    auto part(std::size_t iDepth, std::size_t iWidth) -> decltype(auto)
    {
      return Base_t::parts()[iDepth * nWidthParts() + iWidth];
    }
    auto part(std::size_t iDepth, std::size_t iWidth) const -> decltype(auto)
    {
      return Base_t::parts()[iDepth * nWidthParts() + iWidth];
    }

    virtual Partition_t const* findPart(double w, double d) const override;

    std::vector<double> computeWidthSeps(unsigned int nD, unsigned int nW) const;
    std::vector<double> computeDepthSeps(unsigned int nD, unsigned int nW) const;

    virtual std::string doDescribe(std::string indent, std::string firstIndent) const override;

    template <PartitionBase::AreaRangeMember_t Range, typename BeginIter, typename EndIter>
    static std::vector<double> detectSeparators(BeginIter b,
                                                EndIter e,
                                                std::size_t const nGroups,
                                                std::size_t const startDelta,
                                                std::size_t const stride);

  }; // class GridPartition<>

  //--------------------------------------------------------------------------

} // namespace geo::part

//******************************************************************************
//***  inline implementation
//******************************************************************************
inline std::string geo::part::PartitionBase::describeArea(std::string,
                                                          std::string firstIndent) const
{
  std::ostringstream sstr;
  sstr << firstIndent << "partition covers ";
  dumpArea(sstr);
  return sstr.str();
}

//******************************************************************************
//***  template implementation
//******************************************************************************

namespace geo::part::details {

  //------------------------------------------------------------------------
  //---  some metaprogramming utilities
  //------------------------------------------------------------------------
  template <typename T, typename = void>
  struct is_partition_type : public std::false_type {};

  template <typename Part>
  struct is_partition_type<
    Part,
    std::enable_if_t<std::is_base_of<PartitionBase, std::decay_t<Part>>::value>>
    : public std::true_type {};

  template <typename T>
  constexpr bool is_partition_type_v = is_partition_type<T>();

  //------------------------------------------------------------------------
  template <typename, typename = void>
  struct is_partition_ptr : public std::false_type {};

  template <typename PartPtr>
  struct is_partition_ptr<PartPtr,
                          std::enable_if_t<is_partition_type_v<decltype(*std::declval<PartPtr>())>>>
    : public std::true_type {};

  template <typename T>
  constexpr bool is_partition_ptr_v = is_partition_ptr<T>();

  //------------------------------------------------------------------------
  template <typename, typename = void>
  struct is_partition_ptr_iterator : public std::false_type {};

  template <typename Iter>
  struct is_partition_ptr_iterator<
    Iter,
    std::enable_if_t<is_partition_ptr_v<std::decay_t<typename Iter::value_type>>>>
    : public std::true_type {};

  //------------------------------------------------------------------------
  //---  Sorting objects
  //------------------------------------------------------------------------
  template <AreaOwner::AreaRangeMember_t Range>
  struct RangeLowerBoundExtractor {
    static constexpr auto range = Range;
    using Area_t = AreaOwner::Area_t;

    double operator()(double lower) const { return lower; }
    double operator()(Area_t::Range_t const& r) const { return (*this)(r.lower); }
    double operator()(Area_t const& area) const { return (*this)(area.*range); }
    double operator()(AreaOwner const& area) const { return (*this)(area.area()); }
    double operator()(AreaOwner const* ptr) const { return (*this)(*ptr); }
  };

  //------------------------------------------------------------------------
  template <PartitionBase::AreaRangeMember_t Range>
  struct PartitionRangeLowerBoundExtractor : public RangeLowerBoundExtractor<Range> {
    using Base_t = RangeLowerBoundExtractor<Range>;
    using Area_t = typename Base_t::Area_t;

    using Partition_t = PartitionBase; 

    using Base_t::operator(); // import inherited versions

    auto operator()(Partition_t const& part) { return Base_t::operator()(part.area()); }
    template <typename PartPtr,
              typename = std::enable_if_t<details::is_partition_ptr<PartPtr>::value>>
    auto operator()(PartPtr const& part)
    {
      return operator()(*part);
    }
  };

  //------------------------------------------------------------------------
  template <PartitionBase::AreaRangeMember_t Range>
  struct PartitionSorterByAreaRangeLower {
    using Sorter_t = PartitionSorterByAreaRangeLower<Range>;

    using KeyExtractor_t = PartitionRangeLowerBoundExtractor<Range>;

    static constexpr auto range = Range;

    template <typename T>
    static auto key(T const& obj)
    {
      return KeyExtractor_t()(obj);
    }

    using Key_t = decltype(Sorter_t::key(std::declval<PartitionBase>()));

    static Key_t key(Key_t k) { return k; } // shortcut
    static bool sortKey(Key_t a, Key_t b) { return a < b; }

    template <typename A, typename B>
    bool operator()(A const& a, B const& b) const
    {
      return sortKey(key(a), key(b));
    }
  };

  //------------------------------------------------------------------------

} // namespace geo::part::details

//------------------------------------------------------------------------------
//---  partition data description
//---

template <typename Data>
template <typename Stream>
geo::part::PartitionDataDescriber<Data>::PartitionDataDescriber(Stream&& out,
                                                                Data const* data,
                                                                std::string indent /* = "" */,
                                                                std::string firstIndent /* = "" */
)
{
  out << firstIndent;
  std::string typeName = lar::debug::demangle<Data>();
  if (data) { out << typeName << "[" << ((void*)data) << "]"; }
  else {
    out << "no '" << typeName << "' data";
  }
}

template <typename Stream, typename Data>
void geo::part::describePartitionData(Stream&& out,
                                      Data const* data,
                                      std::string indent /* = "" */,
                                      std::string firstIndent /* = "" */
)
{
  PartitionDataDescriber<Data>(std::forward<Stream>(out), data, indent, firstIndent);
}

//------------------------------------------------------------------------------
//--- geo::part::Partition
//---
template <typename Data>
typename geo::part::Partition<Data>::Subpartitions_t const geo::part::Partition<Data>::NoSubparts;

//------------------------------------------------------------------------------
template <typename Data>
template <typename Pred>
void geo::part::Partition<Data>::walk(Partition_t const* start, Pred&& pred)
{
  if (!start) return;
  pred(*start);

  // recursive implementation
  for (auto const& subPart : start->parts())
    subPart->walk(std::forward<Pred>(pred));
}

//------------------------------------------------------------------------------
//---  geo::part::PartitionWithData
//---
template <typename Data>
std::string geo::part::PartitionWithData<Data>::doDescribe(std::string indent,
                                                           std::string firstIndent) const
{
  std::string msg = Base_t::doDescribe(indent, firstIndent);
  if (data()) {
    std::ostringstream sstr;
    sstr << ": ";
    describePartitionData(sstr, data(), indent);
    msg += sstr.str();
  }
  else {
    msg += " (no data)";
  }
  return msg;
}

//------------------------------------------------------------------------------
//--- geo::part::PartitionContainer
//---
//------------------------------------------------------------------------------
template <typename Data>
auto geo::part::PartitionContainer<Data>::atPoint(double w, double d) const -> Data_t*
{
  if (!Base_t::contains(w, d)) return nullptr; // not our point at all
  // it's ours; see if it belongs to a subpart
  auto part = findPart(w, d);
  return part ? part->atPoint(w, d) : Base_t::data();
}

//------------------------------------------------------------------------------
template <typename Data>
std::string geo::part::PartitionContainer<Data>::doDescribe(std::string indent,
                                                            std::string firstIndent) const
{
  std::string msg = firstIndent + describeIntro();
  if (Base_t::data()) {
    std::ostringstream sstr;
    sstr << ", and ";
    describePartitionData(sstr, Base_t::data(), indent, "");
    msg += sstr.str();
  }

  for (auto const& part : parts()) {
    msg += "\n" + indent + " * ";
    msg += part->describe(indent + "  ", "");
  }

  return msg;
}

//------------------------------------------------------------------------------
template <typename Data>
std::string geo::part::PartitionContainer<Data>::describeIntro() const
{
  return std::to_string(parts().size()) + " subpartitions";
}

//------------------------------------------------------------------------------
//---  geo::part::SortedPartition
//---
template <typename Data, typename Sorter>
auto geo::part::SortedPartition<Data, Sorter>::findPartWithKey(double key) const
  -> Partition_t const*
{
  auto pbegin = Base_t::parts().cbegin();
  auto iPart = std::upper_bound(pbegin, Base_t::parts().cend(), key, sorter);
  return (iPart == pbegin) ? nullptr : (--iPart)->get();
}

//------------------------------------------------------------------------------
template <typename Data, typename Sorter>
void geo::part::SortedPartition<Data, Sorter>::initParts()
{
  /*
   * Initialization tasks:
   * - ensure that the parts are sorted by increasing depth
    */
  std::sort(Base_t::myParts.begin(), Base_t::myParts.end(), sorter);
}

//------------------------------------------------------------------------------
//---  geo::part::DepthPartition
//---
template <typename Data>
std::string geo::part::DepthPartition<Data>::describeIntro() const
{
  std::ostringstream sstr;
  sstr << Base_t::size() << " partitions along depth covering " << Base_t::area();
  return sstr.str();
}

//------------------------------------------------------------------------------
//---  geo::part::WidthPartition
//---
template <typename Data>
std::string geo::part::WidthPartition<Data>::describeIntro() const
{
  std::ostringstream sstr;
  sstr << Base_t::size() << " partitions along width covering " << Base_t::area();
  return sstr.str();
}

//------------------------------------------------------------------------------
//---  geo::part::GridPartition
//---
template <typename Data>
geo::part::GridPartition<Data>::GridPartition(Area_t const& area,
                                              Subpartitions_t&& subpartitions,
                                              unsigned int nDepthPartitions,
                                              unsigned int nWidthPartitions,
                                              Data_t* defData /* = nullptr */
                                              )
  : Base_t(area, std::move(subpartitions), defData)
  , widthSeps(computeWidthSeps(nDepthPartitions, nWidthPartitions))
  , depthSeps(computeDepthSeps(nDepthPartitions, nWidthPartitions))
{
  assert(nWidthPartitions * nDepthPartitions == Base_t::size());
}

template <typename Data>
geo::part::GridPartition<Data>::GridPartition(Area_t const& area,
                                              Subpartitions_t&& subpartitions,
                                              unsigned int nDepthPartitions,
                                              Data_t* defData /* = nullptr */
                                              )
  : GridPartition(area,
                  std::move(subpartitions),
                  nDepthPartitions,
                  (nDepthPartitions ? subpartitions.size() / nDepthPartitions : 0),
                  defData)
{}

//------------------------------------------------------------------------------
template <typename Data>
std::string geo::part::GridPartition<Data>::doDescribe(std::string indent,
                                                       std::string firstIndent) const
{
  std::ostringstream sstr;
  sstr << firstIndent << Base_t::describeIntro() << " in a (WxD) = " << nWidthParts() << " x "
       << nDepthParts() << " grid";
  if (Base_t::data()) {
    sstr << ", and ";
    describePartitionData(sstr, Base_t::data(), indent, "");
  }
  for (std::size_t iDepth = 0; iDepth < nDepthParts(); ++iDepth) {
    for (std::size_t iWidth = 0; iWidth < nWidthParts(); ++iWidth) {
      sstr << "\n"
           << indent << " [" << iDepth << "][" << iWidth << "] "
           << part(iDepth, iWidth)->describe(indent + "  ", "");
    } // for width
  }   // for depth
  return sstr.str();
}

//------------------------------------------------------------------------------
template <typename Data>
auto geo::part::GridPartition<Data>::findPart(double w, double d) const -> Partition_t const*
{
  auto const iWidth = std::upper_bound(widthSeps.cbegin(), widthSeps.cend(), w);
  if (iWidth == widthSeps.cbegin()) return nullptr;
  auto const iDepth = std::upper_bound(depthSeps.cbegin(), depthSeps.cend(), d);
  if (iDepth == depthSeps.cbegin()) return nullptr;
  return part(std::distance(depthSeps.cbegin(), iDepth) - 1U,
              std::distance(widthSeps.cbegin(), iWidth) - 1U)
    .get();
}

//------------------------------------------------------------------------------
template <typename Data>
std::vector<double> geo::part::GridPartition<Data>::computeWidthSeps(unsigned int,
                                                                     unsigned int nW) const
{
  return detectSeparators<&Area_t::width>(
    Base_t::parts().cbegin(), Base_t::parts().cend(), nW, 1U, nW);
}

template <typename Data>
std::vector<double> geo::part::GridPartition<Data>::computeDepthSeps(unsigned int nD,
                                                                     unsigned int nW) const
{
  return detectSeparators<&Area_t::depth>(
    Base_t::parts().cbegin(), Base_t::parts().cend(), nD, nW, 1U);
}

template <typename Data>
template <geo::part::AreaOwner::AreaRangeMember_t Range, typename BeginIter, typename EndIter>
std::vector<double> geo::part::GridPartition<Data>::detectSeparators(BeginIter b,
                                                                     EndIter e,
                                                                     std::size_t const nGroups,
                                                                     std::size_t const startDelta,
                                                                     std::size_t const stride)
{
  /*
   * The iterators are better be random access.
   * The range [b,e[ is considered to be a 2D table stored row after row.
   * This function can operate on rows or columns, given the proper arguments.
   *
   * The full range is split in nGroups "groups" (e.g. rows or columns).
   * Each group g starts at the element g x startDelta.
   * All members of that group are stride elements far from each other.
   *
   * Be the table of size (d x w).
   * To process data row by row:
   * - the group is a row
   * - the number of groups is the number of rows: nGroups = d
   * - each group will have (d x w)/nGroups = w elements
   * - the start offset is the number of group times the size of it:
   *   startDelta = w
   * - the elements are contiguous: stride = 1
   *
   * To process data column by column:
   * - the group is a column
   * - the number of groups is the number of columns: nGroups = w
   * - each group will have (d x w)/nGroups = d elements
   * - the start offset matches the number of group: startDelta = 1
   * - the elements are separated by a full row: stride = w
   *
   */
  // the separator is on the lower bound of selected range of the partition area

  static_assert(details::is_partition_ptr_iterator<BeginIter>(),
                "Begin iterator does not point to a pointer to partition type");

  details::PartitionRangeLowerBoundExtractor<Range> lowerBound;

  std::size_t const nParts = std::distance(b, e);
  std::size_t const nPartsInGroup = nParts / nGroups;

  auto const part = [b](std::size_t index) { return std::next(b, index)->get(); };

  std::vector<double> seps(nGroups);
  for (size_t g = 0; g < nGroups; ++g) {

    double& sep = seps[g];

    // indices of an element in the previous and next group, respectively
    std::size_t index = g * startDelta;
    sep = lowerBound(part(index));

    std::size_t const iend = index + nPartsInGroup * stride;
    while ((index += stride) < iend) {
      double const l = lowerBound(part(index));
      if (sep > l) sep = l;
    } // while (element)

  } // while (groups)
  return seps;
}

//******************************************************************************

#endif // LARCOREALG_GEOMETRY_PARTITIONS_H