spectral_embedding: corrected bug in the computation of the Laplacian

+ minor changes (API, etc.)

Author: Emmanuelle Gouillart

diffusions.py

@@ -170,6 +170,34 @@ def _make_weights_3d(edges, data, beta=130, eps=1.e-6):
     weights = np.exp(- beta*gradients / (10*data.std())) + eps
     return weights
 
+def _make_distances_3d(edges, data):
+    lx, ly, lz = data.shape
+    gradients = np.abs(data[edges[0]/(ly*lz), \
+                                 (edges[0] % (ly*lz))/lz, \
+                                 (edges[0] % (ly*lz))%lz] - \
+                            data[edges[1]/(ly*lz), \
+                                 (edges[1] % (ly*lz))/lz, \
+                                 (edges[1] % (ly*lz)) % lz])
+    return gradients
+
+def _make_adaptive_weights(edges, data):
+    print "adaptive"
+    pixel_nb = len(np.unique(edges.ravel()))
+    gradients = _make_distances_3d(edges, data)
+    i_indices = np.hstack((edges[0], edges[1]))
+    j_indices = np.hstack((edges[1], edges[0]))
+    w = np.hstack((gradients, gradients))
+    nb = np.bincount(i_indices).astype(np.float)
+    total_weight = np.bincount(i_indices, weights=w)
+    sigmas = total_weight / nb
+    sigma_of_edges = np.array([sigmas[edges[0]], sigmas[edges[1]]])
+    return _make_weights_adaptative(gradients, sigma_of_edges)
+
+def _make_weights_adaptative(gradients, sigma_of_edges, eps=1.e-10):
+    sigma_i, sigma_j = sigma_of_edges
+    weights = np.exp(- gradients**2 / (sigma_i * sigma_j)) + eps
+    return weights 
+
 def _make_laplacian_sparse(edges, weights):
     """
     Sparse implementation
@@ -190,19 +218,19 @@ def _make_normed_laplacian(edges, weights):
     """
     Sparse implementation
     """
-    tol = 1.e-8
-    eps = 1.e-5
+    eps = 0
     pixel_nb = len(np.unique(edges.ravel()))
     diag = np.arange(pixel_nb)
     i_indices = np.hstack((edges[0], edges[1]))
     j_indices = np.hstack((edges[1], edges[0]))
     data = np.hstack((-weights, -weights))
     lap = coo_matrix((data, (i_indices, j_indices)), shape=(pixel_nb, pixel_nb))
     w = -np.ravel(lap.sum(axis=1))
+    print w.min(), w.max()
     data *= 1. / (np.sqrt(w[i_indices]*w[j_indices]))
-    data = np.hstack((data, eps*np.ones_like(diag)))
-    i_indices = np.hstack((i_indices, diag))
-    j_indices = np.hstack((j_indices, diag))
+    #data = np.hstack((data, eps*np.ones_like(diag)))
+    #i_indices = np.hstack((i_indices, diag))
+    #j_indices = np.hstack((j_indices, diag))
     lap = coo_matrix((-data, (i_indices, j_indices)),\
             shape=(pixel_nb, pixel_nb))
     return lap.tocsc(), w
@@ -245,7 +273,10 @@ def _trim_edges_weights(edges, weights, mask):
 def _build_laplacian(data, mask=None, normed=False, beta=50):
     lx, ly, lz = data.shape
     edges = _make_edges_3d(lx, ly, lz)
-    weights = _make_weights_3d(edges, data, beta=beta, eps=1.e-10)
+    if beta==None:
+        weights = _make_adaptive_weights(edges, data)
+    else:
+        weights = _make_weights_3d(edges, data, beta=beta, eps=1.e-10)
     if mask is not None:
         edges, weights = _trim_edges_weights(edges, weights, mask)
     if not normed:

spectral_embedding.py

@@ -43,11 +43,13 @@ def spectral_embedding_sparse(adjacency, k_max=14, mode='amg', take_first=True):
         print 'amg'
         sh = adjacency.shape[0]
         adjacency = adjacency.copy()
-        diag = sparse.coo_matrix((diag_weights.ravel(), (range(sh), range(sh))))
+        #diag = sparse.coo_matrix((diag_weights.ravel(), (range(sh), range(sh))))
+        diag = sparse.eye(sh, sh)
         adjacency =  - adjacency + diag
         ml = smoothed_aggregation_solver(adjacency.tocsr())
         X = scipy.rand(adjacency.shape[0], k_max) 
-        X[:, 0] = 1. / np.sqrt(adjacency.shape[0])
+        #X[:, 0] = 1. / np.sqrt(adjacency.shape[0])
+        X[:, 0] = 1. / dd.ravel()
         M = ml.aspreconditioner()
         lambdas, diffusion_map = lobpcg(adjacency, X, M=M, tol=1.e-12, largest=False)
         print lambdas
@@ -145,14 +147,32 @@ def q_score(adjacency, labels):
         #q -= (weights[label == labels].sum()/total_weights)**2
     return 2 * q
 
+def n_cut(adjacency, labels):
+    """ Returns the Q score of a clustering.
+    """
+    q = 0
+    """
+    if isinstance(adjacency, sparse.csc.csc_matrix):
+        adjacency = np.array(adjacency.todense())
+    """
+    weights = adjacency
+    total_weights = 0.5 * weights.sum()
+    for label in np.unique(labels):
+        inds = np.nonzero(labels == label)[0]
+        a = (weights[inds][:, inds]).sum()
+        b = weights[inds].sum()
+        q += (b - a)/b
+    return - q
+
 
 def best_k_means(k, maps, adjacency, n_bst=10):
     from nipy.neurospin.clustering.clustering import _kmeans
     best_score = -np.inf 
     for _ in range(n_bst):
         print "doing kmeans"
         _, labels, _ = _kmeans(maps, nbclusters=k)
-        score = q_score(adjacency, labels)
+        #score = q_score(adjacency, labels)
+        score = n_cut(adjacency, labels)
         if score > best_score:
             best_score = score
             best_labels = labels
@@ -196,12 +216,12 @@ def communities_clustering_sparse(adjacency, k_best=None, k_min=2, k_max=8, n_bs
         this_maps = maps[:k_best - 1].T.copy()
         res, scores = best_k_means(k_best, this_maps, adjacency,
                                  n_bst=4*n_bst)
-        print 'Final : k=%i, score=%s' % (k_best, score)
+        print 'Final : k=%i, score=%s' % (k_best, scores)
     return res, scores
 
 def separate_in_regions(data, mask, k_best=None, k_min=2, k_max=8, \
-                                    center=None, only_connex=True, \
-                                    take_first=True, beta=10, mode='bf'):
+                                center=None, only_connex=True, n_times=4,\
+                                take_first=True, beta=10, mode='bf'):
     labs, nb_labels = ndimage.label(mask)
     if nb_labels > 1:
         if center is None:
@@ -212,10 +232,10 @@ def separate_in_regions(data, mask, k_best=None, k_min=2, k_max=8, \
         mask = labs == ind_max
     lap, w = _build_laplacian(np.atleast_3d(data), mask=np.atleast_3d(mask), \
                 normed=True, beta=beta)
-    return lap, w
     print lap.shape
     res, scores = communities_clustering_sparse(lap, k_best=k_best, \
-                    k_min=k_min, k_max=k_max, take_first=take_first, mode=mode)
+                    k_min=k_min, k_max=k_max, n_bst=n_times, \
+                    take_first=take_first, mode=mode)
     if not only_connex:
         if k_best==None:
             labels = dict()