Skip to content

Commit

Permalink
refactor!: Remove mean reduction from MultiEigenStepperLoop (acts-p…
Browse files Browse the repository at this point in the history 
…roject#3671)

After acts-project#3521 we also want to change the default and remove the weighted mean reduction as it is unstable.

This pull request primarily focuses on removing the `WeightedComponentReducerLoop` struct and its associated functionality from the `MultiEigenStepperLoop` class and its tests. The changes simplify the codebase by eliminating the weighted component reduction logic and replacing it with the `MaxWeightReducerLoop`.

### Removal of `WeightedComponentReducerLoop`:

* [`Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp`](diffhunk://#diff-2ec4ab79520c414c4eb526d00ede45a6a10f8eeba06398d81c315f9be5686e9eL48-L145): Removed the `WeightedComponentReducerLoop` struct and its methods, which were responsible for reducing the multicomponent state by summing weighted values.

### Replacement with `MaxWeightReducerLoop`:

* [`Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp`](diffhunk://#diff-2ec4ab79520c414c4eb526d00ede45a6a10f8eeba06398d81c315f9be5686e9eL243-R145): Updated the `MultiEigenStepperLoop` class to use `MaxWeightReducerLoop` instead of `WeightedComponentReducerLoop` as the default component reducer.

### Test updates:

* [`Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp`](diffhunk://#diff-35cd1eb9d226dec9f6017cd99b5605c2da3bbaa4c3f9ffbe5b0bed7efcfdc0ddL141-L165): Removed the test case for `WeightedComponentReducerLoop`, which validated the position and direction calculations using weighted sums.
  • Loading branch information
@andiwand
andiwand authored Oct 1, 2024
1 parent aa53675 commit dd33d4c
Show file tree
Hide file tree
Showing 3 changed files with 2 additions and 125 deletions.
100 changes: 1 addition & 99 deletions Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,104 +45,6 @@ namespace Acts {

using namespace Acts::UnitLiterals;

/// @brief Reducer struct for the Loop MultiEigenStepper which reduces the
/// multicomponent state to simply by summing the weighted values
///
/// @note Usage is not encouraged, since it can lead to navigation failures
/// as the global position might not be on surface, even if all components
/// are on surface
struct WeightedComponentReducerLoop {
template <typename component_range_t>
static Vector3 toVector3(const component_range_t& comps,
const FreeIndices i) {
return std::accumulate(
comps.begin(), comps.end(), Vector3{Vector3::Zero()},
[i](const auto& sum, const auto& cmp) -> Vector3 {
return sum + cmp.weight * cmp.state.pars.template segment<3>(i);
});
}

template <typename stepper_state_t>
static Vector3 position(const stepper_state_t& s) {
return toVector3(s.components, eFreePos0);
}

template <typename stepper_state_t>
static Vector3 direction(const stepper_state_t& s) {
return toVector3(s.components, eFreeDir0).normalized();
}

// TODO: Maybe we can cache this value and only update it when the parameters
// change
template <typename stepper_state_t>
static ActsScalar qOverP(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * cmp.state.pars[eFreeQOverP];
});
}

template <typename stepper_state_t>
static ActsScalar absoluteMomentum(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[&s](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * s.particleHypothesis.extractMomentum(
cmp.state.pars[eFreeQOverP]);
});
}

template <typename stepper_state_t>
static Vector3 momentum(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), Vector3::Zero().eval(),
[&s](const auto& sum, const auto& cmp) -> Vector3 {
return sum + cmp.weight *
s.particleHypothesis.extractMomentum(
cmp.state.pars[eFreeQOverP]) *
cmp.state.pars.template segment<3>(eFreeDir0);
});
}

template <typename stepper_state_t>
static ActsScalar charge(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[&s](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * s.particleHypothesis.extractCharge(
cmp.state.pars[eFreeQOverP]);
});
}

template <typename stepper_state_t>
static ActsScalar time(const stepper_state_t& s) {
return std::accumulate(
s.components.begin(), s.components.end(), ActsScalar{0.},
[](const auto& sum, const auto& cmp) -> ActsScalar {
return sum + cmp.weight * cmp.state.pars[eFreeTime];
});
}

