Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:12:31

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 Ccl::Label& getCellLabel(Cell2D& cell) {
0096   return cell.label;
0097 }
0098 
0099 bool operator==(const Cell2D& left, const Cell2D& right) {
0100   return left.row == right.row && left.col == right.col;
0101 }
0102 
0103 bool cellComp(const Cell2D& left, const Cell2D& right) {
0104   return (left.row == right.row) ? left.col < right.col : left.row < right.row;
0105 }
0106 
0107 struct Cluster2D {
0108   std::vector<Cell2D> cells;
0109   std::size_t hash{0};
0110 };
0111 
0112 void clusterAddCell(Cluster2D& cl, const Cell2D& cell) {
0113   cl.cells.push_back(cell);
0114 }
0115 
0116 void hash(Cluster2D& cl) {
0117   std::ranges::sort(cl.cells, cellComp);
0118   cl.hash = 0;
0119   for (const Cell2D& c : cl.cells) {
0120     boost::hash_combine(cl.hash, c.col);
0121   }
0122 }
0123 
0124 bool clHashComp(const Cluster2D& left, const Cluster2D& right) {
0125   return left.hash < right.hash;
0126 }
0127 
0128 template <typename RNG>
0129 void genclusterw(int x, int y, int x0, int y0, int x1, int y1,
0130                  std::vector<Cell2D>& cells, RNG& rng, double startp = 0.5,
0131                  double decayp = 0.9) {
0132   std::vector<Cell2D> add;
0133 
0134   auto maybe_add = [&](int x_, int y_) {
0135     Cell2D c(x_, y_);
0136     if (std::uniform_real_distribution<double>()(rng) < startp &&
0137         !rangeContainsValue(cells, c)) {
0138       cells.push_back(c);
0139       add.push_back(c);
0140     }
0141   };
0142 
0143   // NORTH
0144   if (y < y1) {
0145     maybe_add(x, y + 1);
0146   }
0147   // NORTHEAST
0148   if (x < x1 && y < y1) {
0149     maybe_add(x + 1, y + 1);
0150   }
0151   // EAST
0152   if (x < x1) {
0153     maybe_add(x + 1, y);
0154   }
0155   // SOUTHEAST
0156   if (x < x1 && y > y0) {
0157     maybe_add(x + 1, y - 1);
0158   }
0159   // SOUTH
0160   if (y > y0) {
0161     maybe_add(x, y - 1);
0162   }
0163   // SOUTHWEST
0164   if (x > x0 && y > y0) {
0165     maybe_add(x - 1, y - 1);
0166   }
0167   // WEST
0168   if (x > x0) {
0169     maybe_add(x - 1, y);
0170   }
0171   // NORTHWEST
0172   if (x > x0 && y < y1) {
0173     maybe_add(x - 1, y + 1);
0174   }
0175 
0176   for (Cell2D& c : add) {
0177     genclusterw(c.row, c.col, x0, y0, x1, y1, cells, rng, startp * decayp,
0178                 decayp);
0179   }
0180 }
0181 
0182 template <typename RNG>
0183 Cluster2D gencluster(int x0, int y0, int x1, int y1, RNG& rng,
0184                      double startp = 0.5, double decayp = 0.9) {
0185   int x0_ = x0 + 1;
0186   int x1_ = x1 - 1;
0187   int y0_ = y0 + 1;
0188   int y1_ = y1 - 1;
0189 
0190   int x = std::uniform_int_distribution<std::int32_t>(x0_, x1_)(rng);
0191   int y = std::uniform_int_distribution<std::int32_t>(y0_, y1_)(rng);
0192 
0193   std::vector<Cell2D> cells = {Cell2D(x, y)};
0194   genclusterw(x, y, x0_, y0_, x1_, y1_, cells, rng, startp, decayp);
0195 
0196   Cluster2D cl;
0197   cl.cells = std::move(cells);
0198 
0199   return cl;
0200 }
0201 
0202 BOOST_AUTO_TEST_CASE(Grid_2D_rand) {
0203   using Cell = Cell2D;
0204   using CellC = std::vector<Cell>;
0205   using Cluster = Cluster2D;
0206   using ClusterC = std::vector<Cluster>;
0207 
0208   std::size_t sizeX = 1000;
0209   std::size_t sizeY = 1000;
0210   std::size_t startSeed = 71902647;
0211   std::size_t ntries = 100;
0212 
0213   std::cout << "Grid_2D_rand test with parameters: " << std::endl;
0214   std::cout << " sizeX = " << sizeX << std::endl;
0215   std::cout << " sizeY = " << sizeY << std::endl;
0216   std::cout << " startSeed = " << startSeed << std::endl;
0217   std::cout << " ntries = " << ntries << std::endl;
0218 
0219   while (ntries-- > 0) {
0220     std::mt19937_64 rnd(startSeed++);
0221 
0222     std::vector<Cluster> cls;
0223     std::vector<Cell> cells;
0224     for (Rectangle& rect : segment(0, 0, sizeX, sizeY, rnd)) {
0225       auto& [x0, y0, x1, y1] = rect;
0226       Cluster cl = gencluster(x0, y0, x1, y1, rnd);
0227       hash(cl);
0228       cells.insert(cells.end(), cl.cells.begin(), cl.cells.end());
0229       cls.push_back(cl);
0230     }
0231 
0232     std::shuffle(cells.begin(), cells.end(), rnd);
0233 
0234     ClusterC newCls = Ccl::createClusters<CellC, ClusterC>(cells);
0235 
0236     for (Cluster& cl : newCls) {
0237       hash(cl);
0238     }
0239 
0240     std::ranges::sort(cls, clHashComp);
0241     std::ranges::sort(newCls, clHashComp);
0242 
0243     BOOST_CHECK_EQUAL(cls.size(), newCls.size());
0244     for (std::size_t i = 0; i < cls.size(); i++) {
0245       BOOST_CHECK_EQUAL(cls.at(i).hash, newCls.at(i).hash);
0246     }
0247   }
0248 }
0249 
0250 }  // namespace Acts::Test