feat: Add naive component merger for GSF and benchmark #4874

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft

benjaminhuth wants to merge 21 commits into acts-project:main from benjaminhuth:perf/add-cmp-reduction-benchmark

+408 −26

Core/include/Acts/TrackFitting/GsfMixtureReduction.hpp

-Original file line number
+Diff line change
@@ Expand Up @@
                                      std::size_t maxCmpsAfterMerge,
                                      const Surface &surface);
+    /// Naive implementation of component reduction with KL-distance. Recomputes all
+    /// distances in every iteration without any caching or optimization. This
+    /// serves as a baseline for testing and benchmarking the optimized version.
+    /// @param cmpCache the component collection
+    /// @param maxCmpsAfterMerge the number of components we want to reach
+    /// @param surface the surface type on which the components are
+    void reduceMixtureWithKLDistanceNaive(std::vector<GsfComponent> &cmpCache,
+                                          std::size_t maxCmpsAfterMerge,
+                                          const Surface &surface);
     }  // namespace Acts

Core/src/TrackFitting/GsfMixtureReduction.cpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -11,6 +11,10 @@ @@
     #include "Acts/TrackFitting/detail/SymmetricKlDistanceMatrix.hpp"
     #include <algorithm>
+    #include <format>
+    #include <iostream>
+    using namespace Acts;
     template <typename proj_t, typename angle_desc_t>
     void reduceWithKLDistanceImpl(std::vector<Acts::GsfComponent> &cmpCache,
@@ Expand All @@
       while (remainingComponents > maxCmpsAfterMerge) {
         const auto [minI, minJ] = distances.minDistancePair();
+        // auto prevCmpI = cmpCache[minI];
+        // auto prevCmpJ = cmpCache[minJ];
         // Set one component and compute associated distances
         cmpCache[minI] =
             mergeComponents(cmpCache[minI], cmpCache[minJ], proj, desc);
+        /*std::cout << std::format("OPTIM: Merge with phi: {:.3f} + {:.3f} ->
+           {:.3f}, weight: {:.3f} + {:.3f} -> {:.3f}",
+                                 prevCmpI.boundPars[eBoundPhi],
+                                 prevCmpJ.boundPars[eBoundPhi],
+                                 cmpCache[minI].boundPars[eBoundPhi],
+                                 prevCmpI.weight, prevCmpJ.weight,
+                                 cmpCache[minI].weight) << std::endl;*/
         distances.recomputeAssociatedDistances(minI, cmpCache, proj);
         // Set weight of the other component to -1 so we can remove it later and
@@ Expand Down Expand Up @@
       });
     }
+    void reduceMixtureWithKLDistanceNaive(std::vector<Acts::GsfComponent> &cmpCache,
+                                          std::size_t maxCmpsAfterMerge,
+                                          const Surface &surface) {
+      if (cmpCache.size() <= maxCmpsAfterMerge) {
+        return;
+      }
+      auto proj = [](auto &a) -> decltype(auto) { return a; };
+      detail::angleDescriptionSwitch(surface, [&](const auto &desc) {
+        while (cmpCache.size() > maxCmpsAfterMerge) {
+          // Recompute ALL distances every iteration (naive approach)
+          double minDistance = std::numeric_limits<double>::max();
+          std::size_t minI = 0;
+          std::size_t minJ = 0;
+          for (std::size_t i = 0; i < cmpCache.size(); ++i) {
+            for (std::size_t j = i + 1; j < cmpCache.size(); ++j) {
+              double distance = detail::computeSymmetricKlDivergence(
+                  cmpCache[i], cmpCache[j], proj);
+              if (distance < minDistance) {
+                minDistance = distance;
+                minI = i;
+                minJ = j;
+              }
+            }
+          }
+          // auto prevCmpI = cmpCache[minI];
+          // auto prevCmpJ = cmpCache[minJ];
+          // Merge the two closest components
+          cmpCache[minI] =
+              detail::mergeComponents(cmpCache[minI], cmpCache[minJ], proj, desc);
+          /*std::cout << std::format("NAIVE: Merge with phi: {:.3f} + {:.3f} ->
+             {:.3f}, weight: {:.3f} + {:.3f} -> {:.3f}",
+                                   prevCmpI.boundPars[eBoundPhi],
+                                   prevCmpJ.boundPars[eBoundPhi],
+                                   cmpCache[minI].boundPars[eBoundPhi],
+                                   prevCmpI.weight, prevCmpJ.weight,
+                                   cmpCache[minI].weight) << std::endl;*/
+          // Remove the merged component immediately
+          cmpCache.erase(cmpCache.begin() + minJ);
+        }
+      });
+      assert(cmpCache.size() == maxCmpsAfterMerge && "size mismatch");
+    }
     }  // namespace Acts

