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