Condensé matrice de distance à distance complète de la matrice

une matrice de distance condensée telle que retournée par le pdist peut être convertie en une matrice de distance complète en utilisant scipy.spatial.distance.squareform :

>>> import numpy as np
>>> from scipy.spatial.distance import pdist, squareform
>>> points = np.array([[0,1],[1,1],[3,5], [15, 5]])
>>> dist_condensed = pdist(points)
>>> dist_condensed
array([  1.        ,   5.        ,  15.5241747 ,   4.47213595,
        14.56021978,  12.        ])

utiliser squareform pour convertir en matrice complète:

>>> dist = squareform(dist_condensed)
array([[  0.        ,   1.        ,   5.        ,  15.5241747 ],
       [  1.        ,   0.        ,   4.47213595,  14.56021978],
       [  5.        ,   4.47213595,   0.        ,  12.        ],
       [ 15.5241747 ,  14.56021978,  12.        ,   0.        ]])

la Distance entre le point I, j est stockée dans dist [i, j]:

>>> dist[2, 0]
5.0
>>> np.linalg.norm(points[2] - points[0])
5.0

Indices à l'indice condensé

On peut convertir les indices utilisés pour accéder aux éléments de la matrice carrée en indice dans la matrice condensée:

def square_to_condensed(i, j, n):
    assert i != j, "no diagonal elements in condensed matrix"
    if i < j:
        i, j = j, i
    return n*j - j*(j+1)/2 + i - 1 - j

exemple:

>>> square_to_condensed(1, 2, len(points))
3
>>> dist_condensed[3]
4.4721359549995796
>>> dist[1,2]
4.4721359549995796

indice condensé des indices

aussi l'autre direction est possible sans sqaureform qui est meilleur en termes de durée d'exécution et la consommation de mémoire:

import math

def calc_row_idx(k, n):
    return int(math.ceil((1/2.) * (- (-8*k + 4 *n**2 -4*n - 7)**0.5 + 2*n -1) - 1))

def elem_in_i_rows(i, n):
    return i * (n - 1 - i) + (i*(i + 1))/2

def calc_col_idx(k, i, n):
    return int(n - elem_in_i_rows(i + 1, n) + k)

def condensed_to_square(k, n):
    i = calc_row_idx(k, n)
    j = calc_col_idx(k, i, n)
    return i, j

exemple:

>>> condensed_to_square(3, 4)
(1.0, 2.0)

Runtime comparaison avec la forme carrée

>>> import numpy as np
>>> from scipy.spatial.distance import pdist, squareform
>>> points = np.random.random((10**4,3))
>>> %timeit dist_condensed = pdist(points)
1 loops, best of 3: 555 ms per loop

créer la sqaureform s'avère être vraiment lent:

>>> dist_condensed = pdist(points)
>>> %timeit dist = squareform(dist_condensed)
1 loops, best of 3: 2.25 s per loop

si nous recherchons deux points avec une distance maximale, il n'est pas surprenant que la recherche dans la matrice complète soit O(n) alors que sous forme condensée seulement O(n/2):

>>> dist = squareform(dist_condensed)
>>> %timeit dist_condensed.argmax()
10 loops, best of 3: 45.2 ms per loop
>>> %timeit dist.argmax()
10 loops, best of 3: 93.3 ms per loop

obtenir les inideces pour les deux points prend presque pas de temps dans les deux cas, mais bien sûr, il ya un peu de frais généraux pour le calcul de la indice condensé:

>>> idx_flat = dist.argmax()
>>> idx_condensed = dist.argmax()
>>> %timeit list(np.unravel_index(idx_flat, dist.shape))
100000 loops, best of 3: 2.28 µs per loop
>>> %timeit condensed_to_square(idx_condensed, len(points))
100000 loops, best of 3: 14.7 µs per loop

le vecteur de la matrice comprimée correspond à la région triangulaire inférieure de la matrice carrée. Pour convertir un point dans cette région triangulaire, vous devez calculer le nombre de points à la gauche dans le triangle, et le nombre ci-dessus dans la colonne.

, Vous pouvez utiliser la fonction suivante pour convertir:

q = lambda i,j,n: n*j - j*(j+1)/2 + i - 1 - j

Case:

import numpy as np
from scipy.spatial.distance import pdist, squareform
x = np.random.uniform( size = 100 ).reshape( ( 50, 2 ) )
d = pdist( x )
ds = squareform( d )
for i in xrange( 1, 50 ):
    for j in xrange( i ):
        assert ds[ i, j ] == d[ q( i, j, 50 ) ]

C'est le triangle supérieur version ( i < j ), qui doit être intéressant pour certains:

condensed_idx = lambda i,j,n: i*n + j - i*(i+1)/2 - i - 1

C'est très facile à comprendre:

1. avec i*n + j vous allez à la position dans la matrice carrée;
2. avec - i*(i+1)/2 vous supprimez le triangle inférieur (y compris la diagonale) dans toutes les lignes avant i;
3. avec - i vous supprimez les positions de la ligne i avant la diagonale;
4. avec - 1 vous supprimez les positions de la ligne i sur la diagonale.

Case:

import scipy
from scipy.spatial.distance import pdist, squareform
condensed_idx = lambda i,j,n: i*n + j - i*(i+1)/2 - i - 1
n = 50
dim = 2
x = scipy.random.uniform(size = n*dim).reshape((n, dim))
d = pdist(x)
ds = squareform(d)
for i in xrange(1, n-1):
    for j in xrange(i+1, n):
        assert ds[i, j] == d[condensed_idx(i, j, n)]

si vous voulez accéder à l'élément pdist correspondant à l'élément (i,j)-th de la matrice de la distance carrée, le calcul est le suivant: supposons i < j (autrement des indices de retournement) si i == j , la réponse est 0.

X = random((N,m))
dist_matrix = pdist(X)

puis l'élément (i, j)-th est dist_matrix [ind] où

ind = (N - array(range(1,i+1))).sum()  + (j - 1 - i).

si quelqu'un est à la recherche d'une transformation inverse (c.-à-d. si on lui donne un indice d'élément idx , comprendre quel élément (i, j) correspond à lui), voici une solution à vecteur résonnant:

def actual_indices(idx, n):
    n_row_elems = np.cumsum(np.arange(1, n)[::-1])
    ii = (n_row_elems[:, None] - 1 < idx[None, :]).sum(axis=0)
    shifts = np.concatenate([[0], n_row_elems])
    jj = np.arange(1, n)[ii] + idx - shifts[ii]
    return ii, jj

n = 5
k = 10
idx = np.random.randint(0, n, k)
a = pdist(np.random.rand(n, n))
b = squareform(a)

ii, jj = actual_indices(idx, n)]
assert np.allclose(b[ii, jj, a[idx])

Je l'ai utilisé pour calculer les indices des lignes les plus proches dans une matrice.

m = 3  # how many closest
lowest_dist_idx = np.argpartition(-a, -m)[-m:]
ii, jj = actual_indices(lowest_dist_idx, n)  # rows ii[0] and jj[0] are closest