boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Cardinal Quintic B-spline interpolator: Up and running. [CI SKIP]

Browse Source
This commit is contained in:
NAThompson
2019-08-15 11:39:18 -04:00
parent 1955699777
commit 4f9d284e83
2 changed files with 280 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

153
include/boost/math/interpolators/detail/cardinal_quintic_b_spline_detail.hpp
View File

@@ -8,6 +8,7 @@
#define BOOST_MATH_INTERPOLATORS_CARDINAL_QUINTIC_B_SPLINE_DETAIL_HPP
#include <vector>
#include <utility>
#include <boost/math/special_functions/cardinal_b_spline.hpp>
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<Real, Real> left_endpoint_derivatives,
std::pair<Real, Real> right_endpoint_derivatives)
size_t n,
Real t0 /* initial time, left endpoint */,
Real h /*spacing, stepsize*/,
std::pair<Real, Real> left_endpoint_derivatives,
std::pair<Real, Real> 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<Real> 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<Real> diagonal(n+4);
std::vector<Real> diagonal(n+4, 66);
diagonal[0] = 1;
diagonal[1] = 18;
diagonal[n+2] = 18;
diagonal[n+3] = 1;
std::vector<Real> 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<Real>::quiet_NaN();
std::vector<Real> second_superdiagonal(n+4, 1);
second_superdiagonal[0] = 9;
second_superdiagonal[1] = 8;
second_superdiagonal[n+2] = std::numeric_limits<Real>::quiet_NaN();
second_superdiagonal[n+3] = std::numeric_limits<Real>::quiet_NaN();
std::vector<Real> first_subdiagonal(n+4, 26);
first_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
first_subdiagonal[1] = 1;
first_subdiagonal[n+2] = 33;
first_subdiagonal[n+3] = 10;
std::vector<Real> second_subdiagonal(n+4, 1);
second_subdiagonal[0] = std::numeric_limits<Real>::quiet_NaN();
second_subdiagonal[1] = std::numeric_limits<Real>::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<Real>::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<Real>::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<Real>::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<Real>::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:

147
test/cardinal_quintic_b_spline_test.cpp Normal file
View File

@@ -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 <numeric>
#include <utility>
#include <boost/math/interpolators/cardinal_quintic_b_spline.hpp>
using boost::math::interpolators::cardinal_quintic_b_spline;
template<class Real>
void test_constant()
{
Real c = 7.2;
Real t0 = 0;
Real h = Real(1)/Real(16);
size_t n = 513;
std::vector<Real> v(n, c);
std::pair<Real, Real> left_endpoint_derivatives{0, 0};
std::pair<Real, Real> right_endpoint_derivatives{0, 0};
auto qbs = cardinal_quintic_b_spline<Real>(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<Real>::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<Real>::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<Real>::epsilon());
++i;
}
}
template<class Real>
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<Real> y(n);
for (size_t i = 0; i < n; ++i) {
Real t = i*h;
y[i] = m*t + b;
}
std::pair<Real, Real> left_endpoint_derivatives{m, 0};
std::pair<Real, Real> right_endpoint_derivatives{m, 0};
auto qbs = cardinal_quintic_b_spline<Real>(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<Real>::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<Real>::epsilon());
++i;
}
}
template<class Real>
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<Real> 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<Real, Real> left_endpoint_derivatives{b, 2*a};
std::pair<Real, Real> right_endpoint_derivatives{2*a*t_max + b, 2*a};
auto qbs = cardinal_quintic_b_spline<Real>(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<float>();
test_constant<double>();
test_constant<long double>();
test_linear<float>();
test_linear<double>();
test_linear<long double>();
test_quadratic<long double>();
//test_quadratic<long double>();
return boost::math::test::report_errors();
}