Acts/Clusterization/Clusterization.ipp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Clusterization/Clusterization.hpp"
0012
0013 #include <algorithm>
0014 #include <array>
0015 #include <ranges>
0016 #include <vector>
0017
0018 namespace Acts::Ccl {
0019
0020 template <typename Cluster>
0021 void reserve(Cluster& cl, std::size_t n) {
0022   if constexpr (Acts::Ccl::CanReserve<Cluster>) {
0023     clusterReserve(cl, n);
0024   }
0025 }
0026
0027 template <typename Cell, std::size_t GridDim>
0028 struct Compare {
0029   static_assert(GridDim != 1 && GridDim != 2,
0030                 "Only grid dimensions of 1 or 2 are supported");
0031 };
0032
0033 // Comparator function object for cells, column-wise ordering
0034 // Specialization for 1-D grids
0035 template <Acts::Ccl::HasRetrievableColumnInfo Cell>
0036 struct Compare<Cell, 1> {
0037   bool operator()(const Cell& c0, const Cell& c1) const {
0038     int col0 = getCellColumn(c0);
0039     int col1 = getCellColumn(c1);
0040     return col0 < col1;
0041   }
0042 };
0043
0044 // Specialization for 2-D grid
0045 template <typename Cell>
0046   requires(Acts::Ccl::HasRetrievableColumnInfo<Cell> &&
0047            Acts::Ccl::HasRetrievableRowInfo<Cell>)
0048 struct Compare<Cell, 2> {
0049   bool operator()(const Cell& c0, const Cell& c1) const {
0050     int row0 = getCellRow(c0);
0051     int row1 = getCellRow(c1);
0052     int col0 = getCellColumn(c0);
0053     int col1 = getCellColumn(c1);
0054     return (col0 == col1) ? row0 < row1 : col0 < col1;
0055   }
0056 };
0057
0058 template <std::size_t BufSize>
0059 struct ConnectionsBase {
0060   std::size_t nconn{0};
0061   std::array<Label, BufSize> buf;
0062   ConnectionsBase() { std::ranges::fill(buf, NO_LABEL); }
0063 };
0064
0065 template <std::size_t GridDim>
0066 class Connections {};
0067
0068 // On 1-D grid, cells have 1 backward space neighbor, but there can be up to 2
0069 // cells with different times in that space that can connect
0070 template <>
0071 struct Connections<1> : public ConnectionsBase<2> {
0072   using ConnectionsBase::ConnectionsBase;
0073 };
0074
0075 // On a 2-D grid, cells have 4 backward space neighbors, but there can be up to
0076 // 2 cells with different times in each of those spaces that can connect
0077 template <>
0078 struct Connections<2> : public ConnectionsBase<8> {
0079   using ConnectionsBase::ConnectionsBase;
0080 };
0081
0082 // Cell collection logic
0083 template <typename Cell, typename Connect, std::size_t GridDim>
0084 Connections<GridDim> getConnections(std::size_t idx, std::vector<Cell>& cells,
0085                                     std::vector<Label>& labels,
0086                                     Connect&& connect) {
0087   Connections<GridDim> seen;
0088
0089   for (std::size_t i = 0; i < idx; ++i) {
0090     std::size_t idx2 = idx - i - 1;
0091     ConnectResult cr = connect(cells[idx], cells[idx2]);
0092
0093     if (cr == ConnectResult::eDuplicate) {
0094       throw std::invalid_argument(
0095           "Clusterization: input contains duplicate cells");
0096     }
0097     if (cr == ConnectResult::eNoConnStop) {
0098       break;
0099     }
0100     if (cr == ConnectResult::eNoConn) {
0101       continue;
0102     }
0103     if (cr == ConnectResult::eConn) {
0104       seen.buf[seen.nconn] = labels[idx2];
0105       seen.nconn += 1;
0106       if (seen.nconn == seen.buf.size()) {
0107         break;
0108       }
0109     }
0110   }
0111
0112   return seen;
0113 }
0114
0115 template <typename CellCollection, typename ClusterCollection>
0116   requires(Acts::Ccl::CanAcceptCell<typename CellCollection::value_type,
0117                                     typename ClusterCollection::value_type>)
0118 void mergeClusters(Acts::Ccl::ClusteringData& data, const CellCollection& cells,
0119                    ClusterCollection& outv) {
0120   using Cluster = typename ClusterCollection::value_type;
0121
0122   // Accumulate clusters into the output collection
0123   std::size_t previousSize = outv.size();
0124   outv.resize(previousSize + data.nClusters.size());
0125   for (std::size_t i = 0; i < data.nClusters.size(); ++i) {
0126     Acts::Ccl::reserve(outv[previousSize + i], data.nClusters[i]);
0127   }
0128
0129   // Fill clusters with cells
0130   // We are not using enumerate, since that is less optimal than
0131   // this loop
0132   for (std::size_t i = 0; i < cells.size(); ++i) {
0133     Label label = data.labels[i] - 1;
0134     Cluster& cl = outv[previousSize + label];
0135     clusterAddCell(cl, cells[i]);
0136   }
0137
0138   // Due to previous merging, we may have now clusters with
0139   // no cells. We need to remove them
0140   std::size_t invalidClusters = 0ul;
0141   for (std::size_t i = 0; i < data.nClusters.size(); ++i) {
0142     std::size_t idx = data.nClusters.size() - i - 1;
0143     if (data.nClusters[idx] != 0) {
0144       continue;
0145     }
0146     // we have an invalid cluster.
0147     // move them all to the back so that we can remove
0148     // them later
0149     std::swap(outv[previousSize + idx],
0150               outv[outv.size() - invalidClusters - 1]);
0151     ++invalidClusters;
0152   }
0153   outv.resize(outv.size() - invalidClusters);
0154 }
0155
0156 template <typename Cell>
0157   requires(Acts::Ccl::HasRetrievableColumnInfo<Cell> &&
0158            Acts::Ccl::HasRetrievableRowInfo<Cell>)
0159 ConnectResult Connect2D<Cell>::operator()(const Cell& ref,
0160                                           const Cell& iter) const {
0161   int deltaRow = getCellRow(iter) - getCellRow(ref);
0162   int deltaCol = getCellColumn(iter) - getCellColumn(ref);
0163   assert((deltaCol < 0 || (deltaCol == 0 && deltaRow <= 0)) &&
0164          "Not iterating backwards");
0165
0166   switch (deltaCol) {
0167     case 0:
0168       if (deltaRow == 0) {
0169         return ConnectResult::eDuplicate;
0170       } else if (deltaRow == -1) {
0171         return ConnectResult::eConn;
0172       } else {
0173         return ConnectResult::eNoConn;
0174       }
0175     case -1:
0176       if (deltaRow > static_cast<int>(conn8)) {
0177         return ConnectResult::eNoConn;
0178       } else if (deltaRow < -static_cast<int>(conn8)) {
0179         return ConnectResult::eNoConnStop;
0180       } else {
0181         return ConnectResult::eConn;
0182       }
0183     default:
0184       return ConnectResult::eNoConnStop;
0185   }
0186 }
0187
0188 template <Acts::Ccl::HasRetrievableColumnInfo Cell>
0189 ConnectResult Connect1D<Cell>::operator()(const Cell& ref,
0190                                           const Cell& iter) const {
0191   int deltaCol = getCellColumn(iter) - getCellColumn(ref);
0192   assert((deltaCol <= 0) && "Not iterating backwards");
0193
0194   switch (deltaCol) {
0195     case 0:
0196       return ConnectResult::eDuplicate;
0197     case -1:
0198       return ConnectResult::eConn;
0199     default:
0200       return ConnectResult::eNoConnStop;
0201   }
0202 }
0203
0204 template <std::size_t GridDim>
0205 void recordEquivalences(const Connections<GridDim> seen, DisjointSets& ds) {
0206   // Sanity check: first element should always have
0207   // label if nconn > 0
0208   if (seen.nconn > 0 && seen.buf[0] == NO_LABEL) {
0209     throw std::logic_error("seen.nconn > 0 but seen.buf[0] == NO_LABEL");
0210   }
0211   for (std::size_t i = 1; i < seen.nconn; i++) {
0212     // Sanity check: since connection lookup is always backward
0213     // while iteration is forward, all connected cells found here
0214     // should have a label
0215     if (seen.buf[i] == NO_LABEL) {
0216       throw std::logic_error("i < seen.nconn but see.buf[i] == NO_LABEL");
0217     }
0218     // Only record equivalence if needed
0219     if (seen.buf[0] != seen.buf[i]) {
0220       ds.unionSet(seen.buf[0], seen.buf[i]);
0221     }
0222   }
0223 }
0224
0225 template <typename CellCollection, std::size_t GridDim, typename Connect>
0226 void labelClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0227                    Connect&& connect) {
0228   using Cell = typename CellCollection::value_type;
0229
0230   data.labels.resize(cells.size(), NO_LABEL);
0231   // Sort cells by position to enable in-order scan
0232   std::ranges::sort(cells, Acts::Ccl::Compare<Cell, GridDim>());
0233
0234   // First pass: Allocate labels and record equivalences
0235   for (std::size_t nCell(0ul); nCell < cells.size(); ++nCell) {
0236     const Acts::Ccl::Connections<GridDim> seen =
0237         Acts::Ccl::getConnections<Cell, Connect, GridDim>(
0238             nCell, cells, data.labels, std::forward<Connect>(connect));
0239
0240     if (seen.nconn == 0) {
0241       // Allocate new label
0242       data.labels[nCell] = data.ds.makeSet();
0243     } else {
0244       recordEquivalences(seen, data.ds);
0245       // Set label for current cell
0246       data.labels[nCell] = seen.buf[0];
0247     }
0248   }  // loop on cells
0249
0250   // Second pass: Merge labels based on recorded equivalences
0251   int maxNClusters = 0;
0252   for (Label& lbl : data.labels) {
0253     lbl = data.ds.findSet(lbl);
0254     maxNClusters = std::max(maxNClusters, lbl);
0255   }
0256
0257   // Third pass: Keep count of how many cells go in each
0258   // to-be-created clusters
0259   data.nClusters.resize(maxNClusters, 0);
0260   for (const Label label : data.labels) {
0261     ++data.nClusters[label - 1];
0262   }
0263 }
0264
0265 template <typename CellCollection, typename ClusterCollection,
0266           std::size_t GridDim, typename Connect>
0267   requires(GridDim == 1 || GridDim == 2)
0268 void createClusters(Acts::Ccl::ClusteringData& data, CellCollection& cells,
0269                     ClusterCollection& clusters, Connect&& connect) {
0270   if (cells.empty()) {
0271     return;
0272   }
0273   data.clear();
0274
0275   Acts::Ccl::labelClusters<CellCollection, GridDim, Connect>(
0276       data, cells, std::forward<Connect>(connect));
0277   Acts::Ccl::mergeClusters<CellCollection, ClusterCollection>(data, cells,
0278                                                               clusters);
0279 }
0280
0281 }  // namespace Acts::Ccl