]>
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 barycentric_rational | |
8 | ||
9 | #include <random> | |
10 | #include <boost/random/uniform_real_distribution.hpp> | |
11 | #include <boost/type_index.hpp> | |
12 | #include <boost/test/included/unit_test.hpp> | |
13 | #include <boost/test/floating_point_comparison.hpp> | |
14 | #include <boost/math/interpolators/barycentric_rational.hpp> | |
15 | #include <boost/multiprecision/cpp_bin_float.hpp> | |
16 | ||
17 | #ifdef BOOST_HAS_FLOAT128 | |
18 | #include <boost/multiprecision/float128.hpp> | |
19 | #endif | |
20 | ||
21 | using boost::multiprecision::cpp_bin_float_50; | |
22 | ||
23 | template<class Real> | |
24 | void test_interpolation_condition() | |
25 | { | |
26 | std::cout << "Testing interpolation condition for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
27 | std::mt19937 gen(4); | |
28 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
29 | std::vector<Real> x(500); | |
30 | std::vector<Real> y(500); | |
31 | x[0] = dis(gen); | |
32 | y[0] = dis(gen); | |
33 | for (size_t i = 1; i < x.size(); ++i) | |
34 | { | |
35 | x[i] = x[i-1] + dis(gen); | |
36 | y[i] = dis(gen); | |
37 | } | |
38 | ||
39 | boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size()); | |
40 | ||
41 | for (size_t i = 0; i < x.size(); ++i) | |
42 | { | |
43 | auto z = interpolator(x[i]); | |
44 | BOOST_CHECK_CLOSE(z, y[i], 100*std::numeric_limits<Real>::epsilon()); | |
45 | } | |
46 | } | |
47 | ||
48 | template<class Real> | |
49 | void test_interpolation_condition_high_order() | |
50 | { | |
51 | std::cout << "Testing interpolation condition in high order for barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
52 | std::mt19937 gen(5); | |
53 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
54 | std::vector<Real> x(500); | |
55 | std::vector<Real> y(500); | |
56 | x[0] = dis(gen); | |
57 | y[0] = dis(gen); | |
58 | for (size_t i = 1; i < x.size(); ++i) | |
59 | { | |
60 | x[i] = x[i-1] + dis(gen); | |
61 | y[i] = dis(gen); | |
62 | } | |
63 | ||
64 | // Order 5 approximation: | |
65 | boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5); | |
66 | ||
67 | for (size_t i = 0; i < x.size(); ++i) | |
68 | { | |
69 | auto z = interpolator(x[i]); | |
70 | BOOST_CHECK_CLOSE(z, y[i], 100*std::numeric_limits<Real>::epsilon()); | |
71 | } | |
72 | } | |
73 | ||
74 | ||
75 | template<class Real> | |
76 | void test_constant() | |
77 | { | |
78 | std::cout << "Testing that constants are interpolated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
79 | ||
80 | std::mt19937 gen(6); | |
81 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
82 | std::vector<Real> x(500); | |
83 | std::vector<Real> y(500); | |
84 | Real constant = -8; | |
85 | x[0] = dis(gen); | |
86 | y[0] = constant; | |
87 | for (size_t i = 1; i < x.size(); ++i) | |
88 | { | |
89 | x[i] = x[i-1] + dis(gen); | |
90 | y[i] = y[0]; | |
91 | } | |
92 | ||
93 | boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size()); | |
94 | ||
95 | for (size_t i = 0; i < x.size(); ++i) | |
96 | { | |
97 | // Don't evaluate the constant at x[i]; that's already tested in the interpolation condition test. | |
98 | auto z = interpolator(x[i] + dis(gen)); | |
99 | BOOST_CHECK_CLOSE(z, constant, 100*sqrt(std::numeric_limits<Real>::epsilon())); | |
100 | } | |
101 | } | |
102 | ||
103 | template<class Real> | |
104 | void test_constant_high_order() | |
105 | { | |
106 | std::cout << "Testing that constants are interpolated correctly in high order using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
107 | ||
108 | std::mt19937 gen(7); | |
109 | boost::random::uniform_real_distribution<Real> dis(0.1f, 1); | |
110 | std::vector<Real> x(500); | |
111 | std::vector<Real> y(500); | |
112 | Real constant = 5; | |
113 | x[0] = dis(gen); | |
114 | y[0] = constant; | |
115 | for (size_t i = 1; i < x.size(); ++i) | |
116 | { | |
117 | x[i] = x[i-1] + dis(gen); | |
118 | y[i] = y[0]; | |
119 | } | |
120 | ||
121 | // Set interpolation order to 7: | |
122 | boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 7); | |
123 | ||
124 | for (size_t i = 0; i < x.size(); ++i) | |
125 | { | |
126 | auto z = interpolator(x[i] + dis(gen)); | |
127 | BOOST_CHECK_CLOSE(z, constant, 1000*sqrt(std::numeric_limits<Real>::epsilon())); | |
128 | } | |
129 | } | |
130 | ||
131 | ||
132 | template<class Real> | |
133 | void test_runge() | |
134 | { | |
135 | 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"; | |
136 | ||
137 | std::mt19937 gen(8); | |
138 | boost::random::uniform_real_distribution<Real> dis(0.005f, 0.01f); | |
139 | std::vector<Real> x(500); | |
140 | std::vector<Real> y(500); | |
141 | x[0] = -2; | |
142 | y[0] = 1/(1+25*x[0]*x[0]); | |
143 | for (size_t i = 1; i < x.size(); ++i) | |
144 | { | |
145 | x[i] = x[i-1] + dis(gen); | |
146 | y[i] = 1/(1+25*x[i]*x[i]); | |
147 | } | |
148 | ||
149 | boost::math::barycentric_rational<Real> interpolator(x.data(), y.data(), y.size(), 5); | |
150 | ||
151 | for (size_t i = 0; i < x.size(); ++i) | |
152 | { | |
153 | auto t = x[i] + dis(gen); | |
154 | auto z = interpolator(t); | |
155 | BOOST_CHECK_CLOSE(z, 1/(1+25*t*t), 0.03); | |
156 | } | |
157 | } | |
158 | ||
159 | template<class Real> | |
160 | void test_weights() | |
161 | { | |
162 | std::cout << "Testing weights are calculated correctly using barycentric interpolation on type " << boost::typeindex::type_id<Real>().pretty_name() << "\n"; | |
163 | ||
164 | std::mt19937 gen(9); | |
165 | boost::random::uniform_real_distribution<Real> dis(0.005, 0.01); | |
166 | std::vector<Real> x(500); | |
167 | std::vector<Real> y(500); | |
168 | x[0] = -2; | |
169 | y[0] = 1/(1+25*x[0]*x[0]); | |
170 | for (size_t i = 1; i < x.size(); ++i) | |
171 | { | |
172 | x[i] = x[i-1] + dis(gen); | |
173 | y[i] = 1/(1+25*x[i]*x[i]); | |
174 | } | |
175 | ||
176 | boost::math::detail::barycentric_rational_imp<Real> interpolator(x.data(), x.data() + x.size(), y.data(), 0); | |
177 | ||
178 | for (size_t i = 0; i < x.size(); ++i) | |
179 | { | |
180 | auto w = interpolator.weight(i); | |
181 | if (i % 2 == 0) | |
182 | { | |
183 | BOOST_CHECK_CLOSE(w, 1, 0.00001); | |
184 | } | |
185 | else | |
186 | { | |
187 | BOOST_CHECK_CLOSE(w, -1, 0.00001); | |
188 | } | |
189 | } | |
190 | ||
191 | // d = 1: | |
192 | interpolator = boost::math::detail::barycentric_rational_imp<Real>(x.data(), x.data() + x.size(), y.data(), 1); | |
193 | ||
194 | for (size_t i = 1; i < x.size() -1; ++i) | |
195 | { | |
196 | auto w = interpolator.weight(i); | |
197 | auto w_expect = 1/(x[i] - x[i - 1]) + 1/(x[i+1] - x[i]); | |
198 | if (i % 2 == 0) | |
199 | { | |
200 | BOOST_CHECK_CLOSE(w, -w_expect, 0.00001); | |
201 | } | |
202 | else | |
203 | { | |
204 | BOOST_CHECK_CLOSE(w, w_expect, 0.00001); | |
205 | } | |
206 | } | |
207 | ||
208 | } | |
209 | ||
210 | ||
211 | BOOST_AUTO_TEST_CASE(barycentric_rational) | |
212 | { | |
213 | test_weights<double>(); | |
214 | test_constant<float>(); | |
215 | test_constant<double>(); | |
216 | test_constant<long double>(); | |
217 | test_constant<cpp_bin_float_50>(); | |
218 | ||
219 | test_constant_high_order<float>(); | |
220 | test_constant_high_order<double>(); | |
221 | test_constant_high_order<long double>(); | |
222 | test_constant_high_order<cpp_bin_float_50>(); | |
223 | ||
224 | test_interpolation_condition<float>(); | |
225 | test_interpolation_condition<double>(); | |
226 | test_interpolation_condition<long double>(); | |
227 | test_interpolation_condition<cpp_bin_float_50>(); | |
228 | ||
229 | test_interpolation_condition_high_order<float>(); | |
230 | test_interpolation_condition_high_order<double>(); | |
231 | test_interpolation_condition_high_order<long double>(); | |
232 | test_interpolation_condition_high_order<cpp_bin_float_50>(); | |
233 | ||
234 | test_runge<float>(); | |
235 | test_runge<double>(); | |
236 | test_runge<long double>(); | |
237 | test_runge<cpp_bin_float_50>(); | |
238 | ||
239 | #ifdef BOOST_HAS_FLOAT128 | |
240 | test_interpolation_condition<boost::multiprecision::float128>(); | |
241 | test_constant<boost::multiprecision::float128>(); | |
242 | test_constant_high_order<boost::multiprecision::float128>(); | |
243 | test_interpolation_condition_high_order<boost::multiprecision::float128>(); | |
244 | test_runge<boost::multiprecision::float128>(); | |
245 | #endif | |
246 | ||
247 | } |