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