Core/TrackFitting/GsfMixtureReductionTests.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 <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0015 #include "Acts/Surfaces/PlaneSurface.hpp"
0016 #include "Acts/TrackFitting/GsfMixtureReduction.hpp"
0017 #include "Acts/TrackFitting/detail/GsfComponentMerging.hpp"
0018 #include "Acts/Utilities/Zip.hpp"
0019
0020 #include <algorithm>
0021 #include <cstddef>
0022 #include <memory>
0023 #include <numbers>
0024 #include <numeric>
0025 #include <random>
0026 #include <tuple>
0027 #include <vector>
0028
0029 using namespace Acts;
0030 using namespace Acts::UnitLiterals;
0031
0032 namespace {
0033
0034 std::tuple<BoundVector, double> computeMeanAndSumOfWeights(
0035     const std::vector<GsfComponent> &cmps, const Surface &surface) {
0036   if (cmps.empty()) {
0037     return {BoundVector::Zero(), 0.0};
0038   }
0039
0040   const double sumOfWeights = std::accumulate(
0041       cmps.begin(), cmps.end(), 0.0,
0042       [](auto sum, const auto &cmp) { return sum + cmp.weight; });
0043
0044   BoundVector mean;
0045   detail::Gsf::angleDescriptionSwitch(surface, [&](const auto &desc) {
0046     auto [meanTmp, cov] = detail::Gsf::mergeGaussianMixtureMeanCov(
0047         cmps,
0048         [](const auto &c) {
0049           return std::tie(c.weight, c.boundPars, c.boundCov);
0050         },
0051         desc);
0052     mean = meanTmp;
0053   });
0054
0055   return std::make_tuple(mean, sumOfWeights);
0056 }
0057
0058 GsfComponent makeDefaultComponent(double weight) {
0059   GsfComponent cmp;
0060   cmp.boundPars = BoundVector::Zero();
0061   cmp.boundCov = BoundMatrix::Identity();
0062   cmp.weight = weight;
0063   return cmp;
0064 }
0065
0066 void testReductionEquivalence(std::vector<GsfComponent> &cmpsRef,
0067                               std::vector<GsfComponent> &cmpsTest) {
0068   BOOST_REQUIRE(cmpsRef.size() == cmpsTest.size());
0069
0070   // sort by weight since the order of components is not guaranteed to be the
0071   // same after reduction
0072   std::ranges::sort(cmpsRef, {}, [](const auto &c) { return c.weight; });
0073   std::ranges::sort(cmpsTest, {}, [](const auto &c) { return c.weight; });
0074
0075   constexpr static double tol = 1.e-8;
0076
0077   // Compare components
0078   for (const auto &[ref, test] : Acts::zip(cmpsRef, cmpsTest)) {
0079     BOOST_CHECK_CLOSE(ref.weight, test.weight, tol);
0080     BOOST_CHECK_CLOSE(ref.boundPars[eBoundLoc0], test.boundPars[eBoundLoc0],
0081                       tol);
0082     BOOST_CHECK_CLOSE(ref.boundPars[eBoundLoc1], test.boundPars[eBoundLoc1],
0083                       tol);
0084     BOOST_CHECK_CLOSE(ref.boundPars[eBoundPhi], test.boundPars[eBoundPhi], tol);
0085     BOOST_CHECK_CLOSE(ref.boundPars[eBoundTheta], test.boundPars[eBoundTheta],
0086                       tol);
0087     BOOST_CHECK_CLOSE(ref.boundPars[eBoundQOverP], test.boundPars[eBoundQOverP],
0088                       tol);
0089     BOOST_CHECK_CLOSE(ref.boundPars[eBoundTime], test.boundPars[eBoundTime],
0090                       tol);
0091   }
0092 }
0093
0094 }  // namespace
0095
0096 namespace ActsTests {
0097
0098 BOOST_AUTO_TEST_SUITE(TrackFittingSuite)
0099
0100 BOOST_AUTO_TEST_CASE(test_distance_matrix_min_distance) {
0101   std::vector<GsfComponent> cmps = {
0102       {1. / 3., BoundVector::Constant(-2.), BoundMatrix::Identity()},
0103       {1. / 3., BoundVector::Constant(+0.), BoundMatrix::Identity()},
0104       {1. / 3., BoundVector::Constant(+1.), BoundMatrix::Identity()},
0105       {1. / 3., BoundVector::Constant(+4.), BoundMatrix::Identity()}};
0106
0107   detail::Gsf::SymmetricKLDistanceMatrix mat(cmps);
0108
0109   const auto [i, j] = mat.minDistancePair();
0110   BOOST_CHECK_EQUAL(std::min(i, j), 1);
0111   BOOST_CHECK_EQUAL(std::max(i, j), 2);
0112 }
0113
0114 BOOST_AUTO_TEST_CASE(test_distance_matrix_masking) {
0115   std::vector<GsfComponent> cmps = {
0116       {1. / 3., BoundVector::Constant(-2.), BoundMatrix::Identity()},
0117       {1. / 3., BoundVector::Constant(+0.), BoundMatrix::Identity()},
0118       {1. / 3., BoundVector::Constant(+1.), BoundMatrix::Identity()},
0119       {1. / 3., BoundVector::Constant(+4.), BoundMatrix::Identity()}};
0120
0121   const std::size_t cmp_to_mask = 2;
0122
0123   detail::Gsf::SymmetricKLDistanceMatrix mat_full(cmps);
0124   mat_full.maskAssociatedDistances(cmp_to_mask);
0125
0126   cmps.erase(cmps.begin() + cmp_to_mask);
0127   detail::Gsf::SymmetricKLDistanceMatrix mat_small(cmps);
0128
0129   const auto [full_i, full_j] = mat_full.minDistancePair();
0130   const auto [small_i, small_j] = mat_small.minDistancePair();
0131
0132   BOOST_CHECK_EQUAL(std::min(full_i, full_j), 0);
0133   BOOST_CHECK_EQUAL(std::max(full_i, full_j), 1);
0134   BOOST_CHECK_EQUAL(full_i, small_i);
0135   BOOST_CHECK_EQUAL(full_j, small_j);
0136 }
0137
0138 BOOST_AUTO_TEST_CASE(test_distance_matrix_recompute_distance) {
0139   std::vector<GsfComponent> cmps = {
0140       {1. / 3., BoundVector::Constant(-2.), BoundMatrix::Identity()},
0141       {1. / 3., BoundVector::Constant(+0.), BoundMatrix::Identity()},
0142       {1. / 3., BoundVector::Constant(+1.), BoundMatrix::Identity()},
0143       {1. / 3., BoundVector::Constant(+4.), BoundMatrix::Identity()}};
0144
0145   detail::Gsf::SymmetricKLDistanceMatrix mat(cmps);
0146
0147   {
0148     const auto [i, j] = mat.minDistancePair();
0149     BOOST_CHECK_EQUAL(std::min(i, j), 1);
0150     BOOST_CHECK_EQUAL(std::max(i, j), 2);
0151   }
0152
0153   cmps[3].boundPars = BoundVector::Constant(0.1);
0154   mat.recomputeAssociatedDistances(3, cmps);
0155
0156   {
0157     const auto [i, j] = mat.minDistancePair();
0158     BOOST_CHECK_EQUAL(std::min(i, j), 1);
0159     BOOST_CHECK_EQUAL(std::max(i, j), 3);
0160   }
0161
0162   cmps[0].boundPars = BoundVector::Constant(1.01);
0163   mat.recomputeAssociatedDistances(0, cmps);
0164
0165   {
0166     const auto [i, j] = mat.minDistancePair();
0167     BOOST_CHECK_EQUAL(std::min(i, j), 0);
0168     BOOST_CHECK_EQUAL(std::max(i, j), 2);
0169   }
0170 }
0171
0172 BOOST_AUTO_TEST_CASE(test_mixture_reduction) {
0173   // Assume that the components are on a generic plane surface
0174   std::shared_ptr<PlaneSurface> surface =
0175       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0176   const std::size_t NComps = 4;
0177   std::vector<GsfComponent> cmps;
0178
0179   for (auto i = 0ul; i < NComps; ++i) {
0180     GsfComponent a;
0181     a.boundPars = BoundVector::Zero();
0182     a.boundCov = BoundMatrix::Identity();
0183     a.weight = 1.0 / NComps;
0184     cmps.push_back(a);
0185   }
0186
0187   cmps[0].boundPars[eBoundQOverP] = 0.5_GeV;
0188   cmps[1].boundPars[eBoundQOverP] = 1.5_GeV;
0189   cmps[2].boundPars[eBoundQOverP] = 3.5_GeV;
0190   cmps[3].boundPars[eBoundQOverP] = 4.5_GeV;
0191
0192   // Check start properties
0193   const auto [mean0, sumOfWeights0] =
0194       computeMeanAndSumOfWeights(cmps, *surface);
0195
0196   BOOST_CHECK_CLOSE(mean0[eBoundQOverP], 2.5_GeV, 1.e-8);
0197   BOOST_CHECK_CLOSE(sumOfWeights0, 1.0, 1.e-8);
0198
0199   // Reduce by factor of 2 and check if weights and QoP are correct
0200   reduceMixtureWithKLDistance(cmps, 2, *surface);
0201
0202   BOOST_CHECK_EQUAL(cmps.size(), 2);
0203
0204   std::ranges::sort(cmps, {},
0205                     [](const auto &c) { return c.boundPars[eBoundQOverP]; });
0206   BOOST_CHECK_CLOSE(cmps[0].boundPars[eBoundQOverP], 1.0_GeV, 1.e-8);
0207   BOOST_CHECK_CLOSE(cmps[1].boundPars[eBoundQOverP], 4.0_GeV, 1.e-8);
0208
0209   const auto [mean1, sumOfWeights1] =
0210       computeMeanAndSumOfWeights(cmps, *surface);
0211
0212   BOOST_CHECK_CLOSE(mean1[eBoundQOverP], 2.5_GeV, 1.e-8);
0213   BOOST_CHECK_CLOSE(sumOfWeights1, 1.0, 1.e-8);
0214
0215   // Reduce by factor of 2 and check if weights and QoP are correct
0216   reduceMixtureWithKLDistance(cmps, 1, *surface);
0217
0218   BOOST_CHECK_EQUAL(cmps.size(), 1);
0219   BOOST_CHECK_CLOSE(cmps[0].boundPars[eBoundQOverP], 2.5_GeV, 1.e-8);
0220   BOOST_CHECK_CLOSE(cmps[0].weight, 1.0, 1.e-8);
0221 }
0222
0223 BOOST_AUTO_TEST_CASE(test_weight_cut_reduction) {
0224   std::shared_ptr<PlaneSurface> dummy =
0225       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0226   std::vector<GsfComponent> cmps;
0227
0228   // weights do not need to be normalized for this test
0229   for (auto w : {1.0, 2.0, 3.0, 4.0}) {
0230     GsfComponent a;
0231     a.boundPars = BoundVector::Zero();
0232     a.boundCov = BoundMatrix::Identity();
0233     a.weight = w;
0234     cmps.push_back(a);
0235   }
0236
0237   reduceMixtureLargestWeights(cmps, 2, *dummy);
0238
0239   BOOST_CHECK_EQUAL(cmps.size(), 2);
0240   std::ranges::sort(cmps, {}, [](const auto &c) { return c.weight; });
0241
0242   BOOST_CHECK_EQUAL(cmps[0].weight, 3.0);
0243   BOOST_CHECK_EQUAL(cmps[1].weight, 4.0);
0244 }
0245
0246 BOOST_AUTO_TEST_CASE(test_naive_vs_optimized) {
0247   std::shared_ptr<PlaneSurface> surface =
0248       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0249   const std::size_t NComps = 10;
0250
0251   std::mt19937 rng(42);
0252   std::uniform_real_distribution<double> weightDist(0.5, 1.5);
0253   std::uniform_real_distribution<double> loc0Dist(-10.0, 10.0);
0254   std::uniform_real_distribution<double> loc1Dist(-10.0, 10.0);
0255   std::uniform_real_distribution<double> phiDist(-std::numbers::pi,
0256                                                  std::numbers::pi);
0257   std::uniform_real_distribution<double> thetaDist(0.0, std::numbers::pi);
0258   std::uniform_real_distribution<double> qopDist(0.1, 5.0);
0259
0260   std::vector<GsfComponent> cmps;
0261   double weightSum = 0.0;
0262
0263   for (auto i = 0ul; i < NComps; ++i) {
0264     GsfComponent cmp = makeDefaultComponent(weightDist(rng));
0265     cmp.boundPars[eBoundLoc0] = loc0Dist(rng);
0266     cmp.boundPars[eBoundLoc1] = loc1Dist(rng);
0267     cmp.boundPars[eBoundPhi] = phiDist(rng);
0268     cmp.boundPars[eBoundTheta] = thetaDist(rng);
0269     cmp.boundPars[eBoundQOverP] = qopDist(rng);
0270     cmp.boundPars[eBoundTime] = 0.0;
0271     weightSum += cmp.weight;
0272     cmps.push_back(cmp);
0273   }
0274
0275   for (auto &cmp : cmps) {
0276     cmp.weight /= weightSum;
0277   }
0278
0279   for (std::size_t targetSize = 9; targetSize > 0; --targetSize) {
0280     std::vector<GsfComponent> cmpsOptimized = cmps;
0281     std::vector<GsfComponent> cmpsNaive = cmps;
0282
0283     reduceMixtureWithKLDistance(cmpsOptimized, targetSize, *surface);
0284     reduceMixtureWithKLDistanceNaive(cmpsNaive, targetSize, *surface);
0285     BOOST_CHECK_EQUAL(cmpsNaive.size(), targetSize);
0286
0287     testReductionEquivalence(cmpsNaive, cmpsOptimized);
0288   }
0289 }
0290
0291 BOOST_AUTO_TEST_CASE(test_exact_tie_breaking) {
0292   // Four exactly identical components: every pairwise KL distance is exactly
0293   // 0, so every merge step is a tie. The optimized reducer does not attempt
0294   // to match naive's specific tie-break convention (which pair is merged
0295   // first among exact ties is otherwise unspecified) -- but reducing to a
0296   // single component here still must equal naive's result regardless of
0297   // merge order, since merging any subset of exactly-identical Gaussians is
0298   // a commutative, associative operation (the weighted mean/covariance of
0299   // identical components doesn't depend on pairing order). This is a
0300   // regression guard for the underlying merge math, not for tie-break
0301   // ordering.
0302   std::shared_ptr<PlaneSurface> surface =
0303       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0304
0305   std::vector<GsfComponent> cmps;
0306   for (int i = 0; i < 4; ++i) {
0307     GsfComponent cmp = makeDefaultComponent(0.25);
0308     cmp.boundPars[eBoundQOverP] = 1.0;
0309     cmp.boundCov(eBoundQOverP, eBoundQOverP) = 1.0;
0310     cmps.push_back(cmp);
0311   }
0312
0313   std::vector<GsfComponent> cmpsNaive = cmps;
0314   reduceMixtureWithKLDistanceNaive(cmpsNaive, 1, *surface);
0315
0316   std::vector<GsfComponent> cmpsOptimized = cmps;
0317   reduceMixtureWithKLDistance(cmpsOptimized, 1, *surface);
0318
0319   testReductionEquivalence(cmpsNaive, cmpsOptimized);
0320 }
0321
0322 BOOST_AUTO_TEST_CASE(test_naive_vs_optimized_stress) {
0323   std::shared_ptr<PlaneSurface> surface =
0324       CurvilinearSurface(Vector3{0, 0, 0}, Vector3{1, 0, 0}).planeSurface();
0325
0326   std::mt19937 rng(1234);
0327   std::uniform_real_distribution<double> weightDist(0.5, 1.5);
0328   std::uniform_real_distribution<double> qopDist(0.1, 5.0);
0329   std::uniform_real_distribution<double> covDist(0.5, 2.0);
0330
0331   for (int seed = 0; seed < 500; ++seed) {
0332     const std::size_t NComps = 72;
0333     std::vector<GsfComponent> cmps;
0334     double weightSum = 0.0;
0335
0336     for (auto i = 0ul; i < NComps; ++i) {
0337       GsfComponent cmp = makeDefaultComponent(weightDist(rng));
0338       // Force frequent exact ties/duplicates by quantizing qop and cov, and
0339       // by making every third component an exact duplicate of a previous
0340       // one, similar to how real GSF material-mixture components can be
0341       // near/exactly-degenerate.
0342       double qop = qopDist(rng);
0343       if (i % 3 == 0 && !cmps.empty()) {
0344         cmp.boundPars[eBoundQOverP] = cmps[i / 2].boundPars[eBoundQOverP];
0345         cmp.boundCov(eBoundQOverP, eBoundQOverP) =
0346             cmps[i / 2].boundCov(eBoundQOverP, eBoundQOverP);
0347       } else {
0348         cmp.boundPars[eBoundQOverP] = qop;
0349         cmp.boundCov(eBoundQOverP, eBoundQOverP) = covDist(rng);
0350       }
0351       weightSum += cmp.weight;
0352       cmps.push_back(cmp);
0353     }
0354
0355     for (auto &cmp : cmps) {
0356       cmp.weight /= weightSum;
0357     }
0358
0359     for (std::size_t targetSize : {1ul, 2ul, 12ul, 40ul, 71ul}) {
0360       std::vector<GsfComponent> cmpsOptimized = cmps;
0361       std::vector<GsfComponent> cmpsNaive = cmps;
0362
0363       reduceMixtureWithKLDistance(cmpsOptimized, targetSize, *surface);
0364       reduceMixtureWithKLDistanceNaive(cmpsNaive, targetSize, *surface);
0365       BOOST_REQUIRE_EQUAL(cmpsNaive.size(), targetSize);
0366
0367       testReductionEquivalence(cmpsNaive, cmpsOptimized);
0368     }
0369   }
0370 }
0371
0372 BOOST_AUTO_TEST_SUITE_END()
0373
0374 }  // namespace ActsTests