ceph/src/boost/boost/math/interpolators/pchip.hpp

   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