Core/Utilities/AlgebraHelpersTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Utilities/AlgebraHelpers.hpp"
0012 #include "Acts/Utilities/Logger.hpp"
0013
0014 #include <Eigen/Dense>
0015
0016 Acts::Logging::Level logLevel = Acts::Logging::VERBOSE;
0017
0018 namespace Acts::Test {
0019
0020 BOOST_AUTO_TEST_SUITE(AlgebraHelpers)
0021
0022 BOOST_AUTO_TEST_SUITE(SafeInverse)
0023
0024 ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("SafeInverse", logLevel))
0025
0026 BOOST_AUTO_TEST_CASE(SafeInverseSmallMatrix) {
0027   Eigen::Matrix<double, 2, 2> m;
0028   m << 1, 2, 3, 4;
0029
0030   Eigen::Matrix<double, 2, 2> mInvRef;
0031   mInvRef << -2, 1, 1.5, -0.5;
0032
0033   auto mInv = Acts::safeInverse(m);
0034
0035   BOOST_CHECK(mInv);
0036   BOOST_CHECK_EQUAL(*mInv, mInvRef);
0037
0038   Eigen::Matrix<int, 2, 2> identityInt;
0039   identityInt << 1, 0, 0, 1;
0040   auto identityIntInv = Acts::safeInverse(identityInt);
0041
0042   BOOST_CHECK(identityIntInv);
0043   BOOST_CHECK_EQUAL(*identityIntInv, identityInt);
0044 }
0045
0046 BOOST_AUTO_TEST_CASE(safeInverseLargeMatrix) {
0047   const Eigen::Matrix<double, 5, 5> identity{Eigen::MatrixXd::Identity(5, 5)};
0048   auto identityInv = Acts::safeInverse(identity);
0049
0050   BOOST_CHECK(identityInv);
0051   BOOST_CHECK_EQUAL(*identityInv, identity);
0052 }
0053
0054 BOOST_AUTO_TEST_CASE(safeInverseDynamicMatrix) {
0055   Eigen::MatrixXd identity{Eigen::MatrixXd::Identity(2, 2)};
0056
0057   auto identityInv = Acts::safeInverse(identity);
0058
0059   BOOST_CHECK(identityInv);
0060   BOOST_CHECK_EQUAL(*identityInv, identity);
0061 }
0062
0063 BOOST_AUTO_TEST_CASE(SafeInverseBadSmallMatrix) {
0064   Eigen::Matrix<double, 2, 2> m;
0065   m << 1, 1, 2, 2;
0066
0067   auto mInv = Acts::safeInverse(m);
0068
0069   BOOST_CHECK(!mInv);
0070 }
0071
0072 BOOST_AUTO_TEST_CASE(safeInverseBadLargeMatrix) {
0073   const Eigen::Matrix<double, 5, 5> m{Eigen::MatrixXd::Zero(5, 5)};
0074   auto mInv = Acts::safeInverse(m);
0075
0076   BOOST_CHECK(!mInv);
0077 }
0078
0079 BOOST_AUTO_TEST_CASE(SafeInverseFPESmallMatrix) {
0080   Eigen::Matrix<double, 4, 4> m =
0081       Eigen::MatrixXd::Identity(4, 4) * std::numeric_limits<std::size_t>::max();
0082   m(1, 1) = 1;
0083
0084   auto mInv = Acts::safeInverse(m);
0085   BOOST_REQUIRE(mInv.has_value());
0086   auto mInvInv = Acts::safeInverse(*mInv);
0087   BOOST_CHECK(!mInvInv);
0088
0089   ACTS_VERBOSE("Test: SafeInverseFPESmallMatrix" << "\n"
0090                                                  << "m:\n"
0091                                                  << m << "\n"
0092                                                  << "mInv:\n"
0093                                                  << *mInv);
0094 }
0095
0096 BOOST_AUTO_TEST_CASE(SafeInverseFPELargeMatrix) {
0097   Eigen::Matrix<double, 5, 5> m =
0098       Eigen::MatrixXd::Identity(5, 5) * std::numeric_limits<std::size_t>::max();
0099   m(1, 1) = 1;
0100
0101   auto mInv = Acts::safeInverse(m);
0102
0103   BOOST_CHECK(!mInv);
0104
0105   ACTS_VERBOSE("Test: SafeInverseFPELargeMatrix" << "\n"
0106                                                  << "m:\n"
0107                                                  << m);
0108 }
0109
0110 /// This test should not compile
0111 // BOOST_AUTO_TEST_CASE(SafeInverseNonsquareMatrix) {
0112 //   Eigen::Matrix<double, 2, 3> m;
0113 //   m << 1, 2, 3, 4, 5, 6;
0114 //
0115 //   auto mInv = Acts::safeInverse(m);
0116 // }
0117
0118 BOOST_AUTO_TEST_SUITE_END()
0119
0120 BOOST_AUTO_TEST_SUITE(SafeExp)
0121
0122 ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("SafeExp", logLevel))
0123
0124 BOOST_AUTO_TEST_CASE(safeExpDouble) {
0125   using FloatType = double;
0126
0127   // Values within the safe range
0128   BOOST_CHECK_CLOSE(safeExp<FloatType>(0.0), std::exp(0.0), 1e-8);
0129   BOOST_CHECK_CLOSE(safeExp<FloatType>(1.0), std::exp(1.0), 1e-8);
0130   BOOST_CHECK_CLOSE(safeExp<FloatType>(-1.0), std::exp(-1.0), 1e-8);
0131
0132   // Values causing underflow
0133   BOOST_CHECK_EQUAL(safeExp<FloatType>(-600.0), 0.0);
0134
0135   // Values causing overflow
0136   BOOST_CHECK_EQUAL(safeExp<FloatType>(600.0),
0137                     std::numeric_limits<FloatType>::infinity());
0138 }
0139
0140 BOOST_AUTO_TEST_CASE(safeExpFloat) {
0141   using FloatType = float;
0142
0143   // Values within the safe range
0144   BOOST_CHECK_CLOSE(safeExp<FloatType>(0.0f), std::exp(0.0f), 1e-8);
0145   BOOST_CHECK_CLOSE(safeExp<FloatType>(1.0f), std::exp(1.0f), 1e-8);
0146   BOOST_CHECK_CLOSE(safeExp<FloatType>(-1.0f), std::exp(-1.0f), 1e-8);
0147
0148   // Values causing underflow
0149   BOOST_CHECK_EQUAL(safeExp<FloatType>(-60.0f), 0.0f);
0150
0151   // Values causing overflow
0152   BOOST_CHECK_EQUAL(safeExp<FloatType>(60.0f),
0153                     std::numeric_limits<FloatType>::infinity());
0154 }
0155
0156 BOOST_AUTO_TEST_SUITE_END()
0157
0158 BOOST_AUTO_TEST_SUITE(MatrixIdxUnrolling)
0159
0160 template <std::size_t N>
0161 bool testUnrolling()
0162   requires(N > 1)
0163 {
0164   /// Fill a test matrix with unique values
0165   Acts::ActsSquareMatrix<N> testMatrix{Acts::ActsSquareMatrix<N>::Zero()};
0166   std::size_t counter{N};
0167   for (std::size_t i = 0; i < N; ++i) {
0168     for (std::size_t j = 0; j <= i; ++j) {
0169       testMatrix(j, i) = testMatrix(i, j) = counter;
0170       ++counter;
0171     }
0172   }
0173   /// Then unroll the matrix into a vector
0174   std::vector<double> vec{};
0175   for (std::size_t i = 0; i < N; ++i) {
0176     for (std::size_t j = 0; j <= i; ++j) {
0177       vec.push_back(testMatrix(i, j));
0178     }
0179   }
0180
0181   if (vec.size() != Acts::sumUpToN(N)) {
0182     std::cout << "Compressed vector size mismatch: expected "
0183               << Acts::sumUpToN(N) << ", got " << vec.size() << std::endl;
0184     return false;
0185   }
0186   /// Next test whether the indices are correct
0187   for (std::size_t i = 0; i < N; ++i) {
0188     for (std::size_t j = 0; j <= i; ++j) {
0189       const auto idx = Acts::vecIdxFromSymMat<N>(i, j);
0190       if (idx >= vec.size()) {
0191         std::cout << "Index out of bounds: " << idx << " for size "
0192                   << vec.size() << std::endl;
0193         return false;
0194       }
0195       if (idx != Acts::vecIdxFromSymMat<N>(j, i)) {
0196         std::cout << "Index mismatch: expected "
0197                   << Acts::vecIdxFromSymMat<N>(i, j) << ", got " << idx
0198                   << std::endl;
0199         return false;
0200       }
0201       if (vec[idx] != testMatrix(i, j)) {
0202         std::cout << "Value mismatch at index " << idx << ": expected "
0203                   << testMatrix(i, j) << ", got " << vec[idx] << std::endl;
0204         return false;
0205       }
0206       /// Finally, check whether the indices can be recovered
0207       const auto [iBack, jBack] = Acts::symMatIndices<N>(idx);
0208       if (iBack != i || jBack != j) {
0209         std::cout << "Index mismatch: expected (" << i << ", " << j
0210                   << "), got (" << iBack << ", " << jBack << ") for index "
0211                   << idx << std::endl;
0212         return false;
0213       }
0214     }
0215   }
0216   if constexpr (N > 2) {
0217     return testUnrolling<N - 1>();
0218   }
0219   return true;
0220 }
0221 BOOST_AUTO_TEST_CASE(matrixIdxUnrolling) {
0222   BOOST_CHECK_EQUAL(testUnrolling<20>(), true);
0223 }
0224
0225 BOOST_AUTO_TEST_SUITE_END()
0226
0227 BOOST_AUTO_TEST_SUITE_END()
0228
0229 }  // namespace Acts::Test