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

   1 // Copyright Nick Thompson, 2017
   2 // Use, modification and distribution are subject to the
   3 // Boost Software License, Version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #define BOOST_TEST_MODULE barycentric_rational
   8
   9 #include <random>
  10 #include <boost/random/uniform_real_distribution.hpp>
  11 #include <boost/type_index.hpp>
  12 #include <boost/test/included/unit_test.hpp>
  13 #include <boost/test/floating_point_comparison.hpp>
  14 #include <boost/math/interpolators/barycentric_rational.hpp>
  15 #include <boost/multiprecision/cpp_bin_float.hpp>
  16
  17 #ifdef BOOST_HAS_FLOAT128
  18 #include <boost/multiprecision/float128.hpp>
  19 #endif
  20
  21 using boost::multiprecision::cpp_bin_float_50;
  22
  23 template<class Real>
  24 void test_interpolation_condition()
  25 {
  26     std::cout << "Testing interpolation condition for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  27     std::mt19937 gen(4);
  28     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  29     std::vector<Real> x(500);
  30     std::vector<Real> y(500);
  31     x[0] = dis(gen);
  32     y[0] = dis(gen);
  33     for (size_t i = 1; i < x.size(); ++i)
  34     {
  35         x[i] = x[i-1] + dis(gen);
  36         y[i] = dis(gen);
  37     }
  38
  39     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size());
  40
  41     for (size_t i = 0; i < x.size(); ++i)
  42     {
  43         auto z = interpolator(x[i]);
  44         BOOST_CHECK_CLOSE(z, y[i], 100*std::numeric_limits<Real>::epsilon());
  45     }
  46 }
  47
  48 template<class Real>
  49 void test_interpolation_condition_high_order()
  50 {
  51     std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  52     std::mt19937 gen(5);
  53     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  54     std::vector<Real> x(500);
  55     std::vector<Real> y(500);
  56     x[0] = dis(gen);
  57     y[0] = dis(gen);
  58     for (size_t i = 1; i < x.size(); ++i)
  59     {
  60         x[i] = x[i-1] + dis(gen);
  61         y[i] = dis(gen);
  62     }
  63
  64     // Order 5 approximation:
  65     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5);
  66
  67     for (size_t i = 0; i < x.size(); ++i)
  68     {
  69         auto z = interpolator(x[i]);
  70         BOOST_CHECK_CLOSE(z, y[i], 100*std::numeric_limits<Real>::epsilon());
  71     }
  72 }
  73
  74
  75 template<class Real>
  76 void test_constant()
  77 {
  78     std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  79
  80     std::mt19937 gen(6);
  81     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  82     std::vector<Real> x(500);
  83     std::vector<Real> y(500);
  84     Real constant = -8;
  85     x[0] = dis(gen);
  86     y[0] = constant;
  87     for (size_t i = 1; i < x.size(); ++i)
  88     {
  89         x[i] = x[i-1] + dis(gen);
  90         y[i] = y[0];
  91     }
  92
  93     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size());
  94
  95     for (size_t i = 0; i < x.size(); ++i)
  96     {
  97         // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
  98         auto z = interpolator(x[i] + dis(gen));
  99         BOOST_CHECK_CLOSE(z, constant, 100*sqrt(std::numeric_limits<Real>::epsilon()));
 100     }
 101 }
 102
 103 template<class Real>
 104 void test_constant_high_order()
 105 {
 106     std::cout << "Testing that constants are interpolated correctly in high order using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 107
 108     std::mt19937 gen(7);
 109     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 110     std::vector<Real> x(500);
 111     std::vector<Real> y(500);
 112     Real constant = 5;
 113     x[0] = dis(gen);
 114     y[0] = constant;
 115     for (size_t i = 1; i < x.size(); ++i)
 116     {
 117         x[i] = x[i-1] + dis(gen);
 118         y[i] = y[0];
 119     }
 120
 121     // Set interpolation order to 7:
 122     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 7);
 123
 124     for (size_t i = 0; i < x.size(); ++i)
 125     {
 126         auto z = interpolator(x[i] + dis(gen));
 127         BOOST_CHECK_CLOSE(z, constant, 1000*sqrt(std::numeric_limits<Real>::epsilon()));
 128     }
 129 }
 130
 131
 132 template<class Real>
 133 void test_runge()
 134 {
 135     std::cout << "Testing interpolation of Runge's 1/(1+25x^2) function using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 136
 137     std::mt19937 gen(8);
 138     boost::random::uniform_real_distribution<Real> dis(0.005f, 0.01f);
 139     std::vector<Real> x(500);
 140     std::vector<Real> y(500);
 141     x[0] = -2;
 142     y[0] = 1/(1+25*x[0]*x[0]);
 143     for (size_t i = 1; i < x.size(); ++i)
 144     {
 145         x[i] = x[i-1] + dis(gen);
 146         y[i] = 1/(1+25*x[i]*x[i]);
 147     }
 148
 149     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5);
 150
 151     for (size_t i = 0; i < x.size(); ++i)
 152     {
 153         auto t = x[i] + dis(gen);
 154         auto z = interpolator(t);
 155         BOOST_CHECK_CLOSE(z, 1/(1+25*t*t), 0.03);
 156     }
 157 }
 158
 159 template<class Real>
 160 void test_weights()
 161 {
 162     std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 163
 164     std::mt19937 gen(9);
 165     boost::random::uniform_real_distribution<Real> dis(0.005, 0.01);
 166     std::vector<Real> x(500);
 167     std::vector<Real> y(500);
 168     x[0] = -2;
 169     y[0] = 1/(1+25*x[0]*x[0]);
 170     for (size_t i = 1; i < x.size(); ++i)
 171     {
 172         x[i] = x[i-1] + dis(gen);
 173         y[i] = 1/(1+25*x[i]*x[i]);
 174     }
 175
 176     boost::math::detail::barycentric_rational_imp<Real> interpolator(x.data(), x.data() + x.size(), y.data(), 0);
 177
 178     for (size_t i = 0; i < x.size(); ++i)
 179     {
 180         auto w = interpolator.weight(i);
 181         if (i % 2 == 0)
 182         {
 183             BOOST_CHECK_CLOSE(w, 1, 0.00001);
 184         }
 185         else
 186         {
 187             BOOST_CHECK_CLOSE(w, -1, 0.00001);
 188         }
 189     }
 190
 191     // d = 1:
 192     interpolator = boost::math::detail::barycentric_rational_imp<Real>(x.data(), x.data() + x.size(), y.data(), 1);
 193
 194     for (size_t i = 1; i < x.size() -1; ++i)
 195     {
 196         auto w = interpolator.weight(i);
 197         auto w_expect = 1/(x[i] - x[i - 1]) + 1/(x[i+1] - x[i]);
 198         if (i % 2 == 0)
 199         {
 200             BOOST_CHECK_CLOSE(w, -w_expect, 0.00001);
 201         }
 202         else
 203         {
 204             BOOST_CHECK_CLOSE(w, w_expect, 0.00001);
 205         }
 206     }
 207
 208 }
 209
 210
 211 BOOST_AUTO_TEST_CASE(barycentric_rational)
 212 {
 213     test_weights<double>();
 214     test_constant<float>();
 215     test_constant<double>();
 216     test_constant<long double>();
 217     test_constant<cpp_bin_float_50>();
 218
 219     test_constant_high_order<float>();
 220     test_constant_high_order<double>();
 221     test_constant_high_order<long double>();
 222     test_constant_high_order<cpp_bin_float_50>();
 223
 224     test_interpolation_condition<float>();
 225     test_interpolation_condition<double>();
 226     test_interpolation_condition<long double>();
 227     test_interpolation_condition<cpp_bin_float_50>();
 228
 229     test_interpolation_condition_high_order<float>();
 230     test_interpolation_condition_high_order<double>();
 231     test_interpolation_condition_high_order<long double>();
 232     test_interpolation_condition_high_order<cpp_bin_float_50>();
 233
 234     test_runge<float>();
 235     test_runge<double>();
 236     test_runge<long double>();
 237     test_runge<cpp_bin_float_50>();
 238
 239 #ifdef BOOST_HAS_FLOAT128
 240     test_interpolation_condition<boost::multiprecision::float128>();
 241         test_constant<boost::multiprecision::float128>();
 242         test_constant_high_order<boost::multiprecision::float128>();
 243         test_interpolation_condition_high_order<boost::multiprecision::float128>();
 244         test_runge<boost::multiprecision::float128>();
 245 #endif
 246
 247 }