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