# scipy/scipy

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

`Reviewed on PR-266. Closes ticket 1681.`
rgommers committed Jul 17, 2012
2 parents 65e8ad4 + 9095c6d commit 968939976b9c808c048077e813d27eb7f74c7f4c
Showing with 68 additions and 27 deletions.
1. +13 −11 scipy/sparse/csgraph/_laplacian.py
2. +55 −16 scipy/sparse/csgraph/tests/test_graph_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 -
 @@ -1,27 +1,66 @@ # Author: Gael Varoquaux +# Jake Vanderplas # 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()