Add session 0008, artifact censoring, smoothing, and EGG survey script

- Bundle session 0008 EGG data and add BOLD to GitHub release
- Add boolean mask parameter to PLV functions for artifact censoring
- Add artifact_mask_to_volumes helper to map EGG artifacts to volumes
- Update fMRI tutorial: session 0008, 6mm smoothing, artifact censoring, 200 surrogates
- Add survey_egg_sessions.py script for cross-session EGG quality comparison
- Add tests for mask parameter and artifact_mask_to_volumes (208 total)

Author: micahgallen

docs/tutorials/fmri_coupling_real.ipynb

@@ -15,7 +15,7 @@
 import numpy as np
 
 
-def phase_locking_value(phase_a, phase_b):
+def phase_locking_value(phase_a, phase_b, mask=None):
     """Compute the phase-locking value between two phase time series.
 
     PLV measures the consistency of the phase difference between two
@@ -30,6 +30,9 @@ def phase_locking_value(phase_a, phase_b):
     phase_b : array_like, shape (n_timepoints,)
         Reference phase time series in radians. Broadcast against
         columns of ``phase_a`` when ``phase_a`` is 2D.
+    mask : array_like of bool, shape (n_timepoints,), optional
+        Boolean mask where ``True`` = include, ``False`` = exclude.
+        Only the included timepoints contribute to the PLV.
 
     Returns
     -------
@@ -51,10 +54,10 @@ def phase_locking_value(phase_a, phase_b):
     >>> round(plv, 2)
     1.0
     """
-    return np.abs(phase_locking_value_complex(phase_a, phase_b))
+    return np.abs(phase_locking_value_complex(phase_a, phase_b, mask=mask))
 
 
-def phase_locking_value_complex(phase_a, phase_b):
+def phase_locking_value_complex(phase_a, phase_b, mask=None):
     """Compute the complex phase-locking value.
 
     Returns the complex mean of the phase difference, from which both
@@ -66,6 +69,9 @@ def phase_locking_value_complex(phase_a, phase_b):
         Phase time series in radians.
     phase_b : array_like, shape (n_timepoints,)
         Reference phase time series in radians.
+    mask : array_like of bool, shape (n_timepoints,), optional
+        Boolean mask where ``True`` = include, ``False`` = exclude.
+        Only the included timepoints contribute to the complex mean.
 
     Returns
     -------
@@ -87,23 +93,25 @@ def phase_locking_value_complex(phase_a, phase_b):
     phase_a = np.asarray(phase_a, dtype=float)
     phase_b = np.asarray(phase_b, dtype=float)
 
-    if phase_a.ndim == 1 and phase_b.ndim == 1:
+    if (phase_a.ndim == 1 and phase_b.ndim == 1) or phase_a.ndim == 2:
         if phase_a.shape[0] != phase_b.shape[0]:
             raise ValueError(
                 f"phase_a and phase_b must have the same number of timepoints, "
                 f"got {phase_a.shape[0]} and {phase_b.shape[0]}"
             )
-    elif phase_a.ndim == 2:
-        if phase_a.shape[0] != phase_b.shape[0]:
-            raise ValueError(
-                f"phase_a and phase_b must have the same number of timepoints, "
-                f"got {phase_a.shape[0]} and {phase_b.shape[0]}"
-            )
-        # Broadcast phase_b to match: (n_timepoints,) -> (n_timepoints, 1)
-        phase_b = phase_b[:, np.newaxis]
     else:
         raise ValueError(f"phase_a must be 1D or 2D, got {phase_a.ndim}D")
 
+    # Apply mask before computing phase difference
+    if mask is not None:
+        mask = np.asarray(mask, dtype=bool)
+        phase_a = phase_a[mask]
+        phase_b = phase_b[mask]
+
+    # Broadcast phase_b to match 2D phase_a
+    if phase_a.ndim == 2:
+        phase_b = phase_b[:, np.newaxis]
+
     phase_diff = phase_a - phase_b
     cplv = np.mean(np.exp(1j * phase_diff), axis=0)
 

gastropy/coupling/plv.py

@@ -20,7 +20,7 @@
 from .plv import phase_locking_value
 
 
-def surrogate_plv(phase_a, phase_b, buffer_samples=None, n_surrogates=None, stat="median", seed=None):
+def surrogate_plv(phase_a, phase_b, buffer_samples=None, n_surrogates=None, stat="median", seed=None, mask=None):
     """Compute surrogate PLV via circular time-shifting.
 
     Creates a null distribution of PLV by circularly shifting the
