5

给定两个长度相等的数组，一个保存数据，一个保存结果但最初设置为零，例如：

a = numpy.array([1, 0, 0, 1, 0, 1, 0, 0, 1, 1])
b = numpy.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

我想计算 a 中三个相邻元素的所有可能子集的总和。如果和为0或1，则b中对应的三个元素保持不变；只有当总和超过 1 时，b 中对应的三个元素才设为 1，这样计算后 b 变为

array([0, 0, 0, 1, 1, 1, 0, 1, 1, 1])

一个简单的循环将完成此操作：

for x in range(len(a)-2):
    if a[x:x+3].sum() > 1:
        b[x:x+3] = 1

在此之后， b 具有所需的形式。

我必须为大量数据执行此操作，因此速度是一个问题。NumPy 中是否有更快的方法来执行上述操作？

（我知道这类似于卷积，但不完全相同）。

4

3 回答 3

6

您可以从卷积开始，选择超过 1 的值，最后使用“膨胀”：

b = numpy.convolve(a, [1, 1, 1], mode="same") > 1
b = b | numpy.r_[0, b[:-1]] | numpy.r_[b[1:], 0]

由于这避免了 Python 循环，它应该比你的方法更快，但我没有做计时。

另一种方法是使用第二个卷积来扩张：

kernel = [1, 1, 1]
b = numpy.convolve(a, kernel, mode="same") > 1
b = numpy.convolve(b, kernel, mode="same") > 0

如果您有可用的 SciPy，则另一个扩张选项是

b = numpy.convolve(a, [1, 1, 1], mode="same") > 1
b = scipy.ndimage.morphology.binary_dilation(b)

编辑：通过做一些计时，我发现这个解决方案对于大型阵列来说似乎是最快的：

b = numpy.convolve(a, kernel) > 1
b[:-1] |= b[1:]  # Shift and "smearing" to the *left* (smearing with b[1:] |= b[:-1] does not work)
b[:-1] |= b[1:]  # … and again!
b = b[:-2]

对于一百万个条目的数组，它比您在我的机器上使用的原始方法快 200 多倍。正如 EOL 在评论中指出的那样，这个解决方案可能被认为有点脆弱，因为它取决于 NumPy 的实现细节。

于 2012-04-02T11:46:23.993 回答
2

您可以使用以下有效方式计算“卷积”总和：

>>> a0 = a[:-2]
>>> a1 = a[1:-1]
>>> a2 = a[2:]
>>> a_large_sum = a0 + a1 + a2 > 1

b然后可以通过编写意味着“三个相邻值中的至少一个为真”的内容来有效地完成更新a_large_sum：您首先将a_large_sum数组扩展回与a（向右，向左和向右，然后向左）：

>>> a_large_sum_0 = np.hstack([a_large_sum, [False, False]])
>>> a_large_sum_1 = np.hstack([[False], a_large_sum, [False]])
>>> a_large_sum_2 = np.hstack([[False, False], a_large_sum])

然后，您b以有效的方式获得：

>>> b = a_large_sum_0 | a_large_sum_1 | a_large_sum_2

通过利用 NumPy 内部快速循环，这给出了您获得的结果，但以一种非常有效的方式。

PS：这种方法本质上与 Sven 的第一个解决方案相同，但比 Sven 的优雅代码更普通；然而，它同样快。Sven 的第二个解决方案 (double convolve()) 更加优雅，速度提高了一倍。

于 2012-04-02T12:02:31.317 回答
1

您可能还想看看 NumPy 的stride_tricks. 使用 Sven 的时序设置（参见 Sven 的答案中的链接），我发现对于（非常）大的数组，这也是一种快速的方法来做你想做的事情（即你的定义a）：

shape = (len(a)-2,3)
strides = a.strides+a.strides
a_strided = numpy.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)
b = np.r_[numpy.sum(a_strided, axis=-1) > 1, False, False]
b[2:] |= b[1:-1] | b[:-2]

编辑后（见下面的评论）它不再是最快的方式。

这会在您的原始数组上创建一个特别跨步的视图。里面的数据a没有被复制，只是以一种新的方式查看。我们基本上想要创建一个新数组，其中最后一个索引包含我们想要求和的子数组（即您想要求和的三个元素）。这样，我们可以很容易地用最后一个命令求和。

因此，这个新形状的最后一个元素必须是3，第一个元素将是旧的a负 2 的长度（因为我们只能求和到第-2nd 个元素）。

strides 列表包含新数组a_strided到达形状每个维度中的下一个元素所需的步幅（以字节为单位）。如果您将它们设置为相等，则意味着它们a_strided[0,1]都a_strided[1,0]将是a[1]，这正是我们想要的。在普通数组中，情况并非如此（第一个步幅将是“第一维的大小乘以数组第一维的长度（= shape[0]）”），但在这种情况下，我们可以好好利用它。

不确定我是否真的很好地解释了这一切，但只需打印出 a_strided ，您就会看到结果是什么以及这使操作变得多么容易。

于 2012-04-03T20:22:25.620 回答