Python:并行运行嵌套循环，2D移动窗口 [英] Python: Running nested loop, 2D moving window, in Parallel
问题描述
我处理地形数据.对于一个特定问题，我用Python编写了一个函数，该函数使用特定大小的移动窗口在矩阵(高程网格)中滑动.然后，我必须在此窗口上进行分析，并将窗口中央的单元格设置为结果值.
I work with topographic data. For one particular problem, I have written a function in Python which uses a moving window of a particular size to zip through a matrix (grid of elevations). Then I have to perform an analysis on this window and set the cell at the center of the window a resulting value.
我的最终输出是与我的原始矩阵大小相同的矩阵，该矩阵根据我的分析进行了更改.这个问题在一个较小的区域上需要11个小时才能运行，因此我认为并行化内部循环会加快速度. 或者，也可能有一个聪明的矢量化解决方案...
My final output is a matrix the same size as my original matrix which has been altered according to my analysis. This problem takes 11 hours to run on a small area, so I thought parallelizing the inner loop would accelerate things. Alternatively, there may be a clever vectorization solution too...
请参见下面的函数，DEM是2D numpy数组，w是窗口的大小.
See my function below, DEM is a 2D numpy array, w is the size of the window.
def RMSH_det(DEM, w):
import numpy as np
from scipy import signal
[nrows, ncols] = np.shape(DEM)
#create an empty array to store result
rms = DEM*np.nan
# nw=(w*2)**2
# x = np.arange(0,nw)
for i in np.arange(w+1,nrows-w):
for j in np.arange(w+1,ncols-w):
d1 = np.int64(np.arange(i-w,i+w))
d2 = np.int64(np.arange(j-w,j+w))
win = DEM[d1[0]:d1[-1],d2[0]:d2[-1]]
if np.max(np.isnan(win)) == 1:
rms[i,j] = np.nan
else:
win = signal.detrend(win, type = 'linear')
z = np.reshape(win,-1)
nz = np.size(z)
rootms = np.sqrt(1 / (nz - 1) * np.sum((z-np.mean(z))**2))
rms[i,j] = rootms
return(rms)
我一直在SO/SE中寻找问题的解决方案，并遇到了许多嵌套的for循环示例并尝试并行运行它们.我一直努力调整我的代码以匹配示例，并希望获得一些帮助.解决此问题的方法将帮助我使用其他几个移动窗口功能.
I have scoured SO/SE for solutions to my question and come across many examples of nested for loops and trying to run them in parallel. I've struggled to adapt my code to match the examples and would appreciate some help. A solution to this problem would help me work with several other moving window functions I have.
到目前为止，我已经将内部循环移到了它自己的函数中，可以从外部循环内部调用它:
Thus far, I have moved the inner loop into it's own function, which can be called from within the outer loop:
def inLoop(i, w, DEM,rms,ncols):
for j in np.arange(w+1,ncols-w):
d1 = np.int64(np.arange(i-w,i+w))
d2 = np.int64(np.arange(j-w,j+w))
win = DEM[d1[0]:d1[-1],d2[0]:d2[-1]]
if np.max(np.isnan(win)) == 1:
rms[i,j] = np.nan
else:
win = signal.detrend(win, type = 'linear')
z = np.reshape(win,-1)
nz = np.size(z)
rootms = np.sqrt(1 / (nz - 1) * np.sum((z-np.mean(z))**2))
rms[i,j] = rootms
return(rms)
但是我不确定用需要输入到内部循环中的必要变量来编码对Pool的调用的正确方法.参见下面的外循环:
But I wasn't sure of the correct way of coding the call to Pool with the necessary variables that need to be input into the inner loop. See the outerloop below:
for i in np.arange(w+1,nrows-w):
number_of_workers = 8
with Pool(number_of_workers) as p:
#call the pool
p.starmap(inLoop, [i, w, DEM, rms, ncols])
剩余问题:
-
甚至可以通过并行化来优化此代码吗?
Can this code even be optimized by parallelizing?
如何成功存储并行嵌套循环的结果?
How do I successfully store the result of a parallelized nested-for loop?
推荐答案
使用Numba的解决方案
在某些情况下，如果支持您使用的所有功能，则此操作非常容易.在您的代码中，win = signal.detrend(win, type = 'linear')是您必须在Numba中实现的部分，因为不支持此功能.
A solution using Numba
In some cases this is very easy to do, if all functions which you use are supported. In your code win = signal.detrend(win, type = 'linear') is the part you have to implement in Numba, because this function isn't supported.
在Numba中实现下降趋势
如果您查看 detrend的源代码，并为您的问题提取相关部分，如下所示:
If you look at the source-code of detrend, and extract the relevant parts for your problem, it may look like this:
@nb.njit()
def detrend(w):
Npts=w.shape[0]
A=np.empty((Npts,2),dtype=w.dtype)
for i in range(Npts):
A[i,0]=1.*(i+1) / Npts
A[i,1]=1.
coef, resids, rank, s = np.linalg.lstsq(A, w.T)
out=w.T- np.dot(A, coef)
return out.T
我还为np.max(np.isnan(win)) == 1
@nb.njit()
def isnan(win):
for i in range(win.shape[0]):
for j in range(win.shape[1]):
if np.isnan(win[i,j]):
return True
return False
主要功能
正如我在这里使用Numba一样，并行化非常简单，只是外循环上的一个prange和
As I used Numba here, the parallelization is very simple, just a prange on the outer loop and
import numpy as np
import numba as nb
@nb.njit(parallel=True)
def RMSH_det_nb(DEM, w):
[nrows, ncols] = np.shape(DEM)
#create an empty array to store result
rms = DEM*np.nan
for i in nb.prange(w+1,nrows-w):
for j in range(w+1,ncols-w):
win = DEM[i-w:i+w-1,j-w:j+w-1]
if isnan(win):
rms[i,j] = np.nan
else:
win = detrend(win)
z = win.flatten()
nz = z.size
rootms = np.sqrt(1 / (nz - 1) * np.sum((z-np.mean(z))**2))
rms[i,j] = rootms
return rms
时间(小示例)
w = 10
DEM=np.random.rand(100, 100).astype(np.float32)
res1=RMSH_det(DEM, w)
res2=RMSH_det_nb(DEM, w)
print(np.allclose(res1,res2,equal_nan=True))
#True
%timeit res1=RMSH_det(DEM, w)
#1.59 s ± 72 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res2=RMSH_det_nb(DEM, w) #approx. 55 times faster
#29 ms ± 1.85 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
大型数组的时间
w = 10
DEM=np.random.rand(1355, 1165).astype(np.float32)
%timeit res2=RMSH_det_nb(DEM, w)
#6.63 s ± 21.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
[Edit]使用标准方程式的实现
此方法的数值精度较低.尽管此解决方案的速度要快得多.
This method has a lower numerical precision. Although this solution is quite a lot faster.
@nb.njit()
def isnan(win):
for i in range(win.shape[0]):
for j in range(win.shape[1]):
if win[i,j]==np.nan:
return True
return False
@nb.njit()
def detrend(w):
Npts=w.shape[0]
A=np.empty((Npts,2),dtype=w.dtype)
for i in range(Npts):
A[i,0]=1.*(i+1) / Npts
A[i,1]=1.
coef, resids, rank, s = np.linalg.lstsq(A, w.T)
out=w.T- np.dot(A, coef)
return out.T
@nb.njit()
def detrend_2(w,T1,A):
T2=np.dot(A.T,w.T)
coef=np.linalg.solve(T1,T2)
out=w.T- np.dot(A, coef)
return out.T
@nb.njit(parallel=True)
def RMSH_det_nb_normal_eq(DEM,w):
[nrows, ncols] = np.shape(DEM)
#create an empty array to store result
rms = DEM*np.nan
Npts=w*2-1
A=np.empty((Npts,2),dtype=DEM.dtype)
for i in range(Npts):
A[i,0]=1.*(i+1) / Npts
A[i,1]=1.
T1=np.dot(A.T,A)
nz = Npts**2
for i in nb.prange(w+1,nrows-w):
for j in range(w+1,ncols-w):
win = DEM[i-w:i+w-1,j-w:j+w-1]
if isnan(win):
rms[i,j] = np.nan
else:
win = detrend_2(win,T1,A)
rootms = np.sqrt(1 / (nz - 1) * np.sum((win-np.mean(win))**2))
rms[i,j] = rootms
return rms
时间
w = 10
DEM=np.random.rand(100, 100).astype(np.float32)
res1=RMSH_det(DEM, w)
res2=RMSH_det_nb(DEM, w)
print(np.allclose(res1,res2,equal_nan=True))
#True
%timeit res1=RMSH_det(DEM, w)
#1.59 s ± 72 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res2=RMSH_det_nb_normal_eq(DEM,w)
#7.97 ms ± 89.4 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
使用正态方程优化的解决方案
临时数组可以重复使用，以避免昂贵的内存分配，并且使用矩阵乘法的自定义实现.仅在很小的矩阵上才建议这样做，在大多数其他情况下，np.dot(sgeemm)会快得多.
Temporary arrays are reused to avoid costly memory allocations and a custom implementation for matrix multiplication is used. This is only recommendable for very small matrices, in most other cases np.dot (sgeemm) will be a lot faster.
@nb.njit()
def matmult_2(A,B,out):
for j in range(B.shape[1]):
acc1=nb.float32(0)
acc2=nb.float32(0)
for k in range(B.shape[0]):
acc1+=A[0,k]*B[k,j]
acc2+=A[1,k]*B[k,j]
out[0,j]=acc1
out[1,j]=acc2
return out
@nb.njit(fastmath=True)
def matmult_mod(A,B,w,out):
for j in range(B.shape[1]):
for i in range(A.shape[0]):
acc=nb.float32(0)
acc+=A[i,0]*B[0,j]+A[i,1]*B[1,j]
out[j,i]=acc-w[j,i]
return out
@nb.njit()
def detrend_2_opt(w,T1,A,Tempvar_1,Tempvar_2):
T2=matmult_2(A.T,w.T,Tempvar_1)
coef=np.linalg.solve(T1,T2)
return matmult_mod(A, coef,w,Tempvar_2)
@nb.njit(parallel=True)
def RMSH_det_nb_normal_eq_opt(DEM,w):
[nrows, ncols] = np.shape(DEM)
#create an empty array to store result
rms = DEM*np.nan
Npts=w*2-1
A=np.empty((Npts,2),dtype=DEM.dtype)
for i in range(Npts):
A[i,0]=1.*(i+1) / Npts
A[i,1]=1.
T1=np.dot(A.T,A)
nz = Npts**2
for i in nb.prange(w+1,nrows-w):
Tempvar_1=np.empty((2,Npts),dtype=DEM.dtype)
Tempvar_2=np.empty((Npts,Npts),dtype=DEM.dtype)
for j in range(w+1,ncols-w):
win = DEM[i-w:i+w-1,j-w:j+w-1]
if isnan(win):
rms[i,j] = np.nan
else:
win = detrend_2_opt(win,T1,A,Tempvar_1,Tempvar_2)
rootms = np.sqrt(1 / (nz - 1) * np.sum((win-np.mean(win))**2))
rms[i,j] = rootms
return rms
时间
w = 10
DEM=np.random.rand(100, 100).astype(np.float32)
res1=RMSH_det(DEM, w)
res2=RMSH_det_nb_normal_eq_opt(DEM, w)
print(np.allclose(res1,res2,equal_nan=True))
#True
%timeit res1=RMSH_det(DEM, w)
#1.59 s ± 72 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit res2=RMSH_det_nb_normal_eq_opt(DEM,w)
#4.66 ms ± 87.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
isnan的时间
此功能是完全其他的实现.如果NaN刚好位于数组的开头，则速度要快得多，但是无论如何，即使没有加速.我用小数组(大约窗口大小)和@ user3666197建议的大大小作为基准.
This function is a completely other implementation. It is much faster if a NaN is quite at the beginning of the array, but anyway even if not there is some speedup. I benchmarked it with small arrays (approx. window size) and a large size suggested by @user3666197.
case_1=np.full((20,20),np.nan)
case_2=np.full((20,20),0.)
case_2[10,10]=np.nan
case_3=np.full((20,20),0.)
case_4 = np.full( ( int( 1E4 ), int( 1E4 ) ),np.nan)
case_5 = np.ones( ( int( 1E4 ), int( 1E4 ) ) )
%timeit np.any(np.isnan(case_1))
%timeit np.any(np.isnan(case_2))
%timeit np.any(np.isnan(case_3))
%timeit np.any(np.isnan(case_4))
%timeit np.any(np.isnan(case_5))
#2.75 µs ± 73.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
#2.75 µs ± 46.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
#2.76 µs ± 32.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
#81.3 ms ± 2.97 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
#86.7 ms ± 2.26 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit isnan(case_1)
%timeit isnan(case_2)
%timeit isnan(case_3)
%timeit isnan(case_4)
%timeit isnan(case_5)
#244 ns ± 5.02 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
#357 ns ± 1.07 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
#475 ns ± 9.28 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
#235 ns ± 0.933 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
#58.8 ms ± 2.08 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
这篇关于Python:并行运行嵌套循环，2D移动窗口的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
