libs/math/test/cardinal_quintic_b_spline_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2019
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7
   8 #include "math_unit_test.hpp"
   9 #include <numeric>
  10 #include <utility>
  11 #include <boost/math/interpolators/cardinal_quintic_b_spline.hpp>
  12 #ifdef BOOST_HAS_FLOAT128
  13 #include <boost/multiprecision/float128.hpp>
  14 using boost::multiprecision::float128;
  15 #endif
  16 using boost::math::interpolators::cardinal_quintic_b_spline;
  17
  18 template<class Real>
  19 void test_constant()
  20 {
  21     Real c = 7.5;
  22     Real t0 = 0;
  23     Real h = Real(1)/Real(16);
  24     size_t n = 513;
  25     std::vector<Real> v(n, c);
  26     std::pair<Real, Real> left_endpoint_derivatives{0, 0};
  27     std::pair<Real, Real> right_endpoint_derivatives{0, 0};
  28     auto qbs = cardinal_quintic_b_spline<Real>(v.data(), v.size(), t0, h, left_endpoint_derivatives, right_endpoint_derivatives);
  29
  30     size_t i = 0;
  31     while (i < n) {
  32       Real t = t0 + i*h;
  33       CHECK_ULP_CLOSE(c, qbs(t), 3);
  34       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 400*std::numeric_limits<Real>::epsilon());
  35       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 60000*std::numeric_limits<Real>::epsilon());
  36       ++i;
  37     }
  38
  39     i = 0;
  40     while (i < n - 1) {
  41       Real t = t0 + i*h + h/2;
  42       CHECK_ULP_CLOSE(c, qbs(t), 5);
  43       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 600*std::numeric_limits<Real>::epsilon());
  44       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 30000*std::numeric_limits<Real>::epsilon());
  45       t = t0 + i*h + h/4;
  46       CHECK_ULP_CLOSE(c, qbs(t), 4);
  47       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 600*std::numeric_limits<Real>::epsilon());
  48       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 10000*std::numeric_limits<Real>::epsilon());
  49       ++i;
  50     }
  51 }
  52
  53 template<class Real>
  54 void test_constant_estimate_derivatives()
  55 {
  56     Real c = 7.5;
  57     Real t0 = 0;
  58     Real h = Real(1)/Real(16);
  59     size_t n = 513;
  60     std::vector<Real> v(n, c);
  61     auto qbs = cardinal_quintic_b_spline<Real>(v.data(), v.size(), t0, h);
  62
  63     size_t i = 0;
  64     while (i < n) {
  65       Real t = t0 + i*h;
  66       CHECK_ULP_CLOSE(c, qbs(t), 3);
  67       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 1200*std::numeric_limits<Real>::epsilon());
  68       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 200000*std::numeric_limits<Real>::epsilon());
  69       ++i;
  70     }
  71
  72     i = 0;
  73     while (i < n - 1) {
  74       Real t = t0 + i*h + h/2;
  75       CHECK_ULP_CLOSE(c, qbs(t), 8);
  76       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 1200*std::numeric_limits<Real>::epsilon());
  77       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 80000*std::numeric_limits<Real>::epsilon());
  78       t = t0 + i*h + h/4;
  79       CHECK_ULP_CLOSE(c, qbs(t), 5);
  80       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 1200*std::numeric_limits<Real>::epsilon());
  81       CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 38000*std::numeric_limits<Real>::epsilon());
  82       ++i;
  83     }
  84 }
  85
  86
  87 template<class Real>
  88 void test_linear()
  89 {
  90     using std::abs;
  91     Real m = 8.3;
  92     Real b = 7.2;
  93     Real t0 = 0;
  94     Real h = Real(1)/Real(16);
  95     size_t n = 512;
  96     std::vector<Real> y(n);
  97     for (size_t i = 0; i < n; ++i) {
  98       Real t = i*h;
  99       y[i] = m*t + b;
 100     }
 101     std::pair<Real, Real> left_endpoint_derivatives{m, 0};
 102     std::pair<Real, Real> right_endpoint_derivatives{m, 0};
 103     auto qbs = cardinal_quintic_b_spline<Real>(y.data(), y.size(), t0, h, left_endpoint_derivatives, right_endpoint_derivatives);
 104
 105     size_t i = 0;
 106     while (i < n) {
 107       Real t = t0 + i*h;
 108       if (!CHECK_ULP_CLOSE(m*t+b, qbs(t), 3)) {
 109           std::cerr << "  Problem at t = " << t << "\n";
 110       }
 111       if(!CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 100*abs(m*t+b)*std::numeric_limits<Real>::epsilon())) {
 112           std::cerr << "  Problem at t = " << t << "\n";
 113       }
 114       if(!CHECK_MOLLIFIED_CLOSE(0, qbs.double_prime(t), 10000*abs(m*t+b)*std::numeric_limits<Real>::epsilon())) {
 115           std::cerr << "  Problem at t = " << t << "\n";
 116       }
 117       ++i;
 118     }
 119
 120     i = 0;
 121     while (i < n - 1) {
 122       Real t = t0 + i*h + h/2;
 123       if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) {
 124           std::cerr << "  Problem at t = " << t << "\n";
 125       }
 126       CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
 127       t = t0 + i*h + h/4;
 128       if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) {
 129           std::cerr << "  Problem at t = " << t << "\n";
 130       }
 131       CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 3000*std::numeric_limits<Real>::epsilon());
 132       ++i;
 133     }
 134 }
 135
 136 template<class Real>
 137 void test_linear_estimate_derivatives()
 138 {
 139     using std::abs;
 140     Real m = 8.3;
 141     Real b = 7.2;
 142     Real t0 = 0;
 143     Real h = Real(1)/Real(16);
 144     size_t n = 512;
 145     std::vector<Real> y(n);
 146     for (size_t i = 0; i < n; ++i) {
 147       Real t = i*h;
 148       y[i] = m*t + b;
 149     }
 150
 151     auto qbs = cardinal_quintic_b_spline<Real>(y.data(), y.size(), t0, h);
 152
 153     size_t i = 0;
 154     while (i < n) {
 155       Real t = t0 + i*h;
 156       if (!CHECK_ULP_CLOSE(m*t+b, qbs(t), 3)) {
 157           std::cerr << "  Problem at t = " << t << "\n";
 158       }
 159       if(!CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 100*abs(m*t+b)*std::numeric_limits<Real>::epsilon())) {
 160           std::cerr << "  Problem at t = " << t << "\n";
 161       }
 162       if(!CHECK_MOLLIFIED_CLOSE(0, qbs.double_prime(t), 20000*abs(m*t+b)*std::numeric_limits<Real>::epsilon())) {
 163           std::cerr << "  Problem at t = " << t << "\n";
 164       }
 165       ++i;
 166     }
 167
 168     i = 0;
 169     while (i < n - 1) {
 170       Real t = t0 + i*h + h/2;
 171       if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 5)) {
 172           std::cerr << "  Problem at t = " << t << "\n";
 173       }
 174       CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
 175       t = t0 + i*h + h/4;
 176       if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) {
 177           std::cerr << "  Problem at t = " << t << "\n";
 178       }
 179       CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 3000*std::numeric_limits<Real>::epsilon());
 180       ++i;
 181     }
 182 }
 183
 184
 185 template<class Real>
 186 void test_quadratic()
 187 {
 188     Real a = Real(1)/Real(16);
 189     Real b = -3.5;
 190     Real c = -9;
 191     Real t0 = 0;
 192     Real h = Real(1)/Real(16);
 193     size_t n = 513;
 194     std::vector<Real> y(n);
 195     for (size_t i = 0; i < n; ++i) {
 196       Real t = i*h;
 197       y[i] = a*t*t + b*t + c;
 198     }
 199     Real t_max = t0 + (n-1)*h;
 200     std::pair<Real, Real> left_endpoint_derivatives{b, 2*a};
 201     std::pair<Real, Real> right_endpoint_derivatives{2*a*t_max + b, 2*a};
 202
 203     auto qbs = cardinal_quintic_b_spline<Real>(y, t0, h, left_endpoint_derivatives, right_endpoint_derivatives);
 204
 205     size_t i = 0;
 206     while (i < n) {
 207       Real t = t0 + i*h;
 208       CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 3);
 209       ++i;
 210     }
 211
 212     i = 0;
 213     while (i < n -1) {
 214       Real t = t0 + i*h + h/2;
 215       if(!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 5)) {
 216           std::cerr << "  Problem at abscissa t = " << t << "\n";
 217       }
 218
 219       t = t0 + i*h + h/4;
 220       if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 5)) {
 221           std::cerr << "  Problem abscissa t = " << t << "\n";
 222       }
 223       ++i;
 224     }
 225 }
 226
 227
 228 template<class Real>
 229 void test_quadratic_estimate_derivatives()
 230 {
 231     Real a = Real(1)/Real(16);
 232     Real b = -3.5;
 233     Real c = -9;
 234     Real t0 = 0;
 235     Real h = Real(1)/Real(16);
 236     size_t n = 513;
 237     std::vector<Real> y(n);
 238     for (size_t i = 0; i < n; ++i) {
 239       Real t = i*h;
 240       y[i] = a*t*t + b*t + c;
 241     }
 242     auto qbs = cardinal_quintic_b_spline<Real>(y, t0, h);
 243
 244     size_t i = 0;
 245     while (i < n) {
 246       Real t = t0 + i*h;
 247       CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 3);
 248       ++i;
 249     }
 250
 251     i = 0;
 252     while (i < n -1) {
 253       Real t = t0 + i*h + h/2;
 254       if(!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 10)) {
 255           std::cerr << "  Problem at abscissa t = " << t << "\n";
 256       }
 257
 258       t = t0 + i*h + h/4;
 259       if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 6)) {
 260           std::cerr << "  Problem abscissa t = " << t << "\n";
 261       }
 262       ++i;
 263     }
 264 }
 265
 266
 267 int main()
 268 {
 269     test_constant<double>();
 270     test_constant<long double>();
 271
 272     test_constant_estimate_derivatives<double>();
 273     test_constant_estimate_derivatives<long double>();
 274
 275     test_linear<float>();
 276     test_linear<double>();
 277     test_linear<long double>();
 278
 279     test_linear_estimate_derivatives<double>();
 280     test_linear_estimate_derivatives<long double>();
 281
 282     test_quadratic<double>();
 283     test_quadratic<long double>();
 284
 285     test_quadratic_estimate_derivatives<double>();
 286     test_quadratic_estimate_derivatives<long double>();
 287
 288
 289     #ifdef BOOST_HAS_FLOAT128
 290         test_constant<float128>();
 291         test_linear<float128>();
 292         test_linear_estimate_derivatives<float128>();
 293     #endif
 294
 295     return boost::math::test::report_errors();
 296 }