Merge pull request #67 from smoia/fix/degree

Fix computation of degree matrix (from normal matrix to its adjacency)

Author: smoia

nigsp/operations/laplacian.py

@@ -70,8 +70,6 @@ def compute_laplacian(mtx, negval="absolute", selfloops=False):
                 f'Behaviour "{negval}" to deal with negative values is not supported'
             )
 
-    degree = mtx.sum(axis=1)  # This is fixed to across columns
-
     adjacency = deepcopy(mtx)
 
     if selfloops is False:
@@ -89,15 +87,16 @@ def compute_laplacian(mtx, negval="absolute", selfloops=False):
                 f"but specified matrix has {mtx.shape[0]} diagonal elements."
             )
         adjacency[np.diag_indices(adjacency.shape[0])] = selfloops
-        degree = degree + selfloops
     elif selfloops == "degree":
-        adjacency[np.diag_indices(adjacency.shape[0])] = degree
-        degree = degree * 2
+        adjacency[np.diag_indices(adjacency.shape[0])] = 0
+        adjacency[np.diag_indices(adjacency.shape[0])] = adjacency.sum(axis=1)
     else:
         raise NotImplementedError(
             f'Value "{selfloops}" for self-loops settings is not supported'
         )
 
+    degree = adjacency.sum(axis=1)  # This is fixed to across columns
+
     degree_mat = np.zeros_like(mtx)
     degree_mat[np.diag_indices(degree_mat.shape[0])] = degree
 
