interpolators/detail/cardinal_cubic_b_spline_detail.hpp

0001 // Copyright Nick Thompson, 2017
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_CUBIC_B_SPLINE_DETAIL_HPP
0008 #define BOOST_MATH_INTERPOLATORS_CARDINAL_CUBIC_B_SPLINE_DETAIL_HPP
0009
0010 #include <limits>
0011 #include <cmath>
0012 #include <vector>
0013 #include <memory>
0014 #include <boost/math/constants/constants.hpp>
0015 #include <boost/math/special_functions/fpclassify.hpp>
0016 #include <boost/math/special_functions/trunc.hpp>
0017
0018 namespace boost{ namespace math{ namespace interpolators{ namespace detail{
0019
0020
0021 template <class Real>
0022 class cardinal_cubic_b_spline_imp
0023 {
0024 public:
0025     // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them.
0026     // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1).
0027     template <class BidiIterator>
0028     cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
0029                        Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
0030                        Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN());
0031
0032     Real operator()(Real x) const;
0033
0034     Real prime(Real x) const;
0035
0036     Real double_prime(Real x) const;
0037
0038 private:
0039     std::vector<Real> m_beta;
0040     Real m_h_inv;
0041     Real m_a;
0042     Real m_avg;
0043 };
0044
0045
0046
0047 template <class Real>
0048 Real b3_spline(Real x)
0049 {
0050     using std::abs;
0051     Real absx = abs(x);
0052     if (absx < 1)
0053     {
0054         Real y = 2 - absx;
0055         Real z = 1 - absx;
0056         return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z);
0057     }
0058     if (absx < 2)
0059     {
0060         Real y = 2 - absx;
0061         return boost::math::constants::sixth<Real>()*y*y*y;
0062     }
0063     return static_cast<Real>(0);
0064 }
0065
0066 template<class Real>
0067 Real b3_spline_prime(Real x)
0068 {
0069     if (x < 0)
0070     {
0071         return -b3_spline_prime(-x);
0072     }
0073
0074     if (x < 1)
0075     {
0076         return x*(3*boost::math::constants::half<Real>()*x - 2);
0077     }
0078     if (x < 2)
0079     {
0080         return -boost::math::constants::half<Real>()*(2 - x)*(2 - x);
0081     }
0082     return static_cast<Real>(0);
0083 }
0084
0085 template<class Real>
0086 Real b3_spline_double_prime(Real x)
0087 {
0088     if (x < 0)
0089     {
0090         return b3_spline_double_prime(-x);
0091     }
0092
0093     if (x < 1)
0094     {
0095         return 3*x - 2;
0096     }
0097     if (x < 2)
0098     {
0099         return (2 - x);
0100     }
0101     return static_cast<Real>(0);
0102 }
0103
0104
0105 template <class Real>
0106 template <class BidiIterator>
0107 cardinal_cubic_b_spline_imp<Real>::cardinal_cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size,
0108                                              Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0)
0109 {
0110     using boost::math::constants::third;
0111
0112     std::size_t length = end_p - f;
0113
0114     if (length < 5)
0115     {
0116         if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative))
0117         {
0118             throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n");
0119         }
0120         if (length < 3)
0121         {
0122             throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n");
0123         }
0124     }
0125
0126     if (boost::math::isnan(left_endpoint))
0127     {
0128         throw std::logic_error("Left endpoint is NAN; this is disallowed.\n");
0129     }
0130     if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)())
0131     {
0132         throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n");
0133     }
0134     if (step_size <= 0)
0135     {
0136         throw std::logic_error("The step size must be strictly > 0.\n");
0137     }
0138     // Storing the inverse of the stepsize does provide a measurable speedup.
0139     // It's not huge, but nonetheless worthwhile.
0140     m_h_inv = 1/step_size;
0141
0142     // Following Kress's notation, s'(a) = a1, s'(b) = b1
0143     Real a1 = left_endpoint_derivative;
0144     // See the finite-difference table on Wikipedia for reference on how
0145     // to construct high-order estimates for one-sided derivatives:
0146     // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference
0147     // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method.
0148     if (boost::math::isnan(a1))
0149     {
0150         // For simple functions (linear, quadratic, so on)
0151         // almost all the error comes from derivative estimation.
0152         // This does pairwise summation which gives us another digit of accuracy over naive summation.
0153         Real t0 = 4*(f[1] + third<Real>()*f[3]);
0154         Real t1 = -(25*third<Real>()*f[0] + f[4])/4  - 3*f[2];
0155         a1 = m_h_inv*(t0 + t1);
0156     }
0157
0158     Real b1 = right_endpoint_derivative;
0159     if (boost::math::isnan(b1))
0160     {
0161         size_t n = length - 1;
0162         Real t0 = -4*(f[n - 1] + third<Real>()*f[n - 3]);
0163         Real t1 = (25*third<Real>()*f[n] + f[n - 4])/4  + 3*f[n - 2];
0164
0165         b1 = m_h_inv*(t0 + t1);
0166     }
0167
0168     // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h )
0169     // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy.
0170     m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN());
0171
0172     // Since the splines have compact support, they decay to zero very fast outside the endpoints.
0173     // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the
0174     // boundary [a,b] without massive error.
0175     // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average.
0176     // This algorithm for computing the average is recommended in
0177     // http://www.heikohoffmann.de/htmlthesis/node134.html
0178     Real t = 1;
0179     for (size_t i = 0; i < length; ++i)
0180     {
0181         if (boost::math::isnan(f[i]))
0182         {
0183             std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n";
0184             throw std::logic_error(err);
0185         }
0186         m_avg += (f[i] - m_avg) / t;
0187         t += 1;
0188     }
0189
0190
0191     // Now we must solve an almost-tridiagonal system, which requires O(N) operations.
0192     // There are, in fact 5 diagonals, but they only differ from zero on the first and last row,
0193     // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows.
0194     // See Kress, equations 8.41
0195     // The "tridiagonal" matrix is:
0196     // 1  0 -1
0197     // 1  4  1
0198     //    1  4  1
0199     //       1  4  1
0200     //          ....
0201     //          1  4  1
0202     //          1  0 -1
0203     // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good.
0204     std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN());
0205     std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN());
0206
0207     rhs[0] = -2*step_size*a1;
0208     rhs[rhs.size() - 1] = -2*step_size*b1;
0209
0210     super_diagonal[0] = 0;
0211
0212     for(size_t i = 1; i < rhs.size() - 1; ++i)
0213     {
0214         rhs[i] = 6*(f[i - 1] - m_avg);
0215         super_diagonal[i] = 1;
0216     }
0217
0218
0219     // One step of row reduction on the first row to patch up the 5-diagonal problem:
0220     // 1 0 -1 | r0
0221     // 1 4 1  | r1
0222     // mapsto:
0223     // 1 0 -1 | r0
0224     // 0 4 2  | r1 - r0
0225     // mapsto
0226     // 1 0 -1 | r0
0227     // 0 1 1/2| (r1 - r0)/4
0228     super_diagonal[1] = static_cast<Real>(0.5);
0229     rhs[1] = (rhs[1] - rhs[0])/4;
0230
0231     // Now do a tridiagonal row reduction the standard way, until just before the last row:
0232     for (size_t i = 2; i < rhs.size() - 1; ++i)
0233     {
0234         Real diagonal = 4 - super_diagonal[i - 1];
0235         rhs[i] = (rhs[i] - rhs[i - 1])/diagonal;
0236         super_diagonal[i] /= diagonal;
0237     }
0238
0239     // Now the last row, which is in the form
0240     // 1 sd[n-3] 0      | rhs[n-3]
0241     // 0  1     sd[n-2] | rhs[n-2]
0242     // 1  0     -1      | rhs[n-1]
0243     Real final_subdiag = -super_diagonal[rhs.size() - 3];
0244     rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag;
0245     Real final_diag = -1/final_subdiag;
0246     // Now we're here:
0247     // 1 sd[n-3] 0         | rhs[n-3]
0248     // 0  1     sd[n-2]    | rhs[n-2]
0249     // 0  1     final_diag | (rhs[n-1] - rhs[n-3])/diag
0250
0251     final_diag = final_diag - super_diagonal[rhs.size() - 2];
0252     rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2];
0253
0254
0255     // Back substitutions:
0256     m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag;
0257     for(size_t i = rhs.size() - 2; i > 0; --i)
0258     {
0259         m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1];
0260     }
0261     m_beta[0] = m_beta[2] + rhs[0];
0262 }
0263
0264 template<class Real>
0265 Real cardinal_cubic_b_spline_imp<Real>::operator()(Real x) const
0266 {
0267     // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms,
0268     // just the (at most 5) whose support overlaps the argument.
0269     Real z = m_avg;
0270     Real t = m_h_inv*(x - m_a) + 1;
0271
0272     using std::max;
0273     using std::min;
0274     using std::ceil;
0275     using std::floor;
0276
0277     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
0278     size_t k_max = static_cast<size_t>((max)((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), 0l));
0279
0280     for (size_t k = k_min; k <= k_max; ++k)
0281     {
0282         z += m_beta[k]*b3_spline(t - k);
0283     }
0284
0285     return z;
0286 }
0287
0288 template<class Real>
0289 Real cardinal_cubic_b_spline_imp<Real>::prime(Real x) const
0290 {
0291     Real z = 0;
0292     Real t = m_h_inv*(x - m_a) + 1;
0293
0294     using std::max;
0295     using std::min;
0296     using std::ceil;
0297     using std::floor;
0298
0299     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
0300     size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))));
0301
0302     for (size_t k = k_min; k <= k_max; ++k)
0303     {
0304         z += m_beta[k]*b3_spline_prime(t - k);
0305     }
0306     return z*m_h_inv;
0307 }
0308
0309 template<class Real>
0310 Real cardinal_cubic_b_spline_imp<Real>::double_prime(Real x) const
0311 {
0312     Real z = 0;
0313     Real t = m_h_inv*(x - m_a) + 1;
0314
0315     using std::max;
0316     using std::min;
0317     using std::ceil;
0318     using std::floor;
0319
0320     size_t k_min = static_cast<size_t>((max)(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))));
0321     size_t k_max = static_cast<size_t>((min)(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))));
0322
0323     for (size_t k = k_min; k <= k_max; ++k)
0324     {
0325         z += m_beta[k]*b3_spline_double_prime(t - k);
0326     }
0327     return z*m_h_inv*m_h_inv;
0328 }
0329
0330 }}}}
0331 #endif