diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8f4b829ab..6bb84762f 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -22,8 +22,8 @@ jobs: - name: Set up system run: | bash -c "$(curl -fsSL http://neuro.debian.net/_files/neurodebian-travis.sh)" - sudo apt-get update -qq - sudo apt-get install git-annex-standalone + sudo apt-get -qq update + sudo apt-get -qq install git-annex-standalone - name: Configure git for datalad run: | git config --global user.email "runner@github.com" @@ -39,8 +39,29 @@ jobs: python -m pip install tox tox-gh-actions - name: Install AFNI run: | - docker pull afni/afni_make_build - echo "$(pwd)/junifer/api/res/afni" >> $GITHUB_PATH + echo "++ Distro information" + lsb_release -a + + echo "++ Add universe repo" + sudo add-apt-repository -y universe + + echo "++ Add PPAs for R programs" + sudo add-apt-repository -y "ppa:marutter/rrutter4.0" + sudo add-apt-repository -y "ppa:c2d4u.team/c2d4u4.0+" + + echo "++ Update package manager info" + sudo apt-get -qq update + + echo "++ Get main dependencies" + sudo apt-get -qq install -y tcsh libssl-dev gsl-bin libjpeg62 vim curl \ + build-essential libcurl4-openssl-dev libxml2-dev libgfortran-11-dev \ + libgomp1 r-base cmake rsync libxm4 + + sudo ln -s /usr/lib/x86_64-linux-gnu/libgsl.so.27 /usr/lib/x86_64-linux-gnu/libgsl.so.19 + + curl -O https://afni.nimh.nih.gov/pub/dist/bin/misc/@update.afni.binaries + sudo tcsh @update.afni.binaries -package linux_ubuntu_16_64 -bindir /afni + echo "/afni" >> $GITHUB_PATH if: matrix.python-version == 3.10 - name: Check AFNI run: | diff --git a/docs/builtin.rst b/docs/builtin.rst index 84bb10ce0..df91e9924 100644 --- a/docs/builtin.rst +++ b/docs/builtin.rst @@ -166,6 +166,14 @@ Available - Compute root sum of squares of edgewise timeseries - Done - 0.0.1 + * - :class:`junifer.markers.ReHoParcels` + - Calculate regional homogeneity over parcellation + - Done + - 0.0.1 + * - :class:`junifer.markers.ReHoSpheres` + - Calculate regional homogeneity over spheres placed on coordinates + - Done + - 0.0.1 Planned @@ -184,9 +192,6 @@ Planned * - ALFF and (f)ALFF - Detect amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI - :gh:`35` - * - ReHo - - Calculate regional homogeneity - - :gh:`36` * - Permutation entropy, Range entropy, Multiscale entropy and Hurst exponent - Calculate Permutation entropy, Range entropy, Multiscale entropy and Hurst exponent - :gh:`61` diff --git a/docs/changes/latest.inc b/docs/changes/latest.inc index d68363a6f..196986ab5 100644 --- a/docs/changes/latest.inc +++ b/docs/changes/latest.inc @@ -96,6 +96,8 @@ Enhancements - Refactor :class:`junifer.pipeline.PipelineStepMixin` to improve its implementation and validation for pipeline steps (:gh:`152` by `Synchon Mandal`_). +- Implement :class:`junifer.markers.ReHoParcels` and :class:`junifer.markers.ReHoSpheres` markers (:gh:`36` by `Synchon Mandal`_). + Bugs ~~~~ diff --git a/junifer/api/cli.py b/junifer/api/cli.py index ae0c13dfc..55f39f77a 100644 --- a/junifer/api/cli.py +++ b/junifer/api/cli.py @@ -308,6 +308,7 @@ def setup() -> None: # pragma: no cover def afni_docker() -> None: # pragma: no cover """Configure AFNI-Docker wrappers.""" import junifer + pkg_path = Path(junifer.__path__[0]) # type: ignore afni_wrappers_path = pkg_path / "api" / "res" / "afni" msg = f""" diff --git a/junifer/api/res/afni/3dAFNItoNIFTI b/junifer/api/res/afni/3dAFNItoNIFTI new file mode 100755 index 000000000..37e77f583 --- /dev/null +++ b/junifer/api/res/afni/3dAFNItoNIFTI @@ -0,0 +1,3 @@ +#!/bin/bash + +run_afni_docker.sh 3dAFNItoNIFTI "$@" \ No newline at end of file diff --git a/junifer/markers/__init__.py b/junifer/markers/__init__.py index c66f97ae2..2c853aad0 100644 --- a/junifer/markers/__init__.py +++ b/junifer/markers/__init__.py @@ -12,3 +12,4 @@ from .functional_connectivity_parcels import FunctionalConnectivityParcels from .functional_connectivity_spheres import FunctionalConnectivitySpheres from .parcel_aggregation import ParcelAggregation from .sphere_aggregation import SphereAggregation +from .reho import ReHoParcels, ReHoSpheres diff --git a/junifer/markers/reho/__init__.py b/junifer/markers/reho/__init__.py new file mode 100644 index 000000000..e2fb886d0 --- /dev/null +++ b/junifer/markers/reho/__init__.py @@ -0,0 +1,7 @@ +"""Provide imports for reho sub-package.""" + +# Authors: Synchon Mandal +# License: AGPL + +from .reho_parcels import ReHoParcels +from .reho_spheres import ReHoSpheres diff --git a/junifer/markers/reho/reho_base.py b/junifer/markers/reho/reho_base.py new file mode 100644 index 000000000..6aaadbf99 --- /dev/null +++ b/junifer/markers/reho/reho_base.py @@ -0,0 +1,126 @@ +"""Provide base class for regional homogeneity (ReHo).""" + +# Authors: Synchon Mandal +# License: AGPL + + +from typing import TYPE_CHECKING, Any, Dict, List, Optional + +from ...utils import logger, raise_error +from ..base import BaseMarker +from .reho_estimator import ReHoEstimator + + +if TYPE_CHECKING: + from nibabel import Nifti1Image + + +class ReHoBase(BaseMarker): + """Base class for regional homogeneity computation. + + Parameters + ---------- + use_afni : bool, optional + Whether to use AFNI for computing. If None, will use AFNI only + if available (default None). + name : str, optional + The name of the marker. If None, it will use the class name + (default None). + + """ + + _EXT_DEPENDENCIES = [ + { + "name": "afni", + "optional": True, + "commands": ["3dReHo", "3dAFNItoNIFTI"], + }, + ] + + def __init__( + self, + use_afni: Optional[bool] = None, + name: Optional[str] = None, + ) -> None: + super().__init__(on="BOLD", name=name) + self.use_afni = use_afni + + def get_valid_inputs(self) -> List[str]: + """Get valid data types for input. + + Returns + ------- + list of str + The list of data types that can be used as input for this marker. + + """ + return ["BOLD"] + + def get_output_type(self, input_type: str) -> str: + """Get output type. + + Parameters + ---------- + input_type : str + The data type input to the marker. + + Returns + ------- + str + The storage type output by the marker. + + """ + return "table" + + def compute_reho_map( + self, + input: Dict[str, Any], + **reho_params: Any, + ) -> "Nifti1Image": + """Compute. + + Calculates Kendall's W per voxel using neighborhood voxels. + Instead of the time series values themselves, Kendall's W uses the + relative rank ordering of a hood over all time points to evaluate + a parameter W in range 0-1, with 0 reflecting no trend of agreement + between time series and 1 reflecting perfect agreement. For more + information about the method, please check [1]_. + + Parameters + ---------- + input : dict + The BOLD data as dictionary. + **reho_params : dict + Extra keyword arguments for ReHo. + + Returns + ------- + Niimg-like object + + References + ---------- + .. [1] Jiang, L., & Zuo, X. N. (2016). + Regional Homogeneity: A Multimodal, Multiscale Neuroimaging + Marker of the Human Connectome. + The Neuroscientist, Volume 22(5), Pages 486–505. + https://doi.org/10.1177/1073858415595004 + + """ + if self.use_afni is None: + raise_error( + "Parameter `use_afni` must be set to True or False in order " + "to compute this marker. It is currently set to None (default " + "behaviour). This is intended to be for auto-detection. In " + "order for that to happen, please call the `validate` method " + "before calling the `compute` method." + ) + logger.info("Calculating ReHO map.") + # Initialize reho estimator + reho_estimator = ReHoEstimator() + # Fit-transform reho estimator + reho_map = reho_estimator.fit_transform( + use_afni=self.use_afni, + input_data=input, + **reho_params, + ) + return reho_map diff --git a/junifer/markers/reho/reho_estimator.py b/junifer/markers/reho/reho_estimator.py new file mode 100644 index 000000000..69cac2ecd --- /dev/null +++ b/junifer/markers/reho/reho_estimator.py @@ -0,0 +1,510 @@ +"""Provide estimator class for regional homogeneity (ReHo).""" + +# Authors: Synchon Mandal +# Federico Raimondo +# License: AGPL + + +import hashlib +import shutil +import subprocess +import tempfile +from functools import lru_cache +from itertools import product +from pathlib import Path +from typing import TYPE_CHECKING, Any, Dict, List, Optional, cast + +import nibabel as nib +import numpy as np +from nilearn import image as nimg +from nilearn import masking as nmask +from scipy.stats import rankdata + +from ...utils import logger, raise_error +from ..utils import singleton + + +if TYPE_CHECKING: + from nibabel import Nifti1Image + + +@singleton +class ReHoEstimator: + """Estimator class for regional homogeneity. + + This class is a singleton and is used for efficient computation of ReHo, + by caching the ReHo map for a given set of file path and computation + parameters. + + .. warning:: This class can only be used via ReHoBase() and is a deliberate + decision as it serves a specific purpose. + + Attributes + ---------- + temp_dir_path : pathlib.Path + Path to the temporary directory for assets storage. + + """ + + def __init__(self) -> None: + self._file_path = None + # Create temporary directory for intermittent storage of assets during + # computation via afni's 3dReHo + self.temp_dir_path = Path(tempfile.mkdtemp()) + + def __del__(self) -> None: + """Cleanup.""" + # Delete temporary directory and ignore errors for read-only files + shutil.rmtree(self.temp_dir_path, ignore_errors=True) + + def _compute_reho_afni( + self, + data: "Nifti1Image", + nneigh: int = 27, + neigh_rad: Optional[float] = None, + neigh_x: Optional[float] = None, + neigh_y: Optional[float] = None, + neigh_z: Optional[float] = None, + box_rad: Optional[int] = None, + box_x: Optional[int] = None, + box_y: Optional[int] = None, + box_z: Optional[int] = None, + ) -> "Nifti1Image": + """Compute ReHo map via afni's 3dReHo. + + Parameters + ---------- + data : 4D Niimg-like object + Images to process. + nneigh : {7, 19, 27}, optional + Number of voxels in the neighbourhood, inclusive. Can be: + + * 7 : for facewise neighbours only + * 19 : for face- and edge-wise nieghbours + * 27 : for face-, edge-, and node-wise neighbors + + (default 27). + neigh_rad : positive float, optional + The radius of a desired neighbourhood (default None). + neigh_x : positive float, optional + The semi-radius for x-axis of ellipsoidal volumes (default None). + neigh_y : positive float, optional + The semi-radius for y-axis of ellipsoidal volumes (default None). + neigh_z : positive float, optional + The semi-radius for z-axis of ellipsoidal volumes (default None). + box_rad : positive int, optional + The number of voxels outward in a given cardinal direction for a + cubic box centered on a given voxel (default None). + box_x : positive int, optional + The number of voxels for +/- x-axis of cuboidal volumes + (default None). + box_y : positive int, optional + The number of voxels for +/- y-axis of cuboidal volumes + (default None). + box_z : positive int, optional + The number of voxels for +/- z-axis of cuboidal volumes + (default None). + + Returns + ------- + Niimg-like object + + Raises + ------ + RuntimeError + If the 3dReHo command fails due to some issue. + + Notes + ----- + For more information on the publication, please check [1]_ , and for + 3dReHo help check: + https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dReHo.html + + Please note that that you cannot mix ``box_*`` and ``neigh_*`` + arguments. The arguments are prioritized by their order in the function + signature. + + As the process also depends on the conversion of AFNI files to NIFTI + via afni's 3dAFNItoNIFTI, the help for that can be found at: + https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dAFNItoNIFTI.html + + References + ---------- + .. [1] Taylor, P.A., & Saad, Z.S. (2013). + FATCAT: (An Efficient) Functional And Tractographic Connectivity + Analysis Toolbox. + Brain connectivity, Volume 3(5), Pages 523-35. + https://doi.org/10.1089/brain.2013.0154 + + """ + # Save niimg to nii.gz + nifti_in_file_path = self.temp_dir_path / "input.nii" + nib.save(data, nifti_in_file_path) + + # Set 3dReHo command + reho_afni_out_path_prefix = self.temp_dir_path / "reho" + reho_cmd: List[str] = [ + "3dReHo", + f"-prefix {reho_afni_out_path_prefix.resolve()}", + f"-inset {nifti_in_file_path.resolve()}", + ] + # Check ellipsoidal / cuboidal volume arguments + if neigh_rad: + reho_cmd.append(f"-neigh_RAD {neigh_rad}") + elif neigh_x and neigh_y and neigh_z: + reho_cmd.extend( + [ + f"-neigh_X {neigh_x}", + f"-neigh_Y {neigh_y}", + f"-neigh_Z {neigh_z}", + ] + ) + elif box_rad: + reho_cmd.append(f"-box_RAD {box_rad}") + elif box_x and box_y and box_z: + reho_cmd.extend( + [f"-box_X {box_x}", f"-box_Y {box_y}", f"-box_Z {box_z}"] + ) + else: + reho_cmd.append(f"-nneigh {nneigh}") + # Call 3dReHo + reho_cmd_str = " ".join(reho_cmd) + logger.info(f"3dReHo command to be executed: {reho_cmd_str}") + reho_process = subprocess.run( + reho_cmd_str, # string needed with shell=True + stdin=subprocess.DEVNULL, + shell=True, # needed for respecting $PATH + check=False, + ) + if reho_process.returncode == 0: + logger.info( + "3dReHo succeeded with the following output: " + f"{reho_process.stdout}" + ) + else: + raise_error( + msg="3dReHo failed with the following error: " + f"{reho_process.stdout}", + klass=RuntimeError, + ) + + # SHA256 for bypassing memmap + sha256_params = hashlib.sha256(bytes(reho_cmd_str, "utf-8")) + # Convert afni to nifti + reho_afni_to_nifti_out_path = ( + self.temp_dir_path / f"output_{sha256_params.hexdigest()}.nii" + ) + convert_cmd: List[str] = [ + "3dAFNItoNIFTI", + f"-prefix {reho_afni_to_nifti_out_path.resolve()}", + f"{reho_afni_out_path_prefix}+tlrc.BRIK", + ] + # Call 3dAFNItoNIFTI + convert_cmd_str = " ".join(convert_cmd) + logger.info(f"3dAFNItoNIFTI command to be executed: {convert_cmd_str}") + convert_process = subprocess.run( + convert_cmd_str, # string needed with shell=True + stdin=subprocess.DEVNULL, + shell=True, # needed for respecting $PATH + check=False, + ) + if convert_process.returncode == 0: + logger.info( + "3dAFNItoNIFTI succeeded with the following output: " + f"{convert_process.stdout}" + ) + else: + raise_error( + msg="3dAFNItoNIFTI failed with the following error: " + f"{convert_process.stdout}", + klass=RuntimeError, + ) + + # Cleanup intermediate files + for fname in self.temp_dir_path.glob("reho*"): + fname.unlink() + + # Load nifti + output_data = nib.load(reho_afni_to_nifti_out_path) + # Stupid casting + output_data = cast("Nifti1Image", output_data) + return output_data + + def _compute_reho_python( + self, + data: "Nifti1Image", + nneigh: int = 27, + ) -> "Nifti1Image": + """Compute ReHo map. + + Parameters + ---------- + data : 4D Niimg-like object + Images to process. + nneigh : {7, 19, 27, 125}, optional + Number of voxels in the neighbourhood, inclusive. Can be: + + * 7 : for facewise neighbours only + * 19 : for face- and edge-wise nieghbours + * 27 : for face-, edge-, and node-wise neighbors + * 125 : for 5x5 cuboidal volume + + (default 27). + Returns + ------- + Niimg-like object + + Raises + ------ + ValueError + If ``nneigh`` is invalid. + + """ + valid_nneigh = (7, 19, 27, 125) + if nneigh not in valid_nneigh: + raise_error( + f"Invalid value for `nneigh`, should be one of {valid_nneigh}." + ) + + logger.info(f"Computing ReHo map using {nneigh} neighbours.") + # Get scan data + niimg_data = data.get_fdata() + # Get scan dimensions + n_x, n_y, n_z, _ = niimg_data.shape + + # Get rank of every voxel across time series + ranks_niimg_data = rankdata(niimg_data, axis=-1) + + # Initialize 3D array to store tied rank correction for every voxel + tied_rank_corrections = np.zeros((n_x, n_y, n_z), dtype=np.float64) + # Calculate tied rank correction for every voxel + for i_x, i_y, i_z in product(range(n_x), range(n_y), range(n_z)): + # Calculate tied rank count for every voxel across time series + _, tie_count = np.unique( + ranks_niimg_data[i_x, i_y, i_z, :], + return_counts=True, + ) + # Calculate and store tied rank correction for every voxel across + # timeseries + tied_rank_corrections[i_x, i_y, i_z] = np.sum( + tie_count**3 - tie_count + ) + + # Initialize 3D array to store reho map + reho_map = np.ones((n_x, n_y, n_z), dtype=np.float32) + + # Calculate whole brain mask + mni152_whole_brain_mask = nmask.compute_brain_mask( + data, threshold=0.5, mask_type="whole-brain" + ) + # Convert 0 / 1 array to bool + logical_mni152_whole_brain_mask = ( + mni152_whole_brain_mask.get_fdata().astype(bool) + ) + + # Create mask cluster and set start and end indices + if nneigh in (7, 19, 27): + mask_cluster = np.ones((3, 3, 3)) + + if nneigh == 7: + mask_cluster[0, 0, 0] = 0 + mask_cluster[0, 1, 0] = 0 + mask_cluster[0, 2, 0] = 0 + mask_cluster[0, 0, 1] = 0 + mask_cluster[0, 2, 1] = 0 + mask_cluster[0, 0, 2] = 0 + mask_cluster[0, 1, 2] = 0 + mask_cluster[0, 2, 2] = 0 + mask_cluster[1, 0, 0] = 0 + mask_cluster[1, 2, 0] = 0 + mask_cluster[1, 0, 2] = 0 + mask_cluster[1, 2, 2] = 0 + mask_cluster[2, 0, 0] = 0 + mask_cluster[2, 1, 0] = 0 + mask_cluster[2, 2, 0] = 0 + mask_cluster[2, 0, 1] = 0 + mask_cluster[2, 2, 1] = 0 + mask_cluster[2, 0, 2] = 0 + mask_cluster[2, 1, 2] = 0 + mask_cluster[2, 2, 2] = 0 + + elif nneigh == 19: + mask_cluster[0, 0, 0] = 0 + mask_cluster[0, 2, 0] = 0 + mask_cluster[2, 0, 0] = 0 + mask_cluster[2, 2, 0] = 0 + mask_cluster[0, 0, 2] = 0 + mask_cluster[0, 2, 2] = 0 + mask_cluster[2, 0, 2] = 0 + mask_cluster[2, 2, 2] = 0 + + start_idx = 1 + end_idx = 2 + + elif nneigh == 125: + mask_cluster = np.ones((5, 5, 5)) + start_idx = 2 + end_idx = 3 + + # Convert 0 / 1 array to bool + logical_mask_cluster = mask_cluster.astype(bool) + + for i, j, k in product( + range(start_idx, n_x - (end_idx - 1)), + range(start_idx, n_y - (end_idx - 1)), + range(start_idx, n_z - (end_idx - 1)), + ): + # Get mask only for neighbourhood + logical_neighbourhood_mni152_whole_brain_mask = ( + logical_mni152_whole_brain_mask[ + i - start_idx : i + end_idx, + j - start_idx : j + end_idx, + k - start_idx : k + end_idx, + ] + ) + # Perform logical AND to get neighbourhood mask; + # done to take care of brain boundaries + neighbourhood_mask = ( + logical_mask_cluster + & logical_neighbourhood_mni152_whole_brain_mask + ) + # Continue if voxel is restricted by mask + if neighbourhood_mask[1, 1, 1] == 0: + continue + + # Get ranks for the neighbourhood + neighbourhood_ranks = ranks_niimg_data[ + i - start_idx : i + end_idx, + j - start_idx : j + end_idx, + k - start_idx : k + end_idx, + :, + ] + # Get tied ranks corrections for the neighbourhood + neighbourhood_tied_ranks_corrections = tied_rank_corrections[ + i - start_idx : i + end_idx, + j - start_idx : j + end_idx, + k - start_idx : k + end_idx, + ] + # Mask neighbourhood ranks + masked_neighbourhood_ranks = neighbourhood_ranks[ + logical_mask_cluster, : + ] + # Mask tied ranks corrections for the neighbourhood + masked_tied_rank_corrections = ( + neighbourhood_tied_ranks_corrections[logical_mask_cluster] + ) + # Calculate KCC + reho_map[i, j, k] = _kendall_w_reho( + timeseries_ranks=masked_neighbourhood_ranks, + tied_rank_corrections=masked_tied_rank_corrections, + ) + + output = nimg.new_img_like(data, reho_map, copy_header=False) + return output + + @lru_cache(maxsize=None, typed=True) + def _compute( + self, + use_afni: bool, + data: "Nifti1Image", + **reho_params: Any, + ) -> "Nifti1Image": + """Compute the ReHo map with memoization. + + Parameters + ---------- + use_afni : bool + Whether to use afni or not. + data : 4D Niimg-like object + Images to process. + **reho_params : dict + Extra keyword arguments for ReHo. + + Returns + ------- + Niimg-like object + + """ + if use_afni: + output = self._compute_reho_afni(data, **reho_params) + else: + output = self._compute_reho_python(data, **reho_params) + return output + + def fit_transform( + self, + use_afni: bool, + input_data: Dict[str, Any], + **reho_params: Any, + ) -> "Nifti1Image": + """Fit and transform for the estimator. + + Parameters + ---------- + use_afni : bool + Whether to use afni or not. + input_data : dict + The BOLD data as dictionary. + **reho_params : dict + Extra keyword arguments for ReHo. + + Returns + ------- + Niimg-like object + + """ + bold_path = input_data["path"] + bold_data = input_data["data"] + # Clear cache if file path is different from when caching was done + if self._file_path != bold_path: + logger.info(f"Removing ReHo map cache at {self._file_path}.") + # Clear the cache + self._compute.cache_clear() + # Clear temporary directory files + for file_ in self.temp_dir_path.iterdir(): + file_.unlink(missing_ok=True) + # Set the new file path + self._file_path = bold_path + else: + logger.info(f"Using ReHo map cache at {self._file_path}.") + # Compute + return self._compute(use_afni, bold_data, **reho_params) + + +def _kendall_w_reho( + timeseries_ranks: np.ndarray, tied_rank_corrections: np.ndarray +) -> float: + """Calculate Kendall's coefficient of concordance (KCC) for ReHo map. + + ..note:: This function should only be used to calculate KCC for a ReHo map. + For general use, check out ``junifer.stats.kendall_w``. + + Parameters + ---------- + timeseries_matrix : 2D numpy.ndarray + A matrix of ranks of a subset subject's brain voxels. + tied_rank_corrections : 3D numpy.ndarray + A 3D array consisting of the tied rank corrections for the ranks + of a subset subject's brain voxels. + + Returns + ------- + float + Kendall's W (KCC) of the given timeseries matrix. + + """ + m, n = timeseries_ranks.shape # annotators X items + + numerator = (12 * np.sum(np.square(np.sum(timeseries_ranks, axis=0)))) - ( + 3 * m**2 * n * (n + 1) ** 2 + ) + denominator = (m**2 * n * (n**2 - 1)) - ( + m * np.sum(tied_rank_corrections) + ) + + if denominator == 0: + kcc = 1.0 + else: + kcc = numerator / denominator + + return kcc diff --git a/junifer/markers/reho/reho_parcels.py b/junifer/markers/reho/reho_parcels.py new file mode 100644 index 000000000..451ee01a5 --- /dev/null +++ b/junifer/markers/reho/reho_parcels.py @@ -0,0 +1,148 @@ +"""Provide class for regional homogeneity (ReHo) on parcels.""" + +# Authors: Synchon Mandal +# License: AGPL + + +from typing import Any, Dict, Optional + +import numpy as np + +from ...api.decorators import register_marker +from ...utils import logger +from ..parcel_aggregation import ParcelAggregation +from .reho_base import ReHoBase + + +@register_marker +class ReHoParcels(ReHoBase): + """Class for regional homogeneity on parcels. + + Parameters + ---------- + parcellation : str + The name of the parcellation. Check valid options by calling + :func:`junifer.data.parcellations.list_parcellations`. + use_afni : bool, optional + Whether to use AFNI for computing. If None, will use AFNI only + if available (default None). + reho_params : dict, optional + Extra parameters for computing ReHo map as a dictionary (default None). + If ``use_afni = True``, then the valid keys are: + + * ``nneigh`` : {7, 19, 27}, optional (default 27) + Number of voxels in the neighbourhood, inclusive. Can be: + + - 7 : for facewise neighbours only + - 19 : for face- and edge-wise nieghbours + - 27 : for face-, edge-, and node-wise neighbors + + * ``neigh_rad`` : positive float, optional + The radius of a desired neighbourhood (default None). + * ``neigh_x`` : positive float, optional + The semi-radius for x-axis of ellipsoidal volumes (default None). + * ``neigh_y`` : positive float, optional + The semi-radius for y-axis of ellipsoidal volumes (default None). + * ``neigh_z`` : positive float, optional + The semi-radius for z-axis of ellipsoidal volumes (default None). + * ``box_rad`` : positive int, optional + The number of voxels outward in a given cardinal direction for a + cubic box centered on a given voxel (default None). + * ``box_x`` : positive int, optional + The number of voxels for +/- x-axis of cuboidal volumes + (default None). + * ``box_y`` : positive int, optional + The number of voxels for +/- y-axis of cuboidal volumes + (default None). + * ``box_z`` : positive int, optional + The number of voxels for +/- z-axis of cuboidal volumes + (default None). + + else if ``use_afni = False``, then the valid keys are: + + * ``nneigh`` : {7, 19, 27, 125}, optional (default 27) + Number of voxels in the neighbourhood, inclusive. Can be: + + * 7 : for facewise neighbours only + * 19 : for face- and edge-wise nieghbours + * 27 : for face-, edge-, and node-wise neighbors + * 125 : for 5x5 cuboidal volume + + agg_method : str, optional + The method to perform aggregation using. Check valid options in + :func:`junifer.stats.get_aggfunc_by_name` (default "mean"). + agg_method_params : dict, optional + Parameters to pass to the aggregation function. Check valid options in + :func:`junifer.stats.get_aggfunc_by_name` (default None). + mask : str, optional + The name of the mask to apply to regions before extracting signals. + Check valid options by calling :func:`junifer.data.masks.list_masks` + (default None). + name : str, optional + The name of the marker. If None, it will use the class name + (default None). + + """ + + def __init__( + self, + parcellation: str, + use_afni: Optional[bool] = None, + reho_params: Optional[Dict] = None, + agg_method: str = "mean", + agg_method_params: Optional[Dict] = None, + mask: Optional[str] = None, + name: Optional[str] = None, + ) -> None: + self.parcellation = parcellation + self.reho_params = reho_params + self.agg_method = agg_method + self.agg_method_params = agg_method_params + self.mask = mask + super().__init__(use_afni=use_afni, name=name) + + def compute( + self, + input: Dict[str, Any], + extra_input: Optional[Dict[str, Any]] = None, + ) -> Dict[str, Any]: + """Compute. + + Parameters + ---------- + input : dict + The BOLD data as dictionary. + extra_input : dict, optional + The other fields in the pipeline data object (default None). + + Returns + ------- + dict + The computed result as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a 1D numpy.ndarray + * ``columns`` : the column labels for the parcels as a list + * ``row_names`` : ``None`` + + """ + logger.info("Calculating ReHo for parcels.") + # Calculate reho map + if self.reho_params is not None: + reho_map = self.compute_reho_map(input=input, **self.reho_params) + else: + reho_map = self.compute_reho_map(input=input) + # Initialize parcel aggregation + parcel_aggregation = ParcelAggregation( + parcellation=self.parcellation, + method=self.agg_method, + method_params=self.agg_method_params, + mask=self.mask, + on="BOLD", + ) + # Perform aggregation on reho map + parcel_aggregation_input = {"data": reho_map} + output = parcel_aggregation.compute(input=parcel_aggregation_input) + # Only use the first row and expand row dimension + output["data"] = output["data"][0][np.newaxis, :] + return output diff --git a/junifer/markers/reho/reho_spheres.py b/junifer/markers/reho/reho_spheres.py new file mode 100644 index 000000000..420c5e353 --- /dev/null +++ b/junifer/markers/reho/reho_spheres.py @@ -0,0 +1,156 @@ +"""Provide class for regional homogeneity (ReHo) on spheres.""" + +# Authors: Synchon Mandal +# License: AGPL + + +from typing import Any, Dict, Optional + +import numpy as np + +from ...api.decorators import register_marker +from ...utils import logger +from ..sphere_aggregation import SphereAggregation +from .reho_base import ReHoBase + + +@register_marker +class ReHoSpheres(ReHoBase): + """Class for regional homogeneity on spheres. + + Parameters + ---------- + coords : str + The name of the coordinates list to use. See + :func:`junifer.data.coordinates.list_coordinates` for options. + radius : float, optional + The radius of the sphere in millimeters. If None, the signal will be + extracted from a single voxel. See + :class:`nilearn.maskers.NiftiSpheresMasker` for more information + (default None). + use_afni : bool, optional + Whether to use AFNI for computing. If None, will use AFNI only + if available (default None). + reho_params : dict, optional + Extra parameters for computing ReHo map as a dictionary (default None). + If ``use_afni = True``, then the valid keys are: + + * ``nneigh`` : {7, 19, 27}, optional (default 27) + Number of voxels in the neighbourhood, inclusive. Can be: + + - 7 : for facewise neighbours only + - 19 : for face- and edge-wise nieghbours + - 27 : for face-, edge-, and node-wise neighbors + + * ``neigh_rad`` : positive float, optional + The radius of a desired neighbourhood (default None). + * ``neigh_x`` : positive float, optional + The semi-radius for x-axis of ellipsoidal volumes (default None). + * ``neigh_y`` : positive float, optional + The semi-radius for y-axis of ellipsoidal volumes (default None). + * ``neigh_z`` : positive float, optional + The semi-radius for z-axis of ellipsoidal volumes (default None). + * ``box_rad`` : positive int, optional + The number of voxels outward in a given cardinal direction for a + cubic box centered on a given voxel (default None). + * ``box_x`` : positive int, optional + The number of voxels for +/- x-axis of cuboidal volumes + (default None). + * ``box_y`` : positive int, optional + The number of voxels for +/- y-axis of cuboidal volumes + (default None). + * ``box_z`` : positive int, optional + The number of voxels for +/- z-axis of cuboidal volumes + (default None). + + else if ``use_afni = False``, then the valid keys are: + + * ``nneigh`` : {7, 19, 27, 125}, optional (default 27) + Number of voxels in the neighbourhood, inclusive. Can be: + + * 7 : for facewise neighbours only + * 19 : for face- and edge-wise nieghbours + * 27 : for face-, edge-, and node-wise neighbors + * 125 : for 5x5 cuboidal volume + + agg_method : str, optional + The aggregation method to use. + See :func:`junifer.stats.get_aggfunc_by_name` for more information + (default None). + agg_method_params : dict, optional + The parameters to pass to the aggregation method (default None). + mask : str, optional + The name of the mask to apply to regions before extracting signals. + Check valid options by calling :func:`junifer.data.masks.list_masks` + (default None). + name : str, optional + The name of the marker. If None, it will use the class name + (default None). + + """ + + def __init__( + self, + coords: str, + radius: Optional[float] = None, + use_afni: Optional[bool] = None, + reho_params: Optional[Dict] = None, + agg_method: str = "mean", + agg_method_params: Optional[Dict] = None, + mask: Optional[str] = None, + name: Optional[str] = None, + ) -> None: + self.coords = coords + self.radius = radius + self.reho_params = reho_params + self.agg_method = agg_method + self.agg_method_params = agg_method_params + self.mask = mask + super().__init__(use_afni=use_afni, name=name) + + def compute( + self, + input: Dict[str, Any], + extra_input: Optional[Dict[str, Any]] = None, + ) -> Dict[str, Any]: + """Compute. + + Parameters + ---------- + input : dict + The BOLD data as dictionary. + extra_input : dict, optional + The other fields in the pipeline data object (default None). + + Returns + ------- + dict + The computed result as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a 1D numpy.ndarray + * ``columns`` : the column labels for the spheres as a list + * ``rows_col_name`` : ``None`` + + """ + logger.info("Calculating ReHo for spheres.") + # Calculate reho map + if self.reho_params is not None: + reho_map = self.compute_reho_map(input=input, **self.reho_params) + else: + reho_map = self.compute_reho_map(input=input) + # Initialize sphere aggregation + sphere_aggregation = SphereAggregation( + coords=self.coords, + radius=self.radius, + method=self.agg_method, + method_params=self.agg_method_params, + mask=self.mask, + on="BOLD", + ) + # Perform aggregation on reho map + sphere_aggregation_input = {"data": reho_map} + output = sphere_aggregation.compute(input=sphere_aggregation_input) + # Only use the first row and expand row dimension + output["data"] = output["data"][0][np.newaxis, :] + return output diff --git a/junifer/markers/reho/tests/test_reho_estimator.py b/junifer/markers/reho/tests/test_reho_estimator.py new file mode 100644 index 000000000..ad6f3649a --- /dev/null +++ b/junifer/markers/reho/tests/test_reho_estimator.py @@ -0,0 +1,260 @@ +"""Provide tests for ReHo map compute comparison.""" + +# Authors: Synchon Mandal +# License: AGPL + +import time + +import nibabel as nib +import pytest +from scipy.stats import pearsonr + +from junifer.datareader.default import DefaultDataReader +from junifer.markers.reho.reho_estimator import ReHoEstimator +from junifer.pipeline.utils import _check_afni +from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber +from junifer.utils.logging import logger + + +def test_reho_estimator_cache_python() -> None: + """Test that the cache works properly when using Python implementation.""" + # Get subject from datagrabber + with PartlyCloudyTestingDataGrabber() as dg: + subject = dg["sub-01"] + # Read data for subject + subject_data = DefaultDataReader().fit_transform(subject) + # Setup estimator + reho_estimator = ReHoEstimator() + + first_tic = time.time() + reho_map_without_cache = reho_estimator.fit_transform( + use_afni=False, + input_data=subject_data["BOLD"], + nneigh=27, + ) + first_toc = time.time() + logger.info( + f"ReHo estimator in Python without cache: {first_toc - first_tic}" + ) + assert isinstance(reho_map_without_cache, nib.Nifti1Image) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + # Now fit again, should be faster + second_tic = time.time() + reho_map_with_cache = reho_estimator.fit_transform( + use_afni=False, + input_data=subject_data["BOLD"], + nneigh=27, + ) + second_toc = time.time() + logger.info( + f"ReHo estimator in Python with cache: {second_toc - second_tic}" + ) + assert isinstance(reho_map_with_cache, nib.Nifti1Image) + # Check that cache is being used + assert (second_toc - second_tic) < ((first_toc - first_tic) / 1000) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + # Now change a parameter, should compute again, without clearing the + # cache + third_tic = time.time() + reho_map_with_partial_cache = reho_estimator.fit_transform( + use_afni=False, + input_data=subject_data["BOLD"], + nneigh=125, + ) + third_toc = time.time() + logger.info( + f"ReHo estimator in Python with partial cache: {third_toc - third_tic}" + ) + assert isinstance(reho_map_with_partial_cache, nib.Nifti1Image) + # Should require more time + assert (third_toc - third_tic) > ((first_toc - first_tic) / 10) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + # Now fit again with the previous params, should be fast + fourth_tic = time.time() + reho_map_with_new_cache = reho_estimator.fit_transform( + use_afni=False, + input_data=subject_data["BOLD"], + nneigh=125, + ) + fourth_toc = time.time() + logger.info( + f"ReHo estimator in Python with new cache: {fourth_toc - fourth_tic}" + ) + assert isinstance(reho_map_with_new_cache, nib.Nifti1Image) + # Should require less time + assert (fourth_toc - fourth_tic) < ((first_toc - first_tic) / 1000) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + # Now change the data, it should clear the cache + with PartlyCloudyTestingDataGrabber() as dg: + subject = dg["sub-02"] + # Read data for new subject + subject_data = DefaultDataReader().fit_transform(subject) + + fifth_tic = time.time() + reho_map_with_different_cache = reho_estimator.fit_transform( + use_afni=False, + input_data=subject_data["BOLD"], + nneigh=27, + ) + fifth_toc = time.time() + logger.info( + "ReHo estimator in Python with different cache: " + f"{fifth_toc - fifth_tic}" + ) + assert isinstance(reho_map_with_different_cache, nib.Nifti1Image) + # Should take less time + assert (fifth_toc - fifth_tic) > ((first_toc - first_tic) / 10) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +def test_reho_estimator_cache_afni() -> None: + """Test that the cache works properly when using afni.""" + # Get subject from datagrabber + with PartlyCloudyTestingDataGrabber() as dg: + subject = dg["sub-01"] + # Read data for subject + subject_data = DefaultDataReader().fit_transform(subject) + # Setup estimator + reho_estimator = ReHoEstimator() + + first_tic = time.time() + reho_map_without_cache = reho_estimator.fit_transform( + use_afni=True, + input_data=subject_data["BOLD"], + nneigh=19, + ) + first_toc = time.time() + logger.info( + f"ReHo estimator in AFNI without cache: {first_toc - first_tic}" + ) + assert isinstance(reho_map_without_cache, nib.Nifti1Image) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 2 # input + reho + + # Now fit again, should be faster + second_tic = time.time() + reho_map_with_cache = reho_estimator.fit_transform( + use_afni=True, + input_data=subject_data["BOLD"], + nneigh=19, + ) + second_toc = time.time() + logger.info( + f"ReHo estimator in AFNI with cache: {second_toc - second_tic}" + ) + assert isinstance(reho_map_with_cache, nib.Nifti1Image) + assert (second_toc - second_tic) < ((first_toc - first_tic) / 1000) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 2 # input + reho + + # Now change a parameter, should compute again, without clearing the + # cache + third_tic = time.time() + reho_map_with_partial_cache = reho_estimator.fit_transform( + use_afni=True, + input_data=subject_data["BOLD"], + nneigh=27, + ) + third_toc = time.time() + logger.info( + f"ReHo estimator in AFNI with partial cache: {third_toc - third_tic}" + ) + assert isinstance(reho_map_with_partial_cache, nib.Nifti1Image) + # Should require more time + assert (third_toc - third_tic) > ((first_toc - first_tic) / 10) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 3 # input + 2 * reho + + # Now fit again with the previous params, should be fast + fourth_tic = time.time() + reho_map_with_new_cache = reho_estimator.fit_transform( + use_afni=True, + input_data=subject_data["BOLD"], + nneigh=27, + ) + fourth_toc = time.time() + logger.info( + f"ReHo estimator in AFNI with new cache: {fourth_toc - fourth_tic}" + ) + assert isinstance(reho_map_with_new_cache, nib.Nifti1Image) + # Should require less time + assert (fourth_toc - fourth_tic) < ((first_toc - first_tic) / 1000) + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 3 # input + 2 * reho + + # Now change the data, it should clear the cache + with PartlyCloudyTestingDataGrabber() as dg: + subject = dg["sub-02"] + # Read data for new subject + subject_data = DefaultDataReader().fit_transform(subject) + + fifth_tic = time.time() + reho_map_with_different_cache = reho_estimator.fit_transform( + use_afni=True, + input_data=subject_data["BOLD"], + nneigh=27, + ) + fifth_toc = time.time() + logger.info( + f"ReHo estimator in AFNI with different cache: {fifth_toc - fifth_tic}" + ) + assert isinstance(reho_map_with_different_cache, nib.Nifti1Image) + # Should take less time + assert (fifth_toc - fifth_tic) > ((first_toc - first_tic) / 10) + # Count intermediate files + n_files = len([x for x in reho_estimator.temp_dir_path.glob("*")]) + assert n_files == 2 # input + reho + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +def test_reho_estimator_afni_vs_python() -> None: + """Compare afni and Python implementations.""" + # Get subject from datagrabber + with PartlyCloudyTestingDataGrabber() as dg: + subject = dg["sub-01"] + # Read data for subject + subject_data = DefaultDataReader().fit_transform(subject) + # Setup estimator + reho_estimator = ReHoEstimator() + + # Compare using 27 neighbours + reho_map_afni = reho_estimator.fit_transform( + use_afni=True, + input_data=subject_data["BOLD"], + nneigh=27, + ) + reho_map_python = reho_estimator.fit_transform( + use_afni=False, + input_data=subject_data["BOLD"], + nneigh=27, + ) + + # Calculate Pearson correlation coefficient + r, _ = pearsonr( + reho_map_afni.get_fdata().flatten(), + reho_map_python.get_fdata().flatten(), + ) + # Assert good correlation + assert r > 0.70 diff --git a/junifer/markers/reho/tests/test_reho_parcels.py b/junifer/markers/reho/tests/test_reho_parcels.py new file mode 100644 index 000000000..ef3f6bcb3 --- /dev/null +++ b/junifer/markers/reho/tests/test_reho_parcels.py @@ -0,0 +1,118 @@ +"""Provide tests for ReHo on parcels.""" + +# Authors: Synchon Mandal +# License: AGPL + +from pathlib import Path + +import pytest +from nilearn import image as nimg +from scipy.stats import pearsonr + +from junifer.markers.reho.reho_parcels import ReHoParcels +from junifer.pipeline.utils import _check_afni +from junifer.storage.sqlite import SQLiteFeatureStorage +from junifer.testing.datagrabbers import SPMAuditoryTestingDatagrabber + + +PARCELLATION = "Schaefer100x7" + + +def test_reho_parcels_computation() -> None: + """Test ReHoParcels fit-transform.""" + with SPMAuditoryTestingDatagrabber() as dg: + # Use first subject + subject_data = dg["sub001"] + # Load image to memory + fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) + # Initialize marker + reho_parcels_marker = ReHoParcels(parcellation=PARCELLATION) + # Fit transform marker on data + reho_parcels_output = reho_parcels_marker.fit_transform( + {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + ) + # Get BOLD output + reho_parcels_output_bold = reho_parcels_output["BOLD"] + # Assert BOLD output keys + assert "data" in reho_parcels_output_bold + assert "columns" in reho_parcels_output_bold + + reho_parcels_output_bold_data = reho_parcels_output_bold["data"] + # Assert BOLD output data dimension + assert reho_parcels_output_bold_data.ndim == 2 + # Assert BOLD output data is normalized + assert (reho_parcels_output_bold_data > 0).all() and ( + reho_parcels_output_bold_data < 1 + ).all() + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +def test_reho_parcels_computation_comparison() -> None: + """Test ReHoParcels fit-transform implementation comparison..""" + with SPMAuditoryTestingDatagrabber() as dg: + # Use first subject + subject_data = dg["sub001"] + # Load image to memory + fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) + + # Initialize marker with use_afni=False + reho_parcels_marker_python = ReHoParcels( + parcellation=PARCELLATION, use_afni=False + ) + # Fit transform marker on data + reho_parcels_output_python = reho_parcels_marker_python.fit_transform( + {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + ) + # Get BOLD output + reho_parcels_output_bold_python = reho_parcels_output_python["BOLD"] + + # Initialize marker with use_afni=True + reho_parcels_marker_afni = ReHoParcels( + parcellation=PARCELLATION, use_afni=True + ) + # Fit transform marker on data + reho_parcels_output_afni = reho_parcels_marker_afni.fit_transform( + {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + ) + # Get BOLD output + reho_parcels_output_bold_afni = reho_parcels_output_afni["BOLD"] + + # Check for Pearson correlation coefficient + r, _ = pearsonr( + reho_parcels_output_bold_python["data"].flatten(), + reho_parcels_output_bold_afni["data"].flatten(), + ) + assert r >= 0.3 # this is very bad, but they differ... + + +def test_reho_parcels_storage(tmp_path: Path) -> None: + """Test ReHoParcels storage. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + with SPMAuditoryTestingDatagrabber() as dg: + # Use first subject + subject_data = dg["sub001"] + # Load image to memory + fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) + # Initialize marker + reho_parcels_marker = ReHoParcels(parcellation=PARCELLATION) + # Initialize storage + reho_parcels_storage = SQLiteFeatureStorage( + tmp_path / "reho_parcels.sqlite" + ) + # Generate meta + meta = { + "element": {"subject": "sub001"} + } # only requires element key for storing + # Fit transform marker on data with storage + reho_parcels_marker.fit_transform( + input={"BOLD": {"path": "/tmp", "data": fmri_img, "meta": meta}}, + storage=reho_parcels_storage, + ) diff --git a/junifer/markers/reho/tests/test_reho_spheres.py b/junifer/markers/reho/tests/test_reho_spheres.py new file mode 100644 index 000000000..d737cbf23 --- /dev/null +++ b/junifer/markers/reho/tests/test_reho_spheres.py @@ -0,0 +1,118 @@ +"""Provide tests for ReHo on spheres.""" + +# Authors: Synchon Mandal +# License: AGPL + +from pathlib import Path + +import pytest +from nilearn import image as nimg +from scipy.stats import pearsonr + +from junifer.markers.reho.reho_spheres import ReHoSpheres +from junifer.pipeline.utils import _check_afni +from junifer.storage.sqlite import SQLiteFeatureStorage +from junifer.testing.datagrabbers import SPMAuditoryTestingDatagrabber + + +COORDINATES = "DMNBuckner" + + +def test_reho_spheres_computation() -> None: + """Test ReHoSpheres fit-transform.""" + with SPMAuditoryTestingDatagrabber() as dg: + # Use first subject + subject_data = dg["sub001"] + # Load image to memory + fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) + # Initialize marker + reho_spheres_marker = ReHoSpheres(coords=COORDINATES, radius=10.0) + # Fit transform marker on data + reho_spheres_output = reho_spheres_marker.fit_transform( + {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + ) + # Get BOLD output + reho_spheres_output_bold = reho_spheres_output["BOLD"] + # Assert BOLD output keys + assert "data" in reho_spheres_output_bold + assert "columns" in reho_spheres_output_bold + + reho_spheres_output_bold_data = reho_spheres_output_bold["data"] + # Assert BOLD output data dimension + assert reho_spheres_output_bold_data.ndim == 2 + # Assert BOLD output data is normalized + assert (reho_spheres_output_bold_data > 0).all() and ( + reho_spheres_output_bold_data < 1 + ).all() + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +def test_reho_spheres_computation_comparison() -> None: + """Test ReHoSpheres fit-transform implementation comparison..""" + with SPMAuditoryTestingDatagrabber() as dg: + # Use first subject + subject_data = dg["sub001"] + # Load image to memory + fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) + + # Initialize marker with use_afni=False + reho_spheres_marker_python = ReHoSpheres( + coords=COORDINATES, radius=10.0, use_afni=False + ) + # Fit transform marker on data + reho_spheres_output_python = reho_spheres_marker_python.fit_transform( + {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + ) + # Get BOLD output + reho_spheres_output_bold_python = reho_spheres_output_python["BOLD"] + + # Initialize marker with use_afni=True + reho_spheres_marker_afni = ReHoSpheres( + coords=COORDINATES, radius=10.0, use_afni=True + ) + # Fit transform marker on data + reho_spheres_output_afni = reho_spheres_marker_afni.fit_transform( + {"BOLD": {"path": "/tmp", "data": fmri_img, "meta": {}}} + ) + # Get BOLD output + reho_spheres_output_bold_afni = reho_spheres_output_afni["BOLD"] + + # Check for Pearson correlation coefficient + r, _ = pearsonr( + reho_spheres_output_bold_python["data"].flatten(), + reho_spheres_output_bold_afni["data"].flatten(), + ) + assert r >= 0.8 # 0.8 is a loose threshold + + +def test_reho_spheres_storage(tmp_path: Path) -> None: + """Test ReHoSpheres storage. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + + """ + with SPMAuditoryTestingDatagrabber() as dg: + # Use first subject + subject_data = dg["sub001"] + # Load image to memory + fmri_img = nimg.load_img(subject_data["BOLD"]["path"]) + # Initialize marker + reho_spheres_marker = ReHoSpheres(coords=COORDINATES, radius=10.0) + # Initialize storage + reho_spheres_storage = SQLiteFeatureStorage( + tmp_path / "reho_spheres.sqlite" + ) + # Generate meta + meta = { + "element": {"subject": "sub001"} + } # only requires element key for storing + # Fit transform marker on data with storage + reho_spheres_marker.fit_transform( + input={"BOLD": {"path": "/tmp", "data": fmri_img, "meta": meta}}, + storage=reho_spheres_storage, + ) diff --git a/junifer/markers/tests/test_crossparcellation_functional_connectivity.py b/junifer/markers/tests/test_crossparcellation_functional_connectivity.py index 8e9fd53da..485a4900a 100644 --- a/junifer/markers/tests/test_crossparcellation_functional_connectivity.py +++ b/junifer/markers/tests/test_crossparcellation_functional_connectivity.py @@ -71,8 +71,7 @@ def test_store(tmp_path: Path) -> None: crossparcellation.fit_transform(input_dict, storage=storage) features = storage.list_features() assert any( - x["name"] == "BOLD_CrossParcellationFC" - for x in features.values() + x["name"] == "BOLD_CrossParcellationFC" for x in features.values() ) diff --git a/junifer/markers/tests/test_sphere_aggregation.py b/junifer/markers/tests/test_sphere_aggregation.py index 5081fb324..0d373faa5 100644 --- a/junifer/markers/tests/test_sphere_aggregation.py +++ b/junifer/markers/tests/test_sphere_aggregation.py @@ -4,8 +4,8 @@ # License: AGPL import typing -from typing import Dict from pathlib import Path +from typing import Dict import nibabel as nib import pytest diff --git a/junifer/markers/utils.py b/junifer/markers/utils.py index 340cd5d77..4c009b5c7 100644 --- a/junifer/markers/utils.py +++ b/junifer/markers/utils.py @@ -29,7 +29,7 @@ def singleton(cls: Type) -> Type: class The only instance of the class. - """ "" + """ instances: Dict = {} def get_instance(*args: Any, **kwargs: Any) -> Type: @@ -47,7 +47,7 @@ def singleton(cls: Type) -> Type: class The only instance of the class. - """ "" + """ if cls not in instances: instances[cls] = cls(*args, **kwargs) return instances[cls] diff --git a/junifer/pipeline/tests/test_pipeline_step_mixin.py b/junifer/pipeline/tests/test_pipeline_step_mixin.py index 5631baf7a..c517d598f 100644 --- a/junifer/pipeline/tests/test_pipeline_step_mixin.py +++ b/junifer/pipeline/tests/test_pipeline_step_mixin.py @@ -4,10 +4,10 @@ # Synchon Mandal # License: AGPL +import warnings from typing import Dict, List import pytest -import warnings from junifer.pipeline.pipeline_step_mixin import PipelineStepMixin from junifer.pipeline.utils import _check_afni diff --git a/junifer/storage/sqlite.py b/junifer/storage/sqlite.py index 3d9ed05c6..99428da28 100644 --- a/junifer/storage/sqlite.py +++ b/junifer/storage/sqlite.py @@ -4,10 +4,10 @@ # Synchon Mandal # License: AGPL +import json from pathlib import Path from typing import TYPE_CHECKING, Dict, List, Optional, Union -import json import numpy as np import pandas as pd from pandas.core.base import NoNewAttributesMixin diff --git a/tox.ini b/tox.ini index 5be4059a9..47b352ada 100644 --- a/tox.ini +++ b/tox.ini @@ -93,6 +93,10 @@ exclude = __init__.py max-line-length = 79 extend-ignore = + ; Use of `functools.lru_cache` or `functools.cache` on methods can lead to + ; memory leaks. The cache may retain instance references, preventing garbage + ; collection. + B019 ; abstract class with no abstract methods B024 D202