EIC code displayed by LXR master/​include/​Acts/​TrackFitting/​detail/​GsfComponentMerging.hpp	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-11-05 09:57:57

 // 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/.
 
 #pragma once
 
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
 #include "Acts/TrackFitting/GsfComponent.hpp"
 #include "Acts/TrackFitting/GsfOptions.hpp"
 #include "Acts/Utilities/detail/periodic.hpp"
 
 #include <cstddef>
 #include <iosfwd>
 #include <stdexcept>
 #include <tuple>
 #include <utility>
 #include <vector>
 
 namespace Acts::detail::Gsf {
 
 /// Reduce Gaussian mixture into a single parameter vector and covariance, using
 /// the specified method to reduce the mixture.
 ///
 /// @param mixture The mixture iterable
 /// @param surface The surface, on which the bound state is
 /// @param method How to reduce the mixture
 /// @return parameters and covariance as tuple
 std::tuple<BoundVector, BoundMatrix> mergeGaussianMixture(
     std::span<const GsfComponent> mixture, const Surface &surface,
     ComponentMergeMethod method);
 
 /// Merge two Gaussian mixture components into one
 ///
 /// @param a First component
 /// @param b Second component
 /// @param surface The surface, on which the bound state is
 /// @return merged component
 GsfComponent mergeTwoComponents(const GsfComponent &a, const GsfComponent &b,
                                 const Surface &surface);
 
 /// Angle descriptions for the combineBoundGaussianMixture function
 template <BoundIndices Idx>
 struct CyclicAngle {
   constexpr static BoundIndices idx = Idx;
   constexpr static double constant = 1.0;
 };
 
 template <BoundIndices Idx>
 struct CyclicRadiusAngle {
   constexpr static BoundIndices idx = Idx;
   double constant = 1.0;  // the radius
 };
 
 /// A compile time map to provide angle descriptions for different surfaces
 //! [angle description]
 template <Surface::SurfaceType type_t>
 struct AngleDescription {
   using Desc = std::tuple<CyclicAngle<eBoundPhi>>;
 };
 
 template <>
 struct AngleDescription<Surface::Disc> {
   using Desc = std::tuple<CyclicAngle<eBoundLoc1>, CyclicAngle<eBoundPhi>>;
 };
 
 template <>
 struct AngleDescription<Surface::Cylinder> {
   using Desc =
       std::tuple<CyclicRadiusAngle<eBoundLoc0>, CyclicAngle<eBoundPhi>>;
 };
 //! [angle description]
 
 /// Helper function that encapsulates a switch-case to select the correct angle
 /// description dependent on the surface
 /// @param surface The surface, which the bound state is on
 /// @param callable A callable that takes the angle description as argument
 /// @return The return value of the callable
 template <typename Callable>
 auto angleDescriptionSwitch(const Surface &surface, Callable &&callable) {
   switch (surface.type()) {
     case Surface::Cylinder: {
       AngleDescription<Surface::Cylinder>::Desc desc{};
       const auto &bounds =
           static_cast<const CylinderSurface &>(surface).bounds();
       std::get<0>(desc).constant = bounds.get(CylinderBounds::eR);
       return callable(desc);
     }
     case Surface::Disc: {
       AngleDescription<Surface::Disc>::Desc desc{};
       return callable(desc);
     }
     default: {
       AngleDescription<Surface::Plane>::Desc desc{};
       return callable(desc);
     }
   }
 }
 
 /// Combine multiple components into one representative parameter object. The
 /// function takes a range and a projector to allow for arbitrary input to be
 /// combined.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight and parameters from the components
 /// @param angleDesc The angle description object
 /// @return parameters
 template <typename component_range_t, typename projector_t,
           typename angle_desc_t>
 BoundVector mergeGaussianMixtureMean(const component_range_t &cmps,
                                      const projector_t &projector,
                                      const angle_desc_t &angleDesc) {
   // Early return in case of range with length 1
   if (cmps.size() == 1) {
     const auto [weight_l, pars_l] = projector(cmps.front());
     return pars_l;
   }
 
   // Do the (circular) mean with complex arithmetic.
   // For normal values, just keep the real values. For angles, use the complex
   // phase. Weighting then transparently happens by multiplying a real-valued
   // weight.
   //! [circular mean]
   using CVec = Eigen::Matrix<std::complex<double>, eBoundSize, 1>;
   CVec cMean = CVec::Zero();
   double sumOfWeights = 0;
 
   for (const auto &cmp : cmps) {
     const auto [weight_l, pars_l] = projector(cmp);
     CVec cPars_l = pars_l;
 
     const auto setPolar = [&](const auto &desc) {
       cPars_l[desc.idx] = std::polar(1.0, pars_l[desc.idx] / desc.constant);
     };
     std::apply([&](auto... dsc) { (setPolar(dsc), ...); }, angleDesc);
 
     sumOfWeights += weight_l;
     cMean += weight_l * cPars_l;
   }
 
   cMean /= sumOfWeights;
 
   BoundVector mean = cMean.real();
 
   const auto getArg = [&](const auto &desc) {
     mean[desc.idx] = desc.constant * std::arg(cMean[desc.idx]);
   };
   std::apply([&](auto... dsc) { (getArg(dsc), ...); }, angleDesc);
   //! [circular mean]
 
   return mean;
 }
 
 /// Combine multiple components into one representative parameter object. The
 /// function takes a range and a projector to allow for arbitrary input to be
 /// combined.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight and parameters from the components
 /// @return parameters
 template <typename component_range_t, typename projector_t>
 BoundVector mergeGaussianMixtureMaxWeight(const component_range_t &cmps,
                                           const projector_t &projector) {
   const auto maxWeightIt =
       std::ranges::max_element(cmps, {}, [&](const auto &cmp) {
         const auto [weight_l, pars_l] = projector(cmp);
         return weight_l;
       });
   const auto [weight_max, pars_max] = projector(*maxWeightIt);
   return pars_max;
 }
 
 /// Compute the covariance of a Gaussian mixture given the mean. The function
 /// takes a range and a projector to allow for arbitrary input to be combined.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight and parameters from the components
 /// @param mean The precomputed mean of the mixture
 /// @param angleDesc The angle description object
 /// @return covariance
 template <typename component_range_t, typename projector_t,
           typename angle_desc_t>
 BoundMatrix mergeGaussianMixtureCov(const component_range_t &cmps,
                                     const projector_t &projector,
                                     const BoundVector &mean,
                                     const angle_desc_t &angleDesc) {
   if (cmps.size() == 1) {
     const auto [weight_l, pars_l, cov_l] = projector(cmps.front());
     return cov_l;
   }
 
   BoundMatrix cov = BoundMatrix::Zero();
   double sumOfWeights = 0;
 
   for (const auto &cmp : cmps) {
     const auto [weight_l, pars_l, cov_l] = projector(cmp);
 
     cov += weight_l * cov_l;  // MARK: fpeMask(FLTUND, 1, #2347)
 
     BoundVector diff = pars_l - mean;
 
     // Apply corrections for cyclic coordinates
     const auto handleCyclicCov = [&l = pars_l, &m = mean,
                                   &diff = diff](const auto &desc) {
       diff[desc.idx] = difference_periodic(l[desc.idx] / desc.constant,
                                            m[desc.idx] / desc.constant,
                                            2 * std::numbers::pi) *
                        desc.constant;
     };
     std::apply([&](auto... dsc) { (handleCyclicCov(dsc), ...); }, angleDesc);
 
     sumOfWeights += weight_l;
     cov += weight_l * diff * diff.transpose();
   }
 
   cov /= sumOfWeights;
 
   return cov;
 }
 
 /// Combine multiple components into one representative parameter object. The
 /// function takes a range and a projector to allow for arbitrary input to be
 /// combined.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight, parameters and covariance
 ///        from the components
 /// @param angleDesc The angle description object
 /// @return parameters and covariance
 template <typename component_range_t, typename projector_t,
           typename angle_desc_t>
 std::tuple<BoundVector, BoundMatrix> mergeGaussianMixtureMeanCov(
     const component_range_t &cmps, const projector_t &projector,
     const angle_desc_t &angleDesc) {
   if (cmps.size() == 1) {
     const auto [weight_l, pars_l, cov_l] = projector(cmps.front());
     return {pars_l, cov_l};
   }
 
   const BoundVector mean = mergeGaussianMixtureMean(
       cmps,
       [&](const auto &c) {
         const auto [weight_l, pars_l, cov_l] = projector(c);
         return std::tuple(weight_l, pars_l);
       },
       angleDesc);
 
   // MARK: fpeMaskBegin(FLTUND, 1, #2347)
   const BoundMatrix cov =
       mergeGaussianMixtureCov(cmps, projector, mean, angleDesc);
   // MARK: fpeMaskEnd(FLTUND)
 
   return {mean, cov};
 }
 
 /// Reduce Gaussian mixture into a single parameter vector, using the specified
 /// method to reduce the mixture.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight and parameters from the components
 /// @param angleDesc The angle description object
 /// @param method How to reduce the mixture
 /// @return parameters
 template <typename component_range_t, typename projector_t,
           typename angle_desc_t>
 BoundVector mergeGaussianMixtureParams(const component_range_t &cmps,
                                        const projector_t &projector,
                                        const angle_desc_t &angleDesc,
                                        ComponentMergeMethod method) {
   if (method == ComponentMergeMethod::eMean) {
     return mergeGaussianMixtureMean(cmps, projector, angleDesc);
   } else if (method == ComponentMergeMethod::eMaxWeight) {
     return mergeGaussianMixtureMaxWeight(cmps, projector);
   } else {
     throw std::logic_error("Invalid component merge method");
   }
 }
 
 /// Reduce Gaussian mixture into a single parameter vector and covariance, using
 /// the specified method to reduce the mixture.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight, parameters and covariance
 ///        from the components
 /// @param angleDesc The angle description object
 /// @param method How to reduce the mixture
 /// @return parameters and covariance
 template <typename component_range_t, typename projector_t,
           typename angle_desc_t>
 std::tuple<BoundVector, BoundMatrix> mergeGaussianMixture(
     const component_range_t &cmps, const projector_t &projector,
     const angle_desc_t &angleDesc, ComponentMergeMethod method) {
   const auto [mean, cov] =
       mergeGaussianMixtureMeanCov(cmps, projector, angleDesc);
 
   if (method == ComponentMergeMethod::eMean) {
     return {mean, cov};
   } else if (method == ComponentMergeMethod::eMaxWeight) {
     const BoundVector mode =
         mergeGaussianMixtureMaxWeight(cmps, [&](const auto &c) {
           const auto [weight_l, pars_l, cov_l] = projector(c);
           return std::tuple(weight_l, pars_l);
         });
     return {mode, cov};
   } else {
     throw std::logic_error("Invalid component merge method");
   }
 }
 
 /// Reduce Gaussian mixture into a single parameter vector, using the specified
 /// method to reduce the mixture.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight and parameters from the components
 /// @param surface The surface, which the bound state is on
 /// @param method How to reduce the mixture
 /// @return parameters
 template <typename component_range_t, typename projector_t>
 BoundVector mergeGaussianMixtureParams(const component_range_t &cmps,
                                        const projector_t &projector,
                                        const Surface &surface,
                                        ComponentMergeMethod method) {
   return angleDescriptionSwitch(surface, [&](const auto &desc) {
     return mergeGaussianMixtureParams(cmps, projector, desc, method);
   });
 }
 
 /// Computes the Kullback-Leibler distance between two components as shown in
 /// https://arxiv.org/abs/2001.00727v1 but ignoring the weights
 double computeSymmetricKlDivergence(const GsfComponent &a,
                                     const GsfComponent &b);
 
 /// Reduce Gaussian mixture into a single parameter vector and covariance, using
 /// the specified method to reduce the mixture.
 ///
 /// @param cmps The component range to merge
 /// @param projector A projector to extract the weight, parameters and covariance
 ///        from the components
 /// @param surface The surface, which the bound state is on
 /// @param method How to reduce the mixture
 /// @return parameters and covariance
 template <typename component_range_t, typename projector_t>
 std::tuple<BoundVector, BoundMatrix> mergeGaussianMixture(
     const component_range_t &cmps, const projector_t &projector,
     const Surface &surface, ComponentMergeMethod method) {
   return angleDescriptionSwitch(surface, [&](const auto &desc) {
     return mergeGaussianMixture(cmps, projector, desc, method);
   });
 }
 
 /// @brief Class representing a symmetric distance matrix
 ///
 /// Pairwise KL distances are kept in a flat lower-triangular array over the
 /// currently *active* components. Removing a component via
 /// @ref maskAssociatedDistances swaps its row/column with the last active
 /// one and shrinks the active count, so scans only touch present components.
 /// All public methods use original indices; the swap-based compaction is an
 /// internal detail.
 class SymmetricKLDistanceMatrix {
   using Array = Eigen::Array<double, Eigen::Dynamic, 1>;
 
   Array m_distances;
   std::vector<std::pair<std::size_t, std::size_t>> m_mapToPair;
   std::vector<std::size_t> m_activeToOriginal;
   std::vector<std::size_t> m_originalToActive;
   std::size_t m_numberComponents{};
   std::size_t m_numberActive{};
 
   static std::size_t triangularIndex(std::size_t i, std::size_t j) {
     return i * (i - 1) / 2 + j;
   }
 
   template <typename setter_t>
   void setAssociated(std::size_t n, const setter_t &setter) {
     // Rows
     for (std::size_t i = 0ul; i < n; ++i) {
       m_distances[triangularIndex(n, i)] = setter(n, i);
     }
     // Columns
     for (std::size_t i = n + 1; i < m_numberActive; ++i) {
       m_distances[triangularIndex(i, n)] = setter(i, n);
     }
   }
 
  public:
   explicit SymmetricKLDistanceMatrix(std::span<const GsfComponent> cmps);
 
   void recomputeAssociatedDistances(std::size_t originalIdx,
                                     std::span<const GsfComponent> cmps);
 
   void maskAssociatedDistances(std::size_t originalIdx);
 
   std::pair<std::size_t, std::size_t> minDistancePair() const;
 
   std::ostream &toStream(std::ostream &os) const;
 
   friend std::ostream &operator<<(std::ostream &os,
                                   const SymmetricKLDistanceMatrix &m) {
     return m.toStream(os);
   }
 };
 
 }  // namespace Acts::detail::Gsf

This page was automatically generated by the 2.3.7 LXR engine. The LXR team