From fee20ab9325fcbf04d66987dbef4e75a9859d1bf Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Thu, 23 Feb 2017 18:21:06 -0600 Subject: [PATCH] Given a function f, known at evenly spaced samples y_j = f(a + jh), this function constructs an interpolant using compactly supported cubic b splines. The advantage of using splines of compact support over traditional cubic splines is that compact support makes the splines well-conditioned. The interpolant is constructed in O(N) time and can be evaluated in constant time. Its error is O(h^4), and obeys the interpolating condition s(x_j) = f(x_j) for all samples. In addition, f' can be estimated from s', albeit with lower accuracy. This routine is cppcheck clean, and is clean under AddressSanitizer and MemorySanitizer. --- doc/interpolators/cubic_b_spline.qbk | 53 ++++ .../math/interpolators/cubic_b_spline.hpp | 294 ++++++++++++++++++ test/test_cubic_b_spline.cpp | 217 +++++++++++++ 3 files changed, 564 insertions(+) create mode 100644 doc/interpolators/cubic_b_spline.qbk create mode 100644 include/boost/math/interpolators/cubic_b_spline.hpp create mode 100644 test/test_cubic_b_spline.cpp diff --git a/doc/interpolators/cubic_b_spline.qbk b/doc/interpolators/cubic_b_spline.qbk new file mode 100644 index 000000000..862975b89 --- /dev/null +++ b/doc/interpolators/cubic_b_spline.qbk @@ -0,0 +1,53 @@ +[/ +Copyright (c) 2017 Nick Thompson +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) +] + +[section:cubic_b Cubic B-spline interpolation] + +[section:cubic_b_over Overview of Cubic B-Spline Interpolation] + +The cubic b spline routine provided by boost allows fast and accurate interpolation of a function which is known at equally spaced points. +The cubic b-spline interpolation is preferable to traditional cubic spline interpolation as the global support of cubic polynomials makes the interpolation ill-conditioned. + +There are many use cases for interpolating a function at equally spaced points. +One particularly important example is solving ODE's whose coefficients depend on data determined from experiment or numerically. +Since most ODE steppers are adaptive, they must be able to sample the coefficients at arbitrary points; not just at the points we know the values of our function. + +The routine must be provided the following arguments: + +* A vector containing data +* The length of the vector +* The start of the functions domain +* The step size + +Optionally, you may provide two additional arguments to the routine, namely the derivative of the function at the left endpoint, and the derivative at the right endpoint. +If you do not provide these arguments, the routine will estimate them using one-sided finite-difference formulas. +An example of a valid call is + + __boost::math::cubic_b_spline(f.data(), f.size(), t0, h); + +Mathematically, the function we are trying to interpolate has the known values ['f[j] = f(jh+t[sub 0])'] + + +[endsect] + +[section:cubic_b_spline Header] + +[heading Accuracy] + +By far the greatest source of error in this routine is the estimation of the derivative at the endpoint. +If you know the endpoint derivatives, you should certainly provide them; however, in the vast majority of cases, +this is unknown, and the routine does a reasonable job trying to figure it out. + +[heading Testing] + +Since the interpolant obeys ['s(x[sub j]) = f(x[sub j])'] at all interpolation points, +the tests generate random data and evaluate the interpolant at the interpolation points, +validating that equality with the data holds. + +In addition, constant, linear, and quadratic functions are interpolated to ensure that the interpolant behaves as expected. + +[endsect] diff --git a/include/boost/math/interpolators/cubic_b_spline.hpp b/include/boost/math/interpolators/cubic_b_spline.hpp new file mode 100644 index 000000000..809334e6f --- /dev/null +++ b/include/boost/math/interpolators/cubic_b_spline.hpp @@ -0,0 +1,294 @@ +// 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) + +// This implements the compactly supported cubic b spline algorithm described in +// "Numerical Analsis, Graduate Texts in Mathematics 181", by Rainer Kress, section 8.3. +// Splines of compact support are faster to evaluate and are better conditioned than classical cubic splines. + +// Let f be the function we are trying to interpolate, and s be the interpolating spline. +// The routine constructs the interpolant in O(N) time, and evaluating s at a point takes constant time. +// The order of accuracy depends on the regularity of the f, however, assuming f is +// four-times continuously differentiable, the error is of O(h^4). +// In addition, we can differentiate the spline and obtain a good interpolant for f'. +// The main restriction of this method is that the samples of f must be evenly spaced. +// Interpolators for non-evenly sampled data are forthcoming. +// Properties: +// - s(x_j) = f(x_j) +// - All cubic polynomials interpolated exactly + +#ifndef BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP +#define BOOST_MATH_INTERPOLATORS_CUBIC_B_SPLINE_HPP + +#include +#include +#include +#include +#include + +namespace boost{ namespace math{ + +using boost::math::constants::third; +using boost::math::constants::half; + +template +Real b3_spline(Real x) +{ + auto sixth = half()*third(); + using std::abs; + auto absx = abs(x); + if (absx < 1) + { + auto y = 2 - absx; + auto z = 1 - absx; + return sixth*(y*y*y - 4*z*z*z); + } + if (absx < 2) + { + auto y = 2 - absx; + return sixth*y*y*y; + } + return (Real) 0; +} + +template +Real b3_spline_prime(Real x) +{ + if ( x < 0 ) + { + return -b3_spline_prime(-x); + } + + if (x < 1) + { + return x*(3*half()*x - 2); + } + if (x < 2) + { + return -half()*(2 - x)*(2 - x); + } + return (Real) 0; +} + +template +class cubic_b_spline +{ +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). + cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size, + Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), + Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()); + + Real interpolate_at(Real x) const; + + Real interpolate_derivative(Real x) const; + +private: + std::vector m_beta; + Real m_h_inv; + Real m_a; + Real m_avg; + +}; + + +template +cubic_b_spline::cubic_b_spline(const Real* const f, size_t length, Real left_endpoint, Real step_size, + Real left_endpoint_derivative, Real right_endpoint_derivative) : m_a(left_endpoint), m_avg(0) +{ + 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 (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. + auto t0 = 4*(f[1] + third()*f[3]); + auto t1 = -(25*third()*f[0] + f[4])/4 - 3*f[2]; + a1 = m_h_inv*(t0 + t1); + } + + Real b1 = right_endpoint_derivative; + if (boost::math::isnan(b1)) + { + auto n = length - 1; + auto t0 = 4*(f[n-3] + third()*f[n - 1]); + auto t1 = -(25*third()*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::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. + Real compensator = 0; + 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); + } + // Kahan summation: + auto y = f[i] - compensator; + auto t = m_avg + y; + auto diff = t - m_avg; + compensator = diff - y; + m_avg = t; + } + m_avg /= length; + + // 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 rhs(length + 2, std::numeric_limits::quiet_NaN()); + std::vector super_diagonal(length + 2, std::numeric_limits::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) + { + auto 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] + auto final_subdiag = -super_diagonal[rhs.size() - 3]; + rhs[rhs.size() - 1] = (rhs[rhs.size() - 1] - rhs[rhs.size() - 3])/final_subdiag; + auto 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 +Real cubic_b_spline::interpolate_at(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. + auto z = m_avg; + auto t = m_h_inv*(x - m_a) + 1; + + using std::max; + using std::min; + using std::ceil; + + size_t k_min = (size_t) max(static_cast(0), static_cast(ceil(t - 2))); + size_t k_max = (size_t) min(static_cast(m_beta.size() - 1), static_cast(floor(t + 2))); + + for (size_t k = k_min; k <= k_max; ++k) + { + z += m_beta[k]*b3_spline(t - k); + } + + return z; +} + +template +Real cubic_b_spline::interpolate_derivative(Real x) const +{ + Real z = 0; + auto t = m_h_inv*(x - m_a) + 1; + + using std::max; + using std::min; + using std::ceil; + + size_t k_min = (size_t) max(static_cast(0), static_cast(ceil(t - 2))); + size_t k_max = (size_t) min(static_cast(m_beta.size() - 1), static_cast(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 diff --git a/test/test_cubic_b_spline.cpp b/test/test_cubic_b_spline.cpp new file mode 100644 index 000000000..a3f84bf8f --- /dev/null +++ b/test/test_cubic_b_spline.cpp @@ -0,0 +1,217 @@ +#define BOOST_TEST_MODULE test_cubic_b_spline + +#include +#include +#include +#include +#include +#include + + + +using boost::math::constants::third; +using boost::math::constants::half; + +template +void test_b3_spline() +{ + std::cout << "Testing evaluation of spline basis functions on type " << boost::typeindex::type_id().pretty_name() << "\n"; + // Outside the support: + auto eps = std::numeric_limits::epsilon(); + BOOST_CHECK_SMALL(boost::math::b3_spline(2.5), (Real) 0); + BOOST_CHECK_SMALL(boost::math::b3_spline(-2.5), (Real) 0); + BOOST_CHECK_SMALL(boost::math::b3_spline_prime(2.5), (Real) 0); + BOOST_CHECK_SMALL(boost::math::b3_spline_prime(-2.5), (Real) 0); + + // On the boundary of support: + BOOST_CHECK_SMALL(boost::math::b3_spline(2), (Real) 0); + BOOST_CHECK_SMALL(boost::math::b3_spline(-2), (Real) 0); + BOOST_CHECK_SMALL(boost::math::b3_spline_prime(2), (Real) 0); + BOOST_CHECK_SMALL(boost::math::b3_spline_prime(-2), (Real) 0); + + // Special values: + BOOST_CHECK_CLOSE(boost::math::b3_spline(-1), third()*half(), eps); + BOOST_CHECK_CLOSE(boost::math::b3_spline( 1), third()*half(), eps); + BOOST_CHECK_CLOSE(boost::math::b3_spline(0), 2*third(), eps); + + BOOST_CHECK_CLOSE(boost::math::b3_spline_prime(-1), half(), eps); + BOOST_CHECK_CLOSE(boost::math::b3_spline_prime( 1), -half(), eps); + BOOST_CHECK_SMALL(boost::math::b3_spline_prime(0), eps); + + // Properties: B3 is an even function, B3' is an odd function. + for (size_t i = 1; i < 200; ++i) + { + auto arg = i*0.01; + BOOST_CHECK_CLOSE(boost::math::b3_spline(arg), boost::math::b3_spline(arg), eps); + BOOST_CHECK_CLOSE(boost::math::b3_spline_prime(-arg), -boost::math::b3_spline_prime(arg), eps); + } + +} +/* + * This test ensures that the interpolant s(x_j) = f(x_j) at all grid points. + */ +template +void test_interpolation_condition() +{ + std::cout << "Testing interpolation condition for cubic b splines on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(1, 10); + boost::random::uniform_real_distribution step_distribution(0.001, 0.01); + std::vector v(5000); + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = dis(gen); + } + + Real step = step_distribution(gen); + Real a = dis(gen); + boost::math::cubic_b_spline spline(v.data(), v.size(), a, step); + + for (size_t i = 0; i < v.size(); ++i) + { + auto y = spline.interpolate_at(i*step + a); + // This seems like a very large tolerance, but I don't know of any other interpolators + // that will be able to do much better on random data. + BOOST_CHECK_CLOSE(y, v[i], 10000000*std::numeric_limits::epsilon()); + } + +} + +template +void test_constant_function() +{ + std::cout << "Testing that constants are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(-100, 100); + boost::random::uniform_real_distribution step_distribution(0.0001, 0.1); + std::vector v(500); + Real constant = dis(gen); + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = constant; + } + + Real step = step_distribution(gen); + Real a = dis(gen); + boost::math::cubic_b_spline spline(v.data(), v.size(), a, step); + + for (size_t i = 0; i < v.size(); ++i) + { + auto y = spline.interpolate_at(i*step + a); + BOOST_CHECK_CLOSE(y, v[i], 10*std::numeric_limits::epsilon()); + auto y_prime = spline.interpolate_derivative(i*step + a); + BOOST_CHECK_SMALL(y_prime, 5000*std::numeric_limits::epsilon()); + } + + // Test that correctly specified left and right-derivatives work properly: + spline = boost::math::cubic_b_spline(v.data(), v.size(), a, step, 0, 0); + + for (size_t i = 0; i < v.size(); ++i) + { + auto y = spline.interpolate_at(i*step + a); + BOOST_CHECK_CLOSE(y, v[i], std::numeric_limits::epsilon()); + auto y_prime = spline.interpolate_derivative(i*step + a); + BOOST_CHECK_SMALL(y_prime, std::numeric_limits::epsilon()); + } +} + + +template +void test_affine_function() +{ + std::cout << "Testing that affine functions are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(-100, 100); + boost::random::uniform_real_distribution step_distribution(0.0001, 0.01); + std::vector v(500); + Real a = dis(gen); + Real b = dis(gen); + Real step = step_distribution(gen); + Real x0 = dis(gen); + + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = a*(i*step - x0) + b; + } + + + boost::math::cubic_b_spline spline(v.data(), v.size(), x0, step); + + for (size_t i = 0; i < v.size(); ++i) + { + auto y = spline.interpolate_at(i*step + x0); + BOOST_CHECK_CLOSE(y, v[i], 10000*std::numeric_limits::epsilon()); + auto y_prime = spline.interpolate_derivative(i*step + x0); + BOOST_CHECK_CLOSE(y_prime, a, 10000000*std::numeric_limits::epsilon()); + } + + // Test that correctly specified left and right-derivatives work properly: + spline = boost::math::cubic_b_spline(v.data(), v.size(), x0, step, a, a); + + for (size_t i = 0; i < v.size(); ++i) + { + auto y = spline.interpolate_at(i*step + x0); + BOOST_CHECK_CLOSE(y, v[i], 10000*std::numeric_limits::epsilon()); + auto y_prime = spline.interpolate_derivative(i*step + x0); + BOOST_CHECK_CLOSE(y_prime, a, 10000000*std::numeric_limits::epsilon()); + } +} + +template +void test_quadratic_function() +{ + std::cout << "Testing that quadratic functions are interpolated correctly by cubic b splines on type " << boost::typeindex::type_id().pretty_name() << "\n"; + std::random_device rd; + std::mt19937 gen(rd()); + boost::random::uniform_real_distribution dis(1, 5); + boost::random::uniform_real_distribution step_distribution(0.0001, 1); + std::vector v(500); + Real a = dis(gen); + Real b = dis(gen); + Real c = dis(gen); + Real step = step_distribution(gen); + Real x0 = dis(gen); + + for (size_t i = 0; i < v.size(); ++i) + { + v[i] = a*(i*step - x0)*(i*step - x0) + b*(i*step - x0) + c; + } + + boost::math::cubic_b_spline spline(v.data(), v.size(), x0, step); + + for (size_t i = 0; i < v.size(); ++i) + { + auto y = spline.interpolate_at(i*step + x0); + BOOST_CHECK_CLOSE(y, v[i], 100000000*std::numeric_limits::epsilon()); + auto y_prime = spline.interpolate_derivative(i*step + x0); + BOOST_CHECK_CLOSE(y_prime, 2*a*(i*step-x0) + b, 2.0); + } +} + + +BOOST_AUTO_TEST_CASE(test_cubic_b_spline) +{ + test_b3_spline(); + test_b3_spline(); + test_b3_spline(); + + test_interpolation_condition(); + test_interpolation_condition(); + test_interpolation_condition(); + + test_constant_function(); + test_constant_function(); + test_constant_function(); + + test_affine_function(); + test_affine_function(); + test_affine_function(); + + test_quadratic_function(); + test_quadratic_function(); + test_quadratic_function(); + +}