Tests/Benchmarks/CMakeLists.txt

-Original file line number
+Diff line change
@@ Expand Up / @@ -33,3 +33,4 @@ add_benchmark(SympyStepper SympyStepperBenchmark.cpp) @@
     add_benchmark(Stepper StepperBenchmark.cpp)
     add_benchmark(SourceLink SourceLinkBenchmark.cpp)
     add_benchmark(TrackEdm TrackEdmBenchmark.cpp)
+    add_benchmark(GsfComponentReduction GsfComponentReductionBenchmark.cpp)

Tests/Benchmarks/GsfComponentReductionBenchmark.cpp

-Original file line number
+Diff line change
@@ -0,0 +1,98 @@
+    // This file is part of the ACTS project.
+    //
+    // Copyright (C) 2016 CERN for the benefit of the ACTS project
+    //
+    // This Source Code Form is subject to the terms of the Mozilla Public
+    // License, v. 2.0. If a copy of the MPL was not distributed with this
+    // file, You can obtain one at https://mozilla.org/MPL/2.0/.
+    #include "Acts/TrackFitting/GsfMixtureReduction.hpp"
+    #include <boost/program_options.hpp>
+    #include "StepperBenchmarkCommons.hpp"
+    namespace po = boost::program_options;
+    using namespace Acts;
+    using namespace ActsTests;
+    constexpr std::size_t nComponents = 72;  // 12 max components GSF, 6 material
+    template <typename Fun>
+    void test(const std::vector<std::vector<GsfComponent>> &inputData,
+              const Acts::Surface &surface, Fun &&fun) {
+      auto num_runs = 1000;
+      auto res = microBenchmark(
+          [&](const std::vector<GsfComponent> &input) {
+            // Avoid reallocation cost every iteration
+            static std::vector<GsfComponent> inputCopy(nComponents);
+            inputCopy = input;
+            fun(inputCopy, 12, surface);
+            assumeRead(inputCopy);
+          },
+          inputData, num_runs);
+      std::cout << res << std::endl;
+    }
+    int main(int argc, char *argv[]) {
+      bool weight = true;
+      bool kl_optimized = true;
+      bool kl_naive = true;
+      try {
+        po::options_description desc("Allowed options");
+        desc.add_options()("help", "produce help message")(
+            "weight", po::value<bool>(&weight)->default_value(true),
+            "run weight-based reduction (reduceMixtureLargestWeights)")(
+            "kl-optimized", po::value<bool>(&kl_optimized)->default_value(true),
+            "run optimized KL distance reduction (reduceMixtureWithKLDistance)")(
+            "kl-naive", po::value<bool>(&kl_naive)->default_value(true),
+            "run naive KL distance reduction (reduceMixtureWithKLDistanceNaive)");
+        po::variables_map vm;
+        po::store(po::parse_command_line(argc, argv, desc), vm);
+        po::notify(vm);
+        if (vm.contains("help")) {
+          std::cout << desc << std::endl;
+          return 0;
+        }
+      } catch (std::exception &e) {
+        std::cerr << "error: " << e.what() << std::endl;
+        return 1;
+      }
+      std::vector<std::vector<GsfComponent>> data(100);
+      auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+          Acts::Transform3::Identity());
+      for (auto &sample : data) {
+        sample.reserve(nComponents);
+        for (auto i = 0ul; i < nComponents; ++i) {
+          double weight_val = 1.0 / nComponents;
+          auto cov = Acts::BoundMatrix::Identity();
+          auto pars = Acts::BoundVector::Random();
+          sample.push_back({weight_val, pars, cov});
+        }
+      }
+      if (weight) {
+        std::cout << "reduceMixtureLargestWeights" << std::endl;
+        test(data, *surface, reduceMixtureLargestWeights);
+        std::cout << std::endl;
+      }
+      if (kl_optimized) {
+        std::cout << "reduceMixtureWithKLDistance (optimized)" << std::endl;
+        test(data, *surface, reduceMixtureWithKLDistance);
+        std::cout << std::endl;
+      }
+      if (kl_naive) {
+        std::cout << "reduceMixtureWithKLDistanceNaive (baseline)" << std::endl;
+        test(data, *surface, reduceMixtureWithKLDistanceNaive);
+        std::cout << std::endl;
+      }
+    }

Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -379,6 +379,95 @@ BOOST_AUTO_TEST_CASE(test_perigee_surface) { @@
       test_surface(*surface, desc, p, 1.1);
     }
