From 25360863611726b376664a631ce06f7ae2d0b55d Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Wed, 12 Jun 2024 18:53:52 +0300
Subject: [PATCH 01/38] Add physutils dependency placeholder

---
 setup.cfg | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/setup.cfg b/setup.cfg
index 8b89352..e460ed3 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -27,6 +27,8 @@ install_requires =
     pandas
     scipy
     duecredit
+    # TODO: Uncomment this once it's uploaded to PyPI
+    # physutils
 tests_require =
     pytest >=5.3
 test_suite = pytest

From f047587b2327796abe7a2c931661359354f36501 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Wed, 12 Jun 2024 20:36:36 +0300
Subject: [PATCH 02/38] Initial Physio object integration in cardiac.py

---
 phys2denoise/metrics/cardiac.py | 34 +++++++++++++++++++--------------
 1 file changed, 20 insertions(+), 14 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 40ebc47..9c933c9 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -6,9 +6,10 @@
 from .responses import crf
 from .utils import apply_function_in_sliding_window as afsw
 from .utils import convolve_and_rescale
+from physutils import physio
 
 
-def _cardiac_metrics(card, peaks, samplerate, metric, window=6, central_measure="mean"):
+def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
     """
     Compute cardiac metrics.
 
@@ -71,8 +72,11 @@ 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.
     """
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=True)
+
     # 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 +89,14 @@ 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 +113,13 @@ 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
 
 
 @due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
-def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
+def heart_rate(data, window=6, central_measure="mean"):
     """
     Compute average heart rate (HR) in a sliding window.
 
@@ -168,12 +172,12 @@ def heart_rate(card, peaks, samplerate, window=6, central_measure="mean"):
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
     return _cardiac_metrics(
-        card, peaks, samplerate, metric="hrv", window=6, central_measure="mean"
+        data, metric="hrv", window=6, central_measure="mean"
     )
 
 
 @due.dcite(references.PINHERO_ET_AL_2016)
-def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="mean"):
+def heart_rate_variability(data, window=6, central_measure="mean"):
     """
     Compute average heart rate variability (HRV) in a sliding window.
 
@@ -224,12 +228,12 @@ def heart_rate_variability(card, peaks, samplerate, window=6, central_measure="m
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
     return _cardiac_metrics(
-        card, peaks, samplerate, metric="hrv", window=6, central_measure="std"
+        data, metric="hrv", window=6, central_measure="std"
     )
 
 
 @due.dcite(references.CHEN_2020)
-def heart_beat_interval(card, peaks, samplerate, window=6, central_measure="mean"):
+def heart_beat_interval(data, window=6, central_measure="mean"):
     """
     Compute average heart beat interval (HBI) in a sliding window.
 
@@ -273,11 +277,11 @@ def heart_beat_interval(card, peaks, samplerate, window=6, central_measure="mean
         vol. 213, pp. 116707, 2020.
     """
     return _cardiac_metrics(
-        card, peaks, samplerate, metric="hbi", window=6, central_measure="mean"
+        data, metric="hbi", window=6, central_measure="mean"
     )
 
 
-def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
+def cardiac_phase(data, slice_timings, n_scans, t_r):
     """Calculate cardiac phase from cardiac peaks.
 
     Assumes that timing of cardiac events are given in same units
@@ -301,11 +305,13 @@ def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
     phase_card : array_like
         Cardiac phase signal, of shape (n_scans,)
     """
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=True)
     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]

From c7fc2eb373fe5120d56abf78c8b6d56e0b3c6c37 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Wed, 12 Jun 2024 20:37:57 +0300
Subject: [PATCH 03/38] Style fixes

---
 phys2denoise/metrics/cardiac.py | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 9c933c9..646d0b3 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -1,12 +1,12 @@
 """Denoising metrics for cardio recordings."""
 import numpy as np
+from physutils import 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
-from physutils import physio
 
 
 def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
