math/interpolators/pchip.hpp

0001 // Copyright Nick Thompson, 2020
0002 // Use, modification and distribution are subject to the
0003 // Boost Software License, Version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_INTERPOLATORS_PCHIP_HPP
0008 #define BOOST_MATH_INTERPOLATORS_PCHIP_HPP
0009 #include <sstream>
0010 #include <memory>
0011 #include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
0012
0013 namespace boost {
0014 namespace math {
0015 namespace interpolators {
0016
0017 template<class RandomAccessContainer>
0018 class pchip {
0019 public:
0020     using Real = typename RandomAccessContainer::value_type;
0021
0022     pchip(RandomAccessContainer && x, RandomAccessContainer && y,
0023           Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
0024           Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
0025     {
0026         using std::isnan;
0027         if (x.size() < 4)
0028         {
0029             std::ostringstream oss;
0030             oss << __FILE__ << ":" << __LINE__ << ":" << __func__;
0031             oss << " This interpolator requires at least four data points.";
0032             throw std::domain_error(oss.str());
0033         }
0034         RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());
0035         if (isnan(left_endpoint_derivative))
0036         {
0037             // If the derivative is not specified, this seems as good a choice as any.
0038             // In particular, it satisfies the monotonicity constraint 0 <= |y'[0]| < 4Delta_i,
0039             // where Delta_i is the secant slope:
0040             s[0] = (y[1]-y[0])/(x[1]-x[0]);
0041         }
0042         else
0043         {
0044             s[0] = left_endpoint_derivative;
0045         }
0046
0047         for (decltype(s.size()) k = 1; k < s.size()-1; ++k) {
0048             Real hkm1 = x[k] - x[k-1];
0049             Real dkm1 = (y[k] - y[k-1])/hkm1;
0050
0051             Real hk = x[k+1] - x[k];
0052             Real dk = (y[k+1] - y[k])/hk;
0053             Real w1 = 2*hk + hkm1;
0054             Real w2 = hk + 2*hkm1;
0055             if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
0056             {
0057                 s[k] = 0;
0058             }
0059             else
0060             {
0061                 // See here:
0062                 // https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/interp.pdf
0063                 // Un-numbered equation just before Section 3.5:
0064                 s[k] = (w1+w2)/(w1/dkm1 + w2/dk);
0065             }
0066
0067         }
0068         auto n = s.size();
0069         if (isnan(right_endpoint_derivative))
0070         {
0071             s[n-1] = (y[n-1]-y[n-2])/(x[n-1] - x[n-2]);
0072         }
0073         else
0074         {
0075             s[n-1] = right_endpoint_derivative;
0076         }
0077         impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));
0078     }
0079
0080     Real operator()(Real x) const {
0081         return impl_->operator()(x);
0082     }
0083
0084     Real prime(Real x) const {
0085         return impl_->prime(x);
0086     }
0087
0088     friend std::ostream& operator<<(std::ostream & os, const pchip & m)
0089     {
0090         os << *m.impl_;
0091         return os;
0092     }
0093
0094     void push_back(Real x, Real y) {
0095         using std::abs;
0096         using std::isnan;
0097         if (x <= impl_->x_.back()) {
0098              throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
0099         }
0100         impl_->x_.push_back(x);
0101         impl_->y_.push_back(y);
0102         impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());
0103         auto n = impl_->size();
0104         impl_->dydx_[n-1] = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1] - impl_->x_[n-2]);
0105         // Now fix s_[n-2]:
0106         auto k = n-2;
0107         Real hkm1 = impl_->x_[k] - impl_->x_[k-1];
0108         Real dkm1 = (impl_->y_[k] - impl_->y_[k-1])/hkm1;
0109
0110         Real hk = impl_->x_[k+1] - impl_->x_[k];
0111         Real dk = (impl_->y_[k+1] - impl_->y_[k])/hk;
0112         Real w1 = 2*hk + hkm1;
0113         Real w2 = hk + 2*hkm1;
0114         if ( (dk > 0 && dkm1 < 0) || (dk < 0 && dkm1 > 0) || dk == 0 || dkm1 == 0)
0115         {
0116             impl_->dydx_[k] = 0;
0117         }
0118         else
0119         {
0120             impl_->dydx_[k] = (w1+w2)/(w1/dkm1 + w2/dk);
0121         }
0122     }
0123
0124 private:
0125     std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
0126 };
0127
0128 }
0129 }
0130 }
0131 #endif