]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
1 | // Copyright Nick Thompson, 2019 |
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_QUADRATIC_B_SPLINE_DETAIL_HPP | |
8 | #define BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP | |
9 | #include <vector> | |
f67539c2 TL |
10 | #include <cmath> |
11 | #include <stdexcept> | |
92f5a8d4 TL |
12 | |
13 | namespace boost{ namespace math{ namespace interpolators{ namespace detail{ | |
14 | ||
15 | template <class Real> | |
16 | Real b2_spline(Real x) { | |
17 | using std::abs; | |
18 | Real absx = abs(x); | |
19 | if (absx < 1/Real(2)) | |
20 | { | |
21 | Real y = absx - 1/Real(2); | |
22 | Real z = absx + 1/Real(2); | |
23 | return (2-y*y-z*z)/2; | |
24 | } | |
25 | if (absx < Real(3)/Real(2)) | |
26 | { | |
27 | Real y = absx - Real(3)/Real(2); | |
28 | return y*y/2; | |
29 | } | |
1e59de90 | 30 | return static_cast<Real>(0); |
92f5a8d4 TL |
31 | } |
32 | ||
33 | template <class Real> | |
34 | Real b2_spline_prime(Real x) { | |
35 | if (x < 0) { | |
36 | return -b2_spline_prime(-x); | |
37 | } | |
38 | ||
39 | if (x < 1/Real(2)) | |
40 | { | |
41 | return -2*x; | |
42 | } | |
43 | if (x < Real(3)/Real(2)) | |
44 | { | |
45 | return x - Real(3)/Real(2); | |
46 | } | |
1e59de90 | 47 | return static_cast<Real>(0); |
92f5a8d4 TL |
48 | } |
49 | ||
50 | ||
51 | template <class Real> | |
52 | class cardinal_quadratic_b_spline_detail | |
53 | { | |
54 | public: | |
55 | // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them. | |
56 | // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1). | |
57 | ||
58 | cardinal_quadratic_b_spline_detail(const Real* const y, | |
59 | size_t n, | |
60 | Real t0 /* initial time, left endpoint */, | |
61 | Real h /*spacing, stepsize*/, | |
62 | Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(), | |
63 | Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()) | |
64 | { | |
65 | if (h <= 0) { | |
66 | throw std::logic_error("Spacing must be > 0."); | |
67 | } | |
68 | m_inv_h = 1/h; | |
69 | m_t0 = t0; | |
70 | ||
71 | if (n < 3) { | |
72 | throw std::logic_error("The interpolator requires at least 3 points."); | |
73 | } | |
74 | ||
75 | using std::isnan; | |
76 | Real a; | |
77 | if (isnan(left_endpoint_derivative)) { | |
78 | // http://web.media.mit.edu/~crtaylor/calculator.html | |
79 | a = -3*y[0] + 4*y[1] - y[2]; | |
80 | } | |
81 | else { | |
82 | a = 2*h*left_endpoint_derivative; | |
83 | } | |
84 | ||
85 | Real b; | |
86 | if (isnan(right_endpoint_derivative)) { | |
87 | b = 3*y[n-1] - 4*y[n-2] + y[n-3]; | |
88 | } | |
89 | else { | |
90 | b = 2*h*right_endpoint_derivative; | |
91 | } | |
92 | ||
93 | m_alpha.resize(n + 2); | |
94 | ||
95 | // Begin row reduction: | |
96 | std::vector<Real> rhs(n + 2, std::numeric_limits<Real>::quiet_NaN()); | |
97 | std::vector<Real> super_diagonal(n + 2, std::numeric_limits<Real>::quiet_NaN()); | |
98 | ||
99 | rhs[0] = -a; | |
100 | rhs[rhs.size() - 1] = b; | |
101 | ||
102 | super_diagonal[0] = 0; | |
103 | ||
104 | for(size_t i = 1; i < rhs.size() - 1; ++i) { | |
105 | rhs[i] = 8*y[i - 1]; | |
106 | super_diagonal[i] = 1; | |
107 | } | |
108 | ||
109 | // Patch up 5-diagonal problem: | |
110 | rhs[1] = (rhs[1] - rhs[0])/6; | |
111 | super_diagonal[1] = Real(1)/Real(3); | |
112 | // First two rows are now: | |
113 | // 1 0 -1 | -2hy0' | |
114 | // 0 1 1/3| (8y0+2hy0')/6 | |
115 | ||
116 | ||
117 | // Start traditional tridiagonal row reduction: | |
118 | for (size_t i = 2; i < rhs.size() - 1; ++i) { | |
119 | Real diagonal = 6 - super_diagonal[i - 1]; | |
120 | rhs[i] = (rhs[i] - rhs[i - 1])/diagonal; | |
121 | super_diagonal[i] /= diagonal; | |
122 | } | |
123 | ||
124 | // 1 sd[n-1] 0 | rhs[n-1] | |
125 | // 0 1 sd[n] | rhs[n] | |
126 | // -1 0 1 | rhs[n+1] | |
127 | ||
128 | rhs[n+1] = rhs[n+1] + rhs[n-1]; | |
129 | Real bottom_subdiagonal = super_diagonal[n-1]; | |
130 | ||
131 | // We're here: | |
132 | // 1 sd[n-1] 0 | rhs[n-1] | |
133 | // 0 1 sd[n] | rhs[n] | |
134 | // 0 bs 1 | rhs[n+1] | |
135 | ||
136 | rhs[n+1] = (rhs[n+1]-bottom_subdiagonal*rhs[n])/(1-bottom_subdiagonal*super_diagonal[n]); | |
137 | ||
138 | m_alpha[n+1] = rhs[n+1]; | |
139 | for (size_t i = n; i > 0; --i) { | |
140 | m_alpha[i] = rhs[i] - m_alpha[i+1]*super_diagonal[i]; | |
141 | } | |
142 | m_alpha[0] = m_alpha[2] + rhs[0]; | |
143 | } | |
144 | ||
145 | Real operator()(Real t) const { | |
146 | if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { | |
147 | const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; | |
148 | throw std::domain_error(err_msg); | |
149 | } | |
150 | // Let k, gamma be defined via t = t0 + kh + gamma * h. | |
151 | // Now find all j: |k-j+1+gamma|< 3/2, or, in other words | |
152 | // j_min = ceil((t-t0)/h - 1/2) | |
153 | // j_max = floor(t-t0)/h + 5/2) | |
154 | using std::floor; | |
155 | using std::ceil; | |
156 | Real x = (t-m_t0)*m_inv_h; | |
157 | size_t j_min = ceil(x - Real(1)/Real(2)); | |
158 | size_t j_max = ceil(x + Real(5)/Real(2)); | |
159 | if (j_max >= m_alpha.size()) { | |
160 | j_max = m_alpha.size() - 1; | |
161 | } | |
162 | ||
163 | Real y = 0; | |
164 | x += 1; | |
165 | for (size_t j = j_min; j <= j_max; ++j) { | |
166 | y += m_alpha[j]*detail::b2_spline(x - j); | |
167 | } | |
168 | return y; | |
169 | } | |
170 | ||
171 | Real prime(Real t) const { | |
172 | if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { | |
173 | const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; | |
174 | throw std::domain_error(err_msg); | |
175 | } | |
176 | // Let k, gamma be defined via t = t0 + kh + gamma * h. | |
177 | // Now find all j: |k-j+1+gamma|< 3/2, or, in other words | |
178 | // j_min = ceil((t-t0)/h - 1/2) | |
179 | // j_max = floor(t-t0)/h + 5/2) | |
180 | using std::floor; | |
181 | using std::ceil; | |
182 | Real x = (t-m_t0)*m_inv_h; | |
183 | size_t j_min = ceil(x - Real(1)/Real(2)); | |
184 | size_t j_max = ceil(x + Real(5)/Real(2)); | |
185 | if (j_max >= m_alpha.size()) { | |
186 | j_max = m_alpha.size() - 1; | |
187 | } | |
188 | ||
189 | Real y = 0; | |
190 | x += 1; | |
191 | for (size_t j = j_min; j <= j_max; ++j) { | |
192 | y += m_alpha[j]*detail::b2_spline_prime(x - j); | |
193 | } | |
194 | return y*m_inv_h; | |
195 | } | |
196 | ||
197 | Real t_max() const { | |
198 | return m_t0 + (m_alpha.size()-3)/m_inv_h; | |
199 | } | |
200 | ||
201 | private: | |
202 | std::vector<Real> m_alpha; | |
203 | Real m_inv_h; | |
204 | Real m_t0; | |
205 | }; | |
206 | ||
207 | }}}} | |
208 | #endif |