@@ -96,7 +96,10 @@ def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
     card_met = np.empty_like(data)
     for n, i in enumerate(idx_min):
         diff = (
-            np.diff(data.peaks[np.logical_and(data.peaks >= i, data.peaks <= idx_max[n])]) / data.fs
+            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
@@ -171,9 +174,7 @@ def heart_rate(data, 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(
-        data, metric="hrv", window=6, central_measure="mean"
-    )
+    return _cardiac_metrics(data, metric="hrv", window=6, central_measure="mean")
 
 
 @due.dcite(references.PINHERO_ET_AL_2016)
@@ -227,9 +228,7 @@ def heart_rate_variability(data, 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(
-        data, metric="hrv", window=6, central_measure="std"
-    )
+    return _cardiac_metrics(data, metric="hrv", window=6, central_measure="std")
 
 
 @due.dcite(references.CHEN_2020)
@@ -276,9 +275,7 @@ def heart_beat_interval(data, window=6, central_measure="mean"):
     .. [1] J. E. Chen et al., "Resting-state "physiological networks"", Neuroimage,
         vol. 213, pp. 116707, 2020.
     """
-    return _cardiac_metrics(
-        data, metric="hbi", window=6, central_measure="mean"
-    )
+    return _cardiac_metrics(data, metric="hbi", window=6, central_measure="mean")
 
 
 def cardiac_phase(data, slice_timings, n_scans, t_r):

From ec8a3dedff409dca0ddd982cce443a375c37d207 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Wed, 12 Jun 2024 20:46:54 +0300
Subject: [PATCH 04/38] Minor cardiac.py physio initialization fix

---
 phys2denoise/metrics/cardiac.py | 42 ++++++++++-----------------------
 1 file changed, 12 insertions(+), 30 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 646d0b3..7e2e686 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -22,12 +22,8 @@ def _cardiac_metrics(data, metric, 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 : Physio_like
+        Object containing the timeseries of the recorded cardiac signal
     metrics : "hbi", "hr", "hrv", string
         Cardiac metric(s) to calculate.
     window : float, optional
@@ -73,7 +69,7 @@ def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
     # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True)
+    data = physio.check_physio(data, ensure_fs=True, copy=True)
 
     # Convert window to samples, but halves it.
     halfwindow_samples = int(round(window * data.fs / 2))
@@ -132,12 +128,8 @@ def heart_rate(data, 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 : Physio_like
+        Object containing the timeseries of the recorded cardiac signal
     window : float, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -188,12 +180,8 @@ def heart_rate_variability(data, 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 : Physio_like
+        Object containing the timeseries of the recorded cardiac signal
     window : float, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -238,12 +226,8 @@ def heart_beat_interval(data, 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 : Physio_like
+        Object containing the timeseries of the recorded cardiac signal
     window : float, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -286,10 +270,8 @@ def cardiac_phase(data, 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 : Physio_like
+        Object containing the timeseries of the recorded cardiac signal
     slice_timings : 1D array_like
         Slice times, in seconds.
     n_scans : int
@@ -303,7 +285,7 @@ def cardiac_phase(data, slice_timings, n_scans, t_r):
         Cardiac phase signal, of shape (n_scans,)
     """
     # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True)
+    data = physio.check_physio(data, ensure_fs=True, copy=True)
     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))

From 67c4d21d2c0446376d0621072d575ee2cc079cee Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Wed, 12 Jun 2024 21:08:51 +0300
Subject: [PATCH 05/38] Integrate Physio object for respiratory and multimodal
 metrics

---
 phys2denoise/metrics/chest_belt.py | 80 +++++++++++++++---------------
 phys2denoise/metrics/multimodal.py | 20 ++++----
 2 files changed, 49 insertions(+), 51 deletions(-)

diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index eeeb08a..5739276 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 physio
 from scipy.interpolate import interp1d
 from scipy.stats import zscore
 
@@ -12,7 +13,7 @@
 
 
 @due.dcite(references.BIRN_2006)
-def respiratory_variance_time(resp, peaks, troughs, samplerate, lags=(0, 4, 8, 12)):
+def respiratory_variance_time(data, lags=(0, 4, 8, 12)):
     """
     Implement the Respiratory Variance over Time (Birn et al. 2006).
 
@@ -20,14 +21,8 @@ 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 : Physio_like
+        Object containing the timeseries of the recorded respiratory signal
     lags: tuple
         lags in seconds of the RVT output. Default is 0, 4, 8, 12.
 
@@ -42,13 +37,15 @@ 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
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=True)
+    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
@@ -123,15 +120,13 @@ def respiratory_pattern_variability(resp, window):
 
 
 @due.dcite(references.POWER_2020)
-def env(resp, samplerate, window=10):
+def env(data, window=10):
     """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 : Physio_like
+        Object containing the timeseries of the recorded respiratory signal
     window : :obj:`int`, optional
         Size of the sliding window, in seconds.
         Default is 10.
@@ -155,13 +150,16 @@ def env(resp, samplerate, window=10):
        young adults scanned at rest, including systematic changes and 'missed'
        deep breaths," Neuroimage, vol. 204, 2020.
     """
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=True, copy=True)
+
     # 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)
         .rolling(window=window, center=True)
         .apply(respiratory_pattern_variability, args=(window,))
     )
@@ -170,15 +168,13 @@ def env(resp, samplerate, window=10):
 
 
 @due.dcite(references.CHANG_GLOVER_2009)
-def respiratory_variance(resp, samplerate, window=6):
+def respiratory_variance(data, window=6):
     """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 : Physio_like
+        Object containing the timeseries of the recorded respiratory signal
     window : :obj:`int`, optional
         Size of the sliding window, in seconds.
         Default is 6.
@@ -206,27 +202,28 @@ 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.
     """
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=True, copy=True)
+
     # 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
 
 
-def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
+def respiratory_phase(data, n_scans, slice_timings, t_r):
     """Calculate respiratory phase from respiratory signal.
 
     Parameters
     ----------
-    resp : 1D array_like
-        Respiratory signal.
-    sample_rate : float
-        Sample rate of physio, in Hertz.
+    data : Physio_like
+        Object containing the timeseries of the recorded respiratory signal
     n_scans
         Number of volumes in the imaging run.
     slice_timings
@@ -239,16 +236,19 @@ def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
     phase_resp : array_like
         Respiratory phase signal, of shape (n_scans, n_slices).
     """
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=True, copy=True)
+
     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,13 +256,13 @@ 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
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index 7a298b2..3717f78 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -1,6 +1,7 @@
 """These functions compute RETROICOR regressors (Glover et al. 2000)."""
 
 import numpy as np
+from physutils import physio
 
 from .. import references
 from ..due import due
@@ -10,8 +11,7 @@
 
 @due.dcite(references.GLOVER_2000)
 def retroicor(
-    physio,
-    sample_rate,
+    data,
     t_r,
     n_scans,
     slice_timings,
@@ -23,12 +23,8 @@ def retroicor(
 
     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 : Physio_like
+        Object containing the timeseries of the recorded respiratory or cardiac signal
     t_r : float
         Imaging sample rate, in seconds.
     n_scans : int
@@ -59,6 +55,9 @@ def retroicor(
        correction of physiological motion effects in fMRI: RETROICOR“, Magn. Reson. Med.,
        issue 1, vol. 44, pp. 162-167, 2000.
     """
+    # Initialize physio object
+    data = physio.check_physio(data, ensure_fs=True, copy=True)
+
     n_slices = np.shape(slice_timings)  # number of slices
 
     # initialize output variables
@@ -75,7 +74,7 @@ def retroicor(
         # slice sampling times
         if card:
             phase[:, i_slice] = cardiac_phase(
-                physio,
+                data,
                 crslice_timings,
                 n_scans,
                 t_r,
@@ -83,8 +82,7 @@ def retroicor(
 
         if resp:
             phase[:, i_slice] = respiratory_phase(
-                physio,
-                sample_rate,
+                data,
                 n_scans,
                 slice_timings,
                 t_r,

From 2ffa46597ad1cc779c03e4d08e75971105fdb8ad Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Thu, 13 Jun 2024 02:26:02 +0300
Subject: [PATCH 06/38] Add sample support for object and function oriented
 usage of metric functions

---
 phys2denoise/metrics/cardiac.py | 30 ++++++++++++++++++++++++++----
 1 file changed, 26 insertions(+), 4 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 7e2e686..8fb09aa 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -1,6 +1,6 @@
 """Denoising metrics for cardio recordings."""
 import numpy as np
-from physutils import physio
+from physutils import io, physio
 
 from .. import references
 from ..due import due
@@ -262,7 +262,7 @@ def heart_beat_interval(data, window=6, central_measure="mean"):
     return _cardiac_metrics(data, metric="hbi", window=6, central_measure="mean")
 
 
-def cardiac_phase(data, slice_timings, n_scans, t_r):
+def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
     """Calculate cardiac phase from cardiac peaks.
 
     Assumes that timing of cardiac events are given in same units
@@ -284,8 +284,30 @@ def cardiac_phase(data, slice_timings, n_scans, t_r):
     phase_card : array_like
         Cardiac phase signal, of shape (n_scans,)
     """
-    # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True, copy=True)
+    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 not data.peaks:
+        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))

From 9fa3467cd30e88a37a3688c067bdf52c747cb2e6 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 16:08:03 +0300
Subject: [PATCH 07/38] Make metric functions operations

---
 phys2denoise/metrics/cardiac.py    | 4 ++++
 phys2denoise/metrics/chest_belt.py | 5 +++++
 phys2denoise/metrics/multimodal.py | 1 +
 3 files changed, 10 insertions(+)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 8fb09aa..428151d 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -118,6 +118,7 @@ def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
 
 
 @due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
+@physio.make_operation()
 def heart_rate(data, window=6, central_measure="mean"):
     """
     Compute average heart rate (HR) in a sliding window.
@@ -170,6 +171,7 @@ def heart_rate(data, window=6, central_measure="mean"):
 
 
 @due.dcite(references.PINHERO_ET_AL_2016)
+@physio.make_operation()
 def heart_rate_variability(data, window=6, central_measure="mean"):
     """
     Compute average heart rate variability (HRV) in a sliding window.
@@ -220,6 +222,7 @@ def heart_rate_variability(data, window=6, central_measure="mean"):
 
 
 @due.dcite(references.CHEN_2020)
+@physio.make_operation()
 def heart_beat_interval(data, window=6, central_measure="mean"):
     """
     Compute average heart beat interval (HBI) in a sliding window.
@@ -262,6 +265,7 @@ def heart_beat_interval(data, window=6, central_measure="mean"):
     return _cardiac_metrics(data, metric="hbi", window=6, central_measure="mean")
 
 
+@physio.make_operation()
 def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
     """Calculate cardiac phase from cardiac peaks.
 
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 5739276..79467bd 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -13,6 +13,7 @@
 
 
 @due.dcite(references.BIRN_2006)
+@physio.make_operation()
 def respiratory_variance_time(data, lags=(0, 4, 8, 12)):
     """
     Implement the Respiratory Variance over Time (Birn et al. 2006).
@@ -78,6 +79,7 @@ def respiratory_variance_time(data, lags=(0, 4, 8, 12)):
 
 
 @due.dcite(references.POWER_2018)
+@physio.make_operation()
 def respiratory_pattern_variability(resp, window):
     """Calculate respiratory pattern variability.
 
@@ -120,6 +122,7 @@ def respiratory_pattern_variability(resp, window):
 
 
 @due.dcite(references.POWER_2020)
+@physio.make_operation()
 def env(data, window=10):
     """Calculate respiratory pattern variability across a sliding window.
 
@@ -168,6 +171,7 @@ def env(data, window=10):
 
 
 @due.dcite(references.CHANG_GLOVER_2009)
+@physio.make_operation()
 def respiratory_variance(data, window=6):
     """Calculate respiratory variance.
 
@@ -217,6 +221,7 @@ def respiratory_variance(data, window=6):
     return rv_out
 
 
+@physio.make_operation()
 def respiratory_phase(data, n_scans, slice_timings, t_r):
     """Calculate respiratory phase from respiratory signal.
 
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index 3717f78..f8c6968 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -10,6 +10,7 @@
 
 
 @due.dcite(references.GLOVER_2000)
+@physio.make_operation()
 def retroicor(
     data,
     t_r,

From 380ab65739570bed7f5d9cfe91933a792a00ecb7 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 16:20:28 +0300
Subject: [PATCH 08/38] Add Physio and non-Physio compatibility to cardiac.py

---
 phys2denoise/metrics/cardiac.py | 67 ++++++++++++++++++++++++++++-----
 1 file changed, 58 insertions(+), 9 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 428151d..00fe312 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -9,7 +9,9 @@
 from .utils import convolve_and_rescale
 
 
-def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
+def _cardiac_metrics(
+    data, metric, fs=None, peaks=None, window=6, central_measure="mean"
+):
     """
     Compute cardiac metrics.
 
@@ -26,6 +28,10 @@ def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
         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.
@@ -68,8 +74,29 @@ def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
     Annual International Conference of the IEEE Engineering in Medicine and
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
-    # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True, copy=True)
+    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 not data.peaks:
+        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 * data.fs / 2))
@@ -119,7 +146,7 @@ def _cardiac_metrics(data, metric, window=6, central_measure="mean"):
 
 @due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
 @physio.make_operation()
-def heart_rate(data, window=6, central_measure="mean"):
+def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean"):
     """
     Compute average heart rate (HR) in a sliding window.
 
@@ -131,6 +158,10 @@ def heart_rate(data, window=6, central_measure="mean"):
     ----------
     data : Physio_like
         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,12 +198,14 @@ def heart_rate(data, 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(data, metric="hrv", window=6, central_measure="mean")
+    return _cardiac_metrics(
+        data, metric="hrv", fs=fs, peaks=peaks, window=6, central_measure="mean"
+    )
 
 
 @due.dcite(references.PINHERO_ET_AL_2016)
 @physio.make_operation()
-def heart_rate_variability(data, window=6, central_measure="mean"):
+def heart_rate_variability(data, fs=None, peaks=None, window=6, central_measure="mean"):
     """
     Compute average heart rate variability (HRV) in a sliding window.
 
@@ -184,6 +217,10 @@ def heart_rate_variability(data, window=6, central_measure="mean"):
     ----------
     data : Physio_like
         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.
@@ -218,12 +255,14 @@ def heart_rate_variability(data, 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(data, metric="hrv", window=6, central_measure="std")
+    return _cardiac_metrics(
+        data, metric="hrv", fs=None, peaks=None, window=6, central_measure="std"
+    )
 
 
 @due.dcite(references.CHEN_2020)
 @physio.make_operation()
-def heart_beat_interval(data, window=6, central_measure="mean"):
+def heart_beat_interval(data, fs=None, peaks=None, window=6, central_measure="mean"):
     """
     Compute average heart beat interval (HBI) in a sliding window.
 
@@ -231,6 +270,10 @@ def heart_beat_interval(data, window=6, central_measure="mean"):
     ----------
     data : Physio_like
         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.
@@ -262,7 +305,9 @@ def heart_beat_interval(data, window=6, central_measure="mean"):
     .. [1] J. E. Chen et al., "Resting-state "physiological networks"", Neuroimage,
         vol. 213, pp. 116707, 2020.
     """
-    return _cardiac_metrics(data, metric="hbi", window=6, central_measure="mean")
+    return _cardiac_metrics(
+        data, metric="hbi", fs=None, peaks=None, window=6, central_measure="mean"
+    )
 
 
 @physio.make_operation()
@@ -282,6 +327,10 @@ def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
         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
     -------

From 59ee6c8aa11b2eeb43dfa4ea7f37d2b65dbe693f Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 18:27:55 +0300
Subject: [PATCH 09/38] Add Physio and non-Physio compatibility to
 chest_belt.py

---
 phys2denoise/metrics/chest_belt.py | 110 ++++++++++++++++++++++++-----
 1 file changed, 93 insertions(+), 17 deletions(-)

diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 79467bd..e48865a 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -1,7 +1,7 @@
 """Denoising metrics for chest belt recordings."""
 import numpy as np
 import pandas as pd
-from physutils import physio
+from physutils import io, physio
 from scipy.interpolate import interp1d
 from scipy.stats import zscore
 
@@ -14,7 +14,9 @@
 
 @due.dcite(references.BIRN_2006)
 @physio.make_operation()
-def respiratory_variance_time(data, lags=(0, 4, 8, 12)):
+def respiratory_variance_time(
+    data, fs=None, peaks=None, troughs=None, lags=(0, 4, 8, 12)
+):
     """
     Implement the Respiratory Variance over Time (Birn et al. 2006).
 
@@ -24,6 +26,12 @@ def respiratory_variance_time(data, lags=(0, 4, 8, 12)):
     ----------
     data : Physio_like
         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.
 
@@ -38,8 +46,31 @@ def respiratory_variance_time(data, lags=(0, 4, 8, 12)):
        respiratory-variation-related fluctuations from neuronal-activity-related
        fluctuations in fMRI”, NeuroImage, vol.31, pp. 1536-1548, 2006.
     """
-    # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True)
+    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 not data.peaks or not data.troughs:
+        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(data) * timestep, timestep)
@@ -80,13 +111,13 @@ def respiratory_variance_time(data, lags=(0, 4, 8, 12)):
 
 @due.dcite(references.POWER_2018)
 @physio.make_operation()
-def respiratory_pattern_variability(resp, window):
+def respiratory_pattern_variability(data, window):
     """Calculate respiratory pattern variability.
 
     Parameters
     ----------
-    resp : str or 1D numpy.ndarray
-        Tiemseries representing respiration activity.
+    data : Physio_like
+        Object containing the timeseries of the recorded respiratory signal
     window : int
         Window length in samples.
 
@@ -110,8 +141,11 @@ 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=True, 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)
@@ -123,13 +157,15 @@ def respiratory_pattern_variability(resp, window):
 
 @due.dcite(references.POWER_2020)
 @physio.make_operation()
-def env(data, window=10):
+def env(data, fs=None, window=10):
     """Calculate respiratory pattern variability across a sliding window.
 
     Parameters
     ----------
     data : Physio_like
         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.
@@ -153,8 +189,20 @@ def env(data, window=10):
        young adults scanned at rest, including systematic changes and 'missed'
        deep breaths," Neuroimage, vol. 204, 2020.
     """
-    # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True, copy=True)
+    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 * data.fs)
@@ -172,13 +220,15 @@ def env(data, window=10):
 
 @due.dcite(references.CHANG_GLOVER_2009)
 @physio.make_operation()
-def respiratory_variance(data, window=6):
+def respiratory_variance(data, fs=None, window=6):
     """Calculate respiratory variance.
 
     Parameters
     ----------
     data : Physio_like
         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,8 +256,20 @@ def respiratory_variance(data, window=6):
        end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
        issue 4, vol. 47, pp. 1381-1393, 2009.
     """
-    # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True, copy=True)
+    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 * data.fs / 2))
@@ -222,7 +284,7 @@ def respiratory_variance(data, window=6):
 
 
 @physio.make_operation()
-def respiratory_phase(data, n_scans, slice_timings, t_r):
+def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None):
     """Calculate respiratory phase from respiratory signal.
 
     Parameters
@@ -235,14 +297,28 @@ def respiratory_phase(data, n_scans, slice_timings, t_r):
         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).
     """
-    # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True, copy=True)
+    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)

From 56abb720a83ef9fb969b32402dcdd91750d2b19a Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 18:56:55 +0300
Subject: [PATCH 10/38] Add Physio and non-Physio compatibility to retroicor

---
 phys2denoise/metrics/multimodal.py | 50 +++++++++++++++++++++++++-----
 1 file changed, 43 insertions(+), 7 deletions(-)

diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index f8c6968..0f3f076 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -1,7 +1,7 @@
 """These functions compute RETROICOR regressors (Glover et al. 2000)."""
 
 import numpy as np
-from physutils import physio
+from physutils import io, physio
 
 from .. import references
 from ..due import due
@@ -17,8 +17,9 @@ def retroicor(
     n_scans,
     slice_timings,
     n_harmonics,
-    card=False,
-    resp=False,
+    physio_type=None,
+    fs=None,
+    cardiac_peaks=None,
 ):
     """Compute RETROICOR regressors.
 
