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: Derivative estimation. [CI SKIP]

Browse Source
This commit is contained in:
NAThompson
2019-08-18 12:08:34 -04:00
parent 648d2b0022
commit 7364de0090
2 changed files with 153 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -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;

131
test/cardinal_quintic_b_spline_test.cpp
View File

@@ -50,6 +50,39 @@ void test_constant()
}
}
template<class Real>
void test_constant_estimate_derivatives()
{
Real c = 7.5;
Real t0 = 0;
Real h = Real(1)/Real(16);
size_t n = 513;
std::vector<Real> v(n, c);
auto qbs = cardinal_quintic_b_spline<Real>(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<Real>::epsilon());
CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 200000*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), 8);
CHECK_MOLLIFIED_CLOSE(Real(0), qbs.prime(t), 1200*std::numeric_limits<Real>::epsilon());
CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 80000*std::numeric_limits<Real>::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<Real>::epsilon());
CHECK_MOLLIFIED_CLOSE(Real(0), qbs.double_prime(t), 38000*std::numeric_limits<Real>::epsilon());
++i;
}
}
template<class Real>
void test_linear()
@@ -100,6 +133,54 @@ void test_linear()
}
}
template<class Real>
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<Real> 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<Real>(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<Real>::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<Real>::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<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), 3000*std::numeric_limits<Real>::epsilon());
++i;
}
}
template<class Real>
void test_quadratic()
@@ -143,22 +224,72 @@ void test_quadratic()
}
}
template<class Real>
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<Real> 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<Real>(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<double>();
test_constant<long double>();
test_constant_estimate_derivatives<double>();
test_constant_estimate_derivatives<long double>();
test_linear<float>();
test_linear<double>();
test_linear<long double>();
test_linear_estimate_derivatives<double>();
test_linear_estimate_derivatives<long double>();
test_quadratic<double>();
test_quadratic<long double>();
test_quadratic_estimate_derivatives<double>();
test_quadratic_estimate_derivatives<long double>();
#ifdef BOOST_HAS_FLOAT128
test_constant<float128>();
test_linear<float128>();
test_linear_estimate_derivatives<float128>();
#endif
return boost::math::test::report_errors();