]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
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> | |
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); | |
1e59de90 | 56 | tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
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> | |
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> | |
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); | |
1e59de90 | 178 | CHECK_ULP_CLOSE(t, s.prime(t), 138); |
f67539c2 TL |
179 | } |
180 | } | |
181 | ||
182 | template<typename Real> | |
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> | |
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 | ||
1e59de90 | 265 | tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
266 | thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest()); |
267 | } | |
268 | ||
269 | } | |
270 | ||
271 | ||
272 | template<typename Real> | |
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 | ||
1e59de90 | 332 | tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
333 | thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest()); |
334 | } | |
335 | ||
336 | ||
337 | } | |
338 | ||
339 | ||
340 | template<typename Real> | |
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 | ||
1e59de90 | 390 | tlo = boost::math::nextafter(tlo, (std::numeric_limits<Real>::max)()); |
f67539c2 TL |
391 | thi = boost::math::nextafter(thi, std::numeric_limits<Real>::lowest()); |
392 | } | |
393 | } | |
394 | ||
395 | ||
396 | template<typename Real> | |
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 | ||
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 | } |