From e5e39f3f93726952ee23557026a9aba4c8a6bdea Mon Sep 17 00:00:00 2001 From: Fede Date: Thu, 24 Nov 2022 16:47:15 +0100 Subject: [PATCH 01/10] WIP: add fALFF --- junifer/markers/__init__.py | 1 + junifer/markers/falff/__init__.py | 6 + junifer/markers/falff/falff_base.py | 223 +++++++++++++++++++++++++ junifer/markers/falff/falff_parcels.py | 123 ++++++++++++++ junifer/markers/falff/falff_spheres.py | 130 ++++++++++++++ 5 files changed, 483 insertions(+) create mode 100644 junifer/markers/falff/__init__.py create mode 100644 junifer/markers/falff/falff_base.py create mode 100644 junifer/markers/falff/falff_parcels.py create mode 100644 junifer/markers/falff/falff_spheres.py diff --git a/junifer/markers/__init__.py b/junifer/markers/__init__.py index 2c853aad0..2af03f4ea 100644 --- a/junifer/markers/__init__.py +++ b/junifer/markers/__init__.py @@ -13,3 +13,4 @@ from .functional_connectivity_spheres import FunctionalConnectivitySpheres from .parcel_aggregation import ParcelAggregation from .sphere_aggregation import SphereAggregation from .reho import ReHoParcels, ReHoSpheres +from .falff import AmplitudeLowFrequencyFluctuationParcels diff --git a/junifer/markers/falff/__init__.py b/junifer/markers/falff/__init__.py new file mode 100644 index 000000000..dcbb4d18c --- /dev/null +++ b/junifer/markers/falff/__init__.py @@ -0,0 +1,6 @@ +"""Provide imports for falff sub-package.""" + +# Authors: Federico Raimondo +# License: AGPL + +from .falff_parcels import AmplitudeLowFrequencyFluctuationParcels diff --git a/junifer/markers/falff/falff_base.py b/junifer/markers/falff/falff_base.py new file mode 100644 index 000000000..bcf714c88 --- /dev/null +++ b/junifer/markers/falff/falff_base.py @@ -0,0 +1,223 @@ +"""Provide abstract class for computing fALFF.""" +# Authors: Federico Raimondo +# Amir Omidvarnia +# Kaustubh R. Patil +# License: AGPL + +from typing import TYPE_CHECKING, Dict, List, Optional, Tuple + +import numpy as np +from scipy import signal + +from junifer.markers.base import BaseMarker +from junifer.utils import logger + +if TYPE_CHECKING: + from junifer.storage import BaseFeatureStorage + + +class AmplitudeLowFrequencyFluctuationBase(BaseMarker): + """Class for computing fALFF/ALFF. + + Parameters + ---------- + tr : float, optional + The Repetition Time of the BOLD data. If None, will extract + the TR from NIFTI header (default None). + name : str, optional + The name of the marker. If None, will use the class name (default + None). + + Notes + ----- + The `tr` parameter is crucial for the correctness of fALFF/ALFF + computation. If a dataset is correctly preprocessed, the TR should be + extracted from the NIFTI without any issue. However, it has been + reported that some preprocessed data might not have the correct TR in + the NIFTI header. + """ + + def __init__( + self, + highpass: float, + lowpass: float, + order: int, + tr: Optional[float] = None, + name: Optional[str] = None, + ) -> None: + if highpass <= 0: + raise ValueError("Highpass must be positive") + if lowpass <= 0: + raise ValueError("Lowpass must be positive") + if highpass >= lowpass: + raise ValueError("Highpass must be lower than lowpass") + self.highpass = highpass + self.lowpass = lowpass + if order <= 0: + raise ValueError("Order must be positive") + self.order = order + self.tr = tr + + super().__init__(on="BOLD", name=name) + + 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. + + """ + valid = ["BOLD"] + return valid + + def get_output_kind(self, input: List[str]) -> List[str]: + """Get output kind. + + Parameters + ---------- + input : list of str + The kind of data to work on. + + Returns + ------- + list of str + The list of storage kinds. + + """ + return ["table"] + + def get_meta(self, kind: str, fractional: bool) -> Dict: + """Get metadata. + + Parameters + ---------- + kind : str + The kind of pipeline step. + fractional : bool + If true, get the meta for fractional ALFF. + + Returns + ------- + dict + The metadata as a dictionary with the only key 'marker'. + + """ + s_meta = super().get_meta(kind) + suffix = "_fALFF" if fractional else "_ALFF" + s_meta["marker"]["name"] += suffix + return s_meta + + def fit_transform( + self, + input: Dict[str, Dict], + storage: Optional["BaseFeatureStorage"] = None, + ) -> Dict: + """Fit and transform. + + Parameters + ---------- + input : dict + The Junifer Data object. + storage : storage-like, optional + The storage class, for example, SQLiteFeatureStorage. + + Returns + ------- + dict + The processed output as a dictionary. If `storage` is provided, + empty dictionary is returned. + + """ + out = {} + t_meta = input.get("meta", {}).copy() + logger.info("Computing fALFF on BOLD") + t_input = input["BOLD"] + t_meta.update(t_input.get("meta", {})) + + # Compute ALFF and fALFF + alff, falff = self.compute(input=t_input) + + alff_meta = t_meta.copy().update( + self.get_meta("BOLD", fractional=False) + ) + falff_meta = t_meta.copy().update( + self.get_meta("BOLD", fractional=True) + ) + + alff.update(meta=alff_meta) + falff.update(meta=falff_meta) + if storage is not None: + logger.info(f"Storing in {storage}") + self.store(kind="BOLD", out=alff, storage=storage) + self.store(kind="BOLD", out=falff, storage=storage) + else: + logger.info("No storage specified, returning dictionary") + out["BOLD"] = { + "ALFF": alff, + "fALFF": falff, + } + return out + + def compute_falff( + self, timeseries: np.ndarray, labels: List[str], tr: float + ) -> Tuple[Dict, Dict]: + """Compute ALFF and fALFF. + + Parameters + ---------- + timeseries : np.ndarray + The timeseries data. + labels : np.ndarray + The labels for each timeseries. + tr : float + The repetition time. + + Returns + ------- + alff: dict + The computed ALFF as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a numpy.ndarray + * ``columns`` : the column labels for the computed values as a list + + falff: dict + The computed fALFF as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a numpy.ndarray + * ``columns`` : the column labels for the computed values as a list + """ + # bandpass the data within the lowpass and highpass cutoff freqs + Nq = 1 / (2 * tr) + Wn = np.array([self.highpass / Nq, self.lowpass / Nq]) + + b, a = signal.butter(N=self.order, Wn=Wn, btype="bandpass") + ts_filt = signal.filtfilt(b, a, timeseries, axis=0) + + ALFF = np.std(ts_filt, axis=0) + PSD_tot = np.std(timeseries, axis=0) + + fALFF = np.divide(ALFF, PSD_tot) + + out_alff = {"data": ALFF, "columns": labels} + out_falff = {"data": fALFF, "columns": labels} + + return out_alff, out_falff + + def store(self, kind, out, storage): + """Store. + + Parameters + ---------- + kind : {"BOLD"} + The data kind to store. + out : dict + The computed result as a dictionary to store. + storage : storage-like + The storage class, for example, SQLiteFeatureStorage. + """ + logger.debug(f"Storing {kind} in {storage}") + storage.store(kind="table", **out) diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py new file mode 100644 index 000000000..cd90ef1ad --- /dev/null +++ b/junifer/markers/falff/falff_parcels.py @@ -0,0 +1,123 @@ +"""Provide class for computing fALFF on parcels.""" +# Authors: Federico Raimondo +# Amir Omidvarnia +# Kaustubh R. Patil +# License: AGPL + +from typing import Dict, List, Optional, Union, Tuple + + +from ...api.decorators import register_marker +from .falff_base import AmplitudeLowFrequencyFluctuationBase +from .. import ParcelAggregation + + +@register_marker +class AmplitudeLowFrequencyFluctuationParcels( + AmplitudeLowFrequencyFluctuationBase +): + """Class for computing fALFF/ALFF on parcels. + + Parameters + ---------- + parcellation : str or list of str + The name(s) of the parcellation(s). Check valid options by calling + :func:`junifer.data.parcellations.list_parcellations`. + highpass : float + The highpass cutoff frequency for the bandpass filter. + lowpass : float + The lowpass cutoff frequency for the bandpass filter. + order : int + The order of the bandpass filter. + 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). + tr : float, optional + The Repetition Time of the BOLD data. If None, will extract + the TR from NIFTI header (default None). + name : str, optional + The name of the marker. If None, will use the class name (default + None). + + Notes + ----- + The `tr` parameter is crucial for the correctness of fALFF/ALFF + computation. If a dataset is correctly preprocessed, the TR should be + extracted from the NIFTI without any issue. However, it has been + reported that some preprocessed data might not have the correct TR in + the NIFTI header. + + ALFF/fALFF are computed using a bandpass butterworth filter. See + :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more + details. + + + """ + + def __init__( + self, + parcellation: Union[str, List[str]], + highpass: float, + lowpass: float, + order: int, + mask: Optional[str] = None, + tr: Optional[float] = None, + name: Optional[str] = None, + ) -> None: + self.parcellation = parcellation + self.mask = mask + + super().__init__( + highpass=highpass, + lowpass=lowpass, + order=order, + tr=tr, + name=name, + ) + + def compute(self, input: Dict) -> Tuple[Dict, Dict]: + """Compute ALFF and fALFF. + + Parameters + ---------- + input : dict + A single input from the pipeline data object in which to compute + the marker. + extra_input : dict, optional + The other fields in the pipeline data object. Useful for accessing + other data kind that needs to be used in the computation. For + example, the functional connectivity markers can make use of the + confounds if available (default None). + + Returns + ------- + alff: dict + The computed ALFF as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a numpy.ndarray + * ``columns`` : the column labels for the computed values as a list + + falff: dict + The computed fALFF as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a numpy.ndarray + * ``columns`` : the column labels for the computed values as a list + """ + pa = ParcelAggregation( + parcellation=self.parcellation, + method="mean", + mask=self.mask, + on="BOLD", + ) + + # get the 2D timeseries after parcel aggregation + parcels = pa.compute(input) + timeseries = parcels["data"] + labels = parcels["columns"] + tr = self.tr or input["BOLD"].header["pixdim"][4] + + out = super().compute_falff(timeseries, labels, tr) + return out diff --git a/junifer/markers/falff/falff_spheres.py b/junifer/markers/falff/falff_spheres.py new file mode 100644 index 000000000..7f7ee9f7c --- /dev/null +++ b/junifer/markers/falff/falff_spheres.py @@ -0,0 +1,130 @@ +"""Provide class for computing fALFF on spheres.""" +# Authors: Federico Raimondo +# Amir Omidvarnia +# Kaustubh R. Patil +# License: AGPL + +from typing import Dict, Optional, Tuple + + +from ...api.decorators import register_marker +from .falff_base import AmplitudeLowFrequencyFluctuationBase +from .. import SphereAggregation + + +@register_marker +class AmplitudeLowFrequencyFluctuationSpheres( + AmplitudeLowFrequencyFluctuationBase +): + """Class for computing fALFF/ALFF 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 mm. If None, the signal will be extracted + from a single voxel. See :class:`nilearn.maskers.NiftiSpheresMasker` + for more information (default None). + highpass : float + The highpass cutoff frequency for the bandpass filter. + lowpass : float + The lowpass cutoff frequency for the bandpass filter. + order : int + The order of the bandpass filter. + 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). + tr : float, optional + The Repetition Time of the BOLD data. If None, will extract + the TR from NIFTI header (default None). + name : str, optional + The name of the marker. If None, will use the class name (default + None). + + Notes + ----- + The `tr` parameter is crucial for the correctness of fALFF/ALFF + computation. If a dataset is correctly preprocessed, the TR should be + extracted from the NIFTI without any issue. However, it has been + reported that some preprocessed data might not have the correct TR in + the NIFTI header. + + ALFF/fALFF are computed using a bandpass butterworth filter. See + :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more + details. + + + """ + + def __init__( + self, + coords: str, + highpass: float, + lowpass: float, + order: int, + radius: Optional[float] = None, + mask: Optional[str] = None, + tr: Optional[float] = None, + name: Optional[str] = None, + ) -> None: + self.coords = coords + self.radius = radius + self.mask = mask + + super().__init__( + highpass=highpass, + lowpass=lowpass, + order=order, + tr=tr, + name=name, + ) + + def compute(self, input: Dict) -> Tuple[Dict, Dict]: + """Compute ALFF and fALFF. + + Parameters + ---------- + input : dict + A single input from the pipeline data object in which to compute + the marker. + extra_input : dict, optional + The other fields in the pipeline data object. Useful for accessing + other data kind that needs to be used in the computation. For + example, the functional connectivity markers can make use of the + confounds if available (default None). + + Returns + ------- + alff: dict + The computed ALFF as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a numpy.ndarray + * ``columns`` : the column labels for the computed values as a list + + falff: dict + The computed fALFF as dictionary. The dictionary has the following + keys: + + * ``data`` : the actual computed values as a numpy.ndarray + * ``columns`` : the column labels for the computed values as a list + """ + pa = SphereAggregation( + coords=self.coords, + radius=self.radius, + method="mean", + mask=self.mask, + on="BOLD", + ) + + # get the 2D timeseries after parcel aggregation + spheres = pa.compute(input) + timeseries = spheres["data"] + labels = spheres["columns"] + tr = self.tr or input["BOLD"].header["pixdim"][4] + + out = super().compute_falff(timeseries, labels, tr) + return out -- 2.52.0 From 38354bd98a06dc668dff2557b148dcd9fad70b97 Mon Sep 17 00:00:00 2001 From: Fede Date: Thu, 24 Nov 2022 17:14:56 +0100 Subject: [PATCH 02/10] add default values for alff filters --- junifer/markers/falff/falff_parcels.py | 14 +++++++------- junifer/markers/falff/falff_spheres.py | 14 +++++++------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py index cd90ef1ad..6f4ba1d67 100644 --- a/junifer/markers/falff/falff_parcels.py +++ b/junifer/markers/falff/falff_parcels.py @@ -23,12 +23,12 @@ class AmplitudeLowFrequencyFluctuationParcels( parcellation : str or list of str The name(s) of the parcellation(s). Check valid options by calling :func:`junifer.data.parcellations.list_parcellations`. - highpass : float - The highpass cutoff frequency for the bandpass filter. + highpass : float, optional + The highpass cutoff frequency for the bandpass filter (default 0.01). lowpass : float - The lowpass cutoff frequency for the bandpass filter. + The lowpass cutoff frequency for the bandpass filter (default 0.1). order : int - The order of the bandpass filter. + The order of the bandpass filter (default 4). 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` @@ -58,9 +58,9 @@ class AmplitudeLowFrequencyFluctuationParcels( def __init__( self, parcellation: Union[str, List[str]], - highpass: float, - lowpass: float, - order: int, + highpass: float = 0.01, + lowpass: float = 0.1, + order: int = 4, mask: Optional[str] = None, tr: Optional[float] = None, name: Optional[str] = None, diff --git a/junifer/markers/falff/falff_spheres.py b/junifer/markers/falff/falff_spheres.py index 7f7ee9f7c..33b517429 100644 --- a/junifer/markers/falff/falff_spheres.py +++ b/junifer/markers/falff/falff_spheres.py @@ -27,12 +27,12 @@ class AmplitudeLowFrequencyFluctuationSpheres( 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). - highpass : float - The highpass cutoff frequency for the bandpass filter. + highpass : float, optional + The highpass cutoff frequency for the bandpass filter (default 0.01). lowpass : float - The lowpass cutoff frequency for the bandpass filter. + The lowpass cutoff frequency for the bandpass filter (default 0.1). order : int - The order of the bandpass filter. + The order of the bandpass filter (default 4). 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` @@ -62,9 +62,9 @@ class AmplitudeLowFrequencyFluctuationSpheres( def __init__( self, coords: str, - highpass: float, - lowpass: float, - order: int, + highpass: float = 0.01, + lowpass: float = 0.1, + order: int = 4, radius: Optional[float] = None, mask: Optional[str] = None, tr: Optional[float] = None, -- 2.52.0 From 757c084d00bf860c52d9de522dd7e4821f017762 Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Thu, 15 Dec 2022 13:13:44 +0100 Subject: [PATCH 03/10] fALFF wip --- junifer/api/res/afni/3dRSFC | 3 + junifer/markers/falff/falff_base.py | 250 +++++++------- junifer/markers/falff/falff_estimator.py | 319 ++++++++++++++++++ junifer/markers/falff/falff_parcels.py | 50 +-- .../markers/falff/tests/test_falff_parcels.py | 61 ++++ 5 files changed, 525 insertions(+), 158 deletions(-) create mode 100755 junifer/api/res/afni/3dRSFC create mode 100644 junifer/markers/falff/falff_estimator.py create mode 100644 junifer/markers/falff/tests/test_falff_parcels.py diff --git a/junifer/api/res/afni/3dRSFC b/junifer/api/res/afni/3dRSFC new file mode 100755 index 000000000..f3036ac16 --- /dev/null +++ b/junifer/api/res/afni/3dRSFC @@ -0,0 +1,3 @@ +#!/bin/bash + +run_afni_docker.sh 3dRSFC "$@" \ No newline at end of file diff --git a/junifer/markers/falff/falff_base.py b/junifer/markers/falff/falff_base.py index bcf714c88..ba1583daa 100644 --- a/junifer/markers/falff/falff_base.py +++ b/junifer/markers/falff/falff_base.py @@ -6,28 +6,34 @@ from typing import TYPE_CHECKING, Dict, List, Optional, Tuple -import numpy as np -from scipy import signal -from junifer.markers.base import BaseMarker -from junifer.utils import logger - -if TYPE_CHECKING: - from junifer.storage import BaseFeatureStorage +from ..base import BaseMarker +from .falff_estimator import AmplitudeLowFrequencyFluctuationEstimator +from ...utils.logging import raise_error, warn_with_log class AmplitudeLowFrequencyFluctuationBase(BaseMarker): - """Class for computing fALFF/ALFF. + """Base class for (fractional) Amplitude Low Frequency Fluctuation. Parameters ---------- + fractional : bool + Whether to compute fractional ALFF. + highpass : float + Highpass cutoff frequency. + lowpass : float + Lowpass cutoff frequency. + order : int + Order of the filter. tr : float, optional The Repetition Time of the BOLD data. If None, will extract the TR from NIFTI header (default None). + 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, will use the class name (default - None). - + The name of the marker. If None, it will use the class name + (default None). Notes ----- The `tr` parameter is crucial for the correctness of fALFF/ALFF @@ -35,31 +41,83 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): extracted from the NIFTI without any issue. However, it has been reported that some preprocessed data might not have the correct TR in the NIFTI header. + """ + _EXT_DEPENDENCIES = [ + { + "name": "afni", + "optional": True, + "commands": ["3dRSFC", "3dAFNItoNIFTI"], + }, + ] + def __init__( self, + fractional: bool, highpass: float, lowpass: float, order: int, tr: Optional[float] = None, + use_afni: Optional[bool] = None, name: Optional[str] = None, ) -> None: if highpass <= 0: - raise ValueError("Highpass must be positive") + raise_error("Highpass must be positive") if lowpass <= 0: - raise ValueError("Lowpass must be positive") + raise_error("Lowpass must be positive") if highpass >= lowpass: - raise ValueError("Highpass must be lower than lowpass") + raise_error("Highpass must be lower than lowpass") self.highpass = highpass self.lowpass = lowpass - if order <= 0: - raise ValueError("Order must be positive") + if order <= 0 and use_afni is False: + raise_error("Order must be positive") self.order = order self.tr = tr + self.use_afni = use_afni + self.fractional = fractional + # Create a name based on the class name if none is provided + if name is None: + suffix = "_fractional" if fractional else "" + name = f"{self.__class__.__name__}{suffix}" super().__init__(on="BOLD", name=name) + def validate(self, input: List[str]) -> List[str]: + """Validate the the pipeline step. + + Parameters + ---------- + input : list of str + The input to the pipeline step. + + Returns + ------- + list of str + The output of the pipeline step. + + Raises + ------ + ValueError + If the pipeline step object is missing dependencies required for + its working, if the input does not have the required data, or + if AFNI was not used and the order is not positive. + + Warns + ----- + UserWarning + If AFNI is used and the order is not 0. + """ + out = super().validate(input) + if self.use_afni is True and self.order > 0: + warn_with_log( + "AFNI will not consider the order of the filter. Set this " + "parameter to 0 to avoid this warning.") + elif self.use_afni is False and self.order <= 0: + raise_error( + "Order must be positive if AFNI is not used.") + return out + def get_valid_inputs(self) -> List[str]: """Get valid data types for input. @@ -72,152 +130,70 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): valid = ["BOLD"] return valid - def get_output_kind(self, input: List[str]) -> List[str]: - """Get output kind. + def get_output_type(self, input: List[str]) -> List[str]: + """Get output type. Parameters ---------- input : list of str - The kind of data to work on. + The type of data to work on. Returns ------- list of str - The list of storage kinds. + The list of storage types. """ return ["table"] - def get_meta(self, kind: str, fractional: bool) -> Dict: - """Get metadata. - - Parameters - ---------- - kind : str - The kind of pipeline step. - fractional : bool - If true, get the meta for fractional ALFF. - - Returns - ------- - dict - The metadata as a dictionary with the only key 'marker'. - - """ - s_meta = super().get_meta(kind) - suffix = "_fALFF" if fractional else "_ALFF" - s_meta["marker"]["name"] += suffix - return s_meta - - def fit_transform( + def compute( self, input: Dict[str, Dict], - storage: Optional["BaseFeatureStorage"] = None, + extra_input: Optional[Dict] = None, ) -> Dict: - """Fit and transform. + """Compute. Parameters ---------- input : dict - The Junifer Data object. - storage : storage-like, optional - The storage class, for example, SQLiteFeatureStorage. + A single input from the pipeline data object in which to compute + the marker. + extra_input : dict, optional + The other fields in the pipeline data object. Useful for accessing + other data kind that needs to be used in the computation. For + example, the functional connectivity markers can make use of the + confounds if available (default None). Returns ------- dict - The processed output as a dictionary. If `storage` is provided, - empty dictionary is returned. + The computed result as dictionary. This will be either returned + to the user or stored in the storage by calling the store method + with this as a parameter. The dictionary has the following keys: + + * ``data`` : the actual computed values as a numpy.ndarray + * ``columns`` : the column labels for the computed values as a list + * ``row_names`` (if more than one row is present in data): "scan" """ - out = {} - t_meta = input.get("meta", {}).copy() - logger.info("Computing fALFF on BOLD") - t_input = input["BOLD"] - t_meta.update(t_input.get("meta", {})) - - # Compute ALFF and fALFF - alff, falff = self.compute(input=t_input) - - alff_meta = t_meta.copy().update( - self.get_meta("BOLD", fractional=False) - ) - falff_meta = t_meta.copy().update( - self.get_meta("BOLD", fractional=True) + estimator = AmplitudeLowFrequencyFluctuationEstimator( + use_afni=self.use_afni ) - alff.update(meta=alff_meta) - falff.update(meta=falff_meta) - if storage is not None: - logger.info(f"Storing in {storage}") - self.store(kind="BOLD", out=alff, storage=storage) - self.store(kind="BOLD", out=falff, storage=storage) - else: - logger.info("No storage specified, returning dictionary") - out["BOLD"] = { - "ALFF": alff, - "fALFF": falff, - } + alff, falff = estimator.fit_transform( + input_data=input, + highpass=self.highpass, + lowpass=self.lowpass, + order=self.order, + tr=self.tr, + ) + post_data = falff if self.fractional else alff + + post_input = { + 'data': post_data, + 'path': None, + } + + out = self._postprocess(post_input) + return out - - def compute_falff( - self, timeseries: np.ndarray, labels: List[str], tr: float - ) -> Tuple[Dict, Dict]: - """Compute ALFF and fALFF. - - Parameters - ---------- - timeseries : np.ndarray - The timeseries data. - labels : np.ndarray - The labels for each timeseries. - tr : float - The repetition time. - - Returns - ------- - alff: dict - The computed ALFF as dictionary. The dictionary has the following - keys: - - * ``data`` : the actual computed values as a numpy.ndarray - * ``columns`` : the column labels for the computed values as a list - - falff: dict - The computed fALFF as dictionary. The dictionary has the following - keys: - - * ``data`` : the actual computed values as a numpy.ndarray - * ``columns`` : the column labels for the computed values as a list - """ - # bandpass the data within the lowpass and highpass cutoff freqs - Nq = 1 / (2 * tr) - Wn = np.array([self.highpass / Nq, self.lowpass / Nq]) - - b, a = signal.butter(N=self.order, Wn=Wn, btype="bandpass") - ts_filt = signal.filtfilt(b, a, timeseries, axis=0) - - ALFF = np.std(ts_filt, axis=0) - PSD_tot = np.std(timeseries, axis=0) - - fALFF = np.divide(ALFF, PSD_tot) - - out_alff = {"data": ALFF, "columns": labels} - out_falff = {"data": fALFF, "columns": labels} - - return out_alff, out_falff - - def store(self, kind, out, storage): - """Store. - - Parameters - ---------- - kind : {"BOLD"} - The data kind to store. - out : dict - The computed result as a dictionary to store. - storage : storage-like - The storage class, for example, SQLiteFeatureStorage. - """ - logger.debug(f"Storing {kind} in {storage}") - storage.store(kind="table", **out) diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py new file mode 100644 index 000000000..5d5825c2e --- /dev/null +++ b/junifer/markers/falff/falff_estimator.py @@ -0,0 +1,319 @@ +"""Provide estimator class for regional homogeneity (ReHo).""" + +# Authors: Synchon Mandal +# Federico Raimondo +# License: AGPL + +from typing import TYPE_CHECKING, Any, Dict, Tuple, Union, Optional + +import shutil +import subprocess +import tempfile +from functools import lru_cache +from pathlib import Path + +import nibabel as nib +import numpy as np +from scipy import signal + +from nilearn import image as nimg + +from ...utils import logger, raise_error +from ..utils import singleton + + +if TYPE_CHECKING: + from nibabel import Nifti1Image, Nifti2Image + + +@singleton +class AmplitudeLowFrequencyFluctuationEstimator: + """Estimator class for AmplitudeLowFrequencyFluctuationBase. + + This class is a singleton and is used for efficient computation of fALFF, + by caching the voxel-wise ALFF map for a given set of file path and + computation parameters. + + .. warning:: This class can only be used via + :class:`junifer.markers.falff.AmplitudeLowFrequencyFluctuationBase` + as it serves a specific purpose. + + Parameters + ---------- + use_afni : bool + Whether to use afni for computation. If False, will use python. + + """ + + def __init__(self, use_afni: bool) -> None: + self.use_afni = use_afni + 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) + + @staticmethod + def _run_afni_cmd(cmd: str) -> None: + """Run AFNI command. + + Parameters + ---------- + cmd : str + AFNI command to be executed. + + Raises + ------ + RuntimeError + If AFNI command fails. + """ + logger.info(f"AFNI command to be executed: {cmd}") + process = subprocess.run( + cmd, + stdin=subprocess.DEVNULL, + # stdout=subprocess.STDOUT, + # stderr=subprocess.STDOUT, + shell=True, + check=False, + ) + if process.returncode == 0: + logger.info( + "AFNI command succeeded with the following output: " + f"{process.stdout}" + ) + else: + raise_error( + msg="AFNI command failed with the following error: " + f"{process.stdout}", + klass=RuntimeError, + ) + + def _compute_alff_afni( + self, + data: Union["Nifti1Image", "Nifti2Image"], + highpass: float, + lowpass: float, + tr: Optional[float], + ) -> Tuple["Nifti1Image", "Nifti1Image"]: + """Compute ALFF map via afni's commands. + + Parameters + ---------- + data : 4D Niimg-like object + Images to process. + highpass : float + Highpass cutoff frequency. + lowpass : float + Lowpass cutoff frequency. + order : int + Order of the filter. Not used by AFNI. + tr : float, optional + The Repetition Time of the BOLD data. + + Returns + ------- + alff: Niimg-like object + ALFF map. + falff: Niimg-like object + fALFF map. + + Raises + ------ + RuntimeError + If the AFNI commands fails due to some issues + + """ + + # Save niimg to nii.gz + nifti_in_file_path = self.temp_dir_path / "input.nii" + nib.save(data, nifti_in_file_path) + + # Bandpass the data + falff_afni_out_path_prefix = self.temp_dir_path / "temp_falff" + + bp_cmd = ( + "3dRSFC " + f"-prefix {falff_afni_out_path_prefix.resolve()} " + f"-input {nifti_in_file_path.resolve()} " + f"-band {highpass} {lowpass} " + "-no_rsfa " + ) + if tr is not None: + bp_cmd += f"-dt {tr} " + self._run_afni_cmd(bp_cmd) + + alff_fname = self.temp_dir_path / "alff.nii" + + convert_cmd = ( + "3dAFNItoNIFTI " + f"-prefix {alff_fname.resolve()} " + f"{falff_afni_out_path_prefix}_ALFF+tlrc.BRIK " + ) + self._run_afni_cmd(convert_cmd) + + falff_fname = self.temp_dir_path / "falff.nii" + + convert_cmd = ( + "3dAFNItoNIFTI " + f"-prefix {falff_fname.resolve()} " + f"{falff_afni_out_path_prefix}_fALFF+tlrc.BRIK " + ) + self._run_afni_cmd(convert_cmd) + + # Cleanup intermediate files + for fname in self.temp_dir_path.glob("temp_falff*"): + fname.unlink() + + # Load niftis + alff_img = nib.load(alff_fname) + falff_img = nib.load(falff_fname) + + return alff_img, falff_img + + def _compute_alff_python( + self, + data: Union["Nifti1Image", "Nifti2Image"], + highpass: float, + lowpass: float, + order: int, + tr: Optional[float], + ) -> Tuple["Nifti1Image", "Nifti1Image"]: + """Compute (f)ALFF map. + + Parameters + ---------- + data : 4D Niimg-like object + Images to process. + highpass : float + Highpass cutoff frequency. + lowpass : float + Lowpass cutoff frequency. + order : int + Order of the filter. + tr : float, optional + The Repetition Time of the BOLD data. + + Returns + ------- + alff: Niimg-like object + ALFF map. + falff: Niimg-like object + fALFF map. + """ + timeseries = data.get_fdata() + tr = tr or data.header["pixdim"][4] + # bandpass the data within the lowpass and highpass cutoff freqs + Nq = 1 / (2 * tr) + Wn = np.array([highpass / Nq, lowpass / Nq]) + + b, a = signal.butter(N=order, Wn=Wn, btype="bandpass") + ts_filt = signal.filtfilt(b, a, timeseries, axis=0) + + ALFF = np.std(ts_filt, axis=0) + PSD_tot = np.std(timeseries, axis=0) + + fALFF = np.divide(ALFF, PSD_tot) + alff_img = nimg.new_img_like(data, ALFF) + falff_img = nimg.new_img_like(data, fALFF) + return alff_img, falff_img + + @lru_cache(maxsize=None, typed=True) # noqa: B019 + def _compute( + self, + data: Union["Nifti1Image", "Nifti2Image"], + highpass: float, + lowpass: float, + order: int, + tr: Optional[float], + ) -> Tuple["Nifti1Image", "Nifti1Image"]: + """Compute the ALFF map with memorization. + + Parameters + ---------- + data : 4D Niimg-like object + Images to process. + highpass : float + Highpass cutoff frequency. + lowpass : float + Lowpass cutoff frequency. + order : int + Order of the filter. + tr : float, optional + The Repetition Time of the BOLD data. + + Returns + ------- + alff: Niimg-like object + ALFF map. + falff: Niimg-like object + fALFF map. + """ + if self.use_afni: + output = self._compute_alff_afni( + data=data, + highpass=highpass, + lowpass=lowpass, + tr=tr, + ) + else: + output = self._compute_alff_python( + data, highpass=highpass, lowpass=lowpass, order=order, tr=tr + ) + return output + + def fit_transform( + self, + input_data: Dict[str, Any], + highpass: float, + lowpass: float, + order: int, + tr: Optional[float], + ) -> Tuple["Nifti1Image", "Nifti1Image"]: + """Fit and transform for the estimator. + + Parameters + ---------- + input_data : dict + The BOLD data as dictionary. + highpass : float + Highpass cutoff frequency. + lowpass : float + Lowpass cutoff frequency. + order : int + Order of the filter. + tr : float, optional + The Repetition Time of the BOLD data. + + Returns + ------- + alff: Niimg-like object + ALFF map. + falff: Niimg-like object + fALFF map. + """ + 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 fALFF 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 fALFF map cache at {self._file_path}.") + # Compute + return self._compute( + data=bold_data, + highpass=highpass, + lowpass=lowpass, + order=order, + tr=tr, + ) diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py index 6f4ba1d67..d0f04d78d 100644 --- a/junifer/markers/falff/falff_parcels.py +++ b/junifer/markers/falff/falff_parcels.py @@ -23,19 +23,30 @@ class AmplitudeLowFrequencyFluctuationParcels( parcellation : str or list of str The name(s) of the parcellation(s). Check valid options by calling :func:`junifer.data.parcellations.list_parcellations`. + fractional : bool + Whether to compute fractional ALFF. highpass : float, optional The highpass cutoff frequency for the bandpass filter (default 0.01). lowpass : float The lowpass cutoff frequency for the bandpass filter (default 0.1). order : int The order of the bandpass filter (default 4). + tr : float, optional + The Repetition Time of the BOLD data. If None, will extract + the TR from NIFTI header (default None). + use_afni : bool, optional + Whether to use AFNI for computing. If None, will use AFNI only + if available (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). - tr : float, optional - The Repetition Time of the BOLD data. If None, will extract - the TR from NIFTI header (default None). + method : str, optional + The method to perform aggregation using. Check valid options in + :func:`junifer.stats.get_aggfunc_by_name` (default "mean"). + method_params : dict, optional + Parameters to pass to the aggregation function. Check valid options in + :func:`junifer.stats.get_aggfunc_by_name`. name : str, optional The name of the marker. If None, will use the class name (default None). @@ -58,25 +69,32 @@ class AmplitudeLowFrequencyFluctuationParcels( def __init__( self, parcellation: Union[str, List[str]], + fractional: bool, highpass: float = 0.01, lowpass: float = 0.1, order: int = 4, - mask: Optional[str] = None, tr: Optional[float] = None, + use_afni: Optional[bool] = None, + mask: Optional[str] = None, + method: str = "mean", + method_params: Optional[Dict] = None, name: Optional[str] = None, ) -> None: self.parcellation = parcellation self.mask = mask - + self.method = method + self.method_params = method_params super().__init__( + fractional=fractional, highpass=highpass, lowpass=lowpass, order=order, tr=tr, name=name, + use_afni=use_afni, ) - def compute(self, input: Dict) -> Tuple[Dict, Dict]: + def _postprocess(self, input: Dict) -> Dict: """Compute ALFF and fALFF. Parameters @@ -92,32 +110,22 @@ class AmplitudeLowFrequencyFluctuationParcels( Returns ------- - alff: dict + dict The computed ALFF as dictionary. The dictionary has the following keys: - * ``data`` : the actual computed values as a numpy.ndarray - * ``columns`` : the column labels for the computed values as a list - - falff: dict - The computed fALFF as dictionary. The dictionary has the following - keys: - * ``data`` : the actual computed values as a numpy.ndarray * ``columns`` : the column labels for the computed values as a list """ pa = ParcelAggregation( parcellation=self.parcellation, - method="mean", + method=self.method, + method_params=self.method_params, mask=self.mask, - on="BOLD", + on="fALFF", ) # get the 2D timeseries after parcel aggregation - parcels = pa.compute(input) - timeseries = parcels["data"] - labels = parcels["columns"] - tr = self.tr or input["BOLD"].header["pixdim"][4] + out = pa.compute(input) - out = super().compute_falff(timeseries, labels, tr) return out diff --git a/junifer/markers/falff/tests/test_falff_parcels.py b/junifer/markers/falff/tests/test_falff_parcels.py new file mode 100644 index 000000000..887db8368 --- /dev/null +++ b/junifer/markers/falff/tests/test_falff_parcels.py @@ -0,0 +1,61 @@ +"""Provide test for parcel-aggregated (f)ALFF .""" + +# Authors: Federico Raimondo +# Synchon Mandal +# License: AGPL + +from pathlib import Path + +import nibabel as nib +import numpy as np +import pytest +from nilearn import datasets +from nilearn.image import concat_imgs, math_img, new_img_like, resample_to_img +from nilearn.maskers import NiftiLabelsMasker, NiftiMasker +from numpy.testing import assert_array_almost_equal, assert_array_equal +from scipy.stats import trim_mean + +from junifer.datareader import DefaultDataReader +from junifer.markers.falff import AmplitudeLowFrequencyFluctuationParcels +from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber +from junifer.pipeline.utils import _check_afni + + +def test_AmplitudeLowFrequencyFluctuationParcels_python(): + """Test AmplitudeLowFrequencyFluctuationParcels using python.""" + # Get the SPM auditory data: + + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker = AmplitudeLowFrequencyFluctuationParcels( + parcellation="Schaefer100x7", method="mean", use_afni=False, + fractional=False) + jun_values4d = marker.fit_transform(input)["BOLD"] + + assert jun_values4d["data"].ndim == 2 + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +def test_AmplitudeLowFrequencyFluctuationParcels_afni(): + """Test AmplitudeLowFrequencyFluctuationParcels using afni.""" + # Get the SPM auditory data: + + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker = AmplitudeLowFrequencyFluctuationParcels( + parcellation="Schaefer100x7", method="mean", use_afni=True, + fractional=False) + jun_values4d = marker.fit_transform(input)["BOLD"] + + assert jun_values4d["data"].ndim == 2 + + # Again, should be blazing fast + jun_values4d = marker.fit_transform(input)["BOLD"] \ No newline at end of file -- 2.52.0 From 9f3ee356e896545a803a18c411ecea25ac7acb8d Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Thu, 15 Dec 2022 13:26:28 +0100 Subject: [PATCH 04/10] Some more fixes --- junifer/markers/falff/falff_base.py | 19 ++++++++++++++++++- junifer/markers/falff/falff_estimator.py | 6 +++--- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/junifer/markers/falff/falff_base.py b/junifer/markers/falff/falff_base.py index ba1583daa..32d11a9e4 100644 --- a/junifer/markers/falff/falff_base.py +++ b/junifer/markers/falff/falff_base.py @@ -4,7 +4,7 @@ # Kaustubh R. Patil # License: AGPL -from typing import TYPE_CHECKING, Dict, List, Optional, Tuple +from typing import Dict, List, Optional from ..base import BaseMarker @@ -176,6 +176,13 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): * ``row_names`` (if more than one row is present in data): "scan" """ + 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.") estimator = AmplitudeLowFrequencyFluctuationEstimator( use_afni=self.use_afni ) @@ -197,3 +204,13 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): out = self._postprocess(post_input) return out + + def _postprocess(self, input: Dict) -> Dict: + """Postprocess the output of the estimator. + + Parameters + ---------- + input : dict + The output of the estimator. It must have the following + """ + raise NotImplementedError("_postprocess must be implemented") \ No newline at end of file diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py index 5d5825c2e..82c1e255d 100644 --- a/junifer/markers/falff/falff_estimator.py +++ b/junifer/markers/falff/falff_estimator.py @@ -75,8 +75,8 @@ class AmplitudeLowFrequencyFluctuationEstimator: process = subprocess.run( cmd, stdin=subprocess.DEVNULL, - # stdout=subprocess.STDOUT, - # stderr=subprocess.STDOUT, + stdout=subprocess.STDOUT, + stderr=subprocess.STDOUT, shell=True, check=False, ) @@ -205,7 +205,7 @@ class AmplitudeLowFrequencyFluctuationEstimator: fALFF map. """ timeseries = data.get_fdata() - tr = tr or data.header["pixdim"][4] + tr = tr or data.header["pixdim"][4] # type: ignore # bandpass the data within the lowpass and highpass cutoff freqs Nq = 1 / (2 * tr) Wn = np.array([highpass / Nq, lowpass / Nq]) -- 2.52.0 From b3e6b1a0d632604615d385a0e01d23b1f34cb54c Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Thu, 15 Dec 2022 14:05:21 +0100 Subject: [PATCH 05/10] Add suffix to avoid having different parameters on the same file name --- junifer/markers/falff/falff_estimator.py | 8 ++-- junifer/markers/falff/falff_parcels.py | 3 +- junifer/markers/falff/falff_spheres.py | 55 ++++++++++++++---------- 3 files changed, 38 insertions(+), 28 deletions(-) diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py index 82c1e255d..92d7e29d4 100644 --- a/junifer/markers/falff/falff_estimator.py +++ b/junifer/markers/falff/falff_estimator.py @@ -109,8 +109,6 @@ class AmplitudeLowFrequencyFluctuationEstimator: Highpass cutoff frequency. lowpass : float Lowpass cutoff frequency. - order : int - Order of the filter. Not used by AFNI. tr : float, optional The Repetition Time of the BOLD data. @@ -146,7 +144,9 @@ class AmplitudeLowFrequencyFluctuationEstimator: bp_cmd += f"-dt {tr} " self._run_afni_cmd(bp_cmd) - alff_fname = self.temp_dir_path / "alff.nii" + params_suffix = f"_{highpass}_{lowpass}_{tr}" + + alff_fname = self.temp_dir_path / f"alff{params_suffix}.nii" convert_cmd = ( "3dAFNItoNIFTI " @@ -155,7 +155,7 @@ class AmplitudeLowFrequencyFluctuationEstimator: ) self._run_afni_cmd(convert_cmd) - falff_fname = self.temp_dir_path / "falff.nii" + falff_fname = self.temp_dir_path / f"falff{params_suffix}.nii" convert_cmd = ( "3dAFNItoNIFTI " diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py index d0f04d78d..789c0f063 100644 --- a/junifer/markers/falff/falff_parcels.py +++ b/junifer/markers/falff/falff_parcels.py @@ -1,10 +1,11 @@ """Provide class for computing fALFF on parcels.""" + # Authors: Federico Raimondo # Amir Omidvarnia # Kaustubh R. Patil # License: AGPL -from typing import Dict, List, Optional, Union, Tuple +from typing import Dict, List, Optional, Union from ...api.decorators import register_marker diff --git a/junifer/markers/falff/falff_spheres.py b/junifer/markers/falff/falff_spheres.py index 33b517429..64b8da3bf 100644 --- a/junifer/markers/falff/falff_spheres.py +++ b/junifer/markers/falff/falff_spheres.py @@ -1,10 +1,11 @@ """Provide class for computing fALFF on spheres.""" + # Authors: Federico Raimondo # Amir Omidvarnia # Kaustubh R. Patil # License: AGPL -from typing import Dict, Optional, Tuple +from typing import Dict, Optional from ...api.decorators import register_marker @@ -27,19 +28,30 @@ class AmplitudeLowFrequencyFluctuationSpheres( 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). + fractional : bool + Whether to compute fractional ALFF. highpass : float, optional The highpass cutoff frequency for the bandpass filter (default 0.01). lowpass : float The lowpass cutoff frequency for the bandpass filter (default 0.1). order : int The order of the bandpass filter (default 4). + tr : float, optional + The Repetition Time of the BOLD data. If None, will extract + the TR from NIFTI header (default None). + use_afni : bool, optional + Whether to use AFNI for computing. If None, will use AFNI only + if available (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). - tr : float, optional - The Repetition Time of the BOLD data. If None, will extract - the TR from NIFTI header (default None). + method : str, optional + The method to perform aggregation using. Check valid options in + :func:`junifer.stats.get_aggfunc_by_name` (default "mean"). + method_params : dict, optional + Parameters to pass to the aggregation function. Check valid options in + :func:`junifer.stats.get_aggfunc_by_name`. name : str, optional The name of the marker. If None, will use the class name (default None). @@ -62,27 +74,34 @@ class AmplitudeLowFrequencyFluctuationSpheres( def __init__( self, coords: str, + fractional: bool, + radius: Optional[float] = None, highpass: float = 0.01, lowpass: float = 0.1, order: int = 4, - radius: Optional[float] = None, - mask: Optional[str] = None, tr: Optional[float] = None, + use_afni: Optional[bool] = None, + mask: Optional[str] = None, + method: str = "mean", + method_params: Optional[Dict] = None, name: Optional[str] = None, ) -> None: self.coords = coords self.radius = radius self.mask = mask - + self.method = method + self.method_params = method_params super().__init__( + fractional=fractional, highpass=highpass, lowpass=lowpass, order=order, tr=tr, name=name, + use_afni=use_afni, ) - def compute(self, input: Dict) -> Tuple[Dict, Dict]: + def _postprocess(self, input: Dict) -> Dict: """Compute ALFF and fALFF. Parameters @@ -98,33 +117,23 @@ class AmplitudeLowFrequencyFluctuationSpheres( Returns ------- - alff: dict + dict The computed ALFF as dictionary. The dictionary has the following keys: - * ``data`` : the actual computed values as a numpy.ndarray - * ``columns`` : the column labels for the computed values as a list - - falff: dict - The computed fALFF as dictionary. The dictionary has the following - keys: - * ``data`` : the actual computed values as a numpy.ndarray * ``columns`` : the column labels for the computed values as a list """ pa = SphereAggregation( coords=self.coords, radius=self.radius, - method="mean", + method=self.method, + method_params=self.method_params, mask=self.mask, - on="BOLD", + on="fALFF", ) # get the 2D timeseries after parcel aggregation - spheres = pa.compute(input) - timeseries = spheres["data"] - labels = spheres["columns"] - tr = self.tr or input["BOLD"].header["pixdim"][4] + out = pa.compute(input) - out = super().compute_falff(timeseries, labels, tr) return out -- 2.52.0 From c10cf0d8099e716abe3efb9279b9b393bd97e1ce Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Thu, 15 Dec 2022 15:18:37 +0100 Subject: [PATCH 06/10] Fix mixin use_afni --- junifer/markers/falff/falff_estimator.py | 4 ++-- junifer/markers/falff/tests/test_falff_parcels.py | 10 +++++++++- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py index 92d7e29d4..e3ef7c815 100644 --- a/junifer/markers/falff/falff_estimator.py +++ b/junifer/markers/falff/falff_estimator.py @@ -75,8 +75,8 @@ class AmplitudeLowFrequencyFluctuationEstimator: process = subprocess.run( cmd, stdin=subprocess.DEVNULL, - stdout=subprocess.STDOUT, - stderr=subprocess.STDOUT, + # stdout=subprocess.STDOUT, + # stderr=subprocess.STDOUT, shell=True, check=False, ) diff --git a/junifer/markers/falff/tests/test_falff_parcels.py b/junifer/markers/falff/tests/test_falff_parcels.py index 887db8368..545d1527a 100644 --- a/junifer/markers/falff/tests/test_falff_parcels.py +++ b/junifer/markers/falff/tests/test_falff_parcels.py @@ -35,6 +35,7 @@ def test_AmplitudeLowFrequencyFluctuationParcels_python(): fractional=False) jun_values4d = marker.fit_transform(input)["BOLD"] + assert marker.use_afni is False assert jun_values4d["data"].ndim == 2 @@ -53,9 +54,16 @@ def test_AmplitudeLowFrequencyFluctuationParcels_afni(): marker = AmplitudeLowFrequencyFluctuationParcels( parcellation="Schaefer100x7", method="mean", use_afni=True, fractional=False) + assert marker.use_afni is True jun_values4d = marker.fit_transform(input)["BOLD"] assert jun_values4d["data"].ndim == 2 # Again, should be blazing fast - jun_values4d = marker.fit_transform(input)["BOLD"] \ No newline at end of file + marker = AmplitudeLowFrequencyFluctuationParcels( + parcellation="Schaefer100x7", method="mean", + fractional=False) + marker.validate(list(input.keys())) + assert marker.use_afni is True + jun_values4d = marker.fit_transform(input)["BOLD"] + assert jun_values4d["data"].ndim == 2 \ No newline at end of file -- 2.52.0 From ffe5291cb87c29e9a0b0c0ef0b24172666ea480a Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Fri, 16 Dec 2022 15:29:26 +0100 Subject: [PATCH 07/10] ALFF ready for review --- docs/builtin.rst | 12 +- docs/changes/latest.inc | 2 + junifer/markers/__init__.py | 5 +- junifer/markers/falff/__init__.py | 1 + junifer/markers/falff/falff_base.py | 64 +---- junifer/markers/falff/falff_estimator.py | 89 ++++--- junifer/markers/falff/falff_parcels.py | 4 - junifer/markers/falff/falff_spheres.py | 4 - .../falff/tests/test_falff_estimator.py | 241 ++++++++++++++++++ .../markers/falff/tests/test_falff_parcels.py | 137 ++++++++-- .../markers/falff/tests/test_falff_spheres.py | 162 ++++++++++++ 11 files changed, 594 insertions(+), 127 deletions(-) create mode 100644 junifer/markers/falff/tests/test_falff_estimator.py create mode 100644 junifer/markers/falff/tests/test_falff_spheres.py diff --git a/docs/builtin.rst b/docs/builtin.rst index df91e9924..a814a751d 100644 --- a/docs/builtin.rst +++ b/docs/builtin.rst @@ -174,7 +174,14 @@ Available - Calculate regional homogeneity over spheres placed on coordinates - Done - 0.0.1 - + * - :class:`junifer.markers.AmplitudeLowFrequencyFluctuationParcels` + - Calculate (f)ALFF and aggregate using parcellations + - Done + - 0.0.1 + * - :class:`junifer.markers.AmplitudeLowFrequencyFluctuationSpheres` + - Calculate (f)ALFF and aggregate using spheres placed on coordinates + - Done + - 0.0.1 Planned ~~~~~~~ @@ -189,9 +196,6 @@ Planned * - Connectedness - Compute connectedness - :gh:`34` - * - ALFF and (f)ALFF - - Detect amplitude of low-frequency fluctuation (ALFF) for resting-state fMRI - - :gh:`35` * - 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 196986ab5..f8c3d88ed 100644 --- a/docs/changes/latest.inc +++ b/docs/changes/latest.inc @@ -98,6 +98,8 @@ Enhancements - Implement :class:`junifer.markers.ReHoParcels` and :class:`junifer.markers.ReHoSpheres` markers (:gh:`36` by `Synchon Mandal`_). +- Implement :class:`junifer.markers.AmplitudeLowFrequencyFluctuationParcels` and :class:`junifer.markers.AmplitudeLowFrequencyFluctuationSpheres` markers (:gh:`35` by `Fede Raimondo`_). + Bugs ~~~~ diff --git a/junifer/markers/__init__.py b/junifer/markers/__init__.py index 2af03f4ea..bc069e3db 100644 --- a/junifer/markers/__init__.py +++ b/junifer/markers/__init__.py @@ -13,4 +13,7 @@ from .functional_connectivity_spheres import FunctionalConnectivitySpheres from .parcel_aggregation import ParcelAggregation from .sphere_aggregation import SphereAggregation from .reho import ReHoParcels, ReHoSpheres -from .falff import AmplitudeLowFrequencyFluctuationParcels +from .falff import ( + AmplitudeLowFrequencyFluctuationParcels, + AmplitudeLowFrequencyFluctuationSpheres, +) diff --git a/junifer/markers/falff/__init__.py b/junifer/markers/falff/__init__.py index dcbb4d18c..6e1023665 100644 --- a/junifer/markers/falff/__init__.py +++ b/junifer/markers/falff/__init__.py @@ -4,3 +4,4 @@ # License: AGPL from .falff_parcels import AmplitudeLowFrequencyFluctuationParcels +from .falff_spheres import AmplitudeLowFrequencyFluctuationSpheres diff --git a/junifer/markers/falff/falff_base.py b/junifer/markers/falff/falff_base.py index 32d11a9e4..040408308 100644 --- a/junifer/markers/falff/falff_base.py +++ b/junifer/markers/falff/falff_base.py @@ -9,7 +9,7 @@ from typing import Dict, List, Optional from ..base import BaseMarker from .falff_estimator import AmplitudeLowFrequencyFluctuationEstimator -from ...utils.logging import raise_error, warn_with_log +from ...utils.logging import raise_error class AmplitudeLowFrequencyFluctuationBase(BaseMarker): @@ -23,8 +23,6 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): Highpass cutoff frequency. lowpass : float Lowpass cutoff frequency. - order : int - Order of the filter. tr : float, optional The Repetition Time of the BOLD data. If None, will extract the TR from NIFTI header (default None). @@ -57,7 +55,6 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): fractional: bool, highpass: float, lowpass: float, - order: int, tr: Optional[float] = None, use_afni: Optional[bool] = None, name: Optional[str] = None, @@ -70,9 +67,6 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): raise_error("Highpass must be lower than lowpass") self.highpass = highpass self.lowpass = lowpass - if order <= 0 and use_afni is False: - raise_error("Order must be positive") - self.order = order self.tr = tr self.use_afni = use_afni self.fractional = fractional @@ -83,41 +77,6 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): name = f"{self.__class__.__name__}{suffix}" super().__init__(on="BOLD", name=name) - def validate(self, input: List[str]) -> List[str]: - """Validate the the pipeline step. - - Parameters - ---------- - input : list of str - The input to the pipeline step. - - Returns - ------- - list of str - The output of the pipeline step. - - Raises - ------ - ValueError - If the pipeline step object is missing dependencies required for - its working, if the input does not have the required data, or - if AFNI was not used and the order is not positive. - - Warns - ----- - UserWarning - If AFNI is used and the order is not 0. - """ - out = super().validate(input) - if self.use_afni is True and self.order > 0: - warn_with_log( - "AFNI will not consider the order of the filter. Set this " - "parameter to 0 to avoid this warning.") - elif self.use_afni is False and self.order <= 0: - raise_error( - "Order must be positive if AFNI is not used.") - return out - def get_valid_inputs(self) -> List[str]: """Get valid data types for input. @@ -130,21 +89,21 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): valid = ["BOLD"] return valid - def get_output_type(self, input: List[str]) -> List[str]: + def get_output_type(self, input_type: str) -> str: """Get output type. Parameters ---------- - input : list of str - The type of data to work on. + input_type : str + The data type input to the marker. Returns ------- - list of str - The list of storage types. + str + The storage type output by the marker. """ - return ["table"] + return "table" def compute( self, @@ -183,15 +142,14 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): "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.") - estimator = AmplitudeLowFrequencyFluctuationEstimator( - use_afni=self.use_afni - ) + + estimator = AmplitudeLowFrequencyFluctuationEstimator() alff, falff = estimator.fit_transform( + use_afni=self.use_afni, input_data=input, highpass=self.highpass, lowpass=self.lowpass, - order=self.order, tr=self.tr, ) post_data = falff if self.fractional else alff @@ -213,4 +171,4 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): input : dict The output of the estimator. It must have the following """ - raise NotImplementedError("_postprocess must be implemented") \ No newline at end of file + raise NotImplementedError("_postprocess must be implemented") diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py index e3ef7c815..d18fba09f 100644 --- a/junifer/markers/falff/falff_estimator.py +++ b/junifer/markers/falff/falff_estimator.py @@ -4,6 +4,7 @@ # Federico Raimondo # License: AGPL +import typing from typing import TYPE_CHECKING, Any, Dict, Tuple, Union, Optional import shutil @@ -14,7 +15,7 @@ from pathlib import Path import nibabel as nib import numpy as np -from scipy import signal +from scipy.fft import fft, fftfreq from nilearn import image as nimg @@ -45,8 +46,7 @@ class AmplitudeLowFrequencyFluctuationEstimator: """ - def __init__(self, use_afni: bool) -> None: - self.use_afni = use_afni + def __init__(self) -> None: self._file_path = None # Create temporary directory for intermittent storage of assets during # computation via afni's 3dReHo @@ -54,6 +54,7 @@ class AmplitudeLowFrequencyFluctuationEstimator: def __del__(self) -> None: """Cleanup.""" + print("Cleaning up temporary directory...") # Delete temporary directory and ignore errors for read-only files shutil.rmtree(self.temp_dir_path, ignore_errors=True) @@ -72,6 +73,7 @@ class AmplitudeLowFrequencyFluctuationEstimator: If AFNI command fails. """ logger.info(f"AFNI command to be executed: {cmd}") + # TODO: Figure out how to capture stdout and stderr process = subprocess.run( cmd, stdin=subprocess.DEVNULL, @@ -130,7 +132,11 @@ class AmplitudeLowFrequencyFluctuationEstimator: nifti_in_file_path = self.temp_dir_path / "input.nii" nib.save(data, nifti_in_file_path) - # Bandpass the data + params_suffix = f"_{highpass}_{lowpass}_{tr}" + alff_fname = self.temp_dir_path / f"alff{params_suffix}.nii" + falff_fname = self.temp_dir_path / f"falff{params_suffix}.nii" + + # Use afni's 3dRSFC to compute ALFF and fALFF falff_afni_out_path_prefix = self.temp_dir_path / "temp_falff" bp_cmd = ( @@ -138,16 +144,13 @@ class AmplitudeLowFrequencyFluctuationEstimator: f"-prefix {falff_afni_out_path_prefix.resolve()} " f"-input {nifti_in_file_path.resolve()} " f"-band {highpass} {lowpass} " - "-no_rsfa " + "-no_rsfa -nosat -nodetrend " ) if tr is not None: bp_cmd += f"-dt {tr} " self._run_afni_cmd(bp_cmd) - params_suffix = f"_{highpass}_{lowpass}_{tr}" - - alff_fname = self.temp_dir_path / f"alff{params_suffix}.nii" - + # Convert afni's output to nifti convert_cmd = ( "3dAFNItoNIFTI " f"-prefix {alff_fname.resolve()} " @@ -155,8 +158,6 @@ class AmplitudeLowFrequencyFluctuationEstimator: ) self._run_afni_cmd(convert_cmd) - falff_fname = self.temp_dir_path / f"falff{params_suffix}.nii" - convert_cmd = ( "3dAFNItoNIFTI " f"-prefix {falff_fname.resolve()} " @@ -165,7 +166,7 @@ class AmplitudeLowFrequencyFluctuationEstimator: self._run_afni_cmd(convert_cmd) # Cleanup intermediate files - for fname in self.temp_dir_path.glob("temp_falff*"): + for fname in self.temp_dir_path.glob("temp_*"): fname.unlink() # Load niftis @@ -179,7 +180,6 @@ class AmplitudeLowFrequencyFluctuationEstimator: data: Union["Nifti1Image", "Nifti2Image"], highpass: float, lowpass: float, - order: int, tr: Optional[float], ) -> Tuple["Nifti1Image", "Nifti1Image"]: """Compute (f)ALFF map. @@ -192,8 +192,6 @@ class AmplitudeLowFrequencyFluctuationEstimator: Highpass cutoff frequency. lowpass : float Lowpass cutoff frequency. - order : int - Order of the filter. tr : float, optional The Repetition Time of the BOLD data. @@ -204,44 +202,65 @@ class AmplitudeLowFrequencyFluctuationEstimator: falff: Niimg-like object fALFF map. """ - timeseries = data.get_fdata() - tr = tr or data.header["pixdim"][4] # type: ignore + timeseries = data.get_fdata().copy() + if tr is None: + tr = float(data.header["pixdim"][4]) # type: ignore + logger.info(f"TR Not provided, using TR from header = {tr}") # bandpass the data within the lowpass and highpass cutoff freqs - Nq = 1 / (2 * tr) - Wn = np.array([highpass / Nq, lowpass / Nq]) - b, a = signal.butter(N=order, Wn=Wn, btype="bandpass") - ts_filt = signal.filtfilt(b, a, timeseries, axis=0) + ts_fft = fft(timeseries, axis=-1) + ts_fft = typing.cast(np.ndarray, ts_fft) + fft_freqs = np.abs(fftfreq(timeseries.shape[-1], tr)) - ALFF = np.std(ts_filt, axis=0) - PSD_tot = np.std(timeseries, axis=0) + dFreq = fft_freqs[1] - fft_freqs[0] + nyquist = np.max(fft_freqs) + nfft = len(fft_freqs) + logger.info( + f"FFT: nfft = {nfft}, dFreq = {dFreq}, nyquist = {nyquist}" + ) - fALFF = np.divide(ALFF, PSD_tot) - alff_img = nimg.new_img_like(data, ALFF) - falff_img = nimg.new_img_like(data, fALFF) + # First compute the denominator on the broadband signal + all_freq_mask = fft_freqs > 0 + denominator = np.sum(np.abs(ts_fft[..., all_freq_mask]), axis=-1) + + # Compute the numerator on the bandpassed signal + freq_mask = np.logical_and(fft_freqs > highpass, fft_freqs < lowpass) + # Compute ALFF + numerator = np.sum(np.abs(ts_fft[..., freq_mask]), axis=-1) + + # Compute fALFF, but avoid division by zero + denom_mask = denominator <= 0.000001 + denominator[denom_mask] = 1 # set to 1 to avoid division by zero + python_falff = np.divide(numerator, denominator) + # Set the values where denominator is zero to zero + python_falff[denom_mask] = 0 + + python_alff = numerator / np.sqrt(timeseries.shape[-1]) + alff_img = nimg.new_img_like(data, python_alff) + falff_img = nimg.new_img_like(data, python_falff) return alff_img, falff_img @lru_cache(maxsize=None, typed=True) # noqa: B019 def _compute( self, + use_afni: bool, data: Union["Nifti1Image", "Nifti2Image"], highpass: float, lowpass: float, - order: int, tr: Optional[float], ) -> Tuple["Nifti1Image", "Nifti1Image"]: """Compute the ALFF map with memorization. Parameters ---------- + use_afni : bool + Whether to use AFNI for computing. data : 4D Niimg-like object Images to process. highpass : float Highpass cutoff frequency. lowpass : float Lowpass cutoff frequency. - order : int - Order of the filter. tr : float, optional The Repetition Time of the BOLD data. @@ -252,7 +271,7 @@ class AmplitudeLowFrequencyFluctuationEstimator: falff: Niimg-like object fALFF map. """ - if self.use_afni: + if use_afni: output = self._compute_alff_afni( data=data, highpass=highpass, @@ -261,30 +280,30 @@ class AmplitudeLowFrequencyFluctuationEstimator: ) else: output = self._compute_alff_python( - data, highpass=highpass, lowpass=lowpass, order=order, tr=tr + data, highpass=highpass, lowpass=lowpass, tr=tr ) return output def fit_transform( self, + use_afni: bool, input_data: Dict[str, Any], highpass: float, lowpass: float, - order: int, tr: Optional[float], ) -> Tuple["Nifti1Image", "Nifti1Image"]: """Fit and transform for the estimator. Parameters ---------- + use_afni : bool + Whether to use AFNI for computing. input_data : dict The BOLD data as dictionary. highpass : float Highpass cutoff frequency. lowpass : float Lowpass cutoff frequency. - order : int - Order of the filter. tr : float, optional The Repetition Time of the BOLD data. @@ -311,9 +330,9 @@ class AmplitudeLowFrequencyFluctuationEstimator: logger.info(f"Using fALFF map cache at {self._file_path}.") # Compute return self._compute( + use_afni=use_afni, data=bold_data, highpass=highpass, lowpass=lowpass, - order=order, tr=tr, ) diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py index 789c0f063..eec5a4d8f 100644 --- a/junifer/markers/falff/falff_parcels.py +++ b/junifer/markers/falff/falff_parcels.py @@ -30,8 +30,6 @@ class AmplitudeLowFrequencyFluctuationParcels( The highpass cutoff frequency for the bandpass filter (default 0.01). lowpass : float The lowpass cutoff frequency for the bandpass filter (default 0.1). - order : int - The order of the bandpass filter (default 4). tr : float, optional The Repetition Time of the BOLD data. If None, will extract the TR from NIFTI header (default None). @@ -73,7 +71,6 @@ class AmplitudeLowFrequencyFluctuationParcels( fractional: bool, highpass: float = 0.01, lowpass: float = 0.1, - order: int = 4, tr: Optional[float] = None, use_afni: Optional[bool] = None, mask: Optional[str] = None, @@ -89,7 +86,6 @@ class AmplitudeLowFrequencyFluctuationParcels( fractional=fractional, highpass=highpass, lowpass=lowpass, - order=order, tr=tr, name=name, use_afni=use_afni, diff --git a/junifer/markers/falff/falff_spheres.py b/junifer/markers/falff/falff_spheres.py index 64b8da3bf..d25f43ffc 100644 --- a/junifer/markers/falff/falff_spheres.py +++ b/junifer/markers/falff/falff_spheres.py @@ -34,8 +34,6 @@ class AmplitudeLowFrequencyFluctuationSpheres( The highpass cutoff frequency for the bandpass filter (default 0.01). lowpass : float The lowpass cutoff frequency for the bandpass filter (default 0.1). - order : int - The order of the bandpass filter (default 4). tr : float, optional The Repetition Time of the BOLD data. If None, will extract the TR from NIFTI header (default None). @@ -78,7 +76,6 @@ class AmplitudeLowFrequencyFluctuationSpheres( radius: Optional[float] = None, highpass: float = 0.01, lowpass: float = 0.1, - order: int = 4, tr: Optional[float] = None, use_afni: Optional[bool] = None, mask: Optional[str] = None, @@ -95,7 +92,6 @@ class AmplitudeLowFrequencyFluctuationSpheres( fractional=fractional, highpass=highpass, lowpass=lowpass, - order=order, tr=tr, name=name, use_afni=use_afni, diff --git a/junifer/markers/falff/tests/test_falff_estimator.py b/junifer/markers/falff/tests/test_falff_estimator.py new file mode 100644 index 000000000..64c9d4360 --- /dev/null +++ b/junifer/markers/falff/tests/test_falff_estimator.py @@ -0,0 +1,241 @@ +"""Provide test for (f)ALFF estimator.""" + +# Authors: Federico Raimondo +# License: AGPL + +import pytest +import time +from scipy.stats import pearsonr +from nibabel import Nifti1Image + +from junifer.datareader import DefaultDataReader +from junifer.markers.falff.falff_estimator import ( + AmplitudeLowFrequencyFluctuationEstimator, +) +from junifer.testing.datagrabbers import ( + PartlyCloudyTestingDataGrabber, +) +from junifer.pipeline.utils import _check_afni +from junifer.utils import logger + + +def test_AmplitudeLowFrequencyFluctuationEstimator_cache_python() -> None: + """Test that the cache works properly when using python.""" + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + + estimator = AmplitudeLowFrequencyFluctuationEstimator() + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=False, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + first_time = time.time() - start_time + logger.info(f"ALFF Estimator First time: {first_time}") + assert isinstance(alff, Nifti1Image) + assert isinstance(falff, Nifti1Image) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + # Now fit again, should be faster + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=False, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + second_time = time.time() - start_time + logger.info(f"ALFF Estimator Second time: {second_time}") + assert second_time < (first_time / 1000) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + # Now change a parameter, should compute again, without clearing the + # cache + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=False, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.11, + tr=None, + ) + third_time = time.time() - start_time + logger.info(f"ALFF Estimator Third time: {third_time}") + assert third_time > (first_time / 10) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 0 # no files in python + + # Now fit again with the previous params, should be fast + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=False, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + fourth = time.time() - start_time + logger.info(f"ALFF Estimator Fourth time: {fourth}") + assert fourth < (first_time / 1000) + n_files = len([x for x in 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: + input = dg["sub-02"] + + input = DefaultDataReader().fit_transform(input) + + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=False, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + fifth = time.time() - start_time + logger.info(f"ALFF Estimator Fifth time: {fifth}") + assert fifth > (first_time / 10) + n_files = len([x for x in 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_AmplitudeLowFrequencyFluctuationEstimator_cache_afni() -> None: + """Test that the cache works properly when using afni.""" + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + + estimator = AmplitudeLowFrequencyFluctuationEstimator() + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=True, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + first_time = time.time() - start_time + logger.info(f"ALFF Estimator First time: {first_time}") + assert isinstance(alff, Nifti1Image) + assert isinstance(falff, Nifti1Image) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 3 # input + alff + falff + + # Now fit again, should be faster + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=True, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + second_time = time.time() - start_time + logger.info(f"ALFF Estimator Second time: {second_time}") + assert second_time < (first_time / 1000) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 3 # input + alff + falff + + # Now change a parameter, should compute again, without clearing the + # cache + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=True, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.11, + tr=None, + ) + third_time = time.time() - start_time + logger.info(f"ALFF Estimator Third time: {third_time}") + assert third_time > (first_time / 10) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 5 # input + 2 * alff + 2 * falff + + # Now fit again with the previous params, should be fast + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=True, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + fourth = time.time() - start_time + logger.info(f"ALFF Estimator Fourth time: {fourth}") + assert fourth < (first_time / 1000) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 5 # input + 2 * alff + 2 * falff + + # Now change the data, it should clear the cache + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-02"] + + input = DefaultDataReader().fit_transform(input) + + start_time = time.time() + alff, falff = estimator.fit_transform( + use_afni=True, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=None, + ) + fifth = time.time() - start_time + logger.info(f"ALFF Estimator Fifth time: {fifth}") + assert fifth > (first_time / 10) + n_files = len([x for x in estimator.temp_dir_path.glob("*")]) + assert n_files == 3 # input + alff + falff + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +def test_AmplitudeLowFrequencyFluctuationEstimator_afni_vs_python() -> None: + """Test that the cache works properly when using afni.""" + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + estimator = AmplitudeLowFrequencyFluctuationEstimator() + + # Use an arbitrary TR to test the AFNI vs Python implementation + afni_alff, afni_falff = estimator.fit_transform( + use_afni=True, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=2.5, + ) + + python_alff, python_falff = estimator.fit_transform( + use_afni=False, + input_data=input["BOLD"], + highpass=0.01, + lowpass=0.1, + tr=2.5, + ) + + r, _ = pearsonr( + afni_alff.get_fdata().flatten(), python_alff.get_fdata().flatten() + ) + assert r > 0.99 + + r, _ = pearsonr( + afni_falff.get_fdata().flatten(), python_falff.get_fdata().flatten() + ) + assert r > 0.99 diff --git a/junifer/markers/falff/tests/test_falff_parcels.py b/junifer/markers/falff/tests/test_falff_parcels.py index 545d1527a..a4cb63fcc 100644 --- a/junifer/markers/falff/tests/test_falff_parcels.py +++ b/junifer/markers/falff/tests/test_falff_parcels.py @@ -1,27 +1,24 @@ -"""Provide test for parcel-aggregated (f)ALFF .""" +"""Provide test for parcel-aggregated (f)ALFF.""" # Authors: Federico Raimondo # Synchon Mandal # License: AGPL -from pathlib import Path - -import nibabel as nib -import numpy as np import pytest -from nilearn import datasets -from nilearn.image import concat_imgs, math_img, new_img_like, resample_to_img -from nilearn.maskers import NiftiLabelsMasker, NiftiMasker -from numpy.testing import assert_array_almost_equal, assert_array_equal -from scipy.stats import trim_mean + +from pathlib import Path +from numpy.testing import assert_array_equal +from scipy.stats import pearsonr from junifer.datareader import DefaultDataReader from junifer.markers.falff import AmplitudeLowFrequencyFluctuationParcels from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber from junifer.pipeline.utils import _check_afni +from junifer.storage import SQLiteFeatureStorage +from junifer.utils import logger -def test_AmplitudeLowFrequencyFluctuationParcels_python(): +def test_AmplitudeLowFrequencyFluctuationParcels_python() -> None: """Test AmplitudeLowFrequencyFluctuationParcels using python.""" # Get the SPM auditory data: @@ -31,39 +28,127 @@ def test_AmplitudeLowFrequencyFluctuationParcels_python(): input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", method="mean", use_afni=False, - fractional=False) - jun_values4d = marker.fit_transform(input)["BOLD"] + parcellation="Schaefer100x7", + method="mean", + use_afni=False, + fractional=False, + ) + python_values = marker.fit_transform(input)["BOLD"]["data"] assert marker.use_afni is False - assert jun_values4d["data"].ndim == 2 + assert python_values.ndim == 2 + assert python_values.shape == (1, 100) @pytest.mark.skipif( _check_afni() is False, reason="requires afni to be in PATH" ) -def test_AmplitudeLowFrequencyFluctuationParcels_afni(): +def test_AmplitudeLowFrequencyFluctuationParcels_afni() -> None: """Test AmplitudeLowFrequencyFluctuationParcels using afni.""" # Get the SPM auditory data: - with PartlyCloudyTestingDataGrabber() as dg: input = dg["sub-01"] input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", method="mean", use_afni=True, - fractional=False) + parcellation="Schaefer100x7", + method="mean", + use_afni=True, + fractional=False, + ) assert marker.use_afni is True - jun_values4d = marker.fit_transform(input)["BOLD"] + afni_values = marker.fit_transform(input)["BOLD"]["data"] - assert jun_values4d["data"].ndim == 2 + assert afni_values.ndim == 2 + assert afni_values.shape == (1, 100) # Again, should be blazing fast marker = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", method="mean", - fractional=False) - marker.validate(list(input.keys())) + parcellation="Schaefer100x7", method="mean", fractional=False + ) + assert marker.use_afni is None + afni_values2 = marker.fit_transform(input)["BOLD"]["data"] assert marker.use_afni is True - jun_values4d = marker.fit_transform(input)["BOLD"] - assert jun_values4d["data"].ndim == 2 \ No newline at end of file + assert_array_equal(afni_values, afni_values2) + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +@pytest.mark.parametrize( + "fractional", [True, False], ids=["fractional", "non-fractional"] +) +def test_AmplitudeLowFrequencyFluctuationParcels_python_vs_afni( + fractional: bool, +) -> None: + """Test AmplitudeLowFrequencyFluctuationParcels using python. + + Parameters + ---------- + factional : bool + Whether to compute fractional ALFF or not. + """ + + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker_python = AmplitudeLowFrequencyFluctuationParcels( + parcellation="Schaefer100x7", + method="mean", + use_afni=False, + fractional=fractional, + ) + python_values = marker_python.fit_transform(input)["BOLD"]["data"] + + assert marker_python.use_afni is False + assert python_values.ndim == 2 + assert python_values.shape == (1, 100) + + marker_afni = AmplitudeLowFrequencyFluctuationParcels( + parcellation="Schaefer100x7", + method="mean", + use_afni=True, + fractional=fractional, + ) + afni_values = marker_afni.fit_transform(input)["BOLD"]["data"] + + assert marker_afni.use_afni is True + assert afni_values.ndim == 2 + assert afni_values.shape == (1, 100) + + r, p = pearsonr(python_values[0], afni_values[0]) + logger.info(f"Correlation between python and afni: {r} (p={p})") + assert r > 0.99 + + +def test_AmplitudeLowFrequencyFluctuationParcels_storage( + tmp_path: Path, +) -> None: + """Test AmplitudeLowFrequencyFluctuationParcels storage. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + """ + with PartlyCloudyTestingDataGrabber() as dg: + # Use first subject + input = dg["sub-01"] + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker = AmplitudeLowFrequencyFluctuationParcels( + parcellation="Schaefer100x7", + method="mean", + use_afni=False, + fractional=True, + ) + storage = SQLiteFeatureStorage(tmp_path / "alff_parcels.sqlite") + + # Fit transform marker on data with storage + marker.fit_transform( + input=input, + storage=storage, + ) diff --git a/junifer/markers/falff/tests/test_falff_spheres.py b/junifer/markers/falff/tests/test_falff_spheres.py new file mode 100644 index 000000000..a8387367a --- /dev/null +++ b/junifer/markers/falff/tests/test_falff_spheres.py @@ -0,0 +1,162 @@ +"""Provide test for parcel-aggregated (f)ALFF.""" + +# Authors: Federico Raimondo +# Synchon Mandal +# License: AGPL + +import pytest + +from pathlib import Path + +from numpy.testing import assert_array_equal +from scipy.stats import pearsonr + +from junifer.datareader import DefaultDataReader +from junifer.markers.falff import AmplitudeLowFrequencyFluctuationSpheres +from junifer.testing.datagrabbers import PartlyCloudyTestingDataGrabber +from junifer.pipeline.utils import _check_afni +from junifer.storage import SQLiteFeatureStorage +from junifer.utils import logger + + +def test_AmplitudeLowFrequencyFluctuationSpheres_python() -> None: + """Test AmplitudeLowFrequencyFluctuationSpheres using python.""" + # Get the SPM auditory data: + + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker = AmplitudeLowFrequencyFluctuationSpheres( + coords="DMNBuckner", + radius=5, + method="mean", + use_afni=False, + fractional=False, + ) + python_values = marker.fit_transform(input)["BOLD"]["data"] + + assert marker.use_afni is False + assert python_values.ndim == 2 + assert python_values.shape == (1, 6) + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +def test_AmplitudeLowFrequencyFluctuationSpheres_afni() -> None: + """Test AmplitudeLowFrequencyFluctuationSpheres using afni.""" + # Get the SPM auditory data: + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker = AmplitudeLowFrequencyFluctuationSpheres( + coords="DMNBuckner", + radius=5, + method="mean", + use_afni=True, + fractional=False, + ) + assert marker.use_afni is True + afni_values = marker.fit_transform(input)["BOLD"]["data"] + + assert afni_values.ndim == 2 + assert afni_values.shape == (1, 6) + + # Again, should be blazing fast + marker = AmplitudeLowFrequencyFluctuationSpheres( + coords="DMNBuckner", + radius=5, + method="mean", + fractional=False, + ) + assert marker.use_afni is None + afni_values2 = marker.fit_transform(input)["BOLD"]["data"] + assert marker.use_afni is True + assert_array_equal(afni_values, afni_values2) + + +@pytest.mark.skipif( + _check_afni() is False, reason="requires afni to be in PATH" +) +@pytest.mark.parametrize( + "fractional", [True, False], ids=["fractional", "non-fractional"] +) +def test_AmplitudeLowFrequencyFluctuationSpheres_python_vs_afni( + fractional: bool, +) -> None: + """Test AmplitudeLowFrequencyFluctuationSpheres using python. + + Parameters + ---------- + fractional : bool + Whether to compute fractional ALFF or not. + """ + with PartlyCloudyTestingDataGrabber() as dg: + input = dg["sub-01"] + + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker_python = AmplitudeLowFrequencyFluctuationSpheres( + coords="DMNBuckner", + radius=5, + method="mean", + use_afni=False, + fractional=fractional, + ) + python_values = marker_python.fit_transform(input)["BOLD"]["data"] + + assert marker_python.use_afni is False + assert python_values.ndim == 2 + assert python_values.shape == (1, 6) + + marker_afni = AmplitudeLowFrequencyFluctuationSpheres( + coords="DMNBuckner", + radius=5, + method="mean", + use_afni=True, + fractional=fractional, + ) + afni_values = marker_afni.fit_transform(input)["BOLD"]["data"] + + assert marker_afni.use_afni is True + assert afni_values.ndim == 2 + assert afni_values.shape == (1, 6) + + r, p = pearsonr(python_values[0], afni_values[0]) + logger.info(f"Correlation between python and afni: {r} (p={p})") + assert r > 0.99 + + +def test_AmplitudeLowFrequencyFluctuationSpheres_storage( + tmp_path: Path, +) -> None: + """Test AmplitudeLowFrequencyFluctuationSpheres storage. + + Parameters + ---------- + tmp_path : pathlib.Path + The path to the test directory. + """ + with PartlyCloudyTestingDataGrabber() as dg: + # Use first subject + input = dg["sub-01"] + input = DefaultDataReader().fit_transform(input) + # Create ParcelAggregation object + marker = AmplitudeLowFrequencyFluctuationSpheres( + coords="DMNBuckner", + radius=5, + method="mean", + use_afni=False, + fractional=True, + ) + storage = SQLiteFeatureStorage(tmp_path / "alff_parcels.sqlite") + + # Fit transform marker on data with storage + marker.fit_transform( + input=input, + storage=storage, + ) -- 2.52.0 From 666779bd6b2855622a505c063a870db8a54344ad Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Mon, 19 Dec 2022 13:36:55 +0100 Subject: [PATCH 08/10] Rebase + address comments --- junifer/markers/falff/falff_base.py | 24 ++++++++++------- junifer/markers/falff/falff_estimator.py | 26 +++++++++---------- junifer/markers/falff/falff_parcels.py | 10 +++---- junifer/markers/falff/falff_spheres.py | 11 ++++---- .../markers/falff/tests/test_falff_parcels.py | 15 ++++++----- .../markers/falff/tests/test_falff_spheres.py | 17 +++++++----- 6 files changed, 56 insertions(+), 47 deletions(-) diff --git a/junifer/markers/falff/falff_base.py b/junifer/markers/falff/falff_base.py index 040408308..dcf028bc8 100644 --- a/junifer/markers/falff/falff_base.py +++ b/junifer/markers/falff/falff_base.py @@ -1,4 +1,5 @@ """Provide abstract class for computing fALFF.""" + # Authors: Federico Raimondo # Amir Omidvarnia # Kaustubh R. Patil @@ -6,6 +7,7 @@ from typing import Dict, List, Optional +from abc import abstractmethod from ..base import BaseMarker from .falff_estimator import AmplitudeLowFrequencyFluctuationEstimator @@ -19,11 +21,11 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): ---------- fractional : bool Whether to compute fractional ALFF. - highpass : float + highpass : positive float Highpass cutoff frequency. - lowpass : float + lowpass : positive float Lowpass cutoff frequency. - tr : float, optional + tr : positive float, optional The Repetition Time of the BOLD data. If None, will extract the TR from NIFTI header (default None). use_afni : bool, optional @@ -32,6 +34,7 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): name : str, optional The name of the marker. If None, it will use the class name (default None). + Notes ----- The `tr` parameter is crucial for the correctness of fALFF/ALFF @@ -86,8 +89,7 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): The list of data types that can be used as input for this marker. """ - valid = ["BOLD"] - return valid + return ["BOLD"] def get_output_type(self, input_type: str) -> str: """Get output type. @@ -141,7 +143,8 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): "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.") + "before calling the `compute` method." + ) estimator = AmplitudeLowFrequencyFluctuationEstimator() @@ -155,14 +158,15 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): post_data = falff if self.fractional else alff post_input = { - 'data': post_data, - 'path': None, + "data": post_data, + "path": None, } out = self._postprocess(post_input) return out + @abstractmethod def _postprocess(self, input: Dict) -> Dict: """Postprocess the output of the estimator. @@ -171,4 +175,6 @@ class AmplitudeLowFrequencyFluctuationBase(BaseMarker): input : dict The output of the estimator. It must have the following """ - raise NotImplementedError("_postprocess must be implemented") + raise_error( + "_postprocess must be implemented", klass=NotImplementedError + ) diff --git a/junifer/markers/falff/falff_estimator.py b/junifer/markers/falff/falff_estimator.py index d18fba09f..49d0f39d5 100644 --- a/junifer/markers/falff/falff_estimator.py +++ b/junifer/markers/falff/falff_estimator.py @@ -1,4 +1,4 @@ -"""Provide estimator class for regional homogeneity (ReHo).""" +"""Provide estimator class for (f)ALFF.""" # Authors: Synchon Mandal # Federico Raimondo @@ -107,11 +107,11 @@ class AmplitudeLowFrequencyFluctuationEstimator: ---------- data : 4D Niimg-like object Images to process. - highpass : float + highpass : positive float Highpass cutoff frequency. - lowpass : float + lowpass : positive float Lowpass cutoff frequency. - tr : float, optional + tr : positive float, optional The Repetition Time of the BOLD data. Returns @@ -188,11 +188,11 @@ class AmplitudeLowFrequencyFluctuationEstimator: ---------- data : 4D Niimg-like object Images to process. - highpass : float + highpass : positive float Highpass cutoff frequency. - lowpass : float + lowpass : positive float Lowpass cutoff frequency. - tr : float, optional + tr : positive float, optional The Repetition Time of the BOLD data. Returns @@ -257,11 +257,11 @@ class AmplitudeLowFrequencyFluctuationEstimator: Whether to use AFNI for computing. data : 4D Niimg-like object Images to process. - highpass : float + highpass : positive float Highpass cutoff frequency. - lowpass : float + lowpass : positive float Lowpass cutoff frequency. - tr : float, optional + tr : positive float, optional The Repetition Time of the BOLD data. Returns @@ -300,11 +300,11 @@ class AmplitudeLowFrequencyFluctuationEstimator: Whether to use AFNI for computing. input_data : dict The BOLD data as dictionary. - highpass : float + highpass : positive float Highpass cutoff frequency. - lowpass : float + lowpass : positive float Lowpass cutoff frequency. - tr : float, optional + tr : positive float, optional The Repetition Time of the BOLD data. Returns diff --git a/junifer/markers/falff/falff_parcels.py b/junifer/markers/falff/falff_parcels.py index eec5a4d8f..0d68ca34b 100644 --- a/junifer/markers/falff/falff_parcels.py +++ b/junifer/markers/falff/falff_parcels.py @@ -26,11 +26,11 @@ class AmplitudeLowFrequencyFluctuationParcels( :func:`junifer.data.parcellations.list_parcellations`. fractional : bool Whether to compute fractional ALFF. - highpass : float, optional + highpass : positive float, optional The highpass cutoff frequency for the bandpass filter (default 0.01). - lowpass : float + lowpass : positive float, optional The lowpass cutoff frequency for the bandpass filter (default 0.1). - tr : float, optional + tr : positive float, optional The Repetition Time of the BOLD data. If None, will extract the TR from NIFTI header (default None). use_afni : bool, optional @@ -52,7 +52,7 @@ class AmplitudeLowFrequencyFluctuationParcels( Notes ----- - The `tr` parameter is crucial for the correctness of fALFF/ALFF + The ``tr`` parameter is crucial for the correctness of fALFF/ALFF computation. If a dataset is correctly preprocessed, the TR should be extracted from the NIFTI without any issue. However, it has been reported that some preprocessed data might not have the correct TR in @@ -61,8 +61,6 @@ class AmplitudeLowFrequencyFluctuationParcels( ALFF/fALFF are computed using a bandpass butterworth filter. See :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more details. - - """ def __init__( diff --git a/junifer/markers/falff/falff_spheres.py b/junifer/markers/falff/falff_spheres.py index d25f43ffc..38c59f26e 100644 --- a/junifer/markers/falff/falff_spheres.py +++ b/junifer/markers/falff/falff_spheres.py @@ -30,11 +30,11 @@ class AmplitudeLowFrequencyFluctuationSpheres( for more information (default None). fractional : bool Whether to compute fractional ALFF. - highpass : float, optional + highpass : positive float, optional The highpass cutoff frequency for the bandpass filter (default 0.01). - lowpass : float + lowpass : positive float, optional The lowpass cutoff frequency for the bandpass filter (default 0.1). - tr : float, optional + tr : positive float, optional The Repetition Time of the BOLD data. If None, will extract the TR from NIFTI header (default None). use_afni : bool, optional @@ -56,7 +56,7 @@ class AmplitudeLowFrequencyFluctuationSpheres( Notes ----- - The `tr` parameter is crucial for the correctness of fALFF/ALFF + The ``tr`` parameter is crucial for the correctness of fALFF/ALFF computation. If a dataset is correctly preprocessed, the TR should be extracted from the NIFTI without any issue. However, it has been reported that some preprocessed data might not have the correct TR in @@ -65,8 +65,6 @@ class AmplitudeLowFrequencyFluctuationSpheres( ALFF/fALFF are computed using a bandpass butterworth filter. See :func:`scipy.signal.butter` and :func:`scipy.signal.filtfilt` for more details. - - """ def __init__( @@ -119,6 +117,7 @@ class AmplitudeLowFrequencyFluctuationSpheres( * ``data`` : the actual computed values as a numpy.ndarray * ``columns`` : the column labels for the computed values as a list + """ pa = SphereAggregation( coords=self.coords, diff --git a/junifer/markers/falff/tests/test_falff_parcels.py b/junifer/markers/falff/tests/test_falff_parcels.py index a4cb63fcc..aa12c0cdf 100644 --- a/junifer/markers/falff/tests/test_falff_parcels.py +++ b/junifer/markers/falff/tests/test_falff_parcels.py @@ -18,6 +18,9 @@ from junifer.storage import SQLiteFeatureStorage from junifer.utils import logger +_PARCELLATION = "Schaefer100x7" + + def test_AmplitudeLowFrequencyFluctuationParcels_python() -> None: """Test AmplitudeLowFrequencyFluctuationParcels using python.""" # Get the SPM auditory data: @@ -28,7 +31,7 @@ def test_AmplitudeLowFrequencyFluctuationParcels_python() -> None: input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", + parcellation=_PARCELLATION, method="mean", use_afni=False, fractional=False, @@ -52,7 +55,7 @@ def test_AmplitudeLowFrequencyFluctuationParcels_afni() -> None: input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", + parcellation=_PARCELLATION, method="mean", use_afni=True, fractional=False, @@ -65,7 +68,7 @@ def test_AmplitudeLowFrequencyFluctuationParcels_afni() -> None: # Again, should be blazing fast marker = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", method="mean", fractional=False + parcellation=_PARCELLATION, method="mean", fractional=False ) assert marker.use_afni is None afni_values2 = marker.fit_transform(input)["BOLD"]["data"] @@ -96,7 +99,7 @@ def test_AmplitudeLowFrequencyFluctuationParcels_python_vs_afni( input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker_python = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", + parcellation=_PARCELLATION, method="mean", use_afni=False, fractional=fractional, @@ -108,7 +111,7 @@ def test_AmplitudeLowFrequencyFluctuationParcels_python_vs_afni( assert python_values.shape == (1, 100) marker_afni = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", + parcellation=_PARCELLATION, method="mean", use_afni=True, fractional=fractional, @@ -140,7 +143,7 @@ def test_AmplitudeLowFrequencyFluctuationParcels_storage( input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationParcels( - parcellation="Schaefer100x7", + parcellation=_PARCELLATION, method="mean", use_afni=False, fractional=True, diff --git a/junifer/markers/falff/tests/test_falff_spheres.py b/junifer/markers/falff/tests/test_falff_spheres.py index a8387367a..25705c7c4 100644 --- a/junifer/markers/falff/tests/test_falff_spheres.py +++ b/junifer/markers/falff/tests/test_falff_spheres.py @@ -1,4 +1,4 @@ -"""Provide test for parcel-aggregated (f)ALFF.""" +"""Provide test for sphere-aggregated (f)ALFF.""" # Authors: Federico Raimondo # Synchon Mandal @@ -19,6 +19,9 @@ from junifer.storage import SQLiteFeatureStorage from junifer.utils import logger +_COORDINATES = "DMNBuckner" + + def test_AmplitudeLowFrequencyFluctuationSpheres_python() -> None: """Test AmplitudeLowFrequencyFluctuationSpheres using python.""" # Get the SPM auditory data: @@ -29,7 +32,7 @@ def test_AmplitudeLowFrequencyFluctuationSpheres_python() -> None: input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationSpheres( - coords="DMNBuckner", + coords=_COORDINATES, radius=5, method="mean", use_afni=False, @@ -54,7 +57,7 @@ def test_AmplitudeLowFrequencyFluctuationSpheres_afni() -> None: input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationSpheres( - coords="DMNBuckner", + coords=_COORDINATES, radius=5, method="mean", use_afni=True, @@ -68,7 +71,7 @@ def test_AmplitudeLowFrequencyFluctuationSpheres_afni() -> None: # Again, should be blazing fast marker = AmplitudeLowFrequencyFluctuationSpheres( - coords="DMNBuckner", + coords=_COORDINATES, radius=5, method="mean", fractional=False, @@ -101,7 +104,7 @@ def test_AmplitudeLowFrequencyFluctuationSpheres_python_vs_afni( input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker_python = AmplitudeLowFrequencyFluctuationSpheres( - coords="DMNBuckner", + coords=_COORDINATES, radius=5, method="mean", use_afni=False, @@ -114,7 +117,7 @@ def test_AmplitudeLowFrequencyFluctuationSpheres_python_vs_afni( assert python_values.shape == (1, 6) marker_afni = AmplitudeLowFrequencyFluctuationSpheres( - coords="DMNBuckner", + coords=_COORDINATES, radius=5, method="mean", use_afni=True, @@ -147,7 +150,7 @@ def test_AmplitudeLowFrequencyFluctuationSpheres_storage( input = DefaultDataReader().fit_transform(input) # Create ParcelAggregation object marker = AmplitudeLowFrequencyFluctuationSpheres( - coords="DMNBuckner", + coords=_COORDINATES, radius=5, method="mean", use_afni=False, -- 2.52.0 From fb559de90da1160c5a0e84e6b392345ce96628d9 Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Mon, 19 Dec 2022 13:39:52 +0100 Subject: [PATCH 09/10] one more docstring --- junifer/markers/falff/tests/test_falff_spheres.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/junifer/markers/falff/tests/test_falff_spheres.py b/junifer/markers/falff/tests/test_falff_spheres.py index 25705c7c4..6d10ae463 100644 --- a/junifer/markers/falff/tests/test_falff_spheres.py +++ b/junifer/markers/falff/tests/test_falff_spheres.py @@ -91,7 +91,7 @@ def test_AmplitudeLowFrequencyFluctuationSpheres_afni() -> None: def test_AmplitudeLowFrequencyFluctuationSpheres_python_vs_afni( fractional: bool, ) -> None: - """Test AmplitudeLowFrequencyFluctuationSpheres using python. + """Test AmplitudeLowFrequencyFluctuationSpheres python vs afni results. Parameters ---------- -- 2.52.0 From 08154f9b100c71bba1ac7259f2740a7fc0e6e5d0 Mon Sep 17 00:00:00 2001 From: Fede Raimondo Date: Mon, 19 Dec 2022 13:54:38 +0100 Subject: [PATCH 10/10] Fix docs.conf to allow for positive/negative floats --- docs/conf.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/conf.py b/docs/conf.py index f738b1e68..a11f4ed53 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -155,6 +155,8 @@ numpydoc_xref_ignore = { "class", "objects", "Engine", + "positive", + "negative", } # numpydoc_validation_checks = { # "all", -- 2.52.0