I'm working on a project that requires me to build a series of functions that use a moving kernel window to manipulate elevation data stored in a matrix.
My original question of how to optimize two nested for loops was answered here link. The solution involved parallelizing my code with Numba. When I tried to adapt another moving kernel window function to numba, I ran into an issue with numpy.gradient which is currently unsupported by numba. So I broke my function in two: A) a preprocessing function not in Numba and B) the main function written in Numba. I am now trying to get the function run with numba to work properly. I continue to get an error (shown below) and am now pretty frustrated after weeks scouring the stack exchange and internet for a solution.
Functions are shown below:
first function
def DCE_preprocess(DEM, cellsize, w):
[nrows, ncols] = np.shape(DEM)
#initiate an empty array same size as dem
rms = DEM*np.nan
# rms = np.float32(rms)
# #compute the directional cosines
[fx, fy] = np.gradient(DEM, cellsize, cellsize)
grad = np.sqrt(fx**2 + fy**2)
asp = np.arctan2(fy, fx)
grad=np.pi/2-np.arctan(grad) #normal of steepest slope
asp[asp<np.pi]=asp[asp<np.pi]+[np.pi/2]
asp[asp<0]=asp[asp<0]+[2*np.pi]
#spherical to cartesian conversion
r = 1
cy = r * np.cos(grad) * np.sin(asp)
cx = r * np.cos(grad) * np.cos(asp)
cz = r * np.sin(grad)
return(cx,cy,cz)
which returns the input for my second function:
second function
eps = np.finfo(float).eps
@nb.njit(parallel=True)
def DC_eig_par(DEM,w,cx,cy,cz,eps):
[nrows, ncols] = np.shape(DEM)
#
# #initiate an empty array same size as dem
rms = DEM*np.nan
rms.astype(np.float32)
#cycling through the DEM
nw=(w*2)**2
for i in nb.prange(w+1,nrows-w):
for j in range(w+1,(ncols-w)):
d1=nb.int32(np.linspace(i-w,i+w,11))
d2=nb.int32(np.linspace(j-w,j+w,11))
tempx = cx[(d1[0]):(d1[-1]),(d2[0]):(d2[-1])]
tx=np.reshape(tempx,-1)
tempy = cy[(d1[0]):(d1[-1]),(d2[0]):(d2[-1])]
ty=np.reshape(tempy,-1)
tempz = np.empty([10,10], dtype = np.float32)
tempz = cz[(d1[0]):(d1[-1]),(d2[0]):(d2[-1])]
tz=np.reshape(tempz,-1)
if np.max(np.isnan(np.concatenate((tx,ty,tz)))) == 0:
T=np.array([[np.sum(tx**2), np.sum(tx*ty), np.sum(tx*tz)],
[np.sum(ty*tx), np.sum(ty**2), np.sum(ty*tz)],
[np.sum(tz*tx), np.sum(tz*ty), np.sum(tz**2)]])
[Te,_] = np.linalg.eig(T) # this step is a bit different from the matlab version b/c np.eig outputs two values.
l = (Te/nw)
l[l<eps] = 0
rms[i,j] = 1/np.log(l[0]/l[1])
else:
rms[i,j] = np.nan
return(rms)
The second function now gives me the following error:
Traceback (most recent call last):
File "<ipython-input-65-af12df87ecaa>", line 8, in <module>
rp = DC_eig_par(DEM, w,cx,cy,cz,eps)
File "C:\ProgramData\Anaconda3\lib\site-packages\numba\dispatcher.py", line 351, in _compile_for_args
error_rewrite(e, 'typing')
File "C:\ProgramData\Anaconda3\lib\site-packages\numba\dispatcher.py", line 318, in error_rewrite
reraise(type(e), e, None)
File "C:\ProgramData\Anaconda3\lib\site-packages\numba\six.py", line 658, in reraise
raise value.with_traceback(tb)
TypingError: Invalid use of Function(<built-in function getitem>) with argument(s) of type(s): (int32, Literal[int](0))
* parameterized
In definition 0:
All templates rejected with literals.
In definition 1:
All templates rejected without literals.
In definition 2:
All templates rejected with literals.
In definition 3:
All templates rejected without literals.
In definition 4:
All templates rejected with literals.
In definition 5:
All templates rejected without literals.
In definition 6:
All templates rejected with literals.
In definition 7:
All templates rejected without literals.
In definition 8:
All templates rejected with literals.
In definition 9:
All templates rejected without literals.
This error is usually caused by passing an argument of a type that is unsupported by the named function.
[1] During: typing of intrinsic-call at U:\GHP\Projects\NSF - Morris Landslides\Code\Developmemnt\Wavelets\DC_eig_par.py (40)
[2] During: typing of static-get-item at U:\GHP\Projects\NSF - Morris Landslides\Code\Developmemnt\Wavelets\DC_eig_par.py (40)
The code fails on line 40 which is tempx = cx[(d1[0]):(d1[-1]),(d2[0]):(d2[-1])]. I'm not sure what's going on. I've made sure that all variables are of a consistent type throughout, or so I thought. Any insights or help anyone could offer as to a way to resolve this issue would be greatly appreciated.
Running similar functions in parallel cuts operation from 6 minutes for a 1165 x 1355 matrix to seconds. I hope to have a similar speedup with this function. Thanks so much
I was able to successfully resolve this issue after a lot of trial and error and digging more into the TraceBack. My new function reads as follows:
@nb.njit(parallel=True)
def DC_eig_par(DEM,w,cx,cy,cz,eps):
[nrows, ncols] = np.shape(DEM)
#
# #initiate an empty array same size as dem
rms = DEM*np.nan
rms.astype(np.float32)
#Compute RMS cycling through the DEM
nw=(w*2)**2
for i in nb.prange(w+1,nrows-w):
for j in range(w+1,(ncols-w)):
# d1=np.int16(np.linspace(i-w,i+w,11))
# d2=np.int16(np.linspace(j-w,j+w,11))
tempx = cx[i-w:i+w,j-w:j+w]
tx = tempx.flatten()
tempy = cy[i-w:i+w,j-w:j+w]
ty = tempy.flatten()
tempz = cz[i-w:i+w,j-w:j+w]
tz = tempz.flatten()
if (np.isnan(np.concatenate((tx,ty,tz)))).sum() == 0:
T=np.array([[np.sum(tx**2), np.sum(tx*ty), np.sum(tx*tz)],
[np.sum(ty*tx), np.sum(ty**2), np.sum(ty*tz)],
[np.sum(tz*tx), np.sum(tz*ty), np.sum(tz**2)]])
[Te,_] = np.linalg.eig(T) # this step is a bit different from the matlab version b/c np.eig outputs two values.
l = (Te/nw)
l[l<eps] = 0
rms[i,j] = 1/np.log(l[0]/l[1])
else:
rms[i,j] = np.nan
return(rms)
Major changes are that I:
reshape and instead used tempx.flattenmax within my if statement was changed to .sum()The code now runs like a charm. In 6 seconds as compared with the 6 minutes the non-parallelized code takes to run.