在 for 循环中加快性能重写数组答案

【问题标题】：Speed performance rewriting array in for loop在 for 循环中加快性能重写数组
【发布时间】：2018-09-27 21:59:38
【问题描述】：

我有一个带有shape = (500, 500) 的二维数据集。从给定位置(x_0, y_0) 我想将每个元素/像素的距离映射到给定位置。我通过确定与(x_0, y_0) 的所有唯一距离并使用整数映射它们来做到这一点。 6 x 6 数据集的这种映射如下所示：

[9 8 7 6 7 8]
[8 5 4 3 4 5]
[7 4 2 1 2 4]
[6 3 1 0 1 3]
[7 4 2 1 2 4]
[8 5 4 3 4 5]

其中整数对应于存储在以下数组中的唯一距离：

[0.  1.  1.41421356  2.  2.23606798  2.82842712  3.  3.16227766  3.60555128  4.24264069]

确定这些距离的代码如下：

def func(data, (x_0,y_0)):
  y, x = numpy.indices((data.shape))
  r = numpy.sqrt((x - x_0)**2 + (y - y_0)**2)

  float_values = numpy.unique(r.ravel())  # Unique already sorts the result 
  int_values = numpy.arange(float_values.shape[0]).astype(numpy.int) 

  for idx in range(float_values.shape[0])[::-1]:
    r[r == float_values[idx]] = int_values[idx] 

  return float_values, r

for 循环是一个瓶颈。我需要的应用程序需要很长时间。有没有办法加快/提高其性能？或者是否有一种完全不同但更快的方法来获得我需要的输出？

【问题讨论】：

您真的需要整个索引方案，还是只存储实际距离就可以了。当然你会牺牲一点空间，但谁真正在乎呢？
或者，您是否愿意计算一个中心为零的（对称）999x999 数组，然后对其进行索引？
另外，np.unique 有一个 return_inverse 参数，您可以使用它直接获取索引。但是unique 是一项昂贵的操作，您应该尽量避免这样做。

标签： python performance numpy optimization

【解决方案1】：

这是一个使用masking的矢量化方法-

def func_mask_vectorized(data, (x_0, y_0)):
    # Leverage broadcasting with open meshes to create the squared distances/ids
    m,n = data.shape
    Y,X = np.ogrid[:m,:n]
    ids = (X-x_0)**2 + (Y-y_0)**2

    # Setup mask that will help us retrieve the unique "compressed" IDs
    # (similar to what return_inverse does).
    # This is done by setting 1s at ids places and then using that mask to 
    # assign range covered array, in effect setting up the unique compress. IDs.
    mask = np.zeros(ids.max()+1, dtype=bool)
    mask[ids] = 1    
    id_arr = mask.astype(int)
    id_arr[mask] = np.arange(mask.sum())
    r_out = id_arr[ids]

    # Finally extract out the unique ones among the IDs & get their sqrt values
    float_values_out = np.sqrt(np.flatnonzero(mask))
    return float_values_out, r_out

基准测试

使用形状为(500,500) 的数据对建议设置进行计时，使用问题示例中也使用的0-9 的数字范围，并为下面本节中的所有完整解决方案计时 -

In [371]: np.random.seed(0)
     ...: data = np.random.randint(0,10,(500,500))
     ...: x_0 = 2
     ...: y_0 = 3

# Original soln
In [372]: %timeit func(data, (x_0,y_0))
1 loop, best of 3: 6.77 s per loop

# @Daniel's soln
In [373]: %timeit func_return_inverse(data, (x_0,y_0))
10 loops, best of 3: 23.9 ms per loop

# Soln from this post
In [374]: %timeit func_mask_vectorized(data, (x_0,y_0))
100 loops, best of 3: 5.02 ms per loop

对于数字可能扩展到 100 甚至 1000 的情况进行扩展并不会改变这些数字的叠加方式 -

In [397]: np.random.seed(0)
     ...: data = np.random.randint(0,100,(500,500))
     ...: x_0 = 50
     ...: y_0 = 50

In [398]: %timeit func(data, (x_0,y_0))
     ...: %timeit func_return_inverse(data, (x_0,y_0))
     ...: %timeit func_mask_vectorized(data, (x_0,y_0))
1 loop, best of 3: 5.62 s per loop
10 loops, best of 3: 20.7 ms per loop
100 loops, best of 3: 4.28 ms per loop

In [399]: np.random.seed(0)
     ...: data = np.random.randint(0,1000,(500,500))
     ...: x_0 = 500
     ...: y_0 = 500

In [400]: %timeit func(data, (x_0,y_0))
     ...: %timeit func_return_inverse(data, (x_0,y_0))
     ...: %timeit func_mask_vectorized(data, (x_0,y_0))
1 loop, best of 3: 6.87 s per loop
10 loops, best of 3: 21.9 ms per loop
100 loops, best of 3: 5.05 ms per loop

【讨论】：

一如既往的好@Divakar。但是，您介意同时说明定量性能比较数据吗？说使用 from zmq import Stopwatch; aClk = Stopwatch() 并用所述 构建每个 CodeUnderTest sn-p >data 像这样的大小：aClk.start(); func_mask_vectorised(...); aClk.stop() 这将获得 [us] - 用于比较苹果和苹果的性能的精确计时。这很公平，不是吗？
@user3666197 查看计时部分。希望看起来公平。我使用了 IPython 计时工具 - %timeit，因为它更容易且非常准确。
很棒，一如既往。是的，系统地执行速度提高了约 4 倍[us]-精确的时钟分辨率对于测量比获胜的鼻子距离无关紧要:o) 一切顺利，先生。

