]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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 | #ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_DETAIL_HPP | |
8 | #define BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_DETAIL_HPP | |
9 | ||
10 | #include <limits> | |
11 | #include <cmath> | |
12 | #include <vector> | |
13 | #include <memory> | |
14 | #include <boost/math/constants/constants.hpp> | |
15 | #include <boost/math/special_functions/fpclassify.hpp> | |
20effc67 | 16 | #include <boost/math/special_functions/trunc.hpp> |
92f5a8d4 TL |
17 | |
18 | namespace boost{ namespace math{ namespace interpolators{ namespace detail{ | |
19 | ||
20 | ||
21 | template <class Real> | |
22 | class cardinal_cubic_b_spline_imp | |
23 | { | |
24 | public: | |
25 | // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them. | |
26 | // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1). | |
27 | template <class BidiIterator> | |
28 | cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size, | |
29 | Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(), | |
30 | Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()); | |
31 | ||
32 | Real operator()(Real x) const; | |
33 | ||
34 | Real prime(Real x) const; | |
35 | ||
36 | Real double_prime(Real x) const; | |
37 | ||
38 | private: | |
39 | std::vector<Real> m_beta; | |
40 | Real m_h_inv; | |
41 | Real m_a; | |
42 | Real m_avg; | |
43 | }; | |
44 | ||
45 | ||
46 | ||
47 | template <class Real> | |
48 | Real b3_spline(Real x) | |
49 | { | |
50 | using std::abs; | |
51 | Real absx = abs(x); | |
52 | if (absx < 1) | |
53 | { | |
54 | Real y = 2 - absx; | |
55 | Real z = 1 - absx; | |
56 | return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z); | |
57 | } | |
58 | if (absx < 2) | |
59 | { | |
60 | Real y = 2 - absx; | |
61 | return boost::math::constants::sixth<Real>()*y*y*y; | |
62 | } | |
1e59de90 | 63 | return static_cast<Real>(0); |
92f5a8d4 TL |
64 | } |
65 | ||
66 | template<class Real> | |
67 | Real b3_spline_prime(Real x) | |
68 | { | |
69 | if (x < 0) | |
70 | { | |
71 | return -b3_spline_prime(-x); | |
72 | } | |
73 | ||
74 | if (x < 1) | |
75 | { | |
76 | return x*(3*boost::math::constants::half<Real>()*x - 2); | |
77 | } | |
78 | if (x < 2) | |
79 | { | |
80 | return -boost::math::constants::half<Real>()*(2 - x)*(2 - x); | |
81 | } | |
1e59de90 | 82 | return static_cast<Real>(0); |
92f5a8d4 TL |
83 | } |
84 | ||
85 | template<class Real> | |
86 | Real b3_spline_double_prime(Real x) | |
87 | { | |
88 | if (x < 0) | |
89 | { | |
90 | return b3_spline_double_prime(-x); | |
91 | } | |
92 | ||
93 | if (x < 1) | |
94 | { | |
95 | return 3*x - 2; | |
96 | } | |
97 | if (x < 2) | |
98 | { | |
99 | return (2 - x); | |
100 | } | |
1e59de90 | 101 | return static_cast<Real>(0); |
92f5a8d4 TL |
102 | } |
103 | ||
104 | ||
105 | template <class Real> | |
106 | template <class BidiIterator> | |
107 | cardinal_cubic_b_spline_imp<Real>::cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size, | |
108 | Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0) | |
109 | { | |
110 | using boost::math::constants::third; | |
111 | ||
112 | std::size_t length = end_p - f; | |
113 | ||
114 | if (length < 5) | |
115 | { | |
116 | if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative)) | |
117 | { | |
118 | throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n"); | |
119 | } | |
120 | if (length < 3) | |
121 | { | |
122 | throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n"); | |
123 | } | |
124 | } | |
125 | ||
126 | if (boost::math::isnan(left_endpoint)) | |
127 | { | |
128 | throw std::logic_error("Left endpoint is NAN; this is disallowed.\n"); | |
129 | } | |
130 | if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)()) | |
131 | { | |
132 | throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n"); | |
133 | } | |
134 | if (step_size <= 0) | |
135 | { | |
136 | throw std::logic_error("The step size must be strictly > 0.\n"); | |
137 | } | |
138 | // Storing the inverse of the stepsize does provide a measurable speedup. | |
139 | // It's not huge, but nonetheless worthwhile. | |
140 | m_h_inv = 1/step_size; | |
141 | ||
142 | // Following Kress's notation, s'(a) = a1, s'(b) = b1 | |
143 | Real a1 = left_endpoint_derivative; | |
144 | // See the finite-difference table on Wikipedia for reference on how | |
145 | // to construct high-order estimates for one-sided derivatives: | |
146 | // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference | |
147 | // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method. | |
148 | if (boost::math::isnan(a1)) | |
149 | { | |
150 | // For simple functions (linear, quadratic, so on) | |
151 | // almost all the error comes from derivative estimation. | |
152 | // This does pairwise summation which gives us another digit of accuracy over naive summation. | |
153 | Real t0 = 4*(f[1] + third<Real>()*f[3]); | |
154 | Real t1 = -(25*third<Real>()*f[0] + f[4])/4 - 3*f[2]; | |
155 | a1 = m_h_inv*(t0 + t1); | |
156 | } | |
157 | ||
158 | Real b1 = right_endpoint_derivative; | |
159 | if (boost::math::isnan(b1)) | |
160 | { | |
161 | size_t n = length - 1; | |
162 | Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]); | |
163 | Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4 - 3*f[n - 2]; | |
164 | ||
165 | b1 = m_h_inv*(t0 + t1); | |
166 | } | |
167 | ||
168 | // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h ) | |
169 | // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy. | |
170 | m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN()); | |
171 | ||
172 | // Since the splines have compact support, they decay to zero very fast outside the endpoints. | |
173 | // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the | |
174 | // boundary [a,b] without massive error. | |
175 | // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average. | |
176 | // This algorithm for computing the average is recommended in | |
177 | // http://www.heikohoffmann.de/htmlthesis/node134.html | |
178 | Real t = 1; | |
179 | for (size_t i = 0; i < length; ++i) | |
180 | { | |
181 | if (boost::math::isnan(f[i])) | |
182 | { | |
183 | std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n"; | |
184 | throw std::logic_error(err); | |
185 | } | |
186 | m_avg += (f[i] - m_avg) / t; | |
187 | t += 1; | |
188 | } | |
189 | ||
190 | ||
191 | // Now we must solve an almost-tridiagonal system, which requires O(N) operations. | |
192 | // There are, in fact 5 diagonals, but they only differ from zero on the first and last row, | |
193 | // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows. | |
194 | // See Kress, equations 8.41 | |
195 | // The the "tridiagonal" matrix is: | |
196 | // 1 0 -1 | |
197 | // 1 4 1 | |
198 | // 1 4 1 | |
199 | // 1 4 1 | |
200 | // .... | |
201 | // 1 4 1 | |
202 | // 1 0 -1 | |
203 | // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good. | |
204 | std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN()); | |
205 | std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN()); | |
206 | ||
207 | rhs[0] = -2*step_size*a1; | |
208 | rhs[rhs.size() - 1] = -2*step_size*b1; | |
209 | ||
210 | super_diagonal[0] = 0; | |
211 | ||
212 | for(size_t i = 1; i < rhs.size() - 1; ++i) | |
213 | { | |
214 | rhs[i] = 6*(f[i - 1] - m_avg); | |
215 | super_diagonal[i] = 1; | |
216 | } | |
217 | ||
218 | ||
219 | // One step of row reduction on the first row to patch up the 5-diagonal problem: | |
220 | // 1 0 -1 | r0 | |
221 | // 1 4 1 | r1 | |
222 | // mapsto: | |
223 | // 1 0 -1 | r0 | |
224 | // 0 4 2 | r1 - r0 | |
225 | // mapsto | |
226 | // 1 0 -1 | r0 | |
227 | // 0 1 1/2| (r1 - r0)/4 | |
228 | super_diagonal[1] = 0.5; | |
229 | rhs[1] = (rhs[1] - rhs[0])/4; | |
230 | ||
231 | // Now do a tridiagonal row reduction the standard way, until just before the last row: | |
232 | for (size_t i = 2; i < rhs.size() - 1; ++i) | |
233 | { | |
234 | Real diagonal = 4 - super_diagonal[i - 1]; | |
235 | rhs[i] = (rhs[i] - rhs[i - 1])/diagonal; | |
236 | super_diagonal[i] /= diagonal; | |
237 | } | |
238 | ||
239 | // Now the last row, which is in the form | |
240 | // 1 sd[n-3] 0 | rhs[n-3] | |
241 | // 0 1 sd[n-2] | rhs[n-2] | |
242 | // 1 0 -1 | rhs[n-1] | |
243 | Real final_subdiag = -super_diagonal[rhs.size() - 3]; | |
244 | rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag; | |
245 | Real final_diag = -1/final_subdiag; | |
246 | // Now we're here: | |
247 | // 1 sd[n-3] 0 | rhs[n-3] | |
248 | // 0 1 sd[n-2] | rhs[n-2] | |
249 | // 0 1 final_diag | (rhs[n-1] - rhs[n-3])/diag | |
250 | ||
251 | final_diag = final_diag - super_diagonal[rhs.size() - 2]; | |
252 | rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2]; | |
253 | ||
254 | ||
255 | // Back substitutions: | |
256 | m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag; | |
257 | for(size_t i = rhs.size() - 2; i > 0; --i) | |
258 | { | |
259 | m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1]; | |
260 | } | |
261 | m_beta[0] = m_beta[2] + rhs[0]; | |
262 | } | |
263 | ||
264 | template<class Real> | |
265 | Real cardinal_cubic_b_spline_imp<Real>::operator()(Real x) const | |
266 | { | |
267 | // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms, | |
268 | // just the (at most 5) whose support overlaps the argument. | |
269 | Real z = m_avg; | |
270 | Real t = m_h_inv*(x - m_a) + 1; | |
271 | ||
272 | using std::max; | |
273 | using std::min; | |
274 | using std::ceil; | |
275 | using std::floor; | |
276 | ||
1e59de90 TL |
277 | size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)))); |
278 | size_t k_max = static_cast<size_t>((max)((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), 0l)); | |
92f5a8d4 TL |
279 | |
280 | for (size_t k = k_min; k <= k_max; ++k) | |
281 | { | |
282 | z += m_beta[k]*b3_spline(t - k); | |
283 | } | |
284 | ||
285 | return z; | |
286 | } | |
287 | ||
288 | template<class Real> | |
289 | Real cardinal_cubic_b_spline_imp<Real>::prime(Real x) const | |
290 | { | |
291 | Real z = 0; | |
292 | Real t = m_h_inv*(x - m_a) + 1; | |
293 | ||
294 | using std::max; | |
295 | using std::min; | |
296 | using std::ceil; | |
297 | using std::floor; | |
298 | ||
1e59de90 TL |
299 | size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)))); |
300 | size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2)))); | |
92f5a8d4 TL |
301 | |
302 | for (size_t k = k_min; k <= k_max; ++k) | |
303 | { | |
304 | z += m_beta[k]*b3_spline_prime(t - k); | |
305 | } | |
306 | return z*m_h_inv; | |
307 | } | |
308 | ||
309 | template<class Real> | |
310 | Real cardinal_cubic_b_spline_imp<Real>::double_prime(Real x) const | |
311 | { | |
312 | Real z = 0; | |
313 | Real t = m_h_inv*(x - m_a) + 1; | |
314 | ||
315 | using std::max; | |
316 | using std::min; | |
317 | using std::ceil; | |
318 | using std::floor; | |
319 | ||
1e59de90 TL |
320 | size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2)))); |
321 | size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2)))); | |
92f5a8d4 TL |
322 | |
323 | for (size_t k = k_min; k <= k_max; ++k) | |
324 | { | |
325 | z += m_beta[k]*b3_spline_double_prime(t - k); | |
326 | } | |
327 | return z*m_h_inv*m_h_inv; | |
328 | } | |
329 | ||
330 | }}}} | |
331 | #endif |