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

   1 // Copyright Nick Thompson, 2019
   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 vector_barycentric_rational
   8
   9 #include <cmath>
  10 #include <random>
  11 #include <array>
  12 #include <Eigen/Dense>
  13 #include <boost/numeric/ublas/vector.hpp>
  14 #include <boost/random/uniform_real_distribution.hpp>
  15 #include <boost/type_index.hpp>
  16 #include <boost/test/included/unit_test.hpp>
  17 #include <boost/test/tools/floating_point_comparison.hpp>
  18 #include <boost/math/interpolators/barycentric_rational.hpp>
  19 #include <boost/math/interpolators/vector_barycentric_rational.hpp>
  20
  21 using std::sqrt;
  22 using std::abs;
  23 using std::numeric_limits;
  24
  25 template<class Real>
  26 void test_agreement_with_1d()
  27 {
  28     std::cout << "Testing with 1D interpolation on type "
  29               << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  30     std::mt19937 gen(4723);
  31     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  32     std::vector<Real> t(100);
  33     std::vector<Eigen::Vector2d> y(100);
  34     t[0] = dis(gen);
  35     y[0][0] = dis(gen);
  36     y[0][1] = dis(gen);
  37     for (size_t i = 1; i < t.size(); ++i)
  38     {
  39         t[i] = t[i-1] + dis(gen);
  40         y[i][0] = dis(gen);
  41         y[i][1] = dis(gen);
  42     }
  43
  44     std::vector<Eigen::Vector2d> y_copy = y;
  45     std::vector<Real> t_copy = t;
  46     std::vector<Real> t_copy0 = t;
  47     std::vector<Real> t_copy1 = t;
  48
  49     std::vector<Real> y_copy0(y.size());
  50     std::vector<Real> y_copy1(y.size());
  51     for (size_t i = 0; i < y.size(); ++i) {
  52         y_copy0[i] = y[i][0];
  53         y_copy1[i] = y[i][1];
  54     }
  55
  56     boost::random::uniform_real_distribution<Real> dis2(t[0], t[t.size()-1]);
  57     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
  58     boost::math::barycentric_rational<Real> scalar_interpolator0(std::move(t_copy0), std::move(y_copy0));
  59     boost::math::barycentric_rational<Real> scalar_interpolator1(std::move(t_copy1), std::move(y_copy1));
  60
  61
  62     Eigen::Vector2d z;
  63
  64     size_t samples = 0;
  65     while (samples++ < 1000)
  66     {
  67         Real t = dis2(gen);
  68         interpolator(z, t);
  69         BOOST_CHECK_CLOSE(z[0], scalar_interpolator0(t), 10000*numeric_limits<Real>::epsilon());
  70         BOOST_CHECK_CLOSE(z[1], scalar_interpolator1(t), 10000*numeric_limits<Real>::epsilon());
  71     }
  72 }
  73
  74
  75 template<class Real>
  76 void test_interpolation_condition_eigen()
  77 {
  78     std::cout << "Testing interpolation condition for barycentric interpolation on Eigen vectors of type "
  79               << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
  80     std::mt19937 gen(4723);
  81     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
  82     std::vector<Real> t(100);
  83     std::vector<Eigen::Vector2d> y(100);
  84     t[0] = dis(gen);
  85     y[0][0] = dis(gen);
  86     y[0][1] = dis(gen);
  87     for (size_t i = 1; i < t.size(); ++i)
  88     {
  89         t[i] = t[i-1] + dis(gen);
  90         y[i][0] = dis(gen);
  91         y[i][1] = dis(gen);
  92     }
  93
  94     std::vector<Eigen::Vector2d> y_copy = y;
  95     std::vector<Real> t_copy = t;
  96     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
  97
  98     Eigen::Vector2d z;
  99     for (size_t i = 0; i < t_copy.size(); ++i)
 100     {
 101         interpolator(z, t_copy[i]);
 102         BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
 103         BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
 104     }
 105 }
 106
 107 template<class Real>
 108 void test_interpolation_condition_std_array()
 109 {
 110     std::cout << "Testing interpolation condition for barycentric interpolation on std::array vectors of type "
 111               << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 112     std::mt19937 gen(4723);
 113     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 114     std::vector<Real> t(100);
 115     std::vector<std::array<Real, 2>> y(100);
 116     t[0] = dis(gen);
 117     y[0][0] = dis(gen);
 118     y[0][1] = dis(gen);
 119     for (size_t i = 1; i < t.size(); ++i)
 120     {
 121         t[i] = t[i-1] + dis(gen);
 122         y[i][0] = dis(gen);
 123         y[i][1] = dis(gen);
 124     }
 125
 126     std::vector<std::array<Real, 2>> y_copy = y;
 127     std::vector<Real> t_copy = t;
 128     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
 129
 130     std::array<Real, 2> z;
 131     for (size_t i = 0; i < t_copy.size(); ++i)
 132     {
 133         interpolator(z, t_copy[i]);
 134         BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
 135         BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
 136     }
 137 }
 138
 139
 140 template<class Real>
 141 void test_interpolation_condition_ublas()
 142 {
 143     std::cout << "Testing interpolation condition for barycentric interpolation ublas vectors of type "
 144               << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 145     std::mt19937 gen(4723);
 146     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 147     std::vector<Real> t(100);
 148     std::vector<boost::numeric::ublas::vector<Real>> y(100);
 149     t[0] = dis(gen);
 150     y[0].resize(2);
 151     y[0][0] = dis(gen);
 152     y[0][1] = dis(gen);
 153     for (size_t i = 1; i < t.size(); ++i)
 154     {
 155         t[i] = t[i-1] + dis(gen);
 156         y[i].resize(2);
 157         y[i][0] = dis(gen);
 158         y[i][1] = dis(gen);
 159     }
 160
 161     std::vector<Real> t_copy = t;
 162     std::vector<boost::numeric::ublas::vector<Real>> y_copy = y;
 163
 164     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
 165
 166     boost::numeric::ublas::vector<Real> z(2);
 167     for (size_t i = 0; i < t_copy.size(); ++i)
 168     {
 169         interpolator(z, t_copy[i]);
 170         BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
 171         BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
 172     }
 173 }
 174
 175 template<class Real>
 176 void test_interpolation_condition_high_order()
 177 {
 178     std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name()  << "\n";
 179     std::mt19937 gen(5);
 180     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 181     std::vector<Real> t(100);
 182     std::vector<Eigen::Vector2d> y(100);
 183     t[0] = dis(gen);
 184     y[0][0] = dis(gen);
 185     y[0][1] = dis(gen);
 186     for (size_t i = 1; i < t.size(); ++i)
 187     {
 188         t[i] = t[i-1] + dis(gen);
 189         y[i][0] = dis(gen);
 190         y[i][1] = dis(gen);
 191     }
 192
 193     std::vector<Eigen::Vector2d> y_copy = y;
 194     std::vector<Real> t_copy = t;
 195     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
 196
 197     Eigen::Vector2d z;
 198     for (size_t i = 0; i < t_copy.size(); ++i)
 199     {
 200         interpolator(z, t_copy[i]);
 201         BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon());
 202         BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon());
 203     }
 204 }
 205
 206
 207 template<class Real>
 208 void test_constant_eigen()
 209 {
 210     std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on Eigen vectors of type "
 211               << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 212
 213     std::mt19937 gen(6);
 214     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 215     std::vector<Real> t(100);
 216     std::vector<Eigen::Vector2d> y(100);
 217     t[0] = dis(gen);
 218     Real constant0 = dis(gen);
 219     Real constant1 = dis(gen);
 220     y[0][0] = constant0;
 221     y[0][1] = constant1;
 222     for (size_t i = 1; i < t.size(); ++i)
 223     {
 224         t[i] = t[i-1] + dis(gen);
 225         y[i][0] = constant0;
 226         y[i][1] = constant1;
 227     }
 228
 229     std::vector<Eigen::Vector2d> y_copy = y;
 230     std::vector<Real> t_copy = t;
 231     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
 232
 233     Eigen::Vector2d z;
 234     for (size_t i = 0; i < t_copy.size(); ++i)
 235     {
 236         // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
 237         Real t = t_copy[i] + dis(gen);
 238         z = interpolator(t);
 239         BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
 240         BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon()));
 241         Eigen::Vector2d zprime = interpolator.prime(t);
 242         Real zero_0 = zprime[0];
 243         Real zero_1 = zprime[1];
 244         BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon()));
 245         BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
 246     }
 247 }
 248
 249
 250 template<class Real>
 251 void test_constant_std_array()
 252 {
 253     std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on std::array vectors of type "
 254               << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 255
 256     std::mt19937 gen(6);
 257     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 258     std::vector<Real> t(100);
 259     std::vector<std::array<Real, 2>> y(100);
 260     t[0] = dis(gen);
 261     Real constant0 = dis(gen);
 262     Real constant1 = dis(gen);
 263     y[0][0] = constant0;
 264     y[0][1] = constant1;
 265     for (size_t i = 1; i < t.size(); ++i)
 266     {
 267         t[i] = t[i-1] + dis(gen);
 268         y[i][0] = constant0;
 269         y[i][1] = constant1;
 270     }
 271
 272     std::vector<std::array<Real,2>> y_copy = y;
 273     std::vector<Real> t_copy = t;
 274     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
 275
 276     std::array<Real, 2> z;
 277     for (size_t i = 0; i < t_copy.size(); ++i)
 278     {
 279         // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
 280         Real t = t_copy[i] + dis(gen);
 281         z = interpolator(t);
 282         BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
 283         BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon()));
 284         std::array<Real, 2> zprime = interpolator.prime(t);
 285         Real zero_0 = zprime[0];
 286         Real zero_1 = zprime[1];
 287         BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon()));
 288         BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
 289     }
 290 }
 291
 292
 293 template<class Real>
 294 void test_constant_high_order()
 295 {
 296     std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 297
 298     std::mt19937 gen(6);
 299     boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
 300     std::vector<Real> t(100);
 301     std::vector<Eigen::Vector2d> y(100);
 302     t[0] = dis(gen);
 303     Real constant0 = dis(gen);
 304     Real constant1 = dis(gen);
 305     y[0][0] = constant0;
 306     y[0][1] = constant1;
 307     for (size_t i = 1; i < t.size(); ++i)
 308     {
 309         t[i] = t[i-1] + dis(gen);
 310         y[i][0] = constant0;
 311         y[i][1] = constant1;
 312     }
 313
 314     std::vector<Eigen::Vector2d> y_copy = y;
 315     std::vector<Real> t_copy = t;
 316     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
 317
 318     Eigen::Vector2d z;
 319     for (size_t i = 0; i < t_copy.size(); ++i)
 320     {
 321         // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
 322         Real t = t_copy[i] + dis(gen);
 323         z = interpolator(t);
 324         BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon()));
 325         BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon()));
 326         Eigen::Vector2d zprime = interpolator.prime(t);
 327         Real zero_0 = zprime[0];
 328         Real zero_1 = zprime[1];
 329         BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon()));
 330         BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon()));
 331     }
 332 }
 333
 334
 335 template<class Real>
 336 void test_weights()
 337 {
 338     std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 339
 340     std::mt19937 gen(9);
 341     boost::random::uniform_real_distribution<Real> dis(0.005, 0.01);
 342     std::vector<Real> t(100);
 343     std::vector<Eigen::Vector2d> y(100);
 344     t[0] = dis(gen);
 345     y[0][0] = dis(gen);
 346     y[0][1] = dis(gen);
 347     for (size_t i = 1; i < t.size(); ++i)
 348     {
 349         t[i] = t[i-1] + dis(gen);
 350         y[i][0] = dis(gen);
 351         y[i][1] = dis(gen);
 352     }
 353
 354     std::vector<Eigen::Vector2d> y_copy = y;
 355     std::vector<Real> t_copy = t;
 356     boost::math::detail::vector_barycentric_rational_imp<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 1);
 357
 358     for (size_t i = 1; i < t_copy.size() - 1; ++i)
 359     {
 360         Real w = interpolator.weight(i);
 361         Real w_expect = 1/(t_copy[i] - t_copy[i - 1]) + 1/(t_copy[i+1] - t_copy[i]);
 362         if (i % 2 == 0)
 363         {
 364             BOOST_CHECK_CLOSE(w, -w_expect, 0.00001);
 365         }
 366         else
 367         {
 368             BOOST_CHECK_CLOSE(w, w_expect, 0.00001);
 369         }
 370     }
 371 }
 372
 373
 374 BOOST_AUTO_TEST_CASE(vector_barycentric_rational)
 375 {
 376     test_weights<double>();
 377     test_constant_eigen<double>();
 378     test_constant_std_array<double>();
 379     test_constant_high_order<double>();
 380     test_interpolation_condition_eigen<double>();
 381     test_interpolation_condition_ublas<double>();
 382     test_interpolation_condition_std_array<double>();
 383     test_interpolation_condition_high_order<double>();
 384     test_agreement_with_1d<double>();
 385 }