diff --git a/docs/changes/newsfragments/348.change b/docs/changes/newsfragments/348.change new file mode 100644 index 000000000..60b77a9e1 --- /dev/null +++ b/docs/changes/newsfragments/348.change @@ -0,0 +1 @@ +For :class:`.CrossParcellationFC`, ``aggregation_method`` and ``correlation_method`` have been renamed to ``agg_method`` and ``corr_method`` respectively and ``agg_method_params`` has been added; for ``FunctionalConnectivityBase``, :class:`.FunctionalConnectivityParcels`, :class:`.FunctionalConnectivitySpheres`, :class:`.EdgeCentricFCParcels` and :class:`.EdgeCentricFCSpheres`, ``cor_method`` and ``cor_method_params`` have been renamed to ``conn_method`` and ``conn_method_params`` by `Synchon Mandal`_ diff --git a/docs/changes/newsfragments/348.enh b/docs/changes/newsfragments/348.enh new file mode 100644 index 000000000..592390cf9 --- /dev/null +++ b/docs/changes/newsfragments/348.enh @@ -0,0 +1 @@ +``FunctionalConnectivity``-family Markers now use :class:`sklearn.covariance.EmpiricalCovariance` as the default covariance estimator and ``correlation`` as the default connecivity matrix kind by `Synchon Mandal`_ diff --git a/docs/changes/newsfragments/348.feature b/docs/changes/newsfragments/348.feature new file mode 100644 index 000000000..3705f1a86 --- /dev/null +++ b/docs/changes/newsfragments/348.feature @@ -0,0 +1 @@ +Introduce :class:`.JuniferConnectivityMeasure` for customising functional connectivity matrix kinds and measurements by `Synchon Mandal`_ diff --git a/docs/conf.py b/docs/conf.py index bf3ee98e4..bd22517c9 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -170,6 +170,7 @@ numpydoc_xref_ignore = { "Engine", "positive", "negative", + "estimator", } # numpydoc_validation_checks = { # "all", diff --git a/ignore_words.txt b/ignore_words.txt index d0940ab0b..b708b7efe 100644 --- a/ignore_words.txt +++ b/ignore_words.txt @@ -5,3 +5,4 @@ sepulcre arange sinc whit +amin diff --git a/junifer/external/nilearn/__init__.py b/junifer/external/nilearn/__init__.py index 4fdb3c40c..f5806212c 100644 --- a/junifer/external/nilearn/__init__.py +++ b/junifer/external/nilearn/__init__.py @@ -4,6 +4,7 @@ # License: AGPL from .junifer_nifti_spheres_masker import JuniferNiftiSpheresMasker +from .junifer_connectivity_measure import JuniferConnectivityMeasure -__all__ = ["JuniferNiftiSpheresMasker"] +__all__ = ["JuniferNiftiSpheresMasker", "JuniferConnectivityMeasure"] diff --git a/junifer/external/nilearn/junifer_connectivity_measure.py b/junifer/external/nilearn/junifer_connectivity_measure.py new file mode 100644 index 000000000..70b4b6d27 --- /dev/null +++ b/junifer/external/nilearn/junifer_connectivity_measure.py @@ -0,0 +1,483 @@ +"""Provide JuniferConnectivityMeasure class.""" + +# Authors: Synchon Mandal +# License: AGPL + +from typing import Callable, List, Optional + +import numpy as np +from nilearn import signal +from nilearn.connectome import ( + ConnectivityMeasure, + cov_to_corr, + prec_to_partial, + sym_matrix_to_vec, +) +from scipy import linalg +from sklearn.base import clone +from sklearn.covariance import EmpiricalCovariance + +from ...utils import logger, raise_error, warn_with_log + + +__all__ = ["JuniferConnectivityMeasure"] + + +DEFAULT_COV_ESTIMATOR = EmpiricalCovariance(store_precision=False) + + +# New BSD License + +# Copyright (c) The nilearn developers. +# All rights reserved. + + +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: + +# a. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# b. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# c. Neither the name of the nilearn developers nor the names of +# its contributors may be used to endorse or promote products +# derived from this software without specific prior written +# permission. + + +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR +# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +# DAMAGE. + + +def _check_square(matrix: np.ndarray) -> None: + """Raise a ValueError if the input matrix is square. + + Parameters + ---------- + matrix : numpy.ndarray + Input array. + + Raises + ------ + ValueError + If ``matrix`` is not a square matrix. + + """ + if matrix.ndim != 2 or (matrix.shape[0] != matrix.shape[-1]): + raise_error( + f"Expected a square matrix, got array of shape {matrix.shape}." + ) + + +def is_spd(M: np.ndarray, decimal: int = 15) -> bool: # noqa: N803 + """Check that input matrix is symmetric positive definite. + + ``M`` must be symmetric down to specified ``decimal`` places. + The check is performed by checking that all eigenvalues are positive. + + Parameters + ---------- + M : numpy.ndarray + Input matrix to check for symmetric positive definite. + decimal : int, optional + Decimal places to check (default 15). + + Returns + ------- + bool + True if matrix is symmetric positive definite, False otherwise. + + """ + if not np.allclose(M, M.T, atol=0, rtol=10**-decimal): + logger.debug(f"matrix not symmetric to {decimal:d} decimals") + return False + eigvalsh = np.linalg.eigvalsh(M) + ispd = eigvalsh.min() > 0 + if not ispd: + logger.debug(f"matrix has a negative eigenvalue: {eigvalsh.min():.3f}") + return ispd + + +def _check_spd(matrix: np.ndarray) -> None: + """Check ``matrix`` is symmetric positive definite. + + Parameters + ---------- + matrix : numpy.ndarray + Input array. + + Raises + ------ + ValueError + If the input matrix is not symmetric positive definite. + + """ + if not is_spd(matrix, decimal=7): + raise_error("Expected a symmetric positive definite matrix.") + + +def _form_symmetric( + function: Callable[[np.ndarray], np.ndarray], + eigenvalues: np.ndarray, + eigenvectors: np.ndarray, +) -> np.ndarray: + """Return the symmetric matrix. + + Apply ``function`` to ``eigenvalues``, construct symmetric matrix with it + and ``eigenvectors`` and return the constructed symmetric matrix. + + Parameters + ---------- + function : callable (function numpy.ndarray -> numpy.ndarray) + The transform to apply to the eigenvalues. + eigenvalues : numpy.ndarray of shape (n_features, ) + Input argument of the function. + eigenvectors : numpy.ndarray of shape (n_features, n_features) + Unitary matrix. + + Returns + ------- + numpy.ndarray of shape (n_features, n_features) + The symmetric matrix obtained after transforming the eigenvalues, while + keeping the same eigenvectors. + + """ + return np.dot(eigenvectors * function(eigenvalues), eigenvectors.T) + + +def _map_eigenvalues( + function: Callable[[np.ndarray], np.ndarray], symmetric: np.ndarray +) -> np.ndarray: + """Matrix function, for real symmetric matrices. + + The function is applied to the eigenvalues of ``symmetric``. + + Parameters + ---------- + function : callable (function numpy.ndarray -> numpy.ndarray) + The transform to apply to the eigenvalues. + symmetric : numpy.ndarray of shape (n_features, n_features) + The input symmetric matrix. + + Returns + ------- + numpy.ndarray of shape (n_features, n_features) + The new symmetric matrix obtained after transforming the eigenvalues, + while keeping the same eigenvectors. + + Notes + ----- + If input matrix is not real symmetric, no error is reported but result will + be wrong. + + """ + eigenvalues, eigenvectors = linalg.eigh(symmetric) + return _form_symmetric(function, eigenvalues, eigenvectors) + + +def _geometric_mean( + matrices: List[np.ndarray], + init: Optional[np.ndarray] = None, + max_iter: int = 10, + tol: Optional[float] = 1e-7, +) -> np.ndarray: + """Compute the geometric mean of symmetric positive definite matrices. + + The geometric mean of ``n`` positive definite matrices + ``M_1, ..., M_n`` is the minimizer of the sum of squared distances from an + arbitrary matrix to each input matrix ``M_k`` + + .. math:: gmean(M_1, ..., M_n) = argmin_X sum_{k=1}^N dist(X, M_k)^2 + + where the used distance is related to matrices logarithm + + .. math:: dist(X, M_k) = ||log(X^{-1/2} M_k X^{-1/2)}|| + + In case of positive numbers, this mean is the usual geometric mean. + + See Algorithm 3 of [1]_ . + + Parameters + ---------- + matrices : list of numpy.ndarray, all of shape (n_features, n_features) + List of matrices whose geometric mean to compute. Raise an error if the + matrices are not all symmetric positive definite of the same shape. + init : numpy.ndarray of shape (n_features, n_features), optional + Initialization matrix, default to the arithmetic mean of matrices. + Raise an error if the matrix is not symmetric positive definite of the + same shape as the elements of matrices (default None). + max_iter : int, optional + Maximal number of iterations (default 10). + tol : positive float or None, optional + The tolerance to declare convergence: if the gradient norm goes below + this value, the gradient descent is stopped. If None, no check is + performed (default 1e-7). + + Returns + ------- + gmean : numpy.ndarray of shape (n_features, n_features) + Geometric mean of the matrices. + + References + ---------- + .. [1] Fletcher, T., P., Joshi, S. + Riemannian geometry for the statistical analysis of diffusion tensor + data. + Signal Processing, Volume 87, Issue 2, 2007, Pages 250-262 + https://doi.org/10.1016/j.sigpro.2005.12.018. + + """ + # Shape and symmetry positive definiteness checks + n_features = matrices[0].shape[0] + for matrix in matrices: + _check_square(matrix) + if matrix.shape[0] != n_features: + raise_error("Matrices are not of the same shape.") + _check_spd(matrix) + + # Initialization + matrices = np.array(matrices) + if init is None: + gmean = np.mean(matrices, axis=0) + else: + _check_square(init) + if init.shape[0] != n_features: + raise_error("Initialization has incorrect shape.") + _check_spd(init) + gmean = init + + norm_old = np.inf + step = 1.0 + + # Gradient descent + for _ in range(max_iter): + # Computation of the gradient + vals_gmean, vecs_gmean = linalg.eigh(gmean) + gmean_inv_sqrt = _form_symmetric(np.sqrt, 1.0 / vals_gmean, vecs_gmean) + whitened_matrices = [ + gmean_inv_sqrt.dot(matrix).dot(gmean_inv_sqrt) + for matrix in matrices + ] + logs = [_map_eigenvalues(np.log, w_mat) for w_mat in whitened_matrices] + # Covariant derivative is - gmean.dot(logms_mean) + logs_mean = np.mean(logs, axis=0) + if np.any(np.isnan(logs_mean)): + raise_error( + klass=FloatingPointError, + msg="Nan value after logarithm operation.", + ) + + # Norm of the covariant derivative on the tangent space at point gmean + norm = np.linalg.norm(logs_mean) + + # Update of the minimizer + vals_log, vecs_log = linalg.eigh(logs_mean) + gmean_sqrt = _form_symmetric(np.sqrt, vals_gmean, vecs_gmean) + # Move along the geodesic + gmean = gmean_sqrt.dot( + _form_symmetric(np.exp, vals_log * step, vecs_log) + ).dot(gmean_sqrt) + + # Update the norm and the step size + if norm < norm_old: + norm_old = norm + elif norm > norm_old: + step = step / 2.0 + norm = norm_old + if tol is not None and norm / gmean.size < tol: + break + if tol is not None and norm / gmean.size >= tol: + warn_with_log( + f"Maximum number of iterations {max_iter} reached without " + f"getting to the requested tolerance level {tol}." + ) + + return gmean + + +class JuniferConnectivityMeasure(ConnectivityMeasure): + """Class for custom ConnectivityMeasure. + + Differs from :class:`nilearn.connectome.ConnectivityMeasure` in the + following ways: + + * default ``cov_estimator`` is + :class:`sklearn.covariance.EmpiricalCovariance` + * default ``kind`` is ``"correlation"`` + + Parameters + ---------- + cov_estimator : estimator object, optional + The covariance estimator + (default ``EmpiricalCovariance(store_precision=False)``). + kind : {"covariance", "correlation", "partial correlation", \ + "tangent", "precision"}, optional + The matrix kind. For the use of ``"tangent"`` see [1]_ + (default "correlation"). + vectorize : bool, optional + If True, connectivity matrices are reshaped into 1D arrays and only + their flattened lower triangular parts are returned (default False). + discard_diagonal : bool, optional + If True, vectorized connectivity coefficients do not include the + matrices diagonal elements. Used only when vectorize is set to True + (default False). + standardize : bool, optional + If standardize is True, the data are centered and normed: their mean + is put to 0 and their variance is put to 1 in the time dimension + (default True). + + .. note:: + + Added to control passing value to ``standardize`` of + ``signal.clean`` to call new behavior since passing ``"zscore"`` or + True (default) is deprecated. This parameter will be deprecated in + version 0.13 and removed in version 0.15. + + Attributes + ---------- + cov_estimator_ : estimator object + A new covariance estimator with the same parameters as + ``cov_estimator``. + mean_ : numpy.ndarray + The mean connectivity matrix across subjects. For ``"tangent"`` kind, + it is the geometric mean of covariances (a group covariance + matrix that captures information from both correlation and partial + correlation matrices). For other values for ``kind``, it is the + mean of the corresponding matrices. + whitening_ : numpy.ndarray + The inverted square-rooted geometric mean of the covariance matrices. + + References + ---------- + .. [1] Varoquaux, G., Baronnet, F., Kleinschmidt, A. et al. + Detection of brain functional-connectivity difference in + post-stroke patients using group-level covariance modeling. + In Tianzi Jiang, Nassir Navab, Josien P. W. Pluim, and + Max A. Viergever, editors, Medical image computing and + computer-assisted intervention - MICCAI 2010, Lecture notes + in computer science, Pages 200-208. Berlin, Heidelberg, 2010. + Springer. + doi:10/cn2h9c. + + """ + + def __init__( + self, + cov_estimator=DEFAULT_COV_ESTIMATOR, + kind="correlation", + vectorize=False, + discard_diagonal=False, + standardize=True, + ): + super().__init__( + cov_estimator=cov_estimator, + kind=kind, + vectorize=vectorize, + discard_diagonal=discard_diagonal, + standardize=standardize, + ) + + def _fit_transform( + self, + X, # noqa: N803 + do_transform=False, + do_fit=False, + confounds=None, + ): + """Avoid duplication of computation.""" + self._check_input(X, confounds=confounds) + if do_fit: + self.cov_estimator_ = clone(self.cov_estimator) + + # Compute all the matrices, stored in "connectivities" + if self.kind == "correlation": + covariances_std = [ + self.cov_estimator_.fit( + signal.standardize_signal( + x, + detrend=False, + standardize=self.standardize, + ) + ).covariance_ + for x in X + ] + connectivities = [cov_to_corr(cov) for cov in covariances_std] + else: + covariances = [self.cov_estimator_.fit(x).covariance_ for x in X] + if self.kind in ("covariance", "tangent"): + connectivities = covariances + elif self.kind == "precision": + connectivities = [linalg.inv(cov) for cov in covariances] + elif self.kind == "partial correlation": + connectivities = [ + prec_to_partial(linalg.inv(cov)) for cov in covariances + ] + else: + allowed_kinds = ( + "correlation", + "partial correlation", + "tangent", + "covariance", + "precision", + ) + raise_error( + f"Allowed connectivity kinds are {allowed_kinds}. " + f"Got kind {self.kind}." + ) + + # Store the mean + if do_fit: + if self.kind == "tangent": + self.mean_ = _geometric_mean( + covariances, max_iter=30, tol=1e-7 + ) + self.whitening_ = _map_eigenvalues( + lambda x: 1.0 / np.sqrt(x), self.mean_ + ) + else: + self.mean_ = np.mean(connectivities, axis=0) + # Fight numerical instabilities: make symmetric + self.mean_ = self.mean_ + self.mean_.T + self.mean_ *= 0.5 + + # Compute the vector we return on transform + if do_transform: + if self.kind == "tangent": + connectivities = [ + _map_eigenvalues( + np.log, self.whitening_.dot(cov).dot(self.whitening_) + ) + for cov in connectivities + ] + + connectivities = np.array(connectivities) + + if confounds is not None and not self.vectorize: + error_message = ( + "'confounds' are provided but vectorize=False. " + "Confounds are only cleaned on vectorized matrices " + "as second level connectome regression " + "but not on symmetric matrices." + ) + raise_error(error_message) + + if self.vectorize: + connectivities = sym_matrix_to_vec( + connectivities, discard_diagonal=self.discard_diagonal + ) + if confounds is not None: + connectivities = signal.clean( + connectivities, confounds=confounds + ) + + return connectivities diff --git a/junifer/external/nilearn/tests/test_junifer_connectivity_measure.py b/junifer/external/nilearn/tests/test_junifer_connectivity_measure.py new file mode 100644 index 000000000..14dc8a07c --- /dev/null +++ b/junifer/external/nilearn/tests/test_junifer_connectivity_measure.py @@ -0,0 +1,1089 @@ +"""Provide tests for JuniferConnectivityMeasure class.""" + +# Authors: Synchon Mandal +# License: AGPL + +import copy +import warnings +from math import cosh, exp, log, sinh, sqrt +from typing import TYPE_CHECKING, List, Optional, Tuple, Type, Union + +import numpy as np +import pytest +from nilearn.connectome.connectivity_matrices import sym_matrix_to_vec +from nilearn.tests.test_signal import generate_signals +from numpy.testing import assert_array_almost_equal, assert_array_equal +from pandas import DataFrame +from scipy import linalg +from sklearn.covariance import EmpiricalCovariance, LedoitWolf + +from junifer.external.nilearn import JuniferConnectivityMeasure +from junifer.external.nilearn.junifer_connectivity_measure import ( + _check_spd, + _check_square, + _form_symmetric, + _geometric_mean, + _map_eigenvalues, + is_spd, +) + + +if TYPE_CHECKING: + from numpy.typing import ArrayLike + from sklearn.base import BaseEstimator + +# New BSD License + +# Copyright (c) The nilearn developers. +# All rights reserved. + + +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: + +# a. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# b. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# c. Neither the name of the nilearn developers nor the names of +# its contributors may be used to endorse or promote products +# derived from this software without specific prior written +# permission. + + +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR +# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY +# OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH +# DAMAGE. + + +CONNECTIVITY_KINDS = ( + "covariance", + "correlation", + "tangent", + "precision", + "partial correlation", +) + +N_FEATURES = 49 + +N_SUBJECTS = 5 + + +def random_diagonal( + p: int, + v_min: float = 1.0, + v_max: float = 2.0, + random_state: Union[int, np.random.RandomState] = 0, +) -> np.ndarray: + """Generate a random diagonal matrix. + + Parameters + ---------- + p : int + The first dimension of the array. + v_min : float, optional + Minimal element (default 1.). + v_max : float, optional + Maximal element (default 2.). + random_state : int or numpy.random.RandomState instance, optional + random number generator, or seed (default 0). + + Returns + ------- + numpy.ndarray of shape (p, p) + A diagonal matrix with the given minimal and maximal elements. + + """ + random_state = np.random.default_rng(random_state) + diag = random_state.random(p) * (v_max - v_min) + v_min + diag[diag == np.amax(diag)] = v_max + diag[diag == np.amin(diag)] = v_min + return np.diag(diag) + + +def random_spd( + p: int, + eig_min: float, + cond: float, + random_state: Union[int, np.random.RandomState] = 0, +) -> np.ndarray: + """Generate a random symmetric positive definite matrix. + + Parameters + ---------- + p : int + The first dimension of the array. + eig_min : float + Minimal eigenvalue. + cond : float + Condition number, defined as the ratio of the maximum eigenvalue to the + minimum one. + random_state : int or numpy.random.RandomState instance, optional + random number generator, or seed (default 0). + + Returns + ------- + numpy.ndarray of shape (p, p) + A symmetric positive definite matrix with the given minimal eigenvalue + and condition number. + + """ + rand_gen = np.random.default_rng(random_state) + mat = rand_gen.standard_normal((p, p)) + unitary, _ = linalg.qr(mat) + diag = random_diagonal( + p, v_min=eig_min, v_max=cond * eig_min, random_state=random_state + ) + return unitary.dot(diag).dot(unitary.T) + + +def _signals( + n_subjects: int = N_SUBJECTS, +) -> Tuple[List[np.ndarray], np.ndarray]: + """Generate signals and compute covariances while applying confounds. + + Parameters + ---------- + n_subjects : int + Number of subjects. + + Returns + ------- + tuple of list of np.ndarray and np.ndarray. + + """ + n_features = N_FEATURES + signals = [] + for k in range(n_subjects): + n_samples = 200 + k + signal, _, confounds = generate_signals( + n_features=n_features, + n_confounds=5, + length=n_samples, + same_variance=False, + ) + signals.append(signal) + signal -= signal.mean(axis=0) + return signals, confounds + + +@pytest.fixture +def signals() -> List[np.ndarray]: + """Return signals as list of np.ndarray.""" + return _signals(N_SUBJECTS)[0] + + +@pytest.fixture +def signals_and_covariances( + cov_estimator: Union[LedoitWolf, EmpiricalCovariance] +) -> Tuple[List[np.ndarray], List[float]]: + """Return signals and covariances for a covariance estimator. + + Parameters + ---------- + cov_estimator : LedoitWolf instance or EmpiricalCovariance instance + The covariance estimator. + + Returns + ------- + tuple of list of np.ndarray and list of float + + """ + signals, _ = _signals() + emp_covs = [] + ledoit_covs = [] + ledoit_estimator = LedoitWolf() + for k, signal_ in enumerate(signals): + n_samples = 200 + k + signal_ -= signal_.mean(axis=0) + emp_covs.append((signal_.T).dot(signal_) / n_samples) + ledoit_covs.append(ledoit_estimator.fit(signal_).covariance_) + + if isinstance(cov_estimator, LedoitWolf): + return signals, ledoit_covs + elif isinstance(cov_estimator, EmpiricalCovariance): + return signals, emp_covs + + +def test_check_square() -> None: + """Test square matrix assertion.""" + non_square = np.ones((2, 3)) + with pytest.raises(ValueError, match="Expected a square matrix"): + _check_square(non_square) + + +@pytest.mark.parametrize( + "invalid_input", + [ + np.array([[0, 1], [0, 0]]), # non symmetric + np.ones((3, 3)), # non SPD + ], +) +def test_check_spd(invalid_input: np.ndarray) -> None: + """Test matrix is symmetric positive definite. + + Parameters + ---------- + invalid_input : numpy.ndarray + The parametrized invalid input array. + + """ + with pytest.raises( + ValueError, match="Expected a symmetric positive definite matrix." + ): + _check_spd(invalid_input) + + +def test_map_eigenvalues_on_exp_map() -> None: + """Test exponential eigenvalues mapping.""" + sym = np.ones((2, 2)) + sym_exp = exp(1.0) * np.array( + [[cosh(1.0), sinh(1.0)], [sinh(1.0), cosh(1.0)]] + ) + assert_array_almost_equal(_map_eigenvalues(np.exp, sym), sym_exp) + + +def test_map_eigenvalues_on_sqrt_map() -> None: + """Test square-root eigenvalues mapping.""" + spd_sqrt = np.array( + [[2.0, -1.0, 0.0], [-1.0, 2.0, -1.0], [0.0, -1.0, 2.0]] + ) + spd = spd_sqrt.dot(spd_sqrt) + assert_array_almost_equal(_map_eigenvalues(np.sqrt, spd), spd_sqrt) + + +def test_map_eigenvalues_on_log_map() -> None: + """Test logarithmic eigenvalues mapping.""" + spd = np.array([[1.25, 0.75], [0.75, 1.25]]) + spd_log = np.array([[0.0, log(2.0)], [log(2.0), 0.0]]) + assert_array_almost_equal(_map_eigenvalues(np.log, spd), spd_log) + + +def test_geometric_mean_couple() -> None: + """Test geometric mean.""" + n_features = 7 + spd1 = np.ones((n_features, n_features)) + spd1 = spd1.dot(spd1) + n_features * np.eye(n_features) + spd2 = np.tril(np.ones((n_features, n_features))) + spd2 = spd2.dot(spd2.T) + vals_spd2, vecs_spd2 = np.linalg.eigh(spd2) + spd2_sqrt = _form_symmetric(np.sqrt, vals_spd2, vecs_spd2) + spd2_inv_sqrt = _form_symmetric(np.sqrt, 1.0 / vals_spd2, vecs_spd2) + geo = spd2_sqrt.dot( + _map_eigenvalues(np.sqrt, spd2_inv_sqrt.dot(spd1).dot(spd2_inv_sqrt)) + ).dot(spd2_sqrt) + + assert_array_almost_equal(_geometric_mean([spd1, spd2]), geo) + + +def test_geometric_mean_diagonal() -> None: + """Test geometric mean along diagonal.""" + n_matrices = 20 + n_features = 5 + diags = [] + for k in range(n_matrices): + diag = np.eye(n_features) + diag[k % n_features, k % n_features] = 1e4 + k + diag[(n_features - 1) // (k + 1), (n_features - 1) // (k + 1)] = ( + k + 1 + ) * 1e-4 + diags.append(diag) + geo = np.prod(np.array(diags), axis=0) ** (1 / float(len(diags))) + + assert_array_almost_equal(_geometric_mean(diags), geo) + + +def test_geometric_mean_geodesic() -> None: + """Test geometric mean along geodesic.""" + n_matrices = 10 + n_features = 6 + sym = np.arange(n_features) / np.linalg.norm(np.arange(n_features)) + sym = sym * sym[:, np.newaxis] + times = np.arange(n_matrices) + non_singular = np.eye(n_features) + non_singular[1:3, 1:3] = np.array([[-1, -0.5], [-0.5, -1]]) + spds = [ + non_singular.dot(_map_eigenvalues(np.exp, time * sym)).dot( + non_singular.T + ) + for time in times + ] + gmean = non_singular.dot(_map_eigenvalues(np.exp, times.mean() * sym)).dot( + non_singular.T + ) + assert_array_almost_equal(_geometric_mean(spds), gmean) + + +def test_geometric_mean_properties() -> None: + """Test geometric mean properties.""" + n_matrices = 40 + n_features = 15 + spds = [ + random_spd(n_features, eig_min=1.0, cond=10.0, random_state=0) + for _ in range(n_matrices) + ] + input_spds = copy.copy(spds) + + gmean = _geometric_mean(spds) + + # Generic + assert isinstance(spds, list) + for spd, input_spd in zip(spds, input_spds): + assert_array_equal(spd, input_spd) + assert is_spd(gmean, decimal=7) + + +def random_non_singular( + p: int, + sing_min: float = 1.0, + sing_max: float = 2.0, + random_state: Union[int, np.random.RandomState] = 0, +) -> np.ndarray: + """Generate a random nonsingular matrix. + + Parameters + ---------- + p : int + The first dimension of the array. + sing_min : float, optional + Minimal singular value (default 1.). + sing_max : float, optional + Maximal singular value (default 2.). + random_state : int or numpy.random.RandomState instance, optional + random number generator, or seed (default 0). + + Returns + ------- + numpy.ndarray of shape (p, p) + A nonsingular matrix with the given minimal and maximal singular + values. + + """ + rand_gen = np.random.default_rng(random_state) + diag = random_diagonal( + p, v_min=sing_min, v_max=sing_max, random_state=random_state + ) + mat1 = rand_gen.standard_normal((p, p)) + mat2 = rand_gen.standard_normal((p, p)) + unitary1, _ = linalg.qr(mat1) + unitary2, _ = linalg.qr(mat2) + return unitary1.dot(diag).dot(unitary2.T) + + +def test_geometric_mean_properties_check_invariance() -> None: + """Test geometric mean properties' invariance.""" + n_matrices = 40 + n_features = 15 + spds = [ + random_spd(n_features, eig_min=1.0, cond=10.0, random_state=0) + for _ in range(n_matrices) + ] + + gmean = _geometric_mean(spds) + + # Invariance under reordering + spds.reverse() + spds.insert(0, spds[1]) + spds.pop(2) + assert_array_almost_equal(_geometric_mean(spds), gmean) + + # Invariance under congruent transformation + non_singular = random_non_singular(n_features, random_state=0) + spds_cong = [non_singular.dot(spd).dot(non_singular.T) for spd in spds] + assert_array_almost_equal( + _geometric_mean(spds_cong), non_singular.dot(gmean).dot(non_singular.T) + ) + + # Invariance under inversion + spds_inv = [linalg.inv(spd) for spd in spds] + init = linalg.inv(np.mean(spds, axis=0)) + assert_array_almost_equal( + _geometric_mean(spds_inv, init=init), linalg.inv(gmean) + ) + + +def grad_geometric_mean( + mats: "ArrayLike", + init: Optional["ArrayLike"] = None, + max_iter: int = 10, + tol: float = 1e-7, +) -> List[float]: + """Compute gradient of geometric mean. + + Return the norm of the covariant derivative at each iteration step + of geometric_mean. See its docstring for details. + + Norm is intrinsic norm on the tangent space of the manifold of symmetric + positive definite matrices. + + Parameters + ---------- + mats : array-like object + Object that can be converted to np.ndarray. + init : array-like object or None, optional + Initialization matrix (default None). + max_iter : int, optional + Maximum iteration for gradient descent (default 10). + tol : float, optional + Tolerance for norm (default 1e-7). + + Returns + ------- + grad_norm : list of float + Norm of the covariant derivative in the tangent space at each step. + + """ + mats = np.array(mats) + + # Initialization + gmean = init or np.mean(mats, axis=0) + + norm_old = np.inf + step = 1.0 + grad_norm = [] + for _ in range(max_iter): + # Computation of the gradient + vals_gmean, vecs_gmean = linalg.eigh(gmean) + gmean_inv_sqrt = _form_symmetric(np.sqrt, 1.0 / vals_gmean, vecs_gmean) + whitened_mats = [ + gmean_inv_sqrt.dot(mat).dot(gmean_inv_sqrt) for mat in mats + ] + logs = [_map_eigenvalues(np.log, w_mat) for w_mat in whitened_mats] + + # Covariant derivative is - gmean.dot(logs_mean) + logs_mean = np.mean(logs, axis=0) + + # Norm of the covariant derivative on + # the tangent space at point gmean + norm = np.linalg.norm(logs_mean) + + # Update of the minimizer + vals_log, vecs_log = linalg.eigh(logs_mean) + gmean_sqrt = _form_symmetric(np.sqrt, vals_gmean, vecs_gmean) + gmean = gmean_sqrt.dot( + _form_symmetric(np.exp, vals_log * step, vecs_log) + ).dot(gmean_sqrt) + + # Update the norm and the step size + if norm < norm_old: + norm_old = norm + if norm > norm_old: + step = step / 2.0 + norm = norm_old + + grad_norm.append(norm / gmean.size) + if tol is not None and norm / gmean.size < tol: + break + + return grad_norm + + +def test_geometric_mean_properties_check_gradient() -> None: + """Test geometric mean properties' gradient.""" + n_matrices = 40 + n_features = 15 + spds = [ + random_spd(n_features, eig_min=1.0, cond=10.0, random_state=0) + for _ in range(n_matrices) + ] + + grad_norm = grad_geometric_mean(spds, tol=1e-20) + + # Gradient norm is decreasing + difference = np.diff(grad_norm) + assert np.amax(difference) <= 0.0 + + # Check warning if gradient norm in the last step is less than + # tolerance + max_iter = 1 + tol = 1e-20 + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + _geometric_mean(spds, max_iter=max_iter, tol=tol) + assert len(w) == 1 + + grad_norm = grad_geometric_mean(spds, max_iter=max_iter, tol=tol) + + assert len(grad_norm) == max_iter + assert grad_norm[-1] > tol + + +# proportion of badly conditioned matrices +@pytest.mark.parametrize("p", [0.5, 1.0]) +def test_geometric_mean_properties_evaluate_convergence(p: float) -> None: + """Test geometric mean properties' convergence. + + Parameters + ---------- + p : float + Convergence criteria. + + """ + n_matrices = 40 + n_features = 15 + # A warning is printed if tolerance is not reached + spds = [ + random_spd(n_features, eig_min=1e-2, cond=1e6, random_state=0) + for _ in range(int(p * n_matrices)) + ] + spds.extend( + random_spd(n_features, eig_min=1.0, cond=10.0, random_state=0) + for _ in range(int(p * n_matrices), n_matrices) + ) + max_iter = 30 if p < 1 else 60 + + _geometric_mean(spds, max_iter=max_iter, tol=1e-5) + + +def test_geometric_mean_error_non_square_matrix() -> None: + """Test geometric mean error for non-square matrix.""" + n_features = 5 + mat1 = np.ones((n_features, n_features + 1)) + + with pytest.raises(ValueError, match="Expected a square matrix"): + _geometric_mean([mat1]) + + +def test_geometric_mean_error_input_matrices_have_different_shapes() -> None: + """Test geometric mean error for different input matrices shape.""" + n_features = 5 + mat1 = np.eye(n_features) + mat2 = np.ones((n_features + 1, n_features + 1)) + + with pytest.raises( + ValueError, match="Matrices are not of the same shape." + ): + _geometric_mean([mat1, mat2]) + + +def test_geometric_mean_error_non_spd_input_matrix() -> None: + """Test geometric mean error for non SPD input matrix.""" + n_features = 5 + mat2 = np.ones((n_features + 1, n_features + 1)) + + with pytest.raises( + ValueError, match="Expected a symmetric positive definite matrix." + ): + _geometric_mean([mat2]) + + +def test_connectivity_measure_errors(): + """Test errors.""" + # Raising error for input subjects not iterable + conn_measure = JuniferConnectivityMeasure() + + with pytest.raises( + ValueError, match="'subjects' input argument must be an iterable" + ): + conn_measure.fit(1.0) + + # input subjects not 2D numpy.ndarrays + with pytest.raises( + ValueError, match="Each subject must be 2D numpy.ndarray." + ): + conn_measure.fit([np.ones((100, 40)), np.ones((10,))]) + + # input subjects with different number of features + with pytest.raises( + ValueError, match="All subjects must have the same number of features." + ): + conn_measure.fit([np.ones((100, 40)), np.ones((100, 41))]) + + # fit_transform with a single subject and kind=tangent + conn_measure = JuniferConnectivityMeasure(kind="tangent") + + with pytest.raises( + ValueError, + match="Tangent space parametrization .* only be .* group of subjects", + ): + conn_measure.fit_transform([np.ones((100, 40))]) + + +@pytest.mark.parametrize( + "cov_estimator", [EmpiricalCovariance(), LedoitWolf()] +) +@pytest.mark.parametrize("kind", CONNECTIVITY_KINDS) +def test_connectivity_measure_generic( + kind: str, + cov_estimator: Type["BaseEstimator"], + signals_and_covariances: Tuple[List[np.ndarray], List[float]], +) -> None: + """Test generic JuniferConnectivityMeasure. + + Parameters + ---------- + kind : str + The parametrized connectivity matrix kind. + cov_estimator : estimator object + The parametrized covariance estimator. + signals_and_covariances : tuple + The signals and covariances for a covariance estimator. + + """ + signals, covs = signals_and_covariances + + # Check outputs properties + input_covs = copy.copy(covs) + conn_measure = JuniferConnectivityMeasure( + kind=kind, cov_estimator=cov_estimator + ) + connectivities = conn_measure.fit_transform(signals) + + # Generic + assert isinstance(connectivities, np.ndarray) + assert len(connectivities) == len(covs) + + for k, _ in enumerate(connectivities): + assert_array_equal(input_covs[k], covs[k]) + + assert is_spd(covs[k], decimal=7) + + +def _assert_connectivity_tangent(connectivities, conn_measure, covs) -> None: + """Assert tangent connectivity matrix. + + Check output value properties for tangent connectivity measure + that they have the expected relationship + to the input covariance matrices. + + - the geometric mean of the eigenvalues + of the mean covariance matrix is positive-definite + - the whitening matrix (used to transform the data + also produces a positive-definite matrix + + """ + for true_covariance_matrix, estimated_covariance_matrix in zip( + covs, connectivities + ): + assert_array_almost_equal( + estimated_covariance_matrix, estimated_covariance_matrix.T + ) + + assert is_spd(conn_measure.whitening_, decimal=7) + + gmean_sqrt = _map_eigenvalues(np.sqrt, conn_measure.mean_) + assert is_spd(gmean_sqrt, decimal=7) + assert_array_almost_equal( + conn_measure.whitening_.dot(gmean_sqrt), + np.eye(N_FEATURES), + ) + assert_array_almost_equal( + gmean_sqrt.dot( + _map_eigenvalues(np.exp, estimated_covariance_matrix) + ).dot(gmean_sqrt), + true_covariance_matrix, + ) + + +def _assert_connectivity_precision(connectivities, covs) -> None: + """Assert precision connectivity matrix. + + Estimated precision matrix: + - is positive definite + - its product with the true covariance matrix + is close to the identity matrix + + """ + for true_covariance_matrix, estimated_covariance_matrix in zip( + covs, connectivities + ): + assert is_spd(estimated_covariance_matrix, decimal=7) + assert_array_almost_equal( + estimated_covariance_matrix.dot(true_covariance_matrix), + np.eye(N_FEATURES), + ) + + +def _assert_connectivity_correlation( + connectivities, cov_estimator, covs +) -> None: + """Assert correlation connectivity matrix. + + Verify that the estimated covariance matrix: + - is symmetric and positive definite + - has values close to 1 on its diagonal + + If the covariance estimator is EmpiricalCovariance, + the product of: + - the square root of the diagonal of the true covariance matrix + - the estimated covariance matrix + - the square root of the diagonal of the true covariance matrix + + should be close to the true covariance matrix. + + """ + for true_covariance_matrix, estimated_covariance_matrix in zip( + covs, connectivities + ): + assert is_spd(estimated_covariance_matrix, decimal=7) + + assert_array_almost_equal( + np.diag(estimated_covariance_matrix), np.ones(N_FEATURES) + ) + + if cov_estimator == EmpiricalCovariance(): + # square root of the diagonal of the true covariance matrix + d = np.sqrt(np.diag(np.diag(true_covariance_matrix))) + + assert_array_almost_equal( + d.dot(estimated_covariance_matrix).dot(d), + true_covariance_matrix, + ) + + +def _assert_connectivity_partial_correlation(connectivities, covs) -> None: + """Assert partial correlation connectivity matrix.""" + for true_covariance_matrix, estimated_covariance_matrix in zip( + covs, connectivities + ): + precision_matrix = linalg.inv(true_covariance_matrix) + + # square root of the diagonal elements of the precision matrix + d = np.sqrt(np.diag(np.diag(precision_matrix))) + + # normalize the computed partial correlation matrix + # necessary to ensure that the diagonal elements + # of the partial correlation matrix are equal to 1 + normalized_partial_correlation_matrix = d.dot( + estimated_covariance_matrix + ).dot(d) + + # expected value + partial_corrlelation_matrix = -precision_matrix + 2 * np.diag( + np.diag(precision_matrix) + ) + + assert_array_almost_equal( + normalized_partial_correlation_matrix, + partial_corrlelation_matrix, + ) + + +@pytest.mark.parametrize( + "kind", + ["tangent", "precision", "correlation", "partial correlation"], +) +@pytest.mark.parametrize( + "cov_estimator", [EmpiricalCovariance(), LedoitWolf()] +) +def test_connectivity_measure_specific_for_each_kind( + kind: str, + cov_estimator: Type["BaseEstimator"], + signals_and_covariances: Tuple[List[np.ndarray], List[float]], +) -> None: + """Test connectivity matrix for each kind. + + Parameters + ---------- + kind : str + The parametrized connectivity matrix kind. + cov_estimator : estimator object + The parametrized covariance estimator. + signals_and_covariances : tuple + The signals and covariances for a covariance estimator. + + """ + signals, covs = signals_and_covariances + + conn_measure = JuniferConnectivityMeasure( + kind=kind, cov_estimator=cov_estimator + ) + connectivities = conn_measure.fit_transform(signals) + + if kind == "tangent": + _assert_connectivity_tangent(connectivities, conn_measure, covs) + elif kind == "precision": + _assert_connectivity_precision(connectivities, covs) + elif kind == "correlation": + _assert_connectivity_correlation(connectivities, cov_estimator, covs) + elif kind == "partial correlation": + _assert_connectivity_partial_correlation(connectivities, covs) + + +@pytest.mark.parametrize("kind", CONNECTIVITY_KINDS) +def test_connectivity_measure_check_mean( + kind: str, signals: List[np.ndarray] +) -> None: + """Test mean of connectivity matrix for each kind. + + Parameters + ---------- + kind : str + The parametrized connectivity matrix kind. + signals : list of np.ndarray + The input signals. + + """ + conn_measure = JuniferConnectivityMeasure(kind=kind) + conn_measure.fit_transform(signals) + + assert (conn_measure.mean_).shape == (N_FEATURES, N_FEATURES) + + if kind != "tangent": + assert_array_almost_equal( + conn_measure.mean_, + np.mean(conn_measure.transform(signals), axis=0), + ) + + # Check that the mean isn't modified in transform + conn_measure = JuniferConnectivityMeasure(kind="covariance") + conn_measure.fit(signals[:1]) + mean = conn_measure.mean_ + conn_measure.transform(signals[1:]) + + assert_array_equal(mean, conn_measure.mean_) + + +@pytest.mark.parametrize("kind", CONNECTIVITY_KINDS) +def test_connectivity_measure_check_vectorization_option( + kind: str, signals: List[np.ndarray] +) -> None: + """Test vectorization of connectivity matrix for each kind. + + Parameters + ---------- + kind : str + The parametrized connectivity matrix kind. + signals : list of np.ndarray + The input signals. + + """ + conn_measure = JuniferConnectivityMeasure(kind=kind) + connectivities = conn_measure.fit_transform(signals) + conn_measure = JuniferConnectivityMeasure(vectorize=True, kind=kind) + vectorized_connectivities = conn_measure.fit_transform(signals) + + assert_array_almost_equal( + vectorized_connectivities, sym_matrix_to_vec(connectivities) + ) + + # Check not fitted error + with pytest.raises(ValueError, match="has not been fitted. "): + JuniferConnectivityMeasure().inverse_transform( + vectorized_connectivities + ) + + +@pytest.mark.parametrize( + "kind", + ["covariance", "correlation", "precision", "partial correlation"], +) +def test_connectivity_measure_check_inverse_transformation( + kind: str, signals: List[np.ndarray] +) -> None: + """Test inverse transform. + + Parameters + ---------- + kind : str + The parametrized connectivity matrix kind. + signals : list of np.ndarray + The input signals. + + """ + # without vectorization: input matrices are returned with no change + conn_measure = JuniferConnectivityMeasure(kind=kind) + connectivities = conn_measure.fit_transform(signals) + + assert_array_almost_equal( + conn_measure.inverse_transform(connectivities), connectivities + ) + + # with vectorization: input vectors are reshaped into matrices + # if diagonal has not been discarded + conn_measure = JuniferConnectivityMeasure(kind=kind, vectorize=True) + vectorized_connectivities = conn_measure.fit_transform(signals) + + assert_array_almost_equal( + conn_measure.inverse_transform(vectorized_connectivities), + connectivities, + ) + + +@pytest.mark.parametrize( + "kind", + ["covariance", "correlation", "precision", "partial correlation"], +) +def test_connectivity_measure_check_inverse_transformation_discard_diag( + kind: str, signals: List[np.ndarray] +) -> None: + """Test diagonal for inverse transform. + + Parameters + ---------- + kind : str + The parametrized connectivity matrix kind. + signals : list of np.ndarray + The input signals. + + """ + # with vectorization + connectivities = JuniferConnectivityMeasure(kind=kind).fit_transform( + signals + ) + conn_measure = JuniferConnectivityMeasure( + kind=kind, vectorize=True, discard_diagonal=True + ) + vectorized_connectivities = conn_measure.fit_transform(signals) + + if kind in ["correlation", "partial correlation"]: + assert_array_almost_equal( + conn_measure.inverse_transform(vectorized_connectivities), + connectivities, + ) + elif kind in ["covariance", "precision"]: + diagonal = np.array( + [np.diagonal(conn) / sqrt(2) for conn in connectivities] + ) + inverse_transformed = conn_measure.inverse_transform( + vectorized_connectivities, diagonal=diagonal + ) + + assert_array_almost_equal(inverse_transformed, connectivities) + with pytest.raises( + ValueError, match="cannot reconstruct connectivity matrices" + ): + conn_measure.inverse_transform(vectorized_connectivities) + + +def test_connectivity_measure_inverse_transform_tangent( + signals: List[np.ndarray], +) -> None: + """Test that for 'tangent' kind, covariance matrices are reconstructed. + + Parameters + ---------- + signals : list of np.ndarray + The input signals. + + """ + # Without vectorization + tangent_measure = JuniferConnectivityMeasure(kind="tangent") + displacements = tangent_measure.fit_transform(signals) + covariances = JuniferConnectivityMeasure(kind="covariance").fit_transform( + signals + ) + + assert_array_almost_equal( + tangent_measure.inverse_transform(displacements), covariances + ) + + # with vectorization + # when diagonal has not been discarded + tangent_measure = JuniferConnectivityMeasure( + kind="tangent", vectorize=True + ) + vectorized_displacements = tangent_measure.fit_transform(signals) + + assert_array_almost_equal( + tangent_measure.inverse_transform(vectorized_displacements), + covariances, + ) + + # When diagonal has been discarded + tangent_measure = JuniferConnectivityMeasure( + kind="tangent", vectorize=True, discard_diagonal=True + ) + vectorized_displacements = tangent_measure.fit_transform(signals) + + diagonal = np.array( + [np.diagonal(matrix) / sqrt(2) for matrix in displacements] + ) + inverse_transformed = tangent_measure.inverse_transform( + vectorized_displacements, diagonal=diagonal + ) + + assert_array_almost_equal(inverse_transformed, covariances) + with pytest.raises( + ValueError, match="cannot reconstruct connectivity matrices" + ): + tangent_measure.inverse_transform(vectorized_displacements) + + +def test_confounds_connectivity_measure() -> None: + """Test confounds.""" + n_subjects = 10 + + signals, confounds = _signals(n_subjects) + + correlation_measure = JuniferConnectivityMeasure( + kind="correlation", vectorize=True + ) + + # Clean confounds on 10 subjects with confounds filtered to 10 subjects in + # length + cleaned_vectors = correlation_measure.fit_transform( + signals, confounds=confounds[:10] + ) + + zero_matrix = np.zeros((confounds.shape[1], cleaned_vectors.shape[1])) + assert_array_almost_equal( + np.dot(confounds[:10].T, cleaned_vectors), zero_matrix + ) + assert isinstance(cleaned_vectors, np.ndarray) + + # Confounds as pandas DataFrame + confounds_df = DataFrame(confounds[:10]) + correlation_measure.fit_transform(signals, confounds=confounds_df) + + +def test_confounds_connectivity_measure_errors() -> None: + """Test errors for dealing with confounds.""" + # Generate signals and compute covariances and apply confounds while + # computing covariances + signals, confounds = _signals() + + # Raising error for input confounds are not iterable + conn_measure = JuniferConnectivityMeasure(vectorize=True) + msg = "'confounds' input argument must be an iterable" + + with pytest.raises(ValueError, match=msg): + conn_measure._check_input(X=signals, confounds=1.0) + + with pytest.raises(ValueError, match=msg): + conn_measure._fit_transform( + X=signals, do_fit=True, do_transform=True, confounds=1.0 + ) + + with pytest.raises(ValueError, match=msg): + conn_measure.fit_transform(X=signals, y=None, confounds=1.0) + + # Raising error for input confounds are given but not vectorize=True + conn_measure = JuniferConnectivityMeasure(vectorize=False) + with pytest.raises( + ValueError, match="'confounds' are provided but vectorize=False" + ): + conn_measure.fit_transform(signals, None, confounds[:10]) + + +def test_connectivity_measure_standardize( + signals: List[np.ndarray], +) -> None: + """Check warning is raised and then suppressed with setting standardize. + + Parameters + ---------- + signals : list of np.ndarray + The input signals. + + """ + match = "default strategy for standardize" + + with pytest.warns(DeprecationWarning, match=match): + JuniferConnectivityMeasure(kind="correlation").fit_transform(signals) + + with warnings.catch_warnings(record=True) as record: + JuniferConnectivityMeasure( + kind="correlation", standardize="zscore_sample" + ).fit_transform(signals) + for m in record: + assert match not in m.message diff --git a/junifer/markers/functional_connectivity/crossparcellation_functional_connectivity.py b/junifer/markers/functional_connectivity/crossparcellation_functional_connectivity.py index 95fdbd653..0dc154fc7 100644 --- a/junifer/markers/functional_connectivity/crossparcellation_functional_connectivity.py +++ b/junifer/markers/functional_connectivity/crossparcellation_functional_connectivity.py @@ -28,18 +28,24 @@ class CrossParcellationFC(BaseMarker): The name of the first parcellation. parcellation_two : str The name of the second parcellation. - aggregation_method : str, optional - The aggregation method (default "mean"). - correlation_method : str, optional + 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). + corr_method : str, optional Any method that can be passed to - :any:`pandas.DataFrame.corr` (default "pearson"). + :meth:`pandas.DataFrame.corr` (default "pearson"). 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 ` for more details. If None, will not apply any mask (default None). name : str, optional - The name of the marker. If None, will use the class name - (default None). + The name of the marker. If None, will use + ``BOLD_CrossParcellationFC`` (default None). """ @@ -55,8 +61,9 @@ class CrossParcellationFC(BaseMarker): self, parcellation_one: str, parcellation_two: str, - aggregation_method: str = "mean", - correlation_method: str = "pearson", + agg_method: str = "mean", + agg_method_params: Optional[Dict] = None, + corr_method: str = "pearson", masks: Union[str, Dict, List[Union[Dict, str]], None] = None, name: Optional[str] = None, ) -> None: @@ -66,8 +73,9 @@ class CrossParcellationFC(BaseMarker): ) self.parcellation_one = parcellation_one self.parcellation_two = parcellation_two - self.aggregation_method = aggregation_method - self.correlation_method = correlation_method + self.agg_method = agg_method + self.agg_method_params = agg_method_params + self.corr_method = corr_method self.masks = masks super().__init__(on=["BOLD"], name=name) @@ -115,13 +123,17 @@ class CrossParcellationFC(BaseMarker): # Perform aggregation using two parcellations aggregation_parcellation_one = ParcelAggregation( parcellation=self.parcellation_one, - method=self.aggregation_method, + method=self.agg_method, + method_params=self.agg_method_params, masks=self.masks, + on="BOLD", ).compute(input, extra_input=extra_input) aggregation_parcellation_two = ParcelAggregation( parcellation=self.parcellation_two, - method=self.aggregation_method, + method=self.agg_method, + method_params=self.agg_method_params, masks=self.masks, + on="BOLD", ).compute(input, extra_input=extra_input) return { @@ -133,7 +145,7 @@ class CrossParcellationFC(BaseMarker): pd.DataFrame( aggregation_parcellation_two["aggregation"]["data"] ), - method=self.correlation_method, + method=self.corr_method, ).values, # Columns should be named after parcellation 1 "col_names": aggregation_parcellation_one["aggregation"][ diff --git a/junifer/markers/functional_connectivity/edge_functional_connectivity_parcels.py b/junifer/markers/functional_connectivity/edge_functional_connectivity_parcels.py index a638837b6..069f2244a 100644 --- a/junifer/markers/functional_connectivity/edge_functional_connectivity_parcels.py +++ b/junifer/markers/functional_connectivity/edge_functional_connectivity_parcels.py @@ -22,36 +22,40 @@ class EdgeCentricFCParcels(FunctionalConnectivityBase): Parameters ---------- parcellation : str or list of str - The name(s) of the parcellation(s). Check valid options by calling - :func:`.list_parcellations`. + The name(s) of the parcellation(s) to use. + See :func:`.list_parcellations` for options. agg_method : str, optional - The method to perform aggregation of BOLD time series. - Check valid options in :func:`.get_aggfunc_by_name` + 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. Check valid options in - :func:`.get_aggfunc_by_name` (default None). - cor_method : str, optional - The method to perform correlation. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` - (default "covariance"). - cor_method_params : dict, optional - Parameters to pass to the correlation function. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` (default None). + Parameters to pass to the aggregation function. + See :func:`.get_aggfunc_by_name` for options + (default None). + conn_method : str, optional + The method to perform connectivity measure using. + See :class:`.JuniferConnectivityMeasure` for options + (default "correlation"). + conn_method_params : dict, optional + Parameters to pass to :class:`.JuniferConnectivityMeasure`. + If None, ``{"empirical": True}`` will be used, which would mean + :class:`sklearn.covariance.EmpiricalCovariance` is used to compute + covariance. If usage of :class:`sklearn.covariance.LedoitWolf` is + desired, ``{"empirical": False}`` should be passed + (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 ` for more details. If None, will not apply any mask (default None). name : str, optional - The name of the marker. If None, will use the class name (default - None). + The name of the marker. If None, will use + ``BOLD_EdgeCentricFCParcels`` (default None). References ---------- .. [1] Jo et al. (2021) - Subject identification using - edge-centric functional connectivity - doi: https://doi.org/10.1016/j.neuroimage.2021.118204 + Subject identification using edge-centric functional connectivity. + https://doi.org/10.1016/j.neuroimage.2021.118204 """ @@ -60,8 +64,8 @@ class EdgeCentricFCParcels(FunctionalConnectivityBase): parcellation: Union[str, List[str]], agg_method: str = "mean", agg_method_params: Optional[Dict] = None, - cor_method: str = "covariance", - cor_method_params: Optional[Dict] = None, + conn_method: str = "correlation", + conn_method_params: Optional[Dict] = None, masks: Union[str, Dict, List[Union[Dict, str]], None] = None, name: Optional[str] = None, ) -> None: @@ -69,8 +73,8 @@ class EdgeCentricFCParcels(FunctionalConnectivityBase): super().__init__( agg_method=agg_method, agg_method_params=agg_method_params, - cor_method=cor_method, - cor_method_params=cor_method_params, + conn_method=conn_method, + conn_method_params=conn_method_params, masks=masks, name=name, ) diff --git a/junifer/markers/functional_connectivity/edge_functional_connectivity_spheres.py b/junifer/markers/functional_connectivity/edge_functional_connectivity_spheres.py index 169d0fbd4..6d6ebb323 100644 --- a/junifer/markers/functional_connectivity/edge_functional_connectivity_spheres.py +++ b/junifer/markers/functional_connectivity/edge_functional_connectivity_spheres.py @@ -22,42 +22,48 @@ class EdgeCentricFCSpheres(FunctionalConnectivityBase): Parameters ---------- coords : str - The name of the coordinates list to use. See - :func:`.list_coordinates` for options. - radius : float, optional - The radius of the sphere in mm. If None, the signal will be extracted - from a single voxel. See :class:`nilearn.maskers.NiftiSpheresMasker` - for more information (default None). + The name of the coordinates list to use. + See :func:`.list_coordinates` for options. + radius : positive float, optional + The radius of the sphere around each coordinates in millimetres. + If None, the signal will be extracted from a single voxel. + See :class:`.JuniferNiftiSpheresMasker` for more information + (default None). allow_overlap : bool, optional Whether to allow overlapping spheres. If False, an error is raised if - the spheres overlap (default is False). + the spheres overlap (default False). agg_method : str, optional - The aggregation method to use. - See :func:`.get_aggfunc_by_name` for more information - (default None). + The method to perform aggregation using. + See :func:`.get_aggfunc_by_name` for options + (default "mean"). agg_method_params : dict, optional - The parameters to pass to the aggregation method (default None). - cor_method : str, optional - The method to perform correlation using. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` (default "covariance"). - cor_method_params : dict, optional - Parameters to pass to the correlation function. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` (default None). + Parameters to pass to the aggregation function. + See :func:`.get_aggfunc_by_name` for options + (default None). + conn_method : str, optional + The method to perform connectivity measure using. + See :class:`.JuniferConnectivityMeasure` for options + (default "correlation"). + conn_method_params : dict, optional + Parameters to pass to :class:`.JuniferConnectivityMeasure`. + If None, ``{"empirical": True}`` will be used, which would mean + :class:`sklearn.covariance.EmpiricalCovariance` is used to compute + covariance. If usage of :class:`sklearn.covariance.LedoitWolf` is + desired, ``{"empirical": False}`` should be passed + (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 ` for more details. If None, will not apply any mask (default None). name : str, optional - The name of the marker. By default, it will use - KIND_EdgeCentricFCSpheres where KIND is the kind of data it - was applied to (default None). + The name of the marker. If None, will use + ``BOLD_EdgeCentricFCSpheres`` (default None). References ---------- .. [1] Jo et al. (2021) - Subject identification using - edge-centric functional connectivity - doi: https://doi.org/10.1016/j.neuroimage.2021.118204 + Subject identification using edge-centric functional connectivity. + https://doi.org/10.1016/j.neuroimage.2021.118204 """ @@ -68,8 +74,8 @@ class EdgeCentricFCSpheres(FunctionalConnectivityBase): allow_overlap: bool = False, agg_method: str = "mean", agg_method_params: Optional[Dict] = None, - cor_method: str = "covariance", - cor_method_params: Optional[Dict] = None, + conn_method: str = "correlation", + conn_method_params: Optional[Dict] = None, masks: Union[str, Dict, List[Union[Dict, str]], None] = None, name: Optional[str] = None, ) -> None: @@ -81,8 +87,8 @@ class EdgeCentricFCSpheres(FunctionalConnectivityBase): super().__init__( agg_method=agg_method, agg_method_params=agg_method_params, - cor_method=cor_method, - cor_method_params=cor_method_params, + conn_method=conn_method, + conn_method_params=conn_method_params, masks=masks, name=name, ) diff --git a/junifer/markers/functional_connectivity/functional_connectivity_base.py b/junifer/markers/functional_connectivity/functional_connectivity_base.py index 3df63122f..3002fa47e 100644 --- a/junifer/markers/functional_connectivity/functional_connectivity_base.py +++ b/junifer/markers/functional_connectivity/functional_connectivity_base.py @@ -7,9 +7,9 @@ from abc import abstractmethod from typing import Any, ClassVar, Dict, List, Optional, Set, Union -from nilearn.connectome import ConnectivityMeasure -from sklearn.covariance import EmpiricalCovariance +from sklearn.covariance import EmpiricalCovariance, LedoitWolf +from ...external.nilearn import JuniferConnectivityMeasure from ...utils import raise_error from ..base import BaseMarker @@ -23,25 +23,31 @@ class FunctionalConnectivityBase(BaseMarker): Parameters ---------- agg_method : str, optional - The method to perform aggregation using. Check valid options in - :func:`.get_aggfunc_by_name` (default "mean"). + 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). - cor_method : str, optional - The method to perform correlation using. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` - (default "covariance"). - cor_method_params : dict, optional - Parameters to pass to the correlation function. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` (default None). + Parameters to pass to the aggregation function. + Check valid options in :func:`.get_aggfunc_by_name` + (default None). + conn_method : str, optional + The method to perform connectivity measure using. + Check valid options in :class:`.JuniferConnectivityMeasure` + (default "correlation"). + conn_method_params : dict, optional + Parameters to pass to :class:`.JuniferConnectivityMeasure`. + If None, ``{"empirical": True}`` will be used, which would mean + :class:`sklearn.covariance.EmpiricalCovariance` is used to compute + covariance. If usage of :class:`sklearn.covariance.LedoitWolf` is + desired, ``{"empirical": False}`` should be passed + (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 ` for more details. If None, will not apply any mask (default None). name : str, optional - The name of the marker. If None, will use the class name (default - None). + The name of the marker. If None, will use ``BOLD_`` + (default None). """ @@ -57,19 +63,18 @@ class FunctionalConnectivityBase(BaseMarker): self, agg_method: str = "mean", agg_method_params: Optional[Dict] = None, - cor_method: str = "covariance", - cor_method_params: Optional[Dict] = None, + conn_method: str = "correlation", + conn_method_params: Optional[Dict] = None, masks: Union[str, Dict, List[Union[Dict, str]], None] = None, name: Optional[str] = None, ) -> None: self.agg_method = agg_method self.agg_method_params = agg_method_params - self.cor_method = cor_method - self.cor_method_params = cor_method_params or {} - - # default to nilearn behavior - self.cor_method_params["empirical"] = self.cor_method_params.get( - "empirical", False + self.conn_method = conn_method + self.conn_method_params = conn_method_params or {} + # Reverse of nilearn behavior + self.conn_method_params["empirical"] = self.conn_method_params.get( + "empirical", True ) self.masks = masks super().__init__(on="BOLD", name=name) @@ -121,14 +126,21 @@ class FunctionalConnectivityBase(BaseMarker): """ # Perform necessary aggregation aggregation = self.aggregate(input, extra_input=extra_input) - # Compute correlation - if self.cor_method_params["empirical"]: - connectivity = ConnectivityMeasure( - cov_estimator=EmpiricalCovariance(), # type: ignore - kind=self.cor_method, - ) + # Set covariance estimator + if self.conn_method_params["empirical"]: + cov_estimator = EmpiricalCovariance(store_precision=False) else: - connectivity = ConnectivityMeasure(kind=self.cor_method) + cov_estimator = LedoitWolf(store_precision=False) + # Compute correlation + connectivity = JuniferConnectivityMeasure( + cov_estimator=cov_estimator, + kind=self.conn_method, + **{ + k: v + for k, v in self.conn_method_params.items() + if k != "empirical" + }, + ) # Create dictionary for output return { "functional_connectivity": { diff --git a/junifer/markers/functional_connectivity/functional_connectivity_parcels.py b/junifer/markers/functional_connectivity/functional_connectivity_parcels.py index 79a321f01..3daa1af28 100644 --- a/junifer/markers/functional_connectivity/functional_connectivity_parcels.py +++ b/junifer/markers/functional_connectivity/functional_connectivity_parcels.py @@ -22,28 +22,34 @@ class FunctionalConnectivityParcels(FunctionalConnectivityBase): Parameters ---------- parcellation : str or list of str - The name(s) of the parcellation(s). Check valid options by calling - :func:`.list_parcellations`. + The name(s) of the parcellation(s) to use. + See :func:`.list_parcellations` for options. agg_method : str, optional - The method to perform aggregation using. Check valid options in - :func:`.get_aggfunc_by_name` (default "mean"). + 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. Check valid options in - :func:`.get_aggfunc_by_name` (default None). - cor_method : str, optional - The method to perform correlation using. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` - (default "covariance"). - cor_method_params : dict, optional - Parameters to pass to the correlation function. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` (default None). + Parameters to pass to the aggregation function. + See :func:`.get_aggfunc_by_name` for options + (default None). + conn_method : str, optional + The method to perform connectivity measure using. + See :class:`.JuniferConnectivityMeasure` for options + (default "correlation"). + conn_method_params : dict, optional + Parameters to pass to :class:`.JuniferConnectivityMeasure`. + If None, ``{"empirical": True}`` will be used, which would mean + :class:`sklearn.covariance.EmpiricalCovariance` is used to compute + covariance. If usage of :class:`sklearn.covariance.LedoitWolf` is + desired, ``{"empirical": False}`` should be passed + (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 ` for more details. If None, will not apply any mask (default None). name : str, optional - The name of the marker. If None, will use the class name (default - None). + The name of the marker. If None, will use + ``BOLD_FunctionalConnectivityParcels`` (default None). """ @@ -52,8 +58,8 @@ class FunctionalConnectivityParcels(FunctionalConnectivityBase): parcellation: Union[str, List[str]], agg_method: str = "mean", agg_method_params: Optional[Dict] = None, - cor_method: str = "covariance", - cor_method_params: Optional[Dict] = None, + conn_method: str = "correlation", + conn_method_params: Optional[Dict] = None, masks: Union[str, Dict, List[Union[Dict, str]], None] = None, name: Optional[str] = None, ) -> None: @@ -61,8 +67,8 @@ class FunctionalConnectivityParcels(FunctionalConnectivityBase): super().__init__( agg_method=agg_method, agg_method_params=agg_method_params, - cor_method=cor_method, - cor_method_params=cor_method_params, + conn_method=conn_method, + conn_method_params=conn_method_params, masks=masks, name=name, ) diff --git a/junifer/markers/functional_connectivity/functional_connectivity_spheres.py b/junifer/markers/functional_connectivity/functional_connectivity_spheres.py index c4fa36560..e865989a9 100644 --- a/junifer/markers/functional_connectivity/functional_connectivity_spheres.py +++ b/junifer/markers/functional_connectivity/functional_connectivity_spheres.py @@ -23,35 +23,42 @@ class FunctionalConnectivitySpheres(FunctionalConnectivityBase): Parameters ---------- coords : str - The name of the coordinates list to use. See - :func:`.list_coordinates` for options. - radius : float, optional - The radius of the sphere in mm. If None, the signal will be extracted - from a single voxel. See :class:`nilearn.maskers.NiftiSpheresMasker` - for more information (default None). + The name of the coordinates list to use. + See :func:`.list_coordinates` for options. + radius : positive float, optional + The radius of the sphere around each coordinates in millimetres. + If None, the signal will be extracted from a single voxel. + See :class:`.JuniferNiftiSpheresMasker` for more information + (default None). allow_overlap : bool, optional Whether to allow overlapping spheres. If False, an error is raised if - the spheres overlap (default is False). + the spheres overlap (default False). agg_method : str, optional - The aggregation method to use. - See :func:`.get_aggfunc_by_name` for more information - (default None). + The method to perform aggregation using. + See :func:`.get_aggfunc_by_name` for options + (default "mean"). agg_method_params : dict, optional - The parameters to pass to the aggregation method (default None). - cor_method : str, optional - The method to perform correlation using. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` (default "covariance"). - cor_method_params : dict, optional - Parameters to pass to the correlation function. Check valid options in - :class:`nilearn.connectome.ConnectivityMeasure` (default None). + Parameters to pass to the aggregation function. + See :func:`.get_aggfunc_by_name` for options + (default None). + conn_method : str, optional + The method to perform connectivity measure using. + See :class:`.JuniferConnectivityMeasure` for options + (default "correlation"). + conn_method_params : dict, optional + Parameters to pass to :class:`.JuniferConnectivityMeasure`. + If None, ``{"empirical": True}`` will be used, which would mean + :class:`sklearn.covariance.EmpiricalCovariance` is used to compute + covariance. If usage of :class:`sklearn.covariance.LedoitWolf` is + desired, ``{"empirical": False}`` should be passed + (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 ` for more details. If None, will not apply any mask (default None). name : str, optional - The name of the marker. By default, it will use - KIND_FunctionalConnectivitySpheres where KIND is the kind of data it - was applied to (default None). + The name of the marker. If None, will use + ``BOLD_FunctionalConnectivitySpheres`` (default None). """ @@ -62,8 +69,8 @@ class FunctionalConnectivitySpheres(FunctionalConnectivityBase): allow_overlap: bool = False, agg_method: str = "mean", agg_method_params: Optional[Dict] = None, - cor_method: str = "covariance", - cor_method_params: Optional[Dict] = None, + conn_method: str = "correlation", + conn_method_params: Optional[Dict] = None, masks: Union[str, Dict, List[Union[Dict, str]], None] = None, name: Optional[str] = None, ) -> None: @@ -75,8 +82,8 @@ class FunctionalConnectivitySpheres(FunctionalConnectivityBase): super().__init__( agg_method=agg_method, agg_method_params=agg_method_params, - cor_method=cor_method, - cor_method_params=cor_method_params, + conn_method=conn_method, + conn_method_params=conn_method_params, masks=masks, name=name, ) diff --git a/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py b/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py index 2ee35d614..67a046286 100644 --- a/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py +++ b/junifer/markers/functional_connectivity/tests/test_crossparcellation_functional_connectivity.py @@ -27,7 +27,7 @@ def test_init() -> None: CrossParcellationFC( parcellation_one="a", parcellation_two="a", - correlation_method="pearson", + corr_method="pearson", ) @@ -58,7 +58,7 @@ def test_compute(tmp_path: Path) -> None: crossparcellation = CrossParcellationFC( parcellation_one=parcellation_one, parcellation_two=parcellation_two, - correlation_method="spearman", + corr_method="spearman", ) out = crossparcellation.compute(element_data["BOLD"])[ "functional_connectivity" @@ -86,7 +86,7 @@ def test_store(tmp_path: Path) -> None: crossparcellation = CrossParcellationFC( parcellation_one=parcellation_one, parcellation_two=parcellation_two, - correlation_method="spearman", + corr_method="spearman", ) storage = SQLiteFeatureStorage( uri=tmp_path / "test_crossparcellation.sqlite", upsert="ignore" diff --git a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py index 1932718ff..580ea3e08 100644 --- a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py +++ b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_parcels.py @@ -5,6 +5,9 @@ # License: AGPL from pathlib import Path +from typing import Dict + +import pytest from junifer.datareader import DefaultDataReader from junifer.markers.functional_connectivity import EdgeCentricFCParcels @@ -12,20 +15,35 @@ from junifer.storage import SQLiteFeatureStorage from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber -def test_EdgeCentricFCParcels(tmp_path: Path) -> None: +@pytest.mark.parametrize( + "conn_method_params", + [ + {"empirical": False}, + {"empirical": True}, + ], +) +def test_EdgeCentricFCParcels( + tmp_path: Path, + conn_method_params: Dict[str, bool], +) -> None: """Test EdgeCentricFCParcels. Parameters ---------- tmp_path : pathlib.Path The path to the test directory. + conn_method_params : dict + The parametrized parameters to connectivity measure method. """ with PartlyCloudyTestingDataGrabber() as dg: + # Get element data element_data = DefaultDataReader().fit_transform(dg["sub-01"]) + # Setup marker marker = EdgeCentricFCParcels( parcellation="TianxS1x3TxMNInonlinear2009cAsym", - cor_method_params={"empirical": True}, + conn_method="correlation", + conn_method_params=conn_method_params, ) # Check correct output assert "matrix" == marker.get_output_type( @@ -41,8 +59,7 @@ def test_EdgeCentricFCParcels(tmp_path: Path) -> None: assert "data" in edge_fc_bold assert "row_names" in edge_fc_bold assert "col_names" in edge_fc_bold - assert edge_fc_bold["data"].shape[0] == n_edges - assert edge_fc_bold["data"].shape[1] == n_edges + assert edge_fc_bold["data"].shape == (n_edges, n_edges) assert len(set(edge_fc_bold["row_names"])) == n_edges assert len(set(edge_fc_bold["col_names"])) == n_edges diff --git a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py index 112891cec..4ef77a382 100644 --- a/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py +++ b/junifer/markers/functional_connectivity/tests/test_edge_functional_connectivity_spheres.py @@ -5,6 +5,9 @@ # License: AGPL from pathlib import Path +from typing import Dict + +import pytest from junifer.datareader import DefaultDataReader from junifer.markers.functional_connectivity import EdgeCentricFCSpheres @@ -12,19 +15,36 @@ from junifer.storage import SQLiteFeatureStorage from junifer.testing.datagrabbers import SPMAuditoryTestingDataGrabber -def test_EdgeCentricFCSpheres(tmp_path: Path) -> None: +@pytest.mark.parametrize( + "conn_method_params", + [ + {"empirical": False}, + {"empirical": True}, + ], +) +def test_EdgeCentricFCSpheres( + tmp_path: Path, + conn_method_params: Dict[str, bool], +) -> None: """Test EdgeCentricFCSpheres. Parameters ---------- tmp_path : pathlib.Path The path to the test directory. + conn_method_params : dict + The parametrized parameters to connectivity measure method. """ with SPMAuditoryTestingDataGrabber() as dg: + # Get element data element_data = DefaultDataReader().fit_transform(dg["sub001"]) + # Setup marker marker = EdgeCentricFCSpheres( - coords="DMNBuckner", radius=5.0, cor_method="correlation" + coords="DMNBuckner", + radius=5.0, + conn_method="correlation", + conn_method_params=conn_method_params, ) # Check correct output assert "matrix" == marker.get_output_type( @@ -45,13 +65,6 @@ def test_EdgeCentricFCSpheres(tmp_path: Path) -> None: assert len(set(edge_fc_bold["row_names"])) == n_edges assert len(set(edge_fc_bold["col_names"])) == n_edges - # Check empirical correlation method parameters - marker = EdgeCentricFCSpheres( - coords="DMNBuckner", - radius=5.0, - cor_method="correlation", - cor_method_params={"empirical": True}, - ) # Store storage = SQLiteFeatureStorage( uri=tmp_path / "test_edge_fc_spheres.sqlite", upsert="ignore" diff --git a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py index 79847cf18..4ec804b8b 100644 --- a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py +++ b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_parcels.py @@ -6,10 +6,13 @@ # License: AGPL from pathlib import Path +from typing import TYPE_CHECKING, Dict, Type +import pytest from nilearn.connectome import ConnectivityMeasure from nilearn.maskers import NiftiLabelsMasker from numpy.testing import assert_array_almost_equal +from sklearn.covariance import EmpiricalCovariance, LedoitWolf from junifer.data import get_parcellation from junifer.datareader import DefaultDataReader @@ -20,19 +23,42 @@ from junifer.storage import SQLiteFeatureStorage from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber -def test_FunctionalConnectivityParcels(tmp_path: Path) -> None: +if TYPE_CHECKING: + from sklearn.base import BaseEstimator + + +@pytest.mark.parametrize( + "conn_method_params, cov_estimator", + [ + ({"empirical": False}, LedoitWolf(store_precision=False)), + ({"empirical": True}, EmpiricalCovariance(store_precision=False)), + ], +) +def test_FunctionalConnectivityParcels( + tmp_path: Path, + conn_method_params: Dict[str, bool], + cov_estimator: Type["BaseEstimator"], +) -> None: """Test FunctionalConnectivityParcels. Parameters ---------- tmp_path : pathlib.Path The path to the test directory. + conn_method_params : dict + The parametrized parameters to connectivity measure method. + cov_estimator : estimator object + The parametrized covariance estimator. """ with PartlyCloudyTestingDataGrabber() as dg: + # Get element data element_data = DefaultDataReader().fit_transform(dg["sub-01"]) + # Setup marker marker = FunctionalConnectivityParcels( - parcellation="TianxS1x3TxMNInonlinear2009cAsym" + parcellation="TianxS1x3TxMNInonlinear2009cAsym", + conn_method="correlation", + conn_method_params=conn_method_params, ) # Check correct output assert "matrix" == marker.get_output_type( @@ -65,7 +91,7 @@ def test_FunctionalConnectivityParcels(tmp_path: Path) -> None: ) # Compute the connectivity measure connectivity_measure = ConnectivityMeasure( - kind="covariance" + cov_estimator=cov_estimator, kind="correlation" # type: ignore ).fit_transform([extracted_timeseries])[0] # Check that FC are almost equal @@ -73,11 +99,6 @@ def test_FunctionalConnectivityParcels(tmp_path: Path) -> None: connectivity_measure, fc_bold["data"], decimal=3 ) - # Check empirical correlation method parameters - marker = FunctionalConnectivityParcels( - parcellation="TianxS1x3TxMNInonlinear2009cAsym", - cor_method_params={"empirical": True}, - ) # Store storage = SQLiteFeatureStorage( uri=tmp_path / "test_fc_parcels.sqlite", upsert="ignore" diff --git a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py index 40ce9d457..ddf7af192 100644 --- a/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py +++ b/junifer/markers/functional_connectivity/tests/test_functional_connectivity_spheres.py @@ -7,12 +7,13 @@ # License: AGPL from pathlib import Path +from typing import TYPE_CHECKING, Dict, Type import pytest from nilearn.connectome import ConnectivityMeasure from nilearn.maskers import NiftiSpheresMasker from numpy.testing import assert_array_almost_equal -from sklearn.covariance import EmpiricalCovariance +from sklearn.covariance import EmpiricalCovariance, LedoitWolf from junifer.data import get_coordinates from junifer.datareader import DefaultDataReader @@ -23,19 +24,43 @@ from junifer.storage import SQLiteFeatureStorage from junifer.testing.datagrabbers import SPMAuditoryTestingDataGrabber -def test_FunctionalConnectivitySpheres(tmp_path: Path) -> None: +if TYPE_CHECKING: + from sklearn.base import BaseEstimator + + +@pytest.mark.parametrize( + "conn_method_params, cov_estimator", + [ + ({"empirical": False}, LedoitWolf(store_precision=False)), + ({"empirical": True}, EmpiricalCovariance(store_precision=False)), + ], +) +def test_FunctionalConnectivitySpheres( + tmp_path: Path, + conn_method_params: Dict[str, bool], + cov_estimator: Type["BaseEstimator"], +) -> None: """Test FunctionalConnectivitySpheres. Parameters ---------- tmp_path : pathlib.Path The path to the test directory. + conn_method_params : dict + The parametrized parameters to connectivity measure method. + cov_estimator : estimator object + The parametrized covariance estimator. """ with SPMAuditoryTestingDataGrabber() as dg: + # Get element data element_data = DefaultDataReader().fit_transform(dg["sub001"]) + # Setup marker marker = FunctionalConnectivitySpheres( - coords="DMNBuckner", radius=5.0, cor_method="correlation" + coords="DMNBuckner", + radius=5.0, + conn_method="correlation", + conn_method_params=conn_method_params, ) # Check correct output assert "matrix" == marker.get_output_type( @@ -67,7 +92,7 @@ def test_FunctionalConnectivitySpheres(tmp_path: Path) -> None: ) # Compute the connectivity measure connectivity_measure = ConnectivityMeasure( - kind="correlation" + cov_estimator=cov_estimator, kind="correlation" # type: ignore ).fit_transform([extracted_timeseries])[0] # Check that FC are almost equal @@ -88,65 +113,9 @@ def test_FunctionalConnectivitySpheres(tmp_path: Path) -> None: ) -def test_FunctionalConnectivitySpheres_empirical(tmp_path: Path) -> None: - """Test FunctionalConnectivitySpheres with empirical covariance. - - Parameters - ---------- - tmp_path : pathlib.Path - The path to the test directory. - - """ - with SPMAuditoryTestingDataGrabber() as dg: - element_data = DefaultDataReader().fit_transform(dg["sub001"]) - marker = FunctionalConnectivitySpheres( - coords="DMNBuckner", - radius=5.0, - cor_method="correlation", - cor_method_params={"empirical": True}, - ) - # Check correct output - assert "matrix" == marker.get_output_type( - input_type="BOLD", output_feature="functional_connectivity" - ) - - # Fit-transform the data - fc = marker.fit_transform(element_data) - fc_bold = fc["BOLD"]["functional_connectivity"] - - assert "data" in fc_bold - assert "row_names" in fc_bold - assert "col_names" in fc_bold - assert fc_bold["data"].shape == (6, 6) - assert len(set(fc_bold["row_names"])) == 6 - assert len(set(fc_bold["col_names"])) == 6 - - # Compare with nilearn - # Load testing coordinates for the target data - testing_coords, _ = get_coordinates( - coords="DMNBuckner", target_data=element_data["BOLD"] - ) - # Extract timeseries - nifti_spheres_masker = NiftiSpheresMasker( - seeds=testing_coords, radius=5.0 - ) - extracted_timeseries = nifti_spheres_masker.fit_transform( - element_data["BOLD"]["data"] - ) - # Compute the connectivity measure - connectivity_measure = ConnectivityMeasure( - cov_estimator=EmpiricalCovariance(), kind="correlation" # type: ignore - ).fit_transform([extracted_timeseries])[0] - - # Check that FC are almost equal - assert_array_almost_equal( - connectivity_measure, fc_bold["data"], decimal=3 - ) - - def test_FunctionalConnectivitySpheres_error() -> None: """Test FunctionalConnectivitySpheres errors.""" with pytest.raises(ValueError, match="radius should be > 0"): FunctionalConnectivitySpheres( - coords="DMNBuckner", radius=-0.1, cor_method="correlation" + coords="DMNBuckner", radius=-0.1, conn_method="correlation" )