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

Merge pull request #306 from boostorg/quintic_hermite

Browse Source 
Quintic Hermite interpolation
This commit is contained in:
Nick
2020-01-25 08:08:50 +08:00
committed by GitHub
parent cb1b9c102e f1a90ae686
commit 8213c19b30
6 changed files with 474 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

119
doc/interpolators/quintic_hermite.qbk Normal file
View File

@@ -0,0 +1,119 @@
[/
Copyright (c) 2020 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:quintic_hermite Quintic Hermite interpolation]
[heading Synopsis]
``
#include <boost/math/interpolators/quintic_hermite.hpp>
namespace boost::math::interpolators {
template<class RandomAccessContainer>
class quintic_hermite {
public:
using Real = typename RandomAccessContainer::value_type;
quintic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2)
Real operator()(Real x) const;
Real prime(Real x) const;
friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m);
void push_back(Real x, Real y, Real dydx, Real d2ydx2);
};
}
``
[heading Quintic Hermite Interpolation]
The quintic Hermite interpolator takes a list of possibly non-uniformly spaced abscissas, ordinates, and their velocities and accelerations which are used to construct a quintic interpolating polynomial between segments.
The interpolant is /C/[super 2] and its evaluation has [bigo](log(/N/)) complexity.
An example usage is as follows:
std::vector<double> x{1,2,3, 4, 5, 6};
std::vector<double> y{7,8,9,10,11,12};
std::vector<double> dydx{1,1,1,1,1,1};
std::vector<double> d2ydx2{0,0,0,0,0,0};
using boost::math::interpolators::quintic_hermite;
auto spline = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
// evaluate at a point:
double z = spline(3.4);
// evaluate derivative at a point:
double zprime = spline.prime(3.4);
Periodically, it is helpful to see what data the interpolator has.
This can be achieved via
std::cout << spline << "\n";
Note that the interpolator is pimpl'd, so that copying the class is cheap, and hence it can be shared between threads.
(The call operator and `.prime()` are threadsafe.)
The interpolator can be updated in constant time.
Hence we can use `boost::circular_buffer` to do real-time interpolation.
#include <boost/circular_buffer.hpp>
...
boost::circular_buffer<double> initial_x{1,2,3,4};
boost::circular_buffer<double> initial_y{4,5,6,7};
boost::circular_buffer<double> initial_dydx{1,1,1,1};
boost::circular_buffer<double> initial_d2ydx2{0,0,0,0};
auto circular_akima = quintic_hermite(std::move(initial_x), std::move(initial_y), std::move(initial_dydx), std::move(initial_d2ydx2));
// interpolate via call operation:
double y = circular_akima(3.5);
// add new data:
circular_akima.push_back(5, 8, 1, 0);
// interpolate at 4.5:
y = circular_akima(4.5);
[heading Complexity and Performance]
The following google benchmark demonstrates the cost of the call operator for this interpolator:
```
Run on (16 X 4300 MHz CPU s)
CPU Caches:
L1 Data 32K (x8)
L1 Instruction 32K (x8)
L2 Unified 1024K (x8)
L3 Unified 11264K (x1)
Load Average: 0.92, 0.64, 0.35
-----------------------------------------------
Benchmark Time
-----------------------------------------------
BMQuinticHermite<double>/8 8.14 ns
BMQuinticHermite<double>/16 8.69 ns
BMQuinticHermite<double>/32 9.42 ns
BMQuinticHermite<double>/64 9.90 ns
BMQuinticHermite<double>/128 10.4 ns
BMQuinticHermite<double>/256 10.9 ns
BMQuinticHermite<double>/512 11.6 ns
BMQuinticHermite<double>/1024 12.3 ns
BMQuinticHermite<double>/2048 12.8 ns
BMQuinticHermite<double>/4096 13.6 ns
BMQuinticHermite<double>/8192 14.6 ns
BMQuinticHermite<double>/16384 15.5 ns
BMQuinticHermite<double>/32768 17.4 ns
BMQuinticHermite<double>/65536 18.5 ns
BMQuinticHermite<double>/131072 18.8 ns
BMQuinticHermite<double>/262144 19.8 ns
BMQuinticHermite<double>/524288 20.5 ns
BMQuinticHermite<double>/1048576 21.6 ns
BMQuinticHermite<double>/2097152 22.5 ns
BMQuinticHermite<double>/4194304 24.2 ns
BMQuinticHermite<double>/8388608 26.6 ns
BMQuinticHermite<double>/16777216 26.7 ns
BMQuinticHermite<double>_BigO 1.14 lgN
```
[endsect]
[/section:quintic_hermite]

1
doc/math.qbk
View File

@@ -715,6 +715,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22.
[include interpolators/cubic_hermite.qbk]
[include interpolators/makima.qbk]
[include interpolators/pchip.qbk]
[include interpolators/quintic_hermite.qbk]
[endmathpart]
[mathpart quadrature Quadrature and Differentiation]

