From feac68618ac7da66fe811ca42cfe00187bdcacf5 Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Fri, 6 Dec 2024 03:44:46 +0100 Subject: [PATCH] refactor(gx2f): put parts into separate compile units (#3946) This reduces the memory impact of the GX2F during compilation by: - extracting heavy parts to a separate compile unit - removing the templating from the extracted parts (to make them moveable in the first place) ## Impact main: ``` [ 6.03M, max: 1072.02M] [ 9.74s] - Core/src/TrackFitting/GlobalChiSquareFitter.cpp [ 7.05M, max: 3364.18M] [ 52.19s] - Examples/Algorithms/TrackFitting/src/GlobalChiSquareFitterFunction.cpp [ 6.93M, max: 4196.06M] [ 85.36s] - Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp ``` after addMeasurementToGx2fSumsBackend: ``` [ 6.04M, max: 1408.93M] [ 16.37s] - Core/src/TrackFitting/GlobalChiSquareFitter.cpp [ 6.91M, max: 2299.42M] [ 29.93s] - Examples/Algorithms/TrackFitting/src/GlobalChiSquareFitterFunction.cpp [ 6.94M, max: 3269.05M] [ 59.98s] - Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp ``` after computeGx2fDeltaParams: ``` [ 6.40M, max: 1769.83M] [ 20.65s] - Core/src/TrackFitting/GlobalChiSquareFitter.cpp [ 5.94M, max: 1805.96M] [ 24.37s] - Examples/Algorithms/TrackFitting/src/GlobalChiSquareFitterFunction.cpp [ 7.23M, max: 2510.50M] [ 53.56s] - Tests/UnitTests/Core/TrackFitting/Gx2fTests.cpp ``` ## Other discussion on the topic ### Performance of dynamic matrix This still needs investigation. Maybe we can template the functions again on `kMeasDims` ### Untemplate temporary track This was the first attempt, by using `VectorMultiTrajectory` for the internal tracks. This fails in Athena, because there the `Mutable...` is used. This comes in with the calibrator. Currently, it might work with 2 calibrators (one for the fitting Actor and one for the final Actor). Better solution would be to have type-erasure, which might come in the future. ## Summary by CodeRabbit - **New Features** - Introduced a new method for adding measurement data into the fitting system, enhancing the integration process. - Added a function to compute delta parameters for improved clarity in the fitting process. - **Improvements** - Streamlined measurement handling and error management within the fitting framework for better efficiency. - Enhanced the modularity of the fitting process, improving maintainability and clarity. --- .../TrackFitting/GlobalChiSquareFitter.hpp | 128 ++++++------------ .../TrackFitting/GlobalChiSquareFitter.cpp | 95 +++++++++++++ 2 files changed, 133 insertions(+), 90 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index dc0bc221e9a..5e46bc0dba3 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -361,6 +361,29 @@ struct Gx2fSystem { std::size_t m_ndf = 0u; }; +/// @brief Adds a measurement to the GX2F equation system in a modular backend function. +/// +/// This function processes measurement data and integrates it into the GX2F +/// system. +/// +/// @param extendedSystem All parameters of the current equation system to update. +/// @param jacobianFromStart The Jacobian matrix from the start to the current state. +/// @param covarianceMeasurement The covariance matrix of the measurement. +/// @param predicted The predicted state vector based on the track state. +/// @param measurement The measurement vector. +/// @param projector The projection matrix. +/// @param logger A logger instance. +/// +/// @note The dynamic Eigen matrices are suboptimal. We could think of +/// templating again in the future on kMeasDims. We currently use dynamic +/// matrices to reduce the memory during compile time. +void addMeasurementToGx2fSumsBackend( + Gx2fSystem& extendedSystem, + const std::vector& jacobianFromStart, + const Eigen::MatrixXd& covarianceMeasurement, const BoundVector& predicted, + const Eigen::VectorXd& measurement, const Eigen::MatrixXd& projector, + const Logger& logger); + /// @brief Process measurements and fill the aMatrix and bVector /// /// The function processes each measurement for the GX2F Actor fitting process. @@ -370,7 +393,7 @@ struct Gx2fSystem { /// @tparam kMeasDim Number of dimensions of the measurement /// @tparam track_state_t The type of the track state /// -/// @param extendedSystem All parameters of the current equation system +/// @param extendedSystem All parameters of the current equation system to update /// @param jacobianFromStart The Jacobian matrix from start to the current state /// @param trackState The track state to analyse /// @param logger A logger instance @@ -379,44 +402,9 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem, const std::vector& jacobianFromStart, const track_state_t& trackState, const Logger& logger) { - // First we get back the covariance and try to invert it. If the inversion - // fails, we can already abort. const ActsSquareMatrix covarianceMeasurement = trackState.template calibratedCovariance(); - const auto safeInvCovMeasurement = safeInverse(covarianceMeasurement); - if (!safeInvCovMeasurement) { - ACTS_WARNING("addMeasurementToGx2fSums: safeInvCovMeasurement failed."); - ACTS_VERBOSE(" covarianceMeasurement:\n" << covarianceMeasurement); - return; - } - - // Create an extended Jacobian. This one contains only eBoundSize rows, - // because the rest is irrelevant. We fill it in the next steps. - // TODO make dimsExtendedParams template with unrolling - Eigen::MatrixXd extendedJacobian = - Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims()); - - // This part of the Jacobian comes from the material-less propagation - extendedJacobian.topLeftCorner() = - jacobianFromStart[0]; - - // If we have material, loop here over all Jacobians. We add extra columns for - // their phi-theta projections. These parts account for the propagation of the - // scattering angles. - for (std::size_t matSurface = 1; matSurface < jacobianFromStart.size(); - matSurface++) { - const BoundMatrix jac = jacobianFromStart[matSurface]; - - const ActsMatrix jacPhiTheta = - jac * Gx2fConstants::phiThetaProjector; - - // The position, where we need to insert the values in the extended Jacobian - const std::size_t deltaPosition = eBoundSize + 2 * (matSurface - 1); - - extendedJacobian.block(0, deltaPosition) = jacPhiTheta; - } - const BoundVector predicted = trackState.smoothed(); const ActsVector measurement = @@ -425,54 +413,9 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem, const ActsMatrix projector = trackState.template projectorSubspaceHelper().projector(); - const Eigen::MatrixXd projJacobian = projector * extendedJacobian; - - const ActsMatrix projPredicted = projector * predicted; - - const ActsVector residual = measurement - projPredicted; - - // Finally contribute to chi2sum, aMatrix, and bVector - extendedSystem.chi2() += - (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); - - extendedSystem.aMatrix() += - (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval(); - - extendedSystem.bVector() += - (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval() - .transpose(); - - ACTS_VERBOSE( - "Contributions in addMeasurementToGx2fSums:\n" - << " kMeasDim: " << kMeasDim << "\n" - << " predicted: " << predicted.transpose() << "\n" - << " measurement: " << measurement.transpose() << "\n" - << " covarianceMeasurement:\n" - << covarianceMeasurement << "\n" - << " projector:\n" - << projector.eval() << "\n" - << " projJacobian:\n" - << projJacobian.eval() << "\n" - << " projPredicted: " << (projPredicted.transpose()).eval() << "\n" - << " residual: " << (residual.transpose()).eval() << "\n" - << " extendedJacobian:\n" - << extendedJacobian << "\n" - << " aMatrix contribution:\n" - << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval() - << "\n" - << " bVector contribution: " - << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval() - << "\n" - << " chi2sum contribution: " - << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0) - << "\n" - << " safeInvCovMeasurement:\n" - << (*safeInvCovMeasurement)); - - return; + addMeasurementToGx2fSumsBackend(extendedSystem, jacobianFromStart, + covarianceMeasurement, predicted, measurement, + projector, logger); } /// @brief Process material and fill the aMatrix and bVector @@ -694,6 +637,15 @@ std::size_t countMaterialStates( return nMaterialSurfaces; } +/// @brief Solve the gx2f system to get the delta parameters for the update +/// +/// This function computes the delta parameters for the GX2F Actor fitting +/// process by solving the linear equation system [a] * delta = b. It uses the +/// column-pivoting Householder QR decomposition for numerical stability. +/// +/// @param extendedSystem All parameters of the current equation system +Eigen::VectorXd computeGx2fDeltaParams(const Gx2fSystem& extendedSystem); + /// @brief Update parameters (and scattering angles if applicable) /// /// @param params Parameters to be updated @@ -1391,10 +1343,8 @@ class Gx2Fitter { return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements; } - // calculate delta params [a] * delta = b Eigen::VectorXd deltaParamsExtended = - extendedSystem.aMatrix().colPivHouseholderQr().solve( - extendedSystem.bVector()); + computeGx2fDeltaParams(extendedSystem); ACTS_VERBOSE("aMatrix:\n" << extendedSystem.aMatrix() << "\n" @@ -1558,10 +1508,8 @@ class Gx2Fitter { return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements; } - // calculate delta params [a] * delta = b Eigen::VectorXd deltaParamsExtended = - extendedSystem.aMatrix().colPivHouseholderQr().solve( - extendedSystem.bVector()); + computeGx2fDeltaParams(extendedSystem); ACTS_VERBOSE("aMatrix:\n" << extendedSystem.aMatrix() << "\n" diff --git a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp index 50b8b4cab8c..4062124496e 100644 --- a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp +++ b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp @@ -50,3 +50,98 @@ void Acts::Experimental::updateGx2fCovarianceParams( return; } + +void Acts::Experimental::addMeasurementToGx2fSumsBackend( + Gx2fSystem& extendedSystem, + const std::vector& jacobianFromStart, + const Eigen::MatrixXd& covarianceMeasurement, const BoundVector& predicted, + const Eigen::VectorXd& measurement, const Eigen::MatrixXd& projector, + const Logger& logger) { + // First, w try to invert the covariance matrix. If the inversion fails, we + // can already abort. + const auto safeInvCovMeasurement = safeInverse(covarianceMeasurement); + if (!safeInvCovMeasurement) { + ACTS_WARNING("addMeasurementToGx2fSums: safeInvCovMeasurement failed."); + ACTS_VERBOSE(" covarianceMeasurement:\n" << covarianceMeasurement); + return; + } + + // Create an extended Jacobian. This one contains only eBoundSize rows, + // because the rest is irrelevant. We fill it in the next steps. + // TODO make dimsExtendedParams template with unrolling + Eigen::MatrixXd extendedJacobian = + Eigen::MatrixXd::Zero(eBoundSize, extendedSystem.nDims()); + + // This part of the Jacobian comes from the material-less propagation + extendedJacobian.topLeftCorner() = + jacobianFromStart[0]; + + // If we have material, loop here over all Jacobians. We add extra columns for + // their phi-theta projections. These parts account for the propagation of the + // scattering angles. + for (std::size_t matSurface = 1; matSurface < jacobianFromStart.size(); + matSurface++) { + const BoundMatrix jac = jacobianFromStart[matSurface]; + + const ActsMatrix jacPhiTheta = + jac * Gx2fConstants::phiThetaProjector; + + // The position, where we need to insert the values in the extended Jacobian + const std::size_t deltaPosition = eBoundSize + 2 * (matSurface - 1); + + extendedJacobian.template block(0, deltaPosition) = + jacPhiTheta; + } + + const Eigen::MatrixXd projJacobian = projector * extendedJacobian; + + const Eigen::VectorXd projPredicted = projector * predicted; + + const Eigen::VectorXd residual = measurement - projPredicted; + + // Finally contribute to chi2sum, aMatrix, and bVector + extendedSystem.chi2() += + (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); + + extendedSystem.aMatrix() += + (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval(); + + extendedSystem.bVector() += + (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval() + .transpose(); + + ACTS_VERBOSE( + "Contributions in addMeasurementToGx2fSums:\n" + << " predicted: " << predicted.transpose() << "\n" + << " measurement: " << measurement.transpose() << "\n" + << " covarianceMeasurement:\n" + << covarianceMeasurement << "\n" + << " projector:\n" + << projector.eval() << "\n" + << " projJacobian:\n" + << projJacobian.eval() << "\n" + << " projPredicted: " << (projPredicted.transpose()).eval() << "\n" + << " residual: " << (residual.transpose()).eval() << "\n" + << " extendedJacobian:\n" + << extendedJacobian << "\n" + << " aMatrix contribution:\n" + << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval() + << "\n" + << " bVector contribution: " + << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval() + << "\n" + << " chi2sum contribution: " + << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0) + << "\n" + << " safeInvCovMeasurement:\n" + << (*safeInvCovMeasurement)); +} + +Eigen::VectorXd Acts::Experimental::computeGx2fDeltaParams( + const Acts::Experimental::Gx2fSystem& extendedSystem) { + return extendedSystem.aMatrix().colPivHouseholderQr().solve( + extendedSystem.bVector()); +}