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