135
include/boost/math/interpolators/detail/quintic_hermite_detail.hpp Normal file
View File

@@ -0,0 +1,135 @@
#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP
#define BOOST_MATH_INTERPOLATORS_DETAIL_QUINTIC_HERMITE_DETAIL_HPP
#include <algorithm>
#include <stdexcept>
#include <sstream>
#include <cmath>
namespace boost::math::interpolators::detail {
template<class RandomAccessContainer>
class quintic_hermite_detail {
public:
using Real = typename RandomAccessContainer::value_type;
quintic_hermite_detail(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2) : x_{std::move(x)}, y_{std::move(y)}, dydx_{std::move(dydx)}, d2ydx2_{std::move(d2ydx2)}
{
if (x_.size() != y_.size()) {
throw std::domain_error("Number of abscissas must = number of ordinates.");
}
if (x_.size() != dydx_.size()) {
throw std::domain_error("Numbers of derivatives must = number of abscissas.");
}
if (x_.size() != d2ydx2_.size()) {
throw std::domain_error("Number of second derivatives must equal number of abscissas.");
}
if (x_.size() < 2) {
throw std::domain_error("At least 2 abscissas are required.");
}
Real x0 = x_[0];
for (decltype(x_.size()) i = 1; i < x_.size(); ++i) {
Real x1 = x_[i];
if (x1 <= x0)
{
throw std::domain_error("Abscissas must be sorted in strictly increasing order x0 < x1 < ... < x_{n-1}");
}
x0 = x1;
}
}
void push_back(Real x, Real y, Real dydx, Real d2ydx2) {
using std::abs;
using std::isnan;
if (x <= x_.back()) {
throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
}
x_.push_back(x);
y_.push_back(y);
dydx_.push_back(dydx);
d2ydx2_.push_back(d2ydx2);
}
Real operator()(Real x) const {
if (x < x_[0] || x > x_.back()) {
std::ostringstream oss;
oss.precision(std::numeric_limits<Real>::digits10+3);
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
<< x_[0] << ", " << x_.back() << "]";
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 == x_.back()) {
return y_.back();
}
auto it = std::upper_bound(x_.begin(), x_.end(), x);
auto i = std::distance(x_.begin(), it) -1;
Real x0 = *(it-1);
Real x1 = *it;
Real y0 = y_[i];
Real y1 = y_[i+1];
Real v0 = dydx_[i];
Real v1 = dydx_[i+1];
Real a0 = d2ydx2_[i];
Real a1 = d2ydx2_[i+1];
Real dx = (x1-x0);
Real t = (x-x0)/dx;
// See the 'Basis functions' section of:
// https://www.rose-hulman.edu/~finn/CCLI/Notes/day09.pdf
// Also: https://github.com/MrHexxx/QuinticHermiteSpline/blob/master/HermiteSpline.cs
Real y = (1- t*t*t*(10 + t*(-15 + 6*t)))*y0;
y += t*(1+ t*t*(-6 + t*(8 -3*t)))*v0;
y += t*t*(1 + t*(-3 + t*(3-t)))*a0/2;
y += t*t*t*((1 + t*(-2 + t))*a1/2 + (-4 + t*(7 -3*t))*v1 + (10 + t*(-15 + 6*t))*y1);
return y;
}
Real prime(Real x) const {
if (x < x_[0] || x > x_.back()) {
std::ostringstream oss;
oss.precision(std::numeric_limits<Real>::digits10+3);
oss << "Requested abscissa x = " << x << ", which is outside of allowed range ["
<< x_[0] << ", " << x_.back() << "]";
throw std::domain_error(oss.str());
}
if (x == x_.back()) {
return dydx_.back();
}
auto it = std::upper_bound(x_.begin(), x_.end(), x);
auto i = std::distance(x_.begin(), it) -1;
Real x0 = *(it-1);
Real x1 = *it;
Real s0 = dydx_[i];
Real s1 = dydx_[i+1];
// Ridiculous linear interpolation. Fine for now:
Real numerator = s0*(x1-x) + s1*(x-x0);
Real denominator = x1 - x0;
return numerator/denominator;
}
friend std::ostream& operator<<(std::ostream & os, const quintic_hermite_detail & m)
{
os << "(x,y,y') = {";
for (size_t i = 0; i < m.x_.size() - 1; ++i) {
os << "(" << m.x_[i] << ", " << m.y_[i] << ", " << m.dydx_[i] << ", " << m.d2ydx2_[i] << "), ";
}
auto n = m.x_.size()-1;
os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.dydx_[n] << ", " << m.d2ydx2_[n] << ")}";
return os;
}
private:
RandomAccessContainer x_;
RandomAccessContainer y_;
RandomAccessContainer dydx_;
RandomAccessContainer d2ydx2_;
};
}
#endif

