ceph/src/boost/libs/math/test/cardinal_quadratic_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_quadratic_b_spline.hpp>
  12 using boost::math::interpolators::cardinal_quadratic_b_spline;
  13
  14 template<class Real>
  15 void test_constant()
  16 {
  17     Real c = 7.2;
  18     Real t0 = 0;
  19     Real h = Real(1)/Real(16);
  20     size_t n = 512;
  21     std::vector<Real> v(n, c);
  22     auto qbs = cardinal_quadratic_b_spline<Real>(v.data(), v.size(), t0, h);
  23
  24     size_t i = 0;
  25     while (i < n) {
  26       Real t = t0 + i*h;
  27       CHECK_ULP_CLOSE(c, qbs(t), 2);
  28       CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 100*std::numeric_limits<Real>::epsilon());
  29       ++i;
  30     }
  31
  32     i = 0;
  33     while (i < n) {
  34       Real t = t0 + i*h + h/2;
  35       CHECK_ULP_CLOSE(c, qbs(t), 2);
  36       CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 300*std::numeric_limits<Real>::epsilon());
  37       t = t0 + i*h + h/4;
  38       CHECK_ULP_CLOSE(c, qbs(t), 2);
  39       CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 150*std::numeric_limits<Real>::epsilon());
  40       ++i;
  41     }
  42 }
  43
  44 template<class Real>
  45 void test_linear()
  46 {
  47     Real m = 8.3;
  48     Real b = 7.2;
  49     Real t0 = 0;
  50     Real h = Real(1)/Real(16);
  51     size_t n = 512;
  52     std::vector<Real> y(n);
  53     for (size_t i = 0; i < n; ++i) {
  54       Real t = i*h;
  55       y[i] = m*t + b;
  56     }
  57     auto qbs = cardinal_quadratic_b_spline<Real>(y.data(), y.size(), t0, h);
  58
  59     size_t i = 0;
  60     while (i < n) {
  61       Real t = t0 + i*h;
  62       CHECK_ULP_CLOSE(m*t+b, qbs(t), 2);
  63       CHECK_ULP_CLOSE(m, qbs.prime(t), 820);
  64       ++i;
  65     }
  66
  67     i = 0;
  68     while (i < n) {
  69       Real t = t0 + i*h + h/2;
  70       CHECK_ULP_CLOSE(m*t+b, qbs(t), 2);
  71       CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
  72       t = t0 + i*h + h/4;
  73       CHECK_ULP_CLOSE(m*t+b, qbs(t), 3);
  74       CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits<Real>::epsilon());
  75       ++i;
  76     }
  77 }
  78
  79 template<class Real>
  80 void test_quadratic()
  81 {
  82     Real a = 8.2;
  83     Real b = 7.2;
  84     Real c = -9.2;
  85     Real t0 = 0;
  86     Real h = Real(1)/Real(16);
  87     size_t n = 513;
  88     std::vector<Real> y(n);
  89     for (size_t i = 0; i < n; ++i) {
  90       Real t = i*h;
  91       y[i] = a*t*t + b*t + c;
  92     }
  93     Real t_max = t0 + (n-1)*h;
  94     auto qbs = cardinal_quadratic_b_spline<Real>(y, t0, h, b, 2*a*t_max + b);
  95
  96     size_t i = 0;
  97     while (i < n) {
  98       Real t = t0 + i*h;
  99       CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 2);
 100       ++i;
 101     }
 102
 103     i = 0;
 104     while (i < n) {
 105       Real t = t0 + i*h + h/2;
 106       CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 47);
 107
 108       t = t0 + i*h + h/4;
 109       if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 104)) {
 110           std::cerr << "  Problem abscissa t = " << t << "\n";
 111       }
 112       ++i;
 113     }
 114 }
 115
 116 int main()
 117 {
 118     test_constant<float>();
 119     test_constant<double>();
 120     test_constant<long double>();
 121
 122     test_linear<float>();
 123     test_linear<double>();
 124     test_linear<long double>();
 125
 126     test_quadratic<double>();
 127     test_quadratic<long double>();
 128
 129     return boost::math::test::report_errors();
 130 }