From 2a37abd93a985e90b5cc3afe82ec611dad1c0f02 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sat, 18 Jan 2020 14:11:31 -0500 Subject: [PATCH] Cubic Hermite spline: Backend pchip and makima to cubic_hermite. --- doc/interpolators/makima.qbk | 47 +--- doc/interpolators/pchip.qbk | 11 +- .../detail/cubic_hermite_detail.hpp | 5 +- .../interpolators/detail/makima_detail.hpp | 248 ------------------ .../interpolators/detail/pchip_detail.hpp | 199 -------------- include/boost/math/interpolators/makima.hpp | 134 +++++++++- include/boost/math/interpolators/pchip.hpp | 82 +++++- 7 files changed, 216 insertions(+), 510 deletions(-) delete mode 100644 include/boost/math/interpolators/detail/makima_detail.hpp delete mode 100644 include/boost/math/interpolators/detail/pchip_detail.hpp diff --git a/doc/interpolators/makima.qbk b/doc/interpolators/makima.qbk index f0023f8c9..0075cb854 100644 --- a/doc/interpolators/makima.qbk +++ b/doc/interpolators/makima.qbk @@ -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 -`` + + #include 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/8 11.3 ns -BMMakima/16 11.5 ns -BMMakima/32 12.7 ns -BMMakima/64 13.6 ns -BMMakima/128 14.5 ns -BMMakima/256 15.1 ns -BMMakima/512 17.2 ns -BMMakima/1024 18.3 ns -BMMakima/2048 19.6 ns -BMMakima/4096 20.7 ns -BMMakima/8192 22.1 ns -BMMakima/16384 23.1 ns -BMMakima/32768 24.2 ns -BMMakima/65536 25.3 ns -BMMakima/131072 27.1 ns -BMMakima/262144 28.3 ns -BMMakima/524288 29.9 ns -BMMakima/1048576 31.6 ns -BMMakima/2097152 31.8 ns -BMMakima/4194304 33.7 ns -BMMakima/8388608 35.0 ns -BMMakima/16777216 40.0 ns -BMMakima_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] diff --git a/doc/interpolators/pchip.qbk b/doc/interpolators/pchip.qbk index 38938d328..34d47b6fc 100644 --- a/doc/interpolators/pchip.qbk +++ b/doc/interpolators/pchip.qbk @@ -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 -`` + + #include 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] diff --git a/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp b/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp index 595d61297..f7a3b60ae 100644 --- a/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp +++ b/include/boost/math/interpolators/detail/cubic_hermite_detail.hpp @@ -130,7 +130,10 @@ public: return os; } -private: + auto size() const { + return x_.size(); + } + RandomAccessContainer x_; RandomAccessContainer y_; RandomAccessContainer dydx_; diff --git a/include/boost/math/interpolators/detail/makima_detail.hpp b/include/boost/math/interpolators/detail/makima_detail.hpp deleted file mode 100644 index 38dcea747..000000000 --- a/include/boost/math/interpolators/detail/makima_detail.hpp +++ /dev/null @@ -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 -#include -#include -#include -#include -#include - -namespace boost::math::interpolators::detail { - -template -class makima_detail { -public: - using Real = typename RandomAccessContainer::value_type; - - makima_detail(RandomAccessContainer && x, RandomAccessContainer && y, - Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), - Real right_endpoint_derivative = std::numeric_limits::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::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::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::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::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 diff --git a/include/boost/math/interpolators/detail/pchip_detail.hpp b/include/boost/math/interpolators/detail/pchip_detail.hpp deleted file mode 100644 index b1e57816b..000000000 --- a/include/boost/math/interpolators/detail/pchip_detail.hpp +++ /dev/null @@ -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 -#include -#include -#include -#include -#include - -namespace boost::math::interpolators::detail { - -template -class pchip_detail { -public: - using Real = typename RandomAccessContainer::value_type; - - pchip_detail(RandomAccessContainer && x, RandomAccessContainer && y, - Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), - Real right_endpoint_derivative = std::numeric_limits::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::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::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::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::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 diff --git a/include/boost/math/interpolators/makima.hpp b/include/boost/math/interpolators/makima.hpp index c2fc1dad7..4284073ca 100644 --- a/include/boost/math/interpolators/makima.hpp +++ b/include/boost/math/interpolators/makima.hpp @@ -10,7 +10,8 @@ #ifndef BOOST_MATH_INTERPOLATORS_MAKIMA_HPP #define BOOST_MATH_INTERPOLATORS_MAKIMA_HPP #include -#include +#include +#include namespace boost::math::interpolators { @@ -21,8 +22,90 @@ public: makima(RandomAccessContainer && x, RandomAccessContainer && y, Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), - Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()) : impl_(std::make_shared>(std::move(x), std::move(y), left_endpoint_derivative, right_endpoint_derivative)) - {} + Real right_endpoint_derivative = std::numeric_limits::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::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>(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::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> impl_; + std::shared_ptr> impl_; }; } diff --git a/include/boost/math/interpolators/pchip.hpp b/include/boost/math/interpolators/pchip.hpp index 9c22332bc..3a37ddf67 100644 --- a/include/boost/math/interpolators/pchip.hpp +++ b/include/boost/math/interpolators/pchip.hpp @@ -7,7 +7,7 @@ #ifndef BOOST_MATH_INTERPOLATORS_PCHIP_HPP #define BOOST_MATH_INTERPOLATORS_PCHIP_HPP #include -#include +#include namespace boost::math::interpolators { @@ -18,8 +18,54 @@ public: pchip(RandomAccessContainer && x, RandomAccessContainer && y, Real left_endpoint_derivative = std::numeric_limits::quiet_NaN(), - Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()) : impl_(std::make_shared>(std::move(x), std::move(y), left_endpoint_derivative, right_endpoint_derivative)) - {} + Real right_endpoint_derivative = std::numeric_limits::quiet_NaN()) + { + if (x.size() < 4) + { + throw std::domain_error("Must be at least four data points."); + } + RandomAccessContainer s(x.size(), std::numeric_limits::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>(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::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> impl_; + std::shared_ptr> impl_; }; }