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>
test_agreement_with_1d()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]);
57     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
58     boost::math::barycentric_rational<Real> scalar_interpolator0(std::move(t_copy0), std::move(y_copy0));
59     boost::math::barycentric_rational<Real> scalar_interpolator1(std::move(t_copy1), std::move(y_copy1));
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>
test_interpolation_condition_eigen()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;
96     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
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>
test_interpolation_condition_std_array()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;
128     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
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>
test_interpolation_condition_ublas()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 
164     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
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>
test_interpolation_condition_high_order()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;
195     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
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>
test_constant_eigen()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;
231     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
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>
test_constant_std_array()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;
274     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y));
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>
test_constant_high_order()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;
316     boost::math::vector_barycentric_rational<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 5);
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>
test_weights()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;
356     boost::math::detail::vector_barycentric_rational_imp<decltype(t), decltype(y)> interpolator(std::move(t), std::move(y), 1);
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 
BOOST_AUTO_TEST_CASE(vector_barycentric_rational)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 }
386