+    BOOST_AUTO_TEST_CASE(test_phi_angle_wrapping_simple) {
+      // Two components with phi angles spanning the -pi/+pi boundary
+      std::vector<DummyComponent<1>> cmps;
+      // Component 1: phi = +2.8 rad (about 160 degrees), weight = 0.5
+      DummyComponent<1> c1;
+      c1.weight = 0.5;
+      c1.boundPars[0] = 2.8;
+      c1.boundCov = ActsSquareMatrix<1>::Identity();
+      cmps.push_back(c1);
+      // Component 2: phi = -2.8 rad (about -160 degrees), weight = 0.5
+      DummyComponent<1> c2;
+      c2.weight = 0.5;
+      c2.boundPars[0] = -2.8;
+      c2.boundCov = ActsSquareMatrix<1>::Identity();
+      cmps.push_back(c2);
+      // Expected: mean should be +/-pi (going through discontinuity), NOT 0.0
+      // With cyclic angle descriptor
+      const auto desc = std::tuple<detail::CyclicAngle<eBoundLoc0>>{};
+      auto [mean, cov] =
+          detail::gaussianMixtureMeanCov(cmps, std::identity{}, desc);
+      // Check mean is close to +/-pi, not 0
+      BOOST_CHECK(std::abs(std::abs(mean[0]) - std::numbers::pi) < 0.1);
+      // Without cyclic angle descriptor - should give wrong result
+      auto [meanWrong, covWrong] =
+          detail::gaussianMixtureMeanCov(cmps, std::identity{}, std::tuple<>{});
+      // This will likely be close to 0.0 (arithmetic mean), which is wrong
+      BOOST_CHECK(std::abs(meanWrong[0]) < 0.5);
+    }
+    BOOST_AUTO_TEST_CASE(test_mean_shuffle) {
+      std::shared_ptr<PlaneSurface> surface =
+          CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
+      const std::size_t NComps = 10;
+      std::mt19937 rng(43);
+      std::uniform_real_distribution<double> weightDist(0.5, 1.5);
+      std::uniform_real_distribution<double> loc0Dist(-10.0, 10.0);
+      std::uniform_real_distribution<double> loc1Dist(-10.0, 10.0);
+      std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
+                                                     std::numbers::pi);
+      std::uniform_real_distribution<double> thetaDist(0.0, std::numbers::pi);
+      std::uniform_real_distribution<double> qopDist(0.1, 5.0);
+      std::vector<DummyComponent<5>> cmps;
+      double weightSum = 0.0;
+      for (auto i = 0ul; i < NComps; ++i) {
+        DummyComponent<5> cmp;
+        cmp.boundPars[eBoundLoc0] = loc0Dist(rng);
+        cmp.boundPars[eBoundLoc1] = loc1Dist(rng);
+        cmp.boundPars[eBoundPhi] = phiDist(rng);
+        cmp.boundPars[eBoundTheta] = thetaDist(rng);
+        cmp.boundPars[eBoundQOverP] = qopDist(rng);
+        cmp.weight = weightDist(rng);
+        cmp.boundCov = ActsSquareMatrix<5>::Identity();
+        weightSum += cmp.weight;
+        cmps.push_back(cmp);
+      }
+      for (auto &cmp : cmps) {
+        cmp.weight /= weightSum;
+      }
+      std::ranges::sort(cmps, {}, [](const auto &c) { return c.weight; });
+      const auto desc = detail::AngleDescription<Surface::Plane>::Desc{};
+      const auto meanBefore =
+          std::get<0>(detail::gaussianMixtureMeanCov(cmps, std::identity{}, desc));
+      for (int i = 0; i < 20; ++i) {
+        std::cout << "Shuffle iteration " << i + 1 << std::endl;
+        std::shuffle(cmps.begin(), cmps.end(), rng);
+        const auto meanAfter = std::get<0>(
+            detail::gaussianMixtureMeanCov(cmps, std::identity{}, desc));
+        BOOST_CHECK_CLOSE(meanBefore[eBoundLoc0], meanAfter[eBoundLoc0], 1.e-8);
+        BOOST_CHECK_CLOSE(meanBefore[eBoundLoc1], meanAfter[eBoundLoc1], 1.e-8);
+        BOOST_CHECK_CLOSE(meanBefore[eBoundPhi], meanAfter[eBoundPhi], 1.e-8);
+        BOOST_CHECK_CLOSE(meanBefore[eBoundTheta], meanAfter[eBoundTheta], 1.e-8);
+        BOOST_CHECK_CLOSE(meanBefore[eBoundQOverP], meanAfter[eBoundQOverP], 1.e-8);
+      }
+    }
     BOOST_AUTO_TEST_SUITE_END()
     }  // namespace ActsTests

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

feat: Add naive component merger for GSF and benchmark #4874

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!

feat: Add naive component merger for GSF and benchmark #4874

Are you sure you want to change the base?

feat: Add naive component merger for GSF and benchmark #4874

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Uh oh!