interpolators/detail/cardinal_quadratic_b_spline_detail.hpp

0001 // Copyright Nick Thompson, 2019
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP
0008 #define BOOST_MATH_INTERPOLATORS_CARDINAL_QUADRATIC_B_SPLINE_DETAIL_HPP
0009 #include <vector>
0010 #include <cmath>
0011 #include <stdexcept>
0012
0013 namespace boost{ namespace math{ namespace interpolators{ namespace detail{
0014
0015 template <class Real>
0016 Real b2_spline(Real x) {
0017     using std::abs;
0018     Real absx = abs(x);
0019     if (absx < 1/Real(2))
0020     {
0021         Real y = absx - 1/Real(2);
0022         Real z = absx + 1/Real(2);
0023         return (2-y*y-z*z)/2;
0024     }
0025     if (absx < Real(3)/Real(2))
0026     {
0027         Real y = absx - Real(3)/Real(2);
0028         return y*y/2;
0029     }
0030     return static_cast<Real>(0);
0031 }
0032
0033 template <class Real>
0034 Real b2_spline_prime(Real x) {
0035     if (x < 0) {
0036         return -b2_spline_prime(-x);
0037     }
0038
0039     if (x < 1/Real(2))
0040     {
0041         return -2*x;
0042     }
0043     if (x < Real(3)/Real(2))
0044     {
0045         return x - Real(3)/Real(2);
0046     }
0047     return static_cast<Real>(0);
0048 }
0049
0050
0051 template <class Real>
0052 class cardinal_quadratic_b_spline_detail
0053 {
0054 public:
0055     // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
0056     // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).
0057
0058     cardinal_quadratic_b_spline_detail(const Real* const y,
0059                                 size_t n,
0060                                 Real t0 /* initial time, left endpoint */,
0061                                 Real h  /*spacing, stepsize*/,
0062                                 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
0063                                 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
0064     {
0065         if (h <= 0) {
0066             throw std::logic_error("Spacing must be > 0.");
0067         }
0068         m_inv_h = 1/h;
0069         m_t0 = t0;
0070
0071         if (n < 3) {
0072             throw std::logic_error("The interpolator requires at least 3 points.");
0073         }
0074
0075         using std::isnan;
0076         Real a;
0077         if (isnan(left_endpoint_derivative)) {
0078             // http://web.media.mit.edu/~crtaylor/calculator.html
0079             a = -3*y[0] + 4*y[1] - y[2];
0080         }
0081         else {
0082             a = 2*h*left_endpoint_derivative;
0083         }
0084
0085         Real b;
0086         if (isnan(right_endpoint_derivative)) {
0087             b = 3*y[n-1] - 4*y[n-2] + y[n-3];
0088         }
0089         else {
0090             b = 2*h*right_endpoint_derivative;
0091         }
0092
0093         m_alpha.resize(n + 2);
0094
0095         // Begin row reduction:
0096         std::vector<Real> rhs(n + 2, std::numeric_limits<Real>::quiet_NaN());
0097         std::vector<Real> super_diagonal(n + 2, std::numeric_limits<Real>::quiet_NaN());
0098
0099         rhs[0] = -a;
0100         rhs[rhs.size() - 1] = b;
0101
0102         super_diagonal[0] = 0;
0103
0104         for(size_t i = 1; i < rhs.size() - 1; ++i) {
0105             rhs[i] = 8*y[i - 1];
0106             super_diagonal[i] = 1;
0107         }
0108
0109         // Patch up 5-diagonal problem:
0110         rhs[1] = (rhs[1] - rhs[0])/6;
0111         super_diagonal[1] = Real(1)/Real(3);
0112         // First two rows are now:
0113         // 1 0 -1 | -2hy0'
0114         // 0 1 1/3| (8y0+2hy0')/6
0115
0116
0117         // Start traditional tridiagonal row reduction:
0118         for (size_t i = 2; i < rhs.size() - 1; ++i) {
0119             Real diagonal = 6 - super_diagonal[i - 1];
0120             rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
0121             super_diagonal[i] /= diagonal;
0122         }
0123
0124         //  1 sd[n-1] 0     | rhs[n-1]
0125         //  0 1       sd[n] | rhs[n]
0126         // -1 0       1     | rhs[n+1]
0127
0128         rhs[n+1] = rhs[n+1] + rhs[n-1];
0129         Real bottom_subdiagonal = super_diagonal[n-1];
0130
0131         // We're here:
0132         //  1 sd[n-1] 0     | rhs[n-1]
0133         //  0 1       sd[n] | rhs[n]
0134         //  0 bs      1     | rhs[n+1]
0135
0136         rhs[n+1] = (rhs[n+1]-bottom_subdiagonal*rhs[n])/(1-bottom_subdiagonal*super_diagonal[n]);
0137
0138         m_alpha[n+1] = rhs[n+1];
0139         for (size_t i = n; i > 0; --i) {
0140             m_alpha[i] = rhs[i] - m_alpha[i+1]*super_diagonal[i];
0141         }
0142         m_alpha[0] = m_alpha[2] + rhs[0];
0143     }
0144
0145     Real operator()(Real t) const {
0146         if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
0147             const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
0148             throw std::domain_error(err_msg);
0149         }
0150         // Let k, gamma be defined via t = t0 + kh + gamma * h.
0151         // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
0152         // j_min = ceil((t-t0)/h - 1/2)
0153         // j_max = floor(t-t0)/h + 5/2)
0154         using std::floor;
0155         using std::ceil;
0156         Real x = (t-m_t0)*m_inv_h;
0157         auto j_min = static_cast<size_t>(ceil(x - Real(1)/Real(2)));
0158         auto j_max = static_cast<size_t>(ceil(x + Real(5)/Real(2)));
0159         if (j_max >= m_alpha.size()) {
0160             j_max = m_alpha.size() - 1;
0161         }
0162
0163         Real y = 0;
0164         x += 1;
0165         for (size_t j = j_min; j <= j_max; ++j) {
0166             y += m_alpha[j]*detail::b2_spline(x - j);
0167         }
0168         return y;
0169     }
0170
0171     Real prime(Real t) const {
0172         if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) {
0173             const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work.";
0174             throw std::domain_error(err_msg);
0175         }
0176         // Let k, gamma be defined via t = t0 + kh + gamma * h.
0177         // Now find all j: |k-j+1+gamma|< 3/2, or, in other words
0178         // j_min = ceil((t-t0)/h - 1/2)
0179         // j_max = floor(t-t0)/h + 5/2)
0180         using std::floor;
0181         using std::ceil;
0182         Real x = (t-m_t0)*m_inv_h;
0183         auto j_min = static_cast<size_t>(ceil(x - Real(1)/Real(2)));
0184         auto j_max = static_cast<size_t>(ceil(x + Real(5)/Real(2)));
0185         if (j_max >= m_alpha.size()) {
0186             j_max = m_alpha.size() - 1;
0187         }
0188
0189         Real y = 0;
0190         x += 1;
0191         for (size_t j = j_min; j <= j_max; ++j) {
0192             y += m_alpha[j]*detail::b2_spline_prime(x - j);
0193         }
0194         return y*m_inv_h;
0195     }
0196
0197     Real t_max() const {
0198         return m_t0 + (m_alpha.size()-3)/m_inv_h;
0199     }
0200
0201 private:
0202     std::vector<Real> m_alpha;
0203     Real m_inv_h;
0204     Real m_t0;
0205 };
0206
0207 }}}}
0208 #endif