@@ -147,38 +146,38 @@ def normalisation(lapl, degree, norm="symmetric", fix_zeros=True):
     --------
     https://en.wikipedia.org/wiki/Laplacian_matrix
     """
-    degree = deepcopy(degree)
-    if lapl.ndim - degree.ndim > 1:
+    deg = deepcopy(degree)
+    if lapl.ndim - deg.ndim > 1:
         raise NotImplementedError(
-            f"The provided degree matrix is {degree.ndim}D while the "
+            f"The provided degree matrix is {deg.ndim}D while the "
             f"provided laplacian matrix is {lapl.ndim}D."
         )
-    elif lapl.ndim == degree.ndim:
-        if not (degree.diagonal() == degree.sum(axis=1)).all():
+    elif lapl.ndim == deg.ndim:
+        if not (deg.diagonal() == deg.sum(axis=1)).all():
             raise ValueError(
                 "The provided degree matrix is not a diagonal matrix (or a stack of)."
             )
-        degree = degree.diagonal()
+        deg = deepcopy(deg.diagonal())
 
-    if degree.shape != lapl.shape[:-1]:
+    if deg.shape != lapl.shape[:-1]:
         raise ValueError(
-            f"The provided degree matrix has shape {degree.shape} while the "
+            f"The provided degree matrix has shape {deg.shape} while the "
             f"provided matrix has shape {lapl.shape}."
         )
 
     if fix_zeros:
-        degree[degree == 0] = 1
+        deg[deg == 0] = 1
 
     d = np.zeros_like(lapl)
 
     if norm in ["symmetric", "symm"]:
-        d[np.diag_indices(d.shape[0])] = degree ** (-1 / 2)
+        d[np.diag_indices(d.shape[0])] = deg ** (-1 / 2)
         return d @ lapl @ d
     elif norm in ["random walk", "rw", "random walk inflow", "rwi"]:
-        d[np.diag_indices(d.shape[0])] = degree ** (-1)
+        d[np.diag_indices(d.shape[0])] = deg ** (-1)
         return d @ lapl
     elif norm in ["random walk outflow", "rwo"]:
-        d[np.diag_indices(d.shape[0])] = degree ** (-1)
+        d[np.diag_indices(d.shape[0])] = deg ** (-1)
         return lapl @ d
     else:
         raise NotImplementedError(f'Normalisation type "{norm}" is not supported.')

nigsp/tests/test_laplacian.py

@@ -12,41 +12,41 @@
 def test_compute_laplacian():
     def glap(mtx):
         deg = mtx.sum(axis=1)
-        adj = dc(mtx)
-        adj[np.diag_indices(mtx.shape[0])] = 0
 
-        L = np.diag(deg) - adj
+        L = np.diag(deg) - mtx
         return L, deg
 
     mtx = np.random.rand(4, 4)
     mtx = (mtx + mtx.T) / 2
 
-    L, deg = glap(mtx)
+    adj = dc(mtx)
+    adj[np.diag_indices(mtx.shape[0])] = 0
+    L, deg = glap(adj)
     lapl, degree = laplacian.compute_laplacian(mtx)
 
-    assert (lapl == L).all()
-    assert (degree == deg).all()
+    assert np.allclose(lapl, L)
+    assert np.allclose(degree, deg)
 
     lapl, degree = laplacian.compute_laplacian(mtx, selfloops=True)
-    Lsl = np.diag(mtx.sum(axis=1)) - mtx
-    assert (lapl == Lsl).all()
-    assert (degree == deg).all()
+    L, deg = glap(mtx)
+    assert np.allclose(lapl, L)
+    assert np.allclose(degree, deg)
 
     rn_deg = np.random.rand(4)
     lapl, degree = laplacian.compute_laplacian(mtx, selfloops=rn_deg)
-    updeg = deg + rn_deg
     adj = dc(mtx)
     adj[np.diag_indices(mtx.shape[0])] = rn_deg
-    Lsl = np.diag(updeg) - adj
-    assert (lapl == Lsl).all()
-    assert (degree == updeg).all()
+    L, deg = glap(adj)
+    assert np.allclose(lapl, L)
+    assert np.allclose(degree, deg)
 
     lapl, degree = laplacian.compute_laplacian(mtx, selfloops="degree")
-    updeg = deg + deg
+    adj[np.diag_indices(mtx.shape[0])] = 0
+    _, deg = glap(adj)
     adj[np.diag_indices(mtx.shape[0])] = deg
-    Lsl = np.diag(updeg) - adj
-    assert (lapl == Lsl).all()
-    assert (degree == updeg).all()
+    L, deg = glap(adj)
+    assert np.allclose(lapl, L)
+    assert np.allclose(degree, deg)
 
     mtx = mtx - mtx.mean()
 
@@ -56,22 +56,21 @@ def glap(mtx):
     mtx_res = (mtx - mtx.min()) / mtx.max()
 
     L, deg = glap(mtx_abs)
-    lapl, degree = laplacian.compute_laplacian(mtx, negval="absolute")
-    assert (lapl == L).all()
-    assert (degree == deg).all()
+    lapl, degree = laplacian.compute_laplacian(mtx, negval="absolute", selfloops=True)
+    assert np.allclose(lapl, L)
+    assert np.allclose(degree, deg)
 
     L, deg = glap(mtx_rem)
-    lapl, degree = laplacian.compute_laplacian(mtx, negval="remove")
-    assert (lapl == L).all()
-    assert (degree == deg).all()
+    lapl, degree = laplacian.compute_laplacian(mtx, negval="remove", selfloops=True)
+    assert np.allclose(lapl, L)
+    assert np.allclose(degree, deg)
 
     L, deg = glap(mtx_res)
-    lapl, degree = laplacian.compute_laplacian(mtx, negval="rescale")
-    assert (lapl == L).all()
-    assert (degree == deg).all()
+    lapl, degree = laplacian.compute_laplacian(mtx, negval="rescale", selfloops=True)
+    assert np.allclose(lapl, L)
+    assert np.allclose(degree, deg)
 
 
-@mark.xfail
 def test_normalisation():
     L = np.random.rand(4, 4)
     L = (L + L.T) / 2
@@ -80,10 +79,10 @@ def test_normalisation():
 
     lapl_symm = laplacian.normalisation(L, d, norm="symmetric")
 
-    d = np.diag(d)
-    lapl_rwi = laplacian.normalisation(L, d, norm="random walk")
+    deg = np.diag(d)
+    lapl_rwi = laplacian.normalisation(L, deg, norm="random walk")
 
-    lapl_rwo = laplacian.normalisation(L, d, norm="rwo")
+    lapl_rwo = laplacian.normalisation(L, deg, norm="rwo")
 
     d[2] = 1
     d_symm = np.diag(d ** (-1 / 2))