Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SeedFinderCudaTest.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SeedFinderCudaTest.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2023 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
14 #include "Acts/Seeding/Seed.hpp"
18 
19 #include <chrono>
20 #include <fstream>
21 #include <iomanip>
22 #include <iostream>
23 #include <sstream>
24 #include <utility>
25 
26 #include <boost/type_erasure/any_cast.hpp>
27 #include <cuda_profiler_api.h>
28 
29 #include "ATLASCuts.hpp"
30 #include "SpacePoint.hpp"
31 
32 using namespace Acts::UnitLiterals;
33 
34 std::vector<const SpacePoint*> readFile(std::string filename) {
36  int layer;
37  std::vector<const SpacePoint*> readSP;
38 
39  std::ifstream spFile(filename);
40  if (spFile.is_open()) {
41  while (!spFile.eof()) {
42  std::getline(spFile, line);
43  std::stringstream ss(line);
44  std::string linetype;
45  ss >> linetype;
46  float x, y, z, r, varianceR, varianceZ;
47  if (linetype == "lxyz") {
48  ss >> layer >> x >> y >> z >> varianceR >> varianceZ;
49  r = std::hypot(x, y);
50  float f22 = varianceR;
51  float wid = varianceZ;
52  float cov = wid * wid * .08333;
53  if (cov < f22)
54  cov = f22;
55  if (std::abs(z) > 450.) {
56  varianceZ = 9. * cov;
57  varianceR = .06;
58  } else {
59  varianceR = 9. * cov;
60  varianceZ = .06;
61  }
62  SpacePoint* sp =
63  new SpacePoint{x, y, z, r, layer, varianceR, varianceZ};
64  // if(r < 200.){
65  // sp->setClusterList(1,0);
66  // }
67  readSP.push_back(sp);
68  }
69  }
70  }
71  return readSP;
72 }
73 
74 int main(int argc, char** argv) {
75  auto start_pre = std::chrono::system_clock::now();
76 
77  std::string file{"sp.txt"};
78  bool help(false);
79  bool quiet(false);
80  bool allgroup(false);
81  bool do_cpu(true);
82  int nGroupToIterate = 500;
83  int skip = 0;
84  int deviceID = 0;
85  int nTrplPerSpBLimit = 100;
86  int nAvgTrplPerSpBLimit = 2;
87 
88  int opt;
89  while ((opt = getopt(argc, argv, "haf:n:s:d:l:m:qG")) != -1) {
90  switch (opt) {
91  case 'a':
92  allgroup = true;
93  break;
94  case 'f':
95  file = optarg;
96  break;
97  case 'n':
98  nGroupToIterate = atoi(optarg);
99  break;
100  case 's':
101  skip = atoi(optarg);
102  break;
103  case 'd':
104  deviceID = atoi(optarg);
105  break;
106  case 'l':
107  nAvgTrplPerSpBLimit = atoi(optarg);
108  break;
109  case 'm':
110  nTrplPerSpBLimit = atoi(optarg);
111  break;
112  case 'q':
113  quiet = true;
114  break;
115  case 'G':
116  do_cpu = false;
117  break;
118  case 'h':
119  help = true;
120  [[fallthrough]];
121  default: /* '?' */
122  std::cerr << "Usage: " << argv[0] << " [-hq] [-f FILENAME]\n";
123  if (help) {
124  std::cout << " -h : this help" << std::endl;
125  std::cout << " -a ALL : analyze all groups. Default is \""
126  << allgroup << "\"" << std::endl;
127  std::cout
128  << " -f FILE : read spacepoints from FILE. Default is \""
129  << file << "\"" << std::endl;
130  std::cout << " -n NUM : Number of groups to iterate in seed "
131  "finding. Default is "
132  << nGroupToIterate << std::endl;
133  std::cout << " -s SKIP : Number of groups to skip in seed "
134  "finding. Default is "
135  << skip << std::endl;
136  std::cout << " -d DEVID : NVIDIA GPU device ID. Default is "
137  << deviceID << std::endl;
138  std::cout << " -l : A limit on the average number of triplets "
139  "per bottom spacepoint: this is used for determining "
140  "matrix size for triplets per middle space point"
141  << nAvgTrplPerSpBLimit << std::endl;
142  std::cout << " -m : A limit on the number of triplets per "
143  "bottom spacepoint: users do not have to touch this for "
144  "# spacepoints < ~200k"
145  << nTrplPerSpBLimit << std::endl;
146  std::cout << " -q : don't print out all found seeds"
147  << std::endl;
148  std::cout << " -G : only run on GPU, not CPU" << std::endl;
149  }
150 
151  exit(EXIT_FAILURE);
152  }
153  }
154 
155  std::string devName;
156  ACTS_CUDA_ERROR_CHECK(cudaSetDevice(deviceID));
157 
158  std::ifstream f(file);
159  if (!f.good()) {
160  std::cerr << "input file \"" << file << "\" does not exist\n";
161  exit(EXIT_FAILURE);
162  }
163 
164  std::vector<const SpacePoint*> spVec = readFile(file);
165 
166  std::cout << "read " << spVec.size() << " SP from file " << file << std::endl;
167 
168  // Set seed finder configuration
170  // silicon detector max
171  config.rMax = 160._mm;
172  config.deltaRMin = 5._mm;
173  config.deltaRMax = 160._mm;
174  config.deltaRMinTopSP = config.deltaRMin;
175  config.deltaRMinBottomSP = config.deltaRMin;
176  config.deltaRMaxTopSP = config.deltaRMax;
177  config.deltaRMaxBottomSP = config.deltaRMax;
178  config.collisionRegionMin = -250._mm;
179  config.collisionRegionMax = 250._mm;
180  config.zMin = -2800._mm;
181  config.zMax = 2800._mm;
182  config.maxSeedsPerSpM = 5;
183  // 2.7 eta
184  config.cotThetaMax = 7.40627;
185  config.sigmaScattering = 1.00000;
186 
187  config.minPt = 500._MeV;
188  config.impactMax = 10._mm;
189  config = config.toInternalUnits();
190 
192  options.bFieldInZ = 1.99724_T;
193  options.beamPos = {-.5_mm, -.5_mm};
194  options = options.toInternalUnits().calculateDerivedQuantities(config);
195 
196  int numPhiNeighbors = 1;
197 
198  // extent used to store r range for middle spacepoint
199  Acts::Extent rRangeSPExtent;
200 
201  const Acts::Range1D<float> rMiddleSPRange;
202 
203  std::vector<std::pair<int, int>> zBinNeighborsTop;
204  std::vector<std::pair<int, int>> zBinNeighborsBottom;
205 
206  // cuda
207  cudaDeviceProp prop;
208  ACTS_CUDA_ERROR_CHECK(cudaGetDeviceProperties(&prop, deviceID));
209  printf("\n GPU Device %d: \"%s\" with compute capability %d.%d\n\n", deviceID,
210  prop.name, prop.major, prop.minor);
211  config.maxBlockSize = prop.maxThreadsPerBlock;
212  config.nTrplPerSpBLimit = nTrplPerSpBLimit;
213  config.nAvgTrplPerSpBLimit = nAvgTrplPerSpBLimit;
214 
215  // binfinder
216  auto bottomBinFinder = std::make_shared<Acts::BinFinder<SpacePoint>>(
217  Acts::BinFinder<SpacePoint>(zBinNeighborsBottom, numPhiNeighbors));
218  auto topBinFinder = std::make_shared<Acts::BinFinder<SpacePoint>>(
219  Acts::BinFinder<SpacePoint>(zBinNeighborsTop, numPhiNeighbors));
220  Acts::SeedFilterConfig sfconf;
222  config.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
223  Acts::SeedFilter<SpacePoint>(sfconf, &atlasCuts));
224  Acts::SeedFinder<SpacePoint> seedFinder_cpu(config);
225  Acts::SeedFinder<SpacePoint, Acts::Cuda> seedFinder_cuda(config, options);
226 
227  // covariance tool, sets covariances per spacepoint as required
228  auto ct = [=](const SpacePoint& sp, float, float,
229  float) -> std::pair<Acts::Vector3, Acts::Vector2> {
230  Acts::Vector3 position(sp.x(), sp.y(), sp.z());
231  Acts::Vector2 variance(sp.varianceR, sp.varianceZ);
232  return std::make_pair(position, variance);
233  };
234 
235  // setup spacepoint grid config
237  gridConf.minPt = config.minPt;
238  gridConf.rMax = config.rMax;
239  gridConf.zMax = config.zMax;
240  gridConf.zMin = config.zMin;
241  gridConf.deltaRMax = config.deltaRMax;
242  gridConf.cotThetaMax = config.cotThetaMax;
243  // setup spacepoint grid options
245  gridOpts.bFieldInZ = options.bFieldInZ;
246  // create grid with bin sizes according to the configured geometry
247  std::unique_ptr<Acts::SpacePointGrid<SpacePoint>> grid =
248  Acts::SpacePointGridCreator::createGrid<SpacePoint>(gridConf, gridOpts);
249  auto spGroup = Acts::BinnedSPGroup<SpacePoint>(
250  spVec.begin(), spVec.end(), ct, bottomBinFinder, topBinFinder,
251  std::move(grid), rRangeSPExtent, config, options);
252 
253  auto end_pre = std::chrono::system_clock::now();
254  std::chrono::duration<double> elapsec_pre = end_pre - start_pre;
255  double preprocessTime = elapsec_pre.count();
256  std::cout << "Preprocess Time: " << preprocessTime << std::endl;
257 
258  //--------------------------------------------------------------------//
259  // Begin Seed finding //
260  //--------------------------------------------------------------------//
261 
262  auto start_cpu = std::chrono::system_clock::now();
263 
264  int group_count;
265  auto groupIt = Acts::BinnedSPGroupIterator<SpacePoint>(spGroup, skip);
266 
267  //----------- CPU ----------//
268  group_count = 0;
269  std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_cpu;
270 
271  if (do_cpu) {
272  decltype(seedFinder_cpu)::SeedingState state;
273  state.spacePointData.resize(spVec.size());
274  for (; groupIt != spGroup.end(); ++groupIt) {
275  const auto [bottom, middle, top] = *groupIt;
276  seedFinder_cpu.createSeedsForGroup(
277  options, state, spGroup.grid(),
278  std::back_inserter(seedVector_cpu.emplace_back()), bottom, middle,
279  top, rMiddleSPRange);
280  group_count++;
281  if (allgroup == false) {
282  if (group_count >= nGroupToIterate)
283  break;
284  }
285  }
286  // auto timeMetric_cpu = seedFinder_cpu.getTimeMetric();
287  std::cout << "Analyzed " << group_count << " groups for CPU" << std::endl;
288  }
289 
290  auto end_cpu = std::chrono::system_clock::now();
291  std::chrono::duration<double> elapsec_cpu = end_cpu - start_cpu;
292  double cpuTime = elapsec_cpu.count();
293 
294  //----------- CUDA ----------//
295 
296  cudaProfilerStart();
297  auto start_cuda = std::chrono::system_clock::now();
298 
299  group_count = 0;
300  std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_cuda;
301  groupIt = Acts::BinnedSPGroupIterator<SpacePoint>(spGroup, skip);
302 
303  Acts::SpacePointData spacePointData;
304  spacePointData.resize(spVec.size());
305 
306  for (; groupIt != spGroup.end(); ++groupIt) {
307  const auto [bottom, middle, top] = *groupIt;
308  seedVector_cuda.push_back(seedFinder_cuda.createSeedsForGroup(
309  spacePointData, spGroup.grid(), bottom, middle, top));
310  group_count++;
311  if (allgroup == false) {
312  if (group_count >= nGroupToIterate)
313  break;
314  }
315  }
316 
317  auto end_cuda = std::chrono::system_clock::now();
318  std::chrono::duration<double> elapsec_cuda = end_cuda - start_cuda;
319  double cudaTime = elapsec_cuda.count();
320 
321  cudaProfilerStop();
322  std::cout << "Analyzed " << group_count << " groups for CUDA" << std::endl;
323 
324  std::cout << std::endl;
325  std::cout << "----------------------- Time Metric -----------------------"
326  << std::endl;
327  std::cout << " " << (do_cpu ? "CPU" : " ")
328  << " CUDA " << (do_cpu ? "Speedup " : "")
329  << std::endl;
330  std::cout << "Seedfinding_Time " << std::setw(11)
331  << (do_cpu ? std::to_string(cpuTime) : "") << " " << std::setw(11)
332  << cudaTime << " " << std::setw(11)
333  << (do_cpu ? std::to_string(cpuTime / cudaTime) : "") << std::endl;
334  double wallTime_cpu = cpuTime + preprocessTime;
335  double wallTime_cuda = cudaTime + preprocessTime;
336  std::cout << "Wall_time " << std::setw(11)
337  << (do_cpu ? std::to_string(wallTime_cpu) : "") << " "
338  << std::setw(11) << wallTime_cuda << " " << std::setw(11)
339  << (do_cpu ? std::to_string(wallTime_cpu / wallTime_cuda) : "")
340  << std::endl;
341  std::cout << "-----------------------------------------------------------"
342  << std::endl;
343  std::cout << std::endl;
344 
345  int nSeed_cpu = 0;
346  for (auto& outVec : seedVector_cpu) {
347  nSeed_cpu += outVec.size();
348  }
349 
350  int nSeed_cuda = 0;
351  for (auto& outVec : seedVector_cuda) {
352  nSeed_cuda += outVec.size();
353  }
354 
355  std::cout << "Number of Seeds (CPU | CUDA): " << nSeed_cpu << " | "
356  << nSeed_cuda << std::endl;
357 
358  int nMatch = 0;
359 
360  for (size_t i = 0; i < seedVector_cpu.size(); i++) {
361  auto regionVec_cpu = seedVector_cpu[i];
362  auto regionVec_cuda = seedVector_cuda[i];
363 
364  std::vector<std::vector<SpacePoint>> seeds_cpu;
365  std::vector<std::vector<SpacePoint>> seeds_cuda;
366 
367  // for (size_t i_cpu = 0; i_cpu < regionVec_cpu.size(); i_cpu++) {
368  for (auto sd : regionVec_cpu) {
369  std::vector<SpacePoint> seed_cpu;
370  seed_cpu.push_back(*(sd.sp()[0]));
371  seed_cpu.push_back(*(sd.sp()[1]));
372  seed_cpu.push_back(*(sd.sp()[2]));
373 
374  seeds_cpu.push_back(seed_cpu);
375  }
376 
377  for (auto sd : regionVec_cuda) {
378  std::vector<SpacePoint> seed_cuda;
379  seed_cuda.push_back(*(sd.sp()[0]));
380  seed_cuda.push_back(*(sd.sp()[1]));
381  seed_cuda.push_back(*(sd.sp()[2]));
382 
383  seeds_cuda.push_back(seed_cuda);
384  }
385 
386  for (auto seed : seeds_cpu) {
387  for (auto other : seeds_cuda) {
388  if (seed[0] == other[0] && seed[1] == other[1] && seed[2] == other[2]) {
389  nMatch++;
390  break;
391  }
392  }
393  }
394  }
395 
396  if (do_cpu) {
397  std::cout << nMatch << " seeds are matched" << std::endl;
398  std::cout << "Matching rate: " << float(nMatch) / nSeed_cpu * 100 << "%"
399  << std::endl;
400  }
401 
402  if (!quiet) {
403  if (do_cpu) {
404  std::cout << "CPU Seed result:" << std::endl;
405 
406  for (auto& regionVec : seedVector_cpu) {
407  for (size_t i = 0; i < regionVec.size(); i++) {
408  const Acts::Seed<SpacePoint>* seed = &regionVec[i];
409  const SpacePoint* sp = seed->sp()[0];
410  std::cout << " (" << sp->x() << ", " << sp->y() << ", " << sp->z()
411  << ") ";
412  sp = seed->sp()[1];
413  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
414  << sp->z() << ") ";
415  sp = seed->sp()[2];
416  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
417  << sp->z() << ") ";
418  std::cout << std::endl;
419  }
420  }
421 
422  std::cout << std::endl;
423  }
424  std::cout << "CUDA Seed result:" << std::endl;
425 
426  for (auto& regionVec : seedVector_cuda) {
427  for (size_t i = 0; i < regionVec.size(); i++) {
428  const Acts::Seed<SpacePoint>* seed = &regionVec[i];
429  const SpacePoint* sp = seed->sp()[0];
430  std::cout << " (" << sp->x() << ", " << sp->y() << ", " << sp->z()
431  << ") ";
432  sp = seed->sp()[1];
433  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
434  << sp->z() << ") ";
435  sp = seed->sp()[2];
436  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
437  << sp->z() << ") ";
438  std::cout << std::endl;
439  }
440  }
441  }
442 
443  std::cout << std::endl;
444  std::cout << std::endl;
445 
446  return 0;
447 }