python - 如何从一系列点计算从一个点到每个其他点的相对向量

Question

我有一个(x,y)对中的点列表，它代表代理列表的位置。例如，给定 3 个代理，有 3 对点，我将其存储如下:

points = np.array([[x1, y1],
                   [x2, y2],
                   [x3, y3]])

我想计算一个后续数组，即一个代理与其他每个代理的相对位置，但不是其本身。因此，使用上面的数据，我想使用数组 points 生成数组 relative_positions。 point 可以有 N 个位置(我可以同时拥有最多 50-100 个代理)。

因此，使用上述的points，我想生成输出:

relative_positions = [[x2-x1, y2-y1],
                      [x3-x1, y3-y1],
                      [x1-x2, y1-y2],
                      [x3-x2, y3-y2],
                      [x1-x3, y1-y3],
                      [x2-x3, y2-y3]]

例如，给定四个代理位置存储为 numpy 数组:

agent_points = np.array([[10, 1],
                         [30, 3],
                         [25, 10],
                         [5, 5]])

我想生成输出:

relative_positions = [[30-10,  3-1],
                      [25-10, 10-1],
                      [5-10,   5-1],
                      [10-30,  1-3],
                      [25-30, 10-3],
                      [5-30,   5-3],
                      [10-25, 1-10],
                      [30-25, 3-10],
                      [5-25,  5-10],
                      [10-5,   1-5],
                      [30-5,   3-5],
                      [25-5,  10-5]]

我如何有效地做到这一点？我考虑过只计算每个可能的差异，并删除 0 个案例(当它是代理与其自身的相对位置时)，但是我不认为这是一种“纯粹”的方法，因为我可能会意外地删除恰好位于同一点(或非常接近)的代理

Answer 1

方法#1

使用 a 输入数组，您可以执行以下操作 -

d = (a-a[:,None,:])
valid_mask = ~np.eye(len(a),dtype=bool)
out = d[valid_mask]

基本上，我们将 a 扩展到 3D，使第一个轴成为可外部广播，然后我们对其 执行减法>2D 版本，产生 mxmx2 形状的输出，其中 m 为 a.shape[0]。示意性地表示 -

a[:, None, :]    :  4 x 1 x 2
a                :      4 x 2
output           :  4 x 4 x 2

More info .

创建valid_mask的另一种方法是 -

r = np.arange(len(a))
valid_mask = r[:,None] != r

方法#2

我们将利用np.lib.stride_tricks.as_strided来获取3D数组的非对角掩码(沿着前两个轴)，以便我们使用它此处屏蔽差异数组d。此掩码生成的灵感来自发布的 2D 数组问题 here对于 3D 案例看起来像这样 -

def nodiag_view3D(a):
    m = a.shape[0]
    p,q,r = a.strides
    return np.lib.stride_tricks.as_strided(a[:,1:], shape=(m-1,m,2), strides=(p+q,q,r))

要解决我们的问题，那就是 -

d = (a-a[:,None,:])
out = nodiag_view3D(d).reshape(-1,a.shape[1])

<小时/>

展示方法#2 如何改进方法#1 的时机

In [96]: a = np.random.rand(5000,2)

In [97]: d = (a-a[:,None,:])

In [98]: %%timeit
    ...: valid_mask = ~np.eye(len(a),dtype=bool)
    ...: out = d[valid_mask]
1 loop, best of 3: 763 ms per loop

In [99]: %%timeit
    ...: r = np.arange(len(a))
    ...: valid_mask = r[:,None] != r
    ...: out = d[valid_mask]
1 loop, best of 3: 767 ms per loop

In [100]: %timeit nodiag_view3D(d).reshape(-1,a.shape[1])
10 loops, best of 3: 177 ms per loop