根据另一个数组中的所有值查找一个数组的最近索引 - Python / NumPy答案

【问题标题】：Find nearest indices for one array against all values in another array - Python / NumPy根据另一个数组中的所有值查找一个数组的最近索引 - Python / NumPy
【发布时间】：2018-01-03 02:39:29
【问题描述】：

我有一个复数列表，我想在另一个复数列表中找到最接近的值。

我目前使用 numpy 的方法：

import numpy as np

refArray = np.random.random(16);
myArray = np.random.random(1000);


def find_nearest(array, value):
    idx = (np.abs(array-value)).argmin()
    return idx;

for value in np.nditer(myArray):
    index = find_nearest(refArray, value);
    print(index);

不幸的是，这需要很长时间才能获得大量值。是否有更快或更“pythonian”的方式将 myArray 中的每个值与 refArray 中最接近的值匹配？

仅供参考：我的脚本中不一定需要 numpy。

重要提示： myArray 和 refArray 的顺序都很重要，不应更改。如果要应用排序，则应以某种方式保留原始索引。

【问题讨论】：

就时间复杂度而言，滑动窗口可能是最有效的。
我看不出滑动窗口比当前解决方案更有效。据我所知，目前的解决方案是 O(n)，这是最好的希望。然后，为了最小化时间常数，需要做一些权衡。但这取决于你的大案子是否会爆炸你的记忆。如果不是内存问题，使用广播可能会有所收获并使用更多 numpy 计算，但如果 RAM 内存是一个问题，这也可能会减慢您的速度。
@JohanL RAM 不是问题，不幸的是时间是问题。这个简单的循环是我能想到的最简单但也是最好的方法。不幸的是，数组大小为 ref=64 和 values=200,000，匹配需要超过 10 秒，目标将远低于 1 秒...跨度>
你可以试试scipy.spatial.KDTree.query。

标签： python arrays list numpy

【解决方案1】：

这是基于this post 的np.searchsorted 的一种矢量化方法-

def closest_argmin(A, B):
    L = B.size
    sidx_B = B.argsort()
    sorted_B = B[sidx_B]
    sorted_idx = np.searchsorted(sorted_B, A)
    sorted_idx[sorted_idx==L] = L-1
    mask = (sorted_idx > 0) & \
    ((np.abs(A - sorted_B[sorted_idx-1]) < np.abs(A - sorted_B[sorted_idx])) )
    return sidx_B[sorted_idx-mask]

简要说明：

获取左侧位置的排序索引。我们使用 - np.searchsorted(arr1, arr2, side='left') 或只是 np.searchsorted(arr1, arr2) 来执行此操作。现在，searchsorted 期望排序数组作为第一个输入，所以我们需要在那里做一些准备工作。
将左侧位置的值与其紧邻右侧位置的值进行比较 (left + 1) 并查看哪个最接近。我们在计算mask 的步骤中执行此操作。
根据左侧或紧邻的右侧是否最接近，选择相应的。这是通过将mask 值作为转换为ints 的偏移量的索引相减来完成的。

基准测试

原始方法-

def org_app(myArray, refArray):
    out1 = np.empty(myArray.size, dtype=int)
    for i, value in enumerate(myArray):
        # find_nearest from posted question
        index = find_nearest(refArray, value)
        out1[i] = index
    return out1

时间和验证 -

In [188]: refArray = np.random.random(16)
     ...: myArray = np.random.random(1000)
     ...: 

In [189]: %timeit org_app(myArray, refArray)
100 loops, best of 3: 1.95 ms per loop

In [190]: %timeit closest_argmin(myArray, refArray)
10000 loops, best of 3: 36.6 µs per loop

In [191]: np.allclose(closest_argmin(myArray, refArray), org_app(myArray, refArray))
Out[191]: True

50x+ 加快了发布的样本，希望更大的数据集更快！

【讨论】：

你真的需要np.abs吗？我想你也可以使用(A - sorted_B[sorted_idx-1] < sorted_B[sorted_idx] - A) 或(A < (sorted_B[sorted_idx-1] + sorted_B[sorted_idx]) / 2)。
另外，我认为使用np.allclose 比较 index 数组没有意义。
这是一个非常有效的解决方案。与stackoverflow.com/a/52799252/8033585 中的cKDTree 进行比较。
谁不直接使用np.sort对B进行排序？
@Anush 因为我在最后一步需要sidx_B。

【解决方案2】：

一个比@Divakar 短得多的答案，同样使用广播，甚至更快：

abs(myArray[:, None] - refArray[None, :]).argmin(axis=-1)

【讨论】：