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