基于元素距离交错两个numy数组（python）

Question

我有两个数组，看起来像这样：

X = np.array([ 157,  262,  368,  472,  577,  682,  786,  891,  996, 1100, 1204,
       1310, 1415, 1520, 1625, 1731, 1879])

Y = np.array([  30,  135,  240,  345,  450,  555,  660,  765,  870,  975, 1080,
       1185, 1290, 1395, 1500, 1605])

阵列将：

• 将值从开始按升序排序。在
• 有时长度不等。在

我想根据以下内容将这两个元素交叉放入一个新数组Z：

• 每个元素只能使用一次
• 不需要使用所有元素
• 元素Xi只能包含在Z中，前提是Yj在Y中没有其他元素的值差小于abs(Xi - Yj)，并且{}中没有到{}的值距离小于{}的元素。{cd5>中的元素适用于相同的元素

我知道我可以用一堆嵌套的for循环来实现这一点，但我想知道是否有更聪明、更整洁的方法来实现这一点？在

（我意识到，从我提出问题的方式来看，这听起来像是从教科书上剪下来的。不是的。但也许这是一个经典的排序函数，谁知道呢，但对我这个生物学家来说。。。我只能说，我不知道如何高效、简洁地解决问题。）

编辑：不是很好的例子

^{pr2}$

我甚至不完全确定它是否满足所有的标准，但是当重写旧代码时，我只做了一个循环。。。一定还有更好的办法！在

Answer 1

让我们试着为这个例子找出解决方案：

In [1]: import numpy as np

In [5]: X = np.array([1879, 1731])

In [6]: Y = np.array([1481, 1691, 1586, 1796])

我们可以计算X中的值与Y中的值之间的所有距离，如下所示：

^{pr2}$

行对应于X值，列对应于Y值。在

为了找到最接近Y中某个元素的X值，我们将查找 对于X，它对应于dist的列中的最小值 矩阵。每列对应一个特定的Y，因此 列对应于某个X和某个{}之间的最小值。在

从视觉上讲，我们要寻找的是dist中的值，它们是 两者的最小值它们所在的行，和它们所在的列 在。我们称之为“行列最小值”。在

在上面的dist数组中，40是行列的最小值。65是最小值， 但不是行列最小值。在

对于每一列，我们可以找到X索引，该索引通过以下方式将列最小化：

In [6]: idx1 = np.argmin(dist, axis = 0)

In [7]: idx1
Out[7]: array([1, 1, 1, 1])

同样，对于每一行，我们可以这样找到Y索引：

In [8]: idx2 = np.argmin(dist, axis = 1)

In [9]: idx2
Out[9]: array([3, 1])

现在，让我们暂时忘掉这个例子，假设idx1如下所示：

        0,1,2,3,4,5   # the index value 
idx1 = (_,_,_,_,_,2,...)

这意味着在第5列中，第2行的值最小。在

然后，如果第2行、第5列对应于行列最小值，则idx2 必须看起来像这样：

        0,1,2        # index value
idx2 = (_,_,5,...)

我们可以用NumPy来表达这种关系

idx1[idx2] == np.arange(len(X))
idx2[idx1] == np.arange(len(Y))    

所以对应于行列最小值的X，Y值是

X[idx1[idx2] == np.arange(len(X))]

以及

Y[idx2[idx1] == np.arange(len(Y))]

---

import numpy as np
tests = [
    (np.array([1879, 1731]),
     np.array([1481, 1691, 1586, 1806])), 
    (np.array([1879, 1731]),
     np.array([1481, 1691, 1586, 1796])),
    (np.array([ 157,  262,  368,  472,  577,  682,  786,  891,  996, 1100, 1204]),
     np.array([  30,  135,  240,  345,  450,  555,  660,  765,  870,  975])),
    (np.array([ 157, 262, 368, 472, 577, 682, 786, 891, 996, 1100, 1204, 1310,
                1415, 1520, 1625, 1731, 1879]),
     np.array([ 221, 326, 431, 536, 641, 746, 851, 956, 1061, 1166, 1271, 1376,
                1481, 1586, 1691, 1796]))]

def find_close(X,Y):
    new_list = list()
    for i in X:
        delta_i = np.abs(Y - i)
        # print(delta_i)
        delta_reciprocal = np.abs(X - Y[delta_i.argmin()])
        if delta_i.min() == delta_reciprocal.min():
            new_list += sorted([Y[delta_i.argmin()],
                                X[delta_reciprocal.argmin()]])
    Z = np.array(new_list)
    return Z

def alt_find_close(X,Y):
    dist = np.abs(np.subtract.outer(X,Y))
    idx1 = np.argmin(dist, axis = 0)
    idx2 = np.argmin(dist, axis = 1)
    Z = np.r_[X[idx1[idx2] == np.arange(len(X))], Y[idx2[idx1] == np.arange(len(Y))]]
    return Z

for X, Y in tests:
    assert np.allclose(sorted(find_close(X,Y)), sorted(alt_find_close(X,Y)))

---

Timeit结果：

% python -mtimeit -s'import test' 'test.find_close(test.X, test.Y)'
1000 loops, best of 3: 454 usec per loop
% python -mtimeit -s'import test' 'test.alt_find_close(test.X, test.Y)'
10000 loops, best of 3: 40.6 usec per loop

因此alt_find_close明显比find_close快。在

Answer 2

我想你可能想使用scipy.spatial.cKDTree（你可以自己用numpy来构建这样的东西searchsorted，但我看不出有多大意义，除非可以对等距离的问题有更多的控制）。在

但是一般来说你应该小心。您的例子是整数数组，根据发现的情况，可能会有问题（argmin总是找到第一个，所以可能并不重要，但是对于这个例子，如果您有相等的距离，您可能会放弃点）。在

import numpy as np
from scipy.spatial import cKDTree

def find_close_fast(X, Y):
    kX = cKDTree(X[:,None]) # needs to be 2D
    kY = cKDTree(Y[:,None])

    nearest_X = kX.query(Y[:,None], p=1)[1] # might as well use 1-norm

    # Which Y corresponds the other way around?
    nearest_Y = kY.query(X[nearest_X][:,None], p=1)[1]

    w = nearest_Y == np.arange(len(Y))
    result = np.concatenate((X[nearest_X[w]], Y[w]))
    return result

如果数组变大（每个数组可能有几百个左右），这将大大加快。例如：

^{pr2}$

但是我在这里使用浮点数是有原因的，如果你能有相等的距离，结果就不能保证。分类是不同的，没有真正找出原因。在

编辑：实际上，您还可以将两个数组排序为一个数组（并记住哪个数组属于哪个类），然后检查两个不同类的相邻位置。如果一个点有另一个类的两个邻接点，用手挑近一个。这可能更快，而且只使用numpy。在