Handle zero-variance ROIs in correlation matrix

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):

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"):