ある配列のどの要素が他の要素の近くにあるかを見つける最も効率的な方法は何ですか？

Question

Iはnumpy.ndarray 2つの1次元のオブジェクトを有し、第一の配列の要素が第二におけるdxの任意要素内にあるワークアウトしたいです。私は現在持っている何

は

10 loops, best of 3: 16.5 msec per loop 

aとbがあることが保証されている場合は特に、find_all_close機能を最適化するための最良の方法何を次のようにtimeit結果を返します

# setup 
numpy.random.seed(1) 
a = numpy.random.random(1000) # create one array 
numpy.random.seed(2) 
b = numpy.random.random(1000) # create second array 
dx = 1e-4 # close-ness parameter 

# function I want to optimise 
def find_all_close(a, b): 
    # compare one number to all elements of b 
    def _is_coincident(t): 
     return (numpy.abs(b - t) <= dx).any() 
    # vectorize and loop over a 
    is_coincident = numpy.vectorize(_is_coincident) 
    return is_coincident(a).nonzero()[0] 

あるfloat配列がfind_all_closeに渡されたときに昇順に並べ替えられ、おそらくcythonなどで並べ替えられますか？

実際には、10,000から100,000要素（またはそれ以上）の配列で作業しており、この操作全体を数百種類の異なる配列bで実行しています。

Answer 1

これを行う最も簡単な方法は、最初の配列の各要素についてです.2つのバイナリ検索で2番目の配列を検索し、最初の配列の要素の上にdx以下、多くてdxを探します。これはlinearithmic時間です：

left = np.searchsorted(b, a - dx, 'left') 
right = np.searchsorted(b, a + dx, 'right') 
a[left != right] 

線形アルゴリズムは、あなたが最初の配列の要素を反復処理するよう移動ウィンドウを追跡する二番目の配列の中に二つのポインタを持っています。

Answer 2

あなたのアプローチは2次です。ソートされた配列の1パス線形時間アルゴリズムがあります。あなたは、2つの配列を適切なタイミングで実行しなければなりません。

def prox(a,b,dx): 
    ia=ib=ir=0 
    res=zeros(a.size,int32) 
    while ia<a.size and ib<b.size: 
     if abs(a[ia]-b[ib])<dx: 
      res[ir]=ia 
      ir += 1 
      ia += 1 
     elif a[ia]>b[ib] : 
       ib += 1 
     else : ia += 1 
    return res[:ir]  

このコードをNumbaでコンパイルすると、パフォーマンスがさらに向上します。

テスト：それは線形時間の複雑さを持っていないので、

a=rand(1000) 
b=rand(1000) 
a.sort() 
b.sort() 

In [10]: prox(a,b,1e-5) 
Out[10]: 
array([ 35, 90, 159, 165, 174, 252, 276, 380, 383, 467, 508, 515, 641, 
     658, 705, 711, 728, 814, 857, 871, 907, 945]) 

In [11]: %timeit prox(a,b,1e-4) 
100 loops, best of 3: 6.23 ms per loop 
In [12]: prox2=numba.jit(prox) 
In [13]: %timeit prox2(a,b,1e-4) 
10000 loops, best of 3: 19.1 µs per loop 

Answer 3

これは、あなたのデータのソート的な性質を利用しない（私は、実行時にはそれは、キャッシュごとにソートされることから利益を得ない疑いが）、しかしnlogn悪いありえない、確かシンプルさとよくtestednessの面で打つのは難しいです：

from scipy.spatial import cKDTree 
print(cKDTree(a[:, None]).query_ball_point(b[:, None], dx)) 

Answer 4

aとbがソートされた配列である場合は、このソート性質が0でを悪用される可能性が。基本的な考え方は、bのインデックスを左に、aの各要素を配置できるようにして、並べ替えられた順序が維持されるようにします。このようにして、bに特定のビンの右側境界を知り、aの各要素を配置することができます。同じビンの左側の境界を取得するには、以前に見つかったインデックスから1を減算するだけです。次に、それぞれのaとの差を左右の境界と比較して、いずれかがしきい値内にあるかどうかを確認し、そうであればそれを有効なインデックスと見なします。

aの要素が最高bよりも大きいと要素がb最低よりも小さいのとき、すなわちとき、コーナーケースのためのいくつかの追加の作業があるでしょう。 np.searchsortedの'left'検索オプションを使用する場合、適切な境界を見つけるために最低でも1にクリップするだけで、同じインデックスをアレイ全体で一度に使用できます。したがって、実装は次のようになります -

def find_all_close_searchsorted(a, b): 
    lidx = np.searchsorted(b,a,'left').clip(min=1,max=b.size-1) 
    close_mask = (np.abs(b[lidx] - a) <= dx) | (np.abs(b[lidx-1] - a) <= dx) 
    return np.nonzero(close_mask)[0] 

ランタイムテスト -

In [2]: np.random.seed(1) 
    ...: a = np.sort(np.random.random(1000)) # create one array 
    ...: np.random.seed(2) 
    ...: b = np.sort(np.random.random(1000)) # create second array 
    ...: dx = 1e-4 # close-ness parameter 
    ...: 

In [3]: np.allclose(find_all_close_searchsorted(a, b),find_all_close(a, b)) 
Out[3]: True 

In [4]: %timeit find_all_close(a,b) 
100 loops, best of 3: 16 ms per loop 

In [5]: %timeit find_all_close_searchsorted(a,b) 
10000 loops, best of 3: 91.4 µs per loop