时间戳

首页 > 解决方案 > 使用 numpy 快速查找子矩阵

python - 使用 numpy 快速查找子矩阵

问题描述

我正在为 Python编写 Bridson 的泊松圆盘采样 ( https://www.cs.ubc.ca/~rbridson/docs/bridson-siggraph07-poissondisk.pdf ) 的实现。

这种抽样的一个特点是样本之间有一个设定的最小距离，它避免了聚类。新候选点的邻居查找利用了此功能，并使用背景网格来加快搜索速度。

背景网格由M存储True非空单元格和False其他单元格的布尔网格和P存储精确坐标的点网格组成。它们都实现为 n 维numpy数组。

网格大小cellsize的选择方式是每个网格单元中最多有一个样本，然后您只需检查几个最近的行和列。

现在，我使用以下过程来检查该点p是否靠近任何现有点：

    def in_neighborhood(p, n=2):
        indices = (p/cellsize).astype(int)
        indmin = np.maximum(indices-n, np.zeros(ndim, dtype=int))
        indmax = np.minimum(indices+n+1, gridsize)
        if M[tuple(indices)]:
            return True
        for offset in np.ndindex(tuple(indmax - indmin)):
            indcur = tuple([sum(x) for x in zip(indmin,offset)])
            if M[indcur] and squared_distance(p, P[indcur]) < squared_radius:
                return True
        return False

这占用了大部分的执行时间（> 90%）。我怀疑我在for这里实现循环的方式确实不是最理想的，但是当我事先不知道维数时，我真的不知道如何提取子矩阵。任何帮助表示赞赏。

完整的代码和带有一些使用示例的 Jupyter 笔记本位于https://github.com/diregoblin/poisson_disc_sampling

标签： pythonnumpyoptimization

解决方案

切片是这样的：

维度使它有点棘手。最后，将切片准备为数组的开销非常值得。如果您事先知道尺寸，也许可以进行进一步的优化。应该差别不大。作为一个有趣的部分，这使得 M 矩阵几乎是多余的 - 它仅用于检查点本身是否值得尝试，因为切片的其余部分大大加快了乘法的速度。

更新：

我根据评论重新引入了 M 检查，并且还使用了操作np.square而不是np.power提到的操作。

np.power > np.square 转换后的时间：

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
     1    0.196    0.196    5.906    5.906 (Bridson_sampling)

更新后的功能：

def in_neighborhood(p, n=2):
    indices = (p / cellsize).astype(int)
    indmin = np.maximum(indices - n, np.zeros(ndim, dtype=int))
    indmax = np.minimum(indices + n + 1, gridsize)
    if M[tuple(indices)]:
        return True
    a = []
    for i in range(ndim):
        a.append(slice(indmin[i], indmax[i]))
    if np.any(np.logical_and(M[tuple(a)], np.sum(np.square(p - P[tuple(a)]), axis=ndim) < squared_radius)):
        return True

结果似乎是一样的，但是我仍然建议您检查一下——我可能会出错。结果时间快了约 5 倍。仍然不是那么快，但是我认为我现在可以在 6.8 秒内接受约 120K 的调用。

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.207    0.207    6.802    6.802 (Bridson_sampling)
   123689    3.163    0.000    4.818    0.000 (in_neighborhood)
507833/507831    0.224    0.000    2.220    0.000 implement_array_function}
   503091    0.439    0.000    1.610    0.000 (_wrapreduction)
   142170    0.310    0.000    1.469    0.000 (in_limits)

旧部分。

为了完整性和娱乐性，我将留下其余的答案。从失败中学习。

我花了最后几个小时挖掘您的代码无济于事。我设法减少了大约 3% 的执行时间，这可能是分析器的错误，我不高兴。在这个过程中，我对一些事情感到有些困惑，所以我想问问它们，也许它可以帮助一点？老实说，我不知道。

在您在 github 中拥有的算法中，在生成球形样本时，您似乎并没有排除内球体。这可能会为您的算法添加一些不必要的迭代。这是故意的吗？

引用源材料：

当活动列表不为空时，从中选择一个随机索引（比如 i）。从围绕 xi 的半径 r 和 2r 之间的球面环中均匀选择最多 k 个点。对于每个点，依次检查它是否在现有样本的距离 r 内（使用背景网格仅测试附近的样本）。如果一个点距离现有样本足够远，则将其作为下一个样本发出并将其添加到活动列表中。如果在 k 次尝试后没有找到这样的点，则将 i 从活动列表中删除。

另一件事：据我了解，作者建议：

  1. 用一个点初始化列表。
  2. 如果非空，则从列表中随机取一个点
  3. 设置 k = 0
  4. 取一个随机邻居点作为样本。
  5. 如果该点是合适的候选者，则将其添加到列表中
  6. 增加 k
  7. 如果 k = 30 并且没有找到点，则从列表中删除原始点。
  8. 回到 1

这种解释源于：

3 分析步骤 2 精确执行 2N-1 次以产生 N 个样本：每次迭代要么产生一个新样本并将其添加到活动列表中，要么从活动列表中删除一个现有样本。步骤 2 的每次迭代都需要 O(k) 时间，并且由于 k 保持不变（通常非常小），因此算法是线性的。

作者指出，在每次迭代中，您要么添加一个点，要么删除一个点。他没有说一次添加多个点。（现在我想起来了，我不确定这个 2N 如何与他可以不止一次碰到一个点的事实相适应 - 他会的。总的来说，我对算法描述感到困惑。也许这只是不过我累了。）

你所做的略有不同，因为你总是取完整的 30 个样本，全部尝试，然后删除原始元素。从输出来看，它似乎与作者声称的结果相当相似。

另外，为了清楚起见，N 是什么？我假设它是要生成的点数。

  • 结果中生成了约 5K 点
  • 该in_neighborhood函数被调用约 122K 次
  • 这完全在作者声称的范围内 (2N * k)

一个观察：许多计算都是多余的，因为 P[indcur] 通常为零 - 至少在开始时，因为我没有足够的耐心等待以后的打印。所以也许在计算过程中的某种缓存可能会有所帮助。然而，正如我刚刚测试的那样，元组转换相当快，所以维度很可能不是我一开始想的问题。

另一个可能想到的事情是许多浮点运算可能是可能的原因之一。

推荐阅读