libs/math/test/cardinal_b_spline_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 <cmath>
  13 #include <boost/math/special_functions/cardinal_b_spline.hpp>
  14 #include <boost/math/interpolators/detail/cubic_b_spline_detail.hpp>
  15 #ifdef BOOST_HAS_FLOAT128
  16 #include <boost/multiprecision/float128.hpp>
  17 using boost::multiprecision::float128;
  18 #endif
  19
  20 using std::abs;
  21 using boost::math::cardinal_b_spline;
  22 using boost::math::cardinal_b_spline_prime;
  23 using boost::math::forward_cardinal_b_spline;
  24 using boost::math::cardinal_b_spline_double_prime;
  25
  26
  27 template<class Real>
  28 void test_box()
  29 {
  30     Real t = cardinal_b_spline<0>(Real(1.1));
  31     Real expected = 0;
  32     CHECK_ULP_CLOSE(expected, t, 0);
  33     CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0);
  34
  35     t = cardinal_b_spline<0>(Real(-1.1));
  36     expected = 0;
  37     CHECK_ULP_CLOSE(expected, t, 0);
  38
  39     Real h = Real(1)/Real(256);
  40     for (t = -Real(1)/Real(2)+h; t < Real(1)/Real(2); t += h)
  41     {
  42         expected = 1;
  43         CHECK_ULP_CLOSE(expected, cardinal_b_spline<0>(t), 0);
  44         expected = 0;
  45         CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<0>(Real(1.1)), 0);
  46     }
  47
  48     for (t = h; t < 1; t += h)
  49     {
  50         expected = 1;
  51         CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<0>(t), 0);
  52     }
  53 }
  54
  55 template<class Real>
  56 void test_hat()
  57 {
  58     Real t = cardinal_b_spline<1>(Real(2.1));
  59     Real expected = 0;
  60     CHECK_ULP_CLOSE(expected, t, 0);
  61
  62     t = cardinal_b_spline<1>(Real(-2.1));
  63     expected = 0;
  64     CHECK_ULP_CLOSE(expected, t, 0);
  65
  66     Real h = Real(1)/Real(256);
  67     for (t = -1; t <= 1; t += h)
  68     {
  69         expected = 1-abs(t);
  70         if(!CHECK_ULP_CLOSE(expected, cardinal_b_spline<1>(t), 0) )
  71         {
  72             std::cerr << "  Problem at t = " << t << "\n";
  73         }
  74         if (t == -Real(1)) {
  75             if (!CHECK_ULP_CLOSE(Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0)) {
  76                 std::cout << "  Problem at t = " << t << "\n";
  77             }
  78         }
  79         else if (t == Real(1)) {
  80             CHECK_ULP_CLOSE(-Real(1)/Real(2), cardinal_b_spline_prime<1>(t), 0);
  81         }
  82         else if (t < 0) {
  83             CHECK_ULP_CLOSE(Real(1), cardinal_b_spline_prime<1>(t), 0);
  84         }
  85         else if (t == 0) {
  86             CHECK_ULP_CLOSE(Real(0), cardinal_b_spline_prime<1>(t), 0);
  87         }
  88         else if (t > 0) {
  89             CHECK_ULP_CLOSE(Real(-1), cardinal_b_spline_prime<1>(t), 0);
  90         }
  91     }
  92
  93     for (t = 0; t < 2; t += h)
  94     {
  95         expected = 1 - abs(t-1);
  96         CHECK_ULP_CLOSE(expected, forward_cardinal_b_spline<1>(t), 0);
  97     }
  98 }
  99
 100 template<class Real>
 101 void test_quadratic()
 102 {
 103     using std::abs;
 104     auto b2 = [](Real x) {
 105         Real absx = abs(x);
 106         if (absx >= 3/Real(2)) {
 107             return Real(0);
 108         }
 109         if (absx >= 1/Real(2)) {
 110             Real t = absx - 3/Real(2);
 111             return t*t/2;
 112         }
 113         Real t1 = absx - 1/Real(2);
 114         Real t2 = absx + 1/Real(2);
 115         return (2-t1*t1 -t2*t2)/2;
 116     };
 117
 118     auto b2_prime = [&](Real x)->Real {
 119         Real absx = abs(x);
 120         Real signx  = 1;
 121         if (x < 0) {
 122             signx = -1;
 123         }
 124         if (absx >= 3/Real(2)) {
 125             return Real(0);
 126         }
 127         if (absx >= 1/Real(2)) {
 128             return (absx - 3/Real(2))*signx;
 129         }
 130         return -2*absx*signx;
 131     };
 132
 133
 134     Real h = 1/Real(256);
 135     for (Real t = -5; t <= 5; t += h) {
 136         Real expected = b2(t);
 137         CHECK_ULP_CLOSE(expected, cardinal_b_spline<2>(t), 0);
 138         expected = b2_prime(t);
 139
 140         if (!CHECK_ULP_CLOSE(expected, cardinal_b_spline_prime<2>(t), 0))
 141         {
 142             std::cerr << "  Problem at t = " << t << "\n";
 143         }
 144
 145     }
 146 }
 147
 148 template<class Real>
 149 void test_cubic()
 150 {
 151     Real expected = Real(2)/Real(3);
 152     Real computed = cardinal_b_spline<3, Real>(0);
 153     CHECK_ULP_CLOSE(expected, computed, 0);
 154
 155     expected = Real(1)/Real(6);
 156     computed = cardinal_b_spline<3, Real>(1);
 157     CHECK_ULP_CLOSE(expected, computed, 0);
 158
 159     expected = Real(0);
 160     computed = cardinal_b_spline<3, Real>(2);
 161     CHECK_ULP_CLOSE(expected, computed, 0);
 162
 163     Real h = 1/Real(256);
 164     for (Real t = -4; t <= 4; t += h) {
 165         expected = boost::math::detail::b3_spline_prime<Real>(t);
 166         computed = cardinal_b_spline_prime<3>(t);
 167         CHECK_ULP_CLOSE(expected, computed, 0);
 168         expected = boost::math::detail::b3_spline_double_prime<Real>(t);
 169         computed = cardinal_b_spline_double_prime<3>(t);
 170         if (!CHECK_ULP_CLOSE(expected, computed, 0)) {
 171             std::cerr << "  Problem at t = " << t << "\n";
 172         }
 173     }
 174 }
 175
 176 template<class Real>
 177 void test_quintic()
 178 {
 179   Real expected = Real(11)/Real(20);
 180   Real computed = cardinal_b_spline<5, Real>(0);
 181   CHECK_ULP_CLOSE(expected, computed, 0);
 182
 183   expected = Real(13)/Real(60);
 184   computed = cardinal_b_spline<5, Real>(1);
 185   CHECK_ULP_CLOSE(expected, computed, 1);
 186
 187   expected = Real(1)/Real(120);
 188   computed = cardinal_b_spline<5, Real>(2);
 189   CHECK_ULP_CLOSE(expected, computed, 0);
 190
 191   expected = Real(0);
 192   computed = cardinal_b_spline<5, Real>(3);
 193   CHECK_ULP_CLOSE(expected, computed, 0);
 194
 195 }
 196
 197 template<unsigned n, typename Real>
 198 void test_b_spline_derivatives()
 199 {
 200     Real h = 1/Real(256);
 201     Real supp = (n+Real(1))/Real(2);
 202     for (Real t = -supp - 1; t <= supp+1; t+= h)
 203     {
 204         Real expected = cardinal_b_spline<n-1>(t+Real(1)/Real(2)) - cardinal_b_spline<n-1>(t - Real(1)/Real(2));
 205         Real computed = cardinal_b_spline_prime<n>(t);
 206         CHECK_MOLLIFIED_CLOSE(expected, computed, std::numeric_limits<Real>::epsilon());
 207
 208         expected = cardinal_b_spline<n-2>(t+1) - 2*cardinal_b_spline<n-2>(t) + cardinal_b_spline<n-2>(t-1);
 209         computed = cardinal_b_spline_double_prime<n>(t);
 210         CHECK_MOLLIFIED_CLOSE(expected, computed, 2*std::numeric_limits<Real>::epsilon());
 211     }
 212 }
 213
 214 template<unsigned n, typename Real>
 215 void test_partition_of_unity()
 216 {
 217   std::mt19937 gen(323723);
 218   Real supp = (n+1.0)/2.0;
 219   std::uniform_real_distribution<Real> dis(-supp, -supp+1);
 220
 221   for(size_t i = 0; i < 500; ++i) {
 222     Real x = dis(gen);
 223     Real one = 0;
 224     while (x < supp) {
 225         one += cardinal_b_spline<n>(x);
 226         x += 1;
 227     }
 228     if(!CHECK_ULP_CLOSE(Real(1), one, n)) {
 229       std::cerr << "  Partition of unity failure at n = " << n << "\n";
 230     }
 231   }
 232 }
 233
 234
 235 int main()
 236 {
 237     test_box<float>();
 238     test_box<double>();
 239     test_box<long double>();
 240
 241     test_hat<float>();
 242     test_hat<double>();
 243     test_hat<long double>();
 244
 245     test_quadratic<float>();
 246     test_quadratic<double>();
 247     test_quadratic<long double>();
 248
 249     test_cubic<float>();
 250     test_cubic<double>();
 251     test_cubic<long double>();
 252
 253     test_quintic<float>();
 254     test_quintic<double>();
 255     test_quintic<long double>();
 256
 257     test_partition_of_unity<1, double>();
 258     test_partition_of_unity<2, double>();
 259     test_partition_of_unity<3, double>();
 260     test_partition_of_unity<4, double>();
 261     test_partition_of_unity<5, double>();
 262     test_partition_of_unity<6, double>();
 263
 264     test_b_spline_derivatives<3, double>();
 265     test_b_spline_derivatives<4, double>();
 266     test_b_spline_derivatives<5, double>();
 267     test_b_spline_derivatives<6, double>();
 268     test_b_spline_derivatives<7, double>();
 269     test_b_spline_derivatives<8, double>();
 270     test_b_spline_derivatives<9, double>();
 271
 272     test_b_spline_derivatives<3, long double>();
 273     test_b_spline_derivatives<4, long double>();
 274     test_b_spline_derivatives<5, long double>();
 275     test_b_spline_derivatives<6, long double>();
 276     test_b_spline_derivatives<7, long double>();
 277     test_b_spline_derivatives<8, long double>();
 278     test_b_spline_derivatives<9, long double>();
 279
 280
 281 #ifdef BOOST_HAS_FLOAT128
 282     test_box<float128>();
 283     test_hat<float128>();
 284     test_quadratic<float128>();
 285     test_cubic<float128>();
 286     test_quintic<float128>();
 287 #endif
 288
 289     return boost::math::test::report_errors();
 290 }