41
include/boost/math/interpolators/quintic_hermite.hpp Normal file
View File

@@ -0,0 +1,41 @@
#ifndef BOOST_MATH_INTERPOLATORS_QUINTIC_HERMITE_HPP
#define BOOST_MATH_INTERPOLATORS_QUINTIC_HERMITE_HPP
#include <algorithm>
#include <stdexcept>
#include <memory>
#include <boost/math/interpolators/detail/quintic_hermite_detail.hpp>
namespace boost::math::interpolators {
template<class RandomAccessContainer>
class quintic_hermite {
public:
using Real = typename RandomAccessContainer::value_type;
quintic_hermite(RandomAccessContainer && x, RandomAccessContainer && y, RandomAccessContainer && dydx, RandomAccessContainer && d2ydx2)
: impl_(std::make_shared<detail::quintic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2)))
{}
Real operator()(Real x) const {
return impl_->operator()(x);
}
Real prime(Real x) const {
return impl_->prime(x);
}
friend std::ostream& operator<<(std::ostream & os, const quintic_hermite & m)
{
os << *m.impl_;
return os;
}
void push_back(Real x, Real y, Real dydx, Real d2ydx2) {
impl_->push_back(x, y, dydx, d2ydx2);
}
private:
std::shared_ptr<detail::quintic_hermite_detail<RandomAccessContainer>> impl_;
};
}
#endif

1
test/Jamfile.v2
View File

@@ -947,6 +947,7 @@ test-suite misc :
[ 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" : <linkflags>-lquadmath ] ]
[ run makima_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
[ run pchip_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
[ run quintic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
[ run cubic_hermite_test.cpp : : : [ requires cxx17_if_constexpr cxx17_std_apply ] [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : <linkflags>-lquadmath ] ]
[ run catmull_rom_test.cpp ../../test/build//boost_unit_test_framework : : : <define>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 : : : <define>TEST=2 [ requires cxx11_hdr_array cxx11_hdr_initializer_list ] : catmull_rom_test_2 ]

177
test/quintic_hermite_test.cpp Normal file
View File

@@ -0,0 +1,177 @@
/*
* 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 <random>
#include <boost/math/interpolators/quintic_hermite.hpp>
#include <boost/circular_buffer.hpp>
#ifdef BOOST_HAS_FLOAT128
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
#endif
using boost::math::interpolators::quintic_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);
for (auto & t : y) {
t = 7;
}
auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
for (Real t = 0; t <= 81; t += 0.25) {
CHECK_ULP_CLOSE(Real(7), qh(t), 24);
CHECK_ULP_CLOSE(Real(0), qh.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);
auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
for (Real t = 0; t <= 9; t += 0.25) {
CHECK_ULP_CLOSE(Real(t), qh(t), 2);
CHECK_ULP_CLOSE(Real(1), qh.prime(t), 2);
}
}
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);
auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
for (Real t = 0; t <= 9; t += 0.125) {
CHECK_ULP_CLOSE(Real(t*t)/2, qh(t), 2);
CHECK_ULP_CLOSE(t, qh.prime(t), 2);
}
}
template<typename Real>
void test_cubic()
{
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]*x[i]/6;
}
std::vector<Real> dydx(x.size());
for (size_t i = 0; i < y.size(); ++i) {
dydx[i] = x[i]*x[i]/2;
}
std::vector<Real> d2ydx2(x.size());
for (size_t i = 0; i < y.size(); ++i) {
d2ydx2[i] = x[i];
}
auto qh = quintic_hermite(std::move(x), std::move(y), std::move(dydx), std::move(d2ydx2));
for (Real t = 0; t <= 9; t += 0.125) {
CHECK_ULP_CLOSE(Real(t*t*t)/6, qh(t), 5);
//CHECK_ULP_CLOSE(t, qh.prime(t), 2);
}
}
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::default_random_engine rd;
std::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);
}
auto x_copy = x;
auto y_copy = y;
auto dydx_copy = dydx;
auto d2ydx2_copy = d2ydx2;
auto s = quintic_hermite(std::move(x_copy), std::move(y_copy), std::move(dydx_copy), std::move(d2ydx2_copy));
//std::cout << "s = " << s << "\n";
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_interpolation_condition<float>();
test_constant<double>();
test_linear<double>();
test_quadratic<double>();
test_cubic<double>();
test_interpolation_condition<double>();
test_constant<long double>();
test_linear<long double>();
test_quadratic<long double>();
test_cubic<long double>();
test_interpolation_condition<long double>();
#ifdef BOOST_HAS_FLOAT128
test_constant<float128>();
test_linear<float128>();
test_quadratic<float128>();
test_cubic<float128>();
#endif
return boost::math::test::report_errors();
}