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