Added faster, in-place clusterization.

The in-place clusterization could be used instead of the default
clusterization. It clusterizes by sorting a cell container in
such a way that adjacent cells with the same label belong to the
same cluster.

Author: Goetz Gaycken

Core/include/Acts/Clusterization/InPlaceClusterization.hpp

@@ -82,6 +82,9 @@ class DigitizationAlgorithm final : public IAlgorithm {
     /// Minimum number of attempts to derive a valid dgitized measurement when
     /// random numbers are involved.
     std::size_t minMaxRetries = 10;
+    /// Use InPlaceClusterization clusterization instead of the original
+    /// clusterization if doMerge is true
+    bool useInPlaceClusterization = false;
   };
 
   /// Construct the smearing algorithm.

Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp

@@ -40,12 +40,13 @@ class ModuleClusters {
 
   ModuleClusters(Acts::BinUtility segmentation,
                  std::vector<Acts::BoundIndices> geoIndices, bool merge,
-                 double nsigma, bool commonCorner)
+                 double nsigma, bool commonCorner, bool alt = false)
       : m_segmentation(std::move(segmentation)),
         m_geoIndices(std::move(geoIndices)),
         m_merge(merge),
         m_nsigma(nsigma),
-        m_commonCorner(commonCorner) {}
+        m_commonCorner(commonCorner),
+        m_useInPlaceClusterization(alt) {}
 
   void add(DigitizedParameters params, simhit_t simhit);
   std::vector<std::pair<DigitizedParameters, std::set<simhit_t>>>
@@ -58,14 +59,21 @@ class ModuleClusters {
   bool m_merge;
   double m_nsigma;
   bool m_commonCorner;
+  bool m_useInPlaceClusterization;
 
   std::vector<ModuleValue> createCellCollection();
   void merge();
-  ModuleValue squash(std::vector<ModuleValue>& values);
+  void mergeWithInPlaceClusterization();
+  template <typename index_t, unsigned int AXIS>
+  void mergeWithInPlaceClusterizationImpl();
+
+  template <class cell_list_t>
+  ModuleValue squash(cell_list_t& values) const;
   std::vector<std::size_t> nonGeoEntries(
-      std::vector<Acts::BoundIndices>& indices);
+      const std::vector<Acts::BoundIndices>& indices) const;
+  template <class cell_list_t>
   std::vector<std::vector<ModuleValue>> mergeParameters(
-      std::vector<ModuleValue> values);
+      cell_list_t& values) const;
 };
 
 }  // namespace ActsExamples

Examples/Algorithms/Digitization/include/ActsExamples/Digitization/ModuleClusters.hpp

