python - 一个阵列轴的快速插补

Question

我想在 3 维数组中插入一个数据轴。不同值的给定 x 值略有不同，但它们都应映射到相同的 x 值。

由于给定的 x 值不相同，目前我执行以下操作:

import numpy as np
from scipy import interpolate

axes_have = np.ones( ( 2, 72, 2001 ) )
axes_have *= np.linspace( 0, 100, 2001 )[None,None,:]
axes_have += np.linspace( -0.3, 0.3, 144 ).reshape( ( 2, 72 ) )[:,:,None]

arr = np.sin( axes_have )
arr *= np.random.random( (2, 72 ) )[:,:,None]

axis_want = np.linspace( 0, 100, 201 )    

arr_ip = np.zeros( ( 2, 72, 201 ) )
for i in range( arr.shape[0] ):
    for j in range( arr.shape[1] ):
         ip_func = interpolate.PchipInterpolator( axes_have[i,j,:], arr[i,j,:], extrapolate=True )
         arr_ip[i,j,:] = ip_func( axis_want )

使用两个嵌套 for -loops 不出所料地非常慢。

有没有办法提高速度？也许通过做一些 NumPy 数组魔术或并行化。

Answer 1

我已经做了一些初步测试，似乎你的大部分时间都花在了生成实际的插值函数上。似乎只是矢量化不会使速度提高一吨，但它确实使并行化变得容易(这确实提高了产量)。这是一个例子。

import numpy as np
from scipy import interpolate
import timeit
import multiprocessing



def myfunc(arg):
    x, y = arg
    return interpolate.PchipInterpolator(x,
                                         y,
                                         extrapolate=True)

p = multiprocessing.Pool(processes=8)
axes_have = np.ones((2, 72, 2001))
axes_have *= np.linspace(0, 100, 2001)[None, None, :]
axes_have += np.linspace(-0.3, 0.3, 144).reshape((2, 72))[:, :, None]

arr = np.sin(axes_have)
arr *= np.random.random((2, 72))[:, :, None]

axis_want = np.linspace(0, 100, 201)

arr_ip = np.zeros((2, 72, 201))
s_time1 = timeit.default_timer()
for i in range(arr.shape[0]):
    for j in range(arr.shape[1]):
         ip_func = interpolate.PchipInterpolator(axes_have[i, j, :],
                                                 arr[i, j, :],
                                                 extrapolate=True)
         arr_ip[i, j, :] = ip_func(axis_want)
elapsed1 = timeit.default_timer() - s_time1

s_time2 = timeit.default_timer()
flatten_have = [y for x in axes_have for y in x]
flatten_arr = [y for x in arr for y in x]
interp_time = timeit.default_timer()
funcs = p.map(myfunc, zip(flatten_have, flatten_arr))
interp_elapsed = timeit.default_timer() - interp_time
arr_ip = np.asarray([func(axis_want) for func in funcs]).reshape(2, 72, 201)
elapsed2 = timeit.default_timer() - s_time2

print("Elapsed 1: {}".format(elapsed1))
print("Elapsed 2: {}".format(elapsed2))
print("Elapsed interp: {}".format(interp_elapsed))

典型结果(请注意，在没有并行化的情况下，矢量化实现的速度几乎完全相同，而且我有 2 个内核，因此您的运行时间应该大致为(原始时间/内核数)):

Elapsed 1: 10.4025919437
Elapsed 2: 5.11409401894
Elapsed interp: 5.10987687111

不要误会我的意思，可能有一种算法方法可以更快地做到这一点，但这似乎是产生立即加速的最简单方法，因为问题是令人尴尬的并行。