Cuda/Seeding/SeedFinderCudaTest.cpp

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 "Acts/Plugins/Cuda/Seeding/SeedFinder.hpp"
0010 #include "Acts/Seeding/BinnedGroup.hpp"
0011 #include "Acts/Seeding/InternalSeed.hpp"
0012 #include "Acts/Seeding/InternalSpacePoint.hpp"
0013 #include "Acts/Seeding/Seed.hpp"
0014 #include "Acts/Seeding/SeedFilter.hpp"
0015 #include "Acts/Seeding/SeedFinder.hpp"
0016 #include "Acts/Seeding/SpacePointGrid.hpp"
0017 #include "Acts/Utilities/GridBinFinder.hpp"
0018
0019 #include <chrono>
0020 #include <fstream>
0021 #include <iomanip>
0022 #include <iostream>
0023 #include <optional>
0024 #include <sstream>
0025 #include <utility>
0026
0027 #include <boost/type_erasure/any_cast.hpp>
0028 #include <cuda_profiler_api.h>
0029
0030 #include "ATLASCuts.hpp"
0031 #include "SpacePoint.hpp"
0032
0033 using namespace Acts::UnitLiterals;
0034
0035 std::vector<const SpacePoint*> readFile(std::string filename) {
0036   std::string line;
0037   int layer;
0038   std::vector<const SpacePoint*> readSP;
0039
0040   std::ifstream spFile(filename);
0041   if (spFile.is_open()) {
0042     while (!spFile.eof()) {
0043       std::getline(spFile, line);
0044       std::stringstream ss(line);
0045       std::string linetype;
0046       ss >> linetype;
0047       float x, y, z, r, varianceR, varianceZ;
0048       if (linetype == "lxyz") {
0049         ss >> layer >> x >> y >> z >> varianceR >> varianceZ;
0050         r = std::hypot(x, y);
0051         float f22 = varianceR;
0052         float wid = varianceZ;
0053         float cov = wid * wid * .08333;
0054         if (cov < f22)
0055           cov = f22;
0056         if (std::abs(z) > 450.) {
0057           varianceZ = 9. * cov;
0058           varianceR = .06;
0059         } else {
0060           varianceR = 9. * cov;
0061           varianceZ = .06;
0062         }
0063         SpacePoint* sp =
0064             new SpacePoint{x, y, z, r, layer, varianceR, varianceZ};
0065         //     if(r < 200.){
0066         //       sp->setClusterList(1,0);
0067         //     }
0068         readSP.push_back(sp);
0069       }
0070     }
0071   }
0072   return readSP;
0073 }
0074
0075 int main(int argc, char** argv) {
0076   auto start_pre = std::chrono::system_clock::now();
0077
0078   std::string file{"sp.txt"};
0079   bool help(false);
0080   bool quiet(false);
0081   bool allgroup(false);
0082   bool do_cpu(true);
0083   int nGroupToIterate = 500;
0084   int skip = 0;
0085   int deviceID = 0;
0086   int nTrplPerSpBLimit = 100;
0087   int nAvgTrplPerSpBLimit = 2;
0088
0089   int opt;
0090   while ((opt = getopt(argc, argv, "haf:n:s:d:l:m:qG")) != -1) {
0091     switch (opt) {
0092       case 'a':
0093         allgroup = true;
0094         break;
0095       case 'f':
0096         file = optarg;
0097         break;
0098       case 'n':
0099         nGroupToIterate = atoi(optarg);
0100         break;
0101       case 's':
0102         skip = atoi(optarg);
0103         break;
0104       case 'd':
0105         deviceID = atoi(optarg);
0106         break;
0107       case 'l':
0108         nAvgTrplPerSpBLimit = atoi(optarg);
0109         break;
0110       case 'm':
0111         nTrplPerSpBLimit = atoi(optarg);
0112         break;
0113       case 'q':
0114         quiet = true;
0115         break;
0116       case 'G':
0117         do_cpu = false;
0118         break;
0119       case 'h':
0120         help = true;
0121         [[fallthrough]];
0122       default: /* '?' */
0123         std::cerr << "Usage: " << argv[0] << " [-hq] [-f FILENAME]\n";
0124         if (help) {
0125           std::cout << "      -h : this help" << std::endl;
0126           std::cout << "      -a ALL   : analyze all groups. Default is \""
0127                     << allgroup << "\"" << std::endl;
0128           std::cout
0129               << "      -f FILE  : read spacepoints from FILE. Default is \""
0130               << file << "\"" << std::endl;
0131           std::cout << "      -n NUM   : Number of groups to iterate in seed "
0132                        "finding. Default is "
0133                     << nGroupToIterate << std::endl;
0134           std::cout << "      -s SKIP  : Number of groups to skip in seed "
0135                        "finding. Default is "
0136                     << skip << std::endl;
0137           std::cout << "      -d DEVID : NVIDIA GPU device ID. Default is "
0138                     << deviceID << std::endl;
0139           std::cout << "      -l : A limit on the average number of triplets "
0140                        "per bottom spacepoint: this is used for determining "
0141                        "matrix size for triplets per middle space point"
0142                     << nAvgTrplPerSpBLimit << std::endl;
0143           std::cout << "      -m : A limit on the number of triplets per "
0144                        "bottom spacepoint: users do not have to touch this for "
0145                        "# spacepoints < ~200k"
0146                     << nTrplPerSpBLimit << std::endl;
0147           std::cout << "      -q : don't print out all found seeds"
0148                     << std::endl;
0149           std::cout << "      -G : only run on GPU, not CPU" << std::endl;
0150         }
0151
0152         exit(EXIT_FAILURE);
0153     }
0154   }
0155
0156   std::string devName;
0157   ACTS_CUDA_ERROR_CHECK(cudaSetDevice(deviceID));
0158
0159   std::ifstream f(file);
0160   if (!f.good()) {
0161     std::cerr << "input file \"" << file << "\" does not exist\n";
0162     exit(EXIT_FAILURE);
0163   }
0164
0165   std::vector<const SpacePoint*> spVec = readFile(file);
0166
0167   std::cout << "read " << spVec.size() << " SP from file " << file << std::endl;
0168
0169   // Set seed finder configuration
0170   Acts::SeedFinderConfig<SpacePoint> config;
0171   // silicon detector max
0172   config.rMax = 160._mm;
0173   config.deltaRMin = 5._mm;
0174   config.deltaRMax = 160._mm;
0175   config.deltaRMinTopSP = config.deltaRMin;
0176   config.deltaRMinBottomSP = config.deltaRMin;
0177   config.deltaRMaxTopSP = config.deltaRMax;
0178   config.deltaRMaxBottomSP = config.deltaRMax;
0179   config.collisionRegionMin = -250._mm;
0180   config.collisionRegionMax = 250._mm;
0181   config.zMin = -2800._mm;
0182   config.zMax = 2800._mm;
0183   config.maxSeedsPerSpM = 5;
0184   // 2.7 eta
0185   config.cotThetaMax = 7.40627;
0186   config.sigmaScattering = 1.00000;
0187
0188   config.minPt = 500._MeV;
0189   config.impactMax = 10._mm;
0190   config = config.toInternalUnits();
0191
0192   Acts::SeedFinderOptions options;
0193   options.bFieldInZ = 2_T;
0194   options.beamPos = {-.5_mm, -.5_mm};
0195   options = options.toInternalUnits().calculateDerivedQuantities(config);
0196
0197   int numPhiNeighbors = 1;
0198
0199   // extent used to store r range for middle spacepoint
0200   Acts::Extent rRangeSPExtent;
0201
0202   const Acts::Range1D<float> rMiddleSPRange;
0203
0204   std::vector<std::pair<int, int>> zBinNeighborsTop;
0205   std::vector<std::pair<int, int>> zBinNeighborsBottom;
0206
0207   // cuda
0208   cudaDeviceProp prop;
0209   ACTS_CUDA_ERROR_CHECK(cudaGetDeviceProperties(&prop, deviceID));
0210   printf("\n GPU Device %d: \"%s\" with compute capability %d.%d\n\n", deviceID,
0211          prop.name, prop.major, prop.minor);
0212   config.maxBlockSize = prop.maxThreadsPerBlock;
0213   config.nTrplPerSpBLimit = nTrplPerSpBLimit;
0214   config.nAvgTrplPerSpBLimit = nAvgTrplPerSpBLimit;
0215
0216   // binfinder
0217   auto bottomBinFinder = std::make_unique<Acts::GridBinFinder<2ul>>(
0218       Acts::GridBinFinder<2ul>(numPhiNeighbors, zBinNeighborsBottom));
0219   auto topBinFinder = std::make_unique<Acts::GridBinFinder<2ul>>(
0220       Acts::GridBinFinder<2ul>(numPhiNeighbors, zBinNeighborsTop));
0221   Acts::SeedFilterConfig sfconf;
0222   Acts::ATLASCuts<SpacePoint> atlasCuts = Acts::ATLASCuts<SpacePoint>();
0223   config.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
0224       Acts::SeedFilter<SpacePoint>(sfconf, &atlasCuts));
0225   Acts::SeedFinder<SpacePoint, Acts::CylindricalSpacePointGrid<SpacePoint>>
0226       seedFinder_cpu(config);
0227   Acts::SeedFinder<SpacePoint, Acts::Cuda> seedFinder_cuda(config, options);
0228
0229   // covariance tool, sets covariances per spacepoint as required
0230   auto ct = [=](const SpacePoint& sp, float, float, float)
0231       -> std::tuple<Acts::Vector3, Acts::Vector2, std::optional<float>> {
0232     Acts::Vector3 position(sp.x(), sp.y(), sp.z());
0233     Acts::Vector2 variance(sp.varianceR, sp.varianceZ);
0234     return {position, variance, std::nullopt};
0235   };
0236
0237   // setup spacepoint grid config
0238   Acts::CylindricalSpacePointGridConfig gridConf;
0239   gridConf.minPt = config.minPt;
0240   gridConf.rMax = config.rMax;
0241   gridConf.zMax = config.zMax;
0242   gridConf.zMin = config.zMin;
0243   gridConf.deltaRMax = config.deltaRMax;
0244   gridConf.cotThetaMax = config.cotThetaMax;
0245   // setup spacepoint grid options
0246   Acts::CylindricalSpacePointGridOptions gridOpts;
0247   gridOpts.bFieldInZ = options.bFieldInZ;
0248   // create grid with bin sizes according to the configured geometry
0249   Acts::CylindricalSpacePointGrid<SpacePoint> grid =
0250       Acts::CylindricalSpacePointGridCreator::createGrid<SpacePoint>(gridConf,
0251                                                                      gridOpts);
0252   Acts::CylindricalSpacePointGridCreator::fillGrid(
0253       config, options, grid, spVec.begin(), spVec.end(), ct, rRangeSPExtent);
0254
0255   auto spGroup = Acts::CylindricalBinnedGroup<SpacePoint>(
0256       std::move(grid), *bottomBinFinder, *topBinFinder);
0257
0258   auto end_pre = std::chrono::system_clock::now();
0259   std::chrono::duration<double> elapsec_pre = end_pre - start_pre;
0260   double preprocessTime = elapsec_pre.count();
0261   std::cout << "Preprocess Time: " << preprocessTime << std::endl;
0262
0263   //--------------------------------------------------------------------//
0264   //                        Begin Seed finding                          //
0265   //--------------------------------------------------------------------//
0266
0267   auto start_cpu = std::chrono::system_clock::now();
0268
0269   std::array<std::vector<std::size_t>, 2ul> navigation;
0270   navigation[0ul].resize(spGroup.grid().numLocalBins()[0ul]);
0271   navigation[1ul].resize(spGroup.grid().numLocalBins()[1ul]);
0272   std::iota(navigation[0ul].begin(), navigation[0ul].end(), 1ul);
0273   std::iota(navigation[1ul].begin(), navigation[1ul].end(), 1ul);
0274
0275   std::array<std::size_t, 2ul> localPosition =
0276       spGroup.grid().localBinsFromGlobalBin(skip);
0277
0278   int group_count;
0279   auto groupIt = Acts::CylindricalBinnedGroupIterator<SpacePoint>(
0280       spGroup, localPosition, navigation);
0281
0282   //----------- CPU ----------//
0283   group_count = 0;
0284   std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_cpu;
0285
0286   if (do_cpu) {
0287     decltype(seedFinder_cpu)::SeedingState state;
0288     state.spacePointData.resize(spVec.size());
0289     for (; groupIt != spGroup.end(); ++groupIt) {
0290       const auto [bottom, middle, top] = *groupIt;
0291       seedFinder_cpu.createSeedsForGroup(
0292           options, state, spGroup.grid(),
0293           std::back_inserter(seedVector_cpu.emplace_back()), bottom, middle,
0294           top, rMiddleSPRange);
0295       group_count++;
0296       if (allgroup == false) {
0297         if (group_count >= nGroupToIterate)
0298           break;
0299       }
0300     }
0301     // auto timeMetric_cpu = seedFinder_cpu.getTimeMetric();
0302     std::cout << "Analyzed " << group_count << " groups for CPU" << std::endl;
0303   }
0304
0305   auto end_cpu = std::chrono::system_clock::now();
0306   std::chrono::duration<double> elapsec_cpu = end_cpu - start_cpu;
0307   double cpuTime = elapsec_cpu.count();
0308
0309   //----------- CUDA ----------//
0310
0311   cudaProfilerStart();
0312   auto start_cuda = std::chrono::system_clock::now();
0313
0314   group_count = 0;
0315   std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_cuda;
0316   groupIt = Acts::CylindricalBinnedGroupIterator<SpacePoint>(
0317       spGroup, localPosition, navigation);
0318
0319   Acts::SpacePointData spacePointData;
0320   spacePointData.resize(spVec.size());
0321
0322   for (; groupIt != spGroup.end(); ++groupIt) {
0323     const auto [bottom, middle, top] = *groupIt;
0324     seedVector_cuda.push_back(seedFinder_cuda.createSeedsForGroup(
0325         spacePointData, spGroup.grid(), bottom, middle, top));
0326     group_count++;
0327     if (allgroup == false) {
0328       if (group_count >= nGroupToIterate)
0329         break;
0330     }
0331   }
0332
0333   auto end_cuda = std::chrono::system_clock::now();
0334   std::chrono::duration<double> elapsec_cuda = end_cuda - start_cuda;
0335   double cudaTime = elapsec_cuda.count();
0336
0337   cudaProfilerStop();
0338   std::cout << "Analyzed " << group_count << " groups for CUDA" << std::endl;
0339
0340   std::cout << std::endl;
0341   std::cout << "----------------------- Time Metric -----------------------"
0342             << std::endl;
0343   std::cout << "                       " << (do_cpu ? "CPU" : "   ")
0344             << "          CUDA        " << (do_cpu ? "Speedup " : "")
0345             << std::endl;
0346   std::cout << "Seedfinding_Time  " << std::setw(11)
0347             << (do_cpu ? std::to_string(cpuTime) : "") << "  " << std::setw(11)
0348             << cudaTime << "  " << std::setw(11)
0349             << (do_cpu ? std::to_string(cpuTime / cudaTime) : "") << std::endl;
0350   double wallTime_cpu = cpuTime + preprocessTime;
0351   double wallTime_cuda = cudaTime + preprocessTime;
0352   std::cout << "Wall_time         " << std::setw(11)
0353             << (do_cpu ? std::to_string(wallTime_cpu) : "") << "  "
0354             << std::setw(11) << wallTime_cuda << "  " << std::setw(11)
0355             << (do_cpu ? std::to_string(wallTime_cpu / wallTime_cuda) : "")
0356             << std::endl;
0357   std::cout << "-----------------------------------------------------------"
0358             << std::endl;
0359   std::cout << std::endl;
0360
0361   int nSeed_cpu = 0;
0362   for (auto& outVec : seedVector_cpu) {
0363     nSeed_cpu += outVec.size();
0364   }
0365
0366   int nSeed_cuda = 0;
0367   for (auto& outVec : seedVector_cuda) {
0368     nSeed_cuda += outVec.size();
0369   }
0370
0371   std::cout << "Number of Seeds (CPU | CUDA): " << nSeed_cpu << " | "
0372             << nSeed_cuda << std::endl;
0373
0374   int nMatch = 0;
0375
0376   for (std::size_t i = 0; i < seedVector_cpu.size(); i++) {
0377     auto regionVec_cpu = seedVector_cpu[i];
0378     auto regionVec_cuda = seedVector_cuda[i];
0379
0380     std::vector<std::vector<SpacePoint>> seeds_cpu;
0381     std::vector<std::vector<SpacePoint>> seeds_cuda;
0382
0383     // for (std::size_t i_cpu = 0; i_cpu < regionVec_cpu.size(); i_cpu++) {
0384     for (auto sd : regionVec_cpu) {
0385       std::vector<SpacePoint> seed_cpu;
0386       seed_cpu.push_back(*(sd.sp()[0]));
0387       seed_cpu.push_back(*(sd.sp()[1]));
0388       seed_cpu.push_back(*(sd.sp()[2]));
0389
0390       seeds_cpu.push_back(seed_cpu);
0391     }
0392
0393     for (auto sd : regionVec_cuda) {
0394       std::vector<SpacePoint> seed_cuda;
0395       seed_cuda.push_back(*(sd.sp()[0]));
0396       seed_cuda.push_back(*(sd.sp()[1]));
0397       seed_cuda.push_back(*(sd.sp()[2]));
0398
0399       seeds_cuda.push_back(seed_cuda);
0400     }
0401
0402     for (auto seed : seeds_cpu) {
0403       for (auto other : seeds_cuda) {
0404         if (seed[0] == other[0] && seed[1] == other[1] && seed[2] == other[2]) {
0405           nMatch++;
0406           break;
0407         }
0408       }
0409     }
0410   }
0411
0412   if (do_cpu) {
0413     std::cout << nMatch << " seeds are matched" << std::endl;
0414     std::cout << "Matching rate: " << float(nMatch) / nSeed_cpu * 100 << "%"
0415               << std::endl;
0416   }
0417
0418   if (!quiet) {
0419     if (do_cpu) {
0420       std::cout << "CPU Seed result:" << std::endl;
0421
0422       for (auto& regionVec : seedVector_cpu) {
0423         for (std::size_t i = 0; i < regionVec.size(); i++) {
0424           const Acts::Seed<SpacePoint>* seed = &regionVec[i];
0425           const SpacePoint* sp = seed->sp()[0];
0426           std::cout << " (" << sp->x() << ", " << sp->y() << ", " << sp->z()
0427                     << ") ";
0428           sp = seed->sp()[1];
0429           std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
0430                     << sp->z() << ") ";
0431           sp = seed->sp()[2];
0432           std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
0433                     << sp->z() << ") ";
0434           std::cout << std::endl;
0435         }
0436       }
0437
0438       std::cout << std::endl;
0439     }
0440     std::cout << "CUDA Seed result:" << std::endl;
0441
0442     for (auto& regionVec : seedVector_cuda) {
0443       for (std::size_t i = 0; i < regionVec.size(); i++) {
0444         const Acts::Seed<SpacePoint>* seed = &regionVec[i];
0445         const SpacePoint* sp = seed->sp()[0];
0446         std::cout << " (" << sp->x() << ", " << sp->y() << ", " << sp->z()
0447                   << ") ";
0448         sp = seed->sp()[1];
0449         std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
0450                   << sp->z() << ") ";
0451         sp = seed->sp()[2];
0452         std::cout << sp->surface << " (" << sp->x() << ", " << sp->y() << ", "
0453                   << sp->z() << ") ";
0454         std::cout << std::endl;
0455       }
0456     }
0457   }
0458
0459   std::cout << std::endl;
0460   std::cout << std::endl;
0461
0462   return 0;
0463 }