Najszybszy sposób obliczenia odległości między poszczególnymi punktami w pythonie

Question

W moim projekcie muszę obliczyć odległość między euclidianami pomiędzy punktami zapisanymi w tablicy. Macierz wejściowa jest tablicą 2D numpy z 3 kolumnami, które są współrzędnymi (x, y, z), a każdy wiersz definiuje nowy punkt.

Zwykle pracuję z 5000 - 6000 punktów w moich testowych przypadkach.

Mój pierwszy algorytm używa Cythona i mojego drugiego numpy. Uważam, że mój numpy algorytm jest szybszy niż cyton.

edit: z 6000 punktów:

NumPy 1,76 s/Cython 4,36 s

Oto mój kod Cython:

cimport cython 
from libc.math cimport sqrt 
@cython.boundscheck(False) 
@cython.wraparound(False) 
cdef void calcul1(double[::1] M,double[::1] R): 

    cdef int i=0 
    cdef int max = M.shape[0] 
    cdef int x,y 
    cdef int start = 1 

    for x in range(0,max,3): 
    for y in range(start,max,3): 

     R[i]= sqrt((M[y] - M[x])**2 + (M[y+1] - M[x+1])**2 + (M[y+2] - M[x+2])**2) 
     i+=1 

    start += 1 

M jest rzutem pamięci początkowej tablicy wejściowej ale flatten() przez numpy przed wywołaniem funkcji calcul1(), R jest widokiem pamięci macierzy wyjściowej 1D do przechowywania wszystkich wyników.

Oto mój kod Numpy:

def calcul2(M): 

    return np.sqrt(((M[:,:,np.newaxis] - M[:,np.newaxis,:])**2).sum(axis=0)) 

Tutaj M jest początkowy tablicy wejście ale transpose() przez numpy przed wywołaniem funkcji mają współrzędne (x, y, z) jako wiersze i punktów jak kolumny.

Ponadto ta funkcja numpy jest dość wygodna, ponieważ tablica, którą zwraca, jest dobrze zorganizowana. Jest to tablica n na n z liczbą punktów, a każdy punkt ma rząd i kolumnę. Tak na przykład odległości AB jest przechowywany w indeksie przecięcia wiersz A i kolumnie B.

Oto jak ja ich (funkcja Cython) ZAPROSZENIE:

cpdef test(): 

    cdef double[::1] Mf 
    cdef double[::1] out = np.empty(17998000,dtype=np.float64) # (6000² - 6000)/2 

    M = np.arange(6000*3,dtype=np.float64).reshape(6000,3) # Example array with 6000 points 
    Mf = M.flatten() #because my cython algorithm need a 1D array 
    Mt = M.transpose() # because my numpy algorithm need coordinates as rows 

    calcul2(Mt) 

    calcul1(Mf,out) 

robię coś źle tutaj? Mój projekt nie jest wystarczająco szybki.

1: Czy istnieje sposób na poprawienie mojego kodu cythonowego w celu pokonania prędkości numpy?

2: Czy istnieje sposób na poprawę mojego numpy kodu, aby jeszcze szybciej obliczyć?

3: Lub inne rozwiązania, ale musi to być python/cython (jak obliczenia równoległe)?

Dziękuję.

Answer 1

Nie wiesz, gdzie otrzymujesz swoje czasy, ale można użyć scipy.spatial.distance:

M = np.arange(6000*3, dtype=np.float64).reshape(6000,3) 
np_result = calcul2(M) 
sp_result = sd.cdist(M.T, M.T) #Scipy usage 
np.allclose(np_result, sp_result) 
>>> True 

taktowanie:

%timeit calcul2(M) 
1000 loops, best of 3: 313 µs per loop 

%timeit sd.cdist(M.T, M.T) 
10000 loops, best of 3: 86.4 µs per loop 

Co ważne, jej również przydatna okazuje się, że wyjście jest symetryczny:

np.allclose(sp_result, sp_result.T) 
>>> True 

Alternatywą jest obliczanie tylko górnego trójkąta tej tablicy:

%timeit sd.pdist(M.T) 
10000 loops, best of 3: 39.1 µs per loop 

Edycja: Nie jesteś pewien, który indeks chcesz zapiąć, wygląda na to, że możesz to robić w obie strony? Zwijanie drugiego indeksu w celu porównania:

%timeit sd.pdist(M) 
10 loops, best of 3: 135 ms per loop 

Jeszcze około 10 razy szybciej niż obecna implementacja NumPy.