Merge pull request #32 from tsalo/tests

Add smoke tests for metrics and fix miscellaneous bugs

Author: Stefano Moia

phys2denoise/metrics/cardiac.py

@@ -61,7 +61,7 @@ def _crf(t):
     return crf_arr
 
 
-def cardiac_phase(peaks, slice_timings, n_scans, t_r):
+def cardiac_phase(peaks, sample_rate, slice_timings, n_scans, t_r):
     """Calculate cardiac phase from cardiac peaks.
 
     Assumes that timing of cardiac events are given in same units
@@ -71,6 +71,8 @@ def cardiac_phase(peaks, slice_timings, n_scans, t_r):
     ----------
     peaks : 1D array_like
         Cardiac peak times, in seconds.
+    sample_rate : float
+        Sample rate of physio, in Hertz.
     slice_timings : 1D array_like
         Slice times, in seconds.
     n_scans : int
@@ -81,31 +83,37 @@ def cardiac_phase(peaks, slice_timings, n_scans, t_r):
     Returns
     -------
     phase_card : array_like
-        Cardiac phase signal.
+        Cardiac phase signal, of shape (n_scans,)
     """
-    n_slices = np.shape(slice_timings)
-    phase_card = np.zeros(n_scans)
+    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
     for i_slice in range(n_slices):
-        # find previous cardiac peaks
-        previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice]))
-        if np.size(previous_card_peaks) == 0:
-            t1 = 0
-        else:
-            last_peak = previous_card_peaks[0][-1]
-            t1 = peaks[last_peak]
-
-        # find posterior cardiac peaks
-        next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice]))
-        if np.size(next_card_peaks) == 0:
-            t2 = n_scans * t_r
-        else:
-            next_peak = next_card_peaks[0][0]
-            t2 = peaks[next_peak]
-
-        # compute cardiac phase
-        phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / (
-            t2 - t1
-        )
+        # generate slice acquisition timings across all scans
+        times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice]
+        phase_card_crSlice = np.zeros(n_scans)
+        for j_scan in range(n_scans):
+            previous_card_peaks = np.asarray(
+                np.nonzero(card_peaks_sec < times_crSlice[j_scan])
+            )
+            if np.size(previous_card_peaks) == 0:
+                t1 = 0
+            else:
+                last_peak = previous_card_peaks[0][-1]
+                t1 = card_peaks_sec[last_peak]
+            next_card_peaks = np.asarray(
+                np.nonzero(card_peaks_sec > times_crSlice[j_scan])
+            )
+            if np.size(next_card_peaks) == 0:
+                t2 = n_scans * t_r
+            else:
+                next_peak = next_card_peaks[0][0]
+                t2 = card_peaks_sec[next_peak]
+            phase_card_crSlice[j_scan] = (
+                2 * np.math.pi * (times_crSlice[j_scan] - t1)
+            ) / (t2 - t1)
+        phase_card[:, i_slice] = phase_card_crSlice
 
     return phase_card

phys2denoise/metrics/chest_belt.py

@@ -3,7 +3,7 @@
 import numpy as np
 import pandas as pd
 from scipy.ndimage.filters import convolve1d
-from scipy.signal import detrend, resample
+from scipy.signal import detrend
 from scipy.stats import zscore
 
 mpl.use("TkAgg")
@@ -15,17 +15,18 @@
 
 
 @due.dcite(references.POWER_2018)
-def rpv(belt_ts, window):
+def rpv(resp, window=6):
     """Calculate respiratory pattern variability.
 
     Parameters
     ----------
-    belt_ts
-    window
+    resp : 1D array_like
+    window : int
 
     Returns
     -------
-    rpv_arr
+    rpv_val : float
+        Respiratory pattern variability value.
 
     Notes
     -----
@@ -43,35 +44,29 @@ def rpv(belt_ts, window):
        115, pp. 2105-2114, 2018.
     """
     # First, z-score respiratory traces
-    resp_z = zscore(belt_ts)
+    resp_z = zscore(resp)
 
     # Collect upper envelope
     rpv_upper_env = utils.rms_envelope_1d(resp_z, window)
 
     # Calculate standard deviation
