QuadTreeProcessor.h

/**************************************************************************
* BASF2 (Belle Analysis Framework 2)                                     *
* Copyright(C) 2014 - Belle II Collaboration                             *
*                                                                        *
* Author: The Belle II Collaboration                                     *
* Contributors: Viktor Trusov, Thomas Hauth, Nils Braun                  *
*                                                                        *
* This software is provided "as is" without any warranty.                *
**************************************************************************/
#pragma once
 
#include <tracking/trackFindingCDC/legendre/quadtree/QuadTreeNode.h>
#include <tracking/trackFindingCDC/legendre/quadtree/QuadTreeItem.h>
 
#include <tracking/trackFindingCDC/utilities/Algorithms.h>
 
#include <algorithm>
#include <functional>
#include <memory>
#include <map>
#include <vector>
#include <deque>
#include <utility>
 
namespace Belle2 {
  namespace TrackFindingCDC {
 
    template<typename AX, typename AY, class AData>
    class QuadTreeProcessor {
 
    public:
      using Item = QuadTreeItem<AData>;
 
      using QuadTree = QuadTreeNode<AX, AY, Item>;
 
      using XSpan = typename QuadTree::XSpan;
 
      using YSpan = typename QuadTree::YSpan;
 
      using XYSpans = std::pair<XSpan, YSpan>;
 
      using QuadTreeChildren = typename QuadTree::Children;
 
      using CandidateReceiver = std::function<void(const std::vector<AData*>&, QuadTree*)>;
 
    public:
      QuadTreeProcessor(int lastLevel,
                        int seedLevel,
                        const XYSpans& xySpans,
                        bool debugOutput = false)
        : m_quadTree{std::make_unique<QuadTree>(xySpans.first, xySpans.second, 0, nullptr)}
        , m_lastLevel(lastLevel)
        , m_seedLevel(seedLevel)
        , m_debugOutput(debugOutput)
        , m_debugOutputMap()
      {
      }
 
      virtual ~QuadTreeProcessor()
      {
        clear();
      }
 
      void clear()
      {
        m_seededTrees.clear();
        m_quadTree->clearChildren();
        m_quadTree->clearItems();
        m_items.clear();
      }
 
      void seed(std::vector<AData*>& datas)
      {
        // Create the items
        for (AData* data : datas) {
          m_items.emplace_back(data);
        }
 
        // Creating the seed level
        long nSeedBins = pow(2, m_seedLevel);
        m_seededTrees.reserve(nSeedBins * nSeedBins);
 
        // Expand the first levels to the seed sectors
        m_seededTrees.push_back(m_quadTree.get());
        std::vector<QuadTree*> nextSeededTrees;
 
        for (int level = 0; level < m_seedLevel; ++level) {
          for (QuadTree* node : m_seededTrees) {
            if (node->getChildren().empty()) {
              this->createChildren(node, node->getChildren());
            }
            for (QuadTree& child :  node->getChildren()) {
              nextSeededTrees.push_back(&child);
            }
          }
          std::swap(nextSeededTrees, m_seededTrees);
          nextSeededTrees.clear();
        }
 
        // Fill the seed level with the items
        for (QuadTree* seededTree : m_seededTrees) {
          seededTree->reserveItems(m_items.size());
 
          for (Item& item : m_items) {
            if (item.isUsed()) continue;
            if (isInNode(seededTree, item.getPointer())) {
              seededTree->insertItem(&item);
            }
          }
        }
      }
 
    public:
      std::vector<AData*> getAssignedItems()
      {
        std::vector<const CDCWireHit*> result;
        for (QuadTree* seededTree : m_seededTrees) {
          for (Item* item : seededTree->getItems()) {
            result.push_back(item->getPointer());
          }
        }
        std::sort(result.begin(), result.end());
        result.erase(std::unique(result.begin(), result.end()), result.end());
        return result;
      }
 
    public:
      void fill(const CandidateReceiver& candidateReceiver, int nHitsThreshold)
      {
        fill(candidateReceiver, nHitsThreshold, std::numeric_limits<AY>::max());
      }
 
      void fill(const CandidateReceiver& candidateReceiver, int nHitsThreshold, float yLimit)
      {
        std::vector<QuadTree*> quadTrees = m_seededTrees;
        std::sort(quadTrees.begin(), quadTrees.end(), [](QuadTree * quadTree1, QuadTree * quadTree2) {
          return quadTree1->getNItems() > quadTree2->getNItems();
        });
 
        for (QuadTree* tree : quadTrees) {
          erase_remove_if(tree->getItems(), [](Item * hit) { return hit->isUsed(); });
          fillGivenTree(tree, candidateReceiver, nHitsThreshold, yLimit);
        }
      }
 
