Cálculo eficiente de una matriz de distancia euclidiana usando numpy

Tengo un conjunto de puntos en el espacio bidimensional y necesito calcular la distancia de cada punto a cada otro punto.

Tengo un número relativamente pequeño de puntos, tal vez a lo sumo 100. Pero ya que tengo que hacerlo con frecuencia y rapidez para determinar las relaciones entre estos puntos en movimiento, y como soy consciente de que la iteración a través de los puntos podría ser tan mala como complejidad O (n ^ 2), estoy buscando formas de aprovechar la magia de matriz de Numpy (o scipy).

Tal como está en mi código, las coordenadas de cada objeto se almacenan en su clase. Sin embargo, también podría actualizarlos en una matriz numpy cuando actualizo la coordenada de clase.

class Cell(object): """Represents one object in the field.""" def __init__(self,id,x=0,y=0): self.m_id = id self.m_x = x self.m_y = y 

Se me ocurre crear una matriz de distancia euclidiana para evitar la duplicación, pero quizás tenga una estructura de datos más inteligente.

También estoy abierto a punteros a algoritmos ingeniosos.

Además, observo que hay preguntas similares relacionadas con la distancia y el número euclidianos, pero no encontramos ninguna que aborde directamente esta cuestión de poblar de manera eficiente una matriz de distancia completa.

Puedes aprovechar el tipo complex :

 # build a complex array of your cells z = np.array([complex(c.m_x, c.m_y) for c in cells]) 

Primera solucion

 # mesh this array so that you will have all combinations m, n = np.meshgrid(z, z) # get the distance via the norm out = abs(mn) 

Segunda solucion

La malla es la idea principal. Pero numpy es inteligente, por lo que no tiene que generar m & n . Simplemente calcule la diferencia utilizando una versión transpuesta de z . La malla se realiza automáticamente:

 out = abs(z[..., np.newaxis] - z) 

Tercera solucion

Y si z se establece directamente como una matriz bidimensional, puede usar zT lugar de la extraña z[..., np.newaxis] . Así que finalmente, su código se verá así:

 z = np.array([[complex(c.m_x, c.m_y) for c in cells]]) # notice the [[ ... ]] out = abs(zT-z) 

Ejemplo

 >>> z = np.array([[0.+0.j, 2.+1.j, -1.+4.j]]) >>> abs(zT-z) array([[ 0. , 2.23606798, 4.12310563], [ 2.23606798, 0. , 4.24264069], [ 4.12310563, 4.24264069, 0. ]]) 

Como complemento, es posible que desee eliminar los duplicados posteriormente, tomando el triángulo superior:

 >>> np.triu(out) array([[ 0. , 2.23606798, 4.12310563], [ 0. , 0. , 4.24264069], [ 0. , 0. , 0. ]]) 

Algunos puntos de referencia

 >>> timeit.timeit('abs(zT-z)', setup='import numpy as np;z = np.array([[0.+0.j, 2.+1.j, -1.+4.j]])') 4.645645342274779 >>> timeit.timeit('abs(z[..., np.newaxis] - z)', setup='import numpy as np;z = np.array([0.+0.j, 2.+1.j, -1.+4.j])') 5.049334864854522 >>> timeit.timeit('m, n = np.meshgrid(z, z); abs(mn)', setup='import numpy as np;z = np.array([0.+0.j, 2.+1.j, -1.+4.j])') 22.489568296184686 

Aquí es cómo puedes hacerlo usando numpy:

 import numpy as np x = np.array([0,1,2]) y = np.array([2,4,6]) # take advantage of broadcasting, to make a 2dim array of diffs dx = x[..., np.newaxis] - x[np.newaxis, ...] dy = y[..., np.newaxis] - y[np.newaxis, ...] dx => array([[ 0, -1, -2], [ 1, 0, -1], [ 2, 1, 0]]) # stack in one array, to speed up calculations d = np.array([dx,dy]) d.shape => (2, 3, 3) 

Ahora todo lo que queda es calcular la norma L2 a lo largo del eje 0 (como se explica aquí ):

 (d**2).sum(axis=0)**0.5 => array([[ 0. , 2.23606798, 4.47213595], [ 2.23606798, 0. , 2.23606798], [ 4.47213595, 2.23606798, 0. ]]) 

Si no necesita la matriz de distancia completa, será mejor que use kd-tree. Considere scipy.spatial.cKDTree o sklearn.neighbors.KDTree . Esto se debe a que un kd-árbol kan encuentra k-vecinos más cercanos en tiempo O (n log n), y por lo tanto, evita la complejidad O (n ** 2) de calcular todas las distancias n por n.

Jake Vanderplas da este ejemplo utilizando la transmisión en Python Data Science Handbook , que es muy similar a lo que propuso @ shx2.

 import numpy as np rand = random.RandomState(42) X = rand.rand(3, 2) dist_sq = np.sum((X[:, np.newaxis, :] - X[np.newaxis, :, :]) ** 2, axis = -1) dist_sq array([[0. , 0.18543317, 0.81602495], [0.18543317, 0. , 0.22819282], [0.81602495, 0.22819282, 0. ]])