【解决方案2】：

不要乱用“唯一距离”数组。只需预先计算由radicand（平方和）索引的距离。这简直就是

roots = [sqrt(float(i)) for i in range(upper_limit)]
然后，由于像素是连续的，您可以选择从参考点向外循环，只需将roots 的整个适用切片从参考点映射到矩阵边缘。

或者，完全退出循环：让numpy 的矢量化操作为您完成，例如

dist = np.sqrt(dist_matrix)

【讨论】：

您是否介意针对原始代码发布您的代码提案的加速因子（详见 Divakar 的矢量化代码下的评论）？这允许将苹果与苹果进行比较，这很公平，不是吗？

【解决方案3】：

使用return_inverse-参数unique：

def func(data, (x_0,y_0)):
    y, x = numpy.indices(data.shape)
    r = (x - x_0)**2 + (y - y_0)**2
    float_values, r = numpy.unique(r, return_inverse=True)
    return float_values ** 0.5, r.reshape(data.shape)

【讨论】：

您是否介意针对原始代码发布您的代码提案的加速因子（详见 Divakar 的矢量化代码下的评论）？这允许将苹果与苹果进行比较，这很公平，不是吗？

【解决方案4】：

您的索引方案（数据中的整数）与距离的顺序相同。如果总是这样，则可以在没有数据实际内容的情况下生成距离数组。

我将此解决方案基于索引计算，该索引计算使用每个位置到锚位置的 x 和 y 像素偏移量。假设“so”是最小偏移量，“ho”是较大的偏移量，“mo”是任一方向上可能的最大偏移量：

index = ho + (mo+1) * lo - lo * (lo+1) // 2

要计算数组中的距离，我们只需要知道矩阵的维度和锚点像素的位置。

import numpy as np
def distanceArray(x,y,cols,rows):
    maxDx  = max(x,cols-x)
    maxDy  = max(y,rows-y)
    maxD   = max(maxDx,maxDy)
    minD   = min(maxDx,maxDy)
    lo = np.arange(minD)[:,None]
    hi = np.arange(maxD)
    sqs = lo*lo + hi*hi
    unique = np.tri(*sqs.shape,maxD-minD, dtype=bool)[::-1,::-1]
    return np.sqrt(sqs[unique])

如果我们只关注距锚点位置的像素偏移量，我们将获得由数据形状边界（maxDx 和 maxDy）确定的水平和垂直 detlas 范围。

对于距离计算，我们可以忽略垂直/水平方向并创建一个小范围和一个大 (r) 范围。（来自 maxD 和 minD 的 lo 和 hi）

要计算所有平方和，我们可以将两个范围之一转置为一个垂直向量 (lo)，然后在对它们的值进行平方 (hi * hi + lo * lo) 后将其添加到另一个 (hi) .这会产生一个包含所有平方和 (sqs) 组合的二维矩阵。

在该矩阵中，顶部三角形是其对应部分的重复。因此，我们使用三角布尔矩阵屏蔽了重复的距离对。（唯一）屏蔽顶部三角形将确保屏蔽操作得出的平方和的顺序是正确的。

最后，过滤后的 sqs 值以正确的顺序包含我们需要的内容。我们可以仅将代价高昂的平方根函数应用于最终结果。

不将距离计算应用于每个像素应该会带来一些显着的性能提升，因为它允许您仅在需要时使用索引距离。我想将这个 distanceArray 函数的性能与其他解决方案进行比较是不公平的（因为它只做了他们所做的一部分）但是，鉴于不必做某事也是优化的一部分，最终结果可能是更好（在我的非科学测试中大约是 Divakar 的 5 倍）。

请注意，如果您只使用一小部分像素的距离，您可能希望避免所有这些计算，并使用字典作为缓存，根据 dX 和 dY 偏移量“按需”计算距离（键控和有序元组）。这将执行绝对最小数量的计算，并且只会为任何特定的偏移对计算一次距离。您甚至可以继续将该缓存用于其他锚点位置和数据形状，因为无论锚点的位置如何，偏移对总是会产生相同的距离。

[EDIT] 要获得与我用于 distanceArray 相同的索引，您可以使用这个：

def offsets(x,y,cols,rows):
    mo   = max(x,cols-x-1,y,rows-y-1)+1

    dx   = abs(np.arange(cols)-x)
    dy   = abs(np.arange(rows)-y)[:,None]

    mo21 = 2 * mo - 1
    ly = dy*(mo21 - dy )//2  # mo*lo - lo*(lo+1)//2 when dy is lowest
    lx = dx*(mo21 - dx )//2  # mo*lo - lo*(lo+1)//2 when dx is lowest

    return np.maximum(dx,dy) + np.minimum(lx,ly)

offsets(3,3,6,6)

array([[9, 8, 6, 3, 6, 8],
       [8, 7, 5, 2, 5, 7],
       [6, 5, 4, 1, 4, 5],
       [3, 2, 1, 0, 1, 2],
       [6, 5, 4, 1, 4, 5],
       [8, 7, 5, 2, 5, 7]])

【讨论】：