@@ -56,8 +57,43 @@ def retroicor(
        correction of physiological motion effects in fMRI: RETROICOR“, Magn. Reson. Med.,
        issue 1, vol. 44, pp. 162-167, 2000.
     """
-    # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True, copy=True)
+    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
 
@@ -73,7 +109,7 @@ 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(
                 data,
                 crslice_timings,
@@ -81,7 +117,7 @@ def retroicor(
                 t_r,
             )
 
-        if resp:
+        if data.physio_type == "respiratory":
             phase[:, i_slice] = respiratory_phase(
                 data,
                 n_scans,

From 2c19a65a5b0b7795994e86ce127374b5634e1a7a Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 22:26:00 +0300
Subject: [PATCH 11/38] Fix CRF, iCRF and RRF calculation

---
 phys2denoise/metrics/responses.py          |  7 +++++--
 phys2denoise/tests/conftest.py             | 15 +++++++++++++++
 phys2denoise/tests/test_metrics_cardiac.py | 11 +++--------
 3 files changed, 23 insertions(+), 10 deletions(-)

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/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..2206030 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -7,20 +7,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():

From c28914c5eda8f0b9187a58daa58481d3e8a80113 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 23:09:33 +0300
Subject: [PATCH 12/38] Fix cardiac_phase test, return Physio object and metric
 to update history

---
 phys2denoise/metrics/cardiac.py            | 9 +++++----
 phys2denoise/tests/test_metrics_cardiac.py | 9 ++++++---
 2 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 00fe312..3e217d7 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -1,5 +1,6 @@
 """Denoising metrics for cardio recordings."""
 import numpy as np
+from loguru import logger
 from physutils import io, physio
 
 from .. import references
@@ -89,7 +90,7 @@ def _cardiac_metrics(
             and the peak indices separately.
             """
         )
-    if not data.peaks:
+    if data.peaks.size == 0:
         raise ValueError(
             """
             Peaks must be a non-empty list.
@@ -141,7 +142,7 @@ def _cardiac_metrics(
     # Convolve with crf and 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)
@@ -352,7 +353,7 @@ def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
             and the peak indices separately.
             """
         )
-    if not data.peaks:
+    if data.peaks.size == 0:
         raise ValueError(
             """
             Peaks must be a non-empty list.
@@ -392,4 +393,4 @@ def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
             ) / (t2 - t1)
         phase_card[:, i_slice] = phase_card_crSlice
 
-    return phase_card
+    return data, phase_card
diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index 2206030..cf125e2 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -1,5 +1,6 @@
 """Tests for phys2denoise.metrics.cardiac."""
 import numpy as np
+from loguru import logger
 
 from phys2denoise.metrics import cardiac
 
@@ -25,9 +26,11 @@ 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])
-    card_phase = cardiac.cardiac_phase(
-        peaks,
-        sample_rate=sample_rate,
+    data = np.zeros(peaks.shape)
+    _, card_phase = cardiac.cardiac_phase(
+        data,
+        peaks=peaks,
+        fs=sample_rate,
         slice_timings=slice_timings,
         n_scans=n_scans,
         t_r=t_r,

From 625d45f0e9b149c9af31489f02d538772178bbd2 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 23:40:31 +0300
Subject: [PATCH 13/38] Fix chest_belt tests, return Physio object and metric
 to update history with each operation

---
 phys2denoise/metrics/chest_belt.py            | 35 ++++++++++++++-----
 phys2denoise/tests/test_metrics_chest_belt.py | 12 +++----
 2 files changed, 32 insertions(+), 15 deletions(-)

diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index e48865a..e29eccb 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -62,7 +62,7 @@ def respiratory_variance_time(
             and the peak and trough indices separately.
             """
         )
