BoxDivisionHoughTree.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
#pragma once
#include <tracking/trackFindingCDC/hough/trees/WeightedFastHoughTree.h>
#include <tracking/trackFindingCDC/hough/baseelements/SectoredLinearDivision.h>
 
#include <tracking/trackFindingCDC/numerics/LookupTable.h>
 
#include <tracking/trackFindingCDC/utilities/TupleGenerate.h>
 
#include <type_traits>
#include <utility>
#include <tuple>
#include <array>
#include <memory>
 
namespace Belle2 {
  namespace TrackFindingCDC {

    template <class AItemPtr, class AHoughBox, size_t ... divisions>
    class BoxDivisionHoughTree {
 
    public:
      using HoughBox = AHoughBox;

      using BoxDivision = SectoredLinearDivision<HoughBox, divisions...>;

      using HoughTree = WeightedFastHoughTree<AItemPtr, HoughBox, BoxDivision>;

      template <size_t I>
      using Type = typename HoughBox::template Type<I>;

      template <class T>
      using HasType = typename HoughBox::template HasType<T>;

      template <class T>
      using TypeIndex = typename HoughBox::template TypeIndex<T>;

      template <size_t I>
      using Width = typename HoughBox::template Width<I>;

      using Node = typename  HoughTree::Node;
 
    public:
      explicit BoxDivisionHoughTree(int maxLevel, int sectorLevelSkip = 0)
        : m_maxLevel(maxLevel)
        , m_sectorLevelSkip(sectorLevelSkip)
        , m_overlaps((divisions * 0)...)
      {
      }
 
    private:
      template <size_t I>
      using Array = typename Type<I>::Array;

      using Arrays = TupleGenerateN<Array, sizeof...(divisions)>;
 
    public:
      size_t getDivision(size_t i) const
      {
        return m_divisions[i];
      }

      template <size_t I>
      void constructArray(double lowerBound,
                          double upperBound,
                          Width<I> nBinOverlap = 0,
                          Width<I> nBinWidth = 0)
      {
        static_assert(std::is_integral<Width<I> >::value, "Method only applicable to discrete axes");
        const size_t division = getDivision(I);
        const int granularityLevel = m_maxLevel + m_sectorLevelSkip;
        const size_t nBins = std::pow(division, granularityLevel);
 
        if (nBinWidth == 0) {
          nBinWidth = nBinOverlap + 1;
        }
 
        B2ASSERT("Width " << nBinWidth << "is not bigger than overlap " << nBinOverlap,
                 nBinOverlap < nBinWidth);
 
        const auto nPositions = (nBinWidth - nBinOverlap) * nBins + nBinOverlap + 1;
        std::get<I>(m_arrays) = linspace<float>(lowerBound, upperBound, nPositions);
        std::get<I>(m_overlaps) = nBinOverlap;
      }

      template <size_t I>
      void assignArray(Array<I> array, Width<I> overlap = 0)
      {
        std::get<I>(m_arrays) = std::move(array);
        std::get<I>(m_overlaps) = overlap;
 
        // In case of a discrete axes, check whether the size of the array is sensible
        // such that the bin width at the highest granularity level is a whole number given the overlap.
        if (std::is_integral<Width<I> >::value) {
          const int division = getDivision(I);
          const int granularityLevel = m_maxLevel + m_sectorLevelSkip;
          const long int nBins = std::pow(division, granularityLevel);
          const long int nPositions = std::get<I>(m_arrays).size();
          const long int nWidthTimeNBins = nPositions - 1 + (nBins - 1) * overlap;
 
          B2ASSERT("Not enough positions in array to cover all bins.\n"
                   "Expected: positions = " << nBins - (nBins - 1) * overlap + 1 << " at least.\n"
                   "Actual:   positions = " << nPositions << " (overlap = " << overlap << ", bins = " << nBins << ")\n",
                   nWidthTimeNBins >= nBins);
 
          B2ASSERT("Number of positions in array introduces inhomogeneous width at the highest granularity level.\n"
                   "Expected: positions = 'width * bins - (bins - 1) * overlap + 1'\n"
                   "Actual:   positions = " << nPositions << " (overlap = " << overlap << ", bins = " << nBins << ")\n",
                   nWidthTimeNBins % nBins == 0);
        }
      }

      template <class T>
      std::enable_if_t< HasType<T>::value, void>
      assignArray(Array<TypeIndex<T>::value > array, Width<TypeIndex<T>::value > overlap = 0)
      {
        assignArray<TypeIndex<T>::value>(std::move(array), overlap);
      }
 
    public:
      void initialize()
      {
        // Compose the hough space
        HoughBox houghPlane = constructHoughPlane();
        BoxDivision boxDivision(m_overlaps, m_sectorLevelSkip);
        m_houghTree.reset(new HoughTree(std::move(houghPlane), std::move(boxDivision)));
      }

      template <class AItemPtrs>
      void seed(const AItemPtrs& items)
      {
        if (not m_houghTree) initialize();
        m_houghTree->seed(items);
      }

      void fell()
      {
        m_houghTree->fell();
      }

      void raze()
      {
        m_houghTree->raze();
      }
 
    public:
      HoughTree* getTree() const
      {
        return m_houghTree.get();
      }

      int getMaxLevel() const
      {
        return m_maxLevel;
      }

      void setMaxLevel(int maxLevel)
      {
        m_maxLevel = maxLevel;
      }

      int getSectorLevelSkip() const
      {
        return m_sectorLevelSkip;
      }

      void setSectorLevelSkip(int sectorLevelSkip)
      {
        m_sectorLevelSkip = sectorLevelSkip;
      }

      template <size_t I>
      const Array<I>& getArray() const
      {
        return std::get<I>(m_arrays);
      }
 
    private:
      template <size_t... Is>
      HoughBox constructHoughPlaneImpl(const std::index_sequence<Is...>& is __attribute__((unused)))
      {
        return HoughBox(Type<Is>::getRange(std::get<Is>(m_arrays))...);
      }

      HoughBox constructHoughPlane()
      {
        return constructHoughPlaneImpl(std::make_index_sequence<sizeof...(divisions)>());
      }
 
    private:
      const std::array<size_t, sizeof ...(divisions)> m_divisions = {{divisions ...}};

      std::unique_ptr<HoughTree> m_houghTree = nullptr;

      int m_maxLevel;

      int m_sectorLevelSkip;

      typename HoughBox::Delta m_overlaps;

      Arrays m_arrays;
    };
  }
}