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