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/AlgebraHelpers.hpp"
0012 #include "Acts/Utilities/detail/Polynomials.hpp"
0013
0014 #include <algorithm>
0015 #include <numeric>
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 BOOST_AUTO_TEST_CASE(PowerTests) {
0094   for (unsigned p = 0; p <= 15; ++p) {
0095     BOOST_CHECK_EQUAL(std::pow(2., p), Acts::pow(2., p));
0096     BOOST_CHECK_EQUAL(std::pow(0.5, p), Acts::pow(0.5, p));
0097     for (std::size_t k = 1; k <= 15; ++k) {
0098       BOOST_CHECK_EQUAL(std::pow(k, p), Acts::pow(k, p));
0099     }
0100   }
0101   for (int p = 0; p <= 15; ++p) {
0102     BOOST_CHECK_EQUAL(std::pow(2., p), Acts::pow(2., p));
0103     BOOST_CHECK_EQUAL(std::pow(2., -p), Acts::pow(2., -p));
0104     BOOST_CHECK_EQUAL(std::pow(0.5, p), Acts::pow(0.5, p));
0105
0106     BOOST_CHECK_EQUAL(std::pow(0.5, -p), Acts::pow(0.5, -p));
0107   }
0108 }
0109
0110 BOOST_AUTO_TEST_CASE(SumOfIntegers) {
0111   std::array<unsigned, 100> numberSeq{Acts::filledArray<unsigned, 100>(1)};
0112   std::iota(numberSeq.begin(), numberSeq.end(), 1);
0113   for (unsigned i = 1; i <= numberSeq.size(); ++i) {
0114     const unsigned sum =
0115         std::accumulate(numberSeq.begin(), numberSeq.begin() + i, 0);
0116     BOOST_CHECK_EQUAL(sum, Acts::sumUpToN(i));
0117   }
0118 }
0119
0120 BOOST_AUTO_TEST_CASE(BinomialTests) {
0121   BOOST_CHECK_EQUAL(Acts::binomial(1u, 1u), 1u);
0122   for (unsigned n = 2; n <= 10; ++n) {
0123     /// Check that the binomial of (n 1 is always n)
0124     BOOST_CHECK_EQUAL(Acts::binomial(n, 1u), n);
0125     for (unsigned k = 1; k <= n; ++k) {
0126       /// Use recursive formula
0127       ///  n      n -1       n -1
0128       ///     =          +
0129       ///  k      k -1        k
0130       std::cout << "n: " << n << ", k: " << k
0131                 << ", binom(n,k): " << Acts::binomial(n, k)
0132                 << ", binom(n-1, k-1): " << Acts::binomial(n - 1, k - 1)
0133                 << ", binom(n-1,k): " << Acts::binomial(n - 1, k) << std::endl;
0134       BOOST_CHECK_EQUAL(Acts::binomial(n, k), Acts::binomial(n - 1, k - 1) +
0135                                                   Acts::binomial(n - 1, k));
0136       BOOST_CHECK_EQUAL(Acts::binomial(n, k), Acts::binomial(n, n - k));
0137     }
0138   }
0139 }
0140 BOOST_AUTO_TEST_CASE(LegendrePolynomials) {
0141   using namespace Acts::detail::Legendre;
0142   using namespace Acts::detail;
0143   std::cout << "Legdnre coefficients L=0: " << coefficients<0>() << std::endl;
0144   std::cout << "Legdnre coefficients L=1: " << coefficients<1>() << std::endl;
0145   std::cout << "Legdnre coefficients L=2: " << coefficients<2>() << std::endl;
0146   std::cout << "Legdnre coefficients L=3: " << coefficients<3>() << std::endl;
0147   std::cout << "Legdnre coefficients L=4: " << coefficients<4>() << std::endl;
0148   std::cout << "Legdnre coefficients L=5: " << coefficients<5>() << std::endl;
0149   std::cout << "Legdnre coefficients L=6: " << coefficients<6>() << std::endl;
0150   for (unsigned order = 0; order < 10; ++order) {
0151     const double sign = (order % 2 == 0 ? 1. : -1.);
0152     BOOST_CHECK_EQUAL(withinTolerance(legendrePoly(1., order), 1.), true);
0153     BOOST_CHECK_EQUAL(withinTolerance(legendrePoly(-1., order), sign * 1.),
0154                       true);
0155     for (double x = -1.; x <= 1.; x += stepSize) {
0156       const double evalX = legendrePoly(x, order);
0157       const double evalDx = legendrePoly(x, order, 1);
0158       const double evalD2x = legendrePoly(x, order, 2);
0159       /// Check whether the polynomial solves the legendre differental equation
0160       /// (1-x^{2}) d^P_{l}(x) /d^{x} -2x * d^P_{l}(x) / dx + l*(l+1) P_{l}(x) =
0161       /// 0;
0162       const double legendreEq = (1. - Acts::square(x)) * evalD2x -
0163                                 2. * x * evalDx + order * (order + 1) * evalX;
0164       BOOST_CHECK_EQUAL(withinTolerance(legendreEq, 0.), true);
0165       BOOST_CHECK_EQUAL(withinTolerance(evalX, sign * legendrePoly(-x, order)),
0166                         true);
0167       if (order == 0) {
0168         continue;
0169       }
0170       const double evalX_P1 = legendrePoly(x, order + 1);
0171       const double evalX_M1 = legendrePoly(x, order - 1);
0172       /// Recursion formular
0173       /// (n+1) P_{n+1}(x) = (2n+1)*x*P_{n}(x) - n * P_{n-1}(x)
0174       BOOST_CHECK_EQUAL(
0175           withinTolerance((order + 1) * evalX_P1,
0176                           (2 * order + 1) * x * evalX - order * evalX_M1),
0177           true);
0178     }
0179   }
0180 }
0181
0182 BOOST_AUTO_TEST_CASE(ChebychevPolynomials) {
0183   using namespace Acts::detail::Chebychev;
0184   using namespace Acts::detail;
0185
0186   /// Properties of the Chebychev polynomials
0187   ///  T_{n}(x) = (-1)^{n} T_n(-x)
0188   ///  U_{n}(x) = (-1)^{n} U_n(-x)
0189   ///  T_{n}(1) = 1
0190   ///  U_{n}(1) = n+1
0191   ///  T_{n+1}(x) = 2*x*T_{n}(x) - T_{n-1}(x)
0192   ///             = x*T_{n}(x) - (1-x^{2}) U_{n-1}(x)
0193   ///  U_{n+1}(x) = 2*x*U_{n}(x) - U_{n-1}(x)
0194   ///             = x*U_{n}(x) + T_{n+1}(x)
0195   for (unsigned order = 0; order < 10; ++order) {
0196     const double sign = (order % 2 == 0 ? 1. : -1.);
0197     const double T_n1 = chebychevPolyTn(1., order);
0198     const double U_n1 = chebychevPolyUn(1., order);
0199
0200     std::cout << "Order: " << order << " T(1)=" << T_n1 << ", U(1)=" << U_n1
0201               << std::endl;
0202     BOOST_CHECK_EQUAL(withinTolerance(T_n1, 1.), true);
0203     BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyTn(-1., order), sign), true);
0204
0205     BOOST_CHECK_EQUAL(withinTolerance(U_n1, order + 1), true);
0206     BOOST_CHECK_EQUAL(withinTolerance(U_n1, sign * chebychevPolyUn(-1., order)),
0207                       true);
0208     if (order == 0) {
0209       continue;
0210     }
0211     for (double x = -1.; x <= 1.; x += stepSize) {
0212       const double U_np1 = chebychevPolyUn(x, order + 1);
0213       const double U_n = chebychevPolyUn(x, order);
0214       const double U_nm1 = chebychevPolyUn(x, order - 1);
0215
0216       const double T_np1 = chebychevPolyTn(x, order + 1);
0217       const double T_n = chebychevPolyTn(x, order);
0218       const double T_nm1 = chebychevPolyTn(x, order - 1);
0219
0220       BOOST_TEST_MESSAGE(
0221           "Order: " << order << ", x=" << x << ", U_{n+1}(x) = " << U_np1
0222                     << ", 2*x*U_{n}(x) - U_{n-1}=" << (2. * x * U_n - U_nm1)
0223                     << ", x*U_{n}(x) + T_{n+1}(x)=" << (x * U_n + T_np1));
0224
0225       BOOST_CHECK_EQUAL(withinTolerance(U_np1, 2. * x * U_n - U_nm1), true);
0226       BOOST_CHECK_EQUAL(withinTolerance(U_np1, x * U_n + T_np1), true);
0227
0228       BOOST_CHECK_EQUAL(withinTolerance(U_n, sign * chebychevPolyUn(-x, order)),
0229                         true);
0230
0231       BOOST_TEST_MESSAGE("x=" << x << ", T_{n+1}(x) = " << T_np1
0232                               << ", 2*x*T_{n}(x) - T_{n-1}="
0233                               << (2. * x * T_n - T_nm1));
0234       BOOST_CHECK_EQUAL(withinTolerance(T_np1, 2. * x * T_n - T_nm1), true);
0235       BOOST_CHECK_EQUAL(withinTolerance(T_np1, x * T_n - (1 - x * x) * U_nm1),
0236                         true);
0237
0238       BOOST_CHECK_EQUAL(withinTolerance(T_n, sign * chebychevPolyTn(-x, order)),
0239                         true);
0240
0241       /// Check derivative
0242       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyTn(x, order, 1),
0243                                         Chebychev::evalFirstKind(x, order, 1)),
0244                         true);
0245       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyUn(x, order, 1),
0246                                         Chebychev::evalSecondKind(x, order, 1)),
0247                         true);
0248       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyTn(x, order, 2),
0249                                         Chebychev::evalFirstKind(x, order, 2)),
0250                         true);
0251       BOOST_CHECK_EQUAL(withinTolerance(chebychevPolyUn(x, order, 2),
0252                                         Chebychev::evalSecondKind(x, order, 2)),
0253                         true);
0254     }
0255   }
0256 }
0257 BOOST_AUTO_TEST_SUITE_END()
0258
0259 }  // namespace ActsTests