1つの配列軸の高速補間

Question

3次元配列内で1軸のデータを補間したい。指定された異なる値の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) 

を2つの入れ子for -loopsを使用すると、当然非常に遅いです。

速度を向上させる方法はありますか？たぶんNumPy配列のマジックや並列化を行うことによって。

Answer 1

私はいくつかの初期テストを行っており、膨大な時間が実際の補間関数を生成するのに費やされているようです。ベクトル化だけでは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 

私は間違ってはいません。これをより迅速に行うためのアルゴリズム的な方法があるかもしれませんが、問題は恥ずかしく並行しているため、即座にスピードアップするのが最も簡単な方法です。