-    if not data.peaks or not data.troughs:
+    if data.peaks.size == 0 or data.troughs.size == 0:
         raise ValueError(
             """
             Peaks and troughs must be non-empty lists.
@@ -106,7 +106,7 @@ def respiratory_variance_time(
         )
         rvt_lags[:, ind] = temp_rvt
 
-    return rvt_lags
+    return data, rvt_lags
 
 
 @due.dcite(references.POWER_2018)
@@ -142,7 +142,7 @@ def respiratory_pattern_variability(data, window):
        115, pp. 2105-2114, 2018.
     """
     # Initialize physio object
-    data = physio.check_physio(data, ensure_fs=True, copy=True)
+    data = physio.check_physio(data, ensure_fs=False, copy=True)
 
     # First, z-score respiratory traces
     resp_z = zscore(data.data)
@@ -152,7 +152,7 @@ def respiratory_pattern_variability(data, window):
 
     # Calculate standard deviation
     rpv_val = np.std(rpv_upper_env)
-    return rpv_val
+    return data, rpv_val
 
 
 @due.dcite(references.POWER_2020)
@@ -189,6 +189,25 @@ def env(data, fs=None, 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)
@@ -210,12 +229,12 @@ def env(data, fs=None, window=10):
     # Calculate RPV across a rolling window
 
     env_arr = (
-        pd.Series(data)
+        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)
@@ -280,7 +299,7 @@ def respiratory_variance(data, fs=None, window=6):
     # Convolve with rrf
     rv_out = convolve_and_rescale(rv_arr, rrf(data.fs), rescale="zscore")
 
-    return rv_out
+    return data, rv_out
 
 
 @physio.make_operation()
diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py
index 1250d85..27179e4 100644
--- a/phys2denoise/tests/test_metrics_chest_belt.py
+++ b/phys2denoise/tests/test_metrics_chest_belt.py
@@ -14,12 +14,10 @@ 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
 
@@ -33,7 +31,7 @@ def test_respiratory_phase_smoke():
     slice_timings = np.linspace(0, t_r, 22)[1:-1]
     n_samples = int(np.rint((n_scans * t_r) * sample_rate))
     resp = np.random.normal(size=n_samples)
-    resp_phase = chest_belt.respiratory_phase(
+    _, resp_phase = chest_belt.respiratory_phase(
         resp,
         sample_rate=sample_rate,
         slice_timings=slice_timings,
@@ -49,7 +47,7 @@ def test_respiratory_pattern_variability_smoke():
     n_samples = 2000
     resp = np.random.normal(size=n_samples)
     window = 50
-    rpv_val = chest_belt.respiratory_pattern_variability(resp, window)
+    _, rpv_val = chest_belt.respiratory_pattern_variability(resp, window)
     assert isinstance(rpv_val, float)
 
 
@@ -59,7 +57,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 +68,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)

From acf960d792ae64074dae72d680d7bfcd1dfad834 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 24 Jun 2024 23:57:00 +0300
Subject: [PATCH 14/38] Update RVT test

---
 phys2denoise/metrics/chest_belt.py | 2 +-
 phys2denoise/tests/test_rvt.py     | 7 +++++--
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index e29eccb..4e83b7d 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -367,4 +367,4 @@ def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None):
 
         phase_resp[:, i_slice] = phase_resp_crSlice
 
-    return phase_resp
+    return data, phase_resp
diff --git a/phys2denoise/tests/test_rvt.py b/phys2denoise/tests/test_rvt.py
index 13e27d8..67a7bf9 100644
--- a/phys2denoise/tests/test_rvt.py
+++ b/phys2denoise/tests/test_rvt.py
@@ -15,8 +15,11 @@ 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
+    phys, r = respiratory_variance_time(
+        phys.data, fs=phys.fs, peaks=phys.peaks, troughs=phys.troughs
     )
     assert r is not None
     assert len(r) == 18750

From 62e4b6f93fb0c6073736313469cac8a6be6ea0ce Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Tue, 25 Jun 2024 00:00:17 +0300
Subject: [PATCH 15/38] Fix respiratory_phase_smoke test

---
 phys2denoise/tests/test_metrics_chest_belt.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py
index 27179e4..5e5bfd1 100644
--- a/phys2denoise/tests/test_metrics_chest_belt.py
+++ b/phys2denoise/tests/test_metrics_chest_belt.py
@@ -22,7 +22,6 @@ def test_rrf_smoke():
     assert rrf_arr.size == pred_len
 
 
-@mark.xfail
 def test_respiratory_phase_smoke():
     """Basic smoke test for respiratory phase calculation."""
     t_r = 1.0
@@ -33,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,

From e60d49ec41ef0e92b09869a2b28a33380c3fc626 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Tue, 25 Jun 2024 00:07:59 +0300
Subject: [PATCH 16/38] Return physio in retroicor and specify raised error
 type in test_mirrorpad_exception

---
 phys2denoise/metrics/multimodal.py       | 2 +-
 phys2denoise/tests/test_metrics_utils.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index 0f3f076..f5b479e 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -134,4 +134,4 @@ def retroicor(
                 (j_harm + 1) * phase[i_slice]
             )
 
-    return retroicor_regressors, phase
+    return data, retroicor_regressors, phase
diff --git a/phys2denoise/tests/test_metrics_utils.py b/phys2denoise/tests/test_metrics_utils.py
index 916c485..d7e1b46 100644
--- a/phys2denoise/tests/test_metrics_utils.py
+++ b/phys2denoise/tests/test_metrics_utils.py
@@ -20,7 +20,7 @@ def test_mirrorpad_exception(short_arr):
     """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:
+    with pytest.raises(IndexError):
         arr_mirror = mirrorpad_1d(short_arr)
 
 

From 3dbf65612629873114509752e924cbf569646894 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Tue, 25 Jun 2024 00:13:05 +0300
Subject: [PATCH 17/38] Add loguru as a dependency

---
 setup.cfg | 1 +
 1 file changed, 1 insertion(+)

diff --git a/setup.cfg b/setup.cfg
index e460ed3..6f9e718 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -27,6 +27,7 @@ install_requires =
     pandas
     scipy
     duecredit
+    loguru
     # TODO: Uncomment this once it's uploaded to PyPI
     # physutils
 tests_require =

From f7d781116e9948b5403386bae037022a4d4104c5 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Tue, 25 Jun 2024 12:38:58 +0300
Subject: [PATCH 18/38] Add Physio object use cardiac phase unit test

---
 phys2denoise/tests/test_metrics_cardiac.py | 23 ++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index cf125e2..6a7be8c 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -1,6 +1,7 @@
 """Tests for phys2denoise.metrics.cardiac."""
 import numpy as np
 from loguru import logger
+from physutils import physio
 
 from phys2denoise.metrics import cardiac
 
@@ -37,3 +38,25 @@ def test_cardiac_phase_smoke():
     )
     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
+    phys, card_phase = cardiac.cardiac_phase(
+        phys,
+        slice_timings=slice_timings,
+        n_scans=n_scans,
+        t_r=t_r,
+    )
+    logger.debug(f"History: {phys.history}")
+    assert phys.history[0][0] == "phys2denoise.metrics.cardiac.cardiac_phase"
+    assert card_phase.ndim == 2
+    assert card_phase.shape == (n_scans, slice_timings.size)

From 907faba7d592bd5700ed1048687722daa09ef88c Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Tue, 25 Jun 2024 12:39:28 +0300
Subject: [PATCH 19/38] Remove test log

---
 phys2denoise/tests/test_metrics_cardiac.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index 6a7be8c..25faa9a 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -56,7 +56,6 @@ def test_cardiac_phase_smoke_physio_obj():
         n_scans=n_scans,
         t_r=t_r,
     )
-    logger.debug(f"History: {phys.history}")
     assert phys.history[0][0] == "phys2denoise.metrics.cardiac.cardiac_phase"
     assert card_phase.ndim == 2
     assert card_phase.shape == (n_scans, slice_timings.size)

From 2d9fd750197377c14afebda8394551f91f1c080a Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Thu, 18 Jul 2024 15:43:35 +0300
Subject: [PATCH 20/38] Add physutils dependency and peakdet version

---
 setup.cfg | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/setup.cfg b/setup.cfg
index 6f9e718..82a5696 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -28,8 +28,7 @@ install_requires =
     scipy
     duecredit
     loguru
-    # TODO: Uncomment this once it's uploaded to PyPI
-    # physutils
+    physutils
 tests_require =
     pytest >=5.3
 test_suite = pytest
@@ -52,7 +51,7 @@ test =
     %(style)s
     pytest >=5.3
     pytest-cov
-    peakdet
+    peakdet>=0.4.0
     coverage
 devtools =
     pre-commit

From 366175a5a153146f17dec605869046039264aab5 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sat, 20 Jul 2024 00:48:56 +0300
Subject: [PATCH 21/38] Update peakdet dependency

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index 82a5696..562bae4 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -51,7 +51,7 @@ test =
     %(style)s
     pytest >=5.3
     pytest-cov
-    peakdet>=0.4.0
+    peakdet>=0.5.0
     coverage
 devtools =
     pre-commit

From 533c80e055af0fd8de104f99cffeda46874abcef Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 22 Jul 2024 10:27:00 +0300
Subject: [PATCH 22/38] Minor fix

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index 562bae4..23813a9 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -23,7 +23,7 @@ provides =
 python_requires = >=3.6.1
 install_requires =
     numpy >=1.9.3
-    matplotlib >=3.1.1
+    matplotlib>=3.9
     pandas
     scipy
     duecredit

From 3fba78c1ab7752c7f267eee14960f5027af1c4f9 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 22 Jul 2024 10:44:55 +0300
Subject: [PATCH 23/38] circleci: python version upgrade

---
 .circleci/config.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index fc36d30..176c0a7 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -85,7 +85,7 @@ jobs:
 
   style_check:
     docker:
-      - image: cimg/python:3.7
+      - image: cimg/python:3.10
     working_directory: /tmp/src/phys2denoise
     resource_class: small
     steps:

From 43576d0edee92c1678a2c77222fd0630ae431684 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 22 Jul 2024 11:22:25 +0300
Subject: [PATCH 24/38] circleci: python version migration to 3.8 and 3.12

---
 .circleci/config.yml | 24 ++++++++++++------------
 setup.cfg            |  2 +-
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 176c0a7..071888b 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!
+  test38: # 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.8
     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.py38
       - 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.py38
 
-  test310:
+  test312:
     docker:
-      - image: cimg/python:3.10
+      - image: cimg/python:3.12
     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.py312
       - store_artifacts:
           path: /tmp/src/coverage
       - persist_to_workspace:
           root: /tmp
           paths:
-            - src/coverage/.coverage.py310
+            - src/coverage/.coverage.py312
 
   style_check:
     docker:
-      - image: cimg/python:3.10
+      - image: cimg/python:3.8
     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.12
     resource_class: small
     steps:
       - attach_workspace:
@@ -141,5 +141,5 @@ workflows:
             - style_check
       - merge_coverage:
           requires:
-            - test37
-            - test310
+            - test38
+            - test312
diff --git a/setup.cfg b/setup.cfg
index 23813a9..7568f5d 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -20,7 +20,7 @@ provides =
     phys2denoise
 
 [options]
-python_requires = >=3.6.1
+python_requires = >=3.8
 install_requires =
     numpy >=1.9.3
     matplotlib>=3.9

From e2f22419bbde7ea5692bad68352318d4f75fb91e Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 22 Jul 2024 11:23:29 +0300
Subject: [PATCH 25/38] Minor CI fix

---
 .circleci/config.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 071888b..a392ad0 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -133,10 +133,10 @@ workflows:
     # Inside the workflow, you define the jobs you want to run.
     jobs:
       - style_check
-      - test37:
+      - test38:
           requires:
             - style_check
-      - test310:
+      - test312:
           requires:
             - style_check
       - merge_coverage:

From 9efbf5fe1065e9da7c15dbdec52bd20c970b4032 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 22 Jul 2024 11:43:17 +0300
Subject: [PATCH 26/38] CI: Change minimum python version to 3.9

---
 .circleci/config.yml | 14 +++++++-------
 setup.cfg            |  2 +-
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index a392ad0..5508e83 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:
-  test38: # 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.8
+      - 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.py38
+            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,7 +56,7 @@ jobs:
           root: /tmp
           # Must be relative path from root
           paths:
-            - src/coverage/.coverage.py38
+            - src/coverage/.coverage.py39
 
   test312:
     docker:
@@ -85,7 +85,7 @@ jobs:
 
   style_check:
     docker:
-      - image: cimg/python:3.8
+      - image: cimg/python:3.9
     working_directory: /tmp/src/phys2denoise
     resource_class: small
     steps:
@@ -133,7 +133,7 @@ workflows:
     # Inside the workflow, you define the jobs you want to run.
     jobs:
       - style_check
-      - test38:
+      - test39:
           requires:
             - style_check
       - test312:
@@ -141,5 +141,5 @@ workflows:
             - style_check
       - merge_coverage:
           requires:
-            - test38
+            - test39
             - test312
diff --git a/setup.cfg b/setup.cfg
index 7568f5d..11316bc 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -20,7 +20,7 @@ provides =
     phys2denoise
 
 [options]
-python_requires = >=3.8
+python_requires = >=3.9
 install_requires =
     numpy >=1.9.3
     matplotlib>=3.9

From ba5dcbc33db27cd67ce7ab5d6a698e378be44c01 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 22 Jul 2024 11:57:43 +0300
Subject: [PATCH 27/38] CI: Test for python 3.9 and 3.11

---
 .circleci/config.yml | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 5508e83..dd45175 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -58,9 +58,9 @@ jobs:
           paths:
             - src/coverage/.coverage.py39
 
-  test312:
+  test311:
     docker:
-      - image: cimg/python:3.12
+      - image: cimg/python:3.11
     working_directory: /tmp/src/phys2denoise
     resource_class: medium
     steps:
@@ -75,13 +75,13 @@ jobs:
           command: |
             pytest --cov=./phys2denoise
             mkdir /tmp/src/coverage
-            mv ./.coverage /tmp/src/coverage/.coverage.py312
+            mv ./.coverage /tmp/src/coverage/.coverage.py311
       - store_artifacts:
           path: /tmp/src/coverage
       - persist_to_workspace:
           root: /tmp
           paths:
-            - src/coverage/.coverage.py312
+            - src/coverage/.coverage.py311
 
   style_check:
     docker:
@@ -105,7 +105,7 @@ jobs:
   merge_coverage:
     working_directory: /tmp/src/phys2denoise
     docker:
-      - image: cimg/python:3.12
+      - image: cimg/python:3.11
     resource_class: small
     steps:
       - attach_workspace:
@@ -136,10 +136,10 @@ workflows:
       - test39:
           requires:
             - style_check
-      - test312:
+      - test311:
           requires:
             - style_check
       - merge_coverage:
           requires:
             - test39
-            - test312
+            - test311

From 9f5e1629fc3f501c33978594d55e7b5f57013bff Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Mon, 22 Jul 2024 19:57:13 +0300
Subject: [PATCH 28/38] Store computed metrics inside the Physio object

---
 phys2denoise/metrics/cardiac.py    | 19 +++++++++++++++----
 phys2denoise/metrics/chest_belt.py |  9 +++++++++
 phys2denoise/metrics/multimodal.py |  3 +++
 3 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 3e217d7..06d8b50 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -199,9 +199,12 @@ def heart_rate(data, fs=None, peaks=None, 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(
-        data, metric="hrv", fs=fs, peaks=peaks, window=6, central_measure="mean"
+    data, hr = _cardiac_metrics(
+        data, metric="hr", fs=fs, peaks=peaks, window=6, central_measure="mean"
     )
+    data._computed_metrics["heart_rate"] = dict(metric=hr, has_lags=False)
+
+    return data, hr
 
 
 @due.dcite(references.PINHERO_ET_AL_2016)
@@ -256,9 +259,12 @@ def heart_rate_variability(data, fs=None, peaks=None, window=6, central_measure=
     Annual International Conference of the IEEE Engineering in Medicine and
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
-    return _cardiac_metrics(
+    data, hrv = _cardiac_metrics(
         data, metric="hrv", fs=None, peaks=None, window=6, central_measure="std"
     )
+    data._computed_metrics["heart_rate_variability"] = dict(metric=hrv)
+
+    return data, hrv
 
 
 @due.dcite(references.CHEN_2020)
@@ -306,9 +312,12 @@ def heart_beat_interval(data, fs=None, peaks=None, window=6, central_measure="me
     .. [1] J. E. Chen et al., "Resting-state "physiological networks"", Neuroimage,
         vol. 213, pp. 116707, 2020.
     """
-    return _cardiac_metrics(
+    data, hbi = _cardiac_metrics(
         data, metric="hbi", fs=None, peaks=None, window=6, central_measure="mean"
     )
+    data._computed_metrics["heart_beat_interval"] = dict(metric=hbi)
+
+    return data, hbi
 
 
 @physio.make_operation()
@@ -393,4 +402,6 @@ def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
             ) / (t2 - t1)
         phase_card[:, i_slice] = phase_card_crSlice
 
+    data._computed_metrics["cardiac_phase"] = dict(metric=phase_card)
+
     return data, phase_card
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 4e83b7d..1c20aeb 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -106,6 +106,7 @@ def respiratory_variance_time(
         )
         rvt_lags[:, ind] = temp_rvt
 
+    data._computed_metrics["rvt"] = dict(metric=rvt_lags, has_lags=True)
     return data, rvt_lags
 
 
@@ -152,6 +153,8 @@ def respiratory_pattern_variability(data, window):
 
     # Calculate standard deviation
     rpv_val = np.std(rpv_upper_env)
+
+    data._computed_metrics["rpv"] = dict(metric=rpv_val)
     return data, rpv_val
 
 
@@ -234,6 +237,8 @@ def _respiratory_pattern_variability(data, window):
         .apply(_respiratory_pattern_variability, args=(window,))
     )
     env_arr[np.isnan(env_arr)] = 0.0
