PCHIP 犰狳函数

Posted 2023-02-16

【中文标题】PCHIP 犰狳函数

【英文标题】：PCHIP armadillo function

【发布时间】：2020-07-27 22:30:36

【问题描述】：

我正在尝试使用犰狳将一些 MATLAB 代码转换为 C++，其中使用的 MATLAB 函数之一是 interp1。 “很容易”，我虽然，犰狳有 interp1。 “线性应该足够好”，我错了。所以我搜索了 interp1.m 源代码并找到了 Octave 源代码，它使用 pchip.m 我找到了 octave 源代码，但 pchip.m 使用 pchip_deriv.cc ，然后似乎使用了 fortran fcn。

那么在我开始真正深入研究 pchip 的转换之前，是否还有其他包含 pchip 的库或资源可供我使用？

【问题讨论】：

C++ 中的数值配方，作者：Press、Teukolsky、Vetterling 和 Flannery...？

您可以手动实现它，它足够短以适合 lambda。

【参考方案1】：

如果其他人需要这个，我查看here 来找到方程式。对于大型数据集，它似乎无法正常工作，因此使用矩阵公式，我创建了一个滑动窗口 PCHIP interp。仍然不匹配 MATLAB，所以我使用线性插值进行平均，并且非常接近。这是正确的方法吗？可能不会……有用吗？是的！

void pchip_slide(const arma::vec& x, const arma::vec& y, const arma::vec& x_new, arma::vec& y_custom,const int M=4)
vec y_interp=zeros<vec>(size(x_new));
interp1(x,y,x_new,y_interp);
int I=0;
int start_interp=0;
int end_interp=0;
for(int ii=0;ii<x_new.n_elem;ii++)
    
    I=index_min(abs(x-x_new(ii)));
    start_interp=std::max((I-2),0);
    end_interp=std::min((start_interp+M-1),int(numel(x)-1));
    
    vec x_mini=x(span(start_interp,end_interp));
    vec y_mini=y(span(start_interp,end_interp));
    mat x_mat=ones<mat>(x_mini.n_elem,x_mini.n_elem);
    for(int ll=0;ll<x_mini.n_elem-1;ll++)
        x_mat.row(ll)=pow((x_mini.t()),(x_mini.n_elem-ll));
    
    vec c_mini=solve( x_mat, y_mini,solve_opts::fast + solve_opts::no_approx);
    rowvec x_pchip_mini=ones<rowvec>(x_mini.n_elem);
    for (int ll=0;ll<x_mini.n_elem-1;ll++)
        x_pchip_mini(ll)=pow((x_new(ii)),(x_mini.n_elem-ll));
    
    y_custom(ii)=conv_to<double>::from(x_pchip_mini*c_mini);
    if ((x_new(ii)>=x(0))&& (x_new(ii)<=x(x.n_elem-1)))
        y_custom(ii)=(y_custom(ii)*1/M+y_interp(ii)*3/M);
    

return;