]>
Commit | Line | Data |
---|---|---|
f67539c2 TL |
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 | ||
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 | ||
70 | Real operator()(Real x) const { | |
71 | return impl_->operator()(x); | |
72 | } | |
73 | ||
74 | Real prime(Real x) const { | |
75 | return impl_->prime(x); | |
76 | } | |
77 | ||
78 | friend std::ostream& operator<<(std::ostream & os, const pchip & m) | |
79 | { | |
80 | os << *m.impl_; | |
81 | return os; | |
82 | } | |
83 | ||
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 |