math/interpolators/makima.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 // See: https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/
0008 // And: https://doi.org/10.1145/321607.321609
0009
0010 #ifndef BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
0011 #define BOOST_MATH_INTERPOLATORS_MAKIMA_HPP
0012 #include <memory>
0013 #include <cmath>
0014 #include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>
0015
0016 namespace boost {
0017 namespace math {
0018 namespace interpolators {
0019
0020 template<class RandomAccessContainer>
0021 class makima {
0022 public:
0023     using Real = typename RandomAccessContainer::value_type;
0024
0025     makima(RandomAccessContainer && x, RandomAccessContainer && y,
0026            Real left_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN(),
0027            Real right_endpoint_derivative = std::numeric_limits<Real>::quiet_NaN())
0028     {
0029         using std::isnan;
0030         using std::abs;
0031         if (x.size() < 4)
0032         {
0033             throw std::domain_error("Must be at least four data points.");
0034         }
0035         RandomAccessContainer s(x.size(), std::numeric_limits<Real>::quiet_NaN());
0036         Real m2 = (y[3]-y[2])/(x[3]-x[2]);
0037         Real m1 = (y[2]-y[1])/(x[2]-x[1]);
0038         Real m0 = (y[1]-y[0])/(x[1]-x[0]);
0039         // Quadratic extrapolation: m_{-1} = 2m_0 - m_1:
0040         Real mm1 = 2*m0 - m1;
0041         // Quadratic extrapolation: m_{-2} = 2*m_{-1}-m_0:
0042         Real mm2 = 2*mm1 - m0;
0043         Real w1 = abs(m1-m0) + abs(m1+m0)/2;
0044         Real w2 = abs(mm1-mm2) + abs(mm1+mm2)/2;
0045         if (isnan(left_endpoint_derivative))
0046         {
0047             s[0] = (w1*mm1 + w2*m0)/(w1+w2);
0048             if (isnan(s[0]))
0049             {
0050                 s[0] = 0;
0051             }
0052         }
0053         else
0054         {
0055             s[0] = left_endpoint_derivative;
0056         }
0057
0058         w1 = abs(m2-m1) + abs(m2+m1)/2;
0059         w2 = abs(m0-mm1) + abs(m0+mm1)/2;
0060         s[1] = (w1*m0 + w2*m1)/(w1+w2);
0061         if (isnan(s[1])) {
0062             s[1] = 0;
0063         }
0064
0065         for (decltype(s.size()) i = 2; i < s.size()-2; ++i) {
0066             Real mim2 = (y[i-1]-y[i-2])/(x[i-1]-x[i-2]);
0067             Real mim1 = (y[i  ]-y[i-1])/(x[i  ]-x[i-1]);
0068             Real mi   = (y[i+1]-y[i  ])/(x[i+1]-x[i  ]);
0069             Real mip1 = (y[i+2]-y[i+1])/(x[i+2]-x[i+1]);
0070             w1 = abs(mip1-mi) + abs(mip1+mi)/2;
0071             w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
0072             s[i] = (w1*mim1 + w2*mi)/(w1+w2);
0073             if (isnan(s[i])) {
0074                 s[i] = 0;
0075             }
0076         }
0077         // Quadratic extrapolation at the other end:
0078
0079         decltype(s.size()) n = s.size();
0080         Real mnm4 = (y[n-3]-y[n-4])/(x[n-3]-x[n-4]);
0081         Real mnm3 = (y[n-2]-y[n-3])/(x[n-2]-x[n-3]);
0082         Real mnm2 = (y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
0083         Real mnm1 = 2*mnm2 - mnm3;
0084         Real mn = 2*mnm1 - mnm2;
0085         w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
0086         w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
0087
0088         s[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
0089         if (isnan(s[n-2])) {
0090             s[n-2] = 0;
0091         }
0092
0093         w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
0094         w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
0095
0096
0097         if (isnan(right_endpoint_derivative))
0098         {
0099             s[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
0100             if (isnan(s[n-1])) {
0101                 s[n-1] = 0;
0102             }
0103         }
0104         else
0105         {
0106             s[n-1] = right_endpoint_derivative;
0107         }
0108
0109         impl_ = std::make_shared<detail::cubic_hermite_detail<RandomAccessContainer>>(std::move(x), std::move(y), std::move(s));
0110     }
0111
0112     Real operator()(Real x) const {
0113         return impl_->operator()(x);
0114     }
0115
0116     Real prime(Real x) const {
0117         return impl_->prime(x);
0118     }
0119
0120     friend std::ostream& operator<<(std::ostream & os, const makima & m)
0121     {
0122         os << *m.impl_;
0123         return os;
0124     }
0125
0126     void push_back(Real x, Real y) {
0127         using std::abs;
0128         using std::isnan;
0129         if (x <= impl_->x_.back()) {
0130              throw std::domain_error("Calling push_back must preserve the monotonicity of the x's");
0131         }
0132         impl_->x_.push_back(x);
0133         impl_->y_.push_back(y);
0134         impl_->dydx_.push_back(std::numeric_limits<Real>::quiet_NaN());
0135         // dydx_[n-2] was computed by extrapolation. Now dydx_[n-2] -> dydx_[n-3], and it can be computed by the same formula.
0136         decltype(impl_->size()) n = impl_->size();
0137         auto i = n - 3;
0138         Real mim2 = (impl_->y_[i-1]-impl_->y_[i-2])/(impl_->x_[i-1]-impl_->x_[i-2]);
0139         Real mim1 = (impl_->y_[i  ]-impl_->y_[i-1])/(impl_->x_[i  ]-impl_->x_[i-1]);
0140         Real mi   = (impl_->y_[i+1]-impl_->y_[i  ])/(impl_->x_[i+1]-impl_->x_[i  ]);
0141         Real mip1 = (impl_->y_[i+2]-impl_->y_[i+1])/(impl_->x_[i+2]-impl_->x_[i+1]);
0142         Real w1 = abs(mip1-mi) + abs(mip1+mi)/2;
0143         Real w2 = abs(mim1-mim2) + abs(mim1+mim2)/2;
0144         impl_->dydx_[i] = (w1*mim1 + w2*mi)/(w1+w2);
0145         if (isnan(impl_->dydx_[i])) {
0146             impl_->dydx_[i] = 0;
0147         }
0148
0149         Real mnm4 = (impl_->y_[n-3]-impl_->y_[n-4])/(impl_->x_[n-3]-impl_->x_[n-4]);
0150         Real mnm3 = (impl_->y_[n-2]-impl_->y_[n-3])/(impl_->x_[n-2]-impl_->x_[n-3]);
0151         Real mnm2 = (impl_->y_[n-1]-impl_->y_[n-2])/(impl_->x_[n-1]-impl_->x_[n-2]);
0152         Real mnm1 = 2*mnm2 - mnm3;
0153         Real mn = 2*mnm1 - mnm2;
0154         w1 = abs(mnm1 - mnm2) + abs(mnm1+mnm2)/2;
0155         w2 = abs(mnm3 - mnm4) + abs(mnm3+mnm4)/2;
0156
0157         impl_->dydx_[n-2] = (w1*mnm3 + w2*mnm2)/(w1 + w2);
0158         if (isnan(impl_->dydx_[n-2])) {
0159             impl_->dydx_[n-2] = 0;
0160         }
0161
0162         w1 = abs(mn - mnm1) + abs(mn+mnm1)/2;
0163         w2 = abs(mnm2 - mnm3) + abs(mnm2+mnm3)/2;
0164
0165         impl_->dydx_[n-1] = (w1*mnm2 + w2*mnm1)/(w1+w2);
0166         if (isnan(impl_->dydx_[n-1])) {
0167             impl_->dydx_[n-1] = 0;
0168         }
0169     }
0170
0171 private:
0172     std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;
0173 };
0174
0175 }
0176 }
0177 }
0178 #endif