diff --git a/src/rbc/core/metrics/timeseries.py b/src/rbc/core/metrics/timeseries.py index d4d60fbc..b5aa5c7b 100644 --- a/src/rbc/core/metrics/timeseries.py +++ b/src/rbc/core/metrics/timeseries.py @@ -6,11 +6,14 @@ from __future__ import annotations +import logging from pathlib import Path from typing import NamedTuple import numpy as np +logger = logging.getLogger(__name__) + def extract_timeseries( data: np.ndarray, @@ -64,6 +67,10 @@ def extract_timeseries( def correlation_matrix(timeseries: np.ndarray) -> np.ndarray: """Compute pairwise Pearson correlation matrix. + ROIs with zero variance (constant signal) produce ``NaN`` entries + because the Pearson correlation is undefined for constant series. + A warning is logged for any such ROIs. + Args: timeseries: (n_rois, n_timepoints) array. @@ -77,7 +84,27 @@ def correlation_matrix(timeseries: np.ndarray) -> np.ndarray: raise ValueError(f"Need at least 2 ROIs, got {timeseries.shape[0]}") if timeseries.shape[1] < 2: raise ValueError(f"Need at least 2 timepoints, got {timeseries.shape[1]}") - return np.corrcoef(timeseries) + + n_rois = timeseries.shape[0] + std = timeseries.std(axis=1) + valid = std > 0 + + if valid.all(): + return np.corrcoef(timeseries) + + n_zero = int((~valid).sum()) + logger.warning( + "%d of %d ROIs have zero variance (constant signal); " + "their correlation values will be NaN", + n_zero, + n_rois, + ) + + corr = np.full((n_rois, n_rois), np.nan, dtype=np.float64) + np.fill_diagonal(corr, 1.0) + if valid.sum() >= 2: + corr[np.ix_(valid, valid)] = np.corrcoef(timeseries[valid]) + return corr class TimeseriesOutputs(NamedTuple): diff --git a/tests/unit/core/test_timeseries.py b/tests/unit/core/test_timeseries.py index 2445c315..56d89ff2 100644 --- a/tests/unit/core/test_timeseries.py +++ b/tests/unit/core/test_timeseries.py @@ -190,6 +190,41 @@ def test_bounded(self) -> None: assert corr.min() >= -1.0 - 1e-10 assert corr.max() <= 1.0 + 1e-10 + def test_zero_variance_roi_produces_nan(self) -> None: + """Zero-variance ROIs should yield NaN without RuntimeWarning.""" + ts = np.array( + [ + [1.0, 2.0, 3.0, 4.0, 5.0], + [7.0, 7.0, 7.0, 7.0, 7.0], # constant + [5.0, 4.0, 3.0, 2.0, 1.0], + ] + ) + import warnings + + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) + corr = correlation_matrix(ts) + + # Valid ROIs should have real correlation values + assert corr[0, 2] == pytest.approx(-1.0) + # Zero-variance ROI should produce NaN + assert np.isnan(corr[1, 0]) + assert np.isnan(corr[0, 1]) + assert np.isnan(corr[1, 2]) + assert corr[1, 1] == pytest.approx(1.0) + + def test_all_zero_variance_produces_all_nan(self) -> None: + """If all ROIs are constant, the entire matrix should be NaN.""" + ts = np.array( + [ + [5.0, 5.0, 5.0], + [3.0, 3.0, 3.0], + ] + ) + corr = correlation_matrix(ts) + assert np.all(np.isnan(corr[~np.eye(2, dtype=bool)])) + np.testing.assert_allclose(np.diag(corr), 1.0) + def test_rejects_single_roi(self) -> None: """Should reject fewer than 2 ROIs.""" with pytest.raises(ValueError, match="2 ROIs"):