libs/math/test/test_cubic_b_spline.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 test_cubic_b_spline
   8
   9 #include <random>
  10 #include <functional>
  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/floating_point_comparison.hpp>
  15 #include <boost/math/interpolators/cubic_b_spline.hpp>
  16 #include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
  17 #include <boost/multiprecision/cpp_bin_float.hpp>
  18
  19 using boost::multiprecision::cpp_bin_float_50;
  20 using boost::math::constants::third;
  21 using boost::math::constants::half;
  22
  23 template<class Real>
  24 void test_b3_spline()
  25 {
  26     std::cout << "Testing evaluation of spline basis functions on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  27     // Outside the support:
  28     Real eps = std::numeric_limits<Real>::epsilon();
  29     BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(2.5), (Real) 0);
  30     BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(-2.5), (Real) 0);
  31     BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(2.5), (Real) 0);
  32     BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(-2.5), (Real) 0);
  33
  34     // On the boundary of support:
  35     BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(2), (Real) 0);
  36     BOOST_CHECK_SMALL(boost::math::detail::b3_spline<Real>(-2), (Real) 0);
  37     BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(2), (Real) 0);
  38     BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(-2), (Real) 0);
  39
  40     // Special values:
  41     BOOST_CHECK_CLOSE(boost::math::detail::b3_spline<Real>(-1), third<Real>()*half<Real>(), eps);
  42     BOOST_CHECK_CLOSE(boost::math::detail::b3_spline<Real>( 1), third<Real>()*half<Real>(), eps);
  43     BOOST_CHECK_CLOSE(boost::math::detail::b3_spline<Real>(0), 2*third<Real>(), eps);
  44
  45     BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_prime<Real>(-1), half<Real>(), eps);
  46     BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_prime<Real>( 1), -half<Real>(), eps);
  47     BOOST_CHECK_SMALL(boost::math::detail::b3_spline_prime<Real>(0), eps);
  48
  49     // Properties: B3 is an even function, B3' is an odd function.
  50     for (size_t i = 1; i < 200; ++i)
  51     {
  52         Real arg = i*0.01;
  53         BOOST_CHECK_CLOSE(boost::math::detail::b3_spline<Real>(arg), boost::math::detail::b3_spline<Real>(arg), eps);
  54         BOOST_CHECK_CLOSE(boost::math::detail::b3_spline_prime<Real>(-arg), -boost::math::detail::b3_spline_prime<Real>(arg), eps);
  55     }
  56
  57 }
  58 /*
  59  * This test ensures that the interpolant s(x_j) = f(x_j) at all grid points.
  60  */
  61 template<class Real>
  62 void test_interpolation_condition()
  63 {
  64     using std::sqrt;
  65     std::cout << "Testing interpolation condition for cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  66     std::random_device rd;
  67     std::mt19937 gen(rd());
  68     boost::random::uniform_real_distribution<Real> dis(1, 10);
  69     std::vector<Real> v(5000);
  70     for (size_t i = 0; i < v.size(); ++i)
  71     {
  72         v[i] = dis(gen);
  73     }
  74
  75     Real step = 0.01;
  76     Real a = 5;
  77     boost::math::cubic_b_spline<Real> spline(v.data(), v.size(), a, step);
  78
  79     for (size_t i = 0; i < v.size(); ++i)
  80     {
  81         Real y = spline(i*step + a);
  82         // This seems like a very large tolerance, but I don't know of any other interpolators
  83         // that will be able to do much better on random data.
  84         BOOST_CHECK_CLOSE(y, v[i], 10000*sqrt(std::numeric_limits<Real>::epsilon()));
  85     }
  86
  87 }
  88
  89
  90 template<class Real>
  91 void test_constant_function()
  92 {
  93     std::cout << "Testing that constants are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  94     std::vector<Real> v(500);
  95     Real constant = 50.2;
  96     for (size_t i = 0; i < v.size(); ++i)
  97     {
  98         v[i] = 50.2;
  99     }
 100
 101     Real step = 0.02;
 102     Real a = 5;
 103     boost::math::cubic_b_spline<Real> spline(v.data(), v.size(), a, step);
 104
 105     for (size_t i = 0; i < v.size(); ++i)
 106     {
 107         // Do not test at interpolation point; we already know it works there:
 108         Real y = spline(i*step + a + 0.001);
 109         BOOST_CHECK_CLOSE(y, constant, 10*std::numeric_limits<Real>::epsilon());
 110         Real y_prime = spline.prime(i*step + a + 0.002);
 111         BOOST_CHECK_SMALL(y_prime, 5000*std::numeric_limits<Real>::epsilon());
 112     }
 113
 114     // Test that correctly specified left and right-derivatives work properly:
 115     spline = boost::math::cubic_b_spline<Real>(v.data(), v.size(), a, step, 0, 0);
 116
 117     for (size_t i = 0; i < v.size(); ++i)
 118     {
 119         Real y = spline(i*step + a + 0.002);
 120         BOOST_CHECK_CLOSE(y, constant, std::numeric_limits<Real>::epsilon());
 121         Real y_prime = spline.prime(i*step + a + 0.002);
 122         BOOST_CHECK_SMALL(y_prime, std::numeric_limits<Real>::epsilon());
 123     }
 124
 125     //
 126     // Again with iterator constructor:
 127     //
 128     boost::math::cubic_b_spline<Real> spline2(v.begin(), v.end(), a, step);
 129
 130     for (size_t i = 0; i < v.size(); ++i)
 131     {
 132        // Do not test at interpolation point; we already know it works there:
 133        Real y = spline2(i*step + a + 0.001);
 134        BOOST_CHECK_CLOSE(y, constant, 10 * std::numeric_limits<Real>::epsilon());
 135        Real y_prime = spline2.prime(i*step + a + 0.002);
 136        BOOST_CHECK_SMALL(y_prime, 5000 * std::numeric_limits<Real>::epsilon());
 137     }
 138
 139     // Test that correctly specified left and right-derivatives work properly:
 140     spline2 = boost::math::cubic_b_spline<Real>(v.begin(), v.end(), a, step, 0, 0);
 141
 142     for (size_t i = 0; i < v.size(); ++i)
 143     {
 144        Real y = spline2(i*step + a + 0.002);
 145        BOOST_CHECK_CLOSE(y, constant, std::numeric_limits<Real>::epsilon());
 146        Real y_prime = spline2.prime(i*step + a + 0.002);
 147        BOOST_CHECK_SMALL(y_prime, std::numeric_limits<Real>::epsilon());
 148     }
 149 }
 150
 151
 152 template<class Real>
 153 void test_affine_function()
 154 {
 155     using std::sqrt;
 156     std::cout << "Testing that affine functions are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 157     std::vector<Real> v(500);
 158     Real a = 10;
 159     Real b = 8;
 160     Real step = 0.005;
 161
 162     auto f = [a, b](Real x) { return a*x + b; };
 163     for (size_t i = 0; i < v.size(); ++i)
 164     {
 165         v[i] = f(i*step);
 166     }
 167
 168     boost::math::cubic_b_spline<Real> spline(v.data(), v.size(), 0, step);
 169
 170     for (size_t i = 0; i < v.size() - 1; ++i)
 171     {
 172         Real arg = i*step + 0.0001;
 173         Real y = spline(arg);
 174         BOOST_CHECK_CLOSE(y, f(arg), sqrt(std::numeric_limits<Real>::epsilon()));
 175         Real y_prime = spline.prime(arg);
 176         BOOST_CHECK_CLOSE(y_prime, a, 100*sqrt(std::numeric_limits<Real>::epsilon()));
 177     }
 178
 179     // Test that correctly specified left and right-derivatives work properly:
 180     spline = boost::math::cubic_b_spline<Real>(v.data(), v.size(), 0, step, a, a);
 181
 182     for (size_t i = 0; i < v.size() - 1; ++i)
 183     {
 184         Real arg = i*step + 0.0001;
 185         Real y = spline(arg);
 186         BOOST_CHECK_CLOSE(y, f(arg), sqrt(std::numeric_limits<Real>::epsilon()));
 187         Real y_prime = spline.prime(arg);
 188         BOOST_CHECK_CLOSE(y_prime, a, 100*sqrt(std::numeric_limits<Real>::epsilon()));
 189     }
 190 }
 191
 192
 193 template<class Real>
 194 void test_quadratic_function()
 195 {
 196     using std::sqrt;
 197     std::cout << "Testing that quadratic functions are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 198     std::vector<Real> v(500);
 199     Real a = 1.2;
 200     Real b = -3.4;
 201     Real c = -8.6;
 202     Real step = 0.01;
 203
 204     auto f = [a, b, c](Real x) { return a*x*x + b*x + c; };
 205     for (size_t i = 0; i < v.size(); ++i)
 206     {
 207         v[i] = f(i*step);
 208     }
 209
 210     boost::math::cubic_b_spline<Real> spline(v.data(), v.size(), 0, step);
 211
 212     for (size_t i = 0; i < v.size() -1; ++i)
 213     {
 214         Real arg = i*step + 0.001;
 215         Real y = spline(arg);
 216         BOOST_CHECK_CLOSE(y, f(arg), 0.1);
 217         Real y_prime = spline.prime(arg);
 218         BOOST_CHECK_CLOSE(y_prime, 2*a*arg + b, 2.0);
 219     }
 220 }
 221
 222
 223 template<class Real>
 224 void test_trig_function()
 225 {
 226     std::cout << "Testing that sine functions are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 227     std::mt19937 gen;
 228     std::vector<Real> v(500);
 229     Real x0 = 1;
 230     Real step = 0.125;
 231
 232     for (size_t i = 0; i < v.size(); ++i)
 233     {
 234         v[i] = sin(x0 + step * i);
 235     }
 236
 237     boost::math::cubic_b_spline<Real> spline(v.data(), v.size(), x0, step);
 238
 239     boost::random::uniform_real_distribution<Real> absissa(x0, x0 + 499 * step);
 240
 241     for (size_t i = 0; i < v.size(); ++i)
 242     {
 243         Real x = absissa(gen);
 244         Real y = spline(x);
 245         BOOST_CHECK_CLOSE(y, sin(x), 1.0);
 246         auto y_prime = spline.prime(x);
 247         BOOST_CHECK_CLOSE(y_prime, cos(x), 2.0);
 248     }
 249 }
 250
 251 template<class Real>
 252 void test_copy_move()
 253 {
 254     std::cout << "Testing that copy/move operation succeed on cubic b spline\n";
 255     std::vector<Real> v(500);
 256     Real x0 = 1;
 257     Real step = 0.125;
 258
 259     for (size_t i = 0; i < v.size(); ++i)
 260     {
 261         v[i] = sin(x0 + step * i);
 262     }
 263
 264     boost::math::cubic_b_spline<Real> spline(v.data(), v.size(), x0, step);
 265
 266
 267     // Default constructor should compile so that splines can be member variables:
 268     boost::math::cubic_b_spline<Real> d;
 269     d = boost::math::cubic_b_spline<Real>(v.data(), v.size(), x0, step);
 270     BOOST_CHECK_CLOSE(d(x0), sin(x0), 0.01);
 271     // Passing to lambda should compile:
 272     auto f = [=](Real x) { return d(x); };
 273     // Make sure this variable is used.
 274     BOOST_CHECK_CLOSE(f(x0), sin(x0), 0.01);
 275
 276     // Move operations should compile.
 277     auto s = std::move(spline);
 278
 279     // Copy operations should compile:
 280     boost::math::cubic_b_spline<Real> c = d;
 281     BOOST_CHECK_CLOSE(c(x0), sin(x0), 0.01);
 282
 283     // Test with std::bind:
 284     auto h = std::bind(&boost::math::cubic_b_spline<double>::operator(), &s, std::placeholders::_1);
 285     BOOST_CHECK_CLOSE(h(x0), sin(x0), 0.01);
 286 }
 287
 288 template<class Real>
 289 void test_outside_interval()
 290 {
 291     std::cout << "Testing that the spline can be evaluated outside the interpolation interval\n";
 292     std::vector<Real> v(400);
 293     Real x0 = 1;
 294     Real step = 0.125;
 295
 296     for (size_t i = 0; i < v.size(); ++i)
 297     {
 298         v[i] = sin(x0 + step * i);
 299     }
 300
 301     boost::math::cubic_b_spline<Real> spline(v.data(), v.size(), x0, step);
 302
 303     // There's no test here; it simply does it's best to be an extrapolator.
 304     //
 305     std::ostream cnull(0);
 306     cnull << spline(0);
 307     cnull << spline(2000);
 308 }
 309
 310 BOOST_AUTO_TEST_CASE(test_cubic_b_spline)
 311 {
 312     test_b3_spline<float>();
 313     test_b3_spline<double>();
 314     test_b3_spline<long double>();
 315     test_b3_spline<cpp_bin_float_50>();
 316
 317     test_interpolation_condition<float>();
 318     test_interpolation_condition<double>();
 319     test_interpolation_condition<long double>();
 320     test_interpolation_condition<cpp_bin_float_50>();
 321
 322     test_constant_function<float>();
 323     test_constant_function<double>();
 324     test_constant_function<long double>();
 325     test_constant_function<cpp_bin_float_50>();
 326
 327     test_affine_function<float>();
 328     test_affine_function<double>();
 329     test_affine_function<long double>();
 330     test_affine_function<cpp_bin_float_50>();
 331
 332     test_quadratic_function<float>();
 333     test_quadratic_function<double>();
 334     test_quadratic_function<long double>();
 335     test_affine_function<cpp_bin_float_50>();
 336
 337     test_trig_function<float>();
 338     test_trig_function<double>();
 339     test_trig_function<long double>();
 340     test_trig_function<cpp_bin_float_50>();
 341
 342     test_copy_move<double>();
 343     test_outside_interval<double>();
 344 }