我有两个包含浮点值的numpy数组x和y.对于x中的每个值,我想在y中找到最接近的元素,而不重用y中的元素.输出应该是1-1的元素索引到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_x,unsorted_y.在numpy / scipy / Python或使用pandas的最佳方法是什么？谢谢.
编辑：澄清我并不是想要找到所有元素的最佳匹配(例如,不是最小化距离的总和),而是最适合每个元素,并且如果它有时以牺牲其他元素为代价也没关系.我假设y通常比x大得多,与上面的例子相反,因此对于y中的x的每个值通常有很多非常好的拟合,我只想有效地找到它.
有人可以为此展示一个scipy kdtrees的例子吗？文档非常稀疏
kdtree = scipy.spatial.cKDTree([x,y])
kdtree.query([-3]*10) # ?? unsure about what query takes as arg
解决方法:
编辑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
对于较小的阵列,它表现稍差.
如果只是为了丢弃重复项,您将不得不遍历所有项目以实现greedy最近邻居算法.考虑到这一点,这是我能够提出的最快速度：
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
