在Python中寻找两个列表/数组中最近的项

10

我有两个包含浮点数值的numpy数组 xy。对于每个x中的值,我想找到y中最接近的元素,不重复使用y中的元素。输出应该是x元素索引到y元素索引的一一映射。以下是一个不好的做法,它依赖于排序。它会从列表中删除配对的每个元素。不排序会很糟糕,因为匹配将取决于原始输入数组的顺序。

def min_i(values):
    min_index, min_value = min(enumerate(values),
                               key=operator.itemgetter(1))
    return min_index, min_value

# unsorted elements
unsorted_x = randn(10)*10
unsorted_y = randn(10)*10

# sort lists
x = sort(unsorted_x)
y = sort(unsorted_y)

pairs = []
indx_to_search = range(len(y))

for x_indx, x_item in enumerate(x):
    if len(indx_to_search) == 0:
        print "ran out of items to match..."
        break
    # until match is found look for closest item
    possible_values = y[indx_to_search]
    nearest_indx, nearest_item = min_i(possible_values)
    orig_indx = indx_to_search[nearest_indx]
    # remove it
    indx_to_search.remove(orig_indx)
    pairs.append((x_indx, orig_indx))
print "paired items: "
for k,v in pairs:
    print x[k], " paired with ", y[v]

我更喜欢不先对元素进行排序,但如果它们已经排序,则想要获取原始未排序列表unsorted_xunsorted_y中的索引。在numpy/scipy/Python或使用pandas中,最好的方法是什么?谢谢。

编辑:为了澄清,我不是试图找到所有元素的最佳拟合(例如最小化距离总和),而是每个元素的最佳拟合,即使有时会牺牲其他元素也可以。我假设y通常比上述示例中的x大得多,因此对于xy中的每个值,通常有许多非常好的拟合,并且我只想高效地找到其中一个。

有人能展示一下用scipy的kdtrees的例子吗?文档非常稀少。

kdtree = scipy.spatial.cKDTree([x,y])
kdtree.query([-3]*10) # ?? unsure about what query takes as arg

“(5, 10), (6, 0)” 的预期答案是什么?“(10, 5), (6, 0)” 的预期答案是什么? - Robᵩ
我建议使用 scipy.spatial.cKDTree(或旧版 scipy 中的 KDTree)。当然,如果你小心处理,argsorting 也可以实现。 - seberg
@Robᵩ:不应该这样,我错了,但我已经编辑了我的答案。基本上,你可以假设对于y中的每个x值都有多个良好的匹配,但你当然是正确的,会有这些权衡。 - user248237
1
@Jaime,不太确定你的意思,你可以使用它获取查询集外点的k个最近邻居。tree = KDTree(x[:,None]); tree.query(y[:,None], k=1)会找到所有y中最接近的x(基于二次范数,你可以进行更改)。 - seberg
显示剩余9条评论
1个回答

9

编辑2 如果您能选择一定数量的邻居,以保证每个数组项都有唯一的邻居,则使用KDTree解决方案可以表现得非常好。以下是代码:

def nearest_neighbors_kd_tree(x, y, k) :
    x, y = map(np.asarray, (x, y))
    tree =scipy.spatial.cKDTree(y[:, None])    
    ordered_neighbors = tree.query(x[:, None], k)[1]
    nearest_neighbor = np.empty((len(x),), dtype=np.intp)
    nearest_neighbor.fill(-1)
    used_y = set()
    for j, neigh_j in enumerate(ordered_neighbors) :
        for k in neigh_j :
            if k not in used_y :
                nearest_neighbor[j] = k
                used_y.add(k)
                break
    return nearest_neighbor

对于 n=1000 个点的样本,我得到:

In [9]: np.any(nearest_neighbors_kd_tree(x, y, 12) == -1)
Out[9]: True

In [10]: np.any(nearest_neighbors_kd_tree(x, y, 13) == -1)
Out[10]: False

所以最佳值为k=13,然后时间为:
In [11]: %timeit nearest_neighbors_kd_tree(x, y, 13)
100 loops, best of 3: 9.26 ms per loop

但最坏情况下,您可能需要 k=1000,然后:

In [12]: %timeit nearest_neighbors_kd_tree(x, y, 1000)
1 loops, best of 3: 424 ms per loop

哪个比其他选项慢:

In [13]: %timeit nearest_neighbors(x, y)
10 loops, best of 3: 60 ms per loop

In [14]: %timeit nearest_neighbors_sorted(x, y)
10 loops, best of 3: 47.4 ms per loop

编辑 在搜索之前对数组进行排序,适用于超过1000个项目的数组:

def nearest_neighbors_sorted(x, y) :
    x, y = map(np.asarray, (x, y))
    y_idx = np.argsort(y)
    y = y[y_idx]
    nearest_neighbor = np.empty((len(x),), dtype=np.intp)
    for j, xj in enumerate(x) :
        idx = np.searchsorted(y, xj)
        if idx == len(y) or idx != 0 and y[idx] - xj > xj - y[idx-1] :
            idx -= 1
        nearest_neighbor[j] = y_idx[idx]
        y = np.delete(y, idx)
        y_idx = np.delete(y_idx, idx)
    return nearest_neighbor

有一个长度为10000的数组:

In [2]: %timeit nearest_neighbors_sorted(x, y)
1 loops, best of 3: 557 ms per loop

In [3]: %timeit nearest_neighbors(x, y)
1 loops, best of 3: 1.53 s per loop

对于较小的数组,它的性能稍差。


你需要遍历所有项目来实现贪心最近邻算法,即使只是为了丢弃重复项。考虑到这一点,以下是我能够想出的最快方法:
def nearest_neighbors(x, y) :
    x, y = map(np.asarray, (x, y))
    y = y.copy()
    y_idx = np.arange(len(y))
    nearest_neighbor = np.empty((len(x),), dtype=np.intp)
    for j, xj in enumerate(x) :
        idx = np.argmin(np.abs(y - xj))
        nearest_neighbor[j] = y_idx[idx]
        y = np.delete(y, idx)
        y_idx = np.delete(y_idx, idx)

    return nearest_neighbor

现在使用:

n = 1000
x = np.random.rand(n)
y = np.random.rand(2*n)

我得到:
In [11]: %timeit nearest_neighbors(x, y)
10 loops, best of 3: 52.4 ms per loop

有没有一种方法可以使用cKDTree来避免重复项?即使稍微降低性能也可以。 - user248237
另一个问题:有没有办法确保 p.argmin(np.abs(y - xj)) 会忽略像 NaN 这样的缺失值?它是否会在某些情况下选择这些值? - user248237
你需要的是 np.nanargmin - denis
这种方法是否也适用于多维点?因为我总是会收到错误提示:在“tree =scipy.spatial.cKDTree(y[:, None])”这一行中,缓冲区的维数不正确(应该是2,但得到了3)。 - Varlor
1
@jaime 2022年所有3个返回错误: --------------------------------------------------nearest_neighbors_kd_tree: 在scipy.spatial.ckdtree.cKDTree.init()中的ckdtree.pyx中ValueError: 数据必须是2维的--------------------- nearest_neighbors(x, y): ValueError: 形状为(127,) (2,) 的操作数无法广播在一起 ---------------------nearest_neighbors_sorted(x,y): ValueError: 对象对于所需的数组来说太深了 - openSourcerer
1
好的,从KDTree解决方案中删除“[:, None]”后可以无错误运行。然而,没有办法避免重复项而不运行整个过程。在完整数据集上,我遇到了内存问题MemoryError:无法为形状为(144008,144008)且数据类型为float64的数组分配155. GiB的内存。 - openSourcerer

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接