+
+    data._computed_metrics["env"] = dict(metric=env_arr)
     return data, env_arr
 
 
@@ -299,6 +304,8 @@ def respiratory_variance(data, fs=None, window=6):
     # Convolve with rrf
     rv_out = convolve_and_rescale(rv_arr, rrf(data.fs), rescale="zscore")
 
+    data._computed_metrics["respiratory_variance"] = dict(metric=rv_out)
+
     return data, rv_out
 
 
@@ -367,4 +374,6 @@ def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None):
 
         phase_resp[:, i_slice] = phase_resp_crSlice
 
+    data._computed_metrics["respiratory_phase"] = dict(metric=phase_resp)
+
     return data, phase_resp
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index f5b479e..fc48213 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -134,4 +134,7 @@ def retroicor(
                 (j_harm + 1) * phase[i_slice]
             )
 
+    data._computed_metrics["retroicor_regressors"] = dict(
+        metrics=retroicor_regressors, phase=phase
+    )
     return data, retroicor_regressors, phase

From daaed9c86adb650e29e0e5865f9410256786262d Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Tue, 23 Jul 2024 16:07:46 +0300
Subject: [PATCH 29/38] Update .all-contributorsrc

---
 .all-contributorsrc | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/.all-contributorsrc b/.all-contributorsrc