template <typename stepper_state_t>
static FreeVector pars(const stepper_state_t& s) {
return std::accumulate(s.components.begin(), s.components.end(),
FreeVector{FreeVector::Zero()},
[](const auto& sum, const auto& cmp) -> FreeVector {
return sum + cmp.weight * cmp.state.pars;
});
}

template <typename stepper_state_t>
static FreeVector cov(const stepper_state_t& s) {
return std::accumulate(s.components.begin(), s.components.end(),
FreeMatrix{FreeMatrix::Zero()},
[](const auto& sum, const auto& cmp) -> FreeMatrix {
return sum + cmp.weight * cmp.state.cov;
});
}
};

namespace detail {

struct MaxMomentumComponent {
Expand Down Expand Up @@ -240,7 +142,7 @@ using MaxWeightReducerLoop =
/// @tparam small_vector_size A size-hint how much memory should be allocated
/// by the small vector
template <typename extension_t = EigenStepperDefaultExtension,
typename component_reducer_t = WeightedComponentReducerLoop>
typename component_reducer_t = MaxWeightReducerLoop>
class MultiEigenStepperLoop : public EigenStepper<extension_t> {
/// Limits the number of steps after at least one component reached the
/// surface
Expand Down
25 changes: 0 additions & 25 deletions Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,31 +138,6 @@ auto makeDefaultBoundPars(bool cov = true, std::size_t n = 4,
//////////////////////
/// Test the reducers
//////////////////////
BOOST_AUTO_TEST_CASE(test_weighted_reducer) {
// Can use this multistepper since we only care about the state which is
// invariant
using MultiState = typename MultiStepperLoop::State;

constexpr std::size_t N = 4;
const auto multi_pars = makeDefaultBoundPars(false, N);

MultiState state(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize);
SingleStepper singleStepper(defaultBField);

WeightedComponentReducerLoop reducer{};

Acts::Vector3 pos = Acts::Vector3::Zero();
Acts::Vector3 dir = Acts::Vector3::Zero();
for (const auto &[sstate, weight, _] : state.components) {
pos += weight * singleStepper.position(sstate);
dir += weight * singleStepper.direction(sstate);
}
dir.normalize();

BOOST_CHECK_EQUAL(reducer.position(state), pos);
BOOST_CHECK_EQUAL(reducer.direction(state), dir);
}

BOOST_AUTO_TEST_CASE(test_max_weight_reducer) {
// Can use this multistepper since we only care about the state which is
// invariant
Expand Down
2 changes: 1 addition & 1 deletion docs/core/reconstruction/track_fitting.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ Simplified overview of the GSF algorithm.
:::

### The Multi-Stepper
To implement the GSF, a special stepper is needed, that can handle a multi-component state internally: The {class}`Acts::MultiEigenStepperLoop`, which is based on the {class}`Acts::EigenStepper` and thus shares a lot of code with it. It interfaces to the navigation as one aggregate state to limit the navigation overhead, but internally processes a multi-component state. How this aggregation is performed can be configured via a template parameter, by default weighted average is used ({struct}`Acts::WeightedComponentReducerLoop`).
To implement the GSF, a special stepper is needed, that can handle a multi-component state internally: The {class}`Acts::MultiEigenStepperLoop`, which is based on the {class}`Acts::EigenStepper` and thus shares a lot of code with it. It interfaces to the navigation as one aggregate state to limit the navigation overhead, but internally processes a multi-component state. How this aggregation is performed can be configured via a template parameter, by default maximum weight is used ({type}`Acts::MaxWeightReducerLoop`).

Even though the multi-stepper interface exposes only one aggregate state and thus is compatible with most standard tools, there is a special aborter is required to stop the navigation when the surface is reached, the {struct}`Acts::MultiStepperSurfaceReached`. It checks if all components have reached the target surface already and updates their state accordingly. Optionally, it also can stop the propagation when the aggregate state reaches the surface.

Expand Down

0 comments on commit dd33d4c

Please sign in to comment.