From 2d51b2897f50714c56b932612732d1a4d9ed6b24 Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Tue, 18 Aug 2026 14:49:37 +0100
Subject: [PATCH] add nuisance regressor option to regression

---
 src/mesoscopy/io.py                 |  56 +++++++
 src/mesoscopy/process/__init__.py   |  41 ++++-
 src/mesoscopy/process/regression.py |  96 +++++++++++
 tests/test_process.py               | 241 ++++++++++++++++++++++++++++
 4 files changed, 431 insertions(+), 3 deletions(-)

diff --git a/src/mesoscopy/io.py b/src/mesoscopy/io.py
index 6daf35c..8a87360 100644
--- a/src/mesoscopy/io.py
+++ b/src/mesoscopy/io.py
@@ -261,6 +261,62 @@ def _read_hdf5_regressors(path: str) -> tuple[np.ndarray, list[str], np.ndarray]
     return regressors, labels, trial_indices
 
 
+def read_nuisance_regressors(path: str) -> tuple[np.ndarray, list[str], np.ndarray]:
+    """Read an external nuisance regressor file in NPZ or HDF5 format.
+
+    Args:
+        path (str): Path to the nuisance regressor file.
+
+    Returns:
+        tuple[np.ndarray, list[str], np.ndarray]: Nuisance regressor matrix of shape (n_samples, n_regressors),
+        list of regressor labels, and the timestamps (n_samples,) the regressors were recorded at.
+
+    Raises:
+        ValueError: If the file format is unsupported.
+    """
+    if path.endswith(".npz"):
+        return _read_npz_nuisance_regressors(path)
+    if path.endswith(".h5"):
+        return _read_hdf5_nuisance_regressors(path)
+    msg = "Unsupported file format."
+    raise ValueError(msg)
+
+
+def _read_npz_nuisance_regressors(path: str) -> tuple[np.ndarray, list[str], np.ndarray]:
+    """Read a nuisance regressor file in NPZ format.
+
+    Args:
+        path (str): Path to the NPZ file. Every array other than 'timestamps' is treated as a nuisance regressor.
+
+    Returns:
+        tuple[np.ndarray, list[str], np.ndarray]: Nuisance regressor matrix, list of regressor labels, and
+        timestamps.
+    """
+    with np.load(path) as f:
+        labels = sorted(key for key in f.files if key != "timestamps")
+        regressors = np.column_stack([f[label] for label in labels])
+        timestamps = np.array(f["timestamps"])
+    return regressors, labels, timestamps
+
+
+def _read_hdf5_nuisance_regressors(path: str) -> tuple[np.ndarray, list[str], np.ndarray]:
+    """Read a nuisance regressor file in HDF5 format.
+
+    Args:
+        path (str): Path to the HDF5 file. Every dataset other than 'timestamps' is treated as a nuisance
+            regressor.
+
+    Returns:
+        tuple[np.ndarray, list[str], np.ndarray]: Nuisance regressor matrix, list of regressor labels, and
+        timestamps.
+    """
+    with h5py.File(path, "r") as f:
+        labels = sorted(key for key in f if key != "timestamps")
+        regressors = np.column_stack([f[label][:] for label in labels])
+        timestamps = np.array(f["timestamps"][:])
+    return regressors, labels, timestamps
+
+
 def _read_fiji_points(path: str) -> dict[str, tuple[float, float]]:
     """Read a FIJI landmark points file.
 
diff --git a/src/mesoscopy/process/__init__.py b/src/mesoscopy/process/__init__.py
index b9cb209..8366f53 100644
--- a/src/mesoscopy/process/__init__.py
+++ b/src/mesoscopy/process/__init__.py
@@ -62,7 +62,7 @@ def process_cmd(): ...
 )
 def smooth_cmd(path: str, out_dir: str, sigma: int = 2) -> None:
     """Generate a smoothed DeltaF/F recording using a Laplace of Gaussian filter."""
-    if not os.path.exists(out_dir):
+    if not Path(out_dir).exists():
         click.echo(f"Creating output directory {out_dir}...")
         Path(out_dir).mkdir(parents=True)
 
@@ -207,6 +207,19 @@ def regions_cmd(path: str, out_dir: str) -> None:
     default=1.0,
     help="Ridge regularisation strength. Defaults to 1.0.",
 )
+@click.option(
+    "-n",
+    "--nuisance-regressors",
+    "nuisance_regressor_paths",
+    type=click.Path(exists=True),
+    multiple=True,
+    help=(
+        "Path to an external nuisance regressor file (NPZ or HDF5), e.g. behavioural motion energy. Every"
+        " array/dataset in the file other than 'timestamps' is treated as one nuisance regressor and interpolated"
+        " onto the recording's own timestamps before being z-scored and appended to the regressor matrix. May be"
+        " passed multiple times to add nuisance regressors from several files."
+    ),
+)
 @click.option(
     "-f",
     "--fast",
@@ -231,7 +244,13 @@ def regions_cmd(path: str, out_dir: str) -> None:
     help="Save regression results as an HDF5 file. Defaults to False.",
 )
 def regression_cmd(
-    recording_path: str, regressor_path: str, out_dir: str, alpha: float, fast: bool, file_format: str
+    recording_path: str,
+    regressor_path: str,
+    out_dir: str,
+    alpha: float,
+    nuisance_regressor_paths: tuple[str, ...],
+    fast: bool,
+    file_format: str,
 ) -> None:
     """Perform pixel-wise ridge regression on a preprocessed ∆F/F recording."""
     if not Path(out_dir).exists():
@@ -241,10 +260,26 @@ def regression_cmd(
     click.echo(f"Loading preprocessed recording from {recording_path}...")
     # Determine whether we're working with an NWB file
     nwb = bool(recording_path.endswith(".nwb"))
-    session_id, deltaf_series, _ = io.load_deltaf(recording_path, nwb=nwb)
+    session_id, deltaf_series, timestamps = io.load_deltaf(recording_path, nwb=nwb)
 
     click.echo(f"Loading regressors from {regressor_path}...")
     regressors, labels, trial_idx = io.read_regressors(regressor_path)
+    labels = list(labels)
+
+    if nuisance_regressor_paths:
+        # Nuisance regressors are recorded on their own clock (e.g. a behavioural camera), so both series are
+        # anchored to elapsed seconds since their own first sample before interpolating one onto the other -- see
+        # `regr.elapsed_seconds`.
+        target_timestamps = regr.elapsed_seconds(timestamps)
+
+    for nuisance_path in nuisance_regressor_paths:
+        click.echo(f"Loading nuisance regressors from {nuisance_path}...")
+        nuisance, nuisance_labels, nuisance_timestamps = io.read_nuisance_regressors(nuisance_path)
+        nuisance = regr.interpolate_regressors(nuisance, regr.elapsed_seconds(nuisance_timestamps), target_timestamps)
+        if trial_idx is not None:
+            nuisance = nuisance[trial_idx]
+        regressors, labels = regr.append_nuisance_regressors(regressors, labels, nuisance, nuisance_labels)
+        click.echo(f"Added nuisance regressors: {nuisance_labels}")
 
     outpath = out_dir + os.sep + session_id + f"_regression.{file_format}"
 
diff --git a/src/mesoscopy/process/regression.py b/src/mesoscopy/process/regression.py
index 08126e9..7f73fe3 100644
--- a/src/mesoscopy/process/regression.py
+++ b/src/mesoscopy/process/regression.py
@@ -19,6 +19,8 @@
 #  IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 #  SOFTWARE.
 
+from datetime import datetime
+
 import numpy as np
 from sklearn.linear_model import Ridge
 from sklearn.metrics import mean_squared_error
@@ -112,6 +114,100 @@ def ridge_regression_fast(deltaf_series: np.ndarray, regressors: np.ndarray, alp
     return coefficients, r2, mse
 
 
+def elapsed_seconds(timestamps: np.ndarray) -> np.ndarray:
+    """Convert a timestamps array to elapsed seconds since its own first sample.
+
+    Handles both already-numeric timestamps (e.g. from a synthetic/plain-float recording) and the ISO-8601
+    byte/string timestamps written by mesoscopy's own preprocessing stage (the ``/timestamps`` dataset in
+    ``*_preprocessed.h5`` files, which is also what ends up linked into NWB output). This lets two
+    independently-timestamped series (e.g. a dF/F recording and an external nuisance regressor file) be
+    interpolated onto a common timebase, on the assumption that both were started at roughly the same moment
+    (e.g. by a shared trigger) -- each is anchored to its own first sample rather than to any shared absolute
+    clock.
+
+    Args:
+        timestamps (np.ndarray): Timestamps, either numeric (seconds) or ISO-8601 byte/str strings.
+
+    Returns:
+        np.ndarray: Elapsed seconds since the first timestamp, as a float array.
+    """
+    timestamps = np.asarray(timestamps)
+    if timestamps.dtype.kind in "SU":
+        parsed = [datetime.fromisoformat(ts.decode("utf-8") if isinstance(ts, bytes) else ts) for ts in timestamps]
+        return np.array([(ts - parsed[0]).total_seconds() for ts in parsed])
+
+    seconds = timestamps.astype(float)
+    return seconds - seconds[0]
+
+
+def interpolate_regressors(
+    regressors: np.ndarray, source_timestamps: np.ndarray, target_timestamps: np.ndarray
+) -> np.ndarray:
+    """Interpolate a regressor matrix from its own timebase onto a target timebase.
+
+    Used to align external nuisance regressors onto  a recording's dF/F timestamps. Interpolation is linear
+    and independent per column; target timestamps outside the range of ``source_timestamps`` are clamped to the nearest
+    edge value (`numpy.interp`'s default behaviour).
+
+    Args:
+        regressors (np.ndarray): Regressor matrix of shape (n_source_samples, n_regressors).
+        source_timestamps (np.ndarray): Timestamps the regressors were recorded at, shape (n_source_samples,).
+        target_timestamps (np.ndarray): Timestamps to interpolate onto, shape (n_target_samples,).
+
+    Returns:
+        np.ndarray: Interpolated regressor matrix of shape (n_target_samples, n_regressors).
+    """
+    return np.column_stack(
+        [np.interp(target_timestamps, source_timestamps, column) for column in regressors.T],
+    )
+
+
+def append_nuisance_regressors(
+    regressors: np.ndarray,
+    labels: list[str],
+    nuisance: np.ndarray,
+    nuisance_labels: list[str],
+    zscore: bool = True,
+) -> tuple[np.ndarray, list[str]]:
+    """Append external nuisance regressors to an existing regressor matrix.
+
+    Args:
+        regressors (np.ndarray): Existing regressor matrix of shape (n_samples, n_regressors).
+        labels (list[str]): Labels of the existing regressors.
+        nuisance (np.ndarray): Nuisance regressor matrix of shape (n_samples, n_nuisance_regressors), already
+            aligned onto the same timebase as ``regressors`` (see :func:`interpolate_regressors`).
+        nuisance_labels (list[str]): Labels of the nuisance regressors.
+        zscore (bool): Whether to z-score each nuisance column (zero mean, unit variance) before appending, so
+            that ridge regularisation isn't skewed by the nuisance regressors' raw scale. Defaults to True.
+
+    Returns:
+        tuple[np.ndarray, list[str]]: Combined regressor matrix, and combined list of labels.
+
+    Raises:
+        ValueError: If the number of samples doesn't match between ``regressors`` and ``nuisance``, or if any
+            nuisance label duplicates an existing label.
+    """
+    if nuisance.shape[0] != regressors.shape[0]:
+        msg = (
+            f"Nuisance regressors have {nuisance.shape[0]} samples, but existing regressors have "
+            f"{regressors.shape[0]} samples."
+        )
+        raise ValueError(msg)
+
+    duplicate_labels = set(labels) & set(nuisance_labels)
+    if duplicate_labels:
+        msg = f"Nuisance regressor labels {sorted(duplicate_labels)} duplicate existing regressor labels."
+        raise ValueError(msg)
+
+    if zscore:
+        nuisance = (nuisance - nuisance.mean(axis=0)) / nuisance.std(axis=0)
+
+    combined_regressors = np.column_stack([regressors, nuisance])
+    combined_labels = [*labels, *nuisance_labels]
+
+    return combined_regressors, combined_labels
+
+
 def _pixel_ridge_regression(deltaf_series: np.ndarray, regressors: np.ndarray) -> np.ndarray:
     """Perform ridge regression on a single pixel's DeltaF/F series.
 
diff --git a/tests/test_process.py b/tests/test_process.py
index dbf7408..0bece51 100644
--- a/tests/test_process.py
+++ b/tests/test_process.py
@@ -73,6 +73,39 @@ def regressor_h5(tmp_path_factory):
     return str(tmpfile)
 
 
+@pytest.fixture
+def nuisance_regressor_h5(tmp_path_factory):
+    """Create an HDF5 nuisance regressor file on its own (coarser) timebase, spanning the same duration as the
+    (300, 40, 40) preproc_h5 fixture's timestamps (0 to 0.299 seconds)."""
+    tmpfile = tmp_path_factory.mktemp("data") / "nuisance.h5"
+    rng = np.random.default_rng(3)
+    n_samples = 50
+    timestamps = np.linspace(0, 0.3, n_samples)
+    with h5.File(str(tmpfile), "w") as f:
+        f.create_dataset("motion_a", data=rng.normal(10, 2, size=n_samples))
+        f.create_dataset("motion_b", data=rng.normal(-5, 1, size=n_samples))
+        f.create_dataset("timestamps", data=timestamps)
+    return str(tmpfile)
+
+
+@pytest.fixture
+def nuisance_regressor_npz(tmp_path_factory):
+    """Create an NPZ nuisance regressor file on its own (coarser) timebase, spanning the same duration as the
+    (300, 40, 40) preproc_h5 fixture's timestamps (0 to 0.299 seconds). Uses different labels from
+    `nuisance_regressor_h5` so both fixtures can be combined in the same regression run."""
+    tmpfile = tmp_path_factory.mktemp("data") / "nuisance.npz"
+    rng = np.random.default_rng(4)
+    n_samples = 50
+    timestamps = np.linspace(0, 0.3, n_samples)
+    np.savez(
+        tmpfile,
+        pupil_a=rng.normal(10, 2, size=n_samples),
+        pupil_b=rng.normal(-5, 1, size=n_samples),
+        timestamps=timestamps,
+    )
+    return str(tmpfile)
+
+
 @pytest.fixture
 def preproc_h5_bytes_timestamps(tmp_path_factory):
     """Preprocessed HDF5 file with byte-string timestamps, matching the real preprocessing pipeline's
@@ -261,6 +294,135 @@ def test_zero_variance_target_gives_perfect_r2(self):
         assert mse[0, 0] == pytest.approx(0.0, abs=1e-8)
 
 
+# ---------------------------------------------------------------------------
+# regression.elapsed_seconds
+# ---------------------------------------------------------------------------
+
+
+class TestElapsedSeconds:
+    def test_numeric_timestamps_normalized_to_start_at_zero(self):
+        timestamps = np.array([5.0, 6.0, 8.0])
+        result = regr.elapsed_seconds(timestamps)
+        np.testing.assert_allclose(result, [0.0, 1.0, 3.0])
+
+    def test_already_zero_started_numeric_timestamps_unchanged(self):
+        timestamps = np.array([0.0, 0.5, 1.0])
+        result = regr.elapsed_seconds(timestamps)
+        np.testing.assert_allclose(result, timestamps)
+
+    def test_parses_iso_byte_string_timestamps(self):
+        """Matches the /timestamps format written by the real preprocessing pipeline (dtype 'S25')."""
+        timestamps = np.array(
+            [b"2024-01-01T14:00:00.000", b"2024-01-01T14:00:00.500", b"2024-01-01T14:00:01.000"], dtype="S25"
+        )
+        result = regr.elapsed_seconds(timestamps)
+        np.testing.assert_allclose(result, [0.0, 0.5, 1.0])
+
+    def test_parses_iso_str_timestamps(self):
+        timestamps = np.array(["2024-01-01T14:00:00.000", "2024-01-01T14:00:00.500"])
+        result = regr.elapsed_seconds(timestamps)
+        np.testing.assert_allclose(result, [0.0, 0.5])
+
+
+# ---------------------------------------------------------------------------
+# regression.interpolate_regressors
+# ---------------------------------------------------------------------------
+
+
+class TestInterpolateRegressors:
+    def test_output_shape(self):
+        regressors = np.arange(20).reshape(10, 2).astype(float)
+        source_timestamps = np.linspace(0, 1, 10)
+        target_timestamps = np.linspace(0, 1, 25)
+        result = regr.interpolate_regressors(regressors, source_timestamps, target_timestamps)
+        assert result.shape == (25, 2)
+
+    def test_matches_linear_interpolation(self):
+        source_timestamps = np.array([0.0, 1.0, 2.0, 3.0])
+        regressors = np.array([[0.0], [10.0], [20.0], [30.0]])
+        target_timestamps = np.array([0.5, 1.5, 2.5])
+        result = regr.interpolate_regressors(regressors, source_timestamps, target_timestamps)
+        np.testing.assert_allclose(result, [[5.0], [15.0], [25.0]])
+
+    def test_clamps_out_of_range_targets(self):
+        source_timestamps = np.array([1.0, 2.0, 3.0])
+        regressors = np.array([[10.0], [20.0], [30.0]])
+        target_timestamps = np.array([-5.0, 10.0])
+        result = regr.interpolate_regressors(regressors, source_timestamps, target_timestamps)
+        np.testing.assert_allclose(result, [[10.0], [30.0]])
+
+
+# ---------------------------------------------------------------------------
+# regression.append_nuisance_regressors
+# ---------------------------------------------------------------------------
+
+
+class TestAppendNuisanceRegressors:
+    def test_output_shape_and_labels(self):
+        regressors = np.random.default_rng(0).normal(size=(20, 2))
+        nuisance = np.random.default_rng(1).normal(size=(20, 3))
+        combined, labels = regr.append_nuisance_regressors(regressors, ["a", "b"], nuisance, ["n1", "n2", "n3"])
+        assert combined.shape == (20, 5)
+        assert labels == ["a", "b", "n1", "n2", "n3"]
+
+    def test_zscores_nuisance_columns_by_default(self):
+        regressors = np.zeros((20, 1))
+        nuisance = np.random.default_rng(0).normal(loc=100, scale=10, size=(20, 2))
+        combined, _ = regr.append_nuisance_regressors(regressors, ["a"], nuisance, ["n1", "n2"])
+        appended = combined[:, 1:]
+        np.testing.assert_allclose(appended.mean(axis=0), 0, atol=1e-10)
+        np.testing.assert_allclose(appended.std(axis=0), 1, atol=1e-10)
+
+    def test_zscore_false_leaves_nuisance_raw(self):
+        regressors = np.zeros((20, 1))
+        nuisance = np.random.default_rng(0).normal(loc=100, scale=10, size=(20, 2))
+        combined, _ = regr.append_nuisance_regressors(regressors, ["a"], nuisance, ["n1", "n2"], zscore=False)
+        np.testing.assert_allclose(combined[:, 1:], nuisance)
+
+    def test_leaves_original_regressors_untouched(self):
+        regressors = np.random.default_rng(0).normal(size=(20, 1))
+        nuisance = np.random.default_rng(1).normal(loc=100, scale=10, size=(20, 1))
+        combined, _ = regr.append_nuisance_regressors(regressors, ["a"], nuisance, ["n1"])
+        np.testing.assert_allclose(combined[:, :1], regressors)
+
+    def test_raises_on_sample_count_mismatch(self):
+        regressors = np.zeros((20, 1))
+        nuisance = np.zeros((10, 1))
+        with pytest.raises(ValueError, match="samples"):
+            regr.append_nuisance_regressors(regressors, ["a"], nuisance, ["n1"])
+
+    def test_raises_on_duplicate_labels(self):
+        regressors = np.zeros((20, 1))
+        nuisance = np.zeros((20, 1))
+        with pytest.raises(ValueError, match="duplicate"):
+            regr.append_nuisance_regressors(regressors, ["a"], nuisance, ["a"])
+
+
+# ---------------------------------------------------------------------------
+# io.read_nuisance_regressors
+# ---------------------------------------------------------------------------
+
+
+class TestReadNuisanceRegressors:
+    def test_reads_hdf5(self, nuisance_regressor_h5):
+        regressors, labels, timestamps = io.read_nuisance_regressors(nuisance_regressor_h5)
+        assert regressors.shape == (50, 2)
+        assert labels == ["motion_a", "motion_b"]
+        assert timestamps.shape == (50,)
+
+    def test_reads_npz(self, nuisance_regressor_npz):
+        regressors, labels, timestamps = io.read_nuisance_regressors(nuisance_regressor_npz)
+        assert regressors.shape == (50, 2)
+        assert labels == ["pupil_a", "pupil_b"]
+        assert timestamps.shape == (50,)
+
+    def test_raises_on_unsupported_format(self, tmp_path):
+        path = tmp_path / "nuisance.txt"
+        path.touch()
+        with pytest.raises(ValueError, match="Unsupported"):
+            io.read_nuisance_regressors(str(path))
+
+
 # ---------------------------------------------------------------------------
 # CLI commands
 # ---------------------------------------------------------------------------
@@ -377,6 +539,85 @@ def test_regression_cmd_not_fast(preproc_h5, regressor_npz, output_dir):
     assert outpath.is_file()
 
 
+def test_regression_cmd_with_nuisance_regressors(preproc_h5, regressor_npz, nuisance_regressor_h5, output_dir):
+    """Nuisance regressors from an external file should be interpolated, z-scored, and appended."""
+    runner = CliRunner()
+    result = runner.invoke(
+        mesoscopy.cli,
+        args=f"process regression {preproc_h5} {regressor_npz} -o {output_dir} --fast -n {nuisance_regressor_h5}",
+    )
+    assert result.exit_code == 0
+    assert "Added nuisance regressors: ['motion_a', 'motion_b']" in result.output
+
+    outpath = pathlib.Path(output_dir) / "preproc_regression.npz"
+    assert outpath.is_file()
+
+    with np.load(outpath, allow_pickle=True) as f:
+        # 3 task regressors (reg_a/b/c) + 2 nuisance regressors (motion_a/b).
+        assert f["coefficients"].shape == (5, 40, 40)
+        assert list(f["labels"])[-2:] == ["motion_a", "motion_b"]
+
+
+def test_regression_cmd_with_multiple_nuisance_regressor_files(
+    preproc_h5, regressor_npz, nuisance_regressor_h5, nuisance_regressor_npz, output_dir
+):
+    """Passing -n multiple times should combine nuisance regressors from every file."""
+    runner = CliRunner()
+    result = runner.invoke(
+        mesoscopy.cli,
+        args=(
+            f"process regression {preproc_h5} {regressor_npz} -o {output_dir} --fast"
+            f" -n {nuisance_regressor_h5} -n {nuisance_regressor_npz}"
+        ),
+    )
+    assert result.exit_code == 0
+
+    outpath = pathlib.Path(output_dir) / "preproc_regression.npz"
+    with np.load(outpath, allow_pickle=True) as f:
+        # 3 task regressors + 2 nuisance regressors from each of the two files.
+        assert f["coefficients"].shape == (7, 40, 40)
+
+
+def test_regression_cmd_with_nuisance_regressors_and_trial_idx(
+    preproc_h5, regressor_h5, nuisance_regressor_h5, output_dir
+):
+    """Nuisance regressors should be interpolated onto the full recording, then subset by trial_idx to match the
+    (already trial_idx-subset) task regressor matrix."""
+    runner = CliRunner()
+    result = runner.invoke(
+        mesoscopy.cli,
+        args=f"process regression {preproc_h5} {regressor_h5} -o {output_dir} --fast --h5 -n {nuisance_regressor_h5}",
+    )
+    assert result.exit_code == 0
+
+    outpath = pathlib.Path(output_dir) / "preproc_regression.h5"
+    result_h5 = io.read_h5(str(outpath))
+    assert result_h5["/coefficients"].shape == (5, 40, 40)
+    assert result_h5["/trial_idx"].shape == (300,)
+
+
+def test_regression_cmd_with_nuisance_regressors_iso_timestamps(
+    preproc_h5_bytes_timestamps, regressor_npz, nuisance_regressor_h5, output_dir
+):
+    """Real preprocessed recordings store /timestamps as ISO-8601 byte strings, not the plain floats used by the
+    `preproc_h5` fixture -- make sure nuisance regressor alignment handles that format too."""
+    runner = CliRunner()
+    result = runner.invoke(
+        mesoscopy.cli,
+        args=(
+            f"process regression {preproc_h5_bytes_timestamps} {regressor_npz} -o {output_dir} --fast"
+            f" -n {nuisance_regressor_h5}"
+        ),
+    )
+    assert result.exit_code == 0
+
+    outpath = pathlib.Path(output_dir) / "preproc_bytes_regression.npz"
+    assert outpath.is_file()
+
+    with np.load(outpath, allow_pickle=True) as f:
+        assert f["coefficients"].shape == (5, 40, 40)
+
+
 # ---------------------------------------------------------------------------
 # region.extract_region_activity / region.extract_all_regions fixtures
 # ---------------------------------------------------------------------------