]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
1 | /* |
2 | * Copyright Nick Thompson, 2020 | |
3 | * Use, modification and distribution are subject to the | |
4 | * Boost Software License, Version 1.0. (See accompanying file | |
5 | * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
6 | */ | |
7 | ||
8 | #include "math_unit_test.hpp" | |
9 | #include <numeric> | |
10 | #include <utility> | |
11 | #include <random> | |
12 | #include <boost/math/interpolators/pchip.hpp> | |
13 | #include <boost/circular_buffer.hpp> | |
1e59de90 | 14 | #include <boost/assert.hpp> |
f67539c2 TL |
15 | #ifdef BOOST_HAS_FLOAT128 |
16 | #include <boost/multiprecision/float128.hpp> | |
17 | using boost::multiprecision::float128; | |
18 | #endif | |
19 | ||
20 | ||
21 | using boost::math::interpolators::pchip; | |
22 | ||
23 | template<typename Real> | |
24 | void test_constant() | |
25 | { | |
26 | ||
27 | std::vector<Real> x{0,1,2,3, 9, 22, 81}; | |
28 | std::vector<Real> y(x.size()); | |
29 | for (auto & t : y) { | |
30 | t = 7; | |
31 | } | |
32 | ||
33 | auto x_copy = x; | |
34 | auto y_copy = y; | |
35 | auto pchip_spline = pchip(std::move(x_copy), std::move(y_copy)); | |
36 | //std::cout << "Constant value pchip spline = " << pchip_spline << "\n"; | |
37 | ||
38 | for (Real t = x[0]; t <= x.back(); t += 0.25) { | |
39 | CHECK_ULP_CLOSE(Real(7), pchip_spline(t), 2); | |
40 | CHECK_ULP_CLOSE(Real(0), pchip_spline.prime(t), 2); | |
41 | } | |
42 | ||
43 | boost::circular_buffer<Real> x_buf(x.size()); | |
44 | for (auto & t : x) { | |
45 | x_buf.push_back(t); | |
46 | } | |
47 | ||
48 | boost::circular_buffer<Real> y_buf(x.size()); | |
49 | for (auto & t : y) { | |
50 | y_buf.push_back(t); | |
51 | } | |
52 | ||
53 | auto circular_pchip_spline = pchip(std::move(x_buf), std::move(y_buf)); | |
54 | ||
55 | for (Real t = x[0]; t <= x.back(); t += 0.25) { | |
56 | CHECK_ULP_CLOSE(Real(7), circular_pchip_spline(t), 2); | |
57 | CHECK_ULP_CLOSE(Real(0), pchip_spline.prime(t), 2); | |
58 | } | |
59 | ||
60 | circular_pchip_spline.push_back(x.back() + 1, 7); | |
61 | CHECK_ULP_CLOSE(Real(0), circular_pchip_spline.prime(x.back()+1), 2); | |
62 | ||
63 | } | |
64 | ||
65 | template<typename Real> | |
66 | void test_linear() | |
67 | { | |
68 | std::vector<Real> x{0,1,2,3}; | |
69 | std::vector<Real> y{0,1,2,3}; | |
70 | ||
71 | auto x_copy = x; | |
72 | auto y_copy = y; | |
73 | auto pchip_spline = pchip(std::move(x_copy), std::move(y_copy)); | |
74 | ||
75 | CHECK_ULP_CLOSE(y[0], pchip_spline(x[0]), 0); | |
76 | CHECK_ULP_CLOSE(Real(1)/Real(2), pchip_spline(Real(1)/Real(2)), 10); | |
77 | CHECK_ULP_CLOSE(y[1], pchip_spline(x[1]), 0); | |
78 | CHECK_ULP_CLOSE(Real(3)/Real(2), pchip_spline(Real(3)/Real(2)), 10); | |
79 | CHECK_ULP_CLOSE(y[2], pchip_spline(x[2]), 0); | |
80 | CHECK_ULP_CLOSE(Real(5)/Real(2), pchip_spline(Real(5)/Real(2)), 10); | |
81 | CHECK_ULP_CLOSE(y[3], pchip_spline(x[3]), 0); | |
82 | ||
83 | x.resize(45); | |
84 | y.resize(45); | |
85 | for (size_t i = 0; i < x.size(); ++i) { | |
86 | x[i] = i; | |
87 | y[i] = i; | |
88 | } | |
89 | ||
90 | x_copy = x; | |
91 | y_copy = y; | |
92 | pchip_spline = pchip(std::move(x_copy), std::move(y_copy)); | |
93 | for (Real t = 0; t < x.back(); t += 0.5) { | |
94 | CHECK_ULP_CLOSE(t, pchip_spline(t), 0); | |
95 | CHECK_ULP_CLOSE(Real(1), pchip_spline.prime(t), 0); | |
96 | } | |
97 | ||
98 | x_copy = x; | |
99 | y_copy = y; | |
100 | // Test endpoint derivatives: | |
101 | pchip_spline = pchip(std::move(x_copy), std::move(y_copy), Real(1), Real(1)); | |
102 | for (Real t = 0; t < x.back(); t += 0.5) { | |
103 | CHECK_ULP_CLOSE(t, pchip_spline(t), 0); | |
104 | CHECK_ULP_CLOSE(Real(1), pchip_spline.prime(t), 0); | |
105 | } | |
106 | ||
107 | ||
108 | boost::circular_buffer<Real> x_buf(x.size()); | |
109 | for (auto & t : x) { | |
110 | x_buf.push_back(t); | |
111 | } | |
112 | ||
113 | boost::circular_buffer<Real> y_buf(x.size()); | |
114 | for (auto & t : y) { | |
115 | y_buf.push_back(t); | |
116 | } | |
117 | ||
118 | auto circular_pchip_spline = pchip(std::move(x_buf), std::move(y_buf)); | |
119 | ||
120 | for (Real t = x[0]; t <= x.back(); t += 0.25) { | |
121 | CHECK_ULP_CLOSE(t, circular_pchip_spline(t), 2); | |
122 | CHECK_ULP_CLOSE(Real(1), circular_pchip_spline.prime(t), 2); | |
123 | } | |
124 | ||
125 | circular_pchip_spline.push_back(x.back() + 1, y.back()+1); | |
126 | ||
127 | CHECK_ULP_CLOSE(Real(y.back() + 1), circular_pchip_spline(Real(x.back()+1)), 2); | |
128 | CHECK_ULP_CLOSE(Real(1), circular_pchip_spline.prime(Real(x.back()+1)), 2); | |
129 | ||
130 | } | |
131 | ||
132 | template<typename Real> | |
133 | void test_interpolation_condition() | |
134 | { | |
135 | for (size_t n = 4; n < 50; ++n) { | |
136 | std::vector<Real> x(n); | |
137 | std::vector<Real> y(n); | |
138 | std::default_random_engine rd; | |
139 | std::uniform_real_distribution<Real> dis(0,1); | |
140 | Real x0 = dis(rd); | |
141 | x[0] = x0; | |
142 | y[0] = dis(rd); | |
143 | for (size_t i = 1; i < n; ++i) { | |
144 | x[i] = x[i-1] + dis(rd); | |
145 | y[i] = dis(rd); | |
146 | } | |
147 | ||
148 | auto x_copy = x; | |
149 | auto y_copy = y; | |
150 | auto s = pchip(std::move(x_copy), std::move(y_copy)); | |
151 | //std::cout << "s = " << s << "\n"; | |
152 | for (size_t i = 0; i < x.size(); ++i) { | |
153 | CHECK_ULP_CLOSE(y[i], s(x[i]), 2); | |
154 | } | |
155 | ||
156 | x_copy = x; | |
157 | y_copy = y; | |
158 | // The interpolation condition is not affected by the endpoint derivatives, even though these derivatives might be super weird: | |
159 | s = pchip(std::move(x_copy), std::move(y_copy), Real(0), Real(0)); | |
160 | //std::cout << "s = " << s << "\n"; | |
161 | for (size_t i = 0; i < x.size(); ++i) { | |
162 | CHECK_ULP_CLOSE(y[i], s(x[i]), 2); | |
163 | } | |
164 | ||
165 | } | |
166 | } | |
167 | ||
168 | template<typename Real> | |
169 | void test_monotonicity() | |
170 | { | |
171 | for (size_t n = 4; n < 50; ++n) { | |
172 | std::vector<Real> x(n); | |
173 | std::vector<Real> y(n); | |
174 | std::default_random_engine rd; | |
175 | std::uniform_real_distribution<Real> dis(0,1); | |
176 | Real x0 = dis(rd); | |
177 | x[0] = x0; | |
178 | y[0] = dis(rd); | |
179 | // Monotone increasing: | |
180 | for (size_t i = 1; i < n; ++i) { | |
181 | x[i] = x[i-1] + dis(rd); | |
182 | y[i] = y[i-1] + dis(rd); | |
183 | } | |
184 | ||
185 | auto x_copy = x; | |
186 | auto y_copy = y; | |
187 | auto s = pchip(std::move(x_copy), std::move(y_copy)); | |
188 | //std::cout << "s = " << s << "\n"; | |
189 | for (size_t i = 0; i < x.size() - 1; ++i) { | |
190 | Real tmin = x[i]; | |
191 | Real tmax = x[i+1]; | |
192 | Real val = y[i]; | |
193 | CHECK_ULP_CLOSE(y[i], s(x[i]), 2); | |
194 | for (Real t = tmin; t < tmax; t += (tmax-tmin)/16) { | |
195 | Real greater_val = s(t); | |
1e59de90 | 196 | BOOST_ASSERT(val <= greater_val); |
f67539c2 TL |
197 | val = greater_val; |
198 | } | |
199 | } | |
200 | ||
201 | ||
202 | x[0] = dis(rd); | |
203 | y[0] = dis(rd); | |
204 | // Monotone decreasing: | |
205 | for (size_t i = 1; i < n; ++i) { | |
206 | x[i] = x[i-1] + dis(rd); | |
207 | y[i] = y[i-1] - dis(rd); | |
208 | } | |
209 | ||
210 | x_copy = x; | |
211 | y_copy = y; | |
212 | s = pchip(std::move(x_copy), std::move(y_copy)); | |
213 | //std::cout << "s = " << s << "\n"; | |
214 | for (size_t i = 0; i < x.size() - 1; ++i) { | |
215 | Real tmin = x[i]; | |
216 | Real tmax = x[i+1]; | |
217 | Real val = y[i]; | |
218 | CHECK_ULP_CLOSE(y[i], s(x[i]), 2); | |
219 | for (Real t = tmin; t < tmax; t += (tmax-tmin)/16) { | |
220 | Real lesser_val = s(t); | |
1e59de90 | 221 | BOOST_ASSERT(val >= lesser_val); |
f67539c2 TL |
222 | val = lesser_val; |
223 | } | |
224 | } | |
225 | ||
226 | } | |
227 | } | |
228 | ||
229 | ||
230 | int main() | |
231 | { | |
1e59de90 | 232 | #if (__GNUC__ > 7) || defined(_MSC_VER) || defined(__clang__) |
f67539c2 TL |
233 | test_constant<float>(); |
234 | test_linear<float>(); | |
235 | test_interpolation_condition<float>(); | |
236 | test_monotonicity<float>(); | |
237 | ||
238 | test_constant<double>(); | |
239 | test_linear<double>(); | |
240 | test_interpolation_condition<double>(); | |
241 | test_monotonicity<double>(); | |
242 | ||
243 | test_constant<long double>(); | |
244 | test_linear<long double>(); | |
245 | test_interpolation_condition<long double>(); | |
246 | test_monotonicity<long double>(); | |
247 | ||
248 | #ifdef BOOST_HAS_FLOAT128 | |
249 | test_constant<float128>(); | |
250 | test_linear<float128>(); | |
251 | #endif | |
1e59de90 | 252 | #endif |
f67539c2 TL |
253 | return boost::math::test::report_errors(); |
254 | } |