Nội suy nhanh của một trục mảng

Question

Tôi muốn nội suy một trục dữ liệu bên trong mảng 3 chiều. Các giá trị x đã cho cho các giá trị khác nhau hơi khác nhau nhưng tất cả chúng phải được ánh xạ tới cùng một giá trị x.

Kể từ khi x-giá trị nhất định là không giống nhau, hiện nay tôi làm như sau:

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) 

Sử dụng hai lồng nhau for -loops là gì ngạc nhiên rất chậm.

Có cách nào để cải thiện tốc độ không? Có thể bằng cách thực hiện một số phép thuật mảng hoặc song song của NumPy.

Answer 1

Tôi đã thực hiện một số thử nghiệm ban đầu và có vẻ như phần lớn thời gian của bạn được chi tiêu để tạo ra các hàm nội suy thực tế. Nó không có vẻ như chỉ vectorization sẽ tăng tốc độ nó lên một tấn, nhưng nó làm cho nó dễ dàng để song song (mà không tăng tốc độ năng suất). Đây là một ví dụ.

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)) 

Kết quả điển hình (lưu ý rằng việc thực hiện vectorized là khá nhiều chính xác cùng một tốc độ mà không có sự song song và rằng tôi có 2 lõi, vì vậy thời gian chạy của bạn nên được khoảng (thời gian ban đầu/# của lõi)):

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

Đừng làm cho tôi sai, có thể có cách giải thuật nhanh hơn, nhưng điều này có vẻ là cách dễ nhất để tăng tốc ngay lập tức, vì vấn đề là khó hiểu song song.