add nuisance regressor option to regression

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.
 

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}"
 

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.