src/TrackFitting/GsfMixtureReduction.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/TrackFitting/GsfMixtureReduction.hpp"
0010
0011 #include "Acts/TrackFitting/detail/GsfComponentMerging.hpp"
0012
0013 #include <algorithm>
0014 #include <format>
0015 #include <iostream>
0016
0017 using namespace Acts;
0018
0019 namespace Acts::detail::Gsf {
0020
0021 namespace {
0022
0023 void reduceWithKLDistanceImpl(std::vector<GsfComponent> &cmpCache,
0024                               std::size_t maxCmpsAfterMerge,
0025                               const Surface &surface) {
0026   SymmetricKLDistanceMatrix distances(cmpCache);
0027
0028   auto remainingComponents = cmpCache.size();
0029
0030   while (remainingComponents > maxCmpsAfterMerge) {
0031     const auto [minI, minJ] = distances.minDistancePair();
0032
0033     // Set one component and compute associated distances
0034     cmpCache[minI] =
0035         mergeTwoComponents(cmpCache[minI], cmpCache[minJ], surface);
0036     distances.recomputeAssociatedDistances(minI, cmpCache);
0037
0038     // Set weight of the other component to -1 so we can remove it later and
0039     // mask its distances
0040     cmpCache[minJ].weight = -1.0;
0041     distances.maskAssociatedDistances(minJ);
0042
0043     remainingComponents--;
0044   }
0045
0046   // Remove all components which are labeled with weight -1
0047   std::ranges::sort(cmpCache, {}, [&](const auto &c) { return c.weight; });
0048   cmpCache.erase(
0049       std::remove_if(cmpCache.begin(), cmpCache.end(),
0050                      [&](const auto &a) { return a.weight == -1.0; }),
0051       cmpCache.end());
0052
0053   assert(cmpCache.size() == maxCmpsAfterMerge && "size mismatch");
0054 }
0055
0056 }  // namespace
0057
0058 }  // namespace Acts::detail::Gsf
0059
0060 void Acts::reduceMixtureLargestWeights(std::vector<GsfComponent> &cmpCache,
0061                                        std::size_t maxCmpsAfterMerge,
0062                                        const Surface & /*surface*/) {
0063   if (cmpCache.size() <= maxCmpsAfterMerge) {
0064     return;
0065   }
0066
0067   std::nth_element(
0068       cmpCache.begin(), cmpCache.begin() + maxCmpsAfterMerge, cmpCache.end(),
0069       [](const auto &a, const auto &b) { return a.weight > b.weight; });
0070   cmpCache.resize(maxCmpsAfterMerge);
0071 }
0072
0073 void Acts::reduceMixtureWithKLDistance(std::vector<GsfComponent> &cmpCache,
0074                                        std::size_t maxCmpsAfterMerge,
0075                                        const Surface &surface) {
0076   if (cmpCache.size() <= maxCmpsAfterMerge) {
0077     return;
0078   }
0079   detail::Gsf::reduceWithKLDistanceImpl(cmpCache, maxCmpsAfterMerge, surface);
0080 }
0081
0082 void Acts::reduceMixtureWithKLDistanceNaive(
0083     std::vector<Acts::GsfComponent> &cmpCache, std::size_t maxCmpsAfterMerge,
0084     const Surface &surface) {
0085   if (cmpCache.size() <= maxCmpsAfterMerge) {
0086     return;
0087   }
0088
0089   while (cmpCache.size() > maxCmpsAfterMerge) {
0090     // Recompute ALL distances every iteration (naive approach)
0091     double minDistance = std::numeric_limits<double>::max();
0092     std::size_t minI = 0;
0093     std::size_t minJ = 0;
0094
0095     for (std::size_t i = 0; i < cmpCache.size(); ++i) {
0096       for (std::size_t j = i + 1; j < cmpCache.size(); ++j) {
0097         const double distance =
0098             detail::Gsf::computeSymmetricKlDivergence(cmpCache[i], cmpCache[j]);
0099         if (distance < minDistance) {
0100           minDistance = distance;
0101           minI = i;
0102           minJ = j;
0103         }
0104       }
0105     }
0106
0107     // Merge the two closest components
0108     cmpCache[minI] = detail::Gsf::mergeTwoComponents(cmpCache[minI],
0109                                                      cmpCache[minJ], surface);
0110
0111     // Remove the merged component immediately
0112     cmpCache.erase(cmpCache.begin() + minJ);
0113   }
0114
0115   assert(cmpCache.size() == maxCmpsAfterMerge && "size mismatch");
0116 }