From 7364de0090f75871504d74ffa2ee322ea74eec4f Mon Sep 17 00:00:00 2001 From: NAThompson Date: Sun, 18 Aug 2019 12:08:34 -0400 Subject: [PATCH] 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();