ceph/src/boost/libs/math/test/quintic_hermite_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2020
   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 <vector>
  12 #include <array>
  13 #include <boost/random/uniform_real.hpp>
  14 #include <boost/random/mersenne_twister.hpp>
  15 #include <boost/math/interpolators/quintic_hermite.hpp>
  16 #include <boost/math/special_functions/next.hpp>
  17 #include <boost/circular_buffer.hpp>
  18 #ifdef BOOST_HAS_FLOAT128
  19 #include <boost/multiprecision/float128.hpp>
  20 using boost::multiprecision::float128;
  21 #endif
  22
  23
  24 using boost::math::interpolators::quintic_hermite;
  25 using boost::math::interpolators::cardinal_quintic_hermite;
  26 using boost::math::interpolators::cardinal_quintic_hermite_aos;
  27
  28 template<typename Real>
  29 void test_constant()
  30 {
  31     std::vector<Real> x{0,1,2,3, 9, 22, 81};
  32     std::vector<Real> y(x.size());
  33     std::vector<Real> dydx(x.size(), 0);
  34     std::vector<Real> d2ydx2(x.size(), 0);
  35     for (auto & t : y)
  36     {
  37         t = 7;
  38     }
  39
  40     auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
  41     for (Real t = 0; t <= 81; t += 0.25)
  42     {
  43         CHECK_ULP_CLOSE(Real(7), qh(t), 24);
  44         CHECK_ULP_CLOSE(Real(0), qh.prime(t), 24);
  45         CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 24);
  46     }
  47 }
  48
  49
  50 template<typename Real>
  51 void test_linear()
  52 {
  53     std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9};
  54     std::vector<Real> y = x;
  55     std::vector<Real> dydx(x.size(), 1);
  56     std::vector<Real> d2ydx2(x.size(), 0);
  57
  58     auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
  59
  60     for (Real t = 0; t <= 9; t += 0.25)
  61     {
  62         CHECK_ULP_CLOSE(Real(t), qh(t), 2);
  63         CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2);
  64         CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 2);
  65     }
  66
  67     boost::random::mt19937 rng;
  68     boost::random::uniform_real_distribution<Real> dis(0.5,1);
  69     x.resize(512);
  70     x[0] = dis(rng);
  71     Real xmin = x[0];
  72     for (size_t i = 1; i < x.size(); ++i)
  73     {
  74         x[i] = x[i-1] + dis(rng);
  75     }
  76     Real xmax = x.back();
  77
  78     y = x;
  79     dydx.resize(x.size(), 1);
  80     d2ydx2.resize(x.size(), 0);
  81
  82     qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
  83
  84     for (Real t = xmin; t <= xmax; t += 0.125)
  85     {
  86         CHECK_ULP_CLOSE(t, qh(t), 2);
  87         CHECK_ULP_CLOSE(Real(1), qh.prime(t), 100);
  88         CHECK_MOLLIFIED_CLOSE(Real(0), qh.double_prime(t), 200*std::numeric_limits<Real>::epsilon());
  89     }
  90 }
  91
  92 template<typename Real>
  93 void test_quadratic()
  94 {
  95
  96     std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9};
  97     std::vector<Real> y(x.size());
  98     for (size_t i = 0; i < y.size(); ++i)
  99     {
 100         y[i] = x[i]*x[i]/2;
 101     }
 102
 103     std::vector<Real> dydx(x.size());
 104     for (size_t i = 0; i < y.size(); ++i) {
 105         dydx[i] = x[i];
 106     }
 107
 108     std::vector<Real> d2ydx2(x.size(), 1);
 109
 110     auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
 111
 112     for (Real t = 0; t <= 9; t += 0.0078125)
 113     {
 114         CHECK_ULP_CLOSE(Real(t*t)/2, qh(t), 2);
 115         CHECK_ULP_CLOSE(t, qh.prime(t), 12);
 116         CHECK_ULP_CLOSE(Real(1), qh.double_prime(t), 32);
 117     }
 118
 119     boost::random::mt19937 rng;
 120     boost::random::uniform_real_distribution<Real> dis(0.5,1);
 121     x.resize(8);
 122     x[0] = dis(rng);
 123     Real xmin = x[0];
 124     for (size_t i = 1; i < x.size(); ++i)
 125     {
 126         x[i] = x[i-1] + dis(rng);
 127     }
 128     Real xmax = x.back();
 129
 130     y.resize(x.size());
 131     for (size_t i = 0; i < y.size(); ++i)
 132     {
 133         y[i] = x[i]*x[i]/2;
 134     }
 135
 136     dydx.resize(x.size());
 137     for (size_t i = 0; i < y.size(); ++i)
 138     {
 139         dydx[i] = x[i];
 140     }
 141
 142     d2ydx2.resize(x.size(), 1);
 143
 144     qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
 145
 146     for (Real t = xmin; t <= xmax; t += 0.125)
 147     {
 148         CHECK_ULP_CLOSE(Real(t*t)/2, qh(t), 4);
 149         CHECK_ULP_CLOSE(t, qh.prime(t), 53);
 150         CHECK_ULP_CLOSE(Real(1), qh.double_prime(t), 700);
 151     }
 152 }
 153
 154 template<typename Real>
 155 void test_cubic()
 156 {
 157
 158     std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9};
 159     std::vector<Real> y(x.size());
 160     for (size_t i = 0; i < y.size(); ++i)
 161     {
 162         y[i] = x[i]*x[i]*x[i];
 163     }
 164
 165     std::vector<Real> dydx(x.size());
 166     for (size_t i = 0; i < y.size(); ++i) {
 167         dydx[i] = 3*x[i]*x[i];
 168     }
 169
 170     std::vector<Real> d2ydx2(x.size());
 171     for (size_t i = 0; i < y.size(); ++i)
 172     {
 173         d2ydx2[i] = 6*x[i];
 174     }
 175
 176     auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
 177
 178     for (Real t = 0; t <= 9; t += 0.0078125)
 179     {
 180         CHECK_ULP_CLOSE(t*t*t, qh(t), 10);
 181         CHECK_ULP_CLOSE(3*t*t, qh.prime(t), 15);
 182         CHECK_ULP_CLOSE(6*t, qh.double_prime(t), 20);
 183     }
 184 }
 185
 186 template<typename Real>
 187 void test_quartic()
 188 {
 189
 190     std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9, 10, 11};
 191     std::vector<Real> y(x.size());
 192     for (size_t i = 0; i < y.size(); ++i)
 193     {
 194         y[i] = x[i]*x[i]*x[i]*x[i];
 195     }
 196
 197     std::vector<Real> dydx(x.size());
 198     for (size_t i = 0; i < y.size(); ++i)
 199     {
 200         dydx[i] = 4*x[i]*x[i]*x[i];
 201     }
 202
 203     std::vector<Real> d2ydx2(x.size());
 204     for (size_t i = 0; i < y.size(); ++i)
 205     {
 206         d2ydx2[i] = 12*x[i]*x[i];
 207     }
 208
 209     auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
 210
 211     for (Real t = 1; t <= 11; t += 0.0078125)
 212     {
 213         CHECK_ULP_CLOSE(t*t*t*t, qh(t), 100);
 214         CHECK_ULP_CLOSE(4*t*t*t, qh.prime(t), 100);
 215         CHECK_ULP_CLOSE(12*t*t, qh.double_prime(t), 100);
 216     }
 217 }
 218
 219
 220 template<typename Real>
 221 void test_interpolation_condition()
 222 {
 223     for (size_t n = 4; n < 50; ++n) {
 224         std::vector<Real> x(n);
 225         std::vector<Real> y(n);
 226         std::vector<Real> dydx(n);
 227         std::vector<Real> d2ydx2(n);
 228         boost::random::mt19937 rd;
 229         boost::random::uniform_real_distribution<Real> dis(0,1);
 230         Real x0 = dis(rd);
 231         x[0] = x0;
 232         y[0] = dis(rd);
 233         for (size_t i = 1; i < n; ++i) {
 234             x[i] = x[i-1] + dis(rd);
 235             y[i] = dis(rd);
 236             dydx[i] = dis(rd);
 237             d2ydx2[i] = dis(rd);
 238         }
 239
 240         auto x_copy = x;
 241         auto y_copy = y;
 242         auto dydx_copy = dydx;
 243         auto d2ydx2_copy = d2ydx2;
 244         auto s = quintic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy), std::move(d2ydx2_copy));
 245         //std::cout << "s = " << s << "\n";
 246         for (size_t i = 0; i < x.size(); ++i) {
 247             CHECK_ULP_CLOSE(y[i], s(x[i]), 2);
 248             CHECK_ULP_CLOSE(dydx[i], s.prime(x[i]), 2);
 249             CHECK_ULP_CLOSE(d2ydx2[i], s.double_prime(x[i]), 2);
 250         }
 251     }
 252 }
 253
 254 template<typename Real>
 255 void test_cardinal_constant()
 256 {
 257
 258     std::vector<Real> y(25);
 259     std::vector<Real> dydx(y.size(), 0);
 260     std::vector<Real> d2ydx2(y.size(), 0);
 261     for (auto & t : y) {
 262         t = 7;
 263     }
 264     Real x0 = 4;
 265     Real dx = Real(1)/Real(8);
 266
 267     auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx);
 268
 269     for (Real t = x0; t <= x0 + 24*dx; t += 0.25)
 270     {
 271         CHECK_ULP_CLOSE(Real(7), qh(t), 24);
 272         CHECK_ULP_CLOSE(Real(0), qh.prime(t), 24);
 273         CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 24);
 274     }
 275
 276     std::vector<std::array<Real, 3>> data(25);
 277     for (size_t i = 0; i < data.size(); ++i)
 278     {
 279         data[i][0] = 7;
 280         data[i][1] = 0;
 281         data[i][2] = 0;
 282     }
 283
 284     auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx);
 285     for (Real t = x0; t <= x0 + 24*dx; t += 0.25)
 286     {
 287         CHECK_ULP_CLOSE(Real(7), qh_aos(t), 24);
 288         CHECK_ULP_CLOSE(Real(0), qh_aos.prime(t), 24);
 289         CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 24);
 290     }
 291
 292     // Now check the boundaries:
 293     auto [tlo, thi] = qh.domain();
 294     int samples = 5000;
 295     int i = 0;
 296     while (i++ < samples)
 297     {
 298         CHECK_ULP_CLOSE(Real(7), qh(tlo), 2);
 299         CHECK_ULP_CLOSE(Real(7), qh(thi), 2);
 300         CHECK_ULP_CLOSE(Real(7), qh_aos(tlo), 2);
 301         CHECK_ULP_CLOSE(Real(7), qh_aos(thi), 2);
 302         CHECK_ULP_CLOSE(Real(0), qh.prime(tlo), 2);
 303         CHECK_ULP_CLOSE(Real(0), qh.prime(thi), 2);
 304         CHECK_ULP_CLOSE(Real(0), qh_aos.prime(tlo), 2);
 305         CHECK_ULP_CLOSE(Real(0), qh_aos.prime(thi), 2);
 306
 307         tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)());
 308         thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
 309     }
 310 }
 311
 312
 313 template<typename Real>
 314 void test_cardinal_linear()
 315 {
 316     std::vector<Real> y{0,1,2,3,4,5,6,7,8,9};
 317     Real x0 = 0;
 318     Real dx = 1;
 319     std::vector<Real> dydx(y.size(), 1);
 320     std::vector<Real> d2ydx2(y.size(), 0);
 321
 322     auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx);
 323
 324     for (Real t = 0; t <= 9; t += 0.25) {
 325         CHECK_ULP_CLOSE(Real(t), qh(t), 2);
 326         CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2);
 327         CHECK_ULP_CLOSE(Real(0), qh.double_prime(t), 2);
 328     }
 329
 330     std::vector<std::array<Real, 3>> data(10);
 331     for (size_t i = 0; i < data.size(); ++i) {
 332         data[i][0] = i;
 333         data[i][1] = 1;
 334         data[i][2] = 0;
 335     }
 336
 337     auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx);
 338
 339     for (Real t = 0; t <= 9; t += 0.25) {
 340         CHECK_ULP_CLOSE(Real(t), qh_aos(t), 2);
 341         CHECK_ULP_CLOSE(Real(1), qh_aos.prime(t), 2);
 342         CHECK_ULP_CLOSE(Real(0), qh_aos.double_prime(t), 2);
 343     }
 344
 345     // Now check the boundaries:
 346     auto [tlo, thi] = qh.domain();
 347     int samples = 5000;
 348     int i = 0;
 349     while (i++ < samples)
 350     {
 351         CHECK_ULP_CLOSE(Real(tlo), qh(tlo), 2);
 352         CHECK_ULP_CLOSE(Real(thi), qh(thi), 2);
 353         CHECK_ULP_CLOSE(Real(tlo), qh_aos(tlo), 2);
 354         CHECK_ULP_CLOSE(Real(thi), qh_aos(thi), 2);
 355         CHECK_ULP_CLOSE(Real(1), qh.prime(tlo), 2);
 356         CHECK_ULP_CLOSE(Real(1), qh.prime(thi), 128);
 357         CHECK_ULP_CLOSE(Real(1), qh_aos.prime(tlo), 2);
 358         CHECK_ULP_CLOSE(Real(1), qh_aos.prime(thi), 128);
 359
 360         tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)());
 361         thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
 362     }
 363 }
 364
 365 template<typename Real>
 366 void test_cardinal_quadratic()
 367 {
 368     Real x0 = 0;
 369     Real dx = 1;
 370     std::vector<Real> y(10);
 371     for (size_t i = 0; i < y.size(); ++i)
 372     {
 373         y[i] = i*i/Real(2);
 374     }
 375
 376     std::vector<Real> dydx(y.size());
 377     for (size_t i = 0; i < y.size(); ++i) {
 378         dydx[i] = i;
 379     }
 380
 381     std::vector<Real> d2ydx2(y.size(), 1);
 382
 383     auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx);
 384
 385     for (Real t = 0; t <= 9; t += 0.0078125) {
 386         Real computed = qh(t);
 387         CHECK_ULP_CLOSE(Real(t*t)/2, computed, 2);
 388         CHECK_ULP_CLOSE(t, qh.prime(t), 15);
 389         CHECK_ULP_CLOSE(Real(1), qh.double_prime(t), 32);
 390     }
 391
 392     std::vector<std::array<Real, 3>> data(10);
 393     for (size_t i = 0; i < data.size(); ++i) {
 394         data[i][0] = i*i/Real(2);
 395         data[i][1] = i;
 396         data[i][2] = 1;
 397     }
 398     auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx);
 399
 400     for (Real t = 0; t <= 9; t += 0.0078125)
 401     {
 402         Real computed = qh_aos(t);
 403         CHECK_ULP_CLOSE(Real(t*t)/2, computed, 2);
 404         CHECK_ULP_CLOSE(t, qh_aos.prime(t), 12);
 405         CHECK_ULP_CLOSE(Real(1), qh_aos.double_prime(t), 64);
 406     }
 407
 408         // Now check the boundaries:
 409     auto [tlo, thi] = qh.domain();
 410     int samples = 5000;
 411     int i = 0;
 412     while (i++ < samples)
 413     {
 414         CHECK_ULP_CLOSE(tlo*tlo/2, qh(tlo), 16);
 415         CHECK_ULP_CLOSE(thi*thi/2, qh(thi), 16);
 416         CHECK_ULP_CLOSE(tlo*tlo/2, qh_aos(tlo), 16);
 417         CHECK_ULP_CLOSE(thi*thi/2, qh_aos(thi), 16);
 418         CHECK_ULP_CLOSE(tlo, qh.prime(tlo), 16);
 419         CHECK_ULP_CLOSE(thi, qh.prime(thi), 64);
 420         CHECK_ULP_CLOSE(tlo, qh_aos.prime(tlo), 16);
 421         CHECK_ULP_CLOSE(thi, qh_aos.prime(thi), 64);
 422
 423         tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)());
 424         thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
 425     }
 426 }
 427
 428 template<typename Real>
 429 void test_cardinal_cubic()
 430 {
 431     Real x0 = 0;
 432     Real dx = 1;
 433     std::vector<Real> y(10);
 434     for (size_t i = 0; i < y.size(); ++i)
 435     {
 436         y[i] = i*i*i;
 437     }
 438
 439     std::vector<Real> dydx(y.size());
 440     for (size_t i = 0; i < y.size(); ++i) {
 441         dydx[i] = 3*i*i;
 442     }
 443
 444     std::vector<Real> d2ydx2(y.size());
 445     for (size_t i = 0; i < y.size(); ++i) {
 446         d2ydx2[i] = 6*i;
 447     }
 448
 449     auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx);
 450
 451     for (Real t = 0; t <= 9; t += 0.0078125)
 452     {
 453         Real computed = qh(t);
 454         CHECK_ULP_CLOSE(t*t*t, computed, 10);
 455         CHECK_ULP_CLOSE(3*t*t, qh.prime(t), 15);
 456         CHECK_ULP_CLOSE(6*t, qh.double_prime(t), 39);
 457     }
 458
 459     std::vector<std::array<Real, 3>> data(10);
 460     for (size_t i = 0; i < data.size(); ++i) {
 461         data[i][0] = i*i*i;
 462         data[i][1] = 3*i*i;
 463         data[i][2] = 6*i;
 464     }
 465
 466     auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx);
 467     for (Real t = 0; t <= 9; t += 0.0078125)
 468     {
 469         Real computed = qh_aos(t);
 470         CHECK_ULP_CLOSE(t*t*t, computed, 10);
 471         CHECK_ULP_CLOSE(3*t*t, qh_aos.prime(t), 15);
 472         CHECK_ULP_CLOSE(6*t, qh_aos.double_prime(t), 30);
 473     }
 474 }
 475
 476 template<typename Real>
 477 void test_cardinal_quartic()
 478 {
 479     Real x0 = 0;
 480     Real dx = 1;
 481     std::vector<Real> y(7);
 482     for (size_t i = 0; i < y.size(); ++i)
 483     {
 484         y[i] = i*i*i*i;
 485     }
 486
 487     std::vector<Real> dydx(y.size());
 488     for (size_t i = 0; i < y.size(); ++i) {
 489         dydx[i] = 4*i*i*i;
 490     }
 491
 492     std::vector<Real> d2ydx2(y.size());
 493     for (size_t i = 0; i < y.size(); ++i) {
 494         d2ydx2[i] = 12*i*i;
 495     }
 496
 497     auto qh = cardinal_quintic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), x0, dx);
 498
 499     for (Real t = 0; t <= 6; t += 0.0078125)
 500     {
 501         CHECK_ULP_CLOSE(Real(t*t*t*t), qh(t), 250);
 502         CHECK_ULP_CLOSE(4*t*t*t, qh.prime(t), 250);
 503         CHECK_ULP_CLOSE(12*t*t, qh.double_prime(t), 250);
 504     }
 505
 506     std::vector<std::array<Real, 3>> data(7);
 507     for (size_t i = 0; i < data.size(); ++i) {
 508         data[i][0] = i*i*i*i;
 509         data[i][1] = 4*i*i*i;
 510         data[i][2] = 12*i*i;
 511     }
 512
 513     auto qh_aos = cardinal_quintic_hermite_aos(std::move(data), x0, dx);
 514     for (Real t = 0; t <= 6; t += 0.0078125)
 515     {
 516         Real computed = qh_aos(t);
 517         CHECK_ULP_CLOSE(t*t*t*t, computed, 10);
 518         CHECK_ULP_CLOSE(4*t*t*t, qh_aos.prime(t), 64);
 519         CHECK_ULP_CLOSE(12*t*t, qh_aos.double_prime(t), 128);
 520     }
 521 }
 522
 523
 524 int main()
 525 {
 526     test_constant<float>();
 527     test_linear<float>();
 528     test_quadratic<float>();
 529     test_cubic<float>();
 530     test_quartic<float>();
 531     test_interpolation_condition<float>();
 532
 533     test_cardinal_constant<float>();
 534     test_cardinal_linear<float>();
 535     test_cardinal_quadratic<float>();
 536     test_cardinal_cubic<float>();
 537     test_cardinal_quartic<float>();
 538
 539     test_constant<double>();
 540     test_linear<double>();
 541     test_quadratic<double>();
 542     test_cubic<double>();
 543     test_quartic<double>();
 544     test_interpolation_condition<double>();
 545
 546     test_cardinal_constant<double>();
 547     test_cardinal_linear<double>();
 548     test_cardinal_quadratic<double>();
 549     test_cardinal_cubic<double>();
 550     test_cardinal_quartic<double>();
 551
 552     test_constant<long double>();
 553     test_linear<long double>();
 554     test_quadratic<long double>();
 555     test_cubic<long double>();
 556     test_quartic<long double>();
 557     test_interpolation_condition<long double>();
 558
 559     test_cardinal_constant<long double>();
 560     test_cardinal_linear<long double>();
 561     test_cardinal_quadratic<long double>();
 562     test_cardinal_cubic<long double>();
 563     test_cardinal_quartic<long double>();
 564
 565 #ifdef BOOST_HAS_FLOAT128
 566     test_constant<float128>();
 567     //test_linear<float128>();
 568     test_quadratic<float128>();
 569     test_cubic<float128>();
 570     test_quartic<float128>();
 571     test_interpolation_condition<float128>();
 572     test_cardinal_constant<float128>();
 573     test_cardinal_linear<float128>();
 574     test_cardinal_quadratic<float128>();
 575     test_cardinal_cubic<float128>();
 576     test_cardinal_quartic<float128>();
 577 #endif
 578
 579     return boost::math::test::report_errors();
 580 }