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