mirror of
https://github.com/boostorg/math.git
synced 2026-01-19 04:22:09 +00:00
O(1) septic Hermite splines [CI SKIP]
This commit is contained in:
NAThompson
@@ -160,5 +160,225 @@ private:
|
||||
RandomAccessContainer d2ydx2_;
|
||||
RandomAccessContainer d3ydx3_;
|
||||
};
|
||||
|
||||
template<class RandomAccessContainer>
|
||||
class cardinal_septic_hermite_detail {
|
||||
public:
|
||||
using Real = typename RandomAccessContainer::value_type;
|
||||
cardinal_septic_hermite_detail(RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3, Real x0, Real dx)
|
||||
: y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}, d3ydx3_{std::move(d3ydx3)}, x0_{x0}, dx_{dx}
|
||||
{
|
||||
if (y_.size() != dydx_.size()) {
|
||||
throw std::domain_error("Numbers of derivatives must = number of ordinates.");
|
||||
}
|
||||
if (y_.size() != d2ydx2_.size()) {
|
||||
throw std::domain_error("Number of second derivatives must equal number of ordinates.");
|
||||
}
|
||||
if (y_.size() != d3ydx3_.size()) {
|
||||
throw std::domain_error("Number of third derivatives must equal number of ordinates.");
|
||||
}
|
||||
|
||||
if (y_.size() < 2) {
|
||||
throw std::domain_error("At least 2 abscissas are required.");
|
||||
}
|
||||
}
|
||||
|
||||
inline Real operator()(Real x) const {
|
||||
Real xf = x0_ + (y_.size()-1)*dx_;
|
||||
if (x < x0_ || x > xf) {
|
||||
std::ostringstream oss;
|
||||
oss.precision(std::numeric_limits<Real>::digits10+3);
|
||||
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
|
||||
<< x0_ << ", " << xf << "]";
|
||||
throw std::domain_error(oss.str());
|
||||
}
|
||||
// We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work.
|
||||
// Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf.
|
||||
if (x == xf)
|
||||
{
|
||||
return y_.back();
|
||||
}
|
||||
return this->unchecked_evaluation(x);
|
||||
}
|
||||
|
||||
inline Real unchecked_evaluation(Real x) const {
|
||||
using std::floor;
|
||||
auto i = static_cast<decltype(y_.size())>(floor((x-x0_)/dx_));
|
||||
Real xi = x0_ + i*dx_;
|
||||
Real t = (x - xi)/dx_;
|
||||
|
||||
Real y0 = y_[i];
|
||||
Real y1 = y_[i+1];
|
||||
// Velocity:
|
||||
Real v0 = dydx_[i];
|
||||
Real v1 = dydx_[i+1];
|
||||
// Acceleration:
|
||||
Real a0 = d2ydx2_[i];
|
||||
Real a1 = d2ydx2_[i+1];
|
||||
// Jerk:
|
||||
Real j0 = d3ydx3_[i];
|
||||
Real j1 = d3ydx3_[i+1];
|
||||
|
||||
// See:
|
||||
// http://seisweb.usask.ca/classes/GEOL481/2017/Labs/interpolation_utilities_matlab/shermite.m
|
||||
Real t2 = t*t;
|
||||
Real t3 = t2*t;
|
||||
Real t4 = t3*t;
|
||||
Real t5 = t4*t;
|
||||
Real t6 = t5*t;
|
||||
Real t7 = t6*t;
|
||||
Real dx2 = dx_*dx_;
|
||||
Real dx3 = dx2*dx_;
|
||||
|
||||
Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
|
||||
Real z4 = -s;
|
||||
Real z0 = s + 1;
|
||||
Real z1 = 10*t7 - 36*t6 + 45*t5 - 20*t4 + t;
|
||||
Real z2 = 2*t7 - 15*t6/2 + 10*t5 - 5*t4 + t2/2;
|
||||
Real z3 = t7/6 - 2*t6/3 + t5 - 2*t4/3 + t3/6;
|
||||
|
||||
Real z5 = 10*t7 - 34*t6 + 39*t5 - 15*t4;
|
||||
Real z6 = -2*t7 + 13*t6/2 - 7*t5 + 5*t4/2;
|
||||
Real z7 = t7/6 - t6/2 + t5/2 - t4/6;
|
||||
|
||||
return z0*y0 + z1*dx_*v0 + z2*dx2*a0 + z3*dx3*j0 + z4*y1 + z5*dx_*v1 + z6*dx2*a1 + z7*dx3*j1;
|
||||
}
|
||||
|
||||
inline Real prime(Real x) const {
|
||||
Real xf = x0_ + (y_.size()-1)*dx_;
|
||||
if (x < x0_ || x > xf)
|
||||
{
|
||||
std::ostringstream oss;
|
||||
oss.precision(std::numeric_limits<Real>::digits10+3);
|
||||
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
|
||||
<< x0_ << ", " << xf << "]";
|
||||
throw std::domain_error(oss.str());
|
||||
}
|
||||
if (x == xf)
|
||||
{
|
||||
return dydx_.back();
|
||||
}
|
||||
|
||||
return this->unchecked_prime(x);
|
||||
}
|
||||
|
||||
inline Real unchecked_prime(Real x) const {
|
||||
return std::numeric_limits<Real>::quiet_NaN();
|
||||
}
|
||||
|
||||
private:
|
||||
RandomAccessContainer y_;
|
||||
RandomAccessContainer dydx_;
|
||||
RandomAccessContainer d2ydx2_;
|
||||
RandomAccessContainer d3ydx3_;
|
||||
Real x0_;
|
||||
Real dx_;
|
||||
};
|
||||
|
||||
|
||||
template<class RandomAccessContainer>
|
||||
class cardinal_septic_hermite_detail_aos {
|
||||
public:
|
||||
using Point = typename RandomAccessContainer::value_type;
|
||||
using Real = typename Point::value_type;
|
||||
cardinal_septic_hermite_detail_aos(RandomAccessContainer && dat, Real x0, Real dx)
|
||||
: data_{std::move(dat)}, x0_{x0}, dx_{dx}
|
||||
{
|
||||
if (data_.size() < 2) {
|
||||
throw std::domain_error("At least 2 abscissas are required.");
|
||||
}
|
||||
if (data_[0].size() != 4) {
|
||||
throw std::domain_error("There must be 4 data items per struct.");
|
||||
}
|
||||
}
|
||||
|
||||
inline Real operator()(Real x) const {
|
||||
Real xf = x0_ + (data_.size()-1)*dx_;
|
||||
if (x < x0_ || x > xf) {
|
||||
std::ostringstream oss;
|
||||
oss.precision(std::numeric_limits<Real>::digits10+3);
|
||||
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
|
||||
<< x0_ << ", " << xf << "]";
|
||||
throw std::domain_error(oss.str());
|
||||
}
|
||||
// We need t := (x-x_k)/(x_{k+1}-x_k) \in [0,1) for this to work.
|
||||
// Sadly this neccessitates this loathesome check, otherwise we get t = 1 at x = xf.
|
||||
if (x == xf)
|
||||
{
|
||||
return data_.back()[0];
|
||||
}
|
||||
return this->unchecked_evaluation(x);
|
||||
}
|
||||
|
||||
inline Real unchecked_evaluation(Real x) const {
|
||||
using std::floor;
|
||||
auto i = static_cast<decltype(data_.size())>(floor((x-x0_)/dx_));
|
||||
Real xi = x0_ + i*dx_;
|
||||
Real t = (x - xi)/dx_;
|
||||
|
||||
Real y0 = data_[i][0];
|
||||
Real y1 = data_[i+1][0];
|
||||
// Velocity:
|
||||
Real v0 = data_[i][1];
|
||||
Real v1 = data_[i+1][1];
|
||||
// Acceleration:
|
||||
Real a0 = data_[i][2];
|
||||
Real a1 = data_[i+1][2];
|
||||
// Jerk:
|
||||
Real j0 = data_[i][3];
|
||||
Real j1 = data_[i+1][3];
|
||||
|
||||
Real t2 = t*t;
|
||||
Real t3 = t2*t;
|
||||
Real t4 = t3*t;
|
||||
Real t5 = t4*t;
|
||||
Real t6 = t5*t;
|
||||
Real t7 = t6*t;
|
||||
Real dx2 = dx_*dx_;
|
||||
Real dx3 = dx2*dx_;
|
||||
|
||||
Real s = t4*(-35 + t*(84 + t*(-70 + 20*t)));
|
||||
Real z4 = -s;
|
||||
Real z0 = s + 1;
|
||||
Real z1 = 10*t7 - 36*t6 + 45*t5 - 20*t4 + t;
|
||||
Real z2 = 2*t7 - 15*t6/2 + 10*t5 - 5*t4 + t2/2;
|
||||
Real z3 = t7/6 - 2*t6/3 + t5 - 2*t4/3 + t3/6;
|
||||
|
||||
Real z5 = 10*t7 - 34*t6 + 39*t5 - 15*t4;
|
||||
Real z6 = -2*t7 + 13*t6/2 - 7*t5 + 5*t4/2;
|
||||
Real z7 = t7/6 - t6/2 + t5/2 - t4/6;
|
||||
|
||||
return z0*y0 + z1*dx_*v0 + z2*dx2*a0 + z3*dx3*j0 + z4*y1 + z5*dx_*v1 + z6*dx2*a1 + z7*dx3*j1;
|
||||
}
|
||||
|
||||
inline Real prime(Real x) const {
|
||||
Real xf = x0_ + (data_.size()-1)*dx_;
|
||||
if (x < x0_ || x > xf)
|
||||
{
|
||||
std::ostringstream oss;
|
||||
oss.precision(std::numeric_limits<Real>::digits10+3);
|
||||
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
|
||||
<< x0_ << ", " << xf << "]";
|
||||
throw std::domain_error(oss.str());
|
||||
}
|
||||
if (x == xf)
|
||||
{
|
||||
return data_.back()[1];
|
||||
}
|
||||
|
||||
return this->unchecked_prime(x);
|
||||
}
|
||||
|
||||
inline Real unchecked_prime(Real x) const {
|
||||
return std::numeric_limits<Real>::quiet_NaN();
|
||||
}
|
||||
|
||||
private:
|
||||
RandomAccessContainer data_;
|
||||
Real x0_;
|
||||
Real dx_;
|
||||
};
|
||||
|
||||
|
||||
}
|
||||
#endif
|
||||
#endif
|
||||
|
||||
@@ -34,5 +34,28 @@ public:
|
||||
private:
|
||||
std::shared_ptr<detail::septic_hermite_detail<RandomAccessContainer>> impl_;
|
||||
};
|
||||
|
||||
template<class RandomAccessContainer>
|
||||
class cardinal_septic_hermite {
|
||||
public:
|
||||
using Real = typename RandomAccessContainer::value_type;
|
||||
cardinal_septic_hermite(RandomAccessContainer && y, RandomAccessContainer && dydx,
|
||||
RandomAccessContainer && d2ydx2, RandomAccessContainer && d3ydx3, Real x0, Real dx)
|
||||
: impl_(std::make_shared<detail::cardinal_septic_hermite_detail<RandomAccessContainer>>(
|
||||
std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx))
|
||||
{}
|
||||
|
||||
Real operator()(Real x) const {
|
||||
return impl_->operator()(x);
|
||||
}
|
||||
|
||||
Real prime(Real x) const {
|
||||
return impl_->prime(x);
|
||||
}
|
||||
|
||||
private:
|
||||
std::shared_ptr<detail::cardinal_septic_hermite_detail<RandomAccessContainer>> impl_;
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
@@ -13,7 +13,9 @@
|
||||
#include <thread>
|
||||
#include <future>
|
||||
#include <iostream>
|
||||
#if BOOST_HAS_FLOAT128
|
||||
#include <boost/multiprecision/float128.hpp>
|
||||
#endif
|
||||
#include <boost/math/constants/constants.hpp>
|
||||
#include <boost/math/special_functions/detail/daubechies_scaling_integer_grid.hpp>
|
||||
#include <boost/math/filters/daubechies.hpp>
|
||||
@@ -89,7 +91,7 @@ public:
|
||||
|
||||
matched_holder(RandomAccessContainer && y, RandomAccessContainer && dydx, int grid_refinements) : y_{std::move(y)}, dydx_{std::move(dydx)}
|
||||
{
|
||||
h_ = Real(1)/(1<< grid_refinements);
|
||||
inv_h_ = 1<< grid_refinements;
|
||||
}
|
||||
|
||||
inline Real operator()(Real x) const {
|
||||
@@ -97,25 +99,27 @@ public:
|
||||
using std::floor;
|
||||
using std::sqrt;
|
||||
using std::pow;
|
||||
if (x <= 0 || x >= 3) {
|
||||
return 0;
|
||||
}
|
||||
// This is the exact Holder exponent, but it's pessimistic almost everywhere!
|
||||
// It's only exactly right at dyadic rationals.
|
||||
// Some compilers do not evaluate this at compile time:
|
||||
Real const alpha = 2 - log(1+sqrt(Real(3)))/log(Real(2));
|
||||
Real const onemalpham1 = 1/(1-alpha);
|
||||
//const Real alpha = 0.5500156865235042;
|
||||
// So we're gonna make the graph dip a little harder; this will capture more of the self-similar behavior:
|
||||
//Real constexpr const alpha = Real(3)/Real(10);
|
||||
int64_t i = static_cast<int64_t>(floor(x/h_));
|
||||
Real t = (x- i*h_)/h_;
|
||||
Real s = x*inv_h_;
|
||||
Real ii = floor(s);
|
||||
int64_t i = static_cast<int64_t>(ii);
|
||||
Real t = s - ii;
|
||||
Real v = y_[i];
|
||||
Real dphi = dydx_[i+1]*h_;
|
||||
v += (dphi - alpha*(y_[i+1] - y_[i]))*t/(1-alpha);
|
||||
v += (-dphi + y_[i+1] - y_[i])*pow(t, alpha)/(1-alpha);
|
||||
Real dphi = dydx_[i+1]/inv_h_;
|
||||
Real diff = y_[i+1] - y_[i];
|
||||
v += ((dphi - alpha*diff)*t + (-dphi + diff)*pow(t, alpha) )*onemalpham1;
|
||||
return v;
|
||||
}
|
||||
|
||||
private:
|
||||
Real h_;
|
||||
Real inv_h_;
|
||||
RandomAccessContainer y_;
|
||||
RandomAccessContainer dydx_;
|
||||
};
|
||||
@@ -155,7 +159,9 @@ public:
|
||||
daubechies_scaling(int grid_refinements = -1)
|
||||
{
|
||||
static_assert(p <= 15, "Scaling functions only implements up to p = 15");
|
||||
#if BOOST_HAS_FLOAT128
|
||||
using boost::multiprecision::float128;
|
||||
#endif
|
||||
if (grid_refinements < 0)
|
||||
{
|
||||
if (std::is_same_v<Real, float>)
|
||||
@@ -275,13 +281,6 @@ public:
|
||||
auto y = t0.get();
|
||||
auto dydx = t1.get();
|
||||
|
||||
// Storing the vector of x's is unnecessary; it's only because I haven't implemented an equispaced cubic Hermite interpolator:
|
||||
std::vector<Real> x(y.size());
|
||||
Real h = Real(2*p-1)/(x.size()-1);
|
||||
for (size_t i = 0; i < x.size(); ++i) {
|
||||
x[i] = i*h;
|
||||
}
|
||||
|
||||
if constexpr (p==2) {
|
||||
m_mh = std::make_shared<detail::matched_holder<std::vector<Real>>>(std::move(y), std::move(dydx), grid_refinements);
|
||||
}
|
||||
@@ -297,9 +296,9 @@ public:
|
||||
m_qh = std::make_shared<interpolators::detail::cardinal_quintic_hermite_detail<std::vector<Real>>>(std::move(y), std::move(dydx), std::move(d2ydx2), Real(0), dx);
|
||||
}
|
||||
if constexpr (p >= 10) {
|
||||
m_sh = std::make_shared<interpolators::detail::septic_hermite_detail<std::vector<Real>>>(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
Real dx = Real(1)/(1 << grid_refinements);
|
||||
m_sh = std::make_shared<interpolators::detail::cardinal_septic_hermite_detail<std::vector<Real>>>(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), Real(0), dx);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
@@ -320,7 +319,7 @@ public:
|
||||
return m_qh->unchecked_evaluation(x);
|
||||
}
|
||||
if constexpr (p >= 10) {
|
||||
return m_sh->operator()(x);
|
||||
return m_sh->unchecked_evaluation(x);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -339,7 +338,7 @@ public:
|
||||
return m_qh->unchecked_prime(x);
|
||||
}
|
||||
if constexpr (p >= 10) {
|
||||
return m_sh->prime(x);
|
||||
return m_sh->unchecked_prime(x);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -358,7 +357,7 @@ private:
|
||||
// Need this for p = 6,7,8,9:
|
||||
std::shared_ptr<interpolators::detail::cardinal_quintic_hermite_detail<std::vector<Real>>> m_qh;
|
||||
// Need this for p >= 10:
|
||||
std::shared_ptr<interpolators::detail::septic_hermite_detail<std::vector<Real>>> m_sh;
|
||||
std::shared_ptr<interpolators::detail::cardinal_septic_hermite_detail<std::vector<Real>>> m_sh;
|
||||
|
||||
/*Real constant_interpolation(Real x) const {
|
||||
if (x <= 0 || x >= 2*p-1) {
|
||||
|
||||
288
test/septic_hermite_test.cpp
Normal file
288
test/septic_hermite_test.cpp
Normal file
@@ -0,0 +1,288 @@
|
||||
/*
|
||||
* Copyright Nick Thompson, 2020
|
||||
* 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/random/uniform_real.hpp>
|
||||
#include <boost/random/mersenne_twister.hpp>
|
||||
#include <boost/math/interpolators/septic_hermite.hpp>
|
||||
#ifdef BOOST_HAS_FLOAT128
|
||||
#include <boost/multiprecision/float128.hpp>
|
||||
using boost::multiprecision::float128;
|
||||
#endif
|
||||
|
||||
|
||||
using boost::math::interpolators::septic_hermite;
|
||||
using boost::math::interpolators::cardinal_septic_hermite;
|
||||
|
||||
template<typename Real>
|
||||
void test_constant()
|
||||
{
|
||||
|
||||
std::vector<Real> x{0,1,2,3, 9, 22, 81};
|
||||
std::vector<Real> y(x.size());
|
||||
std::vector<Real> dydx(x.size(), 0);
|
||||
std::vector<Real> d2ydx2(x.size(), 0);
|
||||
std::vector<Real> d3ydx3(x.size(), 0);
|
||||
for (auto & t : y) {
|
||||
t = 7;
|
||||
}
|
||||
|
||||
auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
|
||||
for (Real t = 0; t <= 81; t += 0.25) {
|
||||
CHECK_ULP_CLOSE(Real(7), sh(t), 24);
|
||||
CHECK_ULP_CLOSE(Real(0), sh.prime(t), 24);
|
||||
}
|
||||
|
||||
Real x0 = 0;
|
||||
Real dx = 1;
|
||||
y.resize(128, 7);
|
||||
dydx.resize(128, 0);
|
||||
d2ydx2.resize(128, 0);
|
||||
d3ydx3.resize(128, 0);
|
||||
auto csh = cardinal_septic_hermite(std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3), x0, dx);
|
||||
for (Real t = x0; t <= 127; t += 0.25) {
|
||||
CHECK_ULP_CLOSE(Real(7), csh(t), 24);
|
||||
//CHECK_ULP_CLOSE(Real(0), csh.prime(t), 24);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
template<typename Real>
|
||||
void test_linear()
|
||||
{
|
||||
std::vector<Real> x{0,1,2,3, 4,5,6,7,8,9};
|
||||
std::vector<Real> y = x;
|
||||
std::vector<Real> dydx(x.size(), 1);
|
||||
std::vector<Real> d2ydx2(x.size(), 0);
|
||||
std::vector<Real> d3ydx3(x.size(), 0);
|
||||
|
||||
auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
|
||||
for (Real t = 0; t <= 9; t += 0.25) {
|
||||
CHECK_ULP_CLOSE(Real(t), sh(t), 2);
|
||||
//CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2);
|
||||
}
|
||||
|
||||
boost::random::mt19937 rng;
|
||||
boost::random::uniform_real_distribution<Real> dis(0.5,1);
|
||||
x.resize(512);
|
||||
x[0] = dis(rng);
|
||||
Real xmin = x[0];
|
||||
for (size_t i = 1; i < x.size(); ++i) {
|
||||
x[i] = x[i-1] + dis(rng);
|
||||
}
|
||||
Real xmax = x.back();
|
||||
|
||||
y = x;
|
||||
dydx.resize(x.size(), 1);
|
||||
d2ydx2.resize(x.size(), 0);
|
||||
d3ydx3.resize(x.size(), 0);
|
||||
|
||||
sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
|
||||
for (Real t = xmin; t <= xmax; t += 0.125) {
|
||||
CHECK_ULP_CLOSE(t, sh(t), 15);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Real>
|
||||
void test_quadratic()
|
||||
{
|
||||
std::vector<Real> x{0,1,2,3,4,5,6,7,8,9};
|
||||
std::vector<Real> y(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i)
|
||||
{
|
||||
y[i] = x[i]*x[i]/2;
|
||||
}
|
||||
|
||||
std::vector<Real> dydx(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i) {
|
||||
dydx[i] = x[i];
|
||||
}
|
||||
|
||||
std::vector<Real> d2ydx2(x.size(), 1);
|
||||
std::vector<Real> d3ydx3(x.size(), 0);
|
||||
|
||||
auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
|
||||
for (Real t = 0; t <= 9; t += 0.0078125) {
|
||||
CHECK_ULP_CLOSE(t*t/2, sh(t), 100);
|
||||
//CHECK_ULP_CLOSE(t, qh.prime(t), 7);
|
||||
}
|
||||
|
||||
boost::random::mt19937 rng;
|
||||
boost::random::uniform_real_distribution<Real> dis(0.5,1);
|
||||
x.resize(8);
|
||||
x[0] = dis(rng);
|
||||
Real xmin = x[0];
|
||||
for (size_t i = 1; i < x.size(); ++i) {
|
||||
x[i] = x[i-1] + dis(rng);
|
||||
}
|
||||
Real xmax = x.back();
|
||||
|
||||
y.resize(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i)
|
||||
{
|
||||
y[i] = x[i]*x[i]/2;
|
||||
}
|
||||
|
||||
dydx.resize(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i) {
|
||||
dydx[i] = x[i];
|
||||
}
|
||||
|
||||
d2ydx2.resize(x.size(), 1);
|
||||
d3ydx3.resize(x.size(), 0);
|
||||
|
||||
sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
|
||||
for (Real t = xmin; t <= xmax; t += 0.125) {
|
||||
CHECK_ULP_CLOSE(t*t/2, sh(t), 24);
|
||||
//CHECK_ULP_CLOSE(t, qh.prime(t), 36);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Real>
|
||||
void test_cubic()
|
||||
{
|
||||
|
||||
std::vector<Real> x{0,1,2,3,4,5,6,7};
|
||||
Real xmax = x.back();
|
||||
std::vector<Real> y(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i)
|
||||
{
|
||||
y[i] = x[i]*x[i]*x[i];
|
||||
}
|
||||
|
||||
std::vector<Real> dydx(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i) {
|
||||
dydx[i] = 3*x[i]*x[i];
|
||||
}
|
||||
|
||||
std::vector<Real> d2ydx2(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i) {
|
||||
d2ydx2[i] = 6*x[i];
|
||||
}
|
||||
std::vector<Real> d3ydx3(x.size(), 6);
|
||||
|
||||
auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
|
||||
for (Real t = 0; t <= xmax; t += 0.0078125) {
|
||||
CHECK_ULP_CLOSE(t*t*t, sh(t), 151);
|
||||
}
|
||||
}
|
||||
|
||||
template<typename Real>
|
||||
void test_quartic()
|
||||
{
|
||||
|
||||
std::vector<Real> x{0,1,2,3,4,5,6,7,8,9};
|
||||
Real xmax = x.back();
|
||||
std::vector<Real> y(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i)
|
||||
{
|
||||
y[i] = x[i]*x[i]*x[i]*x[i];
|
||||
}
|
||||
|
||||
std::vector<Real> dydx(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i) {
|
||||
dydx[i] = 4*x[i]*x[i]*x[i];
|
||||
}
|
||||
|
||||
std::vector<Real> d2ydx2(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i) {
|
||||
d2ydx2[i] = 12*x[i]*x[i];
|
||||
}
|
||||
|
||||
std::vector<Real> d3ydx3(x.size());
|
||||
for (size_t i = 0; i < y.size(); ++i) {
|
||||
d3ydx3[i] = 24*x[i];
|
||||
}
|
||||
|
||||
auto sh = septic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2), std::move(d3ydx3));
|
||||
|
||||
for (Real t = 1; t <= xmax; t += 0.0078125) {
|
||||
CHECK_ULP_CLOSE(t*t*t*t, sh(t), 117);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
template<typename Real>
|
||||
void test_interpolation_condition()
|
||||
{
|
||||
for (size_t n = 4; n < 50; ++n) {
|
||||
std::vector<Real> x(n);
|
||||
std::vector<Real> y(n);
|
||||
std::vector<Real> dydx(n);
|
||||
std::vector<Real> d2ydx2(n);
|
||||
std::vector<Real> d3ydx3(n);
|
||||
boost::random::mt19937 rd;
|
||||
boost::random::uniform_real_distribution<Real> dis(0,1);
|
||||
Real x0 = dis(rd);
|
||||
x[0] = x0;
|
||||
y[0] = dis(rd);
|
||||
for (size_t i = 1; i < n; ++i) {
|
||||
x[i] = x[i-1] + dis(rd);
|
||||
y[i] = dis(rd);
|
||||
dydx[i] = dis(rd);
|
||||
d2ydx2[i] = dis(rd);
|
||||
d3ydx3[i] = dis(rd);
|
||||
}
|
||||
|
||||
auto x_copy = x;
|
||||
auto y_copy = y;
|
||||
auto dydx_copy = dydx;
|
||||
auto d2ydx2_copy = d2ydx2;
|
||||
auto d3ydx3_copy = d3ydx3;
|
||||
auto s = septic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy), std::move(d2ydx2_copy), std::move(d3ydx3_copy));
|
||||
|
||||
for (size_t i = 0; i < x.size(); ++i) {
|
||||
CHECK_ULP_CLOSE(y[i], s(x[i]), 2);
|
||||
//CHECK_ULP_CLOSE(dydx[i], s.prime(x[i]), 2);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
int main()
|
||||
{
|
||||
test_constant<float>();
|
||||
test_linear<float>();
|
||||
test_quadratic<float>();
|
||||
test_cubic<float>();
|
||||
test_quartic<float>();
|
||||
test_interpolation_condition<float>();
|
||||
|
||||
test_constant<double>();
|
||||
test_linear<double>();
|
||||
test_quadratic<double>();
|
||||
test_cubic<double>();
|
||||
test_quartic<double>();
|
||||
test_interpolation_condition<double>();
|
||||
|
||||
test_constant<long double>();
|
||||
test_linear<long double>();
|
||||
test_quadratic<long double>();
|
||||
test_cubic<long double>();
|
||||
test_quartic<long double>();
|
||||
test_interpolation_condition<long double>();
|
||||
|
||||
#ifdef BOOST_HAS_FLOAT128
|
||||
test_constant<float128>();
|
||||
test_linear<float128>();
|
||||
test_quadratic<float128>();
|
||||
test_cubic<float128>();
|
||||
test_quartic<float128>();
|
||||
test_interpolation_condition<float128>();
|
||||
#endif
|
||||
|
||||
return boost::math::test::report_errors();
|
||||
}
|
||||
Reference in New Issue
Block a user