179 lines · plain
1// Copyright Nick Thompson, 20202// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0.4// (See accompanying file LICENSE_1_0.txt5// 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.3216099 10#ifndef BOOST_MATH_INTERPOLATORS_MAKIMA_HPP11#define BOOST_MATH_INTERPOLATORS_MAKIMA_HPP12#include <memory>13#include <cmath>14#include <boost/math/interpolators/detail/cubic_hermite_detail.hpp>15 16namespace boost {17namespace math {18namespace interpolators {19 20template<class RandomAccessContainer>21class makima {22public: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 else54 {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 else105 {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 171private:172 std::shared_ptr<detail::cubic_hermite_detail<RandomAccessContainer>> impl_;173};174 175}176}177}178#endif179