10

2 つの numpy 配列がxありy、float 値が含まれています。の各値について、 の要素を再利用せずxに、 の最も近い要素を見つけたいと思います。出力は、x の要素のインデックスから y の要素のインデックスへの 1-1 マッピングである必要があります。これは、ソートに依存する悪い方法です。リストからペアになった各要素を削除します。ペアリングは元の入力配列の順序に依存するため、並べ替えがないと、これはうまくいきません。yy

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上記の例とは対照的に、それは一般的にはるかに大きいと想定しているため、 inxの各値には通常多くの非常に適切な適合があり、それを効率的に見つけたいだけです。xy

誰かがこれのために scipy kdtrees の例を示すことができますか? ドキュメントはかなりまばらです

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

1 に答える 1

9

EDIT 2KDTree配列内のすべてのアイテムに対して一意の隣人を持つことを保証するいくつかの隣人を選択できる場合、使用するソリューションは非常にうまく機能します。次のコードを使用します。

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
于 2013-03-12T16:08:55.023 に答える