CellularPathFollower.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/ca/Path.h>
#include <tracking/trackFindingCDC/ca/AutomatonCell.h>
 
#include <tracking/trackFindingCDC/utilities/WeightedRelation.h>
 
#include <cassert>
#include <cmath>
#include <vector>
 
namespace Belle2 {
  namespace TrackFindingCDC {
 
    template<class ACellHolder>
    class CellularPathFollower {
 
    public:
      std::vector<Path<ACellHolder>> followAll(
                                    const std::vector<ACellHolder*>& cellHolders,
                                    const std::vector<WeightedRelation<ACellHolder>>& cellHolderRelations,
                                    Weight minStateToFollow = -INFINITY) const
      {
        B2ASSERT("Expected the relations to be sorted",
                 std::is_sorted(cellHolderRelations.begin(), cellHolderRelations.end()));
 
        // Result
        std::vector<Path<ACellHolder> > paths;
 
        // Stack for back tracking
        Path<ACellHolder> path;
 
        for (ACellHolder* cellHolder : cellHolders) {
          const AutomatonCell& automatonCell = cellHolder->getAutomatonCell();
 
          if (validStartCell(automatonCell, minStateToFollow)) {
            // Cell marks a start point of a path
 
            // Start new path
            path.clear();
 
            // Insert a pointer to the cell into the path
            path.push_back(cellHolder);
 
            // Recursively grow the path
            growAllPaths(path, cellHolderRelations, paths);
            path.pop_back();
          }
        }
        return paths;
      }
 
      Path<ACellHolder> followSingle(ACellHolder* startCellHolder,
                                     const std::vector<WeightedRelation<ACellHolder>>& cellHolderRelations,
                                     Weight minStateToFollow = -INFINITY) const
      {
        assert(std::is_sorted(cellHolderRelations.begin(), cellHolderRelations.end()));
 
        Path<ACellHolder> path;
        if (not startCellHolder) return path;
        const AutomatonCell& startCell = startCellHolder->getAutomatonCell();
        if (not validStartCell(startCell, minStateToFollow)) return path;
 
        // Start new path
        path.reserve(20); // Just a guess
 
        // Insert a pointer to the cell into the path
        path.push_back(startCellHolder);
        bool grew = true;
        while (grew) {
          grew = false;
          ACellHolder* cellHolder = path.back();
 
          auto continuations = asRange(std::equal_range(cellHolderRelations.begin(),
                                                        cellHolderRelations.end(),
                                                        cellHolder));
 
          for (const WeightedRelation<ACellHolder>& relation : continuations) {
            // cppcheck-suppress useStlAlgorithm
            if (isHighestContinuation(relation)) {
              ACellHolder* neighbor = relation.getTo();
              path.push_back(neighbor);
              grew = true;
              break;
            }
          }
        }
        return path;
      }
 
    private:
      void growAllPaths(Path<ACellHolder>& path,
                        const std::vector<WeightedRelation<ACellHolder>>& cellHolderRelations,
                        std::vector<Path<ACellHolder> >& paths) const
      {
        auto growPathByRelation = [&](const WeightedRelation<ACellHolder>& neighborRelation) {
          if (not this->isHighestContinuation(neighborRelation)) return false;
          ACellHolder* neighbor(neighborRelation.getTo());
          path.push_back(neighbor);
          this->growAllPaths(path, cellHolderRelations, paths);
          path.pop_back();
          return true;
        };
 
        ACellHolder* lastCellHolder = path.back();
 
        auto continuations = asRange(std::equal_range(cellHolderRelations.begin(),
                                                      cellHolderRelations.end(),
                                                      lastCellHolder));
        int nRelationsUsed = std::count_if(continuations.begin(),
                                           continuations.end(),
                                           growPathByRelation);
 
        if (nRelationsUsed == 0) {
          // end point of the recursion copy maximal path to the vector.
          paths.push_back(path);
        }
      }
 
      static bool validStartCell(const AutomatonCell& automatonCell,
                                 Weight minStateToFollow)
      {
        return
          automatonCell.hasStartFlag() and
          not automatonCell.hasMaskedFlag() and
          not automatonCell.hasCycleFlag() and
          minStateToFollow <= automatonCell.getCellState();
      }
 
      static bool isHighestContinuation(const WeightedRelation<ACellHolder>& relation)
      {
        ACellHolder* cellHolderPtr(relation.getFrom());
        ACellHolder* neighborCellHolderPtr(relation.getTo());
 
        if (not cellHolderPtr or not neighborCellHolderPtr) return false;
 
        ACellHolder& cellHolder = *cellHolderPtr;
        Weight relationWeight = relation.getWeight();
        ACellHolder& neighborCellHolder = *neighborCellHolderPtr;
 
        return isHighestContinuation(cellHolder, relationWeight, neighborCellHolder);
      }
 
      static bool isHighestContinuation(ACellHolder& cellHolder,
                                        Weight relationWeight,
                                        ACellHolder& neighborCellHolder)
      {
        const AutomatonCell& automatonCell = cellHolder.getAutomatonCell();
        const AutomatonCell& neighborAutomatonCell = neighborCellHolder.getAutomatonCell();
 
        return not neighborAutomatonCell.hasCycleFlag() and not neighborAutomatonCell.hasMaskedFlag() and
               (automatonCell.getCellState() ==
                (neighborAutomatonCell.getCellState() + relationWeight + automatonCell.getCellWeight()));
      }
    };
  }
}