132 lines · plain
1// Copyright Nick Thompson, 20202// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0.4// (See accompanying file LICENSE_1_0.txt5// or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_INTERPOLATORS_PCHIP_HPP8#define BOOST_MATH_INTERPOLATORS_PCHIP_HPP9#include <sstream>10#include <memory>11#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>12 13namespace boost {14namespace math {15namespace interpolators {16 17template<class RandomAccessContainer>18class pchip {19public:20 using Real = typename RandomAccessContainer::value_type;21 22 pchip(RandomAccessContainer && x, RandomAccessContainer && y,23 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),24 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())25 {26 using std::isnan;27 if (x.size() < 4)28 {29 std::ostringstream oss;30 oss << __FILE__ << ":" << __LINE__ << ":" << __func__;31 oss << " This interpolator requires at least four data points.";32 throw std::domain_error(oss.str());33 }34 RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());35 if (isnan(left_endpoint_derivative))36 {37 // If the derivative is not specified, this seems as good a choice as any.38 // In particular, it satisfies the monotonicity constraint 0 <= |y'[0]| < 4Delta_i,39 // where Delta_i is the secant slope:40 s[0] = (y[1]-y[0])/(x[1]-x[0]);41 }42 else43 {44 s[0] = left_endpoint_derivative;45 }46 47 for (decltype(s.size()) k = 1; k < s.size()-1; ++k) {48 Real hkm1 = x[k] - x[k-1];49 Real dkm1 = (y[k] - y[k-1])/hkm1;50 51 Real hk = x[k+1] - x[k];52 Real dk = (y[k+1] - y[k])/hk;53 Real w1 = 2*hk + hkm1;54 Real w2 = hk + 2*hkm1;55 if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)56 {57 s[k] = 0;58 }59 else60 {61 // See here:62 // https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/interp.pdf63 // Un-numbered equation just before Section 3.5:64 s[k] = (w1+w2)/(w1/dkm1 + w2/dk);65 }66 67 }68 auto n = s.size();69 if (isnan(right_endpoint_derivative))70 {71 s[n-1] = (y[n-1]-y[n-2])/(x[n-1] - x[n-2]);72 }73 else74 {75 s[n-1] = right_endpoint_derivative;76 }77 impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));78 }79 80 Real operator()(Real x) const {81 return impl_->operator()(x);82 }83 84 Real prime(Real x) const {85 return impl_->prime(x);86 }87 88 friend std::ostream& operator<<(std::ostream & os, const pchip & m)89 {90 os << *m.impl_;91 return os;92 }93 94 void push_back(Real x, Real y) {95 using std::abs;96 using std::isnan;97 if (x <= impl_->x_.back()) {98 throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");99 }100 impl_->x_.push_back(x);101 impl_->y_.push_back(y);102 impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());103 auto n = impl_->size();104 impl_->dydx_[n-1] = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1] - impl_->x_[n-2]);105 // Now fix s_[n-2]:106 auto k = n-2;107 Real hkm1 = impl_->x_[k] - impl_->x_[k-1];108 Real dkm1 = (impl_->y_[k] - impl_->y_[k-1])/hkm1;109 110 Real hk = impl_->x_[k+1] - impl_->x_[k];111 Real dk = (impl_->y_[k+1] - impl_->y_[k])/hk;112 Real w1 = 2*hk + hkm1;113 Real w2 = hk + 2*hkm1;114 if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)115 {116 impl_->dydx_[k] = 0;117 }118 else119 {120 impl_->dydx_[k] = (w1+w2)/(w1/dkm1 + w2/dk);121 }122 }123 124private:125 std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;126};127 128}129}130}131#endif132