    private:
      void fillGivenTree(QuadTree* node,
                         const CandidateReceiver& candidateReceiver,
                         int nItemsThreshold,
                         AY yLimit)
      {
        if (node->getNItems() < nItemsThreshold) {
          return;
        }
 
        if ((node->getYMin() > yLimit) or (-node->getYMax() > yLimit)) {
          return;
        }
 
        if (isLeaf(node)) {
          callResultFunction(node, candidateReceiver);
          return;
        }
 
        if (node->getChildren().empty()) {
          this->createChildren(node, node->getChildren());
        }
 
        if (!node->checkFilled()) {
          fillChildren(node, node->getItems());
          node->setFilled();
        }
 
        std::vector<QuadTree*> children;
        for (QuadTree& child : node->getChildren()) {
          children.push_back(&child);
        }
        const auto compareNItems = [](const QuadTree * lhs, const QuadTree * rhs) {
          return lhs->getNItems() < rhs->getNItems();
        };
 
        // Explicitly count down the children
        const int nChildren = children.size();
        for (int iChild = 0; iChild < nChildren; ++iChild) {
          auto itHeaviestChild = std::max_element(children.begin(), children.end(), compareNItems);
          QuadTree* heaviestChild = *itHeaviestChild;
          children.erase(itHeaviestChild);
          // After we have processed some children we need to get rid of the already used hits in all the children,
          // because this can change the number of items drastically
          erase_remove_if(heaviestChild->getItems(), [&](Item * hit) { return hit->isUsed(); });
          this->fillGivenTree(heaviestChild, candidateReceiver, nItemsThreshold, yLimit);
        }
      }
 
      void createChildren(QuadTree* node, QuadTreeChildren& m_children) const
      {
        for (int i = 0; i < node->getXNbins(); ++i) {
          for (int j = 0; j < node->getYNbins(); ++j) {
            const XYSpans& xySpans = createChild(node, i, j);
            const XSpan& xSpan = xySpans.first;
            const YSpan& ySpan = xySpans.second;
            m_children.push_back(QuadTree(xSpan, ySpan, node->getLevel() + 1, node));
          }
        }
      }
 
      void fillChildren(QuadTree* node, std::vector<Item*>& items)
      {
        const size_t neededSize = 2 * items.size();
        for (QuadTree& child : node->getChildren()) {
          child.reserveItems(neededSize);
        }
 
        for (Item* item : items) {
          if (item->isUsed()) continue;
 
          for (QuadTree& child : node->getChildren()) {
            if (isInNode(&child, item->getPointer())) {
              child.insertItem(item);
            }
          }
        }
        afterFillDebugHook(node->getChildren());
      }
 
      void callResultFunction(QuadTree* node, const CandidateReceiver& candidateReceiver) const
      {
        const std::vector<Item*>& foundItems = node->getItems();
        std::vector<AData*> candidate;
        candidate.reserve(foundItems.size());
 
        for (Item* item : foundItems) {
          item->setUsedFlag();
          candidate.push_back(item->getPointer());
        }
 
        candidateReceiver(candidate, node);
      }
 
    protected: // Section of specialisable functions
      virtual XYSpans createChild(QuadTree* node, int iX, int iY) const
      {
        AX xMin = node->getXLowerBound(iX);
        AX xMax = node->getXUpperBound(iX);
        AY yMin = node->getYLowerBound(iY);
        AY yMax = node->getYUpperBound(iY);
        return XYSpans({xMin, xMax}, {yMin, yMax});
      }
 
      virtual bool isInNode(QuadTree* node, AData* item) const = 0;
 
      virtual bool isLeaf(QuadTree* node) const
      {
        if (node->getLevel() >= m_lastLevel) {
          return true;
        } else {
          return false;
        }
      }
 
      int getLastLevel() const
      {
        return m_lastLevel;
      }
 
    public: // debug stuff
      virtual void afterFillDebugHook(QuadTreeChildren& children)
      {
        if (not m_debugOutput) return;
        for (QuadTree& childNode : children) {
          if (childNode.getLevel() != getLastLevel()) continue; // Only write the lowest level
          //m_debugOutputMap[ {childNode.getXMean(), childNode.getYMean()}] = childNode.getItems();
        }
      }
 
      const std::map<std::pair<AX, AY>, std::vector<Item*>>& getDebugInformation() const
      {
        return m_debugOutputMap;
      }
 
    protected:
      std::unique_ptr<QuadTree> m_quadTree;
 
      std::deque<Item> m_items;
 
      std::vector<QuadTree*> m_seededTrees;
 
    private:
      int m_lastLevel;
 
      int m_seedLevel;
 
      bool m_debugOutput;
 
      std::map<std::pair<AX, AY>, std::vector<Item*>> m_debugOutputMap;
    };
  }
}