Trouver des voisins avec des coupes efficacement et retour index

Question

J'ai beaucoup de points dans le plan x,y, avec une longueur d'environ 10000, chaque point (x,y) a un rayon intrinsèque r. Ce petit ensemble de données n'est qu'un petit coin de tout mon ensemble de données. J'ai un point intéressé (x1,y1), je veux trouver le point à proximité autour de (x1,y1) dans les 1 et répondre aux critères que la distance entre (x,y) et (x1,y1) est inférieure à r. Je veux retourner l'index de ces bons points, pas les bons points eux-mêmes.

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] 

Dans cette fonction, je peux trouver les bons points autour (x1,y1), mais pas l'indice de ces bons points. TOUTEFOIS, j'ai besoin de l'index ORIGINAL car l'index ORIGINAL est utilisé pour extraire d'autres données associées à la coordonnée (x,y). Comme je l'ai mentionné, l'échantillon de données n'est qu'un petit coin de mon ensemble de données, j'appellerai la fonction ci-dessus environ 1 000 000 fois pour l'ensemble de mes données, donc l'efficacité de la fonction index_finder ci-dessus.

Des réflexions sur une telle tâche?

Answer 1

Approche # 1

On pourrait simplement index dans le premier masque avec son propre masque pour sélectionner les lieux vrais masqués valeurs de la deuxième étape, comme si -

idx[idx] = idx1 

Ainsi, idx aurait les dernières valeurs masquées valides/bonnes places correspondant au tableau original x et y, soit -

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

Ce masque pourrait ensuite être utilisé pour indexer dans d'autres tableaux comme mentionné dans la question.

---

Approche # 2

Comme une autre approche, nous pourrions utiliser deux instructions conditionnelles, créant ainsi deux masques avec eux. Enfin, combinez-les avec AND-ing pour obtenir le masque combiné, qui pourrait être indexé en x et y tableaux pour les sorties finales. Nous n'aurons pas besoin d'obtenir les indices réels de cette façon, donc c'est un avantage de plus.

Par conséquent, la mise en œuvre -

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] 

Si pour une raison quelconque, vous devez toujours les indices correspondants, il suffit de faire -

comb_idx = np.flatnonzero(comb_mask) 

Si vous faites ces opérations pour différentes x1 et y1 paires pour le même ensemble de données x et y, je suggère d'utiliser broadcasting pour le vectoriser à travers tous les x1, y1 paires de données ets, comme indiqué dans this post.

Answer 2

Vous pouvez prendre un masque de vos indices, comme ceci:

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] 

maintenant idx1 est les indices que vous souhaitez extraire.

façon plus rapide en général de le faire est d'utiliser scipy.spatial.KDTree

from scipy.spatial import KDTree 

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

Si vous avez beaucoup de points à la requête sur le même ensemble de données, ce sera beaucoup plus plus vite que d'appeler votre application séquentielle index_finder.

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

AUSSI MAL: si vous avez des distances différentes pour chaque point, vous pouvez faire:

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: Cela devrait fonctionner avec différentes r valeurs pour chaque point

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) 

Ceci est très lent pour une seule paire (x1, y1) car le KDTree prend du temps à remplir. Mais si vous avez des millions de paires, ce sera beaucoup plus rapide.

(je l'ai supposé que vous voulez que l'union de tous les bons points dans les (x, y) données pour tous (x1, y1), si vous voulez séparément, il est également possible en utilisant une méthode similaire, la suppression d'éléments de i[j] selon que d[j] < r[i[j]])

Answer 3

numpy.where semble fait pour trouver les indices

la norme vectorisé calc + np.where() pourrait être plus rapide qu'une boucle

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]]) 

la norme calc peut également inclure le tableau r, la 2ème étape?

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) 

vous pouvez temps le test en 2 étapes vs dont r dans la 1ère norme vectorisé calc?

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)