快速插值单个数组轴

Question

我希望您能够在三维数组中插值一个轴的数据。不同值的给定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循环很慢，这并不奇怪。有没有办法提高速度？也许可以通过一些NumPy数组魔术或并行化来实现。

- leviathan

---

你能添加一个 arr 的示例吗？ - DJK

我的示例中有一个错误，现在已经修复了。现在应该给出 arr。 - leviathan

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

不要误会，也许有一种算法方法可以更快地完成这个任务，但是由于这个问题很容易并行处理，这似乎是立即加速的最简单方法。

- Saedeas