diff options
Diffstat (limited to 'boost/math/interpolators/detail/cubic_b_spline_detail.hpp')
| -rw-r--r-- | boost/math/interpolators/detail/cubic_b_spline_detail.hpp | 287 |
1 files changed, 287 insertions, 0 deletions
diff --git a/boost/math/interpolators/detail/cubic_b_spline_detail.hpp b/boost/math/interpolators/detail/cubic_b_spline_detail.hpp new file mode 100644 index 0000000000..f7b2d6cd29 --- /dev/null +++ b/boost/math/interpolators/detail/cubic_b_spline_detail.hpp @@ -0,0 +1,287 @@ +// Copyright Nick Thompson, 2017 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef CUBIC_B_SPLINE_DETAIL_HPP +#define CUBIC_B_SPLINE_DETAIL_HPP + +#include <limits> +#include <cmath> +#include <vector> +#include <memory> +#include <boost/math/constants/constants.hpp> +#include <boost/math/special_functions/fpclassify.hpp> + +namespace boost{ namespace math{ namespace detail{ + + +template <class Real> +class cubic_b_spline_imp +{ +public: + // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them. + // f[0] = f(a), f[length -1] = b, step_size = (b - a)/(length -1). + template <class BidiIterator> + cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size, + Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(), + Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()); + + Real operator()(Real x) const; + + Real prime(Real x) const; + +private: + std::vector<Real> m_beta; + Real m_h_inv; + Real m_a; + Real m_avg; +}; + + + +template <class Real> +Real b3_spline(Real x) +{ + using std::abs; + Real absx = abs(x); + if (absx < 1) + { + Real y = 2 - absx; + Real z = 1 - absx; + return boost::math::constants::sixth<Real>()*(y*y*y - 4*z*z*z); + } + if (absx < 2) + { + Real y = 2 - absx; + return boost::math::constants::sixth<Real>()*y*y*y; + } + return (Real) 0; +} + +template<class Real> +Real b3_spline_prime(Real x) +{ + if (x < 0) + { + return -b3_spline_prime(-x); + } + + if (x < 1) + { + return x*(3*boost::math::constants::half<Real>()*x - 2); + } + if (x < 2) + { + return -boost::math::constants::half<Real>()*(2 - x)*(2 - x); + } + return (Real) 0; +} + + +template <class Real> +template <class BidiIterator> +cubic_b_spline_imp<Real>::cubic_b_spline_imp(BidiIterator f, BidiIterator end_p, Real left_endpoint, Real step_size, + Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0) +{ + using boost::math::constants::third; + + std::size_t length = end_p - f; + + if (length < 5) + { + if (boost::math::isnan(left_endpoint_derivative) || boost::math::isnan(right_endpoint_derivative)) + { + throw std::logic_error("Interpolation using a cubic b spline with derivatives estimated at the endpoints requires at least 5 points.\n"); + } + if (length < 3) + { + throw std::logic_error("Interpolation using a cubic b spline requires at least 3 points.\n"); + } + } + + if (boost::math::isnan(left_endpoint)) + { + throw std::logic_error("Left endpoint is NAN; this is disallowed.\n"); + } + if (left_endpoint + length*step_size >= (std::numeric_limits<Real>::max)()) + { + throw std::logic_error("Right endpoint overflows the maximum representable number of the specified precision.\n"); + } + if (step_size <= 0) + { + throw std::logic_error("The step size must be strictly > 0.\n"); + } + // Storing the inverse of the stepsize does provide a measurable speedup. + // It's not huge, but nonetheless worthwhile. + m_h_inv = 1/step_size; + + // Following Kress's notation, s'(a) = a1, s'(b) = b1 + Real a1 = left_endpoint_derivative; + // See the finite-difference table on Wikipedia for reference on how + // to construct high-order estimates for one-sided derivatives: + // https://en.wikipedia.org/wiki/Finite_difference_coefficient#Forward_and_backward_finite_difference + // Here, we estimate then to O(h^4), as that is the maximum accuracy we could obtain from this method. + if (boost::math::isnan(a1)) + { + // For simple functions (linear, quadratic, so on) + // almost all the error comes from derivative estimation. + // This does pairwise summation which gives us another digit of accuracy over naive summation. + Real t0 = 4*(f[1] + third<Real>()*f[3]); + Real t1 = -(25*third<Real>()*f[0] + f[4])/4 - 3*f[2]; + a1 = m_h_inv*(t0 + t1); + } + + Real b1 = right_endpoint_derivative; + if (boost::math::isnan(b1)) + { + size_t n = length - 1; + Real t0 = 4*(f[n-3] + third<Real>()*f[n - 1]); + Real t1 = -(25*third<Real>()*f[n - 4] + f[n])/4 - 3*f[n - 2]; + + b1 = m_h_inv*(t0 + t1); + } + + // s(x) = \sum \alpha_i B_{3}( (x- x_i - a)/h ) + // Of course we must reindex from Kress's notation, since he uses negative indices which make C++ unhappy. + m_beta.resize(length + 2, std::numeric_limits<Real>::quiet_NaN()); + + // Since the splines have compact support, they decay to zero very fast outside the endpoints. + // This is often very annoying; we'd like to evaluate the interpolant a little bit outside the + // boundary [a,b] without massive error. + // A simple way to deal with this is just to subtract the DC component off the signal, so we need the average. + // This algorithm for computing the average is recommended in + // http://www.heikohoffmann.de/htmlthesis/node134.html + Real t = 1; + for (size_t i = 0; i < length; ++i) + { + if (boost::math::isnan(f[i])) + { + std::string err = "This function you are trying to interpolate is a nan at index " + std::to_string(i) + "\n"; + throw std::logic_error(err); + } + m_avg += (f[i] - m_avg) / t; + t += 1; + } + + + // Now we must solve an almost-tridiagonal system, which requires O(N) operations. + // There are, in fact 5 diagonals, but they only differ from zero on the first and last row, + // so we can patch up the tridiagonal row reduction algorithm to deal with two special rows. + // See Kress, equations 8.41 + // The the "tridiagonal" matrix is: + // 1 0 -1 + // 1 4 1 + // 1 4 1 + // 1 4 1 + // .... + // 1 4 1 + // 1 0 -1 + // Numerical estimate indicate that as N->Infinity, cond(A) -> 6.9, so this matrix is good. + std::vector<Real> rhs(length + 2, std::numeric_limits<Real>::quiet_NaN()); + std::vector<Real> super_diagonal(length + 2, std::numeric_limits<Real>::quiet_NaN()); + + rhs[0] = -2*step_size*a1; + rhs[rhs.size() - 1] = -2*step_size*b1; + + super_diagonal[0] = 0; + + for(size_t i = 1; i < rhs.size() - 1; ++i) + { + rhs[i] = 6*(f[i - 1] - m_avg); + super_diagonal[i] = 1; + } + + + // One step of row reduction on the first row to patch up the 5-diagonal problem: + // 1 0 -1 | r0 + // 1 4 1 | r1 + // mapsto: + // 1 0 -1 | r0 + // 0 4 2 | r1 - r0 + // mapsto + // 1 0 -1 | r0 + // 0 1 1/2| (r1 - r0)/4 + super_diagonal[1] = 0.5; + rhs[1] = (rhs[1] - rhs[0])/4; + + // Now do a tridiagonal row reduction the standard way, until just before the last row: + for (size_t i = 2; i < rhs.size() - 1; ++i) + { + Real diagonal = 4 - super_diagonal[i - 1]; + rhs[i] = (rhs[i] - rhs[i - 1])/diagonal; + super_diagonal[i] /= diagonal; + } + + // Now the last row, which is in the form + // 1 sd[n-3] 0 | rhs[n-3] + // 0 1 sd[n-2] | rhs[n-2] + // 1 0 -1 | rhs[n-1] + Real final_subdiag = -super_diagonal[rhs.size() - 3]; + rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag; + Real final_diag = -1/final_subdiag; + // Now we're here: + // 1 sd[n-3] 0 | rhs[n-3] + // 0 1 sd[n-2] | rhs[n-2] + // 0 1 final_diag | (rhs[n-1] - rhs[n-3])/diag + + final_diag = final_diag - super_diagonal[rhs.size() - 2]; + rhs[rhs.size() - 1] = rhs[rhs.size() - 1] - rhs[rhs.size() - 2]; + + + // Back substitutions: + m_beta[rhs.size() - 1] = rhs[rhs.size() - 1]/final_diag; + for(size_t i = rhs.size() - 2; i > 0; --i) + { + m_beta[i] = rhs[i] - super_diagonal[i]*m_beta[i + 1]; + } + m_beta[0] = m_beta[2] + rhs[0]; +} + +template<class Real> +Real cubic_b_spline_imp<Real>::operator()(Real x) const +{ + // See Kress, 8.40: Since B3 has compact support, we don't have to sum over all terms, + // just the (at most 5) whose support overlaps the argument. + Real z = m_avg; + Real t = m_h_inv*(x - m_a) + 1; + + using std::max; + using std::min; + using std::ceil; + using std::floor; + + size_t k_min = (size_t) max(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))); + size_t k_max = (size_t) max(min(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))), (long) 0); + for (size_t k = k_min; k <= k_max; ++k) + { + z += m_beta[k]*b3_spline(t - k); + } + + return z; +} + +template<class Real> +Real cubic_b_spline_imp<Real>::prime(Real x) const +{ + Real z = 0; + Real t = m_h_inv*(x - m_a) + 1; + + using std::max; + using std::min; + using std::ceil; + using std::floor; + + size_t k_min = (size_t) max(static_cast<long>(0), boost::math::ltrunc(ceil(t - 2))); + size_t k_max = (size_t) min(static_cast<long>(m_beta.size() - 1), boost::math::ltrunc(floor(t + 2))); + + for (size_t k = k_min; k <= k_max; ++k) + { + z += m_beta[k]*b3_spline_prime(t - k); + } + return z*m_h_inv; +} + +}}} +#endif |