SciPy：在3D网格上插值散乱数据

Question

我的数据是随机分布的 (x, y, z, d)，我需要在 3D 网格上对其进行插值。问题是网格很大：总点数是 nx*ny*nz = 345*188*1501 = 97 354 860.

我尝试使用 Rbf interpolation 同时通过所有点，但出现错误：

numpy.core._exceptions._ArrayMemoryError: Unable to allocate 343. GiB for an array with shape (97354860, 473) and data type float64（这里的473是什么？）

作为我尝试过的例子：

from scipy.interpolate import Rbf
rng = np.random.default_rng()
x, y, z, d = rng.random((4, 50))
rbfi = Rbf(x, y, z, d)

# xi, yi, zi alltogether take about 3 Gb
xi = yi = zi = np.linspace(0, 1, 97354860)
di = rbfi(xi, yi, zi)  # here I get error

但后来我尝试在循环中逐点插值，它似乎工作得足够慢（插值大约需要 20 分钟）：

for i in range(0, xi.shape[0]):
  di[i] = rbfi(xi[i], yi[i], zi[i])

那么问题来了：循环插值和一次插值有什么区别？我认为 Rbf 应该使用循环在每个请求点 (xi, yi, zi) 上内部插入输入数据点。

另外我认为如果我一次通过所有点或者我逐点插值，插值结果应该不会有差异。

Answer 1

正如您提到的，当您使用 NumPy 时，您同时处理所有数据量，这将受到系统内存大小的限制。但是当您使用循环时，由于迭代和直接分配，该工作负载不会应用于 ram。在这方面，我认为使用 chunks 可以非常高效； IE。将某些部分的总大小除以 Chunks。然后，您可以通过在零件上循环并传递 ram 限制来使用数组。因此，建议的方式类似于：

chunk_val = 2000000                  # it is arbitrary and must be choosed based on the system rams size 
chunk = xi.shape[0] // chunk_val
chunk_res = xi.shape[0] % chunk_val

# by array
di = np.array([])
start = 0
for i in range(chunk + 1):
    di = np.append(di, rbfi(xi[start:(i+1) * chunk_val], yi[start:(i+1) * chunk_val], zi[start:(i+1) * chunk_val]))
    start += chunk_val
    if i == chunk:
        di = np.append(di, rbfi(xi[start:xi.shape[0]], yi[start:xi.shape[0]], zi[start:xi.shape[0]]))

# by list
di = []
start = 0
for i in range(chunk + 1):
    di.append(rbfi(xi[start:(i+1) * chunk_val], yi[start:(i+1) * chunk_val], zi[start:(i+1) * chunk_val]))
    start += chunk_val
    if i == chunk:
        di.append(rbfi(xi[start:xi.shape[0]], yi[start:xi.shape[0]], zi[start:xi.shape[0]]))
di = [item for sublist in di for item in sublist]

上述代码中提出了两种方式，一种是数组，一种是列表。我已经在性能方面测试了这种方法，在 Google Collaboratory 上花费了大约 70 s。
如果 Numba 与 rbfi dtype 兼容，它也可以由 Numba 处理。如果是这样，它可能是最快的方法之一，必须对其进行评估。

更新：

我认为唯一的方法是使用块。正如你在rbf source code中看到的那样，rbf中有一系列函数与class相关，它们将是运行并且它们的数据必须存储在系统中的某个位置（rams ) 直到完成所有功能（您可以使用分析器分析 rbf 的每一行，例如 pprofiler 和...以查看源代码的哪一部分消耗的时间最多，ram 或...）。因此，据我所知，当您在迭代中导入数据时，函数会在较小的数据量上运行（可以由系统硬件处理）并且创建的 temporary 数组和...将被释放或在每个循环中被覆盖。因此，迭代不会导致内存泄漏或......

SciPy：在3D网格上插值散乱数据

SciPy: interpolate scattered data on 3D grid

python

interpolation

numpy

scipy

更新：