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/makima.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::makima;
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 akima = makima(std::move(x_copy), std::move(y_copy));
35
36 for (Real t = x[0]; t <= x.back(); t += 0.25) {
37 CHECK_ULP_CLOSE(Real(7), akima(t), 2);
38 CHECK_ULP_CLOSE(Real(0), akima.prime(t), 2);
39 }
40
41 boost::circular_buffer<Real> x_buf(x.size());
42 for (auto & t : x) {
43 x_buf.push_back(t);
44 }
45
46 boost::circular_buffer<Real> y_buf(x.size());
47 for (auto & t : y) {
48 y_buf.push_back(t);
49 }
50
51 auto circular_akima = makima(std::move(x_buf), std::move(y_buf));
52
53 for (Real t = x[0]; t <= x.back(); t += 0.25) {
54 CHECK_ULP_CLOSE(Real(7), circular_akima(t), 2);
55 CHECK_ULP_CLOSE(Real(0), akima.prime(t), 2);
56 }
57
58 circular_akima.push_back(x.back() + 1, 7);
59 CHECK_ULP_CLOSE(Real(0), circular_akima.prime(x.back()+1), 2);
60
61 }
62
63 template<typename Real>
test_linear()64 void test_linear()
65 {
66 std::vector<Real> x{0,1,2,3};
67 std::vector<Real> y{0,1,2,3};
68
69 auto x_copy = x;
70 auto y_copy = y;
71 auto akima = makima(std::move(x_copy), std::move(y_copy));
72
73 CHECK_ULP_CLOSE(y[0], akima(x[0]), 0);
74 CHECK_ULP_CLOSE(Real(1)/Real(2), akima(Real(1)/Real(2)), 10);
75 CHECK_ULP_CLOSE(y[1], akima(x[1]), 0);
76 CHECK_ULP_CLOSE(Real(3)/Real(2), akima(Real(3)/Real(2)), 10);
77 CHECK_ULP_CLOSE(y[2], akima(x[2]), 0);
78 CHECK_ULP_CLOSE(Real(5)/Real(2), akima(Real(5)/Real(2)), 10);
79 CHECK_ULP_CLOSE(y[3], akima(x[3]), 0);
80
81 x.resize(45);
82 y.resize(45);
83 for (size_t i = 0; i < x.size(); ++i) {
84 x[i] = i;
85 y[i] = i;
86 }
87
88 x_copy = x;
89 y_copy = y;
90 akima = makima(std::move(x_copy), std::move(y_copy));
91 for (Real t = 0; t < x.back(); t += 0.5) {
92 CHECK_ULP_CLOSE(t, akima(t), 0);
93 CHECK_ULP_CLOSE(Real(1), akima.prime(t), 0);
94 }
95
96 x_copy = x;
97 y_copy = y;
98 // Test endpoint derivatives:
99 akima = makima(std::move(x_copy), std::move(y_copy), Real(1), Real(1));
100 for (Real t = 0; t < x.back(); t += 0.5) {
101 CHECK_ULP_CLOSE(t, akima(t), 0);
102 CHECK_ULP_CLOSE(Real(1), akima.prime(t), 0);
103 }
104
105
106 boost::circular_buffer<Real> x_buf(x.size());
107 for (auto & t : x) {
108 x_buf.push_back(t);
109 }
110
111 boost::circular_buffer<Real> y_buf(x.size());
112 for (auto & t : y) {
113 y_buf.push_back(t);
114 }
115
116 auto circular_akima = makima(std::move(x_buf), std::move(y_buf));
117
118 for (Real t = x[0]; t <= x.back(); t += 0.25) {
119 CHECK_ULP_CLOSE(t, circular_akima(t), 2);
120 CHECK_ULP_CLOSE(Real(1), circular_akima.prime(t), 2);
121 }
122
123 circular_akima.push_back(x.back() + 1, y.back()+1);
124
125 CHECK_ULP_CLOSE(Real(y.back() + 1), circular_akima(Real(x.back()+1)), 2);
126 CHECK_ULP_CLOSE(Real(1), circular_akima.prime(Real(x.back()+1)), 2);
127
128 }
129
130 template<typename Real>
test_interpolation_condition()131 void test_interpolation_condition()
132 {
133 for (size_t n = 4; n < 50; ++n) {
134 std::vector<Real> x(n);
135 std::vector<Real> y(n);
136 std::default_random_engine rd;
137 std::uniform_real_distribution<Real> dis(0,1);
138 Real x0 = dis(rd);
139 x[0] = x0;
140 y[0] = dis(rd);
141 for (size_t i = 1; i < n; ++i) {
142 x[i] = x[i-1] + dis(rd);
143 y[i] = dis(rd);
144 }
145
146 auto x_copy = x;
147 auto y_copy = y;
148 auto s = makima(std::move(x_copy), std::move(y_copy));
149 //std::cout << "s = " << s << "\n";
150 for (size_t i = 0; i < x.size(); ++i) {
151 CHECK_ULP_CLOSE(y[i], s(x[i]), 2);
152 }
153
154 x_copy = x;
155 y_copy = y;
156 // The interpolation condition is not affected by the endpoint derivatives, even though these derivatives might be super weird:
157 s = makima(std::move(x_copy), std::move(y_copy), Real(0), Real(0));
158 //std::cout << "s = " << s << "\n";
159 for (size_t i = 0; i < x.size(); ++i) {
160 CHECK_ULP_CLOSE(y[i], s(x[i]), 2);
161 }
162
163 }
164 }
165
main()166 int main()
167 {
168 test_constant<float>();
169 test_linear<float>();
170 test_interpolation_condition<float>();
171
172 test_constant<double>();
173 test_linear<double>();
174 test_interpolation_condition<double>();
175
176 test_constant<long double>();
177 test_linear<long double>();
178 test_interpolation_condition<long double>();
179
180 #ifdef BOOST_HAS_FLOAT128
181 test_constant<float128>();
182 test_linear<float128>();
183 #endif
184
185 return boost::math::test::report_errors();
186 }
187