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
0037     // Set weight of the other component to -1 so we can remove it later, and
0038     // remove it from the active set (this must happen before recomputing
0039     // minI's distances since minI's internal slot may move as part of the
0040     // removal's swap-to-last compaction)
0041     cmpCache[minJ].weight = -1.0;
0042     distances.maskAssociatedDistances(minJ);
0043     distances.recomputeAssociatedDistances(minI, cmpCache);
0044
0045     remainingComponents--;
0046   }
0047
0048   // Remove all components which are labeled with weight -1
0049   std::ranges::sort(cmpCache, {}, [&](const auto &c) { return c.weight; });
0050   cmpCache.erase(
0051       std::remove_if(cmpCache.begin(), cmpCache.end(),
0052                      [&](const auto &a) { return a.weight == -1.0; }),
0053       cmpCache.end());
0054
0055   assert(cmpCache.size() == maxCmpsAfterMerge && "size mismatch");
0056 }
0057
0058 }  // namespace
0059
0060 }  // namespace Acts::detail::Gsf
0061
0062 void Acts::reduceMixtureLargestWeights(std::vector<GsfComponent> &cmpCache,
0063                                        std::size_t maxCmpsAfterMerge,
0064                                        const Surface & /*surface*/) {
0065   if (cmpCache.size() <= maxCmpsAfterMerge) {
0066     return;
0067   }
0068
0069   std::nth_element(
0070       cmpCache.begin(), cmpCache.begin() + maxCmpsAfterMerge, cmpCache.end(),
0071       [](const auto &a, const auto &b) { return a.weight > b.weight; });
0072   cmpCache.resize(maxCmpsAfterMerge);
0073 }
0074
0075 void Acts::reduceMixtureWithKLDistance(std::vector<GsfComponent> &cmpCache,
0076                                        std::size_t maxCmpsAfterMerge,
0077                                        const Surface &surface) {
0078   if (cmpCache.size() <= maxCmpsAfterMerge) {
0079     return;
0080   }
0081   detail::Gsf::reduceWithKLDistanceImpl(cmpCache, maxCmpsAfterMerge, surface);
0082 }
0083
0084 void Acts::reduceMixtureWithKLDistanceNaive(
0085     std::vector<Acts::GsfComponent> &cmpCache, std::size_t maxCmpsAfterMerge,
0086     const Surface &surface) {
0087   if (cmpCache.size() <= maxCmpsAfterMerge) {
0088     return;
0089   }
0090
0091   while (cmpCache.size() > maxCmpsAfterMerge) {
0092     // Recompute ALL distances every iteration (naive approach)
0093     double minDistance = std::numeric_limits<double>::max();
0094     std::size_t minI = 0;
0095     std::size_t minJ = 0;
0096
0097     for (std::size_t i = 0; i < cmpCache.size(); ++i) {
0098       for (std::size_t j = i + 1; j < cmpCache.size(); ++j) {
0099         const double distance =
0100             detail::Gsf::computeSymmetricKlDivergence(cmpCache[i], cmpCache[j]);
0101         if (distance < minDistance) {
0102           minDistance = distance;
0103           minI = i;
0104           minJ = j;
0105         }
0106       }
0107     }
0108
0109     // Merge the two closest components
0110     cmpCache[minI] = detail::Gsf::mergeTwoComponents(cmpCache[minI],
0111                                                      cmpCache[minJ], surface);
0112
0113     // Remove the merged component immediately
0114     cmpCache.erase(cmpCache.begin() + minJ);
0115   }
0116
0117   assert(cmpCache.size() == maxCmpsAfterMerge && "size mismatch");
0118 }