Skip to content

Commit

Permalink
fix: fix FPE and refactor GSF component reduction (#1710)
Browse files Browse the repository at this point in the history 
Refactors the KL mixture reduction and fix FPEs occuring there. Also improves the unit testing of these features.
  • Loading branch information
@benjaminhuth
benjaminhuth committed Dec 9, 2022
1 parent 440f092 commit 6045480
Show file tree
Hide file tree
Showing 4 changed files with 245 additions and 124 deletions.
134 changes: 79 additions & 55 deletions Core/include/Acts/TrackFitting/detail/KLMixtureReduction.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,10 @@ namespace Acts {
namespace detail {

/// Computes the Kullback-Leibler distance between two components as shown in
/// https://arxiv.org/abs/2001.00727v1, while ignoring the weights as done in
/// Athena
/// https://arxiv.org/abs/2001.00727v1 but ignoring the weights
template <typename component_t, typename component_projector_t>
auto computeKLDistance(const component_t &a, const component_t &b,
const component_projector_t &proj) {
auto computeSymmetricKlDivergence(const component_t &a, const component_t &b,
const component_projector_t &proj) {
const auto parsA = proj(a).boundPars[eBoundQOverP];
const auto parsB = proj(b).boundPars[eBoundQOverP];
const auto covA = (*proj(a).boundCov)(eBoundQOverP, eBoundQOverP);
Expand Down Expand Up @@ -70,83 +69,111 @@ auto mergeComponents(const component_t &a, const component_t &b,

/// @brief Class representing a symmetric distance matrix
class SymmetricKLDistanceMatrix {
Eigen::VectorXd m_data;
using Array = Eigen::Array<ActsScalar, Eigen::Dynamic, 1>;
using Mask = Eigen::Array<bool, Eigen::Dynamic, 1>;

Array m_distances;
Mask m_mask;
std::vector<std::pair<std::size_t, std::size_t>> m_mapToPair;
std::size_t m_N;
std::size_t m_numberComponents;

template <typename array_t, typename setter_t>
void setAssociated(std::size_t n, array_t &array, setter_t &&setter) {
const auto indexConst = (n - 1) * n / 2;

// Rows
for (auto i = 0ul; i < n; ++i) {
array[indexConst + i] = setter(n, i);
}

// Columns
for (auto i = n + 1; i < m_numberComponents; ++i) {
array[(i - 1) * i / 2 + n] = setter(n, i);
}
}

public:
template <typename component_t, typename projector_t>
SymmetricKLDistanceMatrix(const std::vector<component_t> &cmps,
const projector_t &proj)
: m_data(cmps.size() * (cmps.size() - 1) / 2),
m_mapToPair(m_data.size()),
m_N(cmps.size()) {
for (auto i = 1ul; i < m_N; ++i) {
: m_distances(Array::Zero(cmps.size() * (cmps.size() - 1) / 2)),
m_mask(Mask::Ones(cmps.size() * (cmps.size() - 1) / 2)),
m_mapToPair(m_distances.size()),
m_numberComponents(cmps.size()) {
for (auto i = 1ul; i < m_numberComponents; ++i) {
const auto indexConst = (i - 1) * i / 2;
for (auto j = 0ul; j < i; ++j) {
m_mapToPair.at(indexConst + j) = {i, j};
m_data[indexConst + j] = computeKLDistance(cmps[i], cmps[j], proj);
m_distances[indexConst + j] =
computeSymmetricKlDivergence(cmps[i], cmps[j], proj);
}
}
}

auto at(std::size_t i, std::size_t j) const {
return m_data[i * (i - 1) / 2 + j];
return m_distances[i * (i - 1) / 2 + j];
}

template <typename component_t, typename projector_t>
void recomputeAssociatedDistances(std::size_t n,
const std::vector<component_t> &cmps,
const projector_t &proj) {
const auto indexConst = (n - 1) * n / 2;
assert(cmps.size() == m_numberComponents && "size mismatch");

throw_assert(cmps.size() == m_N, "size mismatch");

// Rows
for (auto i = 0ul; i < n; ++i) {
m_data[indexConst + i] = computeKLDistance(cmps[n], cmps[i], proj);
}

// Columns
for (auto i = n + 1; i < cmps.size(); ++i) {
m_data[(i - 1) * i / 2 + n] = computeKLDistance(cmps[n], cmps[i], proj);
}
setAssociated(n, m_distances, [&](std::size_t i, std::size_t j) {
return computeSymmetricKlDivergence(cmps[i], cmps[j], proj);
});
}

void resetAssociatedDistances(std::size_t n, double value) {
const auto indexConst = (n - 1) * n / 2;

// Rows
for (auto i = 0ul; i < n; ++i) {
m_data[indexConst + i] = value;
}
void maskAssociatedDistances(std::size_t n) {
setAssociated(n, m_mask, [&](std::size_t, std::size_t) { return false; });
}

// Columns
for (auto i = n + 1; i < m_N; ++i) {
m_data[(i - 1) * i / 2 + n] = value;
}
void setAssociatedDistances(std::size_t n, ActsScalar value) {
setAssociated(n, m_mask, [&](std::size_t, std::size_t) { return value; });
}

auto minDistancePair() const {
// TODO Eigen minCoeff does not work for some reason??
// return m_mapToPair.at(m_data.minCoeff());
return m_mapToPair.at(std::distance(
&m_data[0],
std::min_element(m_data.data(), m_data.data() + m_data.size())));
ActsScalar min = std::numeric_limits<ActsScalar>::max();
std::size_t idx = 0;

for (auto i = 0l; i < m_distances.size(); ++i) {
if (auto new_min = std::min(min, m_distances[i]);
m_mask[i] && new_min < min) {
min = new_min;
idx = i;
}
}

return m_mapToPair.at(idx);
}

friend std::ostream &operator<<(std::ostream &os,
const SymmetricKLDistanceMatrix &m) {
for (auto i = 1ul; i < m.m_N; ++i) {
const auto prev_precision = os.precision();
const int width = 8;
const int prec = 2;

os << "\n";
os << std::string(width, ' ') << " | ";
for (auto j = 0ul; j < m.m_numberComponents - 1; ++j) {
os << std::setw(width) << j << " ";
}
os << "\n";
os << std::string((width + 3) + (width + 2) * (m.m_numberComponents - 1),
'-');
os << "\n";

for (auto i = 1ul; i < m.m_numberComponents; ++i) {
const auto indexConst = (i - 1) * i / 2;
Eigen::RowVectorXd vals;
vals.resize(i);
os << std::setw(width) << i << " | ";
for (auto j = 0ul; j < i; ++j) {
vals[j] = m.m_data[indexConst + j];
os << std::setw(width) << std::setprecision(prec)
<< m.m_distances[indexConst + j] << " ";
}
os << vals << "\n";
os << "\n";
}

os << std::setprecision(prev_precision);
return os;
}
};
Expand All @@ -168,20 +195,17 @@ void reduceWithKLDistance(std::vector<component_t> &cmpCache,
while (remainingComponents > maxCmpsAfterMerge) {
const auto [minI, minJ] = distances.minDistancePair();

// Set one component and compute associated distances
cmpCache[minI] =
mergeComponents(cmpCache[minI], cmpCache[minJ], proj, angle_desc);
distances.recomputeAssociatedDistances(minI, cmpCache, proj);
remainingComponents--;

// Reset removed components so that it won't have the shortest distance
// ever, and so that we can sort them by weight in the end to remove them
// Set weight of the other component to -1 so we can remove it later and
// mask its distances
proj(cmpCache[minJ]).weight = -1.0;
proj(cmpCache[minJ]).boundPars[eBoundQOverP] =
std::numeric_limits<double>::max();
(*proj(cmpCache[minJ]).boundCov)(eBoundQOverP, eBoundQOverP) =
std::numeric_limits<double>::max();
distances.resetAssociatedDistances(minJ,
std::numeric_limits<double>::max());
distances.maskAssociatedDistances(minJ);

remainingComponents--;
}

// Remove all components which are labled with weight -1
Expand Down
1 change: 1 addition & 0 deletions Tests/UnitTests/Core/TrackFitting/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,5 @@ add_unittest(GainMatrixUpdater GainMatrixUpdaterTests.cpp)
add_unittest(KalmanFitter KalmanFitterTests.cpp)
add_unittest(Gsf GsfTests.cpp)
add_unittest(GsfComponentMerging GsfComponentMergingTests.cpp)
add_unittest(GsfMixtureReduction GsfMixtureReductionTests.cpp)
add_unittest(Chi2Fitter Chi2FitterTests.cpp)
71 changes: 2 additions & 69 deletions Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,13 +8,14 @@

#include <boost/test/unit_test.hpp>

#include "Acts/Definitions/Units.hpp"
#include "Acts/EventData/detail/TransformationBoundToFree.hpp"
#include "Acts/EventData/detail/TransformationFreeToBound.hpp"
#include "Acts/Surfaces/CylinderSurface.hpp"
#include "Acts/Surfaces/DiscSurface.hpp"
#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Surfaces/PlaneSurface.hpp"
#include "Acts/TrackFitting/detail/KLMixtureReduction.hpp"
#include "Acts/Utilities/detail/gaussian_mixture_helpers.hpp"

#include <random>

Expand Down Expand Up @@ -342,71 +343,3 @@ BOOST_AUTO_TEST_CASE(test_perigee_surface) {
// Here we expect a very bad approximation
test_surface(*surface, desc, p, 1.);
}

BOOST_AUTO_TEST_CASE(test_kl_mixture_reduction) {
auto meanAndSumOfWeights = [](const auto &cmps) {
const auto mean = std::accumulate(
cmps.begin(), cmps.end(), Acts::BoundVector::Zero().eval(),
[](auto sum, const auto &cmp) -> Acts::BoundVector {
return sum + cmp.weight * cmp.boundPars;
});

const double sumOfWeights = std::accumulate(
cmps.begin(), cmps.end(), 0.0,
[](auto sum, const auto &cmp) { return sum + cmp.weight; });

return std::make_tuple(mean, sumOfWeights);
};

// Do not bother with circular angles in this test
const auto desc = std::tuple<>{};

// Need this projection, since we need to write to the lvalue references which
// isn't possible through Identity / std::identity due to perfect forwarding
const auto proj = [](auto &a) -> decltype(auto) { return a; };

const std::size_t NComps = 4;
std::vector<DummyComponent<eBoundSize>> cmps;

for (auto i = 0ul; i < NComps; ++i) {
DummyComponent<eBoundSize> a;
a.boundPars = Acts::BoundVector::Zero();
a.boundCov = Acts::BoundSymMatrix::Identity();
a.weight = 1.0 / NComps;
cmps.push_back(a);
}

cmps[0].boundPars[eBoundQOverP] = 0.5_GeV;
cmps[1].boundPars[eBoundQOverP] = 1.5_GeV;
cmps[2].boundPars[eBoundQOverP] = 3.5_GeV;
cmps[3].boundPars[eBoundQOverP] = 4.5_GeV;

// Check start properties
const auto [mean0, sumOfWeights0] = meanAndSumOfWeights(cmps);

BOOST_CHECK_CLOSE(mean0[eBoundQOverP], 2.5_GeV, 1.e-8);
BOOST_CHECK_CLOSE(sumOfWeights0, 1.0, 1.e-8);

// Reduce by factor of 2 and check if weights and QoP are correct
Acts::detail::reduceWithKLDistance(cmps, 2, proj, desc);

BOOST_CHECK(cmps.size() == 2);

std::sort(cmps.begin(), cmps.end(), [](const auto &a, const auto &b) {
return a.boundPars[eBoundQOverP] < b.boundPars[eBoundQOverP];
});
BOOST_CHECK_CLOSE(cmps[0].boundPars[eBoundQOverP], 1.0_GeV, 1.e-8);
BOOST_CHECK_CLOSE(cmps[1].boundPars[eBoundQOverP], 4.0_GeV, 1.e-8);

const auto [mean1, sumOfWeights1] = meanAndSumOfWeights(cmps);

BOOST_CHECK_CLOSE(mean1[eBoundQOverP], 2.5_GeV, 1.e-8);
BOOST_CHECK_CLOSE(sumOfWeights1, 1.0, 1.e-8);

// Reduce by factor of 2 and check if weights and QoP are correct
Acts::detail::reduceWithKLDistance(cmps, 1, proj, desc);

BOOST_CHECK(cmps.size() == 1);
BOOST_CHECK_CLOSE(cmps[0].boundPars[eBoundQOverP], 2.5_GeV, 1.e-8);
BOOST_CHECK_CLOSE(cmps[0].weight, 1.0, 1.e-8);
}

0 comments on commit 6045480

Please sign in to comment.