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

Cubic Hermite spline: Backend pchip and makima to cubic_hermite.

Browse Source
This commit is contained in:
Nick Thompson
2020-01-18 14:11:31 -05:00
parent d8c2219a23
commit 2a37abd93a
7 changed files with 216 additions and 510 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
doc/interpolators/makima.qbk
View File

@@ -8,9 +8,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
[section:makima Modified Akima interpolation]
[heading Synopsis]
``
#include <boost/math/interpolators/makima.hpp>
``
#include <boost/math/interpolators/makima.hpp>
namespace boost::math::interpolators {
@@ -87,47 +86,7 @@ The modified Akima spline oscillates less than the cubic spline, but has less sm
[heading Complexity and Performance]
The call to the constructor requires [bigo](/N/) operations to compute the weighted slopes.
Each call to the interpolant is [bigo](log(/N/)), where /N/ is the number of points to interpolate.
A google benchmark demonstrating the performance of evaluating the spline is given below:
```
Run on (4 X 2700 MHz CPU s)
CPU Caches:
L1 Data 32K (x2)
L1 Instruction 32K (x2)
L2 Unified 262K (x2)
L3 Unified 3145K (x1)
Load Average: 2.08, 1.83, 1.96
---------------------------------------
Benchmark Time
---------------------------------------
BMMakima<double>/8 11.3 ns
BMMakima<double>/16 11.5 ns
BMMakima<double>/32 12.7 ns
BMMakima<double>/64 13.6 ns
BMMakima<double>/128 14.5 ns
BMMakima<double>/256 15.1 ns
BMMakima<double>/512 17.2 ns
BMMakima<double>/1024 18.3 ns
BMMakima<double>/2048 19.6 ns
BMMakima<double>/4096 20.7 ns
BMMakima<double>/8192 22.1 ns
BMMakima<double>/16384 23.1 ns
BMMakima<double>/32768 24.2 ns
BMMakima<double>/65536 25.3 ns
BMMakima<double>/131072 27.1 ns
BMMakima<double>/262144 28.3 ns
BMMakima<double>/524288 29.9 ns
BMMakima<double>/1048576 31.6 ns
BMMakima<double>/2097152 31.8 ns
BMMakima<double>/4194304 33.7 ns
BMMakima<double>/8388608 35.0 ns
BMMakima<double>/16777216 40.0 ns
BMMakima<double>_BigO 1.63 lgN
```
The complexity and performance is identical to that of the cubic Hermite interpolator, since this object simply constructs derivatives and forwards the data to `cubic_hermite.hpp`.
[endsect]
[/section:makima]

11
doc/interpolators/pchip.qbk
View File

@@ -8,9 +8,8 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
[section:pchip PCHIP interpolation]
[heading Synopsis]
``
#include <boost/math/interpolators/pchip.hpp>
``
#include <boost/math/interpolators/pchip.hpp>
namespace boost::math::interpolators {
@@ -82,11 +81,7 @@ Hence we can use `boost::circular_buffer` to do real-time interpolation:
[heading Complexity and Performance]
The call to the constructor requires [bigo](/N/) operations to compute the weighted slopes.
Each call to the interpolant is [bigo](log(/N/)), where /N/ is the number of points to interpolate.
A google benchmark demonstrating the performance of evaluating the spline is given below:
This interpolator chooses the slopes and forwards data to the cubic Hermite interpolator, so the performance is stated in the documentation for `cubic_hermite.hpp`.
[endsect]

5
include/boost/math/interpolators/detail/cubic_hermite_detail.hpp
View File

@@ -130,7 +130,10 @@ public:
return os;
}
private:
auto size() const {
return x_.size();
}
RandomAccessContainer x_;
RandomAccessContainer y_;
RandomAccessContainer dydx_;

248
include/boost/math/interpolators/detail/makima_detail.hpp
View File

@@ -1,248 +0,0 @@
// 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)
// See: https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/
// And: https://doi.org/10.1145/321607.321609
#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_MAKIMA_DETAIL_HPP
#define BOOST_MATH_INTERPOLATORS_DETAIL_MAKIMA_DETAIL_HPP
#include <stdexcept>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <sstream>
#include <limits>
namespace boost::math::interpolators::detail {
template<class RandomAccessContainer>
class makima_detail {
public:
using Real = typename RandomAccessContainer::value_type;
makima_detail(RandomAccessContainer && x, RandomAccessContainer && y,
Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()) : x_{std::move(x)}, y_{std::move(y)}
{
using std::abs;
using std::isnan;
if (x_.size() != y_.size())
{
throw std::domain_error("There must be the same number of ordinates as abscissas.");
}
if (x_.size() < 4)
{
throw std::domain_error("Must be at least four data points.");
}
Real x0 = x_[0];
for (size_t i = 1; i < x_.size(); ++i) {
Real x1 = x_[i];
if (x1 <= x0) {
throw std::domain_error("Abscissas must be listed in strictly increasing order x0 < x1 < ... < x_{n-1}");
}
x0 = x1;
}
s_.resize(x_.size(), std::numeric_limits<Real>::quiet_NaN());
Real m2 = (y_[3]-y_[2])/(x_[3]-x_[2]);
Real m1 = (y_[2]-y_[1])/(x_[2]-x_[1]);
Real m0 = (y_[1]-y_[0])/(x_[1]-x_[0]);
// Quadratic extrapolation: m_{-1} = 2m_0 - m_1:
Real mm1 = 2*m0 - m1;
// Quadratic extrapolation: m_{-2} = 2*m_{-1}-m_0:
Real mm2 = 2*mm1 - m0;
Real w1 = abs(m1-m0) + abs(m1+m0)/2;
Real w2 = abs(mm1-mm2) + abs(mm1+mm2)/2;
if (isnan(left_endpoint_derivative))
{
s_[0] = (w1*mm1 + w2*m0)/(w1+w2);
if (isnan(s_[0]))
{
s_[0] = 0;
}
}
else
{
s_[0] = left_endpoint_derivative;
}
w1 = abs(m2-m1) + abs(m2+m1)/2;
w2 = abs(m0-mm1) + abs(m0+mm1)/2;
s_[1] = (w1*m0 + w2*m1)/(w1+w2);
if (isnan(s_[1])) {
s_[1] = 0;
}
for (decltype(s_.size()) i = 2; i < s_.size()-2; ++i) {
Real mim2 = (y_[i-1]-y_[i-2])/(x_[i-1]-x_[i-2]);
Real mim1 = (y_[i ]-y_[i-1])/(x_[i ]-x_[i-1]);
Real mi = (y_[i+1]-y_[i ])/(x_[i+1]-x_[i ]);
Real mip1 = (y_[i+2]-y_[i+1])/(x_[i+2]-x_[i+1]);
w1 = abs(mip1-mi) + abs(mip1+mi)/2;
w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
s_[i] = (w1*mim1 + w2*mi)/(w1+w2);
if (isnan(s_[i])) {
s_[i] = 0;
}
}
// Quadratic extrapolation at the other end:
decltype(s_.size()) n = s_.size();
Real mnm4 = (y_[n-3]-y_[n-4])/(x_[n-3]-x_[n-4]);
Real mnm3 = (y_[n-2]-y_[n-3])/(x_[n-2]-x_[n-3]);
Real mnm2 = (y_[n-1]-y_[n-2])/(x_[n-1]-x_[n-2]);
Real mnm1 = 2*mnm2 - mnm3;
Real mn = 2*mnm1 - mnm2;
w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
s_[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
if (isnan(s_[n-2])) {
s_[n-2] = 0;
}
w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
if (isnan(right_endpoint_derivative))
{
s_[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
if (isnan(s_[n-1])) {
s_[n-1] = 0;
}
}
else
{
s_[n-1] = right_endpoint_derivative;
}
}
void push_back(Real x, Real y) {
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);
s_.push_back(std::numeric_limits<Real>::quiet_NaN());
// s_[n-2] was computed by extrapolation. Now s_[n-2] -> s_[n-3], and it can be computed by the same formula.
decltype(s_.size()) n = s_.size();
auto i = n - 3;
Real mim2 = (y_[i-1]-y_[i-2])/(x_[i-1]-x_[i-2]);
Real mim1 = (y_[i ]-y_[i-1])/(x_[i ]-x_[i-1]);
Real mi = (y_[i+1]-y_[i ])/(x_[i+1]-x_[i ]);
Real mip1 = (y_[i+2]-y_[i+1])/(x_[i+2]-x_[i+1]);
Real w1 = abs(mip1-mi) + abs(mip1+mi)/2;
Real w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
s_[i] = (w1*mim1 + w2*mi)/(w1+w2);
if (isnan(s_[i])) {
s_[i] = 0;
}
Real mnm4 = (y_[n-3]-y_[n-4])/(x_[n-3]-x_[n-4]);
Real mnm3 = (y_[n-2]-y_[n-3])/(x_[n-2]-x_[n-3]);
Real mnm2 = (y_[n-1]-y_[n-2])/(x_[n-1]-x_[n-2]);
Real mnm1 = 2*mnm2 - mnm3;
Real mn = 2*mnm1 - mnm2;
w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
s_[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
if (isnan(s_[n-2])) {
s_[n-2] = 0;
}
w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
s_[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
if (isnan(s_[n-1])) {
s_[n-1] = 0;
}
}
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 s0 = s_[i];
Real s1 = s_[i+1];
Real dx = (x1-x0);
Real t = (x-x0)/dx;
// See the section 'Representations' in the page
// https://en.wikipedia.org/wiki/Cubic_Hermite_spline
// This uses the factorized form:
//Real y = y0*(1+2*t)*(1-t)*(1-t) + dx*s0*t*(1-t)*(1-t)
// + y1*t*t*(3-2*t) + dx*s1*t*t*(t-1);
// And then factorized further:
Real y = (1-t)*(1-t)*(y0*(1+2*t) + s0*(x-x0))
+ t*t*(y1*(3-2*t) + dx*s1*(t-1));
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 s_.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 = s_[i];
Real s1 = s_[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 makima_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.s_[i] << "), ";
}
auto n = m.x_.size()-1;
os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.s_[n] << ")}";
return os;
}
private:
RandomAccessContainer x_;
RandomAccessContainer y_;
RandomAccessContainer s_;
};
}
#endif

199
include/boost/math/interpolators/detail/pchip_detail.hpp
View File

@@ -1,199 +0,0 @@
// 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)
// See Fritsch and Carlson: https://doi.org/10.1137/0717021
#ifndef BOOST_MATH_INTERPOLATORS_DETAIL_PCHIP_DETAIL_HPP
#define BOOST_MATH_INTERPOLATORS_DETAIL_PCHIP_DETAIL_HPP
#include <stdexcept>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <sstream>
#include <limits>
namespace boost::math::interpolators::detail {
template<class RandomAccessContainer>
class pchip_detail {
public:
using Real = typename RandomAccessContainer::value_type;
pchip_detail(RandomAccessContainer && x, RandomAccessContainer && y,
Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()) : x_{std::move(x)}, y_{std::move(y)}
{
using std::abs;
using std::isnan;
if (x_.size() != y_.size())
{
throw std::domain_error("There must be the same number of ordinates as abscissas.");
}
if (x_.size() < 4)
{
throw std::domain_error("Must be at least four data points.");
}
Real x0 = x_[0];
for (size_t i = 1; i < x_.size(); ++i) {
Real x1 = x_[i];
if (x1 <= x0) {
throw std::domain_error("Abscissas must be listed in strictly increasing order x0 < x1 < ... < x_{n-1}");
}
x0 = x1;
}
s_.resize(x_.size(), std::numeric_limits<Real>::quiet_NaN());
if (isnan(left_endpoint_derivative))
{
// O(h) finite difference derivative:
s_[0] = (y_[1]-y_[0])/(x_[1]-x_[0]);
}
else
{
s_[0] = left_endpoint_derivative;
}
for (decltype(s_.size()) k = 1; k < s_.size()-1; ++k) {
Real hkm1 = x_[k] - x_[k-1];
Real dkm1 = (y_[k] - y_[k-1])/hkm1;
Real hk = x_[k+1] - x_[k];
Real dk = (y_[k+1] - y_[k])/hk;
Real w1 = 2*hk + hkm1;
Real w2 = hk + 2*hkm1;
if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
{
s_[k] = 0;
}
else
{
s_[k] = (w1+w2)/(w1/dkm1 + w2/dk);
}
}
// Quadratic extrapolation at the other end:
decltype(s_.size()) n = s_.size();
if (isnan(right_endpoint_derivative))
{
s_[n-1] = (y_[n-1]-y_[n-2])/(x_[n-1] - x_[n-2]);
}
else
{
s_[n-1] = right_endpoint_derivative;
}
}
void push_back(Real x, Real y) {
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);
s_.push_back(std::numeric_limits<Real>::quiet_NaN());
decltype(s_.size()) n = s_.size();
s_[n-1] = (y_[n-1]-y_[n-2])/(x_[n-1] - x_[n-2]);
// Now fix s_[n-2]:
auto k = n-2;
Real hkm1 = x_[k] - x_[k-1];
Real dkm1 = (y_[k] - y_[k-1])/hkm1;
Real hk = x_[k+1] - x_[k];
Real dk = (y_[k+1] - y_[k])/hk;
Real w1 = 2*hk + hkm1;
Real w2 = hk + 2*hkm1;
if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
{
s_[k] = 0;
}
else
{
s_[k] = (w1+w2)/(w1/dkm1 + w2/dk);
}
}
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 s0 = s_[i];
Real s1 = s_[i+1];
Real dx = (x1-x0);
Real t = (x-x0)/dx;
// See the section 'Representations' in the page
// https://en.wikipedia.org/wiki/Cubic_Hermite_spline
// This uses the factorized form:
//Real y = y0*(1+2*t)*(1-t)*(1-t) + dx*s0*t*(1-t)*(1-t)
// + y1*t*t*(3-2*t) + dx*s1*t*t*(t-1);
// And then factorized further:
Real y = (1-t)*(1-t)*(y0*(1+2*t) + s0*(x-x0))
+ t*t*(y1*(3-2*t) + dx*s1*(t-1));
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 s_.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 = s_[i];
Real s1 = s_[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 pchip_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.s_[i] << "), ";
}
auto n = m.x_.size()-1;
os << "(" << m.x_[n] << ", " << m.y_[n] << ", " << m.s_[n] << ")}";
return os;
}
private:
RandomAccessContainer x_;
RandomAccessContainer y_;
RandomAccessContainer s_;
};
}
#endif

134
include/boost/math/interpolators/makima.hpp
View File

@@ -10,7 +10,8 @@
#ifndef BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
#define BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
#include <memory>
#include <boost/math/interpolators/detail/makima_detail.hpp>
#include <cmath>
#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
namespace boost::math::interpolators {
@@ -21,8 +22,90 @@ public:
makima(RandomAccessContainer && x, RandomAccessContainer && y,
Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()) : impl_(std::make_shared<detail::makima_detail<RandomAccessContainer>>(std::move(x), std::move(y), left_endpoint_derivative, right_endpoint_derivative))
{}
Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
{
using std::isnan;
using std::abs;
if (x.size() < 4)
{
throw std::domain_error("Must be at least four data points.");
}
RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());
Real m2 = (y[3]-y[2])/(x[3]-x[2]);
Real m1 = (y[2]-y[1])/(x[2]-x[1]);
Real m0 = (y[1]-y[0])/(x[1]-x[0]);
// Quadratic extrapolation: m_{-1} = 2m_0 - m_1:
Real mm1 = 2*m0 - m1;
// Quadratic extrapolation: m_{-2} = 2*m_{-1}-m_0:
Real mm2 = 2*mm1 - m0;
Real w1 = abs(m1-m0) + abs(m1+m0)/2;
Real w2 = abs(mm1-mm2) + abs(mm1+mm2)/2;
if (isnan(left_endpoint_derivative))
{
s[0] = (w1*mm1 + w2*m0)/(w1+w2);
if (isnan(s[0]))
{
s[0] = 0;
}
}
else
{
s[0] = left_endpoint_derivative;
}
w1 = abs(m2-m1) + abs(m2+m1)/2;
w2 = abs(m0-mm1) + abs(m0+mm1)/2;
s[1] = (w1*m0 + w2*m1)/(w1+w2);
if (isnan(s[1])) {
s[1] = 0;
}
for (decltype(s.size()) i = 2; i < s.size()-2; ++i) {
Real mim2 = (y[i-1]-y[i-2])/(x[i-1]-x[i-2]);
Real mim1 = (y[i ]-y[i-1])/(x[i ]-x[i-1]);
Real mi = (y[i+1]-y[i ])/(x[i+1]-x[i ]);
Real mip1 = (y[i+2]-y[i+1])/(x[i+2]-x[i+1]);
w1 = abs(mip1-mi) + abs(mip1+mi)/2;
w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
s[i] = (w1*mim1 + w2*mi)/(w1+w2);
if (isnan(s[i])) {
s[i] = 0;
}
}
// Quadratic extrapolation at the other end:
decltype(s.size()) n = s.size();
Real mnm4 = (y[n-3]-y[n-4])/(x[n-3]-x[n-4]);
Real mnm3 = (y[n-2]-y[n-3])/(x[n-2]-x[n-3]);
Real mnm2 = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
Real mnm1 = 2*mnm2 - mnm3;
Real mn = 2*mnm1 - mnm2;
w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
s[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
if (isnan(s[n-2])) {
s[n-2] = 0;
}
w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
if (isnan(right_endpoint_derivative))
{
s[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
if (isnan(s[n-1])) {
s[n-1] = 0;
}
}
else
{
s[n-1] = right_endpoint_derivative;
}
impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));
}
Real operator()(Real x) const {
return impl_->operator()(x);
@@ -39,11 +122,52 @@ public:
}
void push_back(Real x, Real y) {
impl_->push_back(x, y);
using std::abs;
using std::isnan;
if (x <= impl_->x_.back()) {
throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
}
impl_->x_.push_back(x);
impl_->y_.push_back(y);
impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());
// dydx_[n-2] was computed by extrapolation. Now dydx_[n-2] -> dydx_[n-3], and it can be computed by the same formula.
decltype(impl_->size()) n = impl_->size();
auto i = n - 3;
Real mim2 = (impl_->y_[i-1]-impl_->y_[i-2])/(impl_->x_[i-1]-impl_->x_[i-2]);
Real mim1 = (impl_->y_[i ]-impl_->y_[i-1])/(impl_->x_[i ]-impl_->x_[i-1]);
Real mi = (impl_->y_[i+1]-impl_->y_[i ])/(impl_->x_[i+1]-impl_->x_[i ]);
Real mip1 = (impl_->y_[i+2]-impl_->y_[i+1])/(impl_->x_[i+2]-impl_->x_[i+1]);
Real w1 = abs(mip1-mi) + abs(mip1+mi)/2;
Real w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
impl_->dydx_[i] = (w1*mim1 + w2*mi)/(w1+w2);
if (isnan(impl_->dydx_[i])) {
impl_->dydx_[i] = 0;
}
Real mnm4 = (impl_->y_[n-3]-impl_->y_[n-4])/(impl_->x_[n-3]-impl_->x_[n-4]);
Real mnm3 = (impl_->y_[n-2]-impl_->y_[n-3])/(impl_->x_[n-2]-impl_->x_[n-3]);
Real mnm2 = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1]-impl_->x_[n-2]);
Real mnm1 = 2*mnm2 - mnm3;
Real mn = 2*mnm1 - mnm2;
w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
impl_->dydx_[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
if (isnan(impl_->dydx_[n-2])) {
impl_->dydx_[n-2] = 0;
}
w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
impl_->dydx_[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
if (isnan(impl_->dydx_[n-1])) {
impl_->dydx_[n-1] = 0;
}
}
private:
std::shared_ptr<detail::makima_detail<RandomAccessContainer>> impl_;
std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
};
}

