[MARKER]: Support FunctionalConnectivityLagged

docs/builtin.rst

@@ -292,6 +292,14 @@ Available
        | `Richman et al. (2000) <https://doi.org/10.1152/ajpheart.2000.278.6.H2039>`_
      - Done
      - 0.0.4
+   * - :class:`.FunctionalConnectivityLaggedParcels`
+     - Compute lagged functional connectivity over parcellation
+     - Done
+     - 0.0.6
+   * - :class:`.FunctionalConnectivityLaggedSpheres`
+     - Compute lagged functional connectivity over spheres placed on coordinates
+     - Done
+     - 0.0.6
 
 Planned
 ~~~~~~~

docs/changes/newsfragments/435.feature

@@ -0,0 +1 @@
+Add support for :class:`.FunctionalConnectivityLaggedParcels` and :class:`.FunctionalConnectivityLaggedSpheres` markers by `Synchon Mandal`_

junifer/markers/__init__.pyi

@@ -5,6 +5,8 @@ __all__ = [
     "SphereAggregation",
     "FunctionalConnectivityParcels",
     "FunctionalConnectivitySpheres",
+    "FunctionalConnectivityLaggedParcels",
+    "FunctionalConnectivityLaggedSpheres",
     "CrossParcellationFC",
     "EdgeCentricFCParcels",
     "EdgeCentricFCSpheres",
@@ -24,6 +26,8 @@ from .sphere_aggregation import SphereAggregation
 from .functional_connectivity import (
     FunctionalConnectivityParcels,
     FunctionalConnectivitySpheres,
+    FunctionalConnectivityLaggedParcels,
+    FunctionalConnectivityLaggedSpheres,
     CrossParcellationFC,
     EdgeCentricFCParcels,
     EdgeCentricFCSpheres,

junifer/markers/functional_connectivity/__init__.pyi

@@ -1,13 +1,21 @@
 __all__ = [
     "FunctionalConnectivityParcels",
     "FunctionalConnectivitySpheres",
+    "FunctionalConnectivityLaggedParcels",
+    "FunctionalConnectivityLaggedSpheres",
     "CrossParcellationFC",
     "EdgeCentricFCParcels",
     "EdgeCentricFCSpheres",
 ]
 
 from .functional_connectivity_parcels import FunctionalConnectivityParcels
 from .functional_connectivity_spheres import FunctionalConnectivitySpheres
+from .functional_connectivity_lagged_parcels import (
+    FunctionalConnectivityLaggedParcels,
+)
+from .functional_connectivity_lagged_spheres import (
+    FunctionalConnectivityLaggedSpheres,
+)
 from .crossparcellation_functional_connectivity import CrossParcellationFC
 from .edge_functional_connectivity_parcels import EdgeCentricFCParcels
 from .edge_functional_connectivity_spheres import EdgeCentricFCSpheres

junifer/markers/functional_connectivity/functional_connectivity_lagged_base.py

@@ -0,0 +1,174 @@
+"""Provide abstract base class for lagged functional connectivity (FC)."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+
+from abc import abstractmethod
+from itertools import product
+from typing import Any, ClassVar, Optional, Union
+
+import numpy as np
+from scipy.signal import correlate
+
+from ...typing import Dependencies, MarkerInOutMappings
+from ...utils import raise_error
+from ..base import BaseMarker
+
+
+__all__ = ["FunctionalConnectivityLaggedBase"]
+
+
+class FunctionalConnectivityLaggedBase(BaseMarker):
+    """Abstract base class for lagged functional connectivity markers.
+
+    Parameters
+    ----------
+    max_lag : int
+        The time lag range. The lag ranges from ``-max_lag`` to ``+max_lag``
+        time points.
+    agg_method : str, optional
+        The method to perform aggregation using.
+        Check valid options in :func:`.get_aggfunc_by_name`
+        (default "mean").
+    agg_method_params : dict, optional
+        Parameters to pass to the aggregation function.
+        Check valid options in :func:`.get_aggfunc_by_name`
+        (default None).
+    masks : str, dict or list of dict or str, optional
+        The specification of the masks to apply to regions before extracting
+        signals. Check :ref:`Using Masks <using_masks>` for more details.
+        If None, will not apply any mask (default None).
+    name : str, optional
+        The name of the marker. If None, will use ``BOLD_<class_name>``
+        (default None).
+
+    """
+
+    _DEPENDENCIES: ClassVar[Dependencies] = {"numpy", "scipy"}
+
+    _MARKER_INOUT_MAPPINGS: ClassVar[MarkerInOutMappings] = {
+        "BOLD": {
+            "functional_connectivity": "matrix",
+            "lag": "matrix",
+        },
+    }
+
+    def __init__(
+        self,
+        max_lag: int,
+        agg_method: str = "mean",
+        agg_method_params: Optional[dict] = None,
+        masks: Union[str, dict, list[Union[dict, str]], None] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        self.max_lag = max_lag
+        self.agg_method = agg_method
+        self.agg_method_params = agg_method_params
+        self.masks = masks
+        super().__init__(on="BOLD", name=name)
+
+    @abstractmethod
+    def aggregate(
+        self,
+        input: dict[str, Any],
+        extra_input: Optional[dict[str, Any]] = None,
+    ) -> dict[str, Any]:
+        """Perform aggregation."""
+        raise_error(
+            msg="Concrete classes need to implement aggregate().",
+            klass=NotImplementedError,
+        )
+
+    def compute(
+        self,
+        input: dict[str, Any],
+        extra_input: Optional[dict] = None,
+    ) -> dict:
+        """Compute.
+
+        Parameters
+        ----------
+        input : dict
+            A single input from the pipeline data object in which to compute
+            the marker.
+        extra_input : dict, optional
+            The other fields in the pipeline data object. Useful for accessing
+            other data kind that needs to be used in the computation. For
+            example, the functional connectivity markers can make use of the
+            confounds if available (default None).
+
+        Returns
+        -------
+        dict
+            The computed result as dictionary. This will be either returned
+            to the user or stored in the storage by calling the store method
+            with this as a parameter. The dictionary has the following keys:
+
+            * ``functional_connectivity`` : dictionary with the following keys:
+
+              - ``data`` : functional connectivity matrix as ``numpy.ndarray``
+              - ``row_names`` : ROI labels as list of str
+              - ``col_names`` : ROI labels as list of str
+              - ``matrix_kind`` : the kind of matrix (tril, triu or full)
+
+            * ``lag`` : dictionary with the following keys:
+
+              - ``data`` : lag matrix as ``numpy.ndarray``
+              - ``row_names`` : ROI labels as list of str
+              - ``col_names`` : ROI labels as list of str
+              - ``matrix_kind`` : the kind of matrix (tril, triu or full)
+
+        Notes
+        -----
+        Pearson correlation is used to perform connectivity measure.
+
+        """
+        # Perform necessary aggregation
+        aggregation = self.aggregate(input, extra_input=extra_input)
+        # Initialize variables
+        # transposed to get (n_rois, n_timepoints)
+        fmri_data = aggregation["aggregation"]["data"].T
+        n_rois = fmri_data.shape[0]
+        fc_matrix = np.ones((n_rois, n_rois))
+        lag_matrix = np.zeros((n_rois, n_rois), dtype=int)
+        lags = np.arange(-self.max_lag, self.max_lag + 1)
+        # Compute
+        for i, j in product(range(n_rois), range(n_rois)):
+            if i != j:
+                x = fmri_data[i]
+                y = fmri_data[j]
+                # Compute cross-correlation function (CCF)
+                ccf = correlate(x, y, mode="full", method="auto")
+                # Normalize CCF
+                ccf /= np.sqrt(np.sum(np.abs(x**2)) * np.sum(np.abs(y**2)))
+                # Limit lag range
+                ccf = ccf[
+                    len(ccf) // 2
+                    - self.max_lag : len(ccf) // 2
+                    + self.max_lag
+                    + 1
+                ]
+                # Find peak correlation and corresponding lag
+                peak_idx = np.argmax(np.abs(ccf))  # Peak correlation index
+                peak_corr = ccf[peak_idx]  # Peak correlation value
+                peak_lag = lags[peak_idx]  # Corresponding lag
+                # Store in matrices
+                fc_matrix[i, j] = peak_corr
+                lag_matrix[i, j] = peak_lag
+        # Create dictionary for output
+        roi_labels = aggregation["aggregation"]["col_names"]
+        return {
+            "functional_connectivity": {
+                "data": fc_matrix,
+                "row_names": roi_labels,
+                "col_names": roi_labels,
+                "matrix_kind": "tril",
+            },
+            "lag": {
+                "data": lag_matrix,
+                "row_names": roi_labels,
+                "col_names": roi_labels,
+                "matrix_kind": "full",
+            },
+        }

junifer/markers/functional_connectivity/functional_connectivity_lagged_parcels.py

@@ -0,0 +1,102 @@
+"""Provide class for lagged functional connectivity using parcels."""
+
+# Authors: Synchon Mandal <s.mandal@fz-juelich.de>
+# License: AGPL
+
+
+from typing import Any, Optional, Union
+
+from ...api.decorators import register_marker
+from ..parcel_aggregation import ParcelAggregation
+from .functional_connectivity_lagged_base import (
+    FunctionalConnectivityLaggedBase,
+)
+
+
+__all__ = ["FunctionalConnectivityLaggedParcels"]
+
+
+@register_marker
+class FunctionalConnectivityLaggedParcels(FunctionalConnectivityLaggedBase):
+    """Class for lagged functional connectivity using parcellations.
+
+    Parameters
+    ----------
+    parcellation : str or list of str
+        The name(s) of the parcellation(s) to use.
+        See :func:`.list_data` for options.
+    max_lag : int
+        The time lag range. The lag ranges from ``-max_lag`` to ``+max_lag``
+        time points.
+    agg_method : str, optional
+        The method to perform aggregation using.
+        See :func:`.get_aggfunc_by_name` for options
+        (default "mean").
+    agg_method_params : dict, optional
+        Parameters to pass to the aggregation function.
+        See :func:`.get_aggfunc_by_name` for options
+        (default None).
+    masks : str, dict or list of dict or str, optional
+        The specification of the masks to apply to regions before extracting
+        signals. Check :ref:`Using Masks <using_masks>` for more details.
+        If None, will not apply any mask (default None).
+    name : str, optional
+        The name of the marker. If None, will use
+        ``BOLD_FunctionalConnectivityLaggedParcels`` (default None).
+
+    """
+
+    def __init__(
+        self,
+        parcellation: Union[str, list[str]],
+        max_lag: int,
+        agg_method: str = "mean",
+        agg_method_params: Optional[dict] = None,
+        masks: Union[str, dict, list[Union[dict, str]], None] = None,
+        name: Optional[str] = None,
+    ) -> None:
+        self.parcellation = parcellation
+        super().__init__(
+            max_lag=max_lag,
+            agg_method=agg_method,
+            agg_method_params=agg_method_params,
+            masks=masks,
+            name=name,
+        )
+
+    def aggregate(
+        self, input: dict[str, Any], extra_input: Optional[dict] = None
+    ) -> dict:
+        """Perform parcel aggregation.
+
+        Parameters
+        ----------
+        input : dict
+            A single input from the pipeline data object in which to compute
+            the marker.
+        extra_input : dict, optional
+            The other fields in the pipeline data object. Useful for accessing
+            other data kind that needs to be used in the computation. For
+            example, the functional connectivity markers can make use of the
+            confounds if available (default None).
+
+        Returns
+        -------
+        dict
+            The computed result as dictionary. This will be either returned
+            to the user or stored in the storage by calling the store method
+            with this as a parameter. The dictionary has the following keys:
+
+            * ``aggregation`` : dictionary with the following keys:
+
+                - ``data`` : ROI values as ``numpy.ndarray``
+                - ``col_names`` : ROI labels as list of str
+
+        """
+        return ParcelAggregation(
+            parcellation=self.parcellation,
+            method=self.agg_method,
+            method_params=self.agg_method_params,
+            masks=self.masks,
+            on="BOLD",
+        ).compute(input=input, extra_input=extra_input)