std::vector 的 C++ 高效插值

Question

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

首先，我读取了一个包含 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() function/method，也许它也与 if-else 在结束到 return g 的值。

编辑：额外信息

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

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

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

Answer 1

关于代码的一些提示：

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

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);

和g只用在这里？你为什么在函数的顶部声明它？现在完全不需要了：

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

在调用与迭代器一起使用的标准算法后，要恢复“索引位置”会遇到一些麻烦。你应该只使用迭代器。

None 不过这是整体效率。

Answer 2

输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 = xgdat.begin() - i; //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

std::vector 的 C++ 高效插值

C++ Efficient interpolation of a std::vector

c++

sorting

interpolation

std

stdvector