2017-04-03 1 views
4

제 질문은 여기에 보이는 코드 응답을 확장합니다 : Interpolating a 3d array in Python. How to avoid for loops?. 관련 원액 코드는 다음과 같습니다 :파이썬에서 3D 배열 보간 확장

import numpy as np 
from scipy.interpolate import interp1d 
array = np.random.randint(0, 9, size=(100, 100, 100)) 
x = np.linspace(0, 100, 100) 
x_new = np.linspace(0, 100, 1000) 
new_array = interp1d(x, array, axis=0)(x_new) 
new_array.shape # -> (1000, 100, 100) 

x_new는 상수 1 차원 배열이 때하지만 내가 x_new 는 상수 1 차원 배열 무엇을하지 않은 경우, 대신 인덱스에 따라 접근 방식은 위의 좋은 작품 다른 3 차원 배열의 위도/경도 차원 내 x_new 크기는 355x195x192 (시간 x lat x 길이)이고 현재는 위도와 경도 차원을 반복 할 때 사용합니다. x_new는 각 위도/경도 쌍마다 다르기 때문에 아래에서 볼 수 있듯이 어떻게 루핑을 피할 수 있습니까? 내 루프 프로세스는 불행히도 몇 시간이 걸립니다 ...

x_new=(np.argsort(np.argsort(modell, 0), 0).astype(float) + 1)/np.size(modell, 0) 
## x_new is shape 355x195x192 
## pobs is shape 355x1 
## prism_aligned_tmax_sorted is shape 355x195x192 
interp_func = interpolate.interp1d(pobs, prism_aligned_tmax_sorted,axis=0) 
tmaxmod = np.empty((355, 195, 192,)) 
tmaxmod[:] = np.NAN          
for latt in range(0, 195): 
    for lonn in range(0, 192): 
     temp = interp_func(x_new[:,latt,lonn]) 
     tmaxmod[:,latt,lonn] = temp[:,latt,lonn] 

감사합니다.

답변

1

나는 그 루프를 제거하는 방법을 알고 있지만, 당신은 그것을 좋아하지 않을 것입니다.

문제 interp1d의 사용은 본질적으로는 F를 2 차원 어레이 형상이 x 각 스칼라에 대한 매트릭스 값의 1 차원 영역에 걸쳐 보간 함수, 즉 F(x) 기능 어디를 제공한다는 것이다. 당신이하려고하는 보간법은 오히려 이것입니다 : 각각의 (lat,lon) 쌍에 대한 개별 보간기를 생성하는 것입니다. 이것은 F_(lat,lon)(x)의 줄에 더 있습니다.

이 문제가 발생하는 이유는 각 쿼리 포인트에 대해 행렬 값 F(x)을 계산하는 것이지만 하나의 요소 (요소 이 쌍에 해당하는 쿼리 포인트는 [lat,lon]). 따라서 모든 관련없는 함수 값을 계산하는 불필요한 계산을하고있는 것입니다. 문제는 더 효율적인 방법이 있는지 확신 할 수 없다는 것입니다.

사용 사례는 등 뒤에서 적절한 메모리로 고정시킬 수 있습니다. 루프가 몇 시간 동안 실행된다는 사실은 이것이 실제로 당신의 유스 케이스에서는 가능하지 않다는 것을 암시하지만, 어쨌든이 해결책을 보여줄 것입니다. 아이디어는 3d 배열을 2 차원으로 바꾸고,이 모양으로 보간 한 다음 보간 결과의 유효한 2d 부분 공간을 따라 대각선 요소를 취하는 것입니다. 모든 쿼리 포인트에 대해 관련성이없는 모든 행렬 요소를 계산하지만 배열을 반복하지 않고 단일 인덱싱 단계로 관련 행렬 요소를 추출 할 수 있습니다. 이 비용은 사용 가능한 RAM에 맞지 않는 훨씬 더 큰 보조 어레이를 만드는 것입니다.

import numpy as np 
from scipy.interpolate import interp1d 
arr = np.random.randint(0, 9, size=(3, 4, 5)) 
x = np.linspace(0, 10, 3) 
x_new = np.random.rand(6,4,5)*10 

## x is shape 3 
## arr is shape 3x4x5 
## x_new is shape 6x4x5 

# original, loopy approach 
interp_func = interp1d(x, arr, axis=0) 
res = np.empty((6, 4, 5)) 
for lat in range(res.shape[1]): 
    for lon in range(res.shape[2]): 
     temp = interp_func(x_new[:,lat,lon]) # shape (6,4,5) each iteration 
     res[:,lat,lon] = temp[:,lat,lon] 

# new, vectorized approach 
arr2 = arr.reshape(arr.shape[0],-1) # shape (3,20) 
interp_func2 = interp1d(x,arr2,axis=0) 
x_new2 = x_new.reshape(x_new.shape[0],-1) # shape (6,20) 
temp = interp_func2(x_new2) # shape (6,20,20): 20 larger than original! 
s = x_new2.shape[1] # 20, used for fancy indexing ranges 
res2 = temp[:,range(s),range(s)].reshape(res.shape) # shape (6,20) -> (6,4,5) 

결과 resres2 배열이 정확히 동일하므로 접근 방법은 아마 작동합니다

어쨌든, 여기에 하나와 현재의 접근 방식을 비교, 행동 트릭입니다. 그러나 제가 말씀 드렸듯이, 모양이 (nx,nlat,nlon) 인 질의 배열에 대해서는 대개 메모리가 많이 필요한 모양 인 (nx,nlat*nlon,nlat*nlon)의 보조 배열이 필요합니다.


난 그냥 하나 하나 귀하의 1D 보간을 수행하고 생각할 수있는 유일한 엄격한 대안 : 더블 루프에서 nlat*nlon 보간을 정의. 이것은 인터폴 레이터를 생성하는 데 훨씬 많은 오버 헤드를 가져 오지만, 반면에 당신은 불필요한 작업을 수행하지 않고 보간 된 배열 값을 계산하여 버립니다.

마지막으로, 사용에 따라 다 변수 보간법을 사용하는 것이 좋습니다 (interpolate.interpnd 또는 interpolate.griddata). 여러분의 함수가 위도와 경도의 함수로서 매끄럽다 고 가정하면 더 높은 차원에서 전체 데이터 집합을 보간하는 것이 합리적 일 수 있습니다. 이렇게하면 인터폴 레이터를 한 번만 생성하고 불필요한 보풀이없는 정확한 위치에서 쿼리해야합니다.


현재 구현을 고수하면 결국 보간 축을 마지막 위치로 이동하여 성능을 크게 향상시킬 수 있습니다. 이 방식으로 모든 벡터화 된 연산은 연속적인 메모리 블록 (기본 C 메모리 순서로 가정)에서 작동하며 "1 차원 배열 모음"철학과 잘 어울립니다. 원래 순서를 복원해야 할 경우 그래서 당신은,

arr = arr.transpose(1,2,0) # shape (4,5,3) 
interp_func = interp1d(x, arr, axis=-1) 
... 
for lat ...: 
    for lon ...: 
     res[lat,lon,:] = temp[lat,lon,:] # shape (4,5,6) 

의 라인을 따라 뭔가를해야 마침내 res.transpose(2,0,1)로 다시 순서를 바꾸어 할 수 있습니다.