82
include/boost/math/interpolators/pchip.hpp
View File

@@ -7,7 +7,7 @@
#ifndef BOOST_MATH_INTERPOLATORS_PCHIP_HPP
#define BOOST_MATH_INTERPOLATORS_PCHIP_HPP
#include <memory>
#include <boost/math/interpolators/detail/pchip_detail.hpp>
#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
namespace boost::math::interpolators {
@@ -18,8 +18,54 @@ public:
pchip(RandomAccessContainer && x, RandomAccessContainer && y,
Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()) : impl_(std::make_shared<detail::pchip_detail<RandomAccessContainer>>(std::move(x), std::move(y), left_endpoint_derivative, right_endpoint_derivative))
{}
Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
{
if (x.size() < 4)
{
throw std::domain_error("Must be at least four data points.");
}
RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());
if (isnan(left_endpoint_derivative))
{
// O(h) finite difference derivative:
// This, I believe, is the only derivative guaranteed to be monotonic:
s[0] = (y[1]-y[0])/(x[1]-x[0]);
}
else
{
s[0] = left_endpoint_derivative;
}
for (decltype(s.size()) k = 1; k < s.size()-1; ++k) {
Real hkm1 = x[k] - x[k-1];
Real dkm1 = (y[k] - y[k-1])/hkm1;
Real hk = x[k+1] - x[k];
Real dk = (y[k+1] - y[k])/hk;
Real w1 = 2*hk + hkm1;
Real w2 = hk + 2*hkm1;
if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
{
s[k] = 0;
}
else
{
s[k] = (w1+w2)/(w1/dkm1 + w2/dk);
}
}
// Quadratic extrapolation at the other end:
auto n = s.size();
if (isnan(right_endpoint_derivative))
{
s[n-1] = (y[n-1]-y[n-2])/(x[n-1] - x[n-2]);
}
else
{
s[n-1] = right_endpoint_derivative;
}
impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));
}
Real operator()(Real x) const {
return impl_->operator()(x);
@@ -36,11 +82,37 @@ public:
}
void push_back(Real x, Real y) {
impl_->push_back(x, y);
using std::abs;
using std::isnan;
if (x <= impl_->x_.back()) {
throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
}
impl_->x_.push_back(x);
impl_->y_.push_back(y);
impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());
auto n = impl_->size();
impl_->dydx_[n-1] = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1] - impl_->x_[n-2]);
// Now fix s_[n-2]:
auto k = n-2;
Real hkm1 = impl_->x_[k] - impl_->x_[k-1];
Real dkm1 = (impl_->y_[k] - impl_->y_[k-1])/hkm1;
Real hk = impl_->x_[k+1] - impl_->x_[k];
Real dk = (impl_->y_[k+1] - impl_->y_[k])/hk;
Real w1 = 2*hk + hkm1;
Real w2 = hk + 2*hkm1;
if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
{
impl_->dydx_[k] = 0;
}
else
{
impl_->dydx_[k] = (w1+w2)/(w1/dkm1 + w2/dk);
}
}
private:
std::shared_ptr<detail::pchip_detail<RandomAccessContainer>> impl_;
std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
};
}