根据距离优化点网格（python）

Question

我正面临一个问题，我在微积分时间和约束之间陷入困境......

1. 简单案例：

我在 3D 网格上有点均匀分布的点网格（第一张图片）。

注意：您在那里看不到它，但有一个垂直于网格的第三轴 - 有 125 个点（每边 5 x 5 x 5）

3D grid with dots distributed ~ equally

我的目标是，采用噪声统计数据给出的距离（在这种情况下我们将采用 0.03），用这些条件填充网格：

• 点是随机生成的，并通过距离条件测试发送，这将决定它是否保留它
• 距离条件：如果一个点与其他点的距离为 0.03 +/- 10%，则可以将其放入网格中，否则将其移除。

这是一个填充网格的示例：

Badly filled grid (still in 3D)

问题是，这里有很多我想删除的空白。

更多细节，这里是从一个点到它最近的邻居的所有距离图：

Histogram of distances from dots their closest neighbour

红色矩形是我希望所有蓝色条堆叠的位置，在这种情况下，介于 0.03 - 10% 和 0.03 + 10% 之间。

主要问题是，因为它是随机生成的，所以我必须生成很多点才能希望即使只有一个也能符合约束条件。这是一个相对“简单的案例”，因为之后，我将不得不将其应用于：

1. 硬盒：

   • 点分布不均成超过10个轴
   • 点数最多 1000 个或更多

我不是 python 专家，我只知道基础知识（numpy、matplotlib.pyplot），并且对图表、数组、循环非常熟悉，但是我的代码在计算时间方面并未针对这些进行优化任务类型：

• 生成一个点
• 它与点的距离是针对每个其他点计算的（我不知道还能做些什么来防止扫描整个列表..）
• 如果一个点距离前一个点更近，则成为“最近的一个”
• 如果与“最近的”的距离在约束中，则不会删除。

如果有人有一个简单的想法或一些建议，那就太好了，我也希望我能解释清楚（我尽量只保留最重要的东西，并尝试用正确的英文写:)）。

感谢阅读，

内森。

编辑：这里简化了代码，因为它可能有助于更轻松地回答！ 我直接添加了非均布网格（A，B）。

import numpy as np
import matplotlib.pyplot as plt
plt.close("all")

A = np.random.rand(100)
B = np.random.rand(100)

plt.figure()
plt.plot(A, B,".")
plt.grid()

Nfillers = 10000
SNRD = 0.03

for i in range(Nfillers):
    a, b = np.random.rand(2)
    A = np.hstack((A, np.atleast_1d(a)))
    B = np.hstack((B, np.atleast_1d(b)))
    DCN = 1E+308 # Distance Closest Neighbor
    for j in range(len(A)): # len(A) is updated with each iteration
        D1 = A[j] - a
        D2 = B[j] - a
        distance = np.sqrt(D1**2 + D2**2)
        if distance != 0:
            if distance < DCN:
                DCN = distance
    if DCN < 0.9*SNRD or DCN > 1.1*SNRD:
        A = np.delete(A, len(A)-1)
        B = np.delete(B, len(B)-1)

    print i

plt.figure()
plt.plot(A, B, ".")
plt.grid()

编辑 2：更快的代码：

import numpy as np
import matplotlib.pyplot as plt

A = np.random.rand(100)
B = np.random.rand(100)

plt.figure()
plt.plot(A, B,".")
plt.grid()

Nfillers = 100000
SNRD = 0.03

for i in range(Nfillers):
    a, b = np.random.rand(2)
    DCN = np.min(np.sqrt((A - a)**2 + (B - b)**2))
    if DCN > 0.9*SNRD and DCN < 1.1*SNRD:
        A = np.hstack((A, np.atleast_1d(a)))
        B = np.hstack((B, np.atleast_1d(b)))

    print(i)

plt.figure()
plt.plot(A, B, ".")
plt.grid()

Answer 1

您没有包含任何代码，但这里有一个关于如何加速这类问题（计算多点之间的距离）的粗略想法。

将空间划分为大致所需距离大小的块，d。这可能是块编号到点列表的映射，例如blocks[i,j,k] = list_of_points。那么当你生成一个新的点时，你只需要检查相邻的块。比如：

new_point = random_point()
i, j, k = get_block_indices(new_point)

