Add smoothness metric

nigsp/cli/run.py

@@ -91,6 +91,18 @@ def _get_parser():
         ),
         default=None,
     )
+    opt_metrics.add_argument(
+        "-smooth",
+        "--smoothness",
+        dest="comp_metric",
+        action="append_const",
+        const="smoothness",
+        help=(
+            "Compute the signal smoothness "
+            "(see Shuman et al, 2013, IEEE Signal Proc. Mag.)"
+        ),
+        default=None,
+    )
 
     opt_proc = parser.add_argument_group("Optional Arguments for data processing")
     opt_proc.add_argument(

nigsp/objects.py

@@ -43,6 +43,7 @@ def __init__(
         surr_split=None,
         sdi=None,
         gsdi=None,
+        smoothness=None,
         fc=None,
         fc_split=None,
     ):
@@ -82,6 +83,7 @@ def __init__(
         self.surr_split = deepcopy(surr_split)
         self.sdi = deepcopy(sdi)
         self.gsdi = deepcopy(gsdi)
+        self.smoothness = deepcopy(smoothness)
         self.fc = deepcopy(fc)
         self.fc_split = deepcopy(fc_split)
 
@@ -168,6 +170,11 @@ def compute_gsdi(self, mean=False, keys=None):  # pragma: no cover
         self.gsdi = operations.gsdi(self.ts_split, mean, keys)
         return self
 
+    def compute_smoothness(self, signal):  # pragma: no cover
+        """Implement metrics.smoothness as class method."""
+        self.smoothness = operations.smoothness(self.lapl_mtx, signal)
+        return self
+
     def create_surrogates(self, sc_type="informed", n_surr=1000, seed=None):
         """Implement surrogates.sc_informed and sc_uninformed as class method."""
         sc_args = {"timeseries": self.timeseries, "n_surr": n_surr}

nigsp/operations/__init__.py

@@ -15,7 +15,7 @@
 # Import all operations.
 from .graph import nodestrength, zerocross
 from .laplacian import decomposition, symmetric_normalised_laplacian
-from .metrics import functional_connectivity, gsdi, sdi
+from .metrics import functional_connectivity, gsdi, sdi, smoothness
 from .nifti import apply_atlas, apply_mask, mat_to_vol, unfold_atlas, unmask, vol_to_mat
 from .surrogates import random_sign, sc_informed, sc_uninformed, test_significance
 from .timeseries import (

nigsp/operations/metrics.py

@@ -19,7 +19,7 @@
 
 LGR = logging.getLogger(__name__)
 
-SUPPORTED_METRICS = ["sdi", "dfc", "fc"]
+SUPPORTED_METRICS = ["sdi", "dfc", "fc", "smoothness"]
 
 
 @due.dcite(references.PRETI_2019)
@@ -246,6 +246,57 @@ def _fc(timeseries, mean=False):
     return fc
 
 
+@due.dcite(references.SHUMAN_2013)
+def smoothness(laplacian, signal):
+    """Compute the smoothness of a signal (as defined in [1]) over the graph corresponding to given Laplacian.
+
+    Parameters
+    ----------
+    node_signal : numpy.ndarray
+        any signal defined with one value per node.
+    laplacian : numpy.ndarray
+        graph laplacian to use.
+
+    Returns
+    -------
+    smoothness : float
+        the smoothness of the signal.
+
+    References
+    ----------
+    .. [1] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega and P. Vandergheynst,
+       "The emerging field of signal processing on graphs: Extending high-dimensional
+        data analysis to networks and other irregular domains," in IEEE Signal
+        Processing Magazine, vol. 30, no. 3, pp. 83-98, May 2013
+    """
+    LGR.info("Compute signal smoothness.")
+
+    # Checking shape of signal
+    if signal.ndim == 1:
+        LGR.warning(
+            "2D array required, but signal is a vector. Adding empty dimension."
+        )
+        signal = np.expand_dims(signal, -1)
+    elif signal.ndim > 2:
+        raise ValueError("Signal should be a 2D array.")
+
+    # Checking if laplacian is square
+    if laplacian.shape[0] != laplacian.shape[1]:
+        raise ValueError("Laplacian should be a square matrix.")
+
+    # Checking that the dimensions of signal and laplacian are compatible
+    if signal.shape[0] != laplacian.shape[0]:
+        if signal.shape[1] == laplacian.shape[0]:
+            LGR.warning(
+                "It seems that the signal needs to be transposed to the correct shape."
+            )
+            signal = np.swapaxes(signal, 0, 1)
+        else:
+            raise ValueError("The dimensions of the signal and Laplacian don't match.")
+
+    return np.matmul(signal.T, np.matmul(laplacian, signal))
+
+
 """
 Copyright 2022, Stefano Moia.
 

nigsp/references.py

@@ -1,6 +1,8 @@
 """References to be imported and injected throughout the package."""
 from nigsp.due import Doi
 
+SHUMAN_2013 = Doi("10.1109/MSP.2012.2235192")
+
 PRETI_2019 = Doi("10.1038/s41467-019-12765-7")
 
 GRIFFA_2022 = Doi("10.1016/j.neuroimage.2022.118970")

nigsp/tests/test_metrics.py

@@ -3,7 +3,7 @@
 
 import numpy as np
 from numpy.random import rand
-from pytest import raises
+from pytest import raises, warns
 
 from nigsp.operations import metrics
 from nigsp.utils import prepare_ndim_iteration
@@ -72,6 +72,25 @@ def test_gsdi():
     assert (gsdi_out["beta_over_alpha"] == gsdi_in).all()
 
 
+def test_smoothness():
+    s1 = rand(10)
+    s2 = rand(10, 2)
+    s3 = rand(2, 10)
+    laplacian = rand(10, 10)
+
+    expected_smoothness1 = np.dot(s1.T, np.dot(laplacian, s1))
+    expected_smoothness2 = np.dot(s2.T, np.dot(laplacian, s2))
+    expected_smoothness3 = np.dot(s3, np.dot(laplacian, s3.T))
+
+    computed_smoothness1 = metrics.smoothness(laplacian, s1)
+    computed_smoothness2 = metrics.smoothness(laplacian, s2)
+    computed_smoothness3 = metrics.smoothness(laplacian, s3)
+
+    assert (expected_smoothness1 == computed_smoothness1).all()
+    assert (expected_smoothness2 == computed_smoothness2).all()
+    assert (expected_smoothness3 == computed_smoothness3).all()
+
+
 # ### Break tests
 def test_break_sdi():
     ts1 = np.arange(1, 3)[..., np.newaxis]
@@ -124,3 +143,29 @@ def test_functional_connectivity():
 
     for k in fcd.keys():
         assert (fcd[k] == np.corrcoef(tsd[k])).all()
+
+
+def test_break_smoothness():
+    # shape of signal
+    signal = rand(3, 3, 3)
+    laplacian = rand(3, 3)
+
+    with raises(ValueError) as errorinfo:
+        metrics.smoothness(laplacian, signal)
+    assert "should be a 2D" in str(errorinfo.value)
+
+    # shape of laplacian
+    signal = rand(10, 2)
+    laplacian = rand(10, 9)
+
+    with raises(ValueError) as errorinfo:
+        metrics.smoothness(laplacian, signal)
+    assert "a square matrix" in str(errorinfo.value)
+
+    # shape mismatch
+    signal = rand(10, 2)
+    laplacian = rand(9, 9)
+
+    with raises(ValueError) as errorinfo:
+        metrics.smoothness(laplacian, signal)
+    assert "don't match" in str(errorinfo.value)