在Python中扩展插入3d数组

Question

我的问题扩展了这里看到的代码响应:在Python中插入一个3d数组.如何避免for循环？.相关的原始解决方案代码如下:

import numpy as np
from scipy.interpolate import interp1d
array = np.random.randint(0, 9, size=(100, 100, 100))
x = np.linspace(0, 100, 100)
x_new = np.linspace(0, 100, 1000)
new_array = interp1d(x, array, axis=0)(x_new)
new_array.shape # -> (1000, 100, 100)

当x_new是一个常数1-d数组时,上述方法很有效,但是如果我的x_new 不是常数1-d数组,而是取决于另一个3-d数组中纬度/经度维度的索引.我的x_new大小为355x195x192(时间x纬度x长),现在我正在通过纬度和经度维度进行双循环.由于每个纬度/经度对的x_new不同,如何避免循环,如下所示？不幸的是,我的循环过程需要几个小时......

x_new=(np.argsort(np.argsort(modell, 0), 0).astype(float) + 1) / np.size(modell, 0)
## x_new is shape 355x195x192
## pobs is shape 355x1
## prism_aligned_tmax_sorted  is shape 355x195x192
interp_func = interpolate.interp1d(pobs, prism_aligned_tmax_sorted,axis=0)
tmaxmod = np.empty((355, 195, 192,))
tmaxmod[:] = np.NAN                                    
for latt in range(0, 195):
    for lonn in range(0, 192):
        temp = interp_func(x_new[:,latt,lonn])
        tmaxmod[:,latt,lonn] = temp[:,latt,lonn]

感谢您的所有帮助!

Answer 1

我知道如何摆脱这些循环，但你不会喜欢它。

问题是，这种使用interp1d本质上为您提供了一个在 1d 域上插值的矩阵值函数，F(x)即每个标量x都有一个 2d 数组形状的函数F。您尝试执行的插值实际上是这样的：为(lat,lon)您的每一对创建一个单独的插值器。这更符合F_(lat,lon)(x).

这是一个问题的原因是，对于您的用例，您正在计算每个查询点的矩阵值F(x)，但随后继续丢弃除单个矩阵元素之外的所有矩阵元素（[lat,lon]查询点对应的元素）到这一对）。所以你正在做一堆不必要的计算来计算所有这些不相关的函数值。问题是我不确定是否有更有效的方法。

您的用例可以通过您背后的适当记忆来修复。您的循环运行数小时的事实表明，这对于您的用例来说实际上是不可能的，但无论如何我都会展示这个解决方案。这个想法是将你的 3d 数组变成 2d 数组，用这个形状进行插值，然后沿着插值结果的有效 2d 子空间获取对角线元素。您仍将计算每个查询点的每个不相关的矩阵元素，但您将能够通过单个索引步骤提取相关的矩阵元素，而不是循环数组。这样做的代价是创建一个更大的辅助阵列，这很可能不适合您的可用 RAM。

无论如何，这是实际的技巧，将您当前的方法与以前的方法进行比较：

import numpy as np
from scipy.interpolate import interp1d
arr = np.random.randint(0, 9, size=(3, 4, 5))
x = np.linspace(0, 10, 3)
x_new = np.random.rand(6,4,5)*10

## x is shape 3 
## arr is shape  3x4x5
## x_new is shape  6x4x5

# original, loopy approach
interp_func = interp1d(x, arr, axis=0)
res = np.empty((6, 4, 5))
for lat in range(res.shape[1]):
    for lon in range(res.shape[2]):
        temp = interp_func(x_new[:,lat,lon]) # shape (6,4,5) each iteration
        res[:,lat,lon] = temp[:,lat,lon]

# new, vectorized approach
arr2 = arr.reshape(arr.shape[0],-1) # shape (3,20)
interp_func2 = interp1d(x,arr2,axis=0)
x_new2 = x_new.reshape(x_new.shape[0],-1) # shape (6,20)
temp = interp_func2(x_new2) # shape (6,20,20): 20 larger than original!
s = x_new2.shape[1] # 20, used for fancy indexing ranges
res2 = temp[:,range(s),range(s)].reshape(res.shape) # shape (6,20) -> (6,4,5)

结果res和res2数组完全相同，因此该方法可能有效。但正如我所说，对于 shape 的查询数组，(nx,nlat,nlon)我们需要一个 shape 的辅助数组(nx,nlat*nlon,nlat*nlon)，它通常需要巨大的内存。

我能想到的唯一严格的替代方案就是一一执行一维插值：nlat*nlon在双循环中定义插值器。这将产生更大的创建插值器的开销，但另一方面，您不会做一堆不必要的工作来计算然后丢弃的插值数组值。

最后，根据您的使用情况，您应该考虑使用多元插值（我正在考虑interpolate.interpnd或interpolate.griddata）。假设您的函数作为纬度和经度的函数也是平滑的，那么在更高维度中插入完整数据集可能是有意义的。这样，您需要创建一次插值器，并在您需要的点进行精确查询，而不会出现不必要的麻烦。

如果您最终坚持使用当前的实现，则可以通过将插值轴移动到最后一个位置来极大地提高性能。这样，每个向量化操作都作用于连续的内存块（假设默认的 C 内存顺序），这非常适合“一维数组的集合”哲学。所以你应该做一些类似的事情

arr = arr.transpose(1,2,0) # shape (4,5,3)
interp_func = interp1d(x, arr, axis=-1)
...
for lat ...:
    for lon ...:
        res[lat,lon,:] = temp[lat,lon,:] # shape (4,5,6)

如果需要恢复原来的顺序，最后可以用 调回顺序res.transpose(2,0,1)。