Finden Nachbarn mit Schnitten effizient und Index

Question

Ich habe viele Punkte in der x,y Ebene, mit der Länge um 10000, hat jeder Punkt (x,y) einen Eigenradius r. Dieser kleine Datensatz ist nur eine winzige Ecke meines gesamten Datensatzes. Ich habe einen interessiert Punkt (x1,y1), Ich möchte in der Nähe von Punkt (x1,y1) innerhalb 1 finden und erfüllen die Kriterien, dass der Abstand zwischen (x,y) und (x1,y1) ist weniger als r. Ich möchte den Index dieser guten Punkte zurückgeben, nicht die guten Punkte selbst.

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] 

In dieser Funktion kann ich die guten Punkte um (x1,y1), finden aber nicht den Index dieser guten Punkte. JEDOCH brauche ich den ORIGINAL-Index, weil der ORIGINAL-Index verwendet wird, um andere Daten zu extrahieren, die der Koordinate (x,y) zugeordnet sind. Wie ich bereits erwähnt habe, ist der Beispieldatensatz nur eine winzige Ecke meines gesamten Datensatzes. Ich werde die obige Funktion etwa 1.000.000 Mal für meinen gesamten Datensatz aufrufen, daher ist auch die Effizienz der obigen index_finder Funktion eine Überlegung.

Irgendwelche Gedanken zu einer solchen Aufgabe?

Answer 1

Ansatz # 1

Wir kann einfach Index in die erste Maske mit einer eigenen Maske für die wahren Orte maskierten Werte aus der zweiten Stufe der Auswahl, wie so -

idx[idx] = idx1 

So idx hätte die endgültigen gültigen maskierten Werte/Orte mit guten Werten, die dem ursprünglichen Array x und y entsprechen, dh -

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

Diese Maske könnte dann verwendet werden, um in andere Arrays zu indizieren, wie in der Frage erwähnt.

---

Ansatz # 2

Als weiterer Ansatz, könnten wir zwei bedingte Anweisungen verwenden, wodurch zwei Masken mit ihnen. Schließlich kombinieren Sie sie mit AND-ing, um die kombinierte Maske zu erhalten, die in x und y Arrays für die endgültigen Ausgaben indiziert werden konnte. Wir müssen die tatsächlichen Indizes nicht auf diese Weise erhalten, das ist ein weiterer Vorteil.

Daher ist die Umsetzung -

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] 

aus irgendeinem Grund Wenn man noch die entsprechenden Indizes benötigen, tun nur -

comb_idx = np.flatnonzero(comb_mask) 

Wenn Sie diese Operationen für verschiedene x1 und y1 Paare tun für den gleichen x und y Datensatz, würde ich empfehlen, broadcasting zu verwenden, um es durch all diese x1, y1 gepaarten Daten zu vektorisieren ets, wie in this post gezeigt.

Answer 2

Sie können eine Maske Ihrer Indizes nehmen, etwa so:

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] 

jetzt idx1 die Indizes, die Sie extrahieren möchten.

Schneller Weg im Allgemeinen, dies zu tun ist scipy.spatial.KDTree

from scipy.spatial import KDTree 

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

verwenden Wenn Sie viele Punkte gegen den gleichen Datenbestand abzufragen, wird diese viel schneller als sequentiell Ihre index_finder App aufrufen.

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

auch falsch:, wenn Sie für jeden Punkt unterschiedliche Abstände haben, können Sie tun:

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: Diese mit unterschiedlichen r Werte für jeden Punkt arbeiten sollte

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) 

Dies ist sehr langsam für nur eine (x1, y1) Paar als die KDTree dauert eine Weile, um zu bevölkern. Aber wenn Sie Millionen von Paaren haben, wird dies viel schneller sein.

(Ich habe angenommen, Sie wollen die Vereinigung aller guten Punkte in den (x, y) Daten für alle (x1, y1), wenn man sie separat wollen ist es auch möglich, ein ähnliches Verfahren verwendet wird, Elemente von i[j] Entfernen basierend darauf, ob d[j] < r[i[j]])

Answer 3

numpy.where scheint für die Suche nach den Indizes

die vektorisiert Norm calc + np.where() könnte schneller als eine Schleife

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

gemacht

die Norm calc kann das r Array auch enthalten, der 2. Schritt?

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) 

Sie können den 2-Stufen-Test vs einschließlich r in der 1. vektorisiert Norm berechnet und wollen?

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)