使用移动内核 window 时 Numba 失败:TypingError,定义错误

Numba fails when using moving kernel window: TypingError, definition error

我正在从事一个项目,该项目要求我构建一系列函数,这些函数使用移动内核 window 来处理存储在矩阵中的海拔数据。

我最初提出的如何优化两个嵌套 for 循环的问题在这里 得到了解答。解决方案涉及使用 Numba 并行化我的代码。当我尝试使另一个移动内核 window 函数适应 numba 时,我 运行 遇到了 numpy.gradient 的问题,目前 numba 不支持该功能。所以我把我的函数分成两部分:A) 一个不在 Numba 中的预处理函数和 B) 用 Numba 编写的主要函数。我现在正在尝试使用 numba 使函数 运行 正常工作。我不断收到错误消息(如下所示),在搜索堆栈交换和互联网寻找解决方案数周后,我现在非常沮丧。

功能如下:

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)

其中 returns 我的第二个函数的输入: 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)

第二个函数现在出现以下错误:

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)

代码在第 40 行失败,即 tempx = cx[(d1[0]):(d1[-1]),(d2[0]):(d2[-1])]。我不确定发生了什么。我已经确保所有变量在整个过程中都是一致的类型,或者我是这么认为的。任何人都可以就解决此问题的方法提供任何见解或帮助,我们将不胜感激。

运行 类似的函数将 1165 x 1355 矩阵的并行切割操作从 6 分钟缩短到几秒。我希望通过此功能获得类似的加速。 非常感谢

经过大量的反复试验并深入研究 TraceBack,我成功地解决了这个问题。我的新功能如下:

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

主要变化是我:

  1. 消除了对 reshape 的调用,而是使用 tempx.flatten
  2. 我的 if 语句中对 max 的调用已更改为 .sum()

现在的代码 运行 很有魅力。与 6 分钟相比,非并行代码需要 6 秒 运行.