From 73429a3df597936c17536546a9458d627e14e081 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Sun, 11 Aug 2019 11:55:59 -0400 Subject: [PATCH 1/7] Cardinal Quintic B-spline interpolator: First sketch. [CI SKIP] --- .../cardinal_quintic_b_spline.hpp | 62 +++++++++++++ .../cardinal_quintic_b_spline_detail.hpp | 89 +++++++++++++++++++ 2 files changed, 151 insertions(+) create mode 100644 include/boost/math/interpolators/cardinal_quintic_b_spline.hpp create mode 100644 include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp diff --git a/include/boost/math/interpolators/cardinal_quintic_b_spline.hpp b/include/boost/math/interpolators/cardinal_quintic_b_spline.hpp new file mode 100644 index 000000000..3d72865d9 --- /dev/null +++ b/include/boost/math/interpolators/cardinal_quintic_b_spline.hpp @@ -0,0 +1,62 @@ +// Copyright Nick Thompson, 2019 +// 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 BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_HPP +#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_HPP +#include +#include +#include + + +namespace boost{ namespace math{ namespace interpolators { + +template +class cardinal_quintic_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. + // y[0] = y(a), y[n - 1] = y(b), step_size = (b - a)/(n -1). + cardinal_quintic_b_spline(const Real* const y, + size_t n, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + std::pair left_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}, + std::pair right_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}) + : impl_(std::make_shared>(y, n, t0, h, left_endpoint_derivatives, right_endpoint_derivatives)) + {} + + // Oh the bizarre error messages if we template this on a RandomAccessContainer: + cardinal_quintic_b_spline(std::vector const & y, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + std::pair left_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}, + std::pair right_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limits::quiet_NaN()}) + : impl_(std::make_shared>(y.data(), y.size(), t0, h, left_endpoint_derivatives, right_endpoint_derivatives)) + {} + + + Real operator()(Real t) const { + return impl_->operator()(t); + } + + Real prime(Real t) const { + return impl_->prime(t); + } + + Real double_prime(Real t) const { + return impl_->double_prime(t); + } + + Real t_max() const { + return impl_->t_max(); + } + +private: + std::shared_ptr> impl_; +}; + +}}} +#endif diff --git a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp new file mode 100644 index 000000000..9c172a66b --- /dev/null +++ b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp @@ -0,0 +1,89 @@ +// Copyright Nick Thompson, 2019 +// 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 BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP +#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP +#include +#include + +namespace boost{ namespace math{ namespace interpolators{ namespace detail{ + + +template +class cardinal_quintic_b_spline_detail +{ +public: + // If you don't know the value of the derivative at the endpoints, leave them as nans and the routine will estimate them. + // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1). + + cardinal_quintic_b_spline_detail(const Real* const y, + size_t n, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + std::pair left_endpoint_derivatives, + std::pair right_endpoint_derivatives) + { + if (h <= 0) { + throw std::logic_error("Spacing must be > 0."); + } + m_inv_h = 1/h; + m_t0 = t0; + + if (n < 3) { + throw std::logic_error("The interpolator requires at least 3 points."); + } + + m_alpha.resize(n + 4); + std::vector rhs(n+4); + rhs[0] = 6*h*h*left_endpoint_derivatives.second; + rhs[1] = 24*h*left_endpoint_derivatives.first; + for (size_t i = 2; i < n + 2; ++i) { + rhs[i] = 120*y[i-2]; + } + rhs[n+2] = 24*h*right_endpoint_derivatives.first; + rhs[n+3] = 6*h*h*right_endpoint_derivatives.second; + + std::vector diagonal(n+4); + + } + + Real operator()(Real t) const { + if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + throw std::domain_error(err_msg); + } + return std::numeric_limits::quiet_NaN(); + } + + Real prime(Real t) const { + if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + throw std::domain_error(err_msg); + } + return std::numeric_limits::quiet_NaN(); + } + + Real double_prime(Real t) const { + if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + throw std::domain_error(err_msg); + } + return std::numeric_limits::quiet_NaN(); + } + + + Real t_max() const { + return m_t0 + (m_alpha.size()-3)/m_inv_h; + } + +private: + std::vector m_alpha; + Real m_inv_h; + Real m_t0; +}; + +}}}} +#endif From 1955699777a3e93ba06323b1e21a09a04b60916d Mon Sep 17 00:00:00 2001 From: NAThompson Date: Tue, 13 Aug 2019 09:33:54 -0400 Subject: [PATCH 2/7] Cardinal Quintic B-splines: Documentation. [CI SKIP] --- .../cardinal_quintic_b_spline.qbk | 83 +++++++++++++++++++ doc/math.qbk | 7 +- 2 files changed, 87 insertions(+), 3 deletions(-) create mode 100644 doc/interpolators/cardinal_quintic_b_spline.qbk diff --git a/doc/interpolators/cardinal_quintic_b_spline.qbk b/doc/interpolators/cardinal_quintic_b_spline.qbk new file mode 100644 index 000000000..085482671 --- /dev/null +++ b/doc/interpolators/cardinal_quintic_b_spline.qbk @@ -0,0 +1,83 @@ +[/ +Copyright (c) 2019 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:cardinal_quintic_b Cardinal Quintic B-spline interpolation] + +[heading Synopsis] +`` + #include +`` + + namespace boost{ namespace math{ namespace interpolators { + + template + class cardinal_quintic_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. + // y[0] = y(a), y[n - 1] = y(b), step_size = (b - a)/(n -1). + cardinal_quintic_b_spline(const Real* const y, + size_t n, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + std::pair left_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limit::quiet_NaN()}, + std::pair right_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limit::quiet_NaN()}) + + cardinal_quintic_b_spline(std::vector const & y, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + std::pair left_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limit::quiet_NaN()}, + std::pair right_endpoint_derivatives = {std::numeric_limits::quiet_NaN(), std::numeric_limit::quiet_NaN()}) + + Real operator()(Real t) const; + + Real prime(Real t) const; + + Real double_prime(Real t) const; + + }; + }}} + +[heading Cardinal Quintic B-Spline Interpolation] + +The cardinal quintic B-spline interpolator is very nearly the same as the cubic B-spline interpolator, +with the modification that the basis functions are constructed by convolving a box function with itself five times, +rather than three times as is done with the cubic B-spline. + +The basis functions of the quintic B-spline interpolator are more smooth than the cubic /B/-spline interpolator, and hence this is very useful for computing second derivatives. +For example, the second derivative of the cubic spline interpolator is a piecewise linear function, whereas the second derivative of the quintic /B/-spline is a cubic spline. +The graph of the second derivative of the quintic /B/-spline is therefore more visually appealing, though whether it is in fact more accurate depends on the smoothness of your data. + +And example usage is as follows: + + #include + using boost::math::interpolators::cardinal_quintic_b_spline; + std::vector v(512); + // fill v with data . . . + double t0 = 0; // initial time + double h = 0.125; // spacing + std::pair left_endpoint_derivatives{first_derivative_at_t0, second_derivative_at_t0}; + std::pair right_endpoint_derivatives{first_derivative_at_tf, second_derivative_at_tf}; + auto qs = cardinal_quintic_b_spline(v, t0, h, left_endpoint_derivatives, right_endpoint_derivatives); + + // Evaluate the interpolant at a point: + double y = qs(0.1); + // Evaluate the derivative of the interpolant: + double yp = qs.prime(0.1); + // Evaluate the second derivative of the interpolant: + double ypp = qs.double_prime(0.1); + +This routine will estimate the endpoint derivatives if they are not provided. +/Try to avoid this if possible./ +The endpoint derivatives must be evaluated by finite differences and this is not robust again perturbations in the data. +So if you have some way of knowing the endpoint derivatives, make sure to provide them. + +[heading References] + +Cox, Maurice G. ['Numerical methods for the interpolation and approximation of data by spline functions.] Diss. City, University of London, 1975. + +[endsect] [/section:cardinal_quintic_b] diff --git a/doc/math.qbk b/doc/math.qbk index 06c450071..879dad0d0 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -18,7 +18,7 @@ [block '''''']] [/insert Equation as an image, previous generated with an external tool like Latex.] -[template equation[name] +[template equation[name] [: ''' @@ -35,7 +35,7 @@ [/insert Graph as an SVG or PNG image, previous generated with an external tool like SVG_plot.] -[template graph[name] +[template graph[name] [: ''' @@ -50,7 +50,7 @@ [/insert Indented one-line expression italic and serif font probably using Unicode symbols for Greek and symbols.] [/Example: [expression [sub 1]F[sub 0](a, z) = (1-z)[super -a]]] -[template expression[equation] +[template expression[equation] [: [role serif_italic [equation]] ] @@ -690,6 +690,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [mathpart interpolation Interpolation] [include interpolators/cubic_b_spline.qbk] [include interpolators/cardinal_quadratic_b_spline.qbk] +[include interpolators/cardinal_quintic_b_spline.qbk] [include interpolators/whittaker_shannon.qbk] [include interpolators/barycentric_rational_interpolation.qbk] [include interpolators/vector_barycentric_rational.qbk] From 4f9d284e8333a4b504ca0babee3bd308935c518f Mon Sep 17 00:00:00 2001 From: NAThompson Date: Thu, 15 Aug 2019 11:39:18 -0400 Subject: [PATCH 3/7] Cardinal Quintic B-spline interpolator: Up and running. [CI SKIP] --- .../cardinal_quintic_b_spline_detail.hpp | 153 +++++++++++++++--- test/cardinal_quintic_b_spline_test.cpp | 147 +++++++++++++++++ 2 files changed, 280 insertions(+), 20 deletions(-) create mode 100644 test/cardinal_quintic_b_spline_test.cpp diff --git a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp index 9c172a66b..4c50f28a3 100644 --- a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp +++ b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp @@ -8,6 +8,7 @@ #define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP #include #include +#include namespace boost{ namespace math{ namespace interpolators{ namespace detail{ @@ -20,11 +21,11 @@ public: // y[0] = y(a), y[n -1] = y(b), step_size = (b - a)/(n -1). cardinal_quintic_b_spline_detail(const Real* const y, - size_t n, - Real t0 /* initial time, left endpoint */, - Real h /*spacing, stepsize*/, - std::pair left_endpoint_derivatives, - std::pair right_endpoint_derivatives) + size_t n, + Real t0 /* initial time, left endpoint */, + Real h /*spacing, stepsize*/, + std::pair left_endpoint_derivatives, + std::pair right_endpoint_derivatives) { if (h <= 0) { throw std::logic_error("Spacing must be > 0."); @@ -32,43 +33,155 @@ public: m_inv_h = 1/h; m_t0 = t0; - if (n < 3) { + if (n < 5) { throw std::logic_error("The interpolator requires at least 3 points."); } - m_alpha.resize(n + 4); + if(std::isnan(left_endpoint_derivatives.first) || std::isnan(left_endpoint_derivatives.second) || + std::isnan(right_endpoint_derivatives.first) || std::isnan(right_endpoint_derivatives.second)) { + throw std::logic_error("Derivative estimation is not yet implemented!"); + } + + // This is really challenging my mental limits on by-hand row reduction. + // 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. + std::vector rhs(n+4); - rhs[0] = 6*h*h*left_endpoint_derivatives.second; - rhs[1] = 24*h*left_endpoint_derivatives.first; + rhs[0] = 20*y[0] - 12*h*left_endpoint_derivatives.first + 2*h*h*left_endpoint_derivatives.second; + rhs[1] = 60*y[0] - 12*h*left_endpoint_derivatives.first; for (size_t i = 2; i < n + 2; ++i) { rhs[i] = 120*y[i-2]; } - rhs[n+2] = 24*h*right_endpoint_derivatives.first; - rhs[n+3] = 6*h*h*right_endpoint_derivatives.second; + rhs[n+2] = 60*y[n-1] + 12*h*right_endpoint_derivatives.first; + rhs[n+3] = 20*y[n-1] + 12*h*right_endpoint_derivatives.first + 2*h*h*right_endpoint_derivatives.second; - std::vector diagonal(n+4); + std::vector diagonal(n+4, 66); + diagonal[0] = 1; + diagonal[1] = 18; + diagonal[n+2] = 18; + diagonal[n+3] = 1; + std::vector first_superdiagonal(n+4, 26); + first_superdiagonal[0] = 10; + first_superdiagonal[1] = 33; + first_superdiagonal[n+2] = 1; + // There is one less superdiagonal than diagonal; make sure that if we read it, it shows up as a bug: + first_superdiagonal[n+3] = std::numeric_limits::quiet_NaN(); + + std::vector second_superdiagonal(n+4, 1); + second_superdiagonal[0] = 9; + second_superdiagonal[1] = 8; + second_superdiagonal[n+2] = std::numeric_limits::quiet_NaN(); + second_superdiagonal[n+3] = std::numeric_limits::quiet_NaN(); + + std::vector first_subdiagonal(n+4, 26); + first_subdiagonal[0] = std::numeric_limits::quiet_NaN(); + first_subdiagonal[1] = 1; + first_subdiagonal[n+2] = 33; + first_subdiagonal[n+3] = 10; + + std::vector second_subdiagonal(n+4, 1); + second_subdiagonal[0] = std::numeric_limits::quiet_NaN(); + second_subdiagonal[1] = std::numeric_limits::quiet_NaN(); + second_subdiagonal[n+2] = 8; + second_subdiagonal[n+3] = 9; + + + for (size_t i = 0; i < n+2; ++i) { + Real di = diagonal[i]; + diagonal[i] = 1; + first_superdiagonal[i] /= di; + second_superdiagonal[i] /= di; + rhs[i] /= di; + + // Eliminate first subdiagonal: + Real nfsub = -first_subdiagonal[i+1]; + // Superfluous: + first_subdiagonal[i+1] /= nfsub; + // Not superfluous: + diagonal[i+1] /= nfsub; + first_superdiagonal[i+1] /= nfsub; + second_superdiagonal[i+1] /= nfsub; + rhs[i+1] /= nfsub; + + diagonal[i+1] += first_superdiagonal[i]; + first_superdiagonal[i+1] += second_superdiagonal[i]; + rhs[i+1] += rhs[i]; + // Superfluous, but clarifying: + first_subdiagonal[i+1] = 0; + + // Eliminate second subdiagonal: + Real nssub = -second_subdiagonal[i+2]; + first_subdiagonal[i+2] /= nssub; + diagonal[i+2] /= nssub; + first_superdiagonal[i+2] /= nssub; + second_superdiagonal[i+2] /= nssub; + rhs[i+2] /= nssub; + + first_subdiagonal[i+2] += first_superdiagonal[i]; + diagonal[i+2] += second_superdiagonal[i]; + rhs[i+2] += rhs[i]; + // Superfluous, but clarifying: + second_subdiagonal[i+2] = 0; + } + + // Eliminate last subdiagonal: + Real dnp2 = diagonal[n+2]; + diagonal[n+2] = 1; + first_superdiagonal[n+2] /= dnp2; + rhs[n+2] /= dnp2; + Real nfsubnp3 = -first_subdiagonal[n+3]; + diagonal[n+3] /= nfsubnp3; + rhs[n+3] /= nfsubnp3; + + diagonal[n+3] += first_superdiagonal[n+2]; + rhs[n+3] += rhs[n+2]; + + m_alpha.resize(n + 4, std::numeric_limits::quiet_NaN()); + + m_alpha[n+3] = rhs[n+3]/diagonal[n+3]; + m_alpha[n+2] = rhs[n+2] - first_superdiagonal[n+2]*m_alpha[n+3]; + for (int64_t i = int64_t(n+1); i >= 0; --i) { + m_alpha[i] = rhs[i] - first_superdiagonal[i]*m_alpha[i+1] - second_superdiagonal[i]*m_alpha[i+2]; + } + + + /*std::cout << "alpha = {"; + for (auto & a : m_alpha) { + std::cout << a << ", "; + } + std::cout << "}\n";*/ } Real operator()(Real t) const { - if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { - const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + using boost::math::cardinal_b_spline; + // tf = t0 + (n-1)*h + // alpha.size() = n+4 + if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work."; throw std::domain_error(err_msg); } - return std::numeric_limits::quiet_NaN(); + Real x = (t-m_t0)*m_inv_h; + // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5 + int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1))); + int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) ); + Real s = 0; + for (int64_t j = j_min; j <= j_max; ++j) { + s += m_alpha[j]*cardinal_b_spline<5, Real>(x - j + 2); + } + return s; } Real prime(Real t) const { - if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { - const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work."; throw std::domain_error(err_msg); } return std::numeric_limits::quiet_NaN(); } Real double_prime(Real t) const { - if (t < m_t0 || t > m_t0 + (m_alpha.size()-2)/m_inv_h) { - const char* err_msg = "Tried to evaluate the cardinal quadratic b-spline outside the domain of of interpolation; extrapolation does not work."; + if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) { + const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work."; throw std::domain_error(err_msg); } return std::numeric_limits::quiet_NaN(); @@ -76,7 +189,7 @@ public: Real t_max() const { - return m_t0 + (m_alpha.size()-3)/m_inv_h; + return m_t0 + (m_alpha.size()-5)/m_inv_h; } private: diff --git a/test/cardinal_quintic_b_spline_test.cpp b/test/cardinal_quintic_b_spline_test.cpp new file mode 100644 index 000000000..2de71c87e --- /dev/null +++ b/test/cardinal_quintic_b_spline_test.cpp @@ -0,0 +1,147 @@ +/* + * Copyright Nick Thompson, 2019 + * 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) + */ + +#include "math_unit_test.hpp" +#include +#include +#include +using boost::math::interpolators::cardinal_quintic_b_spline; + +template +void test_constant() +{ + Real c = 7.2; + Real t0 = 0; + Real h = Real(1)/Real(16); + size_t n = 513; + std::vector v(n, c); + std::pair left_endpoint_derivatives{0, 0}; + std::pair right_endpoint_derivatives{0, 0}; + auto qbs = cardinal_quintic_b_spline(v.data(), v.size(), t0, h, left_endpoint_derivatives, right_endpoint_derivatives); + + size_t i = 0; + while (i < n) { + Real t = t0 + i*h; + CHECK_ULP_CLOSE(c, qbs(t), 3); + //CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 100*std::numeric_limits::epsilon()); + ++i; + } + + i = 0; + while (i < n - 1) { + Real t = t0 + i*h + h/2; + CHECK_ULP_CLOSE(c, qbs(t), 4); + //CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 300*std::numeric_limits::epsilon()); + t = t0 + i*h + h/4; + CHECK_ULP_CLOSE(c, qbs(t), 4); + //CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 150*std::numeric_limits::epsilon()); + ++i; + } +} + + +template +void test_linear() +{ + Real m = 8.3; + Real b = 7.2; + Real t0 = 0; + Real h = Real(1)/Real(16); + size_t n = 512; + std::vector y(n); + for (size_t i = 0; i < n; ++i) { + Real t = i*h; + y[i] = m*t + b; + } + std::pair left_endpoint_derivatives{m, 0}; + std::pair right_endpoint_derivatives{m, 0}; + auto qbs = cardinal_quintic_b_spline(y.data(), y.size(), t0, h, left_endpoint_derivatives, right_endpoint_derivatives); + + size_t i = 0; + while (i < n) { + Real t = t0 + i*h; + if (!CHECK_ULP_CLOSE(m*t+b, qbs(t), 3)) { + std::cerr << " Problem at t = " << t << "\n"; + } + //CHECK_ULP_CLOSE(m, qbs.prime(t), 820); + ++i; + } + + i = 0; + while (i < n - 1) { + Real t = t0 + i*h + h/2; + if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) { + std::cerr << " Problem at t = " << t << "\n"; + } + //CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits::epsilon()); + t = t0 + i*h + h/4; + if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) { + std::cerr << " Problem at t = " << t << "\n"; + } + //CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits::epsilon()); + ++i; + } +} + + +template +void test_quadratic() +{ + Real a = Real(1)/Real(16); + Real b = -3.5; + Real c = -9; + Real t0 = 0; + Real h = Real(1)/Real(16); + size_t n = 513; + std::vector y(n); + for (size_t i = 0; i < n; ++i) { + Real t = i*h; + y[i] = a*t*t + b*t + c; + } + Real t_max = t0 + (n-1)*h; + std::pair left_endpoint_derivatives{b, 2*a}; + std::pair right_endpoint_derivatives{2*a*t_max + b, 2*a}; + + auto qbs = cardinal_quintic_b_spline(y, t0, h, left_endpoint_derivatives, right_endpoint_derivatives); + + size_t i = 0; + while (i < n) { + Real t = t0 + i*h; + CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 3); + ++i; + } + + i = 0; + while (i < n -1) { + Real t = t0 + i*h + h/2; + if(!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 5)) { + std::cerr << " Problem at abscissa t = " << t << "\n"; + } + + t = t0 + i*h + h/4; + if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 5)) { + std::cerr << " Problem abscissa t = " << t << "\n"; + } + ++i; + } +} + +int main() +{ + test_constant(); + test_constant(); + test_constant(); + + test_linear(); + test_linear(); + test_linear(); + + test_quadratic(); + //test_quadratic(); + + return boost::math::test::report_errors(); +} From 1483516247169e86b9ae133f789b3066465f1871 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Thu, 15 Aug 2019 14:13:25 -0400 Subject: [PATCH 4/7] Cardinal quintic B-splines: First and second derivatives [CI SKIP] --- .../cardinal_quintic_b_spline_detail.hpp | 35 +++++++++++++++-- test/cardinal_quintic_b_spline_test.cpp | 38 ++++++++++++++----- 2 files changed, 59 insertions(+), 14 deletions(-) diff --git a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp index 4c50f28a3..586c5eefb 100644 --- a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp +++ b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp @@ -6,6 +6,7 @@ #ifndef BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP #define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP +#include #include #include #include @@ -37,8 +38,9 @@ public: throw std::logic_error("The interpolator requires at least 3 points."); } - if(std::isnan(left_endpoint_derivatives.first) || std::isnan(left_endpoint_derivatives.second) || - std::isnan(right_endpoint_derivatives.first) || std::isnan(right_endpoint_derivatives.second)) { + using std::isnan; + if(isnan(left_endpoint_derivatives.first) || isnan(left_endpoint_derivatives.second) || + isnan(right_endpoint_derivatives.first) || isnan(right_endpoint_derivatives.second)) { throw std::logic_error("Derivative estimation is not yet implemented!"); } @@ -153,6 +155,8 @@ public: } Real operator()(Real t) const { + using std::ceil; + using std::floor; using boost::math::cardinal_b_spline; // tf = t0 + (n-1)*h // alpha.size() = n+4 @@ -172,19 +176,42 @@ public: } Real prime(Real t) const { + using std::ceil; + using std::floor; + using boost::math::cardinal_b_spline_prime; if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) { const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work."; throw std::domain_error(err_msg); } - return std::numeric_limits::quiet_NaN(); + Real x = (t-m_t0)*m_inv_h; + // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5 + int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1))); + int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) ); + Real s = 0; + for (int64_t j = j_min; j <= j_max; ++j) { + s += m_alpha[j]*cardinal_b_spline_prime<5, Real>(x - j + 2); + } + return s*m_inv_h; + } Real double_prime(Real t) const { + using std::ceil; + using std::floor; + using boost::math::cardinal_b_spline_double_prime; if (t < m_t0 || t > m_t0 + (m_alpha.size()-5)/m_inv_h) { const char* err_msg = "Tried to evaluate the cardinal quintic b-spline outside the domain of of interpolation; extrapolation does not work."; throw std::domain_error(err_msg); } - return std::numeric_limits::quiet_NaN(); + Real x = (t-m_t0)*m_inv_h; + // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5 + int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1))); + int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) ); + Real s = 0; + for (int64_t j = j_min; j <= j_max; ++j) { + s += m_alpha[j]*cardinal_b_spline_double_prime<5, Real>(x - j + 2); + } + return s*m_inv_h*m_inv_h; } diff --git a/test/cardinal_quintic_b_spline_test.cpp b/test/cardinal_quintic_b_spline_test.cpp index 2de71c87e..10b6b3a1e 100644 --- a/test/cardinal_quintic_b_spline_test.cpp +++ b/test/cardinal_quintic_b_spline_test.cpp @@ -9,12 +9,16 @@ #include #include #include +#ifdef BOOST_HAS_FLOAT128 +#include +using boost::multiprecision::float128; +#endif using boost::math::interpolators::cardinal_quintic_b_spline; template void test_constant() { - Real c = 7.2; + Real c = 7.5; Real t0 = 0; Real h = Real(1)/Real(16); size_t n = 513; @@ -27,18 +31,21 @@ void test_constant() while (i < n) { Real t = t0 + i*h; CHECK_ULP_CLOSE(c, qbs(t), 3); - //CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 100*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 400*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 60000*std::numeric_limits::epsilon()); ++i; } i = 0; while (i < n - 1) { Real t = t0 + i*h + h/2; - CHECK_ULP_CLOSE(c, qbs(t), 4); - //CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 300*std::numeric_limits::epsilon()); + CHECK_ULP_CLOSE(c, qbs(t), 5); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 600*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 30000*std::numeric_limits::epsilon()); t = t0 + i*h + h/4; CHECK_ULP_CLOSE(c, qbs(t), 4); - //CHECK_MOLLIFIED_CLOSE(0, qbs.prime(t), 150*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 600*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 10000*std::numeric_limits::epsilon()); ++i; } } @@ -47,6 +54,7 @@ void test_constant() template void test_linear() { + using std::abs; Real m = 8.3; Real b = 7.2; Real t0 = 0; @@ -67,7 +75,12 @@ void test_linear() if (!CHECK_ULP_CLOSE(m*t+b, qbs(t), 3)) { std::cerr << " Problem at t = " << t << "\n"; } - //CHECK_ULP_CLOSE(m, qbs.prime(t), 820); + if(!CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 100*abs(m*t+b)*std::numeric_limits::epsilon())) { + std::cerr << " Problem at t = " << t << "\n"; + } + if(!CHECK_MOLLIFIED_CLOSE(0, qbs.double_prime(t), 10000*abs(m*t+b)*std::numeric_limits::epsilon())) { + std::cerr << " Problem at t = " << t << "\n"; + } ++i; } @@ -77,12 +90,12 @@ void test_linear() if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) { std::cerr << " Problem at t = " << t << "\n"; } - //CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits::epsilon()); t = t0 + i*h + h/4; if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) { std::cerr << " Problem at t = " << t << "\n"; } - //CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 3000*std::numeric_limits::epsilon()); ++i; } } @@ -132,7 +145,6 @@ void test_quadratic() int main() { - test_constant(); test_constant(); test_constant(); @@ -140,8 +152,14 @@ int main() test_linear(); test_linear(); + test_quadratic(); test_quadratic(); - //test_quadratic(); + + + #ifdef BOOST_HAS_FLOAT128 + test_constant(); + test_linear(); + #endif return boost::math::test::report_errors(); } From 32100ebea498147c3b4296a39aba538cd90bba1f Mon Sep 17 00:00:00 2001 From: NAThompson Date: Sun, 18 Aug 2019 11:09:05 -0400 Subject: [PATCH 5/7] Cardinal quintic B-spline interpolation: Add tests to jamfile and kick off build. --- test/Jamfile.v2 | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 53905d8d1..937f3abcf 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -941,6 +941,7 @@ test-suite misc : [ run cardinal_b_spline_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run whittaker_shannon_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] ] [ run cardinal_quadratic_b_spline_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] ] + [ run cardinal_quintic_b_spline_test.cpp : : : [ requires cxx11_auto_declarations cxx11_constexpr cxx11_smart_ptr cxx11_defaulted_functions ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : -lquadmath ] ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=1 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_1 ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=2 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_2 ] [ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : TEST=3 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_3 ] From 7364de0090f75871504d74ffa2ee322ea74eec4f Mon Sep 17 00:00:00 2001 From: NAThompson Date: Sun, 18 Aug 2019 12:08:34 -0400 Subject: [PATCH 6/7] Cardinal Quintic B-spline: Derivative estimation. [CI SKIP] --- .../cardinal_quintic_b_spline_detail.hpp | 34 +++-- test/cardinal_quintic_b_spline_test.cpp | 131 ++++++++++++++++++ 2 files changed, 153 insertions(+), 12 deletions(-) diff --git a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp index 586c5eefb..9663ff061 100644 --- a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp +++ b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp @@ -34,14 +34,28 @@ public: m_inv_h = 1/h; m_t0 = t0; - if (n < 5) { - throw std::logic_error("The interpolator requires at least 3 points."); + if (n < 8) { + throw std::logic_error("The quntic B-spline interpolator requires at least 8 points."); } using std::isnan; - if(isnan(left_endpoint_derivatives.first) || isnan(left_endpoint_derivatives.second) || - isnan(right_endpoint_derivatives.first) || isnan(right_endpoint_derivatives.second)) { - throw std::logic_error("Derivative estimation is not yet implemented!"); + // This interpolator has error of order h^6, so the derivatives should be estimated with the same error. + // See: https://en.wikipedia.org/wiki/Finite_difference_coefficient + if (isnan(left_endpoint_derivatives.first)) { + 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; + left_endpoint_derivatives.first = tmp/h; + } + if (isnan(right_endpoint_derivatives.first)) { + 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; + right_endpoint_derivatives.first = tmp/h; + } + if(isnan(left_endpoint_derivatives.second)) { + 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; + left_endpoint_derivatives.second = tmp/(h*h); + } + if (isnan(right_endpoint_derivatives.second)) { + 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; + right_endpoint_derivatives.second = tmp/(h*h); } // This is really challenging my mental limits on by-hand row reduction. @@ -146,12 +160,6 @@ public: m_alpha[i] = rhs[i] - first_superdiagonal[i]*m_alpha[i+1] - second_superdiagonal[i]*m_alpha[i+2]; } - - /*std::cout << "alpha = {"; - for (auto & a : m_alpha) { - std::cout << a << ", "; - } - std::cout << "}\n";*/ } Real operator()(Real t) const { @@ -165,11 +173,13 @@ public: throw std::domain_error(err_msg); } Real x = (t-m_t0)*m_inv_h; - // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5 + // Support of B_5 is [-3, 3]. So -3 < x - j + 2 < 3, so x-1 < j < x+5. + // TODO: Zero pad m_alpha so that only the domain check is necessary. int64_t j_min = std::max(int64_t(0), int64_t(ceil(x-1))); int64_t j_max = std::min(int64_t(m_alpha.size() - 1), int64_t(floor(x+5)) ); Real s = 0; for (int64_t j = j_min; j <= j_max; ++j) { + // TODO: Use Cox 1972 to generate all integer translates of B5 simultaneously. s += m_alpha[j]*cardinal_b_spline<5, Real>(x - j + 2); } return s; diff --git a/test/cardinal_quintic_b_spline_test.cpp b/test/cardinal_quintic_b_spline_test.cpp index 10b6b3a1e..31b9fd910 100644 --- a/test/cardinal_quintic_b_spline_test.cpp +++ b/test/cardinal_quintic_b_spline_test.cpp @@ -50,6 +50,39 @@ void test_constant() } } +template +void test_constant_estimate_derivatives() +{ + Real c = 7.5; + Real t0 = 0; + Real h = Real(1)/Real(16); + size_t n = 513; + std::vector v(n, c); + auto qbs = cardinal_quintic_b_spline(v.data(), v.size(), t0, h); + + size_t i = 0; + while (i < n) { + Real t = t0 + i*h; + CHECK_ULP_CLOSE(c, qbs(t), 3); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 1200*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 200000*std::numeric_limits::epsilon()); + ++i; + } + + i = 0; + while (i < n - 1) { + Real t = t0 + i*h + h/2; + CHECK_ULP_CLOSE(c, qbs(t), 8); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 1200*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 80000*std::numeric_limits::epsilon()); + t = t0 + i*h + h/4; + CHECK_ULP_CLOSE(c, qbs(t), 5); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 1200*std::numeric_limits::epsilon()); + CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 38000*std::numeric_limits::epsilon()); + ++i; + } +} + template void test_linear() @@ -100,6 +133,54 @@ void test_linear() } } +template +void test_linear_estimate_derivatives() +{ + using std::abs; + Real m = 8.3; + Real b = 7.2; + Real t0 = 0; + Real h = Real(1)/Real(16); + size_t n = 512; + std::vector y(n); + for (size_t i = 0; i < n; ++i) { + Real t = i*h; + y[i] = m*t + b; + } + + auto qbs = cardinal_quintic_b_spline(y.data(), y.size(), t0, h); + + size_t i = 0; + while (i < n) { + Real t = t0 + i*h; + if (!CHECK_ULP_CLOSE(m*t+b, qbs(t), 3)) { + std::cerr << " Problem at t = " << t << "\n"; + } + if(!CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 100*abs(m*t+b)*std::numeric_limits::epsilon())) { + std::cerr << " Problem at t = " << t << "\n"; + } + if(!CHECK_MOLLIFIED_CLOSE(0, qbs.double_prime(t), 20000*abs(m*t+b)*std::numeric_limits::epsilon())) { + std::cerr << " Problem at t = " << t << "\n"; + } + ++i; + } + + i = 0; + while (i < n - 1) { + Real t = t0 + i*h + h/2; + if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 5)) { + std::cerr << " Problem at t = " << t << "\n"; + } + CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 1500*std::numeric_limits::epsilon()); + t = t0 + i*h + h/4; + if(!CHECK_ULP_CLOSE(m*t+b, qbs(t), 4)) { + std::cerr << " Problem at t = " << t << "\n"; + } + CHECK_MOLLIFIED_CLOSE(m, qbs.prime(t), 3000*std::numeric_limits::epsilon()); + ++i; + } +} + template void test_quadratic() @@ -143,22 +224,72 @@ void test_quadratic() } } + +template +void test_quadratic_estimate_derivatives() +{ + Real a = Real(1)/Real(16); + Real b = -3.5; + Real c = -9; + Real t0 = 0; + Real h = Real(1)/Real(16); + size_t n = 513; + std::vector y(n); + for (size_t i = 0; i < n; ++i) { + Real t = i*h; + y[i] = a*t*t + b*t + c; + } + auto qbs = cardinal_quintic_b_spline(y, t0, h); + + size_t i = 0; + while (i < n) { + Real t = t0 + i*h; + CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 3); + ++i; + } + + i = 0; + while (i < n -1) { + Real t = t0 + i*h + h/2; + if(!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 10)) { + std::cerr << " Problem at abscissa t = " << t << "\n"; + } + + t = t0 + i*h + h/4; + if (!CHECK_ULP_CLOSE(a*t*t + b*t + c, qbs(t), 6)) { + std::cerr << " Problem abscissa t = " << t << "\n"; + } + ++i; + } +} + + int main() { test_constant(); test_constant(); + test_constant_estimate_derivatives(); + test_constant_estimate_derivatives(); + test_linear(); test_linear(); test_linear(); + test_linear_estimate_derivatives(); + test_linear_estimate_derivatives(); + test_quadratic(); test_quadratic(); + test_quadratic_estimate_derivatives(); + test_quadratic_estimate_derivatives(); + #ifdef BOOST_HAS_FLOAT128 test_constant(); test_linear(); + test_linear_estimate_derivatives(); #endif return boost::math::test::report_errors(); From e8b9c128864ac8bc367a4526cc157379ace22d01 Mon Sep 17 00:00:00 2001 From: NAThompson Date: Thu, 29 Aug 2019 09:02:20 -0400 Subject: [PATCH 7/7] Make sure not to allow instantiation on integer types. [CI SKIP] --- .../interpolators/detail/cardinal_quintic_b_spline_detail.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp index 9663ff061..09ce63d82 100644 --- a/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp +++ b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp @@ -28,6 +28,7 @@ public: std::pair left_endpoint_derivatives, std::pair right_endpoint_derivatives) { + static_assert(!std::is_integral::value, "The quintic B-spline interpolator only works with floating point types."); if (h <= 0) { throw std::logic_error("Spacing must be > 0."); }