Skip to content

Commit

Permalink
refactor(gx2f): put parts into separate compile units (#3946)
Browse files Browse the repository at this point in the history 
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.

<!-- This is an auto-generated comment: release notes by coderabbit.ai -->
## 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.
<!-- end of auto-generated comment: release notes by coderabbit.ai -->
  • Loading branch information
@AJPfleger
AJPfleger authored Dec 6, 2024
1 parent f57e260 commit feac686
Show file tree
Hide file tree
Showing 2 changed files with 133 additions and 90 deletions.
128 changes: 38 additions & 90 deletions Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<BoundMatrix>& 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.
Expand All @@ -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
Expand All @@ -379,44 +402,9 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem,
const std::vector<BoundMatrix>& 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<kMeasDim> covarianceMeasurement =
trackState.template calibratedCovariance<kMeasDim>();

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<eBoundSize, eBoundSize>() =
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<eBoundSize, 2> 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<eBoundSize, 2>(0, deltaPosition) = jacPhiTheta;
}

const BoundVector predicted = trackState.smoothed();

const ActsVector<kMeasDim> measurement =
Expand All @@ -425,54 +413,9 @@ void addMeasurementToGx2fSums(Gx2fSystem& extendedSystem,
const ActsMatrix<kMeasDim, eBoundSize> projector =
trackState.template projectorSubspaceHelper<kMeasDim>().projector();

const Eigen::MatrixXd projJacobian = projector * extendedJacobian;

const ActsMatrix<kMeasDim, 1> projPredicted = projector * predicted;

const ActsVector<kMeasDim> 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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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"
Expand Down
95 changes: 95 additions & 0 deletions Core/src/TrackFitting/GlobalChiSquareFitter.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,3 +50,98 @@ void Acts::Experimental::updateGx2fCovarianceParams(

return;
}

void Acts::Experimental::addMeasurementToGx2fSumsBackend(
Gx2fSystem& extendedSystem,
const std::vector<BoundMatrix>& 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<eBoundSize, eBoundSize>() =
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<eBoundSize, 2> 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<eBoundSize, 2>(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());
}

0 comments on commit feac686

Please sign in to comment.