1 // Copyright Nick Thompson, 2020 2 // Use, modification and distribution are subject to the 3 // Boost Software License, Version 1.0. 4 // (See accompanying file LICENSE_1_0.txt 5 // or copy at http://www.boost.org/LICENSE_1_0.txt) 6 7 #ifndef BOOST_MATH_INTERPOLATORS_PCHIP_HPP 8 #define BOOST_MATH_INTERPOLATORS_PCHIP_HPP 9 #include <memory> 10 #include <boost/math/interpolators/detail/cubic_hermite_detail.hpp> 11 12 namespace boost { 13 namespace math { 14 namespace interpolators { 15 16 template<class RandomAccessContainer> 17 class pchip { 18 public: 19 using Real = typename RandomAccessContainer::value_type; 20 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 ())21 pchip(RandomAccessContainer && x, RandomAccessContainer && y, 22 Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(), 23 Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN()) 24 { 25 using std::isnan; 26 if (x.size() < 4) 27 { 28 throw std::domain_error("Must be at least four data points."); 29 } 30 RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN()); 31 if (isnan(left_endpoint_derivative)) 32 { 33 // O(h) finite difference derivative: 34 // This, I believe, is the only derivative guaranteed to be monotonic: 35 s[0] = (y[1]-y[0])/(x[1]-x[0]); 36 } 37 else 38 { 39 s[0] = left_endpoint_derivative; 40 } 41 42 for (decltype(s.size()) k = 1; k < s.size()-1; ++k) { 43 Real hkm1 = x[k] - x[k-1]; 44 Real dkm1 = (y[k] - y[k-1])/hkm1; 45 46 Real hk = x[k+1] - x[k]; 47 Real dk = (y[k+1] - y[k])/hk; 48 Real w1 = 2*hk + hkm1; 49 Real w2 = hk + 2*hkm1; 50 if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0) 51 { 52 s[k] = 0; 53 } 54 else 55 { 56 s[k] = (w1+w2)/(w1/dkm1 + w2/dk); 57 } 58 59 } 60 // Quadratic extrapolation at the other end: 61 auto n = s.size(); 62 if (isnan(right_endpoint_derivative)) 63 { 64 s[n-1] = (y[n-1]-y[n-2])/(x[n-1] - x[n-2]); 65 } 66 else 67 { 68 s[n-1] = right_endpoint_derivative; 69 } 70 impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s)); 71 } 72 operator ()(Real x) const73 Real operator()(Real x) const { 74 return impl_->operator()(x); 75 } 76 prime(Real x) const77 Real prime(Real x) const { 78 return impl_->prime(x); 79 } 80 operator <<(std::ostream & os,const pchip & m)81 friend std::ostream& operator<<(std::ostream & os, const pchip & m) 82 { 83 os << *m.impl_; 84 return os; 85 } 86 push_back(Real x,Real y)87 void push_back(Real x, Real y) { 88 using std::abs; 89 using std::isnan; 90 if (x <= impl_->x_.back()) { 91 throw std::domain_error("Calling push_back must preserve the monotonicity of the x's"); 92 } 93 impl_->x_.push_back(x); 94 impl_->y_.push_back(y); 95 impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN()); 96 auto n = impl_->size(); 97 impl_->dydx_[n-1] = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1] - impl_->x_[n-2]); 98 // Now fix s_[n-2]: 99 auto k = n-2; 100 Real hkm1 = impl_->x_[k] - impl_->x_[k-1]; 101 Real dkm1 = (impl_->y_[k] - impl_->y_[k-1])/hkm1; 102 103 Real hk = impl_->x_[k+1] - impl_->x_[k]; 104 Real dk = (impl_->y_[k+1] - impl_->y_[k])/hk; 105 Real w1 = 2*hk + hkm1; 106 Real w2 = hk + 2*hkm1; 107 if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0) 108 { 109 impl_->dydx_[k] = 0; 110 } 111 else 112 { 113 impl_->dydx_[k] = (w1+w2)/(w1/dkm1 + w2/dk); 114 } 115 } 116 117 private: 118 std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_; 119 }; 120 121 } 122 } 123 } 124 #endif 125