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 7e59adf48..ae87ddaf3 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -692,6 +692,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] 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..09ce63d82 --- /dev/null +++ b/include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp @@ -0,0 +1,240 @@ +// 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 +#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) + { + 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."); + } + m_inv_h = 1/h; + m_t0 = t0; + + if (n < 8) { + throw std::logic_error("The quntic B-spline interpolator requires at least 8 points."); + } + + using std::isnan; + // 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. + // 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] = 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] = 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, 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]; + } + + } + + 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 + 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); + } + 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. + // 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; + } + + 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); + } + 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); + } + 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; + } + + + Real t_max() const { + return m_t0 + (m_alpha.size()-5)/m_inv_h; + } + +private: + std::vector m_alpha; + Real m_inv_h; + Real m_t0; +}; + +}}}} +#endif diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 39eb99678..e5f3122ab 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -943,6 +943,7 @@ test-suite misc : [ run gegenbauer_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 ] diff --git a/test/cardinal_quintic_b_spline_test.cpp b/test/cardinal_quintic_b_spline_test.cpp new file mode 100644 index 000000000..31b9fd910 --- /dev/null +++ b/test/cardinal_quintic_b_spline_test.cpp @@ -0,0 +1,296 @@ +/* + * 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 +#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.5; + 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(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), 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(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; + } +} + +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() +{ + 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; + } + 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"; + } + 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; + } + + 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), 3000*std::numeric_limits::epsilon()); + ++i; + } +} + +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() +{ + 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; + } +} + + +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(); +}