python - Numpy:考虑项目的邻居及其在数组中的位置的快速计算

Question

我有 4 个二维 numpy 数组，称为 a、b、c、d，每个数组由 n 行和 m 列组成.我需要做的是给 b 和 d 的每个元素一个计算如下的值(伪代码):

min_coords = min_of_neighbors_coords(x, y)
b[x,y] = a[x,y] * a[min_coords];
d[x,y] = c[min_coords];

其中 min_of_neighbors_coords 是一个函数，它在给定数组元素的坐标的情况下返回具有较低值的“邻居”元素的坐标。即，考虑数组:

1, 2, 5
3, 7, 2
2, 3, 6

min_of_neighbors_coords(1, 1) 将引用值为 7 的中心元素，并返回元组 (0, 0):数字1的坐标。

我设法使用 for 循环(每个元素一个元素)来做到这一点，但该算法非常慢，我正在寻找一种改进它的方法，避免循环并要求对 numpy 进行计算。

这可能吗？

Answer 1

编辑 我将原来的答案保留在底部。正如 Paul 在评论中指出的那样，原始答案并没有真正回答 OP 的问题，使用 ndimage 过滤器可以更轻松地实现。以下更繁琐的功能应该做正确的事情。它需要两个数组，a 和 c，并返回 a 的窗口最小值和 c 中的值a 中窗口最小值的位置:

def neighbor_min(a, c):
    ac = np.concatenate((a[None], c[None]))
    rows, cols = ac.shape[1:]
    ret = np.empty_like(ac)

    # Fill in the center
    win_ac = as_strided(ac, shape=(2, rows-2, cols, 3),
                        strides=ac.strides+ac.strides[1:2])
    win_ac = win_ac[np.ogrid[:2, :rows-2, :cols] +
                    [np.argmin(win_ac[0], axis=2)]]
    win_ac = as_strided(win_ac, shape=(2, rows-2, cols-2, 3),
                        strides=win_ac.strides+win_ac.strides[2:3])
    ret[:, 1:-1, 1:-1] =  win_ac[np.ogrid[:2, :rows-2, :cols-2] +
                                 [np.argmin(win_ac[0], axis=2)]]

    # Fill the top, bottom, left and right borders
    win_ac = as_strided(ac[:, :2, :], shape=(2, 2, cols-2, 3),
                        strides=ac.strides+ac.strides[2:3])
    win_ac = win_ac[np.ogrid[:2, :2, :cols-2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, 0, 1:-1] = win_ac[:, np.argmin(win_ac[0], axis=0),
                             np.ogrid[:cols-2]]
    win_ac = as_strided(ac[:, -2:, :], shape=(2, 2, cols-2, 3),
                        strides=ac.strides+ac.strides[2:3])
    win_ac = win_ac[np.ogrid[:2, :2, :cols-2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, -1, 1:-1] = win_ac[:, np.argmin(win_ac[0], axis=0),
                             np.ogrid[:cols-2]]
    win_ac = as_strided(ac[:, :, :2], shape=(2, rows-2, 2, 3),
                        strides=ac.strides+ac.strides[1:2])
    win_ac = win_ac[np.ogrid[:2, :rows-2, :2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, 1:-1, 0] = win_ac[:, np.ogrid[:rows-2],
                             np.argmin(win_ac[0], axis=1)]
    win_ac = as_strided(ac[:, :, -2:], shape=(2, rows-2, 2, 3),
                        strides=ac.strides+ac.strides[1:2])
    win_ac = win_ac[np.ogrid[:2, :rows-2, :2] +
                    [np.argmin(win_ac[0], axis=2)]]
    ret[:, 1:-1, -1] = win_ac[:, np.ogrid[:rows-2],
                             np.argmin(win_ac[0], axis=1)]
    # Fill the corners
    win_ac = ac[:, :2, :2]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, 0, 0] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
    win_ac = ac[:, :2, -2:]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, 0, -1] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
    win_ac = ac[:, -2:, -2:]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, -1, -1] = win_ac[:, np.argmin(win_ac[0], axis=-1)]
    win_ac = ac[:, -2:, :2]
    win_ac = win_ac[:, np.ogrid[:2],
                    np.argmin(win_ac[0], axis=-1)]
    ret[:, -1, 0] = win_ac[:, np.argmin(win_ac[0], axis=-1)]

    return ret

返回的是一个 (2, rows, cols) 数组，可以解包成两个数组:

>>> a = np.random.randint(100, size=(5,5))
>>> c = np.random.randint(100, size=(5,5))
>>> a
array([[42, 54, 18, 88, 26],
       [80, 65, 83, 31,  4],
       [51, 52, 18, 88, 52],
       [ 1, 70,  5,  0, 89],
       [47, 34, 27, 67, 68]])
