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