Core/Utilities/PolynomialTests.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/detail/Polynomials.hpp"
0012
0013 #include <algorithm>
0014 #include <format>
0015
0016 namespace {
0017 constexpr double stepSize = 1.e-4;
0018
0019 constexpr bool withinTolerance(const double value, const double expect) {
0020   constexpr double tolerance = 1.e-10;
0021   return std::abs(value - expect) < tolerance;
0022 }
0023 template <typename T, std::size_t D>
0024 std::ostream& operator<<(std::ostream& ostr,
0025
0026                          const std::array<T, D>& arr
0027
0028 ) {
0029   ostr << "[";
0030   for (std::size_t t = 0; t < D; ++t) {
0031     ostr << arr[t];
0032     if (t + 1 != D) {
0033       ostr << ", ";
0034     }
0035   }
0036   ostr << "]";
0037   return ostr;
0038 }
0039
0040 template <std::size_t O, std::size_t D>
0041 bool checkDerivative(const std::array<double, O>& unitArray,
0042                      const std::array<double, (O - D + 1)>& recentDeriv) {
0043   const auto nthDerivative = Acts::detail::derivativeCoefficients<D>(unitArray);
0044   const auto firstDeriv = Acts::detail::derivativeCoefficients<1>(recentDeriv);
0045   std::cout << D << "-th derivative: " << nthDerivative << std::endl;
0046   if (nthDerivative.size() != firstDeriv.size()) {
0047     std::cout << __func__ << "<" << D << ">" << ":" << __LINE__
0048               << " - nthDerivative: " << nthDerivative.size()
0049               << " vs. firstDeriv: " << firstDeriv.size() << std::endl;
0050     return false;
0051   }
0052   bool good{true};
0053   for (std::size_t i = 0; i < nthDerivative.size(); ++i) {
0054     if (!withinTolerance(firstDeriv[i], nthDerivative[i])) {
0055       std::cout << __func__ << "<" << D << ">" << ":" << __LINE__
0056                 << " - nthDerivative: " << nthDerivative[i]
0057                 << " vs. firstDeriv: " << firstDeriv[i] << std::endl;
0058       good = false;
0059     }
0060   }
0061   for (std::size_t i = 0; i < firstDeriv.size(); ++i) {
0062     if (!withinTolerance(firstDeriv[i], (i + 1) * recentDeriv[i + 1])) {
0063       std::cout << __func__ << "<" << D << ">" << ":" << __LINE__
0064                 << " - firstDeriv: " << firstDeriv[i]
0065                 << " vs. expect: " << ((i + 1) * recentDeriv[i + 1])
0066                 << std::endl;
0067       good = false;
0068     }
0069   }
0070   if (!good) {
0071     return false;
0072   }
0073   if constexpr (D + 1 < O) {
0074     if (!checkDerivative<O, D + 1>(unitArray, firstDeriv)) {
0075       return false;
0076     }
0077   }
0078   return good;
0079 }
0080
0081 }  // namespace
0082
0083 namespace ActsTests {
0084
0085 BOOST_AUTO_TEST_SUITE(UtilitiesSuite)
0086
0087 BOOST_AUTO_TEST_CASE(DerivativeCoeffs) {
0088   constexpr std::size_t order = 20;
0089   std::array<double, order> unitCoeffs{Acts::filledArray<double, order>(1)};
0090   const bool result = checkDerivative<order, 1>(unitCoeffs, unitCoeffs);
0091   BOOST_CHECK_EQUAL(result, true);
0092 }
0093
0094 BOOST_AUTO_TEST_CASE(LegendrePolynomials) {
0095   using namespace Acts::detail::Legendre;
0096   using namespace Acts::detail;
0097   std::cout << "Legendre coefficients L=0: " << coefficients<0>() << std::endl;
0098   std::cout << "Legendre coefficients L=1: " << coefficients<1>() << std::endl;
0099   std::cout << "Legendre coefficients L=2: " << coefficients<2>() << std::endl;
0100   std::cout << "Legendre coefficients L=3: " << coefficients<3>() << std::endl;
0101   std::cout << "Legendre coefficients L=4: " << coefficients<4>() << std::endl;
0102   std::cout << "Legendre coefficients L=5: " << coefficients<5>() << std::endl;
0103   std::cout << "Legendre coefficients L=6: " << coefficients<6>() << std::endl;
0104   for (unsigned order = 0; order < 10; ++order) {
0105     const double sign = (order % 2 == 0 ? 1. : -1.);
0106     BOOST_CHECK_EQUAL(withinTolerance(legendrePoly(1., order), 1.), true);
0107     BOOST_CHECK_EQUAL(withinTolerance(legendrePoly(-1., order), sign * 1.),
0108                       true);
0109     for (double x = -1.; x <= 1.; x += stepSize) {
0110       const double evalX = legendrePoly(x, order);
0111       const double evalDx = legendrePoly(x, order, 1);
0112       const double evalD2x = legendrePoly(x, order, 2);
0113       /// Check whether the polynomial solves the legendre differental equation
0114       /// (1-x^{2}) d^P_{l}(x) /d^{x} -2x * d^P_{l}(x) / dx + l*(l+1) P_{l}(x) =
0115       /// 0;
0116       const double legendreEq = (1. - Acts::square(x)) * evalD2x -
0117                                 2. * x * evalDx + order * (order + 1) * evalX;
0118
0119       BOOST_CHECK_EQUAL(withinTolerance(legendreEq, 0.), true);
0120       BOOST_CHECK_EQUAL(withinTolerance(evalX, sign * legendrePoly(-x, order)),
0121                         true);
0122       if (!withinTolerance(legendreEq, 0.)) {
0123         std::cout << std::format(
0124                          "x: {1:.3f}, P_{0}(x)={2:.3f}, d/dx P_{0}(x)={3:.3f}, "
0125                          "d^2/d^2x P_{0}(x)={4:.3f} --> legendreEq: {5:.3f}",
0126                          order, x, evalX, evalDx, evalD2x, legendreEq)
0127                   << std::endl;
0128       }
0129       if (order == 0) {
0130         continue;
0131       }
0132       const double evalX_P1 = legendrePoly(x, order + 1);
0133       const double evalX_M1 = legendrePoly(x, order - 1);
0134       /// Recursion formula
0135       /// (n+1) P_{n+1}(x) = (2n+1)*x*P_{n}(x) - n * P_{n-1}(x)
0136       BOOST_CHECK_EQUAL(
0137           withinTolerance((order + 1) * evalX_P1,
0138                           (2 * order + 1) * x * evalX - order * evalX_M1),
0139           true);
0140       for (unsigned d = 0; d <= std::min(2u, order); ++d) {
0141         const double constExpEval = legendrePoly(x, order, d);
0142         const double runTimeEval = Legendre::evaluate(x, order, d);
0143         BOOST_CHECK_EQUAL(withinTolerance(constExpEval, runTimeEval), true);
0144         if (!withinTolerance(constExpEval, runTimeEval)) {
0145           std::cout << "Compile time evaluation (" << constExpEval
0146                     << ") vs. run time evaluation (" << runTimeEval
0147                     << ") differ - order: " << order << ", derivative: " << d
0148                     << std::endl;
0149         }
0150       }
0151     }
0152   }
0153 }
0154
0155 BOOST_AUTO_TEST_CASE(ChebychevPolynomials) {
0156   using namespace Acts::detail::Chebychev;
0157   using namespace Acts::detail;
0158
0159   /// Properties of the Chebychev polynomials
0160   ///  T_{n}(x) = (-1)^{n} T_n(-x)
0161   ///  U_{n}(x) = (-1)^{n} U_n(-x)
0162   ///  T_{n}(1) = 1
0163   ///  U_{n}(1) = n+1
0164   ///  T_{n+1}(x) = 2*x*T_{n}(x) - T_{n-1}(x)
0165   ///             = x*T_{n}(x) - (1-x^{2}) U_{n-1}(x)
0166   ///  U_{n+1}(x) = 2*x*U_{n}(x) - U_{n-1}(x)
0167   ///             = x*U_{n}(x) + T_{n+1}(x)
0168   for (unsigned order = 0; order < 10; ++order) {
0169     const double sign = (order % 2 == 0 ? 1. : -1.);
0170     const double T_n1 = chebychevPolyTn(1., order);
0171     const double U_n1 = chebychevPolyUn(1., order);
0172
0173     std::cout << "Order: " << order << " T(1)=" << T_n1 << ", U(1)=" << U_n1
0174               << std::endl;
0175     BOOST_CHECK_EQUAL(withinTolerance(T_n1, 1.), true);
0176     BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyTn(-1., order), sign), true);
0177
0178     BOOST_CHECK_EQUAL(withinTolerance(U_n1, order + 1), true);
0179     BOOST_CHECK_EQUAL(withinTolerance(U_n1, sign * chebychevPolyUn(-1., order)),
0180                       true);
0181     if (order == 0) {
0182       continue;
0183     }
0184     for (double x = -1.; x <= 1.; x += stepSize) {
0185       const double U_np1 = chebychevPolyUn(x, order + 1);
0186       const double U_n = chebychevPolyUn(x, order);
0187       const double U_nm1 = chebychevPolyUn(x, order - 1);
0188
0189       const double T_np1 = chebychevPolyTn(x, order + 1);
0190       const double T_n = chebychevPolyTn(x, order);
0191       const double T_nm1 = chebychevPolyTn(x, order - 1);
0192
0193       BOOST_TEST_MESSAGE(
0194           "Order: " << order << ", x=" << x << ", U_{n+1}(x) = " << U_np1
0195                     << ", 2*x*U_{n}(x) - U_{n-1}=" << (2. * x * U_n - U_nm1)
0196                     << ", x*U_{n}(x) + T_{n+1}(x)=" << (x * U_n + T_np1));
0197
0198       BOOST_CHECK_EQUAL(withinTolerance(U_np1, 2. * x * U_n - U_nm1), true);
0199       BOOST_CHECK_EQUAL(withinTolerance(U_np1, x * U_n + T_np1), true);
0200
0201       BOOST_CHECK_EQUAL(withinTolerance(U_n, sign * chebychevPolyUn(-x, order)),
0202                         true);
0203
0204       BOOST_TEST_MESSAGE("x=" << x << ", T_{n+1}(x) = " << T_np1
0205                               << ", 2*x*T_{n}(x) - T_{n-1}="
0206                               << (2. * x * T_n - T_nm1));
0207       BOOST_CHECK_EQUAL(withinTolerance(T_np1, 2. * x * T_n - T_nm1), true);
0208       BOOST_CHECK_EQUAL(withinTolerance(T_np1, x * T_n - (1 - x * x) * U_nm1),
0209                         true);
0210
0211       BOOST_CHECK_EQUAL(withinTolerance(T_n, sign * chebychevPolyTn(-x, order)),
0212                         true);
0213
0214       /// Check derivative
0215       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyTn(x, order, 1),
0216                                         Chebychev::evalFirstKind(x, order, 1)),
0217                         true);
0218       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyUn(x, order, 1),
0219                                         Chebychev::evalSecondKind(x, order, 1)),
0220                         true);
0221       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyTn(x, order, 2),
0222                                         Chebychev::evalFirstKind(x, order, 2)),
0223                         true);
0224       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyUn(x, order, 2),
0225                                         Chebychev::evalSecondKind(x, order, 2)),
0226                         true);
0227     }
0228   }
0229 }
0230 BOOST_AUTO_TEST_SUITE_END()
0231
0232 }  // namespace ActsTests