index 46f9a58..e8b53b1 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -110,6 +110,18 @@
         "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",
+        "bug",
+        "review"
+      ]
     }
   ],
   "contributorsPerLine": 8,

From bc300f046c69e8315668234bd82a4e66107512a8 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sat, 10 Aug 2024 20:40:18 +0300
Subject: [PATCH 30/38] Add wrapper to determine whether to return a metric or
 a physio object

---
 phys2denoise/metrics/cardiac.py               | 40 ++++++++++++++-----
 phys2denoise/metrics/chest_belt.py            | 14 +++----
 phys2denoise/metrics/utils.py                 | 38 ++++++++++++++++++
 phys2denoise/tests/test_metrics_cardiac.py    |  2 +-
 phys2denoise/tests/test_metrics_chest_belt.py |  8 ++--
 phys2denoise/tests/test_rvt.py                |  4 +-
 6 files changed, 81 insertions(+), 25 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 06d8b50..72c7e8d 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -7,11 +7,17 @@
 from ..due import due
 from .responses import crf
 from .utils import apply_function_in_sliding_window as afsw
-from .utils import convolve_and_rescale
+from .utils import convolve_and_rescale, return_physio_or_metric
 
 
 def _cardiac_metrics(
-    data, metric, fs=None, peaks=None, window=6, central_measure="mean"
+    data,
+    metric,
+    fs=None,
+    peaks=None,
+    window=6,
+    central_measure="mean",
+    return_physio=False,
 ):
     """
     Compute cardiac metrics.
@@ -146,6 +152,7 @@ def _cardiac_metrics(
 
 
 @due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
+@return_physio_or_metric()
 @physio.make_operation()
 def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean"):
     """
@@ -200,14 +207,19 @@ def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean"):
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
     data, hr = _cardiac_metrics(
-        data, metric="hr", fs=fs, peaks=peaks, window=6, central_measure="mean"
+        data,
+        metric="hr",
+        fs=fs,
+        peaks=peaks,
+        window=window,
+        central_measure=central_measure,
     )
-    data._computed_metrics["heart_rate"] = dict(metric=hr, has_lags=False)
 
     return data, hr
 
 
 @due.dcite(references.PINHERO_ET_AL_2016)
+@return_physio_or_metric()
 @physio.make_operation()
 def heart_rate_variability(data, fs=None, peaks=None, window=6, central_measure="mean"):
     """
@@ -260,14 +272,19 @@ def heart_rate_variability(data, fs=None, peaks=None, window=6, central_measure=
     Biology Society (EMBC), doi: 10.1109/EMBC.2016.7591347.
     """
     data, hrv = _cardiac_metrics(
-        data, metric="hrv", fs=None, peaks=None, window=6, central_measure="std"
+        data,
+        metric="hrv",
+        fs=fs,
+        peaks=peaks,
+        window=window,
+        central_measure=central_measure,
     )
-    data._computed_metrics["heart_rate_variability"] = dict(metric=hrv)
 
     return data, hrv
 
 
 @due.dcite(references.CHEN_2020)
+@return_physio_or_metric()
 @physio.make_operation()
 def heart_beat_interval(data, fs=None, peaks=None, window=6, central_measure="mean"):
     """
@@ -313,13 +330,18 @@ def heart_beat_interval(data, fs=None, peaks=None, window=6, central_measure="me
         vol. 213, pp. 116707, 2020.
     """
     data, hbi = _cardiac_metrics(
-        data, metric="hbi", fs=None, peaks=None, window=6, central_measure="mean"
+        data,
+        metric="hbi",
+        fs=fs,
+        peaks=peaks,
+        window=window,
+        central_measure=central_measure,
     )
-    data._computed_metrics["heart_beat_interval"] = dict(metric=hbi)
 
     return data, hbi
 
 
+@return_physio_or_metric()
 @physio.make_operation()
 def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
     """Calculate cardiac phase from cardiac peaks.
@@ -402,6 +424,4 @@ def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
             ) / (t2 - t1)
         phase_card[:, i_slice] = phase_card_crSlice
 
-    data._computed_metrics["cardiac_phase"] = dict(metric=phase_card)
-
     return data, phase_card
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 1c20aeb..0a7b800 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -9,10 +9,11 @@
 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)
+@return_physio_or_metric()
 @physio.make_operation()
 def respiratory_variance_time(
     data, fs=None, peaks=None, troughs=None, lags=(0, 4, 8, 12)
@@ -106,11 +107,11 @@ def respiratory_variance_time(
         )
         rvt_lags[:, ind] = temp_rvt
 
-    data._computed_metrics["rvt"] = dict(metric=rvt_lags, has_lags=True)
     return data, rvt_lags
 
 
 @due.dcite(references.POWER_2018)
+@return_physio_or_metric()
 @physio.make_operation()
 def respiratory_pattern_variability(data, window):
     """Calculate respiratory pattern variability.
@@ -154,11 +155,11 @@ def respiratory_pattern_variability(data, window):
     # Calculate standard deviation
     rpv_val = np.std(rpv_upper_env)
 
-    data._computed_metrics["rpv"] = dict(metric=rpv_val)
     return data, rpv_val
 
 
 @due.dcite(references.POWER_2020)
+@return_physio_or_metric()
 @physio.make_operation()
 def env(data, fs=None, window=10):
     """Calculate respiratory pattern variability across a sliding window.
@@ -238,11 +239,11 @@ def _respiratory_pattern_variability(data, window):
     )
     env_arr[np.isnan(env_arr)] = 0.0
 
-    data._computed_metrics["env"] = dict(metric=env_arr)
     return data, env_arr
 
 
 @due.dcite(references.CHANG_GLOVER_2009)
+@return_physio_or_metric()
 @physio.make_operation()
 def respiratory_variance(data, fs=None, window=6):
     """Calculate respiratory variance.
@@ -304,11 +305,10 @@ def respiratory_variance(data, fs=None, window=6):
     # Convolve with rrf
     rv_out = convolve_and_rescale(rv_arr, rrf(data.fs), rescale="zscore")
 
-    data._computed_metrics["respiratory_variance"] = dict(metric=rv_out)
-
     return data, rv_out
 
 
+@return_physio_or_metric()
 @physio.make_operation()
 def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None):
     """Calculate respiratory phase from respiratory signal.
@@ -374,6 +374,4 @@ def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None):
 
         phase_resp[:, i_slice] = phase_resp_crSlice
 
-    data._computed_metrics["respiratory_phase"] = dict(metric=phase_resp)
-
     return data, phase_resp
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index 34fb21b..4378866 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -1,8 +1,11 @@
 """Miscellaneous utility functions for metric calculation."""
+import functools
 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
 
@@ -332,3 +335,38 @@ 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):
+        @functools.wraps(func)
+        def wrapper(*args, **kwargs):
+            physio, metric = func(*args, **kwargs)
+            if isinstance(args[0], Physio):
+                physio._computed_metrics[func.__name__] = dict(
+                    metric=metric, args=kwargs
+                )
+                if return_physio:
+                    return physio, metric
+                else:
+                    return metric
+            else:
+                return metric
+
+        return wrapper
+
+    return determine_return_type
diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index 25faa9a..f9b69ad 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -28,7 +28,7 @@ def test_cardiac_phase_smoke():
     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(
+    card_phase = cardiac.cardiac_phase(
         data,
         peaks=peaks,
         fs=sample_rate,
diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py
index 5e5bfd1..a4257ec 100644
--- a/phys2denoise/tests/test_metrics_chest_belt.py
+++ b/phys2denoise/tests/test_metrics_chest_belt.py
@@ -30,7 +30,7 @@ def test_respiratory_phase_smoke():
     slice_timings = np.linspace(0, t_r, 22)[1:-1]
     n_samples = int(np.rint((n_scans * t_r) * sample_rate))
     resp = np.random.normal(size=n_samples)
-    _, resp_phase = chest_belt.respiratory_phase(
+    resp_phase = chest_belt.respiratory_phase(
         resp,
         fs=sample_rate,
         slice_timings=slice_timings,
@@ -46,7 +46,7 @@ def test_respiratory_pattern_variability_smoke():
     n_samples = 2000
     resp = np.random.normal(size=n_samples)
     window = 50
-    _, rpv_val = chest_belt.respiratory_pattern_variability(resp, window)
+    rpv_val = chest_belt.respiratory_pattern_variability(resp, window)
     assert isinstance(rpv_val, float)
 
 
@@ -56,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, fs=samplerate, window=window)
+    env_arr = chest_belt.env(resp, fs=samplerate, window=window)
     assert env_arr.ndim == 1
     assert env_arr.shape == (n_samples,)
 
@@ -67,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, fs=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_rvt.py b/phys2denoise/tests/test_rvt.py
index 67a7bf9..20d1690 100644
--- a/phys2denoise/tests/test_rvt.py
+++ b/phys2denoise/tests/test_rvt.py
@@ -17,8 +17,8 @@ def test_respiratory_variance_time(fake_phys):
     phys = peakdet.operations.peakfind_physio(phys)
 
     # TODO: Change to a simpler call once physutils are
-    # integrated to peakdet
-    phys, r = respiratory_variance_time(
+    # integrated to peakdet/prep4phys
+    r = respiratory_variance_time(
         phys.data, fs=phys.fs, peaks=phys.peaks, troughs=phys.troughs
     )
     assert r is not None

From 9883ce296f7824faa757d9feeeb3c01c912374ea Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sat, 10 Aug 2024 20:40:48 +0300
Subject: [PATCH 31/38] [deps] Cap numpy to <2.0

---
 setup.cfg | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/setup.cfg b/setup.cfg
index 9f26a68..f1a08fd 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -21,8 +21,8 @@ provides =
 [options]
 python_requires = >=3.9
 install_requires =
-    numpy >=1.9.3
-    matplotlib>=3.9
+    numpy >=1.9.3, <2
+    matplotlib
     pandas
     scipy
     duecredit

From 181304fe10d6b785b4577cf8ac1883b276bd37c8 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sat, 10 Aug 2024 20:51:09 +0300
Subject: [PATCH 32/38] [docs]: Update computing metrics

---
 docs/user_guide/metrics.rst | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

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.

From 43107db3d956f316255e44a934d28edc7eba5a0e Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sun, 11 Aug 2024 21:35:19 +0300
Subject: [PATCH 33/38] Metric or Physio returning optimization

---
 phys2denoise/metrics/cardiac.py            | 17 ++++++++++++-----
 phys2denoise/metrics/chest_belt.py         | 10 +++++-----
 phys2denoise/metrics/multimodal.py         |  1 +
 phys2denoise/metrics/utils.py              |  5 +++--
 phys2denoise/tests/test_metrics_cardiac.py | 18 +++++++++++++++++-
 5 files changed, 38 insertions(+), 13 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 72c7e8d..63f1475 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -17,7 +17,7 @@ def _cardiac_metrics(
     peaks=None,
     window=6,
     central_measure="mean",
-    return_physio=False,
+    **kwargs,
 ):
     """
     Compute cardiac metrics.
