1 /*
2  * Copyright Nick Thompson, 2020
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 <array>
13 #include <vector>
14 #include <boost/math/interpolators/cubic_hermite.hpp>
15 #include <boost/math/special_functions/next.hpp>
16 #include <boost/circular_buffer.hpp>
17 #ifdef BOOST_HAS_FLOAT128
18 #include <boost/multiprecision/float128.hpp>
19 using boost::multiprecision::float128;
20 #endif
21 
22 
23 using boost::math::interpolators::cubic_hermite;
24 using boost::math::interpolators::cardinal_cubic_hermite;
25 using boost::math::interpolators::cardinal_cubic_hermite_aos;
26 
27 
28 template<typename Real>
test_constant()29 void test_constant()
30 {
31     Real x0 = 0;
32     std::vector<Real> x{x0,1,2,3, 9, 22, 81};
33     std::vector<Real> y(x.size());
34     for (auto & t : y)
35     {
36         t = 7;
37     }
38 
39     std::vector<Real> dydx(x.size(), Real(0));
40     auto x_copy = x;
41     auto y_copy = y;
42     auto dydx_copy = dydx;
43     auto hermite_spline = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
44 
45     // Now check the boundaries:
46     Real tlo = x.front();
47     Real thi = x.back();
48     int samples = 5000;
49     int i = 0;
50     while (i++ < samples)
51     {
52         CHECK_ULP_CLOSE(Real(7), hermite_spline(tlo), 2);
53         CHECK_ULP_CLOSE(Real(7), hermite_spline(thi), 2);
54         CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(tlo), 2);
55         CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(thi), 2);
56         tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
57         thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
58     }
59 
60     boost::circular_buffer<Real> x_buf(x.size());
61     for (auto & t : x) {
62         x_buf.push_back(t);
63     }
64 
65     boost::circular_buffer<Real> y_buf(x.size());
66     for (auto & t : y) {
67         y_buf.push_back(t);
68     }
69 
70     boost::circular_buffer<Real> dydx_buf(x.size());
71     for (auto & t : dydx) {
72         dydx_buf.push_back(t);
73     }
74 
75     auto circular_hermite_spline = cubic_hermite(std::move(x_buf), std::move(y_buf), std::move(dydx_buf));
76 
77     for (Real t = x[0]; t <= x.back(); t += 0.25) {
78         CHECK_ULP_CLOSE(Real(7), circular_hermite_spline(t), 2);
79         CHECK_ULP_CLOSE(Real(0), circular_hermite_spline.prime(t), 2);
80     }
81 
82     circular_hermite_spline.push_back(x.back() + 1, 7, 0);
83     CHECK_ULP_CLOSE(Real(0), circular_hermite_spline.prime(x.back()+1), 2);
84 
85 }
86 
87 template<typename Real>
test_linear()88 void test_linear()
89 {
90     std::vector<Real> x{0,1,2,3};
91     std::vector<Real> y{0,1,2,3};
92     std::vector<Real> dydx{1,1,1,1};
93 
94     auto x_copy = x;
95     auto y_copy = y;
96     auto dydx_copy = dydx;
97     auto hermite_spline = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
98 
99     CHECK_ULP_CLOSE(y[0], hermite_spline(x[0]), 0);
100     CHECK_ULP_CLOSE(Real(1)/Real(2), hermite_spline(Real(1)/Real(2)), 10);
101     CHECK_ULP_CLOSE(y[1], hermite_spline(x[1]), 0);
102     CHECK_ULP_CLOSE(Real(3)/Real(2), hermite_spline(Real(3)/Real(2)), 10);
103     CHECK_ULP_CLOSE(y[2], hermite_spline(x[2]), 0);
104     CHECK_ULP_CLOSE(Real(5)/Real(2), hermite_spline(Real(5)/Real(2)), 10);
105     CHECK_ULP_CLOSE(y[3], hermite_spline(x[3]), 0);
106 
107     x.resize(45);
108     y.resize(45);
109     dydx.resize(45);
110     for (size_t i = 0; i < x.size(); ++i) {
111         x[i] = i;
112         y[i] = i;
113         dydx[i] = 1;
114     }
115 
116     x_copy = x;
117     y_copy = y;
118     dydx_copy = dydx;
119     hermite_spline = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
120     for (Real t = 0; t < x.back(); t += 0.5) {
121         CHECK_ULP_CLOSE(t, hermite_spline(t), 0);
122         CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(t), 0);
123     }
124 
125     boost::circular_buffer<Real> x_buf(x.size());
126     for (auto & t : x) {
127         x_buf.push_back(t);
128     }
129 
130     boost::circular_buffer<Real> y_buf(x.size());
131     for (auto & t : y) {
132         y_buf.push_back(t);
133     }
134 
135     boost::circular_buffer<Real> dydx_buf(x.size());
136     for (auto & t : dydx) {
137         dydx_buf.push_back(t);
138     }
139 
140     auto circular_hermite_spline = cubic_hermite(std::move(x_buf), std::move(y_buf), std::move(dydx_buf));
141 
142     for (Real t = x[0]; t <= x.back(); t += 0.25) {
143         CHECK_ULP_CLOSE(t, circular_hermite_spline(t), 2);
144         CHECK_ULP_CLOSE(Real(1), circular_hermite_spline.prime(t), 2);
145     }
146 
147     circular_hermite_spline.push_back(x.back() + 1, y.back()+1, 1);
148 
149     CHECK_ULP_CLOSE(Real(y.back() + 1), circular_hermite_spline(Real(x.back()+1)), 2);
150     CHECK_ULP_CLOSE(Real(1), circular_hermite_spline.prime(Real(x.back()+1)), 2);
151 
152 }
153 
154 template<typename Real>
test_quadratic()155 void test_quadratic()
156 {
157     std::vector<Real> x(50);
158     std::default_random_engine rd;
159     std::uniform_real_distribution<Real> dis(0.1,1);
160     Real x0 = dis(rd);
161     x[0] = x0;
162     for (size_t i = 1; i < x.size(); ++i) {
163         x[i] = x[i-1] + dis(rd);
164     }
165     Real xmax = x.back();
166 
167     std::vector<Real> y(x.size());
168     std::vector<Real> dydx(x.size());
169     for (size_t i = 0; i < x.size(); ++i) {
170         y[i] = x[i]*x[i]/2;
171         dydx[i] = x[i];
172     }
173 
174     auto s = cubic_hermite(std::move(x), std::move(y), std::move(dydx));
175     for (Real t = x0; t <= xmax; t+= 0.0125)
176     {
177         CHECK_ULP_CLOSE(t*t/2, s(t), 5);
178         CHECK_ULP_CLOSE(t, s.prime(t), 65);
179     }
180 }
181 
182 template<typename Real>
test_interpolation_condition()183 void test_interpolation_condition()
184 {
185     for (size_t n = 4; n < 50; ++n) {
186         std::vector<Real> x(n);
187         std::vector<Real> y(n);
188         std::vector<Real> dydx(n);
189         std::default_random_engine rd;
190         std::uniform_real_distribution<Real> dis(0,1);
191         Real x0 = dis(rd);
192         x[0] = x0;
193         y[0] = dis(rd);
194         for (size_t i = 1; i < n; ++i) {
195             x[i] = x[i-1] + dis(rd);
196             y[i] = dis(rd);
197             dydx[i] = dis(rd);
198         }
199 
200         auto x_copy = x;
201         auto y_copy = y;
202         auto dydx_copy = dydx;
203         auto s = cubic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy));
204         //std::cout << "s = " << s << "\n";
205         for (size_t i = 0; i < x.size(); ++i) {
206             CHECK_ULP_CLOSE(y[i], s(x[i]), 2);
207             CHECK_ULP_CLOSE(dydx[i], s.prime(x[i]), 2);
208         }
209     }
210 }
211 
212 template<typename Real>
test_cardinal_constant()213 void test_cardinal_constant()
214 {
215     Real x0 = 0;
216     Real dx = 2;
217     std::vector<Real> y(25);
218     for (auto & t : y) {
219         t = 7;
220     }
221 
222     std::vector<Real> dydx(y.size(), Real(0));
223 
224     auto hermite_spline = cardinal_cubic_hermite(std::move(y), std::move(dydx), x0, dx);
225 
226     for (Real t = x0; t <= x0 + 24*dx; t += 0.25) {
227         CHECK_ULP_CLOSE(Real(7), hermite_spline(t), 2);
228         CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(t), 2);
229     }
230 
231     // Array of structs:
232 
233     std::vector<std::array<Real, 2>> data(25);
234     for (auto & t : data) {
235         t[0] = 7;
236         t[1] = 0;
237     }
238     auto hermite_spline_aos = cardinal_cubic_hermite_aos(std::move(data), x0, dx);
239 
240     for (Real t = x0; t <= x0 + 24*dx; t += 0.25) {
241         if (!CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(t), 2)) {
242             std::cerr << "  Wrong evaluation at t = " << t << "\n";
243         }
244         if (!CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(t), 2)) {
245             std::cerr << "  Wrong evaluation at t = " << t << "\n";
246         }
247     }
248 
249     // Now check the boundaries:
250     Real tlo = x0;
251     Real thi = x0 + (25-1)*dx;
252     int samples = 5000;
253     int i = 0;
254     while (i++ < samples)
255     {
256         CHECK_ULP_CLOSE(Real(7), hermite_spline(tlo), 2);
257         CHECK_ULP_CLOSE(Real(7), hermite_spline(thi), 2);
258         CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(tlo), 2);
259         CHECK_ULP_CLOSE(Real(7), hermite_spline_aos(thi), 2);
260         CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(tlo), 2);
261         CHECK_ULP_CLOSE(Real(0), hermite_spline.prime(thi), 2);
262         CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(tlo), 2);
263         CHECK_ULP_CLOSE(Real(0), hermite_spline_aos.prime(thi), 2);
264 
265         tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
266         thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
267     }
268 
269 }
270 
271 
272 template<typename Real>
test_cardinal_linear()273 void test_cardinal_linear()
274 {
275     Real x0 = 0;
276     Real dx = 1;
277     std::vector<Real> y{0,1,2,3};
278     std::vector<Real> dydx{1,1,1,1};
279     auto y_copy = y;
280     auto dydx_copy = dydx;
281     auto hermite_spline = cardinal_cubic_hermite(std::move(y_copy), std::move(dydx_copy), x0, dx);
282 
283     CHECK_ULP_CLOSE(y[0], hermite_spline(0), 0);
284     CHECK_ULP_CLOSE(Real(1)/Real(2), hermite_spline(Real(1)/Real(2)), 10);
285     CHECK_ULP_CLOSE(y[1], hermite_spline(1), 0);
286     CHECK_ULP_CLOSE(Real(3)/Real(2), hermite_spline(Real(3)/Real(2)), 10);
287     CHECK_ULP_CLOSE(y[2], hermite_spline(2), 0);
288     CHECK_ULP_CLOSE(Real(5)/Real(2), hermite_spline(Real(5)/Real(2)), 10);
289     CHECK_ULP_CLOSE(y[3], hermite_spline(3), 0);
290 
291 
292     y.resize(45);
293     dydx.resize(45);
294     for (size_t i = 0; i < y.size(); ++i) {
295         y[i] = i;
296         dydx[i] = 1;
297     }
298 
299     hermite_spline = cardinal_cubic_hermite(std::move(y), std::move(dydx), x0, dx);
300     for (Real t = 0; t < 44; t += 0.5) {
301         CHECK_ULP_CLOSE(t, hermite_spline(t), 0);
302         CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(t), 0);
303     }
304 
305     std::vector<std::array<Real, 2>> data(45);
306     for (size_t i = 0; i < data.size(); ++i) {
307         data[i][0] = i;
308         data[i][1] = 1;
309     }
310 
311     auto hermite_spline_aos = cardinal_cubic_hermite_aos(std::move(data), x0, dx);
312     for (Real t = 0; t < 44; t += 0.5) {
313         CHECK_ULP_CLOSE(t, hermite_spline_aos(t), 0);
314         CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(t), 0);
315     }
316 
317     Real tlo = x0;
318     Real thi = x0 + (45-1)*dx;
319     int samples = 5000;
320     int i = 0;
321     while (i++ < samples)
322     {
323         CHECK_ULP_CLOSE(Real(tlo), hermite_spline(tlo), 2);
324         CHECK_ULP_CLOSE(Real(thi), hermite_spline(thi), 2);
325         CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(tlo), 2);
326         CHECK_ULP_CLOSE(Real(1), hermite_spline.prime(thi), 2);
327         CHECK_ULP_CLOSE(Real(tlo), hermite_spline_aos(tlo), 2);
328         CHECK_ULP_CLOSE(Real(thi), hermite_spline_aos(thi), 2);
329         CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(tlo), 2);
330         CHECK_ULP_CLOSE(Real(1), hermite_spline_aos.prime(thi), 2);
331 
332         tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
333         thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
334     }
335 
336 
337 }
338 
339 
340 template<typename Real>
test_cardinal_quadratic()341 void test_cardinal_quadratic()
342 {
343     Real x0 = -1;
344     Real dx = Real(1)/Real(256);
345 
346     std::vector<Real> y(50);
347     std::vector<Real> dydx(y.size());
348     for (size_t i = 0; i < y.size(); ++i) {
349         Real x = x0 + i*dx;
350         y[i] = x*x/2;
351         dydx[i] = x;
352     }
353 
354     auto s = cardinal_cubic_hermite(std::move(y), std::move(dydx), x0, dx);
355     for (Real t = x0; t <= x0 + 49*dx; t+= 0.0125)
356     {
357         CHECK_ULP_CLOSE(t*t/2, s(t), 12);
358         CHECK_ULP_CLOSE(t, s.prime(t), 70);
359     }
360 
361     std::vector<std::array<Real, 2>> data(50);
362     for (size_t i = 0; i < data.size(); ++i) {
363         Real x = x0 + i*dx;
364         data[i][0] = x*x/2;
365         data[i][1] = x;
366     }
367 
368 
369     auto saos = cardinal_cubic_hermite_aos(std::move(data), x0, dx);
370     for (Real t = x0; t <= x0 + 49*dx; t+= 0.0125)
371     {
372         CHECK_ULP_CLOSE(t*t/2, saos(t), 12);
373         CHECK_ULP_CLOSE(t, saos.prime(t), 70);
374     }
375 
376     auto [tlo, thi] = s.domain();
377     int samples = 5000;
378     int i = 0;
379     while (i++ < samples)
380     {
381         CHECK_ULP_CLOSE(Real(tlo*tlo/2), s(tlo), 3);
382         CHECK_ULP_CLOSE(Real(thi*thi/2), s(thi), 3);
383         CHECK_ULP_CLOSE(Real(tlo), s.prime(tlo), 3);
384         CHECK_ULP_CLOSE(Real(thi), s.prime(thi), 3);
385         CHECK_ULP_CLOSE(Real(tlo*tlo/2), saos(tlo), 3);
386         CHECK_ULP_CLOSE(Real(thi*thi/2), saos(thi), 3);
387         CHECK_ULP_CLOSE(Real(tlo), saos.prime(tlo), 3);
388         CHECK_ULP_CLOSE(Real(thi), saos.prime(thi), 3);
389 
390         tlo = boost::math::nextafter(tlo, std::numeric_limits<Real>::max());
391         thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest());
392     }
393 }
394 
395 
396 template<typename Real>
test_cardinal_interpolation_condition()397 void test_cardinal_interpolation_condition()
398 {
399     for (size_t n = 4; n < 50; ++n) {
400         std::vector<Real> y(n);
401         std::vector<Real> dydx(n);
402         std::default_random_engine rd;
403         std::uniform_real_distribution<Real> dis(0.1,1);
404         Real x0 = Real(2);
405         Real dx = Real(1)/Real(128);
406         for (size_t i = 0; i < n; ++i) {
407             y[i] = dis(rd);
408             dydx[i] = dis(rd);
409         }
410 
411         auto y_copy = y;
412         auto dydx_copy = dydx;
413         auto s = cardinal_cubic_hermite(std::move(y_copy), std::move(dydx_copy), x0, dx);
414         for (size_t i = 0; i < y.size(); ++i) {
415             CHECK_ULP_CLOSE(y[i], s(x0 + i*dx), 2);
416             CHECK_ULP_CLOSE(dydx[i], s.prime(x0 + i*dx), 2);
417         }
418     }
419 }
420 
421 
422 
main()423 int main()
424 {
425     test_constant<float>();
426     test_linear<float>();
427     test_quadratic<float>();
428     test_interpolation_condition<float>();
429     test_cardinal_constant<float>();
430     test_cardinal_linear<float>();
431     test_cardinal_quadratic<float>();
432     test_cardinal_interpolation_condition<float>();
433 
434     test_constant<double>();
435     test_linear<double>();
436     test_quadratic<double>();
437     test_interpolation_condition<double>();
438     test_cardinal_constant<double>();
439     test_cardinal_linear<double>();
440     test_cardinal_quadratic<double>();
441     test_cardinal_interpolation_condition<double>();
442 
443     test_constant<long double>();
444     test_linear<long double>();
445     test_quadratic<long double>();
446     test_interpolation_condition<long double>();
447     test_cardinal_constant<long double>();
448     test_cardinal_linear<long double>();
449     test_cardinal_quadratic<long double>();
450     test_cardinal_interpolation_condition<long double>();
451 
452 
453 #ifdef BOOST_HAS_FLOAT128
454     test_constant<float128>();
455     test_linear<float128>();
456     test_cardinal_constant<float128>();
457     test_cardinal_linear<float128>();
458 #endif
459 
460     return boost::math::test::report_errors();
461 }
462