Sycl/Seeding/SeedFinderSyclTest.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/EventData/SpacePointData.hpp"
0010 #include "Acts/Plugins/Sycl/Seeding/SeedFinder.hpp"
0011 #include "Acts/Plugins/Sycl/Utilities/QueueWrapper.hpp"
0012 #include "Acts/Seeding/BinnedGroup.hpp"
0013 #include "Acts/Seeding/InternalSeed.hpp"
0014 #include "Acts/Seeding/InternalSpacePoint.hpp"
0015 #include "Acts/Seeding/Seed.hpp"
0016 #include "Acts/Seeding/SeedFilter.hpp"
0017 #include "Acts/Seeding/SeedFinder.hpp"
0018 #include "Acts/Seeding/SeedFinderConfig.hpp"
0019 #include "Acts/Seeding/SpacePointGrid.hpp"
0020 #include "Acts/Utilities/GridBinFinder.hpp"
0021 #include "Acts/Utilities/Logger.hpp"
0022
0023 #include <chrono>
0024 #include <cmath>
0025 #include <ctime>
0026 #include <fstream>
0027 #include <iomanip>
0028 #include <iostream>
0029 #include <limits>
0030 #include <memory>
0031 #include <optional>
0032 #include <sstream>
0033 #include <string>
0034
0035 #include <boost/type_erasure/any_cast.hpp>
0036
0037 #include "ATLASCuts.hpp"
0038 #include "CommandLineArguments.h"
0039 #include "SpacePoint.hpp"
0040 #include "vecmem/memory/sycl/device_memory_resource.hpp"
0041 #include "vecmem/memory/sycl/host_memory_resource.hpp"
0042
0043 using namespace Acts::UnitLiterals;
0044
0045 auto readFile(const std::string& filename) -> std::vector<const SpacePoint*> {
0046   std::string line;
0047   std::vector<const SpacePoint*> readSP;
0048
0049   std::ifstream spFile(filename);
0050   if (spFile.is_open()) {
0051     int id = 0;
0052     while (std::getline(spFile, line)) {
0053       std::stringstream ss(line);
0054       std::string linetype;
0055       ss >> linetype;
0056       if (linetype == "lxyz") {
0057         float x;
0058         float y;
0059         float z;
0060         float r;
0061         float varianceR;
0062         float varianceZ;
0063         int layer;
0064         ss >> layer >> x >> y >> z >> varianceR >> varianceZ;
0065         r = std::hypot(x, y);
0066         float f22 = varianceR;
0067         float wid = varianceZ;
0068         float cov = wid * wid * .08333;
0069         if (cov < f22)
0070           cov = f22;
0071         if (std::abs(z) > 450.) {
0072           varianceZ = 9. * cov;
0073           varianceR = .06;
0074         } else {
0075           varianceR = 9. * cov;
0076           varianceZ = .06;
0077         }
0078         readSP.emplace_back(
0079             new SpacePoint(x, y, z, r, layer, varianceR, varianceZ, id));
0080         ++id;
0081       }
0082     }
0083   }
0084   return readSP;
0085 }
0086
0087 template <typename external_spacepoint_t>
0088 auto setupSeedFinderConfiguration()
0089     -> Acts::SeedFinderConfig<external_spacepoint_t> {
0090   Acts::SeedFinderConfig<SpacePoint> config;
0091   // silicon detector max
0092   config.rMax = 160._mm;
0093   config.deltaRMin = 5._mm;
0094   config.deltaRMax = 160._mm;
0095   config.deltaRMinTopSP = config.deltaRMin;
0096   config.deltaRMinBottomSP = config.deltaRMin;
0097   config.deltaRMaxTopSP = config.deltaRMax;
0098   config.deltaRMaxBottomSP = config.deltaRMax;
0099   config.collisionRegionMin = -250._mm;
0100   config.collisionRegionMax = 250._mm;
0101   config.zMin = -2800._mm;
0102   config.zMax = 2800._mm;
0103   config.maxSeedsPerSpM = 5;
0104   // 2.7 eta
0105   config.cotThetaMax = 7.40627;
0106   config.sigmaScattering = 1.00000;
0107   config.minPt = 500._MeV;
0108   config.impactMax = 10._mm;
0109   // for sycl
0110   config.nTrplPerSpBLimit = 100;
0111   config.nAvgTrplPerSpBLimit = 6;
0112   return config;
0113 }
0114
0115 auto setupSeedFinderOptions() {
0116   Acts::SeedFinderOptions options;
0117   options.bFieldInZ = 2_T;
0118   options.beamPos = {-.5_mm, -.5_mm};
0119   return options;
0120 }
0121
0122 template <typename external_spacepoint_t>
0123 auto setupSpacePointGridConfig(
0124     const Acts::SeedFinderConfig<external_spacepoint_t>& config,
0125     const Acts::SeedFinderOptions& options)
0126     -> std::pair<Acts::CylindricalSpacePointGridConfig,
0127                  Acts::CylindricalSpacePointGridOptions> {
0128   Acts::CylindricalSpacePointGridConfig gridConf{};
0129   gridConf.minPt = config.minPt;
0130   gridConf.rMax = config.rMax;
0131   gridConf.zMax = config.zMax;
0132   gridConf.zMin = config.zMin;
0133   gridConf.deltaRMax = config.deltaRMax;
0134   gridConf.cotThetaMax = config.cotThetaMax;
0135
0136   Acts::CylindricalSpacePointGridOptions gridOpts{};
0137   gridOpts.bFieldInZ = options.bFieldInZ;
0138   return std::make_pair(gridConf, gridOpts);
0139 }
0140
0141 auto main(int argc, char** argv) -> int {
0142   auto start_prep = std::chrono::system_clock::now();
0143
0144   CommandLineArguments cmdlTool;
0145   cmdlTool.parse(argc, argv);
0146
0147   if (!cmdlTool.inpFileName.empty() && !cmdlTool.inpFileExists) {
0148     std::cerr << "Input file not found\n";
0149     return -1;
0150   }
0151
0152   auto spVec = readFile(cmdlTool.inpFileName);
0153
0154   int numPhiNeighbors = 1;
0155
0156   // extent used to store r range for middle spacepoint
0157   Acts::Extent rRangeSPExtent;
0158
0159   const Acts::Range1D<float> rMiddleSPRange;
0160
0161   std::vector<std::pair<int, int>> zBinNeighborsTop;
0162   std::vector<std::pair<int, int>> zBinNeighborsBottom;
0163
0164   auto bottomBinFinder = std::make_unique<Acts::GridBinFinder<2ul>>(
0165       Acts::GridBinFinder<2ul>(numPhiNeighbors, zBinNeighborsBottom));
0166   auto topBinFinder = std::make_unique<Acts::GridBinFinder<2ul>>(
0167       Acts::GridBinFinder<2ul>(numPhiNeighbors, zBinNeighborsTop));
0168   auto config = setupSeedFinderConfiguration<SpacePoint>();
0169   config = config.toInternalUnits().calculateDerivedQuantities();
0170   auto options = setupSeedFinderOptions();
0171   options = options.toInternalUnits().calculateDerivedQuantities(config);
0172
0173   Acts::ATLASCuts<SpacePoint> atlasCuts = Acts::ATLASCuts<SpacePoint>();
0174   Acts::Sycl::DeviceExperimentCuts deviceAtlasCuts;
0175   config.seedFilter = std::make_unique<Acts::SeedFilter<SpacePoint>>(
0176       Acts::SeedFilter<SpacePoint>(Acts::SeedFilterConfig(), &atlasCuts));
0177
0178   const Acts::Logging::Level logLvl =
0179       cmdlTool.csvFormat ? Acts::Logging::WARNING : Acts::Logging::INFO;
0180   Acts::Sycl::QueueWrapper queue(
0181       cmdlTool.deviceName,
0182       Acts::getDefaultLogger("Sycl::QueueWrapper", logLvl));
0183   vecmem::sycl::host_memory_resource resource(queue.getQueue());
0184   vecmem::sycl::device_memory_resource device_resource(queue.getQueue());
0185   Acts::Sycl::SeedFinder<SpacePoint> syclSeedFinder(
0186       config, options, deviceAtlasCuts, queue, resource, &device_resource);
0187   Acts::SeedFinder<SpacePoint, Acts::CylindricalSpacePointGrid<SpacePoint>>
0188       normalSeedFinder(config);
0189   auto globalTool = [=](const SpacePoint& sp, float /*unused*/,
0190                         float /*unused*/, float /*unused*/)
0191       -> std::tuple<Acts::Vector3, Acts::Vector2, std::optional<float>> {
0192     Acts::Vector3 position(sp.x(), sp.y(), sp.z());
0193     Acts::Vector2 covariance(sp.varianceR, sp.varianceZ);
0194     return std::make_tuple(position, covariance, std::nullopt);
0195   };
0196   auto [gridConfig, gridOpts] = setupSpacePointGridConfig(config, options);
0197   gridConfig = gridConfig.toInternalUnits();
0198   gridOpts = gridOpts.toInternalUnits();
0199   Acts::CylindricalSpacePointGrid<SpacePoint> grid =
0200       Acts::CylindricalSpacePointGridCreator::createGrid<SpacePoint>(gridConfig,
0201                                                                      gridOpts);
0202   Acts::CylindricalSpacePointGridCreator::fillGrid(config, options, grid,
0203                                                    spVec.begin(), spVec.end(),
0204                                                    globalTool, rRangeSPExtent);
0205
0206   std::array<std::vector<std::size_t>, 2ul> navigation;
0207   auto spGroup = Acts::CylindricalBinnedGroup<SpacePoint>(
0208       std::move(grid), *bottomBinFinder, *topBinFinder, std::move(navigation));
0209
0210   auto end_prep = std::chrono::system_clock::now();
0211
0212   std::chrono::duration<double> elapsec_prep = end_prep - start_prep;
0213   double prepTime = elapsec_prep.count();
0214
0215   if (!cmdlTool.csvFormat) {
0216     std::cout << "read " << spVec.size() << " SP from file "
0217               << cmdlTool.inpFileName << std::endl;
0218     std::cout << "Preparation time: " << std::to_string(prepTime) << std::endl;
0219   }
0220
0221   // -------------------------------------- //
0222   // ----------- EXECUTE ON CPU ----------- //
0223   // -------------------------------------- //
0224
0225   auto start_cpu = std::chrono::system_clock::now();
0226   uint group_count = 0;
0227   std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_cpu;
0228
0229   if (!cmdlTool.onlyGpu) {
0230     decltype(normalSeedFinder)::SeedingState state;
0231     for (auto [bottom, middle, top] : spGroup) {
0232       normalSeedFinder.createSeedsForGroup(
0233           options, state, spGroup.grid(),
0234           std::back_inserter(seedVector_cpu.emplace_back()), bottom, middle,
0235           top, rMiddleSPRange);
0236       group_count++;
0237       if (!cmdlTool.allgroup && group_count >= cmdlTool.groups) {
0238         break;
0239       }
0240     }
0241   }
0242
0243   auto end_cpu = std::chrono::system_clock::now();
0244
0245   if (!cmdlTool.csvFormat) {
0246     std::cout << "Analyzed " << group_count << " groups for CPU" << std::endl;
0247   }
0248
0249   // -------------------------------------- //
0250   // -------- EXECUTE ON GPU - SYCL ------- //
0251   // -------------------------------------- //
0252
0253   auto start_sycl = std::chrono::system_clock::now();
0254
0255   group_count = 0;
0256   std::vector<std::vector<Acts::Seed<SpacePoint>>> seedVector_sycl;
0257
0258   Acts::SpacePointData spacePointData;
0259   spacePointData.resize(spVec.size());
0260
0261   for (auto [bottom, middle, top] : spGroup) {
0262     seedVector_sycl.push_back(syclSeedFinder.createSeedsForGroup(
0263         spacePointData, spGroup.grid(), bottom, middle, top));
0264     group_count++;
0265     if (!cmdlTool.allgroup && group_count >= cmdlTool.groups) {
0266       break;
0267     }
0268   }
0269   auto end_sycl = std::chrono::system_clock::now();
0270
0271   if (!cmdlTool.csvFormat) {
0272     std::cout << "Analyzed " << group_count << " groups for SYCL" << std::endl;
0273   }
0274
0275   std::chrono::duration<double> elapsec_cpu = end_cpu - start_cpu;
0276   double cpuTime = elapsec_cpu.count();
0277
0278   std::chrono::duration<double> elapsec_sycl = end_sycl - start_sycl;
0279   double syclTime = elapsec_sycl.count();
0280
0281   auto textWidth = 20;
0282   auto numWidth = 11;
0283
0284   int nSeed_cpu = 0;
0285   int nSeed_sycl = 0;
0286   int nMatch = 0;
0287
0288   if (cmdlTool.matches && !cmdlTool.onlyGpu) {
0289     for (auto& outVec : seedVector_cpu) {
0290       nSeed_cpu += outVec.size();
0291     }
0292
0293     for (auto& outVec : seedVector_sycl) {
0294       nSeed_sycl += outVec.size();
0295     }
0296
0297     for (std::size_t i = 0; i < seedVector_cpu.size(); i++) {
0298       auto regionVec_cpu = seedVector_cpu[i];
0299       auto regionVec_sycl = seedVector_sycl[i];
0300
0301       std::vector<std::vector<SpacePoint>> seeds_cpu;
0302       std::vector<std::vector<SpacePoint>> seeds_sycl;
0303
0304       for (const auto& sd : regionVec_cpu) {
0305         std::vector<SpacePoint> seed_cpu;
0306         seed_cpu.push_back(*(sd.sp()[0]));
0307         seed_cpu.push_back(*(sd.sp()[1]));
0308         seed_cpu.push_back(*(sd.sp()[2]));
0309         seeds_cpu.push_back(seed_cpu);
0310       }
0311       for (const auto& sd : regionVec_sycl) {
0312         std::vector<SpacePoint> seed_sycl;
0313         seed_sycl.push_back(*(sd.sp()[0]));
0314         seed_sycl.push_back(*(sd.sp()[1]));
0315         seed_sycl.push_back(*(sd.sp()[2]));
0316         seeds_sycl.push_back(seed_sycl);
0317       }
0318
0319       for (auto seed : seeds_cpu) {
0320         for (auto other : seeds_sycl) {
0321           if (seed[0] == other[0] && seed[1] == other[1] &&
0322               seed[2] == other[2]) {
0323             nMatch++;
0324             break;
0325           }
0326         }
0327       }
0328     }
0329   }
0330
0331   if (!cmdlTool.csvFormat) {
0332     std::cout << std::endl;
0333     std::cout
0334         << "------------------------- Time Metric -------------------------"
0335         << std::endl;
0336     std::cout << std::setw(textWidth) << " Device:";
0337     std::cout << std::setw(numWidth) << "CPU";
0338     std::cout << std::setw(numWidth) << "SYCL";
0339     std::cout << std::setw(textWidth) << "Speedup/ Agreement" << std::endl;
0340     std::cout << std::setw(textWidth) << " Time (s):";
0341     std::cout << std::setw(numWidth) << std::to_string(cpuTime);
0342     std::cout << std::setw(numWidth) << std::to_string(syclTime);
0343     std::cout << std::setw(textWidth) << std::to_string(cpuTime / syclTime);
0344     std::cout << std::endl;
0345
0346     if (cmdlTool.matches && !cmdlTool.onlyGpu) {
0347       std::cout << std::setw(textWidth) << " Seeds found:";
0348       std::cout << std::setw(numWidth) << std::to_string(nSeed_cpu);
0349       std::cout << std::setw(numWidth) << std::to_string(nSeed_sycl);
0350       std::cout << std::setw(textWidth)
0351                 << std::to_string(float(nMatch) / float(nSeed_cpu) * 100);
0352       std::cout << std::endl;
0353     }
0354
0355     std::cout
0356         << "---------------------------------------------------------------"
0357         << std::endl;
0358     std::cout << std::endl;
0359   } else {
0360     std::cout << cpuTime << ',' << syclTime << ',' << cpuTime / syclTime << ','
0361               << nSeed_cpu << ',' << nSeed_sycl << ',' << nMatch << '\n';
0362   }
0363
0364   for (const auto* S : spVec) {
0365     delete[] S;
0366   }
0367
0368   return 0;
0369 }