@@ -154,7 +154,7 @@ def _cardiac_metrics(
 @due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
 @return_physio_or_metric()
 @physio.make_operation()
-def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean"):
+def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean", **kwargs):
     """
     Compute average heart rate (HR) in a sliding window.
 
@@ -213,6 +213,7 @@ def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean"):
         peaks=peaks,
         window=window,
         central_measure=central_measure,
+        **kwargs,
     )
 
     return data, hr
@@ -221,7 +222,9 @@ def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean"):
 @due.dcite(references.PINHERO_ET_AL_2016)
 @return_physio_or_metric()
 @physio.make_operation()
-def heart_rate_variability(data, fs=None, peaks=None, window=6, central_measure="mean"):
+def heart_rate_variability(
+    data, fs=None, peaks=None, window=6, central_measure="mean", **kwargs
+):
     """
     Compute average heart rate variability (HRV) in a sliding window.
 
@@ -278,6 +281,7 @@ def heart_rate_variability(data, fs=None, peaks=None, window=6, central_measure=
         peaks=peaks,
         window=window,
         central_measure=central_measure,
+        **kwargs,
     )
 
     return data, hrv
@@ -286,7 +290,9 @@ def heart_rate_variability(data, fs=None, peaks=None, window=6, central_measure=
 @due.dcite(references.CHEN_2020)
 @return_physio_or_metric()
 @physio.make_operation()
-def heart_beat_interval(data, fs=None, peaks=None, window=6, central_measure="mean"):
+def heart_beat_interval(
+    data, fs=None, peaks=None, window=6, central_measure="mean", **kwargs
+):
     """
     Compute average heart beat interval (HBI) in a sliding window.
 
@@ -336,6 +342,7 @@ def heart_beat_interval(data, fs=None, peaks=None, window=6, central_measure="me
         peaks=peaks,
         window=window,
         central_measure=central_measure,
+        **kwargs,
     )
 
     return data, hbi
@@ -343,7 +350,7 @@ def heart_beat_interval(data, fs=None, peaks=None, window=6, central_measure="me
 
 @return_physio_or_metric()
 @physio.make_operation()
-def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None):
+def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None, **kwargs):
     """Calculate cardiac phase from cardiac peaks.
 
     Assumes that timing of cardiac events are given in same units
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 0a7b800..deae46f 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -16,7 +16,7 @@
 @return_physio_or_metric()
 @physio.make_operation()
 def respiratory_variance_time(
-    data, fs=None, peaks=None, troughs=None, lags=(0, 4, 8, 12)
+    data, fs=None, peaks=None, troughs=None, lags=(0, 4, 8, 12), **kwargs
 ):
     """
     Implement the Respiratory Variance over Time (Birn et al. 2006).
@@ -113,7 +113,7 @@ def respiratory_variance_time(
 @due.dcite(references.POWER_2018)
 @return_physio_or_metric()
 @physio.make_operation()
-def respiratory_pattern_variability(data, window):
+def respiratory_pattern_variability(data, window, **kwargs):
     """Calculate respiratory pattern variability.
 
     Parameters
@@ -161,7 +161,7 @@ def respiratory_pattern_variability(data, window):
 @due.dcite(references.POWER_2020)
 @return_physio_or_metric()
 @physio.make_operation()
-def env(data, fs=None, window=10):
+def env(data, fs=None, window=10, **kwargs):
     """Calculate respiratory pattern variability across a sliding window.
 
     Parameters
@@ -245,7 +245,7 @@ def _respiratory_pattern_variability(data, window):
 @due.dcite(references.CHANG_GLOVER_2009)
 @return_physio_or_metric()
 @physio.make_operation()
-def respiratory_variance(data, fs=None, window=6):
+def respiratory_variance(data, fs=None, window=6, **kwargs):
     """Calculate respiratory variance.
 
     Parameters
@@ -310,7 +310,7 @@ def respiratory_variance(data, fs=None, window=6):
 
 @return_physio_or_metric()
 @physio.make_operation()
-def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None):
+def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None, **kwargs):
     """Calculate respiratory phase from respiratory signal.
 
     Parameters
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index fc48213..85f168f 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -20,6 +20,7 @@ def retroicor(
     physio_type=None,
     fs=None,
     cardiac_peaks=None,
+    **kwargs,
 ):
     """Compute RETROICOR regressors.
 
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index 4378866..edd66ec 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -360,8 +360,9 @@ def wrapper(*args, **kwargs):
                 physio._computed_metrics[func.__name__] = dict(
                     metric=metric, args=kwargs
                 )
-                if return_physio:
-                    return physio, metric
+                return_physio_value = kwargs.get("return_physio", return_physio)
+                if return_physio_value:
+                    return physio
                 else:
                     return metric
             else:
diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index f9b69ad..839b040 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -50,12 +50,28 @@ def test_cardiac_phase_smoke_physio_obj():
     data = np.zeros(peaks.shape)
     phys = physio.Physio(data, sample_rate, physio_type="cardiac")
     phys._metadata["peaks"] = peaks
