interpolators/detail/cardinal_quintic_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_QUINTIC_B_SPLINE_DETAIL_HPP
0008 #define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
0009 #include <cmath>
0010 #include <vector>
0011 #include <utility>
0012 #include <boost/math/special_functions/cardinal_b_spline.hpp>
0013
0014 namespace boost{ namespace math{ namespace interpolators{ namespace detail{
0015
0016
0017 template <class Real>
0018 class cardinal_quintic_b_spline_detail
0019 {
0020 public:
0021     // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
0022     // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1).
0023
0024     cardinal_quintic_b_spline_detail(const Real* const y,
0025                                      size_t n,
0026                                      Real t0 /* initial time, left endpoint */,
0027                                      Real h  /*spacing, stepsize*/,
0028                                      std::pair<Real, Real> left_endpoint_derivatives,
0029                                      std::pair<Real, Real> right_endpoint_derivatives)
0030     {
0031         static_assert(!std::is_integral<Real>::value, "The quintic B-spline interpolator only works with floating point types.");
0032         if (h <= 0) {
0033             throw std::logic_error("Spacing must be > 0.");
0034         }
0035         m_inv_h = 1/h;
0036         m_t0 = t0;
0037
0038         if (n < 8) {
0039             throw std::logic_error("The quintic B-spline interpolator requires at least 8 points.");
0040         }
0041
0042         using std::isnan;
0043         // This interpolator has error of order h^6, so the derivatives should be estimated with the same error.
0044         // See: https://en.wikipedia.org/wiki/Finite_difference_coefficient
0045         if (isnan(left_endpoint_derivatives.first)) {
0046             Real tmp = -49*y[0]/20 + 6*y[1] - 15*y[2]/2 + 20*y[3]/3 - 15*y[4]/4 + 6*y[5]/5 - y[6]/6;
0047             left_endpoint_derivatives.first = tmp/h;
0048         }
0049         if (isnan(right_endpoint_derivatives.first)) {
0050             Real tmp = 49*y[n-1]/20 - 6*y[n-2] + 15*y[n-3]/2 - 20*y[n-4]/3 + 15*y[n-5]/4 - 6*y[n-6]/5 + y[n-7]/6;
0051             right_endpoint_derivatives.first = tmp/h;
0052         }
0053         if(isnan(left_endpoint_derivatives.second)) {
0054             Real tmp = 469*y[0]/90 - 223*y[1]/10 + 879*y[2]/20 - 949*y[3]/18 + 41*y[4] - 201*y[5]/10 + 1019*y[6]/180 - 7*y[7]/10;
0055             left_endpoint_derivatives.second = tmp/(h*h);
0056         }
0057         if (isnan(right_endpoint_derivatives.second)) {
0058             Real tmp = 469*y[n-1]/90 - 223*y[n-2]/10 + 879*y[n-3]/20 - 949*y[n-4]/18 + 41*y[n-5] - 201*y[n-6]/10 + 1019*y[n-7]/180 - 7*y[n-8]/10;
0059             right_endpoint_derivatives.second = tmp/(h*h);
0060         }
0061
0062         // This is really challenging my mental limits on by-hand row reduction.
0063         // I debated bringing in a dependency on a sparse linear solver, but given that that would cause much agony for users I decided against it.
0064
0065         std::vector<Real> rhs(n+4);
0066         rhs[0] = 20*y[0] - 12*h*left_endpoint_derivatives.first +  2*h*h*left_endpoint_derivatives.second;
0067         rhs[1] = 60*y[0] - 12*h*left_endpoint_derivatives.first;
0068         for (size_t i = 2; i < n + 2; ++i) {
0069             rhs[i] = 120*y[i-2];
0070         }
0071         rhs[n+2] = 60*y[n-1] + 12*h*right_endpoint_derivatives.first;
0072         rhs[n+3] = 20*y[n-1] + 12*h*right_endpoint_derivatives.first +  2*h*h*right_endpoint_derivatives.second;
0073
0074         std::vector<Real> diagonal(n+4, 66);
0075         diagonal[0] = 1;
0076         diagonal[1] = 18;
0077         diagonal[n+2] = 18;
0078         diagonal[n+3] = 1;
0079
0080         std::vector<Real> first_superdiagonal(n+4, 26);
0081         first_superdiagonal[0] = 10;
0082         first_superdiagonal[1] = 33;
0083         first_superdiagonal[n+2] = 1;
0084         // There is one less superdiagonal than diagonal; make sure that if we read it, it shows up as a bug:
0085         first_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();
0086
0087         std::vector<Real> second_superdiagonal(n+4, 1);
0088         second_superdiagonal[0] = 9;
0089         second_superdiagonal[1] = 8;
0090         second_superdiagonal[n+2] = std::numeric_limits<Real>::quiet_NaN();
0091         second_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();
0092
0093         std::vector<Real> first_subdiagonal(n+4, 26);
0094         first_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
0095         first_subdiagonal[1] = 1;
0096         first_subdiagonal[n+2] = 33;
0097         first_subdiagonal[n+3] = 10;
0098
0099         std::vector<Real> second_subdiagonal(n+4, 1);
0100         second_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
0101         second_subdiagonal[1] = std::numeric_limits<Real>::quiet_NaN();
0102         second_subdiagonal[n+2] = 8;
0103         second_subdiagonal[n+3] = 9;
0104
0105
0106         for (size_t i = 0; i < n+2; ++i) {
0107             Real di = diagonal[i];
0108             diagonal[i] = 1;
0109             first_superdiagonal[i] /= di;
0110             second_superdiagonal[i] /= di;
0111             rhs[i] /= di;
0112
0113             // Eliminate first subdiagonal:
0114             Real nfsub = -first_subdiagonal[i+1];
0115             // Superfluous:
0116             first_subdiagonal[i+1] /= nfsub;
0117             // Not superfluous:
0118             diagonal[i+1] /= nfsub;
0119             first_superdiagonal[i+1] /= nfsub;
0120             second_superdiagonal[i+1] /= nfsub;
0121             rhs[i+1] /= nfsub;
0122
0123             diagonal[i+1] += first_superdiagonal[i];
0124             first_superdiagonal[i+1] += second_superdiagonal[i];
0125             rhs[i+1] += rhs[i];
0126             // Superfluous, but clarifying:
0127             first_subdiagonal[i+1] = 0;
0128
0129             // Eliminate second subdiagonal:
0130             Real nssub = -second_subdiagonal[i+2];
0131             first_subdiagonal[i+2] /= nssub;
0132             diagonal[i+2] /= nssub;
0133             first_superdiagonal[i+2] /= nssub;
0134             second_superdiagonal[i+2] /= nssub;
0135             rhs[i+2] /= nssub;
0136
0137             first_subdiagonal[i+2] += first_superdiagonal[i];
0138             diagonal[i+2] += second_superdiagonal[i];
0139             rhs[i+2] += rhs[i];
0140             // Superfluous, but clarifying:
0141             second_subdiagonal[i+2] = 0;
0142         }
0143
0144         // Eliminate last subdiagonal:
0145         Real dnp2 = diagonal[n+2];
0146         diagonal[n+2] = 1;
0147         first_superdiagonal[n+2] /= dnp2;
0148         rhs[n+2] /= dnp2;
0149         Real nfsubnp3 = -first_subdiagonal[n+3];
0150         diagonal[n+3] /= nfsubnp3;
0151         rhs[n+3] /= nfsubnp3;
0152
0153         diagonal[n+3] += first_superdiagonal[n+2];
0154         rhs[n+3] += rhs[n+2];
0155
0156         m_alpha.resize(n + 4, std::numeric_limits<Real>::quiet_NaN());
0157
0158         m_alpha[n+3] = rhs[n+3]/diagonal[n+3];
0159         m_alpha[n+2] = rhs[n+2] - first_superdiagonal[n+2]*m_alpha[n+3];
0160         for (int64_t i = int64_t(n+1); i >= 0; --i) {
0161             m_alpha[i] = rhs[i] - first_superdiagonal[i]*m_alpha[i+1] - second_superdiagonal[i]*m_alpha[i+2];
0162         }
0163
0164     }
0165
0166     Real operator()(Real t) const {
0167         using std::ceil;
0168         using std::floor;
0169         using boost::math::cardinal_b_spline;
0170         // tf = t0 + (n-1)*h
0171         // alpha.size() = n+4
0172         if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
0173             const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
0174             throw std::domain_error(err_msg);
0175         }
0176         Real x = (t-m_t0)*m_inv_h;
0177         // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5.
0178         // TODO: Zero pad m_alpha so that only the domain check is necessary.
0179         int64_t j_min = (std::max)(int64_t(0), int64_t(ceil(x-1)));
0180         int64_t j_max = (std::min)(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
0181         Real s = 0;
0182         for (int64_t j = j_min; j <= j_max; ++j) {
0183             // TODO: Use Cox 1972 to generate all integer translates of B5 simultaneously.
0184             s += m_alpha[j]*cardinal_b_spline<5, Real>(x - j + 2);
0185         }
0186         return s;
0187     }
0188
0189     Real prime(Real t) const {
0190         using std::ceil;
0191         using std::floor;
0192         using boost::math::cardinal_b_spline_prime;
0193         if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
0194             const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
0195             throw std::domain_error(err_msg);
0196         }
0197         Real x = (t-m_t0)*m_inv_h;
0198         // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
0199         int64_t j_min = (std::max)(int64_t(0), int64_t(ceil(x-1)));
0200         int64_t j_max = (std::min)(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
0201         Real s = 0;
0202         for (int64_t j = j_min; j <= j_max; ++j) {
0203             s += m_alpha[j]*cardinal_b_spline_prime<5, Real>(x - j + 2);
0204         }
0205         return s*m_inv_h;
0206
0207     }
0208
0209     Real double_prime(Real t) const {
0210         using std::ceil;
0211         using std::floor;
0212         using boost::math::cardinal_b_spline_double_prime;
0213         if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) {
0214             const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work.";
0215             throw std::domain_error(err_msg);
0216         }
0217         Real x = (t-m_t0)*m_inv_h;
0218         // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5
0219         int64_t j_min = (std::max)(int64_t(0), int64_t(ceil(x-1)));
0220         int64_t j_max = (std::min)(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) );
0221         Real s = 0;
0222         for (int64_t j = j_min; j <= j_max; ++j) {
0223             s += m_alpha[j]*cardinal_b_spline_double_prime<5, Real>(x - j + 2);
0224         }
0225         return s*m_inv_h*m_inv_h;
0226     }
0227
0228
0229     Real t_max() const {
0230         return m_t0 + (m_alpha.size()-5)/m_inv_h;
0231     }
0232
0233 private:
0234     std::vector<Real> m_alpha;
0235     Real m_inv_h;
0236     Real m_t0;
0237 };
0238
0239 }}}}
0240 #endif