나는 그 루프를 제거하는 방법을 알고 있지만, 당신은 그것을 좋아하지 않을 것입니다.
문제 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)
결과 res
및 res2
배열이 정확히 동일하므로 접근 방법은 아마 작동합니다
어쨌든, 여기에 하나와 현재의 접근 방식을 비교, 행동 트릭입니다. 그러나 제가 말씀 드렸듯이, 모양이 (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)
로 다시 순서를 바꾸어 할 수 있습니다.