有啥方法可以比 for 循环更快地遍历数组吗？

Posted 2023-04-18

【中文标题】有啥方法可以比 for 循环更快地遍历数组吗？

【英文标题】：Is there any way I can iterate through an array faster than a for loop?有什么方法可以比 for 循环更快地遍历数组吗？

【发布时间】：2019-03-21 04:33:27

【问题描述】：

我正在编写一个代码，将天文地图上的像素通量与另一幅地图上的相应区域进行比较。两张地图都是 numpy 数据数组。

为此，我需要将第一张地图 (Av) 上的像素索引转换为它们在天空坐标上的等效值，然后将这些天空坐标转换为它们在第二张地图 (CO) 上等效的像素索引。然后，我缩放第二张地图的通量以匹配第一张地图的值。在那之后，我必须继续处理数据。

问题在于，第一张地图上有数千个像素，代码需要很长时间才能完成它应该做的事情，这很麻烦故障排除。我发现这部分代码中最慢的是for循环。

有什么方法可以遍历一个 numpy 数组，能够使用索引并计算每个像素的数据，比 for 循环更快？有没有更好的方法来做到这一点有吗？

在伪代码中，我的代码是这样的：

for pixel i,j in 1st map:
    sky_x1,sky_y1 = pixel_2_skycoord(i,j)
    i2,j2 = skycoord_2_pixel(sky_x1,sky_y1)

    Avmap.append(Avflux[i,j])
    COmap.append(COflux[i2,j2]*scale)

实际代码是：

for i in xrange(0,sAv_y-1):
    for j in xrange(0,sAv_x-1):
        if not np.isnan(Avdata[i,j]):           

            y,x=wcs.utils.skycoord_to_pixel(wcs.utils.pixel_to_skycoord(i,j,wAv,0),wcs=wCO)

            x=x.astype(int)+0 #the zero is because i don't understand the problem with numpy but it fixes it anyway
            y=y.astype(int)+0 #i couldn't get the number from an array with 1 value but adding zero resolves it somehow
            COflux=COdata[x,y]

            ylist.append(Avdata[i,j])
            xlist.append(COflux*(AvArea/COArea))

【问题讨论】：

顺便说一句：这是不正确的：if (Avdata[i,j] is not np.nan) 你必须使用类似np.isnan的东西

【参考方案1】：

这里的罪魁祸首是两个 for 循环。 Numpy 有许多功能可以防止使用 for 循环来允许快速编译代码。诀窍是矢量化您的代码。

您可以查看 numpy 的 meshgrid 函数将此数据转换为矢量化形式，然后您可以使用 this SO question 之类的东西将任意函数应用于该向量。

类似的东西：

x_width = 15
y_width = 10

x, y = np.meshgrid(range(x_width), range(y_width))

def translate(x, y, x_o, y_o):
    x_new = x + x_o
    y_new = y + y_o
    return x_new, y_new

x_new, y_new = translate(x, y, 3, 3)

x_new[4,5], y[4,5]

(8, 4)

【讨论】：

这个，以及你分享的问题，真的很有帮助。虽然这不是一个解决方案，因为现在我不是遍历数组而是在整个内存块上应用一个函数，我的 8GB RAM 和更多的交换内存不能胜任这项任务。我发现当您可以腾出 RAM 以节省时间时，此解决方案有效。我以后肯定会用这个。

当然，这是您大部分时间的权衡，CPU 与 RAM。虽然映射到屏幕上像素的东西会填满你的内存，但我很惊讶。可能还有其他优化（可能在原地工作而不是复制数组）。我很高兴这有帮助。

【参考方案2】：

您必须避免循环，并在底层 C 代码、Numpy 或 Astropy 中进行大量计算以进行天空/像素转换。使用astropy.wcs 有几个选项可以做到这一点。

第一个是SkyCoord。让我们首先为您的像素索引创建一个值网格：

In [30]: xx, yy = np.mgrid[:5, :5] 
    ...: xx, yy 
Out[30]: 
(array([[0, 0, 0, 0, 0],
        [1, 1, 1, 1, 1],
        [2, 2, 2, 2, 2],
        [3, 3, 3, 3, 3],
        [4, 4, 4, 4, 4]]), array([[0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4],
        [0, 1, 2, 3, 4]]))

现在我们可以从像素索引创建SkyCoord 对象（这是一个 Numpy 数组子类），并使用 wcs：

In [33]: from astropy.coordinates import SkyCoord 
    ...: sky = SkyCoord.from_pixel(xx, yy, wcs) 
    ...: sky                                                                                   
Out[33]: 
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
    [[(53.17127889, -27.78771333), (53.17127889, -27.78765778),
      (53.17127889, -27.78760222), (53.17127889, -27.78754667),
      (53.17127889, -27.78749111)],
      ....

请注意，这里使用的是wcs.utils.skycoord_to_pixel。该对象还具有使用 wcs 投影到像素的方法。出于实际目的，我将在这里相同：

In [34]: sky.to_pixel(wcs)                                                                     
Out[34]: 
(array([[ 0.00000000e+00, -1.11022302e-16, -2.22044605e-16,
         -3.33066907e-16,  1.13149046e-10],
        ...
        [ 4.00000000e+00,  4.00000000e+00,  4.00000000e+00,
          4.00000000e+00,  4.00000000e+00]]),
 array([[-6.31503738e-11,  1.00000000e+00,  2.00000000e+00,
          3.00000000e+00,  4.00000000e+00],
        ...
        [-1.11457732e-10,  1.00000000e+00,  2.00000000e+00,
          3.00000000e+00,  4.00000000e+00]]))

我们得到一个新的 x 和 y 索引的浮点值元组。因此，您需要将这些值四舍五入并转换为 int 以将其用作数组索引。

第二种选择是使用较低级别的函数，例如wcs.pixel_to_world_values 和 wcs.world_to_pixel_values，它们接受 Nx2 数组并返回：

In [37]: wcs.pixel_to_world_values(np.array([xx.ravel(), yy.ravel()]).T)                       
Out[37]: 
array([[ 53.17127889, -27.78771333],
       [ 53.17127889, -27.78765778],
       [ 53.17127889, -27.78760222],
       [ 53.17127889, -27.78754667],
       ...

【讨论】：

谢谢，我发现这些软件包很难理解，因为在快速使用它们之前我无法花时间研究它们的工作原理。