@@ -50,6 +50,9 @@ def surrogate_plv(phase_a, phase_b, buffer_samples=None, n_surrogates=None, stat
         of surrogate PLV values.
     seed : int or np.random.Generator, optional
         Random seed for reproducibility when ``n_surrogates`` is set.
+    mask : array_like of bool, shape (n_timepoints,), optional
+        Boolean mask where ``True`` = include. Passed through to
+        ``phase_locking_value`` for each surrogate shift.
 
     Returns
     -------
@@ -101,7 +104,7 @@ def surrogate_plv(phase_a, phase_b, buffer_samples=None, n_surrogates=None, stat
     surr_plvs = []
     for shift in shifts:
         shifted_b = np.roll(phase_b, int(shift))
-        plv = phase_locking_value(phase_a, shifted_b)
+        plv = phase_locking_value(phase_a, shifted_b, mask=mask)
         surr_plvs.append(plv)
 
     surr_plvs = np.array(surr_plvs)

gastropy/coupling/surrogate.py

@@ -24,7 +24,7 @@
 
 _DATA_DIR = Path(__file__).parent
 
-_FMRI_SESSIONS = ("0001", "0003", "0004")
+_FMRI_SESSIONS = ("0001", "0003", "0004", "0008")
 
 
 def _load_npz(filename):
@@ -44,8 +44,8 @@ def load_fmri_egg(session="0001"):
     Parameters
     ----------
     session : str
-        Session identifier. Available: ``"0001"``, ``"0003"``, ``"0004"``
-        (three baseline sessions from the semi_precision study).
+        Session identifier. Available: ``"0001"``, ``"0003"``, ``"0004"``,
+        ``"0008"`` (from the semi_precision study).
 
     Returns
     -------
@@ -176,6 +176,20 @@ def list_datasets():
             "sha256:7736cee95003ceb8f55a763be94f8130ab356eddbcaef18a8451954a95de53d7",
         ),
     },
+    "0008": {
+        "bold": (
+            "sub-01_ses-20241101_task-rest_run-01_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz",
+            "sha256:f17edfd40e5770bad7e195168fd0009b1d8fb00e80306042b4c0aad2c894cdc8",
+        ),
+        "mask": (
+            "sub-01_ses-20241101_task-rest_run-01_space-MNI152NLin2009cAsym_desc-brain_mask.nii.gz",
+            "sha256:73f538f88fe28515012cadb5e19ee0c1feaf5da8f40de77794d028a218f820fd",
+        ),
+        "confounds": (
+            "sub-01_ses-20241101_task-rest_run-01_desc-confounds_timeseries.tsv",
+            "sha256:d8f6880b41b996f2a68485226d425629d0d3217c65d1aaeeb0cdfd3023346ae9",
+        ),
+    },
 }
 
 
@@ -192,7 +206,7 @@ def fetch_fmri_bold(session="0001", data_dir=None):
     Parameters
     ----------
     session : str
-        Session identifier. Available: ``"0001"``.
+        Session identifier. Available: ``"0001"``, ``"0008"``.
     data_dir : str or Path, optional
         Directory to store downloaded files. Default uses
         ``pooch``'s OS-appropriate cache directory.

gastropy/data/__init__.py

@@ -315,7 +315,7 @@ def bold_voxelwise_phases(
     return phases
 
 
-def compute_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indices=None):
+def compute_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indices=None, artifact_mask=None):
     """Compute a voxelwise PLV map between EGG and BOLD phases.
 
     Parameters
@@ -332,6 +332,9 @@ def compute_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indices=None):
     mask_indices : array_like, optional
         Boolean mask or integer indices mapping voxels back to the
         3D volume. Required if ``vol_shape`` is provided.
+    artifact_mask : array_like of bool, shape (n_timepoints,), optional
+        Boolean mask where ``True`` = clean volume, ``False`` = artifact.
+        Artifact volumes are excluded from the PLV computation.
 
     Returns
     -------
@@ -349,7 +352,7 @@ def compute_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indices=None):
         )
 
     # PLV: (n_timepoints, n_voxels) vs (n_timepoints,)
-    plv = phase_locking_value(bold_phases.T, egg_phase)
+    plv = phase_locking_value(bold_phases.T, egg_phase, mask=artifact_mask)
 
     if vol_shape is not None and mask_indices is not None:
         vol = np.zeros(vol_shape, dtype=float)
@@ -359,7 +362,7 @@ def compute_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indices=None):
     return plv
 
 
-def compute_surrogate_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indices=None, **kwargs):
+def compute_surrogate_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indices=None, artifact_mask=None, **kwargs):
     """Compute a surrogate PLV map via circular time-shifting.
 
     Same as ``compute_plv_map`` but uses surrogate PLV to generate
