1# Author: Gael Varoquaux <gael.varoquaux@normalesup.org>
2#         Jake Vanderplas <vanderplas@astro.washington.edu>
3# License: BSD
4import numpy as np
5from numpy.testing import assert_allclose, assert_array_almost_equal
6from pytest import raises as assert_raises
7from scipy import sparse
8
9from scipy.sparse import csgraph
10
11
12def _explicit_laplacian(x, normed=False):
13    if sparse.issparse(x):
14        x = x.todense()
15    x = np.asarray(x)
16    y = -1.0 * x
17    for j in range(y.shape[0]):
18        y[j,j] = x[j,j+1:].sum() + x[j,:j].sum()
19    if normed:
20        d = np.diag(y).copy()
21        d[d == 0] = 1.0
22        y /= d[:,None]**.5
23        y /= d[None,:]**.5
24    return y
25
26
27def _check_symmetric_graph_laplacian(mat, normed):
28    if not hasattr(mat, 'shape'):
29        mat = eval(mat, dict(np=np, sparse=sparse))
30
31    if sparse.issparse(mat):
32        sp_mat = mat
33        mat = sp_mat.todense()
34    else:
35        sp_mat = sparse.csr_matrix(mat)
36
37    laplacian = csgraph.laplacian(mat, normed=normed)
38    n_nodes = mat.shape[0]
39    if not normed:
40        assert_array_almost_equal(laplacian.sum(axis=0), np.zeros(n_nodes))
41    assert_array_almost_equal(laplacian.T, laplacian)
42    assert_array_almost_equal(laplacian,
43            csgraph.laplacian(sp_mat, normed=normed).todense())
44
45    assert_array_almost_equal(laplacian,
46            _explicit_laplacian(mat, normed=normed))
47
48
49def test_laplacian_value_error():
50    for t in int, float, complex:
51        for m in ([1, 1],
52                  [[[1]]],
53                  [[1, 2, 3], [4, 5, 6]],
54                  [[1, 2], [3, 4], [5, 5]]):
55            A = np.array(m, dtype=t)
56            assert_raises(ValueError, csgraph.laplacian, A)
57
58
59def test_symmetric_graph_laplacian():
60    symmetric_mats = ('np.arange(10) * np.arange(10)[:, np.newaxis]',
61            'np.ones((7, 7))',
62            'np.eye(19)',
63            'sparse.diags([1, 1], [-1, 1], shape=(4,4))',
64            'sparse.diags([1, 1], [-1, 1], shape=(4,4)).todense()',
65            'np.asarray(sparse.diags([1, 1], [-1, 1], shape=(4,4)).todense())',
66            'np.vander(np.arange(4)) + np.vander(np.arange(4)).T')
67    for mat_str in symmetric_mats:
68        for normed in True, False:
69            _check_symmetric_graph_laplacian(mat_str, normed)
70
71
72def _assert_allclose_sparse(a, b, **kwargs):
73    # helper function that can deal with sparse matrices
74    if sparse.issparse(a):
75        a = a.toarray()
76    if sparse.issparse(b):
77        b = a.toarray()
78    assert_allclose(a, b, **kwargs)
79
80
81def _check_laplacian(A, desired_L, desired_d, normed, use_out_degree):
82    for arr_type in np.array, sparse.csr_matrix, sparse.coo_matrix:
83        for t in int, float, complex:
84            adj = arr_type(A, dtype=t)
85            L = csgraph.laplacian(adj, normed=normed, return_diag=False,
86                                  use_out_degree=use_out_degree)
87            _assert_allclose_sparse(L, desired_L, atol=1e-12)
88            L, d = csgraph.laplacian(adj, normed=normed, return_diag=True,
89                                  use_out_degree=use_out_degree)
90            _assert_allclose_sparse(L, desired_L, atol=1e-12)
91            _assert_allclose_sparse(d, desired_d, atol=1e-12)
92
93
94def test_asymmetric_laplacian():
95    # adjacency matrix
96    A = [[0, 1, 0],
97         [4, 2, 0],
98         [0, 0, 0]]
99
100    # Laplacian matrix using out-degree
101    L = [[1, -1, 0],
102         [-4, 4, 0],
103         [0, 0, 0]]
104    d = [1, 4, 0]
105    _check_laplacian(A, L, d, normed=False, use_out_degree=True)
106
107    # normalized Laplacian matrix using out-degree
108    L = [[1, -0.5, 0],
109         [-2, 1, 0],
110         [0, 0, 0]]
111    d = [1, 2, 1]
112    _check_laplacian(A, L, d, normed=True, use_out_degree=True)
113
114    # Laplacian matrix using in-degree
115    L = [[4, -1, 0],
116         [-4, 1, 0],
117         [0, 0, 0]]
118    d = [4, 1, 0]
119    _check_laplacian(A, L, d, normed=False, use_out_degree=False)
120
121    # normalized Laplacian matrix using in-degree
122    L = [[1, -0.5, 0],
123         [-2, 1, 0],
124         [0, 0, 0]]
125    d = [2, 1, 1]
126    _check_laplacian(A, L, d, normed=True, use_out_degree=False)
127
128
129def test_sparse_formats():
130    for fmt in ('csr', 'csc', 'coo', 'lil', 'dok', 'dia', 'bsr'):
131        mat = sparse.diags([1, 1], [-1, 1], shape=(4,4), format=fmt)
132        for normed in True, False:
133            _check_symmetric_graph_laplacian(mat, normed)
134
135