diff --git a/.all-contributorsrc b/.all-contributorsrc
index 46f9a58..1e8e6aa 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -110,6 +110,20 @@
"review",
"userTesting"
]
+ },
+ {
+ "login": "maestroque",
+ "name": "George Kikas",
+ "avatar_url": "https://avatars.githubusercontent.com/u/74024609?v=4",
+ "profile": "https://github.com/maestroque",
+ "contributions": [
+ "code",
+ "ideas",
+ "infra",
+ "docs",
+ "bug",
+ "review"
+ ]
}
],
"contributorsPerLine": 8,
diff --git a/.circleci/config.yml b/.circleci/config.yml
index fc36d30..dd45175 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -15,14 +15,14 @@ orbs:
# Define a job to be invoked later in a workflow.
# See: https://circleci.com/docs/2.0/configuration-reference/#jobs
jobs:
- test37: # This is the name of the job, feel free to change it to better match what you're trying to do!
+ test39: # This is the name of the job, feel free to change it to better match what you're trying to do!
# These next lines defines a Docker executors: https://circleci.com/docs/2.0/executor-types/
# You can specify an image from Dockerhub or use one of the convenience images from CircleCI's Developer Hub
# A list of available CircleCI Docker convenience images are available here: https://circleci.com/developer/images/image/cimg/python
# The executor is the environment in which the steps below will be executed - below will use a python 3.6.14 container
# Change the version below to your required version of python
docker:
- - image: cimg/python:3.7
+ - image: cimg/python:3.9
working_directory: /tmp/src/phys2denoise
resource_class: medium
# Checkout the code as the first step. This is a dedicated CircleCI step.
@@ -46,7 +46,7 @@ jobs:
command: |
pytest --cov=./phys2denoise
mkdir /tmp/src/coverage
- mv ./.coverage /tmp/src/coverage/.coverage.py37
+ mv ./.coverage /tmp/src/coverage/.coverage.py39
- store_artifacts:
path: /tmp/src/coverage
# Persist the specified paths (workspace/echo-output) into the workspace for use in downstream job.
@@ -56,11 +56,11 @@ jobs:
root: /tmp
# Must be relative path from root
paths:
- - src/coverage/.coverage.py37
+ - src/coverage/.coverage.py39
- test310:
+ test311:
docker:
- - image: cimg/python:3.10
+ - image: cimg/python:3.11
working_directory: /tmp/src/phys2denoise
resource_class: medium
steps:
@@ -75,17 +75,17 @@ jobs:
command: |
pytest --cov=./phys2denoise
mkdir /tmp/src/coverage
- mv ./.coverage /tmp/src/coverage/.coverage.py310
+ mv ./.coverage /tmp/src/coverage/.coverage.py311
- store_artifacts:
path: /tmp/src/coverage
- persist_to_workspace:
root: /tmp
paths:
- - src/coverage/.coverage.py310
+ - src/coverage/.coverage.py311
style_check:
docker:
- - image: cimg/python:3.7
+ - image: cimg/python:3.9
working_directory: /tmp/src/phys2denoise
resource_class: small
steps:
@@ -105,7 +105,7 @@ jobs:
merge_coverage:
working_directory: /tmp/src/phys2denoise
docker:
- - image: cimg/python:3.10
+ - image: cimg/python:3.11
resource_class: small
steps:
- attach_workspace:
@@ -133,13 +133,13 @@ workflows:
# Inside the workflow, you define the jobs you want to run.
jobs:
- style_check
- - test37:
+ - test39:
requires:
- style_check
- - test310:
+ - test311:
requires:
- style_check
- merge_coverage:
requires:
- - test37
- - test310
+ - test39
+ - test311
diff --git a/docs/user_guide/metrics.rst b/docs/user_guide/metrics.rst
index e387600..779f31b 100644
--- a/docs/user_guide/metrics.rst
+++ b/docs/user_guide/metrics.rst
@@ -71,6 +71,10 @@ example shows how to compute the heart rate and the heart rate variability using
from phys2denoise.metrics.chest_belt import respiratory_variance_time
- # Given that the respiratory signal is stored in `data`, the peaks in `peaks`, the troughs in `troughs`
+ # Given that the respiratory signal is stored in `data` (which is not a physio.Physio instance), the peaks in `peaks`, the troughs in `troughs`
# and the sample rate in `sample_rate`
- _, rvt = respiratory_variance_time(data, peaks, troughs, sample_rate)
+ rvt = respiratory_variance_time(data, peaks, troughs, sample_rate)
+
+The computed respiratory variance time is stored in the variable ``rvt``. An internal check is performed to verify if the input data is a Physio object or not
+determining the appropriate output format. If the input is not a ``Physio`` object, the output will be a numpy array only containing the computed metric. Otherwise,
+the output will be a tuple with the updated Physio object and the computed metric.
diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 40ebc47..5b7588d 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -1,14 +1,24 @@
"""Denoising metrics for cardio recordings."""
import numpy as np
+from loguru import logger
+from physutils import io, physio
from .. import references
from ..due import due
from .responses import crf
from .utils import apply_function_in_sliding_window as afsw
-from .utils import convolve_and_rescale
-
-
-def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure="mean"):
+from .utils import convolve_and_rescale, return_physio_or_metric
+
+
+def _cardiac_metrics(
+ data,
+ metric,
+ peaks=None,
+ fs=None,
+ window=6,
+ central_measure="mean",
+ **kwargs,
+):
"""
Compute cardiac metrics.
@@ -21,14 +31,14 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
Parameters
----------
- card : list or 1D numpy.ndarray
- Timeseries of recorded cardiac signal
- peaks : list or 1D numpy.ndarray
- array of peak indexes for card.
- samplerate : float
- Sampling rate for card, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded cardiac signal
metrics : "hbi", "hr", "hrv", string
Cardiac metric(s) to calculate.
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
+ peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+ Indices of peaks in `data`
window : float, optional
Size of the sliding window, in seconds.
Default is 6.
@@ -71,8 +81,32 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
Annual International Conference of the IEEE Engineering in Medicine and
Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
"""
+ if isinstance(data, physio.Physio):
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=True, copy=True)
+ elif fs is not None and peaks is not None:
+ data = io.load_physio(data, fs=fs)
+ data._metadata["peaks"] = peaks
+ else:
+ raise ValueError(
+ """
+ To use this function you should either provide a Physio object
+ with existing peaks metadata (e.g. using the peakdet module), or
+ by providing the physiological data timeseries, the sampling frequency,
+ and the peak indices separately.
+ """
+ )
+ if data.peaks.size == 0:
+ raise ValueError(
+ """
+ Peaks must be a non-empty list.
+ Make sure to run peak detection on your physiological data first,
+ using the peakdet module, or other software of your choice.
+ """
+ )
+
# Convert window to samples, but halves it.
- halfwindow_samples = int(round(window * samplerate / 2))
+ halfwindow_samples = int(round(window * data.fs / 2))
if central_measure in ["mean", "average", "avg"]:
central_measure_operator = np.mean
@@ -85,14 +119,17 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
f" {central_measure} is not a supported metric of centrality."
)
- idx_arr = np.arange(len(card))
+ idx_arr = np.arange(len(data))
idx_min = afsw(idx_arr, np.min, halfwindow_samples)
idx_max = afsw(idx_arr, np.max, halfwindow_samples)
- card_met = np.empty_like(card)
+ card_met = np.empty_like(data)
for n, i in enumerate(idx_min):
diff = (
- np.diff(peaks[np.logical_and(peaks >= i, peaks <= idx_max[n])]) / samplerate
+ np.diff(
+ data.peaks[np.logical_and(data.peaks >= i, data.peaks <= idx_max[n])]
+ )
+ / data.fs
)
if metric == "hbi":
card_met[n] = central_measure_operator(diff) if diff.size > 0 else 0
@@ -109,13 +146,15 @@ def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure=
card_met[np.isnan(card_met)] = 0.0
# Convolve with crf and rescale
- card_met = convolve_and_rescale(card_met, crf(samplerate), rescale="rescale")
+ card_met = convolve_and_rescale(card_met, crf(data.fs), rescale="rescale")
- return card_met
+ return data, card_met
@due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
-def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
+@return_physio_or_metric()
+@physio.make_operation()
+def heart_rate(data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs):
"""
Compute average heart rate (HR) in a sliding window.
@@ -125,12 +164,12 @@ def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
Parameters
----------
- card : list or 1D numpy.ndarray
- Timeseries of recorded cardiac signal
- peaks : list or 1D numpy.ndarray
- array of peak indexes for card.
- samplerate : float
- Sampling rate for card, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded cardiac signal
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
+ peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+ Indices of peaks in `data`
window : float, optional
Size of the sliding window, in seconds.
Default is 6.
@@ -167,13 +206,25 @@ def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
Annual International Conference of the IEEE Engineering in Medicine and
Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
"""
- return _cardiac_metrics(
- card, peaks, samplerate, metric="hrv", window=6, central_measure="mean"
+ data, hr = _cardiac_metrics(
+ data,
+ metric="hr",
+ peaks=peaks,
+ fs=fs,
+ window=window,
+ central_measure=central_measure,
+ **kwargs,
)
+ return data, hr
+
@due.dcite(references.PINHERO_ET_AL_2016)
-def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="mean"):
+@return_physio_or_metric()
+@physio.make_operation()
+def heart_rate_variability(
+ data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs
+):
"""
Compute average heart rate variability (HRV) in a sliding window.
@@ -183,12 +234,12 @@ def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="m
Parameters
----------
- card : list or 1D numpy.ndarray
- Timeseries of recorded cardiac signal
- peaks : list or 1D numpy.ndarray
- array of peak indexes for card.
- samplerate : float
- Sampling rate for card, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded cardiac signal
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
+ peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+ Indices of peaks in `data`
window : float, optional
Size of the sliding window, in seconds.
Default is 6.
@@ -223,24 +274,36 @@ def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="m
Annual International Conference of the IEEE Engineering in Medicine and
Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
"""
- return _cardiac_metrics(
- card, peaks, samplerate, metric="hrv", window=6, central_measure="std"
+ data, hrv = _cardiac_metrics(
+ data,
+ metric="hrv",
+ peaks=peaks,
+ fs=fs,
+ window=window,
+ central_measure=central_measure,
+ **kwargs,
)
+ return data, hrv
+
@due.dcite(references.CHEN_2020)
-def heart_beat_interval(card, peaks, samplerate, window=6, central_measure="mean"):
+@return_physio_or_metric()
+@physio.make_operation()
+def heart_beat_interval(
+ data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs
+):
"""
Compute average heart beat interval (HBI) in a sliding window.
Parameters
----------
- card : list or 1D numpy.ndarray
- Timeseries of recorded cardiac signal
- peaks : list or 1D numpy.ndarray
- array of peak indexes for card.
- samplerate : float
- Sampling rate for card, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded cardiac signal
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
+ peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+ Indices of peaks in `data`
window : float, optional
Size of the sliding window, in seconds.
Default is 6.
@@ -272,12 +335,22 @@ def heart_beat_interval(card, peaks, samplerate, window=6, central_measure="mean
.. [1] J. E. Chen et al., "Resting-state "physiological networks"", Neuroimage,
vol. 213, pp. 116707, 2020.
"""
- return _cardiac_metrics(
- card, peaks, samplerate, metric="hbi", window=6, central_measure="mean"
+ data, hbi = _cardiac_metrics(
+ data,
+ metric="hbi",
+ peaks=peaks,
+ fs=fs,
+ window=window,
+ central_measure=central_measure,
+ **kwargs,
)
+ return data, hbi
-def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
+
+@return_physio_or_metric()
+@physio.make_operation()
+def cardiac_phase(data, slice_timings, n_scans, t_r, peaks=None, fs=None, **kwargs):
"""Calculate cardiac phase from cardiac peaks.
Assumes that timing of cardiac events are given in same units
@@ -285,27 +358,53 @@ def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
Parameters
----------
- peaks : 1D array_like
- Cardiac peak times, in seconds.
- sample_rate : float
- Sample rate of physio, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded cardiac signal
slice_timings : 1D array_like
Slice times, in seconds.
n_scans : int
Number of volumes in the imaging run.
t_r : float
Sampling rate of the imaging run, in seconds.
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
+ peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+ Indices of peaks in `data`
Returns
-------
phase_card : array_like
Cardiac phase signal, of shape (n_scans,)
"""
+ if isinstance(data, physio.Physio):
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=True, copy=True)
+ elif fs is not None and peaks is not None:
+ data = io.load_physio(data, fs=fs)
+ data._metadata["peaks"] = peaks
+ else:
+ raise ValueError(
+ """
+ To use this function you should either provide a Physio object
+ with existing peaks metadata (e.g. using the peakdet module), or
+ by providing the physiological data timeseries, the sampling frequency,
+ and the peak indices separately.
+ """
+ )
+ if data.peaks.size == 0:
+ raise ValueError(
+ """
+ Peaks must be a non-empty list.
+ Make sure to run peak detection on your physiological data first,
+ using the peakdet module, or other software of your choice.
+ """
+ )
+
assert slice_timings.ndim == 1, "Slice times must be a 1D array"
n_slices = np.size(slice_timings)
phase_card = np.zeros((n_scans, n_slices))
- card_peaks_sec = peaks / sample_rate
+ card_peaks_sec = data.peaks / data.fs
for i_slice in range(n_slices):
# generate slice acquisition timings across all scans
times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice]
@@ -332,4 +431,4 @@ def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
) / (t2 - t1)
phase_card[:, i_slice] = phase_card_crSlice
- return phase_card
+ return data, phase_card
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index eeeb08a..4e080fc 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -1,6 +1,7 @@
"""Denoising metrics for chest belt recordings."""
import numpy as np
import pandas as pd
+from physutils import io, physio
from scipy.interpolate import interp1d
from scipy.stats import zscore
@@ -8,11 +9,15 @@
from ..due import due
from .responses import rrf
from .utils import apply_function_in_sliding_window as afsw
-from .utils import convolve_and_rescale, rms_envelope_1d
+from .utils import convolve_and_rescale, return_physio_or_metric, rms_envelope_1d
@due.dcite(references.BIRN_2006)
-def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 12)):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_variance_time(
+ data, peaks=None, troughs=None, fs=None, lags=(0, 4, 8, 12), **kwargs
+):
"""
Implement the Respiratory Variance over Time (Birn et al. 2006).
@@ -20,14 +25,14 @@ def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 1
Parameters
----------
- resp: array_like
- respiratory belt data - samples x 1
- peaks: array_like
- peaks found by peakdet algorithm
- troughs: array_like
- troughs found by peakdet algorithm
- samplerate: float
- sample rate in hz of respiratory belt
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded respiratory signal
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
+ peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+ Indices of peaks in `data`
+ troughs : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
+ Indices of troughs in `data`
lags: tuple
lags in seconds of the RVT output. Default is 0, 4, 8, 12.
@@ -42,13 +47,38 @@ def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 1
respiratory-variation-related fluctuations from neuronal-activity-related
fluctuations in fMRI”, NeuroImage, vol.31, pp. 1536-1548, 2006.
"""
- timestep = 1 / samplerate
+ if isinstance(data, physio.Physio):
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=True, copy=True)
+ elif fs is not None and peaks is not None and troughs is not None:
+ data = io.load_physio(data, fs=fs)
+ data._metadata["peaks"] = peaks
+ data._metadata["troughs"] = troughs
+ else:
+ raise ValueError(
+ """
+ To use this function you should either provide a Physio object
+ with existing peaks and troughs metadata (e.g. using the peakdet module), or
+ by providing the physiological data timeseries, the sampling frequency,
+ and the peak and trough indices separately.
+ """
+ )
+ if data.peaks.size == 0 or data.troughs.size == 0:
+ raise ValueError(
+ """
+ Peaks and troughs must be non-empty lists.
+ Make sure to run peak/trough detection on your physiological data first,
+ using the peakdet module, or other software of your choice.
+ """
+ )
+
+ timestep = 1 / data.fs
# respiration belt timing
- time = np.arange(0, len(resp) * timestep, timestep)
- peak_vals = resp[peaks]
- trough_vals = resp[troughs]
- peak_time = time[peaks]
- trough_time = time[troughs]
+ time = np.arange(0, len(data) * timestep, timestep)
+ peak_vals = data[data.peaks]
+ trough_vals = data[data.troughs]
+ peak_time = time[data.peaks]
+ trough_time = time[data.troughs]
mid_peak_time = (peak_time[:-1] + peak_time[1:]) / 2
period = np.diff(peak_time)
# interpolate peak values over all timepoints
@@ -77,17 +107,19 @@ def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 1
)
rvt_lags[:, ind] = temp_rvt
- return rvt_lags
+ return data, rvt_lags
@due.dcite(references.POWER_2018)
-def respiratory_pattern_variability(resp, window):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_pattern_variability(data, window, **kwargs):
"""Calculate respiratory pattern variability.
Parameters
----------
- resp : str or 1D numpy.ndarray
- Tiemseries representing respiration activity.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded respiratory signal
window : int
Window length in samples.
@@ -111,27 +143,33 @@ def respiratory_pattern_variability(resp, window):
data," Proceedings of the National Academy of Sciences, issue 9, vol.
115, pp. 2105-2114, 2018.
"""
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=False, copy=True)
+
# First, z-score respiratory traces
- resp_z = zscore(resp)
+ resp_z = zscore(data.data)
# Collect upper envelope
rpv_upper_env = rms_envelope_1d(resp_z, window)
# Calculate standard deviation
rpv_val = np.std(rpv_upper_env)
- return rpv_val
+
+ return data, rpv_val
@due.dcite(references.POWER_2020)
-def env(resp, samplerate, window=10):
+@return_physio_or_metric()
+@physio.make_operation()
+def env(data, fs=None, window=10, **kwargs):
"""Calculate respiratory pattern variability across a sliding window.
Parameters
----------
- resp : (X,) :obj:`numpy.ndarray`
- A 1D array with the respiratory belt time series.
- samplerate : :obj:`float`
- Sampling rate for resp, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded respiratory signal
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
window : :obj:`int`, optional
Size of the sliding window, in seconds.
Default is 10.
@@ -155,30 +193,67 @@ def env(resp, samplerate, window=10):
young adults scanned at rest, including systematic changes and 'missed'
deep breaths," Neuroimage, vol. 204, 2020.
"""
+
+ def _respiratory_pattern_variability(data, window):
+ """
+ Respiratory pattern variability function without utilizing
+ the physutils.Physio object history, only to be used within
+ the context of this function. This is done to only store the
+ chest_belt.env operation call to the history and not all the
+ subsequent sub-operations
+ """
+ # First, z-score respiratory traces
+ resp_z = zscore(data)
+
+ # Collect upper envelope
+ rpv_upper_env = rms_envelope_1d(resp_z, window)
+
+ # Calculate standard deviation
+ rpv_val = np.std(rpv_upper_env)
+ return rpv_val
+
+ if isinstance(data, physio.Physio):
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=True, copy=True)
+ elif fs is not None:
+ data = io.load_physio(data, fs=fs)
+ else:
+ raise ValueError(
+ """
+ To use this function you should either provide a Physio object
+ with the sampling frequency encapsulated, or
+ by providing the physiological data timeseries and the sampling
+ frequency separately.
+ """
+ )
+
# Convert window to Hertz
- window = int(window * samplerate)
+ window = int(window * data.fs)
# Calculate RPV across a rolling window
env_arr = (
- pd.Series(resp)
+ pd.Series(data.data)
.rolling(window=window, center=True)
- .apply(respiratory_pattern_variability, args=(window,))
+ .apply(_respiratory_pattern_variability, args=(window,))
)
env_arr[np.isnan(env_arr)] = 0.0
- return env_arr
+
+ return data, env_arr
@due.dcite(references.CHANG_GLOVER_2009)
-def respiratory_variance(resp, samplerate, window=6):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_variance(data, fs=None, window=6, **kwargs):
"""Calculate respiratory variance.
Parameters
----------
- resp : (X,) :obj:`numpy.ndarray`
- A 1D array with the respiratory belt time series.
- samplerate : :obj:`float`
- Sampling rate for resp, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded respiratory signal
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
window : :obj:`int`, optional
Size of the sliding window, in seconds.
Default is 6.
@@ -206,49 +281,81 @@ def respiratory_variance(resp, samplerate, window=6):
end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
issue 4, vol. 47, pp. 1381-1393, 2009.
"""
+ if isinstance(data, physio.Physio):
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=True, copy=True)
+ elif fs is not None:
+ data = io.load_physio(data, fs=fs)
+ else:
+ raise ValueError(
+ """
+ To use this function you should either provide a Physio object
+ with the sampling frequency encapsulated, or
+ by providing the physiological data timeseries and the sampling
+ frequency separately.
+ """
+ )
+
# Convert window to Hertz
- halfwindow_samples = int(round(window * samplerate / 2))
+ halfwindow_samples = int(round(window * data.fs / 2))
# Raw respiratory variance
- rv_arr = afsw(resp, np.std, halfwindow_samples)
+ rv_arr = afsw(data, np.std, halfwindow_samples)
# Convolve with rrf
- rv_out = convolve_and_rescale(rv_arr, rrf(samplerate), rescale="zscore")
+ rv_out = convolve_and_rescale(rv_arr, rrf(data.fs), rescale="zscore")
- return rv_out
+ return data, rv_out
-def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
+@return_physio_or_metric()
+@physio.make_operation()
+def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None, **kwargs):
"""Calculate respiratory phase from respiratory signal.
Parameters
----------
- resp : 1D array_like
- Respiratory signal.
- sample_rate : float
- Sample rate of physio, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded respiratory signal
n_scans
Number of volumes in the imaging run.
slice_timings
Slice times, in seconds.
t_r
Sample rate of the imaging run, in seconds.
+ fs : float, optional if data is a physutils.Physio
+ Sampling rate of `data` in Hz
Returns
-------
phase_resp : array_like
Respiratory phase signal, of shape (n_scans, n_slices).
"""
+ if isinstance(data, physio.Physio):
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=True, copy=True)
+ elif fs is not None:
+ data = io.load_physio(data, fs=fs)
+ else:
+ raise ValueError(
+ """
+ To use this function you should either provide a Physio object
+ with the sampling frequency encapsulated, or
+ by providing the physiological data timeseries and the sampling
+ frequency separately.
+ """
+ )
+
assert slice_timings.ndim == 1, "Slice times must be a 1D array"
n_slices = np.size(slice_timings)
phase_resp = np.zeros((n_scans, n_slices))
# generate histogram from respiratory signal
# TODO: Replace with numpy.histogram
- resp_hist, resp_hist_bins = np.histogram(resp, bins=100)
+ resp_hist, resp_hist_bins = np.histogram(data, bins=100)
# first compute derivative of respiration signal
- resp_diff = np.diff(resp, n=1)
+ resp_diff = np.diff(data, n=1)
for i_slice in range(n_slices):
# generate slice acquisition timings across all scans
@@ -256,15 +363,15 @@ def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
phase_resp_crSlice = np.zeros(n_scans)
for j_scan in range(n_scans):
iphys = int(
- max([1, round(times_crSlice[j_scan] * sample_rate)])
+ max([1, round(times_crSlice[j_scan] * data.fs)])
) # closest idx in resp waveform
iphys = min([iphys, len(resp_diff)]) # cannot be longer than resp_diff
- thisBin = np.argmin(abs(resp[iphys] - resp_hist_bins))
+ thisBin = np.argmin(abs(data[iphys] - resp_hist_bins))
numerator = np.sum(resp_hist[0:thisBin])
phase_resp_crSlice[j_scan] = (
- np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(resp))
+ np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(data))
)
phase_resp[:, i_slice] = phase_resp_crSlice
- return phase_resp
+ return data, phase_resp
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index 7a298b2..fa28bae 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -1,34 +1,35 @@
"""These functions compute RETROICOR regressors (Glover et al. 2000)."""
import numpy as np
+from physutils import io, physio
from .. import references
from ..due import due
from .cardiac import cardiac_phase
from .chest_belt import respiratory_phase
+from .utils import return_physio_or_metric
@due.dcite(references.GLOVER_2000)
+@return_physio_or_metric()
+@physio.make_operation()
def retroicor(
- physio,
- sample_rate,
+ data,
t_r,
n_scans,
slice_timings,
n_harmonics,
- card=False,
- resp=False,
+ physio_type=None,
+ fs=None,
+ cardiac_peaks=None,
+ **kwargs,
):
"""Compute RETROICOR regressors.
Parameters
----------
- physio : array_like
- 1D array, whether cardiac or respiratory signal.
- If cardiac, the array is a set of peaks in seconds.
- If respiratory, the array is the actual respiratory signal.
- sample_rate : float
- Physio sample rate, in Hertz.
+ data : physutils.Physio, np.ndarray, or array-like object
+ Object containing the timeseries of the recorded respiratory or cardiac signal
t_r : float
Imaging sample rate, in seconds.
n_scans : int
@@ -59,6 +60,44 @@ def retroicor(
correction of physiological motion effects in fMRI: RETROICOR“, Magn. Reson. Med.,
issue 1, vol. 44, pp. 162-167, 2000.
"""
+ if isinstance(data, physio.Physio):
+ # Initialize physio object
+ data = physio.check_physio(data, ensure_fs=True, copy=True)
+ if data.physio_type is None and physio_type is not None:
+ data._physio_type = physio_type
+ elif data.physio_type is None and physio_type is None:
+ raise ValueError(
+ """
+ Since the provided Physio object does not specify a `physio_type`,
+ this function's `physio_type` parameter must be specified as a
+ value from {'cardiac', 'respiratory'}
+ """
+ )
+
+ elif fs is not None and physio_type is not None:
+ data = io.load_physio(data, fs=fs)
+ data._physio_type = physio_type
+ if data.physio_type == "cardiac":
+ data._metadata["peaks"] = cardiac_peaks
+ else:
+ raise ValueError(
+ """
+ To use this function you should either provide a Physio object
+ with existing peaks metadata if it describes a cardiac signal
+ (e.g. using the peakdet module), or
+ by providing the physiological data timeseries, the sampling frequency,
+ the physio_type and the peak indices separately.
+ """
+ )
+ if not data.peaks and data.physio_type == "cardiac":
+ raise ValueError(
+ """
+ Peaks must be a non-empty list for cardiac data.
+ Make sure to run peak detection on your cardiac data first,
+ using the peakdet module, or other software of your choice.
+ """
+ )
+
n_slices = np.shape(slice_timings) # number of slices
# initialize output variables
@@ -73,18 +112,17 @@ def retroicor(
# Compute physiological phases using the timings of physio events (e.g. peaks)
# slice sampling times
- if card:
+ if data.physio_type == "cardiac":
phase[:, i_slice] = cardiac_phase(
- physio,
+ data,
crslice_timings,
n_scans,
t_r,
)
- if resp:
+ if data.physio_type == "respiratory":
phase[:, i_slice] = respiratory_phase(
- physio,
- sample_rate,
+ data,
n_scans,
slice_timings,
t_r,
@@ -99,4 +137,9 @@ def retroicor(
(j_harm + 1) * phase[i_slice]
)
- return retroicor_regressors, phase
+ data._computed_metrics["retroicor_regressors"] = dict(
+ metrics=retroicor_regressors, phase=phase
+ )
+ retroicor_regressors = dict(metrics=retroicor_regressors, phase=phase)
+
+ return data, retroicor_regressors
diff --git a/phys2denoise/metrics/responses.py b/phys2denoise/metrics/responses.py
index aabcec9..1a2aeee 100644
--- a/phys2denoise/metrics/responses.py
+++ b/phys2denoise/metrics/responses.py
@@ -2,6 +2,7 @@
import logging
import numpy as np
+from loguru import logger
from .. import references
from ..due import due
@@ -56,10 +57,12 @@ def _crf(t):
) * np.exp(-0.5 * (((t - 12) ** 2) / 9))
return rf
- time_stamps = np.arange(0, time_length, 1 / samplerate)
+ time_stamps = np.arange(0, time_length, samplerate)
time_stamps -= onset
+ logger.debug(f"Time stamps: {time_stamps}")
crf_arr = _crf(time_stamps)
crf_arr = crf_arr / max(abs(crf_arr))
+ logger.debug(f"CRF: {crf_arr}")
if inverse:
return -crf_arr
@@ -135,7 +138,7 @@ def _rrf(t):
rf = 0.6 * t**2.1 * np.exp(-t / 1.6) - 0.0023 * t**3.54 * np.exp(-t / 4.25)
return rf
- time_stamps = np.arange(0, time_length, 1 / samplerate)
+ time_stamps = np.arange(0, time_length, samplerate)
time_stamps -= onset
rrf_arr = _rrf(time_stamps)
rrf_arr = rrf_arr / max(abs(rrf_arr))
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index 34fb21b..e21a085 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -1,8 +1,12 @@
"""Miscellaneous utility functions for metric calculation."""
+import functools
+import inspect
import logging
import numpy as np
+from loguru import logger
from numpy.lib.stride_tricks import sliding_window_view as swv
+from physutils.physio import Physio
from scipy.interpolate import interp1d
from scipy.stats import zscore
@@ -33,7 +37,7 @@ def print_metric_call(metric, args):
msg = f"{msg}\n"
- LGR.info(msg)
+ logger.info(msg)
def mirrorpad_1d(arr, buffer=250):
@@ -58,7 +62,11 @@ def mirrorpad_1d(arr, buffer=250):
post_mirror = np.take(mirror, idx, axis=0)
except IndexError:
len(arr)
- LGR.warning(
+ # LGR.warning(
+ # f"Requested buffer size ({buffer}) is longer than input array length "
+ # f"({len(arr)}). Fixing buffer size to array length."
+ # )
+ logger.warning(
f"Requested buffer size ({buffer}) is longer than input array length "
f"({len(arr)}). Fixing buffer size to array length."
)
@@ -332,3 +340,117 @@ def export_metric(
)
return fileprefix
+
+
+def return_physio_or_metric(*, return_physio=True):
+ """
+ Decorator to check if the input is a Physio object.
+
+ Parameters
+ ----------
+ func : function
+ Function to be decorated
+
+ Returns
+ -------
+ function
+ Decorated function
+ """
+
+ def determine_return_type(func):
+ convolved_metrics = [
+ "respiratory_variance",
+ "heart_rate",
+ "heart_rate_variability",
+ "heart_beat_interval",
+ ]
+
+ @functools.wraps(func)
+ def wrapper(*args, **kwargs):
+ physio, metric = func(*args, **kwargs)
+ default_args = get_default_args(func)
+ if isinstance(args[0], Physio):
+ if "lags" in kwargs:
+ has_lags = True if len(kwargs["lags"]) > 1 else False
+ elif "lags" in default_args:
+ has_lags = True if len(default_args["lags"]) > 1 else False
+ else:
+ has_lags = False
+
+ is_convolved = True if func.__name__ in convolved_metrics else False
+
+ physio._computed_metrics[func.__name__] = Metric(
+ func.__name__,
+ metric,
+ kwargs,
+ has_lags=has_lags,
+ is_convolved=is_convolved,
+ )
+
+ return_physio_value = kwargs.get("return_physio", return_physio)
+ if return_physio_value:
+ return physio
+ else:
+ return metric
+ else:
+ return metric
+
+ return wrapper
+
+ return determine_return_type
+
+
+def get_default_args(func):
+ # Get the signature of the function
+ sig = inspect.signature(func)
+
+ # Extract default values for each parameter
+ defaults = {
+ k: v.default
+ for k, v in sig.parameters.items()
+ if v.default is not inspect.Parameter.empty
+ }
+
+ return defaults
+
+
+class Metric:
+ def __init__(self, name, data, args, has_lags=False, is_convolved=False):
+ self.name = name
+ self._data = data
+ self._args = args
+ self._has_lags = has_lags
+ self._is_convolved = is_convolved
+
+ def __array__(self):
+ return self.data
+
+ def __getitem__(self, slicer):
+ return self.data[slicer]
+
+ def __len__(self):
+ return len(self.data)
+
+ @property
+ def ndim(self):
+ return self.data.ndim
+
+ @property
+ def shape(self):
+ return self.data.shape
+
+ @property
+ def data(self):
+ return self._data
+
+ @property
+ def args(self):
+ return self._args
+
+ @property
+ def has_lags(self):
+ return self._has_lags
+
+ @property
+ def is_convolved(self):
+ return self._is_convolved
diff --git a/phys2denoise/tests/conftest.py b/phys2denoise/tests/conftest.py
index fcf05f3..757a071 100644
--- a/phys2denoise/tests/conftest.py
+++ b/phys2denoise/tests/conftest.py
@@ -1,5 +1,20 @@
import numpy as np
import pytest
+from _pytest.logging import LogCaptureFixture
+from loguru import logger
+
+
+@pytest.fixture
+def caplog(caplog: LogCaptureFixture):
+ handler_id = logger.add(
+ caplog.handler,
+ format="{message}",
+ level=20,
+ filter=lambda record: record["level"].no >= caplog.handler.level,
+ enqueue=False,
+ )
+ yield caplog
+ logger.remove(handler_id)
@pytest.fixture(scope="module")
diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index bb10d76..2d776cf 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -1,5 +1,7 @@
"""Tests for phys2denoise.metrics.cardiac."""
import numpy as np
+from loguru import logger
+from physutils import physio
from phys2denoise.metrics import cardiac
@@ -7,20 +9,15 @@
def test_crf_smoke():
"""Basic smoke test for CRF calculation."""
samplerate = 0.01 # in seconds
- oversampling = 20
time_length = 20
onset = 0
tr = 0.72
crf_arr = cardiac.crf(
- samplerate,
- oversampling=oversampling,
- time_length=time_length,
- onset=onset,
- tr=tr,
+ samplerate, time_length=time_length, onset=onset, inverse=False
)
- pred_len = np.rint(time_length / (tr / oversampling))
+ pred_len = np.rint(time_length * 1 / samplerate)
assert crf_arr.ndim == 1
- assert crf_arr.shape == pred_len
+ assert len(crf_arr) == pred_len
def test_cardiac_phase_smoke():
@@ -30,12 +27,55 @@ def test_cardiac_phase_smoke():
sample_rate = 1 / 0.01
slice_timings = np.linspace(0, t_r, 22)[1:-1]
peaks = np.array([0.534, 0.577, 10.45, 20.66, 50.55, 90.22])
+ data = np.zeros(peaks.shape)
card_phase = cardiac.cardiac_phase(
- peaks,
- sample_rate=sample_rate,
+ data,
+ peaks=peaks,
+ fs=sample_rate,
slice_timings=slice_timings,
n_scans=n_scans,
t_r=t_r,
)
assert card_phase.ndim == 2
assert card_phase.shape == (n_scans, slice_timings.size)
+
+
+def test_cardiac_phase_smoke_physio_obj():
+ """Basic smoke test for cardiac phase calculation."""
+ t_r = 1.0
+ n_scans = 200
+ sample_rate = 1 / 0.01
+ slice_timings = np.linspace(0, t_r, 22)[1:-1]
+ peaks = np.array([0.534, 0.577, 10.45, 20.66, 50.55, 90.22])
+ data = np.zeros(peaks.shape)
+ phys = physio.Physio(data, sample_rate, physio_type="cardiac")
+ phys._metadata["peaks"] = peaks
+
+ # Test where the physio object is returned
+ phys = cardiac.cardiac_phase(
+ phys,
+ slice_timings=slice_timings,
+ n_scans=n_scans,
+ t_r=t_r,
+ )
+
+ assert phys.history[0][0] == "phys2denoise.metrics.cardiac.cardiac_phase"
+ assert phys.computed_metrics["cardiac_phase"].ndim == 2
+ assert phys.computed_metrics["cardiac_phase"].shape == (
+ n_scans,
+ slice_timings.size,
+ )
+ assert phys.computed_metrics["cardiac_phase"].args["slice_timings"] is not None
+ assert phys.computed_metrics["cardiac_phase"].args["n_scans"] is not None
+ assert phys.computed_metrics["cardiac_phase"].args["t_r"] is not None
+
+ # Test where the metric is returned
+ card_phase = cardiac.cardiac_phase(
+ phys,
+ slice_timings=slice_timings,
+ n_scans=n_scans,
+ t_r=t_r,
+ return_physio=False,
+ )
+ assert card_phase.ndim == 2
+ assert card_phase.shape == (n_scans, slice_timings.size)
diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py
index 1250d85..a4257ec 100644
--- a/phys2denoise/tests/test_metrics_chest_belt.py
+++ b/phys2denoise/tests/test_metrics_chest_belt.py
@@ -14,17 +14,14 @@ def test_rrf_smoke():
tr = 0.72
rrf_arr = chest_belt.rrf(
samplerate,
- oversampling=oversampling,
time_length=time_length,
onset=onset,
- tr=tr,
)
- pred_len = int(np.rint(time_length / (tr / oversampling)))
+ pred_len = int(np.rint(time_length * (1 / samplerate)))
assert rrf_arr.ndim == 1
assert rrf_arr.size == pred_len
-@mark.xfail
def test_respiratory_phase_smoke():
"""Basic smoke test for respiratory phase calculation."""
t_r = 1.0
@@ -35,7 +32,7 @@ def test_respiratory_phase_smoke():
resp = np.random.normal(size=n_samples)
resp_phase = chest_belt.respiratory_phase(
resp,
- sample_rate=sample_rate,
+ fs=sample_rate,
slice_timings=slice_timings,
n_scans=n_scans,
t_r=t_r,
@@ -59,7 +56,7 @@ def test_env_smoke():
resp = np.random.normal(size=n_samples)
samplerate = 1 / 0.01
window = 6
- env_arr = chest_belt.env(resp, samplerate=samplerate, window=window)
+ env_arr = chest_belt.env(resp, fs=samplerate, window=window)
assert env_arr.ndim == 1
assert env_arr.shape == (n_samples,)
@@ -70,6 +67,6 @@ def test_respiratory_variance_smoke():
resp = np.random.normal(size=n_samples)
samplerate = 1 / 0.01
window = 6
- rv_arr = chest_belt.respiratory_variance(resp, samplerate=samplerate, window=window)
+ rv_arr = chest_belt.respiratory_variance(resp, fs=samplerate, window=window)
assert rv_arr.ndim == 2
assert rv_arr.shape == (n_samples, 2)
diff --git a/phys2denoise/tests/test_metrics_utils.py b/phys2denoise/tests/test_metrics_utils.py
index 916c485..1995aaf 100644
--- a/phys2denoise/tests/test_metrics_utils.py
+++ b/phys2denoise/tests/test_metrics_utils.py
@@ -16,12 +16,12 @@ def test_mirrorpad(short_arr):
assert all(arr_mirror == expected_arr_mirror)
-def test_mirrorpad_exception(short_arr):
+def test_mirrorpad_exception(short_arr, caplog):
"""When passing array that is too short to perform mirrorpadding, the
function should give an error."""
arr = np.array(short_arr)
- with pytest.raises(Exception) as e_info:
- arr_mirror = mirrorpad_1d(short_arr)
+ arr_mirror = mirrorpad_1d(arr)
+ assert caplog.text.count("Fixing buffer size to array length.") > 0
def test_rms_envelope():
diff --git a/phys2denoise/tests/test_rvt.py b/phys2denoise/tests/test_rvt.py
index 13e27d8..0bc4066 100644
--- a/phys2denoise/tests/test_rvt.py
+++ b/phys2denoise/tests/test_rvt.py
@@ -1,4 +1,5 @@
import peakdet
+from physutils import io, physio
from phys2denoise.metrics.chest_belt import respiratory_variance_time
@@ -15,8 +16,49 @@ def test_respiratory_variance_time(fake_phys):
phys = peakdet.Physio(fake_phys, fs=62.5)
phys = peakdet.operations.filter_physio(phys, cutoffs=3, method="lowpass")
phys = peakdet.operations.peakfind_physio(phys)
- r = test_respiratory_variance_time(
- phys.data, phys.peaks, phys.troughs, samplerate=phys.fs
+
+ # TODO: Change to a simpler call once physutils are
+ # integrated to peakdet/prep4phys
+ r = respiratory_variance_time(
+ phys.data, fs=phys.fs, peaks=phys.peaks, troughs=phys.troughs
+ )
+ assert r is not None
+ assert len(r) == 18750
+
+
+def test_respiratory_variance_time(fake_phys):
+ phys = peakdet.Physio(fake_phys, fs=62.5)
+ phys = peakdet.operations.filter_physio(phys, cutoffs=3, method="lowpass")
+ phys = peakdet.operations.peakfind_physio(phys)
+
+ # TODO: Change to a simpler call once physutils are
+ # integrated to peakdet/prep4phys
+ r = respiratory_variance_time(
+ phys.data, fs=phys.fs, peaks=phys.peaks, troughs=phys.troughs
)
assert r is not None
assert len(r) == 18750
+
+
+def test_respiratory_variance_time_physio_obj(fake_phys):
+ phys = peakdet.Physio(fake_phys, fs=62.5)
+ phys = peakdet.operations.filter_physio(phys, cutoffs=3, method="lowpass")
+ phys = peakdet.operations.peakfind_physio(phys)
+
+ # TODO: Change to a simpler call once physutils are
+ # integrated to peakdet/prep4phys
+ physio_obj = physio.Physio(phys.data, phys.fs)
+ physio_obj._metadata["peaks"] = phys.peaks
+ physio_obj._metadata["troughs"] = phys.troughs
+ physio_obj = respiratory_variance_time(physio_obj)
+
+ assert (
+ physio_obj.history[0][0]
+ == "phys2denoise.metrics.chest_belt.respiratory_variance_time"
+ )
+ assert physio_obj.computed_metrics["respiratory_variance_time"].ndim == 2
+ assert physio_obj.computed_metrics["respiratory_variance_time"].shape == (18750, 4)
+ assert physio_obj.computed_metrics["respiratory_variance_time"].has_lags == True
+
+ # assert r is not None
+ # assert len(r) == 18750
diff --git a/setup.cfg b/setup.cfg
index a3f8ff9..81268a1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -19,13 +19,15 @@ provides =
phys2denoise
[options]
-python_requires = >=3.6.1
+python_requires = >=3.9
install_requires =
- numpy >=1.9.3
- matplotlib >=3.1.1
+ numpy >=1.9.3, <2
+ matplotlib
pandas
scipy
duecredit
+ loguru
+ physutils >=0.2.1
tests_require =
pytest >=5.3
test_suite = pytest
@@ -48,7 +50,7 @@ test =
%(style)s
pytest >=5.3
pytest-cov
- peakdet
+ peakdet>=0.5.0
coverage
devtools =
pre-commit