python – 更有效的循环方式？

2023年6月2日 198次阅读

我从一个更大的脚本中获得了一小段代码.我发现当调用函数t_area时,它负责大部分运行时间.我自己测试了这个功能,它并不慢,它需要花费很多时间,因为它必须运行的次数我相信.以下是调用函数的代码：

tri_area = np.zeros((numx,numy),dtype=float)
for jj in range(0,numy-1):
    for ii in range(0,numx-1):
      xp = x[ii,jj]
      yp = y[ii,jj]
      zp = surface[ii,jj]
      ap = np.array((xp,yp,zp))

      xp = xp+dx
      zp = surface[ii+1,jj]
      bp = np.array((xp,yp,zp))

      yp = yp+dx
      zp = surface[ii+1,jj+1]
      dp = np.array((xp,yp,zp))

      xp = xp-dx
      zp = surface[ii,jj+1]
      cp = np.array((xp,yp,zp))

      tri_area[ii,jj] = t_area(ap,bp,cp,dp)

这里使用的数组大小为216 x 217,x和y的值也是如此.我对python编码很新,我过去使用过MATLAB.所以我的问题是,有没有办法绕过这两个for循环,或者更有效的方式来运行这段代码？寻找任何帮助加快这一点！谢谢！

编辑：

感谢大家的帮助,这已经清除了很多混乱.我被问到循环中使用的函数t_area,下面的代码如下：

def t_area(a,b,c,d):
ab=b-a
ac=c-a
tri_area_a = 0.5*linalg.norm(np.cross(ab,ac))

db=b-d
dc=c-d
tri_area_d = 0.5*linalg.norm(np.cross(db,dc))

ba=a-b
bd=d-b
tri_area_b = 0.5*linalg.norm(np.cross(ba,bd))

ca=a-c
cd=d-c
tri_area_c = 0.5*linalg.norm(np.cross(ca,cd))

av_area = (tri_area_a + tri_area_b + tri_area_c + tri_area_d)*0.5
return(av_area)

对于令人困惑的记谱法感到抱歉,当时它有意义,现在回想起来我可能会改变它.谢谢！

最佳答案 在我们开始之前的一个警告.范围(0,numy-1),等于范围(numy-1),产生从0到numy-2的数字,不包括numy-1.那是因为你有从0到numy-2的numy-1值.虽然MATLAB具有基于1的索引,但Python基于0,因此在转换中对索引进行一些小心.考虑到你有tri_area = np.zeros((numx,numy),dtype = float),tri_area [ii,jj]永远不会以你设置循环的方式访问最后一行或列.因此,我怀疑正确的意图是写范围(numy).

由于功能t_area()是可矢量化的,因此您可以完全取消循环.矢量化意味着numpy通过处理引擎盖下的循环来同时对整个阵列应用一些操作,在那里它们将更快.

首先,我们在(m,n,3)数组中堆叠每个(i,j)元素的所有aps,其中(m,n)是x的大小.如果我们取两个(m,n,3)数组的叉积,默认情况下操作将应用于最后一个轴.这意味着np.cross(a,b)将对每个元素(i,j)采用[i,j]和b [i,j]中的3个数的叉积.类似地,np.linalg.norm(a,axis = 2)将为每个元素(i,j)计算a [i,j]中3个数的范数.这也将有效地减少我们的数组大小(m,n).这里有点谨慎,因为我们需要明确说明我们希望在第二轴上完成此操作.

请注意,在以下示例中,我的索引关系可能与您的索引关系不对应.使这项工作的最低限度是表面从x和y有一个额外的行和列.

import numpy as np

def _t_area(a, b, c):
    ab = b - a
    ac = c - a
    return 0.5 * np.linalg.norm(np.cross(ab, ac), axis=2)

def t_area(x, y, surface, dx):
    a = np.zeros((x.shape[0], y.shape[0], 3), dtype=float)
    b = np.zeros_like(a)
    c = np.zeros_like(a)
    d = np.zeros_like(a)

    a[...,0] = x
    a[...,1] = y
    a[...,2] = surface[:-1,:-1]

    b[...,0] = x + dx
    b[...,1] = y
    b[...,2] = surface[1:,:-1]

    c[...,0] = x
    c[...,1] = y + dx
    c[...,2] = surface[:-1,1:]

    d[...,0] = bp[...,0]
    d[...,1] = cp[...,1]
    d[...,2] = surface[1:,1:]

    # are you sure you didn't mean 0.25???
    return 0.5 * (_t_area(a, b, c) + _t_area(d, b, c) + _t_area(b, a, d) + _t_area(c, a, d))

nx, ny = 250, 250

dx = np.random.random()
x = np.random.random((nx, ny))
y = np.random.random((nx, ny))
surface = np.random.random((nx+1, ny+1))

tri_area = t_area(x, y, surface, dx)

此示例中的x支持索引0-249,而表面0-250. surface [： – 1],表面[0：-1]的简写,将返回从0开始直到最后一行的所有行,但不包括它. -1提供相同的功能并在MATLAB中结束.因此,surface [： – 1]将返回索引0-249的行.类似地,surface [1：]将返回索引1-250的行,这与表面[ii 1]的行相同.

注意：在知道t_area()可以完全矢量化之前,我已经写过这一节.因此,虽然这个答案的目的已经过时了,但我将把它留作遗产来表明如果函数不是可矢量化的话可以进行哪些优化.

不要为每个昂贵的元素调用函数,而应该将它传递给x,y,surface和dx并在内部迭代.这意味着只有一个函数调用和更少的开销.

此外,您不应该为每个循环创建ap,bp,cp和dp的数组,这又会增加开销.一旦在循环外部分配它们,只需更新它们的值.

最后一个改变应该是循环的顺序.默认情况下,Numpy数组是行主要的(而MATLAB是列专业),因此ii作为外部循环表现得更好.您不会注意到您的大小数组的差异,但是,嘿,为什么不呢？

总的来说,修改后的功能应如下所示.

def t_area(x, y, surface, dx):
    # I assume numx == x.shape[0]. If not, pass it as an extra argument.
    tri_area = np.zeros(x.shape, dtype=float)

    ap = np.zeros((3,), dtype=float)
    bp = np.zeros_like(ap)
    cp = np.zeros_like(ap)
    dp = np.zeros_like(ap)

    for ii in range(x.shape[0]-1): # do you really want range(numx-1) or just range(numx)?
        for jj in range(x.shape[1]-1):
            xp = x[ii,jj]
            yp = y[ii,jj]
            zp = surface[ii,jj]
            ap[:] = (xp, yp, zp)

            # get `bp`, `cp` and `dp` in a similar manner and compute `tri_area[ii,jj]`
点赞