ceph/src/boost/libs/math/test/whittaker_shannon_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 <random>
  12 #include <boost/core/demangle.hpp>
  13 #include <boost/math/interpolators/whittaker_shannon.hpp>
  14 #ifdef BOOST_HAS_FLOAT128
  15 #include <boost/multiprecision/float128.hpp>
  16 using boost::multiprecision::float128;
  17 #endif
  18
  19 using boost::math::interpolators::whittaker_shannon;
  20
  21 template<class Real>
  22 void test_trivial()
  23 {
  24     Real t0 = 0;
  25     Real h = Real(1)/Real(16);
  26     std::vector<Real> v{1.5};
  27     std::vector<Real> v_copy = v;
  28     auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
  29
  30
  31     Real expected = 0;
  32     if(!CHECK_MOLLIFIED_CLOSE(expected, ws.prime(0), 10*std::numeric_limits<Real>::epsilon())) {
  33         std::cerr << "  Problem occurred at abscissa " << 0 << "\n";
  34     }
  35
  36     expected = -v_copy[0]/h;
  37     if(!CHECK_MOLLIFIED_CLOSE(expected, ws.prime(h), 10*std::numeric_limits<Real>::epsilon())) {
  38         std::cerr << "  Problem occurred at abscissa " << 0 << "\n";
  39     }
  40 }
  41
  42 template<class Real>
  43 void test_knots()
  44 {
  45     Real t0 = 0;
  46     Real h = Real(1)/Real(16);
  47     size_t n = 512;
  48     std::vector<Real> v(n);
  49     std::mt19937 gen(323723);
  50     std::uniform_real_distribution<Real> dis(1.0, 2.0);
  51
  52     for(size_t i = 0;  i < n; ++i) {
  53       v[i] = static_cast<Real>(dis(gen));
  54     }
  55     auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
  56
  57     size_t i = 0;
  58     while (i < n) {
  59       Real t = t0 + i*h;
  60       Real expected = ws[i];
  61       Real computed = ws(t);
  62       CHECK_ULP_CLOSE(expected, computed, 16);
  63       ++i;
  64     }
  65 }
  66
  67 template<class Real>
  68 void test_bump()
  69 {
  70     using std::exp;
  71     using std::abs;
  72     using std::sqrt;
  73     auto bump = [](Real x) { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };
  74
  75     auto bump_prime = [&bump](Real x) { Real z = 1-x*x; return -2*x*bump(x)/(z*z); };
  76
  77     Real t0 = -1;
  78     size_t n = 2049;
  79     Real h = Real(2)/Real(n-1);
  80
  81     std::vector<Real> v(n);
  82     for(size_t i = 0; i < n; ++i) {
  83         Real t = t0 + i*h;
  84         v[i] = bump(t);
  85     }
  86
  87
  88     std::vector<Real> v_copy = v;
  89     auto ws = whittaker_shannon<decltype(v)>(std::move(v), t0, h);
  90
  91     // Test the knots:
  92     for(size_t i = v_copy.size()/4; i < 3*v_copy.size()/4; ++i) {
  93         Real t = t0 + i*h;
  94         Real expected = v_copy[i];
  95         Real computed = ws(t);
  96         if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 10*std::numeric_limits<Real>::epsilon())) {
  97             std::cerr << "  Problem occurred at abscissa " << t << "\n";
  98         }
  99
 100         Real expected_prime = bump_prime(t);
 101         Real computed_prime = ws.prime(t);
 102         if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, 1000*std::numeric_limits<Real>::epsilon())) {
 103             std::cerr << "  Problem occurred at abscissa " << t << "\n";
 104         }
 105
 106     }
 107
 108     std::mt19937 gen(323723);
 109     std::uniform_real_distribution<long double> dis(-0.85, 0.85);
 110
 111     size_t i = 0;
 112     while (i++ < 1000)
 113     {
 114         Real t = static_cast<Real>(dis(gen));
 115         Real expected = bump(t);
 116         Real computed = ws(t);
 117         if(!CHECK_MOLLIFIED_CLOSE(expected, computed, 10*std::numeric_limits<Real>::epsilon())) {
 118             std::cerr << "  Problem occurred at abscissa " << t << "\n";
 119         }
 120
 121         Real expected_prime = bump_prime(t);
 122         Real computed_prime = ws.prime(t);
 123         if(!CHECK_MOLLIFIED_CLOSE(expected_prime, computed_prime, sqrt(std::numeric_limits<Real>::epsilon()))) {
 124             std::cerr << "  Problem occurred at abscissa " << t << "\n";
 125         }
 126     }
 127 }
 128
 129
 130 int main()
 131 {
 132     test_knots<float>();
 133     test_knots<double>();
 134     test_knots<long double>();
 135
 136     test_bump<double>();
 137     test_bump<long double>();
 138
 139     test_trivial<float>();
 140     test_trivial<double>();
 141     return boost::math::test::report_errors();
 142 }