-    rpv_arr = np.std(rpv_upper_env)
-    return rpv_arr
+    rpv_val = np.std(rpv_upper_env)
+    return rpv_val
 
 
 @due.dcite(references.POWER_2020)
-def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
+def env(resp, samplerate, window=10):
     """Calculate respiratory pattern variability across a sliding window.
 
     Parameters
     ----------
-    belt_ts : (X,) :obj:`numpy.ndarray`
+    resp : (X,) :obj:`numpy.ndarray`
         A 1D array with the respiratory belt time series.
     samplerate : :obj:`float`
-        Sampling rate for belt_ts, in Hertz.
-    out_samplerate : :obj:`float`
-        Sampling rate for the output time series in seconds.
-        Corresponds to TR in fMRI data.
+        Sampling rate for resp, in Hertz.
     window : :obj:`int`, optional
-        Size of the sliding window, in the same units as out_samplerate.
-        Default is 6.
-    lags : (Y,) :obj:`tuple` of :obj:`int`, optional
-        List of lags to apply to the rv estimate. In the same units as
-        out_samplerate.
+        Size of the sliding window, in seconds.
+        Default is 10.
 
     Returns
     -------
@@ -92,42 +87,37 @@ def env(belt_ts, samplerate, out_samplerate, window=10, lags=(0,)):
        young adults scanned at rest, including systematic changes and 'missed'
        deep breaths," Neuroimage, vol. 204, 2020.
     """
-    window = window * samplerate / out_samplerate
+    # Convert window to Hertz
+    window = int(window * samplerate)
+
     # Calculate RPV across a rolling window
     env_arr = (
-        pd.Series(belt_ts).rolling(window=window, center=True).apply(rpv, window=window)
+        pd.Series(resp).rolling(window=window, center=True).apply(rpv, args=(window,))
     )
     env_arr[np.isnan(env_arr)] = 0.0
     return env_arr
 
 
 @due.dcite(references.CHANG_GLOVER_2009)
-def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
+def rv(resp, samplerate, window=6, lags=(0,)):
     """Calculate respiratory variance.
 
     Parameters
     ----------
-    belt_ts : (X,) :obj:`numpy.ndarray`
+    resp : (X,) :obj:`numpy.ndarray`
         A 1D array with the respiratory belt time series.
     samplerate : :obj:`float`
-        Sampling rate for belt_ts, in Hertz.
-    out_samplerate : :obj:`float`
-        Sampling rate for the output time series in seconds.
-        Corresponds to TR in fMRI data.
+        Sampling rate for resp, in Hertz.
     window : :obj:`int`, optional
-        Size of the sliding window, in the same units as out_samplerate.
+        Size of the sliding window, in seconds.
         Default is 6.
-    lags : (Y,) :obj:`tuple` of :obj:`int`, optional
-        List of lags to apply to the rv estimate. In the same units as
-        out_samplerate.
 
     Returns
     -------
-    rv_out : (Z, 2Y) :obj:`numpy.ndarray`
-        Respiratory variance values, with lags applied, downsampled to
-        out_samplerate, convolved with an RRF, and detrended/normalized.
-        The first Y columns are not convolved with the RRF, while the second Y
-        columns are.
+    rv_out : (X, 2) :obj:`numpy.ndarray`
+        Respiratory variance values.
+        The first column is raw RV values, after detrending/normalization.
+        The second column is RV values convolved with the RRF, after detrending/normalization.
 
     Notes
     -----
@@ -145,25 +135,19 @@ def rv(belt_ts, samplerate, out_samplerate, window=6, lags=(0,)):
        end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
        issue 4, vol. 47, pp. 1381-1393, 2009.
     """
+    # Convert window to Hertz
+    window = int(window * samplerate)
+
     # Raw respiratory variance
-    rv_arr = pd.Series(belt_ts).rolling(window=window, center=True).std()
+    rv_arr = pd.Series(resp).rolling(window=window, center=True).std()
     rv_arr[np.isnan(rv_arr)] = 0.0
 
-    # Apply lags
-    n_out_samples = int((belt_ts.shape[0] / samplerate) / out_samplerate)
-    # convert lags from out_samplerate to samplerate
-    delays = [abs(int(lag * samplerate)) for lag in lags]
-    rv_with_lags = utils.apply_lags(rv_arr, lags=delays)
-
-    # Downsample to out_samplerate
-    rv_with_lags = resample(rv_with_lags, num=n_out_samples, axis=0)
-
     # Convolve with rrf
