Skip to content

Commit 53cc5be

Browse files
committed
Handle zero-variance ROIs in correlation matrix computation
Skip zero-variance ROIs when computing np.corrcoef to avoid the RuntimeWarning about division by zero. Affected off-diagonal entries are set to NaN; the diagonal is always 1.0. A warning is logged identifying how many ROIs were constant. Closes #261
1 parent 8837200 commit 53cc5be

2 files changed

Lines changed: 63 additions & 1 deletion

File tree

src/rbc/core/metrics/timeseries.py

Lines changed: 28 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,11 +6,14 @@
66

77
from __future__ import annotations
88

9+
import logging
910
from pathlib import Path
1011
from typing import NamedTuple
1112

1213
import numpy as np
1314

15+
logger = logging.getLogger(__name__)
16+
1417

1518
def extract_timeseries(
1619
data: np.ndarray,
@@ -64,6 +67,10 @@ def extract_timeseries(
6467
def correlation_matrix(timeseries: np.ndarray) -> np.ndarray:
6568
"""Compute pairwise Pearson correlation matrix.
6669
70+
ROIs with zero variance (constant signal) produce ``NaN`` entries
71+
because the Pearson correlation is undefined for constant series.
72+
A warning is logged for any such ROIs.
73+
6774
Args:
6875
timeseries: (n_rois, n_timepoints) array.
6976
@@ -77,7 +84,27 @@ def correlation_matrix(timeseries: np.ndarray) -> np.ndarray:
7784
raise ValueError(f"Need at least 2 ROIs, got {timeseries.shape[0]}")
7885
if timeseries.shape[1] < 2:
7986
raise ValueError(f"Need at least 2 timepoints, got {timeseries.shape[1]}")
80-
return np.corrcoef(timeseries)
87+
88+
n_rois = timeseries.shape[0]
89+
std = timeseries.std(axis=1)
90+
valid = std > 0
91+
92+
if valid.all():
93+
return np.corrcoef(timeseries)
94+
95+
n_zero = int((~valid).sum())
96+
logger.warning(
97+
"%d of %d ROIs have zero variance (constant signal); "
98+
"their correlation values will be NaN",
99+
n_zero,
100+
n_rois,
101+
)
102+
103+
corr = np.full((n_rois, n_rois), np.nan, dtype=np.float64)
104+
np.fill_diagonal(corr, 1.0)
105+
if valid.sum() >= 2:
106+
corr[np.ix_(valid, valid)] = np.corrcoef(timeseries[valid])
107+
return corr
81108

82109

83110
class TimeseriesOutputs(NamedTuple):

tests/unit/core/test_timeseries.py

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -190,6 +190,41 @@ def test_bounded(self) -> None:
190190
assert corr.min() >= -1.0 - 1e-10
191191
assert corr.max() <= 1.0 + 1e-10
192192

193+
def test_zero_variance_roi_produces_nan(self) -> None:
194+
"""Zero-variance ROIs should yield NaN without RuntimeWarning."""
195+
ts = np.array(
196+
[
197+
[1.0, 2.0, 3.0, 4.0, 5.0],
198+
[7.0, 7.0, 7.0, 7.0, 7.0], # constant
199+
[5.0, 4.0, 3.0, 2.0, 1.0],
200+
]
201+
)
202+
import warnings
203+
204+
with warnings.catch_warnings():
205+
warnings.simplefilter("error", RuntimeWarning)
206+
corr = correlation_matrix(ts)
207+
208+
# Valid ROIs should have real correlation values
209+
assert corr[0, 2] == pytest.approx(-1.0)
210+
# Zero-variance ROI should produce NaN
211+
assert np.isnan(corr[1, 0])
212+
assert np.isnan(corr[0, 1])
213+
assert np.isnan(corr[1, 2])
214+
assert corr[1, 1] == pytest.approx(1.0)
215+
216+
def test_all_zero_variance_produces_all_nan(self) -> None:
217+
"""If all ROIs are constant, the entire matrix should be NaN."""
218+
ts = np.array(
219+
[
220+
[5.0, 5.0, 5.0],
221+
[3.0, 3.0, 3.0],
222+
]
223+
)
224+
corr = correlation_matrix(ts)
225+
assert np.all(np.isnan(corr[~np.eye(2, dtype=bool)]))
226+
np.testing.assert_allclose(np.diag(corr), 1.0)
227+
193228
def test_rejects_single_roi(self) -> None:
194229
"""Should reject fewer than 2 ROIs."""
195230
with pytest.raises(ValueError, match="2 ROIs"):

0 commit comments

Comments
 (0)