@@ -375,6 +378,9 @@ def compute_surrogate_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indic
         3D volume dimensions for reshaping output.
     mask_indices : array_like, optional
         Mask indices for 3D volume reconstruction.
+    artifact_mask : array_like of bool, shape (n_timepoints,), optional
+        Boolean mask where ``True`` = clean volume, ``False`` = artifact.
+        Artifact volumes are excluded from the surrogate PLV computation.
     **kwargs
         Additional arguments passed to ``surrogate_plv``
         (e.g., ``buffer_samples``, ``n_surrogates``, ``stat``, ``seed``).
@@ -398,7 +404,7 @@ def compute_surrogate_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indic
             f"bold_phases has {bold_phases.shape[1]}. They must match."
         )
 
-    surr = _surrogate_plv(bold_phases.T, egg_phase, **kwargs)
+    surr = _surrogate_plv(bold_phases.T, egg_phase, mask=artifact_mask, **kwargs)
 
     if vol_shape is not None and mask_indices is not None and surr.ndim == 1:
         vol = np.zeros(vol_shape, dtype=float)
@@ -408,6 +414,63 @@ def compute_surrogate_plv_map(egg_phase, bold_phases, vol_shape=None, mask_indic
     return surr
 
 
+def artifact_mask_to_volumes(sample_mask, trigger_times, sfreq, tr, begin_cut=21, end_cut=21):
+    """Map a sample-level artifact mask to volume-level.
+
+    Converts the per-sample boolean artifact mask (from
+    ``detect_phase_artifacts``) to a per-volume boolean mask
+    suitable for censoring bad volumes in PLV computation.
+
+    A volume is marked as artifact (``False``) if **any** EGG sample
+    within its time window is flagged in the sample mask.
+
+    Parameters
+    ----------
+    sample_mask : array_like of bool
+        Per-sample artifact mask where ``True`` = artifact sample.
+        This is the ``artifact_mask`` field from
+        ``detect_phase_artifacts``, at the EGG sampling rate.
+    trigger_times : array_like
+        Scanner trigger onset times in seconds (from
+        ``find_scanner_triggers``).
+    sfreq : float
+        EGG sampling frequency in Hz (e.g., 10.0).
+    tr : float
+        fMRI repetition time in seconds.
+    begin_cut : int, optional
+        Volumes to remove from the start. Default is 21.
+    end_cut : int, optional
+        Volumes to remove from the end. Default is 21.
+
+    Returns
+    -------
+    vol_mask : np.ndarray of bool
+        Per-volume mask where ``True`` = clean volume (include),
+        ``False`` = artifact volume (exclude). Length equals
+        ``n_volumes - begin_cut - end_cut``.
+    """
+    sample_mask = np.asarray(sample_mask, dtype=bool)
+    trigger_times = np.asarray(trigger_times, dtype=float)
+    n_samples = len(sample_mask)
+    n_volumes = len(trigger_times)
+
+    vol_clean = np.ones(n_volumes, dtype=bool)
+    for k in range(n_volumes):
+        t_start = trigger_times[k]
+        t_end = trigger_times[k + 1] if k + 1 < n_volumes else t_start + tr
+        s = min(int(round(t_start * sfreq)), n_samples)
+        e = min(int(round(t_end * sfreq)), n_samples)
+        if s < e and np.any(sample_mask[s:e]):
+            vol_clean[k] = False
+
+    # Apply volume cuts
+    if begin_cut + end_cut >= n_volumes:
+        return np.array([], dtype=bool)
+    if end_cut == 0:
+        return vol_clean[begin_cut:]
+    return vol_clean[begin_cut:-end_cut]
+
+
 def load_bold(bold_path, mask_path):
     """Load fMRIPrep BOLD and brain mask as a masked 2D array.
 
@@ -529,6 +592,7 @@ def to_nifti(data_3d, affine):
     "create_volume_windows",
     "phase_per_volume",
     "apply_volume_cuts",
+    "artifact_mask_to_volumes",
     "regress_confounds",
     "bold_voxelwise_phases",
     "compute_plv_map",

gastropy/data/fmri_egg_session_0008.npz

@@ -29,7 +29,7 @@
 
 # --- Config ---
 TARGET_SFREQ = 10.0  # Hz — standard for gastropy
-SEMI_PRECISION_SESSIONS = ["0001", "0003", "0004"]  # 3 baseline sessions
+SEMI_PRECISION_SESSIONS = ["0001", "0003", "0004", "0008"]
 TR = 1.856  # seconds