]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | // Copyright Nick Thompson, 2019 |
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 vector_barycentric_rational | |
8 | ||
9 | #include <cmath> | |
10 | #include <random> | |
11 | #include <array> | |
12 | #include <Eigen/Dense> | |
13 | #include <boost/numeric/ublas/vector.hpp> | |
14 | #include <boost/random/uniform_real_distribution.hpp> | |
15 | #include <boost/type_index.hpp> | |
16 | #include <boost/test/included/unit_test.hpp> | |
17 | #include <boost/test/tools/floating_point_comparison.hpp> | |
18 | #include <boost/math/interpolators/barycentric_rational.hpp> | |
19 | #include <boost/math/interpolators/vector_barycentric_rational.hpp> | |
20 | ||
21 | using std::sqrt; | |
22 | using std::abs; | |
23 | using std::numeric_limits; | |
24 | ||
25 | template<class Real> | |
26 | void test_agreement_with_1d() | |
27 | { | |
28 | std::cout << "Testing with 1D interpolation on type " | |
29 | << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
30 | std::mt19937 gen(4723); | |
31 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
32 | std::vector<Real> t(100); | |
33 | std::vector<Eigen::Vector2d> y(100); | |
34 | t[0] = dis(gen); | |
35 | y[0][0] = dis(gen); | |
36 | y[0][1] = dis(gen); | |
37 | for (size_t i = 1; i < t.size(); ++i) | |
38 | { | |
39 | t[i] = t[i-1] + dis(gen); | |
40 | y[i][0] = dis(gen); | |
41 | y[i][1] = dis(gen); | |
42 | } | |
43 | ||
44 | std::vector<Eigen::Vector2d> y_copy = y; | |
45 | std::vector<Real> t_copy = t; | |
46 | std::vector<Real> t_copy0 = t; | |
47 | std::vector<Real> t_copy1 = t; | |
48 | ||
49 | std::vector<Real> y_copy0(y.size()); | |
50 | std::vector<Real> y_copy1(y.size()); | |
51 | for (size_t i = 0; i < y.size(); ++i) { | |
52 | y_copy0[i] = y[i][0]; | |
53 | y_copy1[i] = y[i][1]; | |
54 | } | |
55 | ||
56 | boost::random::uniform_real_distribution<Real> dis2(t[0], t[t.size()-1]); | |
1e59de90 TL |
57 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y)); |
58 | boost::math::interpolators::barycentric_rational<Real> scalar_interpolator0(std::move(t_copy0), std::move(y_copy0)); | |
59 | boost::math::interpolators::barycentric_rational<Real> scalar_interpolator1(std::move(t_copy1), std::move(y_copy1)); | |
92f5a8d4 TL |
60 | |
61 | ||
62 | Eigen::Vector2d z; | |
63 | ||
64 | size_t samples = 0; | |
65 | while (samples++ < 1000) | |
66 | { | |
67 | Real t = dis2(gen); | |
68 | interpolator(z, t); | |
69 | BOOST_CHECK_CLOSE(z[0], scalar_interpolator0(t), 10000*numeric_limits<Real>::epsilon()); | |
70 | BOOST_CHECK_CLOSE(z[1], scalar_interpolator1(t), 10000*numeric_limits<Real>::epsilon()); | |
71 | } | |
72 | } | |
73 | ||
74 | ||
75 | template<class Real> | |
76 | void test_interpolation_condition_eigen() | |
77 | { | |
78 | std::cout << "Testing interpolation condition for barycentric interpolation on Eigen vectors of type " | |
79 | << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
80 | std::mt19937 gen(4723); | |
81 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
82 | std::vector<Real> t(100); | |
83 | std::vector<Eigen::Vector2d> y(100); | |
84 | t[0] = dis(gen); | |
85 | y[0][0] = dis(gen); | |
86 | y[0][1] = dis(gen); | |
87 | for (size_t i = 1; i < t.size(); ++i) | |
88 | { | |
89 | t[i] = t[i-1] + dis(gen); | |
90 | y[i][0] = dis(gen); | |
91 | y[i][1] = dis(gen); | |
92 | } | |
93 | ||
94 | std::vector<Eigen::Vector2d> y_copy = y; | |
95 | std::vector<Real> t_copy = t; | |
1e59de90 | 96 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y)); |
92f5a8d4 TL |
97 | |
98 | Eigen::Vector2d z; | |
99 | for (size_t i = 0; i < t_copy.size(); ++i) | |
100 | { | |
101 | interpolator(z, t_copy[i]); | |
102 | BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon()); | |
103 | BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon()); | |
104 | } | |
105 | } | |
106 | ||
107 | template<class Real> | |
108 | void test_interpolation_condition_std_array() | |
109 | { | |
110 | std::cout << "Testing interpolation condition for barycentric interpolation on std::array vectors of type " | |
111 | << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
112 | std::mt19937 gen(4723); | |
113 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
114 | std::vector<Real> t(100); | |
115 | std::vector<std::array<Real, 2>> y(100); | |
116 | t[0] = dis(gen); | |
117 | y[0][0] = dis(gen); | |
118 | y[0][1] = dis(gen); | |
119 | for (size_t i = 1; i < t.size(); ++i) | |
120 | { | |
121 | t[i] = t[i-1] + dis(gen); | |
122 | y[i][0] = dis(gen); | |
123 | y[i][1] = dis(gen); | |
124 | } | |
125 | ||
126 | std::vector<std::array<Real, 2>> y_copy = y; | |
127 | std::vector<Real> t_copy = t; | |
1e59de90 | 128 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y)); |
92f5a8d4 TL |
129 | |
130 | std::array<Real, 2> z; | |
131 | for (size_t i = 0; i < t_copy.size(); ++i) | |
132 | { | |
133 | interpolator(z, t_copy[i]); | |
134 | BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon()); | |
135 | BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon()); | |
136 | } | |
137 | } | |
138 | ||
139 | ||
140 | template<class Real> | |
141 | void test_interpolation_condition_ublas() | |
142 | { | |
143 | std::cout << "Testing interpolation condition for barycentric interpolation ublas vectors of type " | |
144 | << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
145 | std::mt19937 gen(4723); | |
146 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
147 | std::vector<Real> t(100); | |
148 | std::vector<boost::numeric::ublas::vector<Real>> y(100); | |
149 | t[0] = dis(gen); | |
150 | y[0].resize(2); | |
151 | y[0][0] = dis(gen); | |
152 | y[0][1] = dis(gen); | |
153 | for (size_t i = 1; i < t.size(); ++i) | |
154 | { | |
155 | t[i] = t[i-1] + dis(gen); | |
156 | y[i].resize(2); | |
157 | y[i][0] = dis(gen); | |
158 | y[i][1] = dis(gen); | |
159 | } | |
160 | ||
161 | std::vector<Real> t_copy = t; | |
162 | std::vector<boost::numeric::ublas::vector<Real>> y_copy = y; | |
163 | ||
1e59de90 | 164 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y)); |
92f5a8d4 TL |
165 | |
166 | boost::numeric::ublas::vector<Real> z(2); | |
167 | for (size_t i = 0; i < t_copy.size(); ++i) | |
168 | { | |
169 | interpolator(z, t_copy[i]); | |
170 | BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon()); | |
171 | BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon()); | |
172 | } | |
173 | } | |
174 | ||
175 | template<class Real> | |
176 | void test_interpolation_condition_high_order() | |
177 | { | |
178 | std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
179 | std::mt19937 gen(5); | |
180 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
181 | std::vector<Real> t(100); | |
182 | std::vector<Eigen::Vector2d> y(100); | |
183 | t[0] = dis(gen); | |
184 | y[0][0] = dis(gen); | |
185 | y[0][1] = dis(gen); | |
186 | for (size_t i = 1; i < t.size(); ++i) | |
187 | { | |
188 | t[i] = t[i-1] + dis(gen); | |
189 | y[i][0] = dis(gen); | |
190 | y[i][1] = dis(gen); | |
191 | } | |
192 | ||
193 | std::vector<Eigen::Vector2d> y_copy = y; | |
194 | std::vector<Real> t_copy = t; | |
1e59de90 | 195 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5); |
92f5a8d4 TL |
196 | |
197 | Eigen::Vector2d z; | |
198 | for (size_t i = 0; i < t_copy.size(); ++i) | |
199 | { | |
200 | interpolator(z, t_copy[i]); | |
201 | BOOST_CHECK_CLOSE(z[0], y_copy[i][0], 100*numeric_limits<Real>::epsilon()); | |
202 | BOOST_CHECK_CLOSE(z[1], y_copy[i][1], 100*numeric_limits<Real>::epsilon()); | |
203 | } | |
204 | } | |
205 | ||
206 | ||
207 | template<class Real> | |
208 | void test_constant_eigen() | |
209 | { | |
210 | std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on Eigen vectors of type " | |
211 | << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
212 | ||
213 | std::mt19937 gen(6); | |
214 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
215 | std::vector<Real> t(100); | |
216 | std::vector<Eigen::Vector2d> y(100); | |
217 | t[0] = dis(gen); | |
218 | Real constant0 = dis(gen); | |
219 | Real constant1 = dis(gen); | |
220 | y[0][0] = constant0; | |
221 | y[0][1] = constant1; | |
222 | for (size_t i = 1; i < t.size(); ++i) | |
223 | { | |
224 | t[i] = t[i-1] + dis(gen); | |
225 | y[i][0] = constant0; | |
226 | y[i][1] = constant1; | |
227 | } | |
228 | ||
229 | std::vector<Eigen::Vector2d> y_copy = y; | |
230 | std::vector<Real> t_copy = t; | |
1e59de90 | 231 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y)); |
92f5a8d4 TL |
232 | |
233 | Eigen::Vector2d z; | |
234 | for (size_t i = 0; i < t_copy.size(); ++i) | |
235 | { | |
236 | // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. | |
237 | Real t = t_copy[i] + dis(gen); | |
238 | z = interpolator(t); | |
239 | BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon())); | |
240 | BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon())); | |
241 | Eigen::Vector2d zprime = interpolator.prime(t); | |
242 | Real zero_0 = zprime[0]; | |
243 | Real zero_1 = zprime[1]; | |
244 | BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon())); | |
245 | BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon())); | |
246 | } | |
247 | } | |
248 | ||
249 | ||
250 | template<class Real> | |
251 | void test_constant_std_array() | |
252 | { | |
253 | std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on std::array vectors of type " | |
254 | << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
255 | ||
256 | std::mt19937 gen(6); | |
257 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
258 | std::vector<Real> t(100); | |
259 | std::vector<std::array<Real, 2>> y(100); | |
260 | t[0] = dis(gen); | |
261 | Real constant0 = dis(gen); | |
262 | Real constant1 = dis(gen); | |
263 | y[0][0] = constant0; | |
264 | y[0][1] = constant1; | |
265 | for (size_t i = 1; i < t.size(); ++i) | |
266 | { | |
267 | t[i] = t[i-1] + dis(gen); | |
268 | y[i][0] = constant0; | |
269 | y[i][1] = constant1; | |
270 | } | |
271 | ||
272 | std::vector<std::array<Real,2>> y_copy = y; | |
273 | std::vector<Real> t_copy = t; | |
1e59de90 | 274 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y)); |
92f5a8d4 TL |
275 | |
276 | std::array<Real, 2> z; | |
277 | for (size_t i = 0; i < t_copy.size(); ++i) | |
278 | { | |
279 | // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. | |
280 | Real t = t_copy[i] + dis(gen); | |
281 | z = interpolator(t); | |
282 | BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon())); | |
283 | BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon())); | |
284 | std::array<Real, 2> zprime = interpolator.prime(t); | |
285 | Real zero_0 = zprime[0]; | |
286 | Real zero_1 = zprime[1]; | |
287 | BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon())); | |
288 | BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon())); | |
289 | } | |
290 | } | |
291 | ||
292 | ||
293 | template<class Real> | |
294 | void test_constant_high_order() | |
295 | { | |
296 | std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
297 | ||
298 | std::mt19937 gen(6); | |
299 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
300 | std::vector<Real> t(100); | |
301 | std::vector<Eigen::Vector2d> y(100); | |
302 | t[0] = dis(gen); | |
303 | Real constant0 = dis(gen); | |
304 | Real constant1 = dis(gen); | |
305 | y[0][0] = constant0; | |
306 | y[0][1] = constant1; | |
307 | for (size_t i = 1; i < t.size(); ++i) | |
308 | { | |
309 | t[i] = t[i-1] + dis(gen); | |
310 | y[i][0] = constant0; | |
311 | y[i][1] = constant1; | |
312 | } | |
313 | ||
314 | std::vector<Eigen::Vector2d> y_copy = y; | |
315 | std::vector<Real> t_copy = t; | |
1e59de90 | 316 | boost::math::interpolators::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5); |
92f5a8d4 TL |
317 | |
318 | Eigen::Vector2d z; | |
319 | for (size_t i = 0; i < t_copy.size(); ++i) | |
320 | { | |
321 | // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. | |
322 | Real t = t_copy[i] + dis(gen); | |
323 | z = interpolator(t); | |
324 | BOOST_CHECK_CLOSE(z[0], constant0, 100*sqrt(numeric_limits<Real>::epsilon())); | |
325 | BOOST_CHECK_CLOSE(z[1], constant1, 100*sqrt(numeric_limits<Real>::epsilon())); | |
326 | Eigen::Vector2d zprime = interpolator.prime(t); | |
327 | Real zero_0 = zprime[0]; | |
328 | Real zero_1 = zprime[1]; | |
329 | BOOST_CHECK_SMALL(zero_0, sqrt(numeric_limits<Real>::epsilon())); | |
330 | BOOST_CHECK_SMALL(zero_1, sqrt(numeric_limits<Real>::epsilon())); | |
331 | } | |
332 | } | |
333 | ||
334 | ||
335 | template<class Real> | |
336 | void test_weights() | |
337 | { | |
338 | std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
339 | ||
340 | std::mt19937 gen(9); | |
341 | boost::random::uniform_real_distribution<Real> dis(0.005, 0.01); | |
342 | std::vector<Real> t(100); | |
343 | std::vector<Eigen::Vector2d> y(100); | |
344 | t[0] = dis(gen); | |
345 | y[0][0] = dis(gen); | |
346 | y[0][1] = dis(gen); | |
347 | for (size_t i = 1; i < t.size(); ++i) | |
348 | { | |
349 | t[i] = t[i-1] + dis(gen); | |
350 | y[i][0] = dis(gen); | |
351 | y[i][1] = dis(gen); | |
352 | } | |
353 | ||
354 | std::vector<Eigen::Vector2d> y_copy = y; | |
355 | std::vector<Real> t_copy = t; | |
1e59de90 | 356 | boost::math::interpolators::detail::vector_barycentric_rational_imp<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 1); |
92f5a8d4 TL |
357 | |
358 | for (size_t i = 1; i < t_copy.size() - 1; ++i) | |
359 | { | |
360 | Real w = interpolator.weight(i); | |
361 | Real w_expect = 1/(t_copy[i] - t_copy[i - 1]) + 1/(t_copy[i+1] - t_copy[i]); | |
362 | if (i % 2 == 0) | |
363 | { | |
364 | BOOST_CHECK_CLOSE(w, -w_expect, 0.00001); | |
365 | } | |
366 | else | |
367 | { | |
368 | BOOST_CHECK_CLOSE(w, w_expect, 0.00001); | |
369 | } | |
370 | } | |
371 | } | |
372 | ||
373 | ||
374 | BOOST_AUTO_TEST_CASE(vector_barycentric_rational) | |
375 | { | |
376 | test_weights<double>(); | |
377 | test_constant_eigen<double>(); | |
378 | test_constant_std_array<double>(); | |
379 | test_constant_high_order<double>(); | |
380 | test_interpolation_condition_eigen<double>(); | |
381 | test_interpolation_condition_ublas<double>(); | |
382 | test_interpolation_condition_std_array<double>(); | |
383 | test_interpolation_condition_high_order<double>(); | |
384 | test_agreement_with_1d<double>(); | |
385 | } |