하나의 배열 축을 빠르게 보간합니다.

Question

3 차원 배열 내부에 한 축의 데이터를 보간하고 싶습니다. 다른 값에 대해 지정된 x 값은 약간 다르지만 모두 동일한 x 값에 매핑되어야합니다. 주어진 x 값 이후

현재 내가 다음을 수행, 동일하지 않은 두 개의 중첩 된 for -loops를 사용

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) 

는 당연히 매우 느립니다.

속도를 개선 할 방법이 있습니까? 어쩌면 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 

오해하지 마세요.이 작업을보다 신속하게 처리 할 수있는 알고리즘 방법이있을 수 있지만 문제가 당황 스럽기 때문에 즉시 속도를 낼 수있는 가장 쉬운 방법 인 것 같습니다.