@@ -196,11 +196,13 @@ ProcessCode DigitizationAlgorithm::execute(const AlgorithmContext& ctx) const {
 
     // Run the digitizer. Iterate over the hits for this surface inside the
     // visitor so we do not need to lookup the variant object per-hit.
+
     std::visit(
         [&](const auto& digitizer) {
           ModuleClusters moduleClusters(
               digitizer.geometric.segmentation, digitizer.geometric.indices,
-              m_cfg.doMerge, m_cfg.mergeNsigma, m_cfg.mergeCommonCorner);
+              m_cfg.doMerge, m_cfg.mergeNsigma, m_cfg.mergeCommonCorner,
+              m_cfg.useInPlaceClusterization);
 
           for (auto h = moduleSimHits.begin(); h != moduleSimHits.end(); ++h) {
             const auto& simHit = *h;

Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp

@@ -9,6 +9,7 @@
 #include "ActsExamples/Digitization/ModuleClusters.hpp"
 
 #include "Acts/Clusterization/Clusterization.hpp"
+#include "Acts/Clusterization/InPlaceClusterization.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 #include "ActsExamples/Digitization/MeasurementCreation.hpp"
 #include "ActsFatras/Digitization/Channelizer.hpp"
@@ -19,6 +20,7 @@
 #include <cstdlib>
 #include <limits>
 #include <memory>
+#include <ranges>
 #include <stdexcept>
 #include <type_traits>
 
@@ -48,7 +50,12 @@ void ModuleClusters::add(DigitizedParameters params, simhit_t simhit) {
 std::vector<std::pair<DigitizedParameters, std::set<ModuleClusters::simhit_t>>>
 ModuleClusters::digitizedParameters() {
   if (m_merge) {  // (re-)build the clusters
-    merge();
+    if (m_useInPlaceClusterization && !m_moduleValues.empty() &&
+        std::holds_alternative<Cluster::Cell>(m_moduleValues.front().value)) {
+      mergeWithInPlaceClusterization();
+    } else {
+      merge();
+    }
   }
   std::vector<std::pair<DigitizedParameters, std::set<simhit_t>>> retv;
   for (ModuleValue& mval : m_moduleValues) {
@@ -96,6 +103,145 @@ std::vector<ModuleValue> ModuleClusters::createCellCollection() {
   }
   return cells;
 }
+void ModuleClusters::mergeWithInPlaceClusterization() {
+  if (m_moduleValues.size() < std::numeric_limits<std::uint16_t>::max()) {
+    mergeWithInPlaceClusterizationImpl<std::uint16_t, 0>();
+  } else {
+    mergeWithInPlaceClusterizationImpl<std::uint32_t, 0>();
+  }
+}
+
+// Helper to adapt an in-place clustered cell collection to imitate a
+// ModuleValue collection of a cluster cell_t must have an m_srcIndex member
+// which refers to a ModuleValue
+template <typename cell_t>
+struct ClusterModuleValues {
+  // iterator to iterate over the ModuleValues of a cluster
+  struct iterator {
+    ModuleValue& operator*() {
+      assert(m_moduleValues && m_iter->m_srcIndex < m_moduleValues->size());
+      return (*m_moduleValues)[m_iter->m_srcIndex];
+    }
+
+    // Test whether two iterators refer to the same ModuleValue
+    // Assumes that all ModuleValues belong to the same container, and
+    // that the underlying in-place clustered cell containers are the same.
+    bool operator==(const iterator& other) const {
+      assert(m_moduleValues == other.m_moduleValues);
+      return other.m_iter == m_iter;
+    }
+    iterator& operator++() {
+      ++m_iter;
+      return *this;
+    }
+
+    std::vector<ModuleValue>* m_moduleValues;
+    std::span<cell_t>::iterator m_iter;
+  };
+
+  // Iterator to refer to the first ModuleValue of a range
+  iterator begin() { return iterator{m_moduleValues, m_cells.begin()}; }
+  // Iterator to refer to the ModuleValue after the last one of a range
+  iterator end() { return iterator{m_moduleValues, m_cells.end()}; }
+
+  // Get the number of ModuleValues in the range
+  std::size_t size() const { return m_cells.size(); }
+
+  // Return true if the ModuleValue range is empty
+  bool empty() const { return m_cells.empty(); }
+  ModuleValue& operator[](unsigned int idx) {
+    assert(idx < m_cells.size() && m_moduleValues &&
+           m_cells[idx].m_srcIndex < m_moduleValues->size());
+    return (*m_moduleValues)[m_cells[idx].m_srcIndex];
+  }
+
+  // return the ModuleValue at a certain position in the range (range checked)
+  // @param idx the position in the range where 0 refers to the first element in the range
+  // May throw a range_error if idx is not a valid index in the range.
+  ModuleValue& at(unsigned int idx) {
+    assert(idx < m_cells.size() && m_moduleValues &&
+           m_cells[idx].m_srcIndex < m_moduleValues->size());
+    return (*m_moduleValues).at(m_cells[idx].m_srcIndex);
+  }
+
+  // return the ModuleValue at a certain position in the range (no range check)
+  // @param idx the position in the range where 0 refers to the first element in the range
+  const ModuleValue& operator[](unsigned int idx) const {
+    assert(idx < m_cells.size() && m_moduleValues &&
+           m_cells[idx].m_srcIndex < m_moduleValues->size());
+    return (*m_moduleValues)[m_cells[idx]];
+  }
+
+  std::vector<ModuleValue>* m_moduleValues;
+  std::span<cell_t> m_cells;  // A range of cells which refer to a ModuleValue
+                              // in the ModuleValue collection
+};
+
+template <typename index_t, unsigned int AXIS>
+void ModuleClusters::mergeWithInPlaceClusterizationImpl() {
+  // create a cell collection which contains the coordinates of the ModuleValues
+  // and which refer back to the source ModuleValue.
+  using Cell = Acts::InPlaceClusterization::Cell<std::int16_t, 2, index_t>;
+  std::vector<Cell> cells;
+  cells.reserve(m_moduleValues.size());
+  {
+    index_t idx = 0;
+    for (const ModuleValue& mval : m_moduleValues) {
+      if (std::holds_alternative<Cluster::Cell>(mval.value)) {
+        cells.push_back(Cell(
+            std::array<std::int16_t, 2>{
+                static_cast<std::int16_t>(
+                    std::get<ActsExamples::Cluster::Cell>(mval.value).bin[0]),
+                static_cast<std::int16_t>(
+                    std::get<ActsExamples::Cluster::Cell>(mval.value).bin[1])},
+            idx));
+      }
+      ++idx;
+    }
+  }
+  // Clusterize the cell collection in-place assuming a 4 or 8 fold connectivity
+  if (m_commonCorner) {
+    Acts::InPlaceClusterization::clusterize<
+        AXIS, std::vector<Cell>, index_t,
+        Acts::InPlaceClusterization::EConnectionType::CommonEdgeOrCorner>(
+        cells);
+  } else {
+    Acts::InPlaceClusterization::clusterize<
+        AXIS, std::vector<Cell>, index_t,
+        Acts::InPlaceClusterization::EConnectionType::CommonEdge>(cells);
+  }
+
+  std::vector<ModuleValue> newVals;
+  newVals.reserve(m_moduleValues.size());
+  // Helper to create a cluster module value for the given cell range
+  // the ModuleValues the cell range refers to will be squashed and then
+  // added to a new Cluster.
+  auto addClusters = [this, &newVals](std::vector<Cell>& all_cells,
+                                      unsigned int idx_begin,
+                                      unsigned int idx_end) {
+    auto cluster_cells =
+        std::span(all_cells.begin() + idx_begin, all_cells.begin() + idx_end);
+    if (cluster_cells.size() == 1) {
+      // no need to merge parameters
+      ClusterModuleValues<Cell> temp{&m_moduleValues, cluster_cells};
+      newVals.push_back(squash(temp));
+    } else {
+      ClusterModuleValues<Cell> temp{&m_moduleValues, cluster_cells};
+      std::ranges::sort(cluster_cells, [](const Cell& a, const Cell& b) {
+        return Acts::InPlaceClusterization::traits::getCellCoordinate(a, AXIS) <
+               Acts::InPlaceClusterization::traits::getCellCoordinate(b, AXIS);
+      });
+
+      for (std::vector<ModuleValue>& remerged : mergeParameters(temp)) {
+        newVals.push_back(squash(remerged));
+      }
+    }
+  };
+
+  // Create ModuleValues for each cell range which forms a cluster.
+  for_each_cluster(cells, addClusters);
+  m_moduleValues = std::move(newVals);
+}
 
 void ModuleClusters::merge() {
   std::vector<ModuleValue> cells = createCellCollection();
@@ -117,7 +263,6 @@ void ModuleClusters::merge() {
       // only. Still have to check if they match based on the other
       // indices (a good example of this would a for a timing
       // detector).
-
       for (std::vector<ModuleValue>& remerged : mergeParameters(cellv)) {
         newVals.push_back(squash(remerged));
       }
@@ -134,7 +279,7 @@ void ModuleClusters::merge() {
 
 // ATTN: returns vector of index into `indices'
 std::vector<std::size_t> ModuleClusters::nonGeoEntries(
-    std::vector<Acts::BoundIndices>& indices) {
+    const std::vector<Acts::BoundIndices>& indices) const {
   std::vector<std::size_t> retv;
   for (std::size_t i = 0; i < indices.size(); i++) {
     auto idx = indices.at(i);
@@ -146,8 +291,9 @@ std::vector<std::size_t> ModuleClusters::nonGeoEntries(
 }
 
 // Merging based on parameters
+template <class cell_list_t>
 std::vector<std::vector<ModuleValue>> ModuleClusters::mergeParameters(
-    std::vector<ModuleValue> values) {
+    cell_list_t& values) const {
   std::vector<std::vector<ModuleValue>> retv;
 
   std::vector<bool> used(values.size(), false);
@@ -181,7 +327,7 @@ std::vector<std::vector<ModuleValue>> ModuleClusters::mergeParameters(
       // The cluster to be merged into can have more than one
       // associated value at this point, so we have to consider them
       // all
-      for (ModuleValue& thisval : thisvec) {
+      for (const ModuleValue& thisval : thisvec) {
         // Loop over non-geometric dimensions
         for (auto k : nonGeoEntries(thisval.paramIndices)) {
           double p_i = thisval.paramValues.at(k);
@@ -221,14 +367,15 @@ std::vector<std::vector<ModuleValue>> ModuleClusters::mergeParameters(
   return retv;
 }
 
-ModuleValue ModuleClusters::squash(std::vector<ModuleValue>& values) {
+template <class cell_list_t>
+ModuleValue ModuleClusters::squash(cell_list_t& values) const {
   ModuleValue mval;
   double tot = 0;
   double tot2 = 0;
   std::vector<double> weights;
 
   // First, start by computing cell weights
-  for (ModuleValue& other : values) {
+  for (const ModuleValue& other : values) {
     if (std::holds_alternative<Cluster::Cell>(other.value)) {
       weights.push_back(std::get<Cluster::Cell>(other.value).activation);
     } else {
@@ -240,7 +387,7 @@ ModuleValue ModuleClusters::squash(std::vector<ModuleValue>& values) {
 
   // Now, go over the non-geometric indices
   for (std::size_t i = 0; i < values.size(); i++) {
-    ModuleValue& other = values.at(i);
+    const ModuleValue& other = values.at(i);
     for (std::size_t j = 0; j < other.paramIndices.size(); j++) {
       auto idx = other.paramIndices.at(j);
       if (!rangeContainsValue(m_geoIndices, idx)) {
@@ -272,7 +419,7 @@ ModuleValue ModuleClusters::squash(std::vector<ModuleValue>& values) {
   std::size_t b1max = 0;
 
   for (std::size_t i = 0; i < values.size(); i++) {
-    ModuleValue& other = values.at(i);
+    const ModuleValue& other = values.at(i);
     if (!std::holds_alternative<Cluster::Cell>(other.value)) {
       continue;
     }

Examples/Algorithms/Digitization/src/ModuleClusters.cpp

@@ -746,6 +746,7 @@ def addDigitization(
     rnd: Optional[acts.examples.RandomNumbers] = None,
     doMerge: Optional[bool] = None,
     mergeCommonCorner: Optional[bool] = None,
+    useInPlaceClusterization: Optional[bool] = None,
     minEnergyDeposit: Optional[float] = None,
     logLevel: Optional[acts.logging.Level] = None,
 ) -> acts.examples.Sequencer:
@@ -786,6 +787,7 @@ def addDigitization(
         **acts.examples.defaultKWArgs(
             doMerge=doMerge,
             mergeCommonCorner=mergeCommonCorner,
+            useInPlaceClusterization=useInPlaceClusterization,
         ),
     )
 

Examples/Python/python/acts/examples/simulation.py

@@ -56,7 +56,7 @@ void addDigitization(Context& ctx) {
         outputParticleMeasurementsMap, outputSimHitMeasurementsMap,
         surfaceByIdentifier, randomNumbers, doOutputCells, doClusterization,
         doMerge, mergeCommonCorner, minEnergyDeposit, digitizationConfigs,
-        minMaxRetries);
+        minMaxRetries, useInPlaceClusterization);
 
     c.def_readonly("mergeNsigma", &Config::mergeNsigma);