効率よく切れ目のある近隣を見つけてインデックスを返す

Question

x,y平面に多くの点があり、長さは10000であり、各点(x,y)の固有半径はrです。この小さなデータセットは、私の全データセットのわずか1つです。私は興味のあるポイント(x1,y1)を持っています、(x1,y1)の範囲内で近くのポイントを見つけて、(x,y)と(x1,y1)の間の距離がrより小さいという基準を満たしたいと思います。私は良いポイントそのものではなく、それらの良いポイントのインデックスを返したい。この機能で

import numpy as np 
np.random.seed(2000) 
x = 20.*np.random.rand(10000) 
y = 20.*np.random.rand(10000) 
r = 0.3*np.random.rand(10000) 
x1 = 10. ### (x1,y1) is an interest point 
y1 = 12. 
def index_finder(x,y,r,x1,y1): 
    idx = (abs(x - x1) < 1.) & (abs(y - y1) < 1.) ### This cut will probably cut 90% of the data 
    x_temp = x[idx] ### but if I do like this, then I lose the track of the original index 
    y_temp = y[idx] 
    dis_square = (x_temp - x1)*(x_temp - x1) + (y_temp - y1)*(y_temp - y1) 
    idx1 = dis_square < r*r ### after this cut, there are only a few left 
    x_good = x_temp[idx1] 
    y_good = y_temp[idx1] 

、私はそれらの良い点のインデックス(x1,y1)周りの良い点を見つけることはできませんが。ただし、ORIGINALインデックスは座標(x,y)に関連する他のデータを抽出するために使用されるため、ORIGINALインデックスが必要です。私が言及したように、サンプル・データ・セットは私の全データ・セットのほんの小さなコーナーであり、私は上記の関数を私の全データ・セットに対して約1,000,000回呼び出すので、上記index_finder関数の効率も考慮に入れます。

このようなタスクに関する考えはありますか？我々は、単に真の場所を選択するための独自のマスクを第一のマスクへのインデックスがそうのように、第二段階から値をマスクすることができ

Answer 1

アプローチ＃1

- したがって

idx[idx] = idx1 

、idx元の配列xとyに対応する最終的な有効なマスクされた値/良好な値の場所を持っています。つまり、

x_good = x[idx] 
y_good = y[idx] 

このマスクを使用して、質問に記載されているように他の配列にインデックスを付けることができます。

---

アプローチ＃2

別のアプローチとして、我々はこのように彼らと二つのマスクを作成し、2つの条件文を使用することができます。最後に、AND-ingと組み合わせて組み合わせマスクを取得します。最終的な出力はxとyになります。私たちは実際の指数をそのように得る必要はないので、それがもう一つの利点です。したがって

、実装 -

X = x-x1 
Y = y-y1 
mask1 = (np.abs(X) < 1.) & (np.abs(Y) < 1.) 
mask2 = X**2 + Y*2 < r**2 
comb_mask = mask1 & mask2 

x_good = x[comb_mask] 
y_good = y[comb_mask] 

何らかの理由で、あなたはまだ対応するインデックスが必要な場合は、単に行う -

comb_idx = np.flatnonzero(comb_mask) 

違うx1とy1ペアのためにこれらの操作を行っている場合同じxとyのデータセットの場合、私はbroadcastingを使用して、すべてx1、y1のペアのデータをベクトル化することをお勧めしますets、this postに示すように。

Answer 2

ので、同様にあなたは、あなたのインデックスのマスクを取ることができます：今

def index_finder(x,y,r,x1,y1): 
    idx = np.nonzero((abs(x - x1) < 1.) & (abs(y - y1) < 1.)) #numerical, not boolean 
    mask = (x[idx] - x1)*(x[idx] - x1) + (y[idx] - y1)*(y[idx] - y1) < r*r 
    idx1 = [i[mask] for i in idx] 
    x_good = x_temp[idx1] 
    y_good = y_temp[idx1] 

idx1はあなたが抽出したい指標です。一般的にこれを行うには

より高速な方法は、あなたが同じデータセットに対して照会するための多くのポイントを持っている場合、これは順次ごindex_finderアプリを呼び出すよりも高速ずっとなりますscipy.spatial.KDTree

from scipy.spatial import KDTree 

xy = np.stack((x,y)) 
kdt = KDTree(xy) 
kdt.query_ball_point([x1, y1], r) 

を使用することです。

x1y1 = np.stack((x1, y1)) #`x1` and `y1` are arrays of coordinates. 
kdt.query_ball_point(x1y1, r) 

ALSO WRONG：あなたは、各点に対して異なる距離を持っている場合、あなたが行うことができます：

def query_variable_ball(kdtree, x, y, r): 
    out = [] 
    for x_, y_, r_ in zip(x, y, r): 
     out.append(kdt.query_ball_point([x_, y_], r_) 
    return out 

xy = np.stack((x,y)) 
kdt = KDTree(xy) 
query_variable_ball(kdt, x1, y1, r) 

EDITを2：これは、各ポイント

に異なる r値で動作するはずです

from scipy.spatial import KDTree 

def index_finder_kd(x, y, r, x1, y1): # all arrays 
    xy = np.stack((x,y), axis = -1) 
    x1y1 = np.stack((x1, y1), axis = -1) 
    xytree = KDTree(xy) 
    d, i = xytree.query(x1y1, k = None, distance_upper_bound = 1.) 
    good_idx = np.zeros(x.size, dtype = bool) 
    for idx, dist in zip(i, d): 
     good_idx[idx] |= r[idx] > dist 
    x_good = x[good_idx] 
    y_good = y[good_idx] 
    return x_good, y_good, np.flatnonzero(good_idx) 

非常にはKDTreeのように1つだけ(x1, y1)のペアとしてゆっくりと入力されます。しかし、何百万というペアがあれば、これははるかに高速になります。

（私はあなたがそれらを個別に、それはどうかd[j] < r[i[j]]に基づいてi[j]の要素を削除、同様の方法を使用しても可能ですしたい場合は、すべての(x1, y1)ため(x, y)データ内のすべての良い点の労働組合をしたいと仮定しました）

Answer 3

numpy.whereはインデックス

を見つけるために作られたようです

ベクトル化ノルムカルク+ np.where()ループ

sq_norm = (x - x1)**2 + (y - y1)**2 # no need to take 10000 sqrt 
idcs = np.where(sq_norm < 1.) 

len(idcs[0]) 
Out[193]: 69 

np.stack((idcs[0], x[idcs], y[idcs]), axis=1)[:5] 
Out[194]: 
array([[ 38.  , 9.47165956, 11.94250173], 
     [ 39.  , 9.6966941 , 11.67505453], 
     [ 276.  , 10.68835317, 12.11589316], 
     [ 288.  , 9.93632584, 11.07624915], 
     [ 344.  , 9.48644057, 12.04911857]]) 

より速いかもしれません

ノルムcalcにはr配列も含めることができます。第2ステップですか？

r_sq_norm = (x[idcs] - x1)**2 + (y[idcs] - y1)**2 - r[idcs]**2 
r_idcs = np.where(r_sq_norm < 0.) 

idcs[0][r_idcs] 
Out[11]: array([1575, 3476, 3709], dtype=int64) 

あなたが時間に第一ベクトル化ノルムCALCでrを含めた対2段階のテストをしたいのでしょうか？

sq_norm = (x - x1)**2 + (y - y1)**2 - r**2 
idcs = np.where(sq_norm < 0.) 

idcs[0] 
Out[13]: array([1575, 3476, 3709], dtype=int64)