Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-15 08:15:07

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 #include <boost/test/unit_test.hpp>
0010 
0011 #include "Acts/Clusterization/Clusterization.hpp"
0012 #include "Acts/Utilities/Helpers.hpp"
0013 
0014 #include <algorithm>
0015 #include <array>
0016 #include <cstdlib>
0017 #include <iostream>
0018 #include <iterator>
0019 #include <memory>
0020 #include <random>
0021 #include <stdexcept>
0022 #include <utility>
0023 #include <vector>
0024 
0025 #include <boost/functional/hash.hpp>
0026 
0027 namespace Acts::Test {
0028 
0029 using Rectangle = std::array<int, 4>;
0030 
0031 std::vector<Rectangle> concat(std::vector<std::vector<Rectangle>> vecs) {
0032   std::vector<Rectangle> flat;
0033   for (std::vector<Rectangle>& v : vecs) {
0034     if (flat.empty()) {
0035       flat = std::move(v);
0036     } else {
0037       flat.reserve(flat.size() + v.size());
0038       std::move(v.begin(), v.end(), std::back_inserter(flat));
0039       v.clear();
0040     }
0041   }
0042   return flat;
0043 }
0044 
0045 template <typename RNG>
0046 std::vector<Rectangle> segment(int x0, int y0, int x1, int y1, RNG& rng) {
0047   // Leave enough space for one-cell buffer around clusters
0048   int xmin = x0 + 3;
0049   int xmax = x1 - 3;
0050   int ymin = y0 + 3;
0051   int ymax = y1 - 3;
0052 
0053   // terminal case 1
0054   if (xmax < xmin || ymax < ymin) {
0055     return {{x0, y0, x1, y1}};
0056   }
0057 
0058   std::bernoulli_distribution cointoss;
0059   bool splitx = cointoss(rng);
0060   bool splity = cointoss(rng);
0061 
0062   // terminal case 2
0063   if (!(splitx || splity)) {
0064     return {{x0, y0, x1, y1}};
0065   }
0066 
0067   int x_ = std::uniform_int_distribution<std::int32_t>(xmin, xmax)(rng);
0068   int y_ = std::uniform_int_distribution<std::int32_t>(ymin, ymax)(rng);
0069 
0070   if (splitx && !splity) {
0071     return concat({segment(x0, y0, x_, y1, rng), segment(x_, y0, x1, y1, rng)});
0072   } else if (!splitx && splity) {
0073     return concat({segment(x0, y0, x1, y_, rng), segment(x0, y_, x1, y1, rng)});
0074   } else if (splitx && splity) {
0075     return concat({segment(x0, y0, x_, y_, rng), segment(x_, y0, x1, y_, rng),
0076                    segment(x0, y_, x_, y1, rng), segment(x_, y_, x1, y1, rng)});
0077   }
0078   throw std::runtime_error("unreachable");
0079 }
0080 
0081 struct Cell2D {
0082   Cell2D(int rowv, int colv) : row(rowv), col(colv) {}
0083   int row, col;
0084   Ccl::Label label{Ccl::NO_LABEL};
0085 };
0086 
0087 int getCellRow(const Cell2D& cell) {
0088   return cell.row;
0089 }
0090 
0091 int getCellColumn(const Cell2D& cell) {
0092   return cell.col;
0093 }
0094 
0095 bool operator==(const Cell2D& left, const Cell2D& right) {
0096   return left.row == right.row && left.col == right.col;
0097 }
0098 
0099 bool cellComp(const Cell2D& left, const Cell2D& right) {
0100   return (left.row == right.row) ? left.col < right.col : left.row < right.row;
0101 }
0102 
0103 struct Cluster2D {
0104   std::vector<Cell2D> cells;
0105   std::size_t hash{0};
0106 };
0107 
0108 void clusterAddCell(Cluster2D& cl, const Cell2D& cell) {
0109   cl.cells.push_back(cell);
0110 }
0111 
0112 void hash(Cluster2D& cl) {
0113   std::ranges::sort(cl.cells, cellComp);
0114   cl.hash = 0;
0115   for (const Cell2D& c : cl.cells) {
0116     boost::hash_combine(cl.hash, c.col);
0117   }
0118 }
0119 
0120 bool clHashComp(const Cluster2D& left, const Cluster2D& right) {
0121   return left.hash < right.hash;
0122 }
0123 
0124 template <typename RNG>
0125 void genclusterw(int x, int y, int x0, int y0, int x1, int y1,
0126                  std::vector<Cell2D>& cells, RNG& rng, double startp = 0.5,
0127                  double decayp = 0.9) {
0128   std::vector<Cell2D> add;
0129 
0130   auto maybe_add = [&](int x_, int y_) {
0131     Cell2D c(x_, y_);
0132     if (std::uniform_real_distribution<double>()(rng) < startp &&
0133         !rangeContainsValue(cells, c)) {
0134       cells.push_back(c);
0135       add.push_back(c);
0136     }
0137   };
0138 
0139   // NORTH
0140   if (y < y1) {
0141     maybe_add(x, y + 1);
0142   }
0143   // NORTHEAST
0144   if (x < x1 && y < y1) {
0145     maybe_add(x + 1, y + 1);
0146   }
0147   // EAST
0148   if (x < x1) {
0149     maybe_add(x + 1, y);
0150   }
0151   // SOUTHEAST
0152   if (x < x1 && y > y0) {
0153     maybe_add(x + 1, y - 1);
0154   }
0155   // SOUTH
0156   if (y > y0) {
0157     maybe_add(x, y - 1);
0158   }
0159   // SOUTHWEST
0160   if (x > x0 && y > y0) {
0161     maybe_add(x - 1, y - 1);
0162   }
0163   // WEST
0164   if (x > x0) {
0165     maybe_add(x - 1, y);
0166   }
0167   // NORTHWEST
0168   if (x > x0 && y < y1) {
0169     maybe_add(x - 1, y + 1);
0170   }
0171 
0172   for (Cell2D& c : add) {
0173     genclusterw(c.row, c.col, x0, y0, x1, y1, cells, rng, startp * decayp,
0174                 decayp);
0175   }
0176 }
0177 
0178 template <typename RNG>
0179 Cluster2D gencluster(int x0, int y0, int x1, int y1, RNG& rng,
0180                      double startp = 0.5, double decayp = 0.9) {
0181   int x0_ = x0 + 1;
0182   int x1_ = x1 - 1;
0183   int y0_ = y0 + 1;
0184   int y1_ = y1 - 1;
0185 
0186   int x = std::uniform_int_distribution<std::int32_t>(x0_, x1_)(rng);
0187   int y = std::uniform_int_distribution<std::int32_t>(y0_, y1_)(rng);
0188 
0189   std::vector<Cell2D> cells = {Cell2D(x, y)};
0190   genclusterw(x, y, x0_, y0_, x1_, y1_, cells, rng, startp, decayp);
0191 
0192   Cluster2D cl;
0193   cl.cells = std::move(cells);
0194 
0195   return cl;
0196 }
0197 
0198 BOOST_AUTO_TEST_CASE(Grid_2D_rand) {
0199   using Cell = Cell2D;
0200   using CellC = std::vector<Cell>;
0201   using Cluster = Cluster2D;
0202   using ClusterC = std::vector<Cluster>;
0203 
0204   std::size_t sizeX = 1000;
0205   std::size_t sizeY = 1000;
0206   std::size_t startSeed = 71902647;
0207   std::size_t ntries = 100;
0208 
0209   std::cout << "Grid_2D_rand test with parameters: " << std::endl;
0210   std::cout << " sizeX = " << sizeX << std::endl;
0211   std::cout << " sizeY = " << sizeY << std::endl;
0212   std::cout << " startSeed = " << startSeed << std::endl;
0213   std::cout << " ntries = " << ntries << std::endl;
0214 
0215   while (ntries-- > 0) {
0216     std::mt19937_64 rnd(startSeed++);
0217 
0218     std::vector<Cluster> cls;
0219     std::vector<Cell> cells;
0220     for (Rectangle& rect : segment(0, 0, sizeX, sizeY, rnd)) {
0221       auto& [x0, y0, x1, y1] = rect;
0222       Cluster cl = gencluster(x0, y0, x1, y1, rnd);
0223       hash(cl);
0224       cells.insert(cells.end(), cl.cells.begin(), cl.cells.end());
0225       cls.push_back(cl);
0226     }
0227 
0228     std::shuffle(cells.begin(), cells.end(), rnd);
0229 
0230     Acts::Ccl::ClusteringData data;
0231     ClusterC newCls;
0232     Ccl::createClusters<CellC, ClusterC>(data, cells, newCls);
0233 
0234     for (Cluster& cl : newCls) {
0235       hash(cl);
0236     }
0237 
0238     std::ranges::sort(cls, clHashComp);
0239     std::ranges::sort(newCls, clHashComp);
0240 
0241     BOOST_CHECK_EQUAL(cls.size(), newCls.size());
0242     for (std::size_t i = 0; i < cls.size(); i++) {
0243       BOOST_CHECK_EQUAL(cls.at(i).hash, newCls.at(i).hash);
0244     }
0245   }
0246 }
0247 
0248 }  // namespace Acts::Test