C++标准::向量的高效插值

C++ Efficient interpolation of a std::vector

本文关键字：高效 插值 向量 标准 C++      更新时间：2024-09-27

我需要通过插值找到未知函数的值。问题是我创建的这个太低效了。

首先，我读取了一个数据文件，其中包含y=g(T(和T，但形式是离散的。我将它们的值存储在std::vector<double>中。

在此之后，我将T(std::vector<double> Tgdat(转换为x(std::vector<double> xgdat(。这将是伴随y轴的x轴(std::vector<double> gdat(。

然后，我创建一个函数来插值我的向量std::vector<double> gdat，这样，给定某个x(它的值在向量std::vector<double> xgdat的两个元素之间(，程序可以为g(x(吐出一些值。这个函数通过引用接收矢量，不是因为我想修改它们(这就是为什么我也把它们作为const传递(，而是为了让计算机不必创建它的副本

double geffx (double x, const std::vector<double> &gdat, const std::vector<double> &xgdat)
{
//Local variables
double g;
int k,l;
//Find the index of the element of xgdat that is nearest to x
auto i = min_element(xgdat.begin(), xgdat.end(),
[x] (double a, double b)
{
return abs(x-a)<abs(x-b);
});
k = std::distance(xgdat.begin(), i); //Nearest index
//Find the index of the element of xgdat that is nearest to x
//and it is not the same index as before
auto j = min_element(xgdat.begin(), xgdat.end(),
[x,&xgdat,k] (double a, double b)
{
if (a!=xgdat[k]) return abs(x-a)<abs(x-b);
else return false;
});
l = std::distance(xgdat.begin(), j); //Second nearest index
//Interpolation:
if(xgdat[k]<xgdat[l]) 
g = gdat[k]+(x-xgdat[k])*(gdat[l]-gdat[k])/(xgdat[l]-xgdat[k]);
else 
g = gdat[l]+(x-xgdat[l])*(gdat[k]-gdat[l])/(xgdat[k]-xgdat[l]);
return g;
}

这似乎效率很低，但我无法用更有效的方式来解决同样的问题。我尝试过const的方法，也通过引用传递，但我想最大的问题是min_element()函数/方法，也许它也与最后返回g值的if-else有关。

编辑：额外信息

我使用g++作为编译器，元素的数量是275。

由于该函数是EDO求解器的一部分，因此在每一步(1e4步(中都会多次调用它，直到它收敛。我需要4个插值器，每个插值器被调用不止一次进行评估，所以我认为函数需要被访问超过1e6次。

当我用一个常数(不需要插值器(代替g(x(时，执行时间大约是1-10秒。现在是45分钟-1小时。(非常糟糕，我知道，这就是我需要帮助的原因(

失去min_element，失去绝对距离比较。它是一个凸变换，可以有效地搜索，但不能被C++标准库中存在的任何函数搜索。

无论如何，你不想要最接近的两个点，你想把你的评价放在上面和下面。("最接近的两个"总是给出括号对的唯一情况是当样本间隔均匀时，如果这是真的，你根本不需要搜索，你可以直接使用间隔来计算索引(。

使用lower_bound在排序数组中对您关心的x进行有效的二进制搜索。这是你括号的一边，另一个指数低一个。

最后，代码的顶部看起来像：

//Find the index of the element of xgdat that is nearest-above to x
auto i = lower_bound(xgdat.begin(), xgdat.end(), x); 
//If the vector values are in decreasing order use:
//auto i = lower_bound(xgdat.rbegin(), xgdat.rend(), x);
k = i - xgdat.begin(); //Nearest index
if (i == xgdat.end())
--k;  // extrapolating above
else if (*i == x)
return gdat[k];
l = k? k - 1: 1; //nearest-below index, except when extrapolating downward
// proceed with linear interpolation/extrapolation using l and k

关于代码的一些提示：

在需要的地方声明变量，而不是全部在顶部。

if(xgdat[k]<xgdat[l]) 
g = gdat[k]+(x-xgdat[k])*(gdat[l]-gdat[k])/(xgdat[l]-xgdat[k]);
else 
g = gdat[l]+(x-xgdat[l])*(gdat[k]-gdat[l])/(xgdat[k]-xgdat[l]);

除了交换k和l之外，这两行看起来是一样的。所以不要重复这行：只交换k和l！

if(xgdat[k]>=xgdat[l]) std::swap(k,l);

CCD_ 14只在这里使用？为什么要在函数的顶部声明它？现在根本不需要它：

return gdat[k]+(x-xgdat[k])*(gdat[l]-gdat[k])/(xgdat[l]-xgdat[k]);

你会遇到一些麻烦得到"；索引位置"；调用使用迭代器的标准算法后退出。您应该只使用迭代器。

但这些都不是整体效率。

相关文章：