Accelerazione distanza tra tutte le possibili coppie in una matrice

Question

Ho un array di x, y, z coordinate diversi (~ 10^10) punti (solo 5 illustrati)

a= [[ 34.45 14.13 2.17] 
    [ 32.38 24.43 23.12] 
    [ 33.19 3.28 39.02] 
    [ 36.34 27.17 31.61] 
    [ 37.81 29.17 29.94]] 

voglio fare una nuova matrice con solo quei punti che sono almeno una certa distanza d lontano da tutti gli altri punti nella lista. Ho scritto un codice utilizzando while ciclo,

import numpy as np 
from scipy.spatial import distance 

d=0.1 #or some distance 
i=0 
selected_points=[] 
while i < len(a): 
      interdist=[] 
      j=i+1 
      while j<len(a): 
       interdist.append(distance.euclidean(a[i],a[j])) 
       j+=1 

      if all(dis >= d for dis in interdist): 
       np.array(selected_points.append(a[i])) 
      i+=1 

Questo funziona, ma sta prendendo molto lungo per eseguire questo calcolo. Ho letto da qualche parte che i loop while sono molto lenti.

Mi chiedevo se qualcuno ha qualche suggerimento su come velocizzare questo calcolo.

EDIT: Mentre il mio obiettivo di trovare le particelle che sono almeno una certa distanza da tutti gli altri rimane la stessa, ho capito che c'è un grave difetto nel mio codice, diciamo che ho 3 particelle, il mio codice fa quanto segue, per la prima iterazione di , calcola le distanze 1->2, 1->3, diciamo che 1->2 è inferiore alla distanza di soglia d, quindi il codice getta via la particella 1. Per la successiva iterazione di i, fa solo 2->3, e diciamo che trova che è maggiore di d, quindi mantiene la particella 2, ma questo è sbagliato! dal 2 dovrebbe anche essere scartato con la particella 1. La soluzione di @svohara è quella giusta!

Answer 1

Per grandi insiemi di dati e bassi dimensioni (come i dati 3-dimensionali), talvolta si trova un grande vantaggio di utilizzare un metodo di indicizzazione spaziale. Una scelta popolare per i dati a bassa dimensionalità è l'albero k-d.

La strategia è indicizzare il set di dati. Quindi interrogare l'indice usando lo stesso set di dati, per restituire i 2 vicini più vicini per ogni punto. Il primo vicino più prossimo è sempre il punto stesso (con dist = 0), quindi vogliamo davvero sapere quanto è lontano il prossimo punto più vicino (secondo vicino più prossimo). Per quei punti in cui il 2-NN è> soglia, hai il risultato.

from scipy.spatial import cKDTree as KDTree 
import numpy as np 

#a is the big data as numpy array N rows by 3 cols 
a = np.random.randn(10**8, 3).astype('float32') 

# This will create the index, prepare to wait... 
# NOTE: took 7 minutes on my mac laptop with 10^8 rand 3-d numbers 
# there are some parameters that could be tweaked for faster indexing, 
# and there are implementations (not in scipy) that can construct 
# the kd-tree using parallel computing strategies (GPUs, e.g.) 
k = KDTree(a) 

#ask for the 2-nearest neighbors by querying the index with the 
# same points 
(dists, idxs) = k.query(a, 2) 
# (dists, idxs) = k.query(a, 2, n_jobs=4) # to use more CPUs on query... 

#Note: 9 minutes for query on my laptop, 2 minutes with n_jobs=6 
# So less than 10 minutes total for 10^8 points. 

# If the second NN is > thresh distance, then there is no other point 
# in the data set closer. 
thresh_d = 0.1 #some threshold, equiv to 'd' in O.P.'s code 
d_slice = dists[:, 1] #distances to second NN for each point 
res = np.flatnonzero(d_slice >= thresh_d) 

Answer 2

1. Eliminare l'append, deve essere molto lento. Puoi avere un vettore statico di distanze e usare [] per mettere il numero nella giusta posizione.

2. Utilizzare min anziché tutti. Hai solo bisogno di controllare se la distanza minima è maggiore di x.

3. In realtà, è possibile interrompere l'aggiunta nel momento in cui si trova una distanza inferiore al limite e quindi è possibile eliminare entrambi i punti. In questo modo non devi nemmeno salvare alcuna distanza (a meno che non ne abbiate bisogno in seguito).

   4. Poiché d (a, b) = d (b, a) è possibile eseguire il ciclo interno solo per i seguenti punti, dimenticare le distanze già calcolate. Se ne hai bisogno, puoi scegliere il più veloce dall'array.

Dal tuo commento, credo che questo avrebbe fatto, se non si hanno punti di ripetute.

selected_points = [] 
for p1 in a: 
    save_point = True 
    for p2 in a: 
     if p1!=p2 and distance.euclidean(p1,p2)<d: 
      save_point = False 
      break 
    if save_point: 
     selected_points.append(p1) 

return selected_points 

Alla fine ho controllare a, b e b, a, perché non si dovrebbe modificare un elenco durante l'elaborazione di esso, ma si può essere più intelligente utilizzando alcune variabili adizionale.

Answer 3

Ecco un approccio vettorializzare utilizzando distance.pdist -

# Store number of pts (number of rows in a) 
m = a.shape[0] 

# Get the first of pairwise indices formed with the pairs of rows from a 
# Simpler version, but a bit slow : idx1,_ = np.triu_indices(m,1) 
shifts_arr = np.zeros(m*(m-1)/2,dtype=int) 
shifts_arr[np.arange(m-1,1,-1).cumsum()] = 1 
idx1 = shifts_arr.cumsum() 

# Get the IDs of pairs of rows that are more than "d" apart and thus select 
# the rest of the rows using a boolean mask created with np.in1d for the 
# entire range of number of rows in a. Index into a to get the selected points. 
selected_pts = a[~np.in1d(np.arange(m),idx1[distance.pdist(a) < d])] 

Per un enorme insieme di dati come 10e10, potremmo avere per eseguire le operazioni in blocchi in base alla memoria di sistema disponibile.

Answer 4

l'algoritmo è quadratico (10^20 operazioni), Ecco un approccio lineare se la distribuzione è quasi casuale. Divide il tuo spazio in scatole di dimensioni d/sqrt(3)^3. Metti ogni punto nella sua scatola.

Poi, per ogni casella,

• se v'è un solo punto, non resta che calcolare la distanza con i punti in un piccolo quartiere.

• altro non c'è niente da fare.