boost/math/interpolators/detail/cardinal_quadratic_b_spline_detail.hpp

   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