Merge branch 'pull-266-laplacian' into master.

Reviewed on PR-266. Closes ticket 1681.

scipy/sparse/csgraph/_laplacian.py

@@ -8,7 +8,8 @@
 # License: BSD
 
 import numpy as np
-from scipy.sparse import isspmatrix
+from scipy.sparse import isspmatrix, coo_matrix
+
 
 ###############################################################################
 # Graph laplacian
@@ -31,7 +32,7 @@ def laplacian(csgraph, normed=False, return_diag=False):
     -------
     lap: ndarray
         The N x N laplacian matrix of graph.
-    
+
     diag: ndarray
         The length-N diagonal of the laplacian matrix.
         diag is returned only if return_diag is True.
@@ -70,7 +71,7 @@ def laplacian(csgraph, normed=False, return_diag=False):
                    or np.issubdtype(csgraph.dtype, np.uint)):
         csgraph = csgraph.astype(np.float)
 
-    if isspmatrix(csgraph):            
+    if isspmatrix(csgraph):
         return _laplacian_sparse(csgraph, normed=normed,
                                  return_diag=return_diag)
     else:
@@ -89,14 +90,14 @@ def _laplacian_sparse(graph, normed=False, return_diag=False):
         # The sparsity pattern of the matrix has holes on the diagonal,
         # we need to fix that
         diag_idx = lap.row[diag_mask]
-
-        lap = lap.tolil()
-
         diagonal_holes = list(set(range(n_nodes)).difference(
                                 diag_idx))
-        lap[diagonal_holes, diagonal_holes] = 1
-        lap = lap.tocoo()
+        new_data = np.concatenate([lap.data, np.ones(len(diagonal_holes))])
+        new_row = np.concatenate([lap.row, diagonal_holes])
+        new_col = np.concatenate([lap.col, diagonal_holes])
+        lap = coo_matrix((new_data, (new_row, new_col)), shape=lap.shape)
         diag_mask = (lap.row == lap.col)
+
     lap.data[diag_mask] = 0
     w = -np.asarray(lap.sum(axis=1)).squeeze()
     if normed:
@@ -105,17 +106,18 @@ def _laplacian_sparse(graph, normed=False, return_diag=False):
         w[w_zeros] = 1
         lap.data /= w[lap.row]
         lap.data /= w[lap.col]
-        lap.data[diag_mask] = (1 - w_zeros).astype(lap.data.dtype)
+        lap.data[diag_mask] = (1 - w_zeros[lap.row[diag_mask]]).astype(lap.data.dtype)
     else:
         lap.data[diag_mask] = w[lap.row[diag_mask]]
+
     if return_diag:
         return lap, w
     return lap
 
 
 def _laplacian_dense(graph, normed=False, return_diag=False):
     n_nodes = graph.shape[0]
-    lap = -graph.copy()
+    lap = -np.asarray(graph)  # minus sign leads to a copy
 
     # set diagonal to zero
     lap.flat[::n_nodes + 1] = 0
@@ -129,7 +131,7 @@ def _laplacian_dense(graph, normed=False, return_diag=False):
         lap.flat[::n_nodes + 1] = 1 - w_zeros
     else:
         lap.flat[::n_nodes + 1] = w
+
     if return_diag:
         return lap, w
     return lap
-

scipy/sparse/csgraph/tests/test_graph_laplacian.py

@@ -1,27 +1,66 @@
 # Author: Gael Varoquaux <gael.varoquaux@normalesup.org>
+#         Jake Vanderplas <vanderplas@astro.washington.edu>
 # License: BSD
 
 import numpy as np
 from scipy import sparse
 
 from scipy.sparse import csgraph
 
+def _explicit_laplacian(x, normed=False):
+    if sparse.issparse(x):
+        x = x.todense()
+    x = np.asarray(x)
+    y = -1.0 * x
+    for j in range(y.shape[0]):
+        y[j,j] = x[j,j+1:].sum() + x[j,:j].sum()
+    if normed:
+        d = np.diag(y).copy()
+        d[d == 0] = 1.0
+        y /= d[:,None]**.5
+        y /= d[None,:]**.5
+    return y
 
-def test_graph_laplacian():
-    for mat in (np.arange(10) * np.arange(10)[:, np.newaxis],
-                np.ones((7, 7)),
-                np.eye(19),
-                np.vander(np.arange(4)) + np.vander(np.arange(4)).T,
-               ):
+def _check_graph_laplacian(mat, normed):
+    if not hasattr(mat, 'shape'):
+        mat = eval(mat, dict(np=np, sparse=sparse))
+
+    if sparse.issparse(mat):
+        sp_mat = mat
+        mat = sp_mat.todense()
+    else:
         sp_mat = sparse.csr_matrix(mat)
+
+    laplacian = csgraph.laplacian(mat, normed=normed)
+    n_nodes = mat.shape[0]
+    if not normed:
+        np.testing.assert_array_almost_equal(laplacian.sum(axis=0),
+                                             np.zeros(n_nodes))
+    np.testing.assert_array_almost_equal(laplacian.T,
+                                         laplacian)
+    np.testing.assert_array_almost_equal(\
+        laplacian,
+        csgraph.laplacian(sp_mat, normed=normed).todense())
+
+    np.testing.assert_array_almost_equal(
+        laplacian,
+        _explicit_laplacian(mat, normed=normed))
+
+def test_graph_laplacian():
+    mats = ('np.arange(10) * np.arange(10)[:, np.newaxis]',
+            'np.ones((7, 7))',
+            'np.eye(19)',
+            'sparse.diags([1, 1], [-1, 1], shape=(4,4))',
+            'sparse.diags([1, 1], [-1, 1], shape=(4,4)).todense()',
+            'np.asarray(sparse.diags([1, 1], [-1, 1], shape=(4,4)).todense())',
+            'np.vander(np.arange(4)) + np.vander(np.arange(4)).T',
+            )
+
+    for mat_str in mats:
         for normed in (True, False):
-            laplacian = csgraph.laplacian(mat, normed=normed)
-            n_nodes = mat.shape[0]
-            if not normed:
-                np.testing.assert_array_almost_equal(laplacian.sum(axis=0),
-                                                     np.zeros(n_nodes))
-            np.testing.assert_array_almost_equal(laplacian.T,
-                                                 laplacian)
-            np.testing.assert_array_almost_equal(\
-                laplacian,
-                csgraph.laplacian(sp_mat, normed=normed).todense())
+            yield _check_graph_laplacian, mat_str, normed
+
+
+if __name__ == '__main__':
+    import nose
+    nose.runmodule()