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 <cmath>
  10 #include <random>
  11 #include <boost/random/uniform_real_distribution.hpp>
  12 #include <boost/type_index.hpp>
  13 #include <boost/test/included/unit_test.hpp>
  14 #include <boost/test/tools/floating_point_comparison.hpp>
  15 #include <boost/math/interpolators/barycentric_rational.hpp>
  16 #include <boost/multiprecision/cpp_bin_float.hpp>
  17
  18 #ifdef BOOST_HAS_FLOAT128
  19 #include <boost/multiprecision/float128.hpp>
  20 #endif
  21
  22 using std::sqrt;
  23 using std::abs;
  24 using std::numeric_limits;
  25 using boost::multiprecision::cpp_bin_float_50;
  26
  27 template<class Real>
  28 void test_interpolation_condition()
  29 {
  30     std::cout << "Testing interpolation condition for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  31     std::mt19937 gen(4);
  32     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  33     std::vector<Real> x(100);
  34     std::vector<Real> y(100);
  35     x[0] = dis(gen);
  36     y[0] = dis(gen);
  37     for (size_t i = 1; i < x.size(); ++i)
  38     {
  39         x[i] = x[i-1] + dis(gen);
  40         y[i] = dis(gen);
  41     }
  42
  43     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size());
  44
  45     for (size_t i = 0; i < x.size(); ++i)
  46     {
  47         Real z = interpolator(x[i]);
  48         BOOST_CHECK_CLOSE(z, y[i], 100*numeric_limits<Real>::epsilon());
  49     }
  50
  51     // Make sure that the move constructor does the same thing:
  52     std::vector<Real> x_copy = x;
  53     std::vector<Real> y_copy = y;
  54     boost::math::barycentric_rational<Real> move_interpolator(std::move(x), std::move(y));
  55
  56     for (size_t i = 0; i < x_copy.size(); ++i)
  57     {
  58         Real z = move_interpolator(x_copy[i]);
  59         BOOST_CHECK_CLOSE(z, y_copy[i], 100*numeric_limits<Real>::epsilon());
  60     }
  61 }
  62
  63 template<class Real>
  64 void test_interpolation_condition_high_order()
  65 {
  66     std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  67     std::mt19937 gen(5);
  68     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  69     std::vector<Real> x(100);
  70     std::vector<Real> y(100);
  71     x[0] = dis(gen);
  72     y[0] = dis(gen);
  73     for (size_t i = 1; i < x.size(); ++i)
  74     {
  75         x[i] = x[i-1] + dis(gen);
  76         y[i] = dis(gen);
  77     }
  78
  79     // Order 5 approximation:
  80     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5);
  81
  82     for (size_t i = 0; i < x.size(); ++i)
  83     {
  84         Real z = interpolator(x[i]);
  85         BOOST_CHECK_CLOSE(z, y[i], 100*numeric_limits<Real>::epsilon());
  86     }
  87 }
  88
  89
  90 template<class Real>
  91 void test_constant()
  92 {
  93     std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  94
  95     std::mt19937 gen(6);
  96     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  97     std::vector<Real> x(100);
  98     std::vector<Real> y(100);
  99     Real constant = -8;
 100     x[0] = dis(gen);
 101     y[0] = constant;
 102     for (size_t i = 1; i < x.size(); ++i)
 103     {
 104         x[i] = x[i-1] + dis(gen);
 105         y[i] = y[0];
 106     }
 107
 108     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size());
 109
 110     for (size_t i = 0; i < x.size(); ++i)
 111     {
 112         // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
 113         Real t = x[i] + dis(gen);
 114         Real z = interpolator(t);
 115         BOOST_CHECK_CLOSE(z, constant, 100*sqrt(numeric_limits<Real>::epsilon()));
 116         BOOST_CHECK_SMALL(interpolator.prime(t), sqrt(numeric_limits<Real>::epsilon()));
 117     }
 118 }
 119
 120 template<class Real>
 121 void test_constant_high_order()
 122 {
 123     std::cout << "Testing that constants are interpolated correctly in high order using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 124
 125     std::mt19937 gen(7);
 126     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 127     std::vector<Real> x(100);
 128     std::vector<Real> y(100);
 129     Real constant = 5;
 130     x[0] = dis(gen);
 131     y[0] = constant;
 132     for (size_t i = 1; i < x.size(); ++i)
 133     {
 134         x[i] = x[i-1] + dis(gen);
 135         y[i] = y[0];
 136     }
 137
 138     // Set interpolation order to 7:
 139     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 7);
 140
 141     for (size_t i = 0; i < x.size(); ++i)
 142     {
 143         Real t = x[i] + dis(gen);
 144         Real z = interpolator(t);
 145         BOOST_CHECK_CLOSE(z, constant, 1000*sqrt(numeric_limits<Real>::epsilon()));
 146         BOOST_CHECK_SMALL(interpolator.prime(t), 100*sqrt(numeric_limits<Real>::epsilon()));
 147     }
 148 }
 149
 150
 151 template<class Real>
 152 void test_runge()
 153 {
 154     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";
 155
 156     std::mt19937 gen(8);
 157     boost::random::uniform_real_distribution<Real> dis(0.005f, 0.01f);
 158     std::vector<Real> x(100);
 159     std::vector<Real> y(100);
 160     x[0] = -2;
 161     y[0] = 1/(1+25*x[0]*x[0]);
 162     for (size_t i = 1; i < x.size(); ++i)
 163     {
 164         x[i] = x[i-1] + dis(gen);
 165         y[i] = 1/(1+25*x[i]*x[i]);
 166     }
 167
 168     boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5);
 169
 170     for (size_t i = 0; i < x.size(); ++i)
 171     {
 172         Real t = x[i];
 173         Real z = interpolator(t);
 174         BOOST_CHECK_CLOSE(z, y[i], 0.03);
 175         Real z_prime = interpolator.prime(t);
 176         Real num = -50*t;
 177         Real denom = (1+25*t*t)*(1+25*t*t);
 178         if (abs(num/denom) > 0.00001)
 179         {
 180             BOOST_CHECK_CLOSE_FRACTION(z_prime, num/denom, 0.03);
 181         }
 182     }
 183
 184
 185     Real tol = 0.0001;
 186     for (size_t i = 0; i < x.size(); ++i)
 187     {
 188         Real t = x[i] + dis(gen);
 189         Real z = interpolator(t);
 190         BOOST_CHECK_CLOSE(z, 1/(1+25*t*t), tol);
 191         Real z_prime = interpolator.prime(t);
 192         Real num = -50*t;
 193         Real denom = (1+25*t*t)*(1+25*t*t);
 194         Real runge_prime = num/denom;
 195
 196         if (abs(runge_prime) > 0 && abs(z_prime - runge_prime)/abs(runge_prime) > tol)
 197         {
 198             std::cout << "Error too high for t = " << t << " which is a distance " << t - x[i] << " from node " << i << "/" << x.size() << " associated with data (" << x[i] << ", " << y[i] << ")\n";
 199             BOOST_CHECK_CLOSE_FRACTION(z_prime, runge_prime, tol);
 200         }
 201     }
 202 }
 203
 204 template<class Real>
 205 void test_weights()
 206 {
 207     std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 208
 209     std::mt19937 gen(9);
 210     boost::random::uniform_real_distribution<Real> dis(0.005, 0.01);
 211     std::vector<Real> x(100);
 212     std::vector<Real> y(100);
 213     x[0] = -2;
 214     y[0] = 1/(1+25*x[0]*x[0]);
 215     for (size_t i = 1; i < x.size(); ++i)
 216     {
 217         x[i] = x[i-1] + dis(gen);
 218         y[i] = 1/(1+25*x[i]*x[i]);
 219     }
 220
 221     boost::math::detail::barycentric_rational_imp<Real> interpolator(x.data(), x.data() + x.size(), y.data(), 0);
 222
 223     for (size_t i = 0; i < x.size(); ++i)
 224     {
 225         Real w = interpolator.weight(i);
 226         if (i % 2 == 0)
 227         {
 228             BOOST_CHECK_CLOSE(w, 1, 0.00001);
 229         }
 230         else
 231         {
 232             BOOST_CHECK_CLOSE(w, -1, 0.00001);
 233         }
 234     }
 235
 236     // d = 1:
 237     interpolator = boost::math::detail::barycentric_rational_imp<Real>(x.data(), x.data() + x.size(), y.data(), 1);
 238
 239     for (size_t i = 1; i < x.size() -1; ++i)
 240     {
 241         Real w = interpolator.weight(i);
 242         Real w_expect = 1/(x[i] - x[i - 1]) + 1/(x[i+1] - x[i]);
 243         if (i % 2 == 0)
 244         {
 245             BOOST_CHECK_CLOSE(w, -w_expect, 0.00001);
 246         }
 247         else
 248         {
 249             BOOST_CHECK_CLOSE(w, w_expect, 0.00001);
 250         }
 251     }
 252
 253 }
 254
 255
 256 BOOST_AUTO_TEST_CASE(barycentric_rational)
 257 {
 258     // The tests took too long at the higher precisions.
 259     // They still pass, but the CI system is starting to time out,
 260     // so I figured it'd be polite to comment out the most expensive tests.
 261     test_weights<double>();
 262
 263     test_constant<float>();
 264     //test_constant<double>();
 265     test_constant<long double>();
 266     //test_constant<cpp_bin_float_50>();
 267
 268     //test_constant_high_order<float>();
 269     test_constant_high_order<double>();
 270     //test_constant_high_order<long double>();
 271     //test_constant_high_order<cpp_bin_float_50>();
 272
 273     test_interpolation_condition<float>();
 274     test_interpolation_condition<double>();
 275     //test_interpolation_condition<long double>();
 276     //test_interpolation_condition<cpp_bin_float_50>();
 277
 278     //test_interpolation_condition_high_order<float>();
 279     test_interpolation_condition_high_order<double>();
 280     //test_interpolation_condition_high_order<long double>();
 281     //test_interpolation_condition_high_order<cpp_bin_float_50>();
 282
 283     test_runge<double>();
 284     //test_runge<long double>();
 285     //test_runge<cpp_bin_float_50>();
 286
 287 #ifdef BOOST_HAS_FLOAT128
 288     //test_interpolation_condition<boost::multiprecision::float128>();
 289     //test_constant<boost::multiprecision::float128>();
 290     //test_constant_high_order<boost::multiprecision::float128>();
 291     //test_interpolation_condition_high_order<boost::multiprecision::float128>();
 292     //test_runge<boost::multiprecision::float128>();
 293 #endif
 294
 295 }