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 barycentric_rational
8
9 #include <cmath>
10 #include <random>
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/tools/floating_point_comparison.hpp>
15 #include <boost/math/interpolators/barycentric_rational.hpp>
16 #include <boost/multiprecision/cpp_bin_float.hpp>
17
18 #ifdef BOOST_HAS_FLOAT128
19 #include <boost/multiprecision/float128.hpp>
20 #endif
21
22 using std::sqrt;
23 using std::abs;
24 using std::numeric_limits;
25 using boost::multiprecision::cpp_bin_float_50;
26
27 template<class Real>
test_interpolation_condition()28 void test_interpolation_condition()
29 {
30 std::cout << "Testing interpolation condition for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
31 std::mt19937 gen(4);
32 boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
33 std::vector<Real> x(100);
34 std::vector<Real> y(100);
35 x[0] = dis(gen);
36 y[0] = dis(gen);
37 for (size_t i = 1; i < x.size(); ++i)
38 {
39 x[i] = x[i-1] + dis(gen);
40 y[i] = dis(gen);
41 }
42
43 boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size());
44
45 for (size_t i = 0; i < x.size(); ++i)
46 {
47 Real z = interpolator(x[i]);
48 BOOST_CHECK_CLOSE(z, y[i], 100*numeric_limits<Real>::epsilon());
49 }
50
51 // Make sure that the move constructor does the same thing:
52 std::vector<Real> x_copy = x;
53 std::vector<Real> y_copy = y;
54 boost::math::barycentric_rational<Real> move_interpolator(std::move(x), std::move(y));
55
56 for (size_t i = 0; i < x_copy.size(); ++i)
57 {
58 Real z = move_interpolator(x_copy[i]);
59 BOOST_CHECK_CLOSE(z, y_copy[i], 100*numeric_limits<Real>::epsilon());
60 }
61 }
62
63 template<class Real>
test_interpolation_condition_high_order()64 void test_interpolation_condition_high_order()
65 {
66 std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
67 std::mt19937 gen(5);
68 boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
69 std::vector<Real> x(100);
70 std::vector<Real> y(100);
71 x[0] = dis(gen);
72 y[0] = dis(gen);
73 for (size_t i = 1; i < x.size(); ++i)
74 {
75 x[i] = x[i-1] + dis(gen);
76 y[i] = dis(gen);
77 }
78
79 // Order 5 approximation:
80 boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5);
81
82 for (size_t i = 0; i < x.size(); ++i)
83 {
84 Real z = interpolator(x[i]);
85 BOOST_CHECK_CLOSE(z, y[i], 100*numeric_limits<Real>::epsilon());
86 }
87 }
88
89
90 template<class Real>
test_constant()91 void test_constant()
92 {
93 std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
94
95 std::mt19937 gen(6);
96 boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
97 std::vector<Real> x(100);
98 std::vector<Real> y(100);
99 Real constant = -8;
100 x[0] = dis(gen);
101 y[0] = constant;
102 for (size_t i = 1; i < x.size(); ++i)
103 {
104 x[i] = x[i-1] + dis(gen);
105 y[i] = y[0];
106 }
107
108 boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size());
109
110 for (size_t i = 0; i < x.size(); ++i)
111 {
112 // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test.
113 Real t = x[i] + dis(gen);
114 Real z = interpolator(t);
115 BOOST_CHECK_CLOSE(z, constant, 100*sqrt(numeric_limits<Real>::epsilon()));
116 BOOST_CHECK_SMALL(interpolator.prime(t), sqrt(numeric_limits<Real>::epsilon()));
117 }
118 }
119
120 template<class Real>
test_constant_high_order()121 void test_constant_high_order()
122 {
123 std::cout << "Testing that constants are interpolated correctly in high order using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
124
125 std::mt19937 gen(7);
126 boost::random::uniform_real_distribution<Real> dis(0.1f, 1);
127 std::vector<Real> x(100);
128 std::vector<Real> y(100);
129 Real constant = 5;
130 x[0] = dis(gen);
131 y[0] = constant;
132 for (size_t i = 1; i < x.size(); ++i)
133 {
134 x[i] = x[i-1] + dis(gen);
135 y[i] = y[0];
136 }
137
138 // Set interpolation order to 7:
139 boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 7);
140
141 for (size_t i = 0; i < x.size(); ++i)
142 {
143 Real t = x[i] + dis(gen);
144 Real z = interpolator(t);
145 BOOST_CHECK_CLOSE(z, constant, 1000*sqrt(numeric_limits<Real>::epsilon()));
146 BOOST_CHECK_SMALL(interpolator.prime(t), 100*sqrt(numeric_limits<Real>::epsilon()));
147 }
148 }
149
150
151 template<class Real>
test_runge()152 void test_runge()
153 {
154 std::cout << "Testing interpolation of Runge's 1/(1+25x^2) function using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
155
156 std::mt19937 gen(8);
157 boost::random::uniform_real_distribution<Real> dis(0.005f, 0.01f);
158 std::vector<Real> x(100);
159 std::vector<Real> y(100);
160 x[0] = -2;
161 y[0] = 1/(1+25*x[0]*x[0]);
162 for (size_t i = 1; i < x.size(); ++i)
163 {
164 x[i] = x[i-1] + dis(gen);
165 y[i] = 1/(1+25*x[i]*x[i]);
166 }
167
168 boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5);
169
170 for (size_t i = 0; i < x.size(); ++i)
171 {
172 Real t = x[i];
173 Real z = interpolator(t);
174 BOOST_CHECK_CLOSE(z, y[i], 0.03);
175 Real z_prime = interpolator.prime(t);
176 Real num = -50*t;
177 Real denom = (1+25*t*t)*(1+25*t*t);
178 if (abs(num/denom) > 0.00001)
179 {
180 BOOST_CHECK_CLOSE_FRACTION(z_prime, num/denom, 0.03);
181 }
182 }
183
184
185 Real tol = 0.0001;
186 for (size_t i = 0; i < x.size(); ++i)
187 {
188 Real t = x[i] + dis(gen);
189 Real z = interpolator(t);
190 BOOST_CHECK_CLOSE(z, 1/(1+25*t*t), tol);
191 Real z_prime = interpolator.prime(t);
192 Real num = -50*t;
193 Real denom = (1+25*t*t)*(1+25*t*t);
194 Real runge_prime = num/denom;
195
196 if (abs(runge_prime) > 0 && abs(z_prime - runge_prime)/abs(runge_prime) > tol)
197 {
198 std::cout << "Error too high for t = " << t << " which is a distance " << t - x[i] << " from node " << i << "/" << x.size() << " associated with data (" << x[i] << ", " << y[i] << ")\n";
199 BOOST_CHECK_CLOSE_FRACTION(z_prime, runge_prime, tol);
200 }
201 }
202 }
203
204 template<class Real>
test_weights()205 void test_weights()
206 {
207 std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n";
208
209 std::mt19937 gen(9);
210 boost::random::uniform_real_distribution<Real> dis(0.005, 0.01);
211 std::vector<Real> x(100);
212 std::vector<Real> y(100);
213 x[0] = -2;
214 y[0] = 1/(1+25*x[0]*x[0]);
215 for (size_t i = 1; i < x.size(); ++i)
216 {
217 x[i] = x[i-1] + dis(gen);
218 y[i] = 1/(1+25*x[i]*x[i]);
219 }
220
221 boost::math::detail::barycentric_rational_imp<Real> interpolator(x.data(), x.data() + x.size(), y.data(), 0);
222
223 for (size_t i = 0; i < x.size(); ++i)
224 {
225 Real w = interpolator.weight(i);
226 if (i % 2 == 0)
227 {
228 BOOST_CHECK_CLOSE(w, 1, 0.00001);
229 }
230 else
231 {
232 BOOST_CHECK_CLOSE(w, -1, 0.00001);
233 }
234 }
235
236 // d = 1:
237 interpolator = boost::math::detail::barycentric_rational_imp<Real>(x.data(), x.data() + x.size(), y.data(), 1);
238
239 for (size_t i = 1; i < x.size() -1; ++i)
240 {
241 Real w = interpolator.weight(i);
242 Real w_expect = 1/(x[i] - x[i - 1]) + 1/(x[i+1] - x[i]);
243 if (i % 2 == 0)
244 {
245 BOOST_CHECK_CLOSE(w, -w_expect, 0.00001);
246 }
247 else
248 {
249 BOOST_CHECK_CLOSE(w, w_expect, 0.00001);
250 }
251 }
252
253 }
254
255
BOOST_AUTO_TEST_CASE(barycentric_rational)256 BOOST_AUTO_TEST_CASE(barycentric_rational)
257 {
258 // The tests took too long at the higher precisions.
259 // They still pass, but the CI system is starting to time out,
260 // so I figured it'd be polite to comment out the most expensive tests.
261 test_weights<double>();
262
263 test_constant<float>();
264 //test_constant<double>();
265 test_constant<long double>();
266 //test_constant<cpp_bin_float_50>();
267
268 //test_constant_high_order<float>();
269 test_constant_high_order<double>();
270 //test_constant_high_order<long double>();
271 //test_constant_high_order<cpp_bin_float_50>();
272
273 test_interpolation_condition<float>();
274 test_interpolation_condition<double>();
275 //test_interpolation_condition<long double>();
276 //test_interpolation_condition<cpp_bin_float_50>();
277
278 //test_interpolation_condition_high_order<float>();
279 test_interpolation_condition_high_order<double>();
280 //test_interpolation_condition_high_order<long double>();
281 //test_interpolation_condition_high_order<cpp_bin_float_50>();
282
283 test_runge<double>();
284 //test_runge<long double>();
285 //test_runge<cpp_bin_float_50>();
286
287 #ifdef BOOST_HAS_FLOAT128
288 //test_interpolation_condition<boost::multiprecision::float128>();
289 //test_constant<boost::multiprecision::float128>();
290 //test_constant_high_order<boost::multiprecision::float128>();
291 //test_interpolation_condition_high_order<boost::multiprecision::float128>();
292 //test_runge<boost::multiprecision::float128>();
293 #endif
294
295 }
296