>>> c
array([[94, 94, 29,  6, 76],
       [81, 47, 67, 21, 26],
       [44, 92, 20, 32, 90],
       [81, 25, 32, 68, 25],
       [49, 43, 71, 79, 77]])
>>> neighbor_min(a, c)
array([[[42, 18, 18,  4,  4],
        [42, 18, 18,  4,  4],
        [ 1,  1,  0,  0,  0],
        [ 1,  1,  0,  0,  0],
        [ 1,  1,  0,  0,  0]],

       [[94, 29, 29, 26, 26],
        [94, 29, 29, 26, 26],
        [81, 81, 68, 68, 68],
        [81, 81, 68, 68, 68],
        [81, 81, 68, 68, 68]]])

OP 的案例可以这样解决:

def bd_from_ac(a, c):
    b,d = neighbor_min(a, c)
    return a*b, d

虽然性能受到严重影响，但速度仍然非常快:

In [3]: a = np.random.rand(1000, 1000)

In [4]: c = np.random.rand(1000, 1000)

In [5]: %timeit bd_from_ac(a, c)
1 loops, best of 3: 570 ms per loop

---

您实际上并没有将最小相邻元素的坐标用于获取它以外的任何其他用途，因此您不妨跳过该部分并创建一个 min_neighbor 函数。如果您不想求助于 cython 进行快速循环，您将不得不使用滚动窗口 View ，例如 Paul 的链接中概述的那样。这通常会将您的 (m, n) 数组转换为相同数据的 (m-2, n-2, 3, 3) View ，然后您将在最后两个轴上应用 np.min。

不幸的是，您必须一次在一个轴上应用它，因此您将必须创建数据的 (m-2, n-2, 3) 副本。幸运的是，您可以分两步计算最小值，首先沿一个轴开窗和最小化，然后沿另一个轴开窗和最小化，并获得相同的结果。所以最多你将拥有输入大小的中间存储。如果需要，您甚至可以将输出数组重新用作中间存储并避免内存分配，但那是 left as exercise ...

下面的函数就是这样做的。它有点冗长，因为它不仅要处理中心区域，还要处理四个边缘和四个角的特殊情况。除此之外，它是一个非常紧凑的实现:

def neighbor_min(a):
    rows, cols = a.shape
    ret = np.empty_like(a)

    # Fill in the center
    win_a = as_strided(a, shape=(m-2, n, 3),
                       strides=a.strides+a.strides[:1])
    win_a = win_a.min(axis=2)
    win_a = as_strided(win_a, shape=(m-2, n-2, 3),
                       strides=win_a.strides+win_a.strides[1:])
    ret[1:-1, 1:-1] = win_a.min(axis=2)

    # Fill the top, bottom, left and right borders
    win_a = as_strided(a[:2, :], shape=(2, cols-2, 3),
                       strides=a.strides+a.strides[1:])
    ret[0, 1:-1] = win_a.min(axis=2).min(axis=0)
    win_a = as_strided(a[-2:, :], shape=(2, cols-2, 3),
                       strides=a.strides+a.strides[1:])
    ret[-1, 1:-1] = win_a.min(axis=2).min(axis=0)
    win_a = as_strided(a[:, :2], shape=(rows-2, 2, 3),
                       strides=a.strides+a.strides[:1])
    ret[1:-1, 0] = win_a.min(axis=2).min(axis=1)
    win_a = as_strided(a[:, -2:], shape=(rows-2, 2, 3),
                       strides=a.strides+a.strides[:1])
    ret[1:-1, -1] = win_a.min(axis=2).min(axis=1)

    # Fill the corners
    ret[0, 0] = a[:2, :2].min()
    ret[0, -1] = a[:2, -2:].min()
    ret[-1, -1] = a[-2:, -2:].min()
    ret[-1, 0] = a[-2:, :2].min()

    return ret

您现在可以执行以下操作:

>>> a = np.random.randint(10, size=(5, 5))
>>> a
array([[0, 3, 1, 8, 9],
       [7, 2, 7, 5, 7],
       [4, 2, 6, 1, 9],
       [2, 8, 1, 2, 3],
       [7, 7, 6, 8, 0]])
>>> neighbor_min(a)
array([[0, 0, 1, 1, 5],
       [0, 0, 1, 1, 1],
       [2, 1, 1, 1, 1],
       [2, 1, 1, 0, 0],
       [2, 1, 1, 0, 0]])

你原来的问题可以这样解决:

def bd_from_ac(a, c):
    return a*neighbor_min(a), neighbor_min(c)

作为性能基准:

In [2]: m, n = 1000, 1000

In [3]: a = np.random.rand(m, n)

In [4]: c = np.random.rand(m, n)

In [5]: %timeit bd_from_ac(a, c)
1 loops, best of 3: 123 ms per loop