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