EIC 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) 2020 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 
9 #include "Acts/Plugins/Cuda/Seeding/Seedfinder.hpp"
10 #include "Acts/Seeding/BinFinder.hpp"
11 #include "Acts/Seeding/BinnedSPGroup.hpp"
12 #include "Acts/Seeding/InternalSeed.hpp"
13 #include "Acts/Seeding/InternalSpacePoint.hpp"
14 #include "Acts/Seeding/Seed.hpp"
15 #include "Acts/Seeding/SeedFilter.hpp"
16 #include "Acts/Seeding/Seedfinder.hpp"
17 #include "Acts/Seeding/SpacePointGrid.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 std::vector<const SpacePoint*> readFile(std::string filename) {
33  std::string line;
34  int layer;
35  std::vector<const SpacePoint*> readSP;
36 
37  std::ifstream spFile(filename);
38  if (spFile.is_open()) {
39  while (!spFile.eof()) {
40  std::getline(spFile, line);
41  std::stringstream ss(line);
42  std::string linetype;
43  ss >> linetype;
44  float x, y, z, r, varianceR, varianceZ;
45  if (linetype == "lxyz") {
46  ss >> layer >> x >> y >> z >> varianceR >> varianceZ;
47  r = std::sqrt(x * x + y * y);
48  float f22 = varianceR;
49  float wid = varianceZ;
50  float cov = wid * wid * .08333;
51  if (cov < f22)
52  cov = f22;
53  if (std::abs(z) > 450.) {
54  varianceZ = 9. * cov;
55  varianceR = .06;
56  } else {
57  varianceR = 9. * cov;
58  varianceZ = .06;
59  }
60  SpacePoint* sp =
61  new SpacePoint{x, y, z, r, layer, varianceR, varianceZ};
62  // if(r < 200.){
63  // sp->setClusterList(1,0);
64  // }
65  readSP.push_back(sp);
66  }
67  }
68  }
69  return readSP;
70 }
71 
72 int main(int argc, char** argv) {
73  auto start_pre = std::chrono::system_clock::now();
74 
75  std::string file{"sp.txt"};
76  bool help(false);
77  bool quiet(false);
78  bool allgroup(false);
79  bool do_cpu(true);
80  int nGroupToIterate = 500;
81  int skip = 0;
82  int deviceID = 0;
83  int nTrplPerSpBLimit = 100;
84  int nAvgTrplPerSpBLimit = 2;
85 
86  int opt;
87  while ((opt = getopt(argc, argv, "haf:n:s:d:l:m:qG")) != -1) {
88  switch (opt) {
89  case 'a':
90  allgroup = true;
91  break;
92  case 'f':
93  file = optarg;
94  break;
95  case 'n':
96  nGroupToIterate = atoi(optarg);
97  break;
98  case 's':
99  skip = atoi(optarg);
100  break;
101  case 'd':
102  deviceID = atoi(optarg);
103  break;
104  case 'l':
105  nAvgTrplPerSpBLimit = atoi(optarg);
106  break;
107  case 'm':
108  nTrplPerSpBLimit = atoi(optarg);
109  break;
110  case 'q':
111  quiet = true;
112  break;
113  case 'G':
114  do_cpu = false;
115  break;
116  case 'h':
117  help = true;
118  [[fallthrough]];
119  default: /* '?' */
120  std::cerr << "Usage: " << argv[0] << " [-hq] [-f FILENAME]\n";
121  if (help) {
122  std::cout << " -h : this help" << std::endl;
123  std::cout << " -a ALL : analyze all groups. Default is \""
124  << allgroup << "\"" << std::endl;
125  std::cout
126  << " -f FILE : read spacepoints from FILE. Default is \""
127  << file << "\"" << std::endl;
128  std::cout << " -n NUM : Number of groups to iterate in seed "
129  "finding. Default is "
130  << nGroupToIterate << std::endl;
131  std::cout << " -s SKIP : Number of groups to skip in seed "
132  "finding. Default is "
133  << skip << std::endl;
134  std::cout << " -d DEVID : NVIDIA GPU device ID. Default is "
135  << deviceID << std::endl;
136  std::cout << " -l : A limit on the average number of triplets "
137  "per bottom spacepoint: this is used for determining "
138  "matrix size for triplets per middle space point"
139  << nAvgTrplPerSpBLimit << std::endl;
140  std::cout << " -m : A limit on the number of triplets per "
141  "bottom spacepoint: users do not have to touch this for "
142  "# spacepoints < ~200k"
143  << nTrplPerSpBLimit << std::endl;
144  std::cout << " -q : don't print out all found seeds"
145  << std::endl;
146  std::cout << " -G : only run on GPU, not CPU" << std::endl;
147  }
148 
149  exit(EXIT_FAILURE);
150  }
151  }
152 
153  std::string devName;
154  ACTS_CUDA_ERROR_CHECK(cudaSetDevice(deviceID));
155 
156  std::ifstream f(file);
157  if (!f.good()) {
158  std::cerr << "input file \"" << file << "\" does not exist\n";
159  exit(EXIT_FAILURE);
160  }
161 
162  std::vector<const SpacePoint*> spVec = readFile(file);
163 
164  std::cout << "read " << spVec.size() << " SP from file " << file << std::endl;
165 
166  // Set seed finder configuration
167  Acts::SeedfinderConfig<SpacePoint> config;
168  // silicon detector max
169  config.rMax = 160.;
170  config.deltaRMin = 5.;
171  config.deltaRMax = 160.;
172  config.collisionRegionMin = -250.;
173  config.collisionRegionMax = 250.;
174  config.zMin = -2800.;
175  config.zMax = 2800.;
176  config.maxSeedsPerSpM = 5;
177  // 2.7 eta
178  config.cotThetaMax = 7.40627;
179  config.sigmaScattering = 1.00000;
180 
181  config.minPt = 500.;
182  config.bFieldInZ = 0.00199724;
183 
184  config.beamPos = {-.5, -.5};
185  config.impactMax = 10.;
186 
187  // cuda
188  cudaDeviceProp prop;
189  ACTS_CUDA_ERROR_CHECK(cudaGetDeviceProperties(&prop, deviceID));
190  printf("\n GPU Device %d: \"%s\" with compute capability %d.%d\n\n", deviceID,
191  prop.name, prop.major, prop.minor);
192  config.maxBlockSize = prop.maxThreadsPerBlock;
193  config.nTrplPerSpBLimit = nTrplPerSpBLimit;
194  config.nAvgTrplPerSpBLimit = nAvgTrplPerSpBLimit;
195 
196  // binfinder
197  auto bottomBinFinder = std::make_shared<Acts::BinFinder<SpacePoint>>(
198  Acts::BinFinder<SpacePoint>());
199  auto topBinFinder = std::make_shared<Acts::BinFinder<SpacePoint>>(
200  Acts::BinFinder<SpacePoint>());
201  Acts::SeedFilterConfig sfconf;
202  Acts::ATLASCuts<SpacePoint> atlasCuts = Acts::ATLASCuts<SpacePoint>();
203  config.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
204  Acts::SeedFilter<SpacePoint>(sfconf, &atlasCuts));
205  Acts::Seedfinder<SpacePoint> seedfinder_cpu(config);
206  Acts::Seedfinder<SpacePoint, Acts::Cuda> seedfinder_cuda(config);
207 
208  // covariance tool, sets covariances per spacepoint as required
209  auto ct = [=](const SpacePoint& sp, float, float, float) -> Acts::Vector2D {
210  return {sp.varianceR, sp.varianceZ};
211  };
212 
213  // setup spacepoint grid config
214  Acts::SpacePointGridConfig gridConf;
215  gridConf.bFieldInZ = config.bFieldInZ;
216  gridConf.minPt = config.minPt;
217  gridConf.rMax = config.rMax;
218  gridConf.zMax = config.zMax;
219  gridConf.zMin = config.zMin;
220  gridConf.deltaRMax = config.deltaRMax;
221  gridConf.cotThetaMax = config.cotThetaMax;
222  // create grid with bin sizes according to the configured geometry
223  std::unique_ptr<Acts::SpacePointGrid<SpacePoint>> grid =
224  Acts::SpacePointGridCreator::createGrid<SpacePoint>(gridConf);
225  auto spGroup = Acts::BinnedSPGroup<SpacePoint>(spVec.begin(), spVec.end(), ct,
226  bottomBinFinder, topBinFinder,
227  std::move(grid), config);
228 
229  auto end_pre = std::chrono::system_clock::now();
230  std::chrono::duration<double> elapsec_pre = end_pre - start_pre;
231  double preprocessTime = elapsec_pre.count();
232  std::cout << "Preprocess Time: " << preprocessTime << std::endl;
233 
234  //--------------------------------------------------------------------//
235  // Begin Seed finding //
236  //--------------------------------------------------------------------//
237 
238  auto start_cpu = std::chrono::system_clock::now();
239 
240  int group_count;
241  auto groupIt = spGroup.begin();
242 
243  //----------- CPU ----------//
244  group_count = 0;
245  std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_cpu;
246  groupIt = spGroup.begin();
247 
248  if (do_cpu) {
249  for (int i_s = 0; i_s < skip; i_s++)
250  ++groupIt;
251  for (; !(groupIt == spGroup.end()); ++groupIt) {
252  seedVector_cpu.push_back(seedfinder_cpu.createSeedsForGroup(
253  groupIt.bottom(), groupIt.middle(), groupIt.top()));
254  group_count++;
255  if (allgroup == false) {
256  if (group_count >= nGroupToIterate)
257  break;
258  }
259  }
260  // auto timeMetric_cpu = seedfinder_cpu.getTimeMetric();
261  std::cout << "Analyzed " << group_count << " groups for CPU" << std::endl;
262  }
263 
264  auto end_cpu = std::chrono::system_clock::now();
265  std::chrono::duration<double> elapsec_cpu = end_cpu - start_cpu;
266  double cpuTime = elapsec_cpu.count();
267 
268  //----------- CUDA ----------//
269 
270  cudaProfilerStart();
271  auto start_cuda = std::chrono::system_clock::now();
272 
273  group_count = 0;
274  std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_cuda;
275  groupIt = spGroup.begin();
276 
277  for (int i_s = 0; i_s < skip; i_s++)
278  ++groupIt;
279  for (; !(groupIt == spGroup.end()); ++groupIt) {
280  seedVector_cuda.push_back(seedfinder_cuda.createSeedsForGroup(
281  groupIt.bottom(), groupIt.middle(), groupIt.top()));
282  group_count++;
283  if (allgroup == false) {
284  if (group_count >= nGroupToIterate)
285  break;
286  }
287  }
288  auto end_cuda = std::chrono::system_clock::now();
289  std::chrono::duration<double> elapsec_cuda = end_cuda - start_cuda;
290  double cudaTime = elapsec_cuda.count();
291 
292  cudaProfilerStop();
293  std::cout << "Analyzed " << group_count << " groups for CUDA" << std::endl;
294 
295  std::cout << std::endl;
296  std::cout << "----------------------- Time Metric -----------------------"
297  << std::endl;
298  std::cout << " " << (do_cpu ? "CPU" : " ")
299  << " CUDA " << (do_cpu ? "Speedup " : "")
300  << std::endl;
301  std::cout << "Seedfinding_Time " << std::setw(11)
302  << (do_cpu ? std::to_string(cpuTime) : "") << " " << std::setw(11)
303  << cudaTime << " " << std::setw(11)
304  << (do_cpu ? std::to_string(cpuTime / cudaTime) : "") << std::endl;
305  double wallTime_cpu = cpuTime + preprocessTime;
306  double wallTime_cuda = cudaTime + preprocessTime;
307  std::cout << "Wall_time " << std::setw(11)
308  << (do_cpu ? std::to_string(wallTime_cpu) : "") << " "
309  << std::setw(11) << wallTime_cuda << " " << std::setw(11)
310  << (do_cpu ? std::to_string(wallTime_cpu / wallTime_cuda) : "")
311  << std::endl;
312  std::cout << "-----------------------------------------------------------"
313  << std::endl;
314  std::cout << std::endl;
315 
316  int nSeed_cpu = 0;
317  for (auto& outVec : seedVector_cpu) {
318  nSeed_cpu += outVec.size();
319  }
320 
321  int nSeed_cuda = 0;
322  for (auto& outVec : seedVector_cuda) {
323  nSeed_cuda += outVec.size();
324  }
325 
326  std::cout << "Number of Seeds (CPU | CUDA): " << nSeed_cpu << " | "
327  << nSeed_cuda << std::endl;
328 
329  int nMatch = 0;
330 
331  for (size_t i = 0; i < seedVector_cpu.size(); i++) {
332  auto regionVec_cpu = seedVector_cpu[i];
333  auto regionVec_cuda = seedVector_cuda[i];
334 
335  std::vector<std::vector<SpacePoint>> seeds_cpu;
336  std::vector<std::vector<SpacePoint>> seeds_cuda;
337 
338  // for (size_t i_cpu = 0; i_cpu < regionVec_cpu.size(); i_cpu++) {
339  for (auto sd : regionVec_cpu) {
340  std::vector<SpacePoint> seed_cpu;
341  seed_cpu.push_back(*(sd.sp()[0]));
342  seed_cpu.push_back(*(sd.sp()[1]));
343  seed_cpu.push_back(*(sd.sp()[2]));
344 
345  seeds_cpu.push_back(seed_cpu);
346  }
347 
348  for (auto sd : regionVec_cuda) {
349  std::vector<SpacePoint> seed_cuda;
350  seed_cuda.push_back(*(sd.sp()[0]));
351  seed_cuda.push_back(*(sd.sp()[1]));
352  seed_cuda.push_back(*(sd.sp()[2]));
353 
354  seeds_cuda.push_back(seed_cuda);
355  }
356 
357  for (auto seed : seeds_cpu) {
358  for (auto other : seeds_cuda) {
359  if (seed[0] == other[0] && seed[1] == other[1] && seed[2] == other[2]) {
360  nMatch++;
361  break;
362  }
363  }
364  }
365  }
366 
367  if (do_cpu) {
368  std::cout << nMatch << " seeds are matched" << std::endl;
369  std::cout << "Matching rate: " << float(nMatch) / nSeed_cpu * 100 << "%"
370  << std::endl;
371  }
372 
373  if (!quiet) {
374  if (do_cpu) {
375  std::cout << "CPU Seed result:" << std::endl;
376 
377  for (auto& regionVec : seedVector_cpu) {
378  for (size_t i = 0; i < regionVec.size(); i++) {
379  const Acts::Seed<SpacePoint>* seed = &regionVec[i];
380  const SpacePoint* sp = seed->sp()[0];
381  std::cout << " (" << sp->x() << ", " << sp->y() << ", " << sp->z()
382  << ") ";
383  sp = seed->sp()[1];
384  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
385  << sp->z() << ") ";
386  sp = seed->sp()[2];
387  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
388  << sp->z() << ") ";
389  std::cout << std::endl;
390  }
391  }
392 
393  std::cout << std::endl;
394  }
395  std::cout << "CUDA Seed result:" << std::endl;
396 
397  for (auto& regionVec : seedVector_cuda) {
398  for (size_t i = 0; i < regionVec.size(); i++) {
399  const Acts::Seed<SpacePoint>* seed = &regionVec[i];
400  const SpacePoint* sp = seed->sp()[0];
401  std::cout << " (" << sp->x() << ", " << sp->y() << ", " << sp->z()
402  << ") ";
403  sp = seed->sp()[1];
404  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
405  << sp->z() << ") ";
406  sp = seed->sp()[2];
407  std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
408  << sp->z() << ") ";
409  std::cout << std::endl;
410  }
411  }
412  }
413 
414  std::cout << std::endl;
415  std::cout << std::endl;
416 
417  return 0;
418 }