-    rrf_arr = rrf(out_samplerate, oversampling=1)
-    rv_convolved = convolve1d(rv_with_lags, rrf_arr, axis=0)
+    rrf_arr = rrf(samplerate, oversampling=1)
+    rv_convolved = convolve1d(rv_arr, rrf_arr, axis=0)
 
     # Concatenate the raw and convolved versions
-    rv_combined = np.hstack((rv_with_lags, rv_convolved))
+    rv_combined = np.stack((rv_arr, rv_convolved), axis=-1)
 
     # Detrend and normalize
     rv_combined = rv_combined - np.mean(rv_combined, axis=0)
@@ -240,14 +224,15 @@ def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r):
     Returns
     -------
     phase_resp : array_like
-        Respiratory phase signal.
+        Respiratory phase signal, of shape (n_scans, n_slices).
     """
-    n_slices = np.shape(slice_timings)
+    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 = plt.hist(resp, bins=100)
+    resp_hist, resp_hist_bins, _ = plt.hist(resp, bins=100)
 
     # first compute derivative of respiration signal
     resp_diff = np.diff(resp, n=1)

phys2denoise/tests/test_metrics_cardiac.py

@@ -0,0 +1,41 @@
+"""Tests for phys2denoise.metrics.cardiac."""
+import numpy as np
+
+from phys2denoise.metrics import cardiac
+
+
+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,
+    )
+    pred_len = np.rint(time_length / (tr / oversampling))
+    assert crf_arr.ndim == 1
+    assert crf_arr.shape == pred_len
+
+
+def test_cardiac_phase_smoke():
+    """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])
+    card_phase = cardiac.cardiac_phase(
+        peaks,
+        sample_rate=sample_rate,
+        slice_timings=slice_timings,
+        n_scans=n_scans,
+        t_r=t_r,
+    )
+    assert card_phase.ndim == 2
+    assert card_phase.shape == (n_scans, slice_timings.size)

phys2denoise/tests/test_metrics_chest_belt.py

@@ -0,0 +1,73 @@
+"""Tests for phys2denoise.metrics.chest_belt."""
+import numpy as np
+
+from phys2denoise.metrics import chest_belt
+
+
+def test_rrf_smoke():
+    """Basic smoke test for RRF calculation."""
+    samplerate = 0.01  # in seconds
+    oversampling = 20
+    time_length = 20
+    onset = 0
+    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)))
+    assert rrf_arr.ndim == 1
+    assert rrf_arr.size == pred_len
+
+
+def test_respiratory_phase_smoke():
+    """Basic smoke test for respiratory phase calculation."""
+    t_r = 1.0
+    n_scans = 200
+    sample_rate = 1 / 0.01
+    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,
+        sample_rate=sample_rate,
+        slice_timings=slice_timings,
+        n_scans=n_scans,
+        t_r=t_r,
+    )
+    assert resp_phase.ndim == 2
+    assert resp_phase.shape == (n_scans, slice_timings.size)
+
+
+def test_rpv_smoke():
+    """Basic smoke test for respiratory pattern variability calculation."""
+    n_samples = 2000
+    resp = np.random.normal(size=n_samples)
+    window = 50
+    rpv_val = chest_belt.rpv(resp, window)
+    assert isinstance(rpv_val, float)
+
+
+def test_env_smoke():
+    """Basic smoke test for ENV calculation."""
+    n_samples = 2000
+    resp = np.random.normal(size=n_samples)
+    samplerate = 1 / 0.01
+    window = 6
+    env_arr = chest_belt.env(resp, samplerate=samplerate, window=window)
+    assert env_arr.ndim == 1
+    assert env_arr.shape == (n_samples,)
+
+
+def test_rv_smoke():
+    """Basic smoke test for respiratory variance calculation."""
+    n_samples = 2000
+    resp = np.random.normal(size=n_samples)
+    samplerate = 1 / 0.01
+    window = 6
+    rv_arr = chest_belt.rv(resp, samplerate=samplerate, window=window)
+    assert rv_arr.ndim == 2
+    assert rv_arr.shape == (n_samples, 2)