for ai in (i-1, i, i+1):
    for aj in (j-1, j, j+1):
        for ak in (k-1, k, k+1):
            for other_point in blocks[ai, aj, ak]:
                # check whether other_point is too close
                found_close_other_point = check(new_point, other_point, d)
                if found_close_other_point:
                    return  # start again

# only reach if we haven't found a close point
try:
    blocks[i, j, k].append(new_point)
except KeyError:
    # create new block
    blocks[i, j, k] = [new_point]

在您想要存储多少块和需要检查多少个邻居之间进行权衡。补充几点：

• 您不想存储每个块，而是在该块尚不存在 (KeyError) 时将其处理为空。在最坏的情况下，如果你像这样懒惰地这样做，你将需要与点一样多的块。

• 您可以使用itertools.product 来避免嵌套的 for 循环。

• 如果您使用块的大小，您还可以引入额外的跳过条件，例如如果您知道您可能无法在一个块中容纳超过n 个点，那么您只需检查len(block[i, j, k]) > n。

• 纯 python 实现可能不是最快的，你可以使用http://numba.pydata.org/ 来加快速度，这对于这类数字问题来说非常容易和快速，同时仍然像 python 一样阅读。

---

编辑：一个实现。

import itertools
import numpy as np
from math import ceil
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


def compute_min_distances(ps):
    """Compute the minimum distances between all points ``ps``,
    with shape (num_points, num_dims).
    """
    ndim = ps.shape[-1]
    p1, p2 = ps.reshape(-1, 1, ndim), ps.reshape(1, -1, ndim)
    dd = np.sum((p1 - p2)**2, -1)**0.5
    # want to ignore distance to self
    np.fill_diagonal(dd, np.inf)
    return np.min(dd, axis=0)


def get_block_indices(p, B):
    """Given B blocks, get the indices of ``p``.
    """
    return [int(c * B) for c in p]


def init_blocks(ps, B):
    """Set up the initial block mapping.
    """
    blocks = {}
    for p in ps:
        i, j, k = get_block_indices(p, B)
        try:
            blocks[i, j, k].append(p)
        except KeyError:
            blocks[i, j, k] = [p]
    return blocks


def distance(p1, p2):
    """Distance bewtween two points.
    """
    return sum((a - b)**2 for a, b in zip(p1, p2))**0.5


def gen_close_points(p0s, D=0.03, N=1000):
    """Takes ``p0s`` and generates ``N`` new points with each of 
    which has a closest neighbour +- 10% of ``D``.
    """
    # number of blocks required so that only neigbours matter
    B = ceil(1 / D)

    # mapping of blocks to points
    blocks = init_blocks(p0s, B)

    n = 0
    while n < N:
        new_p = np.random.rand(3)
        i, j, k = get_block_indices(new_p, B)

        neighbour_inds = [(b - 1, b, b + 1) for b in (i, j, k)]

        def gen_neighbours():
            for ai, aj, ak in itertools.product(*neighbour_inds):
                try:
                    # try yielding each point
                    yield from blocks[ai, aj, ak]
                except KeyError:
                    # skip if empty
                    continue

        ds = (distance(new_p, p) for p in gen_neighbours())
        min_d = min(ds, default=0)

        if 0.9 * D < min_d < 1.1 * D:
            try:
                blocks[i, j, k].append(new_p)
            except KeyError:
                blocks[i, j, k] = [new_p]
            n += 1

    # combine into array of all points
    return np.concatenate(tuple(itertools.chain(blocks.values())))

让我们用一些初始点来试一试：

p0s = np.random.rand(100, 3)
d0s = compute_min_distances(p0s)
plt.hist(d);

现在让我们再生成 10000 个点，每个点都指定最近距离：

ps = gen_close_points(p0s, 0.03, 10000)

并检查距离：

plt.hist(compute_min_distances(ps));

看起来不错！最后让我们将旧点和新点绘制在一起：

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter3D(ps[:, 0], ps[:, 1], ps[:, 2], alpha=0.1)
ax.scatter3D(p0s[:, 0], p0s[:, 1], p0s[:, 2])

目前假设数据是范围为 (0, 1) 的 3 维数据，但应该是对其进行概括的一个不错的起点。