-    phys, card_phase = cardiac.cardiac_phase(
+
+    # 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"]["metric"].ndim == 2
+    assert phys.computed_metrics["cardiac_phase"]["metric"].shape == (
+        n_scans,
+        slice_timings.size,
+    )
+
+    # 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)

From b0f411017b646e421be11fefd010e796b35d73c9 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sun, 11 Aug 2024 22:06:47 +0300
Subject: [PATCH 34/38] Add Metric class for easier properties accessing

---
 phys2denoise/metrics/utils.py              | 62 +++++++++++++++++++++-
 phys2denoise/tests/test_metrics_cardiac.py |  8 ++-
 2 files changed, 66 insertions(+), 4 deletions(-)

diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index edd66ec..eab3e7a 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -353,13 +353,29 @@ def return_physio_or_metric(*, return_physio=True):
     """
 
     def determine_return_type(func):
+        metrics_with_lags = ["respiratory_variance_time"]
+        convolved_metrics = [
+            "respiratory_variance",
+            "heart_rate",
+            "heart_rate_variability",
+            "heart_beat_interval",
+        ]
+
         @functools.wraps(func)
         def wrapper(*args, **kwargs):
             physio, metric = func(*args, **kwargs)
             if isinstance(args[0], Physio):
-                physio._computed_metrics[func.__name__] = dict(
-                    metric=metric, args=kwargs
+                has_lags = True if func.__name__ in metrics_with_lags else 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
@@ -371,3 +387,45 @@ def wrapper(*args, **kwargs):
         return wrapper
 
     return determine_return_type
+
+
+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/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py
index 839b040..2d776cf 100644
--- a/phys2denoise/tests/test_metrics_cardiac.py
+++ b/phys2denoise/tests/test_metrics_cardiac.py
@@ -58,12 +58,16 @@ def test_cardiac_phase_smoke_physio_obj():
         n_scans=n_scans,
         t_r=t_r,
     )
+
     assert phys.history[0][0] == "phys2denoise.metrics.cardiac.cardiac_phase"
-    assert phys.computed_metrics["cardiac_phase"]["metric"].ndim == 2
-    assert phys.computed_metrics["cardiac_phase"]["metric"].shape == (
+    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(

From 10509bd43b843480d19121e60e16340300d498e6 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Wed, 21 Aug 2024 17:13:49 +0200
Subject: [PATCH 35/38] More unit tests, physio/metric wrapper optimization,
 doc fixes

---
 .all-contributorsrc                |  2 ++
 phys2denoise/metrics/cardiac.py    | 10 ++++----
 phys2denoise/metrics/chest_belt.py | 10 ++++----
 phys2denoise/metrics/multimodal.py |  8 ++++--
 phys2denoise/metrics/utils.py      | 25 +++++++++++++++++--
 phys2denoise/tests/test_rvt.py     | 39 ++++++++++++++++++++++++++++++
 6 files changed, 80 insertions(+), 14 deletions(-)

diff --git a/.all-contributorsrc b/.all-contributorsrc
index e8b53b1..1e8e6aa 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -119,6 +119,8 @@
       "contributions": [
         "code",
         "ideas",
+        "infra",
+        "docs",
         "bug",
         "review"
       ]
diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 63f1475..711aacb 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -31,7 +31,7 @@ def _cardiac_metrics(
 
     Parameters
     ----------
-    data : Physio_like
+    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.
@@ -164,7 +164,7 @@ def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean", **kw
 
     Parameters
     ----------
-    data : Physio_like
+    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
@@ -234,7 +234,7 @@ def heart_rate_variability(
 
     Parameters
     ----------
-    data : Physio_like
+    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
@@ -298,7 +298,7 @@ def heart_beat_interval(
 
     Parameters
     ----------
-    data : Physio_like
+    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
@@ -358,7 +358,7 @@ def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None, **kwar
 
     Parameters
     ----------
-    data : Physio_like
+    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.
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index deae46f..4f7fc2c 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -25,7 +25,7 @@ def respiratory_variance_time(
 
     Parameters
     ----------
-    data : Physio_like
+    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
@@ -118,7 +118,7 @@ def respiratory_pattern_variability(data, window, **kwargs):
 
     Parameters
     ----------
-    data : Physio_like
+    data : physutils.Physio, np.ndarray, or array-like object
         Object containing the timeseries of the recorded respiratory signal
     window : int
         Window length in samples.
@@ -166,7 +166,7 @@ def env(data, fs=None, window=10, **kwargs):
 
     Parameters
     ----------
-    data : Physio_like
+    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
@@ -250,7 +250,7 @@ def respiratory_variance(data, fs=None, window=6, **kwargs):
 
     Parameters
     ----------
-    data : Physio_like
+    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
@@ -315,7 +315,7 @@ def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None, **kwargs):
 
     Parameters
     ----------
-    data : Physio_like
+    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.
diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py
index 85f168f..fa28bae 100644
--- a/phys2denoise/metrics/multimodal.py
+++ b/phys2denoise/metrics/multimodal.py
@@ -7,9 +7,11 @@
 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(
     data,
@@ -26,7 +28,7 @@ def retroicor(
 
     Parameters
     ----------
-    data : Physio_like
+    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.
@@ -138,4 +140,6 @@ def retroicor(
     data._computed_metrics["retroicor_regressors"] = dict(
         metrics=retroicor_regressors, phase=phase
     )
-    return data, retroicor_regressors, phase
+    retroicor_regressors = dict(metrics=retroicor_regressors, phase=phase)
+
+    return data, retroicor_regressors
diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index eab3e7a..cf14ef6 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -1,5 +1,6 @@
 """Miscellaneous utility functions for metric calculation."""
 import functools
+import inspect
 import logging
 
 import numpy as np
@@ -353,7 +354,6 @@ def return_physio_or_metric(*, return_physio=True):
     """
 
     def determine_return_type(func):
-        metrics_with_lags = ["respiratory_variance_time"]
         convolved_metrics = [
             "respiratory_variance",
             "heart_rate",
@@ -364,8 +364,15 @@ def determine_return_type(func):
         @functools.wraps(func)
         def wrapper(*args, **kwargs):
             physio, metric = func(*args, **kwargs)
+            default_args = get_default_args(func)
             if isinstance(args[0], Physio):
-                has_lags = True if func.__name__ in metrics_with_lags else False
+                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(
@@ -389,6 +396,20 @@ def wrapper(*args, **kwargs):
     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
diff --git a/phys2denoise/tests/test_rvt.py b/phys2denoise/tests/test_rvt.py
index 20d1690..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
 
@@ -23,3 +24,41 @@ def test_respiratory_variance_time(fake_phys):
     )
     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

From 612c616f3ba26d3048f6d80fbe25a17b27322d25 Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Wed, 21 Aug 2024 17:48:09 +0200
Subject: [PATCH 36/38] Minor fix

---
 phys2denoise/metrics/cardiac.py    | 16 ++++++++--------
 phys2denoise/metrics/chest_belt.py |  2 +-
 2 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py
index 711aacb..5b7588d 100644
--- a/phys2denoise/metrics/cardiac.py
+++ b/phys2denoise/metrics/cardiac.py
@@ -13,8 +13,8 @@
 def _cardiac_metrics(
     data,
     metric,
-    fs=None,
     peaks=None,
+    fs=None,
     window=6,
     central_measure="mean",
     **kwargs,
@@ -154,7 +154,7 @@ def _cardiac_metrics(
 @due.dcite(references.CHANG_CUNNINGHAM_GLOVER_2009)
 @return_physio_or_metric()
 @physio.make_operation()
-def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean", **kwargs):
+def heart_rate(data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs):
     """
     Compute average heart rate (HR) in a sliding window.
 
@@ -209,8 +209,8 @@ def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean", **kw
     data, hr = _cardiac_metrics(
         data,
         metric="hr",
-        fs=fs,
         peaks=peaks,
+        fs=fs,
         window=window,
         central_measure=central_measure,
         **kwargs,
@@ -223,7 +223,7 @@ def heart_rate(data, fs=None, peaks=None, window=6, central_measure="mean", **kw
 @return_physio_or_metric()
 @physio.make_operation()
 def heart_rate_variability(
-    data, fs=None, peaks=None, window=6, central_measure="mean", **kwargs
+    data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs
 ):
     """
     Compute average heart rate variability (HRV) in a sliding window.
@@ -277,8 +277,8 @@ def heart_rate_variability(
     data, hrv = _cardiac_metrics(
         data,
         metric="hrv",
-        fs=fs,
         peaks=peaks,
+        fs=fs,
         window=window,
         central_measure=central_measure,
         **kwargs,
@@ -291,7 +291,7 @@ def heart_rate_variability(
 @return_physio_or_metric()
 @physio.make_operation()
 def heart_beat_interval(
-    data, fs=None, peaks=None, window=6, central_measure="mean", **kwargs
+    data, peaks=None, fs=None, window=6, central_measure="mean", **kwargs
 ):
     """
     Compute average heart beat interval (HBI) in a sliding window.
@@ -338,8 +338,8 @@ def heart_beat_interval(
     data, hbi = _cardiac_metrics(
         data,
         metric="hbi",
-        fs=fs,
         peaks=peaks,
+        fs=fs,
         window=window,
         central_measure=central_measure,
         **kwargs,
@@ -350,7 +350,7 @@ def heart_beat_interval(
 
 @return_physio_or_metric()
 @physio.make_operation()
-def cardiac_phase(data, slice_timings, n_scans, t_r, fs=None, peaks=None, **kwargs):
+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
diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py
index 4f7fc2c..4e080fc 100644
--- a/phys2denoise/metrics/chest_belt.py
+++ b/phys2denoise/metrics/chest_belt.py
@@ -16,7 +16,7 @@
 @return_physio_or_metric()
 @physio.make_operation()
 def respiratory_variance_time(
-    data, fs=None, peaks=None, troughs=None, lags=(0, 4, 8, 12), **kwargs
+    data, peaks=None, troughs=None, fs=None, lags=(0, 4, 8, 12), **kwargs
 ):
     """
     Implement the Respiratory Variance over Time (Birn et al. 2006).

From e528d216883dde68447837ab25cdb18713e06f7f Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sun, 25 Aug 2024 15:06:04 +0200
Subject: [PATCH 37/38] [test] Fix failing warning assertion

---
 phys2denoise/metrics/utils.py            | 8 ++++++--
 phys2denoise/tests/test_metrics_utils.py | 6 +++---
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py
index cf14ef6..e21a085 100644
--- a/phys2denoise/metrics/utils.py
+++ b/phys2denoise/metrics/utils.py
@@ -37,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):
@@ -62,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."
         )
diff --git a/phys2denoise/tests/test_metrics_utils.py b/phys2denoise/tests/test_metrics_utils.py
index d7e1b46..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(IndexError):
-        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():

From 6266b6f78950d3c924a97249272fb4acb6c1264d Mon Sep 17 00:00:00 2001
From: maestroque <kikasgeorg@ece.auth.gr>
Date: Sun, 25 Aug 2024 15:36:26 +0200
Subject: [PATCH 38/38] Add updated physutils version

---
 setup.cfg | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/setup.cfg b/setup.cfg
index f1a08fd..81268a1 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -27,7 +27,7 @@ install_requires =
     scipy
     duecredit
     loguru
-    physutils
+    physutils >=0.2.1
 tests_require =
     pytest >=5.3
 test_suite = pytest