ceph/src/boost/boost/math/interpolators/makima.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 // See: https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/
   8 // And: https://doi.org/10.1145/321607.321609
   9
  10 #ifndef BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
  11 #define BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
  12 #include <memory>
  13 #include <cmath>
  14 #include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
  15
  16 namespace boost {
  17 namespace math {
  18 namespace interpolators {
  19
  20 template<class RandomAccessContainer>
  21 class makima {
  22 public:
  23     using Real = typename RandomAccessContainer::value_type;
  24
  25     makima(RandomAccessContainer && x, RandomAccessContainer && y,
  26            Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
  27            Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
  28     {
  29         using std::isnan;
  30         using std::abs;
  31         if (x.size() < 4)
  32         {
  33             throw std::domain_error("Must be at least four data points.");
  34         }
  35         RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());
  36         Real m2 = (y[3]-y[2])/(x[3]-x[2]);
  37         Real m1 = (y[2]-y[1])/(x[2]-x[1]);
  38         Real m0 = (y[1]-y[0])/(x[1]-x[0]);
  39         // Quadratic extrapolation: m_{-1} = 2m_0 - m_1:
  40         Real mm1 = 2*m0 - m1;
  41         // Quadratic extrapolation: m_{-2} = 2*m_{-1}-m_0:
  42         Real mm2 = 2*mm1 - m0;
  43         Real w1 = abs(m1-m0) + abs(m1+m0)/2;
  44         Real w2 = abs(mm1-mm2) + abs(mm1+mm2)/2;
  45         if (isnan(left_endpoint_derivative))
  46         {
  47             s[0] = (w1*mm1 + w2*m0)/(w1+w2);
  48             if (isnan(s[0]))
  49             {
  50                 s[0] = 0;
  51             }
  52         }
  53         else
  54         {
  55             s[0] = left_endpoint_derivative;
  56         }
  57
  58         w1 = abs(m2-m1) + abs(m2+m1)/2;
  59         w2 = abs(m0-mm1) + abs(m0+mm1)/2;
  60         s[1] = (w1*m0 + w2*m1)/(w1+w2);
  61         if (isnan(s[1])) {
  62             s[1] = 0;
  63         }
  64
  65         for (decltype(s.size()) i = 2; i < s.size()-2; ++i) {
  66             Real mim2 = (y[i-1]-y[i-2])/(x[i-1]-x[i-2]);
  67             Real mim1 = (y[i  ]-y[i-1])/(x[i  ]-x[i-1]);
  68             Real mi   = (y[i+1]-y[i  ])/(x[i+1]-x[i  ]);
  69             Real mip1 = (y[i+2]-y[i+1])/(x[i+2]-x[i+1]);
  70             w1 = abs(mip1-mi) + abs(mip1+mi)/2;
  71             w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
  72             s[i] = (w1*mim1 + w2*mi)/(w1+w2);
  73             if (isnan(s[i])) {
  74                 s[i] = 0;
  75             }
  76         }
  77         // Quadratic extrapolation at the other end:
  78
  79         decltype(s.size()) n = s.size();
  80         Real mnm4 = (y[n-3]-y[n-4])/(x[n-3]-x[n-4]);
  81         Real mnm3 = (y[n-2]-y[n-3])/(x[n-2]-x[n-3]);
  82         Real mnm2 = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
  83         Real mnm1 = 2*mnm2 - mnm3;
  84         Real mn = 2*mnm1 - mnm2;
  85         w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
  86         w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
  87
  88         s[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
  89         if (isnan(s[n-2])) {
  90             s[n-2] = 0;
  91         }
  92
  93         w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
  94         w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
  95
  96
  97         if (isnan(right_endpoint_derivative))
  98         {
  99             s[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
 100             if (isnan(s[n-1])) {
 101                 s[n-1] = 0;
 102             }
 103         }
 104         else
 105         {
 106             s[n-1] = right_endpoint_derivative;
 107         }
 108
 109         impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));
 110     }
 111
 112     Real operator()(Real x) const {
 113         return impl_->operator()(x);
 114     }
 115
 116     Real prime(Real x) const {
 117         return impl_->prime(x);
 118     }
 119
 120     friend std::ostream& operator<<(std::ostream & os, const makima & m)
 121     {
 122         os << *m.impl_;
 123         return os;
 124     }
 125
 126     void push_back(Real x, Real y) {
 127         using std::abs;
 128         using std::isnan;
 129         if (x <= impl_->x_.back()) {
 130              throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
 131         }
 132         impl_->x_.push_back(x);
 133         impl_->y_.push_back(y);
 134         impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());
 135         // dydx_[n-2] was computed by extrapolation. Now dydx_[n-2] -> dydx_[n-3], and it can be computed by the same formula.
 136         decltype(impl_->size()) n = impl_->size();
 137         auto i = n - 3;
 138         Real mim2 = (impl_->y_[i-1]-impl_->y_[i-2])/(impl_->x_[i-1]-impl_->x_[i-2]);
 139         Real mim1 = (impl_->y_[i  ]-impl_->y_[i-1])/(impl_->x_[i  ]-impl_->x_[i-1]);
 140         Real mi   = (impl_->y_[i+1]-impl_->y_[i  ])/(impl_->x_[i+1]-impl_->x_[i  ]);
 141         Real mip1 = (impl_->y_[i+2]-impl_->y_[i+1])/(impl_->x_[i+2]-impl_->x_[i+1]);
 142         Real w1 = abs(mip1-mi) + abs(mip1+mi)/2;
 143         Real w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
 144         impl_->dydx_[i] = (w1*mim1 + w2*mi)/(w1+w2);
 145         if (isnan(impl_->dydx_[i])) {
 146             impl_->dydx_[i] = 0;
 147         }
 148
 149         Real mnm4 = (impl_->y_[n-3]-impl_->y_[n-4])/(impl_->x_[n-3]-impl_->x_[n-4]);
 150         Real mnm3 = (impl_->y_[n-2]-impl_->y_[n-3])/(impl_->x_[n-2]-impl_->x_[n-3]);
 151         Real mnm2 = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1]-impl_->x_[n-2]);
 152         Real mnm1 = 2*mnm2 - mnm3;
 153         Real mn = 2*mnm1 - mnm2;
 154         w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
 155         w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
 156
 157         impl_->dydx_[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
 158         if (isnan(impl_->dydx_[n-2])) {
 159             impl_->dydx_[n-2] = 0;
 160         }
 161
 162         w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
 163         w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
 164
 165         impl_->dydx_[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
 166         if (isnan(impl_->dydx_[n-1])) {
 167             impl_->dydx_[n-1] = 0;
 168         }
 169     }
 170
 171 private:
 172     std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
 173 };
 174
 175 }
 176 }
 177 }
 178 #endif