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