2013-03-12 14 views
7

float値を含むnumpy配列が2つあります。xyです。 xの各値について、yの要素を再利用することなく、最も近い要素をyに配置したいと考えています。出力は、xからyまでの要素のインデックスまでの要素のインデックスの1-1のマッピングでなければなりません。ソートに頼ることは悪いことです。これは、リストからペアになった各要素を削除します。ソートしなければ、ペアリングは元の入力配列の順序に依存するため、これは悪いことになります。Pythonの2つのリスト/配列に最も近いアイテムを見つける

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でこれを行う最善の方法は何ですか?ありがとう。

編集:私はすべてのelemets全体でベストフィット(例えば距離の和を最小化することではない)のではなく、各要素のための最良のフィットを見つけようとしていないよ、それは費用で時々だ場合、それは大丈夫だと明確にします他の要素の。私はyが上記の例とは逆に一般にxよりもかなり大きいと仮定しているので、xの各値には通常多くの良い適合があり、yにあります。

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

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

インデックスを検索するためのバイナリ検索の並べ替えはおそらく最善の策だと思います。 – mgilson

+0

@mgilton:scipy/numpyにバイナリサーチアルゴスが組み込まれていますか? – user248237dfsf

+0

はい:[numpy.searchsorted](http://docs.scipy.org/doc/numpy/reference/generated/numpy.searchsorted.html) – mgilson

答えて

6

EDIT 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 

他のオプションより遅いです:

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付:検索は1000個の以上のアイテムの配列のために報わ前に、配列をソート

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 

EDIT要素長配列:

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 
+0

ありがとうございます。 'cKDTree'を使って重複することなくそれを行う方法はありますか?わずかなパフォーマンスでもヒット? – user248237dfsf

+0

別の質問: 'p.argmin(np.abs(y-xj))'がNaNのような欠損値を無視するようにする方法はありますか?それがそれを選ぶケースはいつですか? – user248237dfsf

+0

[np.nanargmin](http://docs.scipy.org/doc/numpy/reference/generated/numpy.nanargmin.html)はあなたが望むものです。 – denis

-1

このくらい単純化されたコード完全にうまくいった。

N=12 
M=15 

X = [np.random.random() for i in range(N)] 
Y = [np.random.random() for i in range(M)] 

pair = [] 

for x in X: 
    t = [abs(x-y) for y in Y] 
    ind = t.index(min(t)) 
    pair.append((x,Y[ind])) 
    X.remove(x) 
    Y.remove(Y[ind]) 

print(pair) 
+1

これは悪い考えです。まず、コードを反復しながらXから要素を取り除くので、コードは機能しません!さらに、元のポスターのすべての説明を本当に読みましたか?あなたは本当に彼/彼女の完全な質問に答えているようではありません。 –

関連する問題