Skip to content

Commit f496f9d

Browse files
committed
[FIX] ASR: round half away from zero for window/bad-chan counts
numpy's np.round (banker's rounding) and int() truncation diverge from MATLAB's round-half-away-from-zero at exact-half inputs, giving +-1 counts for the bad-channel count and window length N in clean_windows and asr_calibrate. Add a _round_half_away helper and apply it at those sites.
1 parent 107cd1b commit f496f9d

2 files changed

Lines changed: 41 additions & 4 deletions

File tree

meegkit/asr.py

Lines changed: 12 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -270,6 +270,15 @@ def transform(self, X, y=None, **kwargs):
270270
return out
271271

272272

273+
def _round_half_away(x):
274+
"""Round half away from zero (MATLAB `round` semantics).
275+
276+
numpy's ``np.round`` rounds half to even (banker's rounding); this
277+
matches MATLAB/Octave which round halves away from zero.
278+
"""
279+
return np.floor(np.abs(x) + 0.5) * np.sign(x)
280+
281+
273282
def clean_windows(X, sfreq, max_bad_chans=0.2, zthresholds=[-3.5, 5],
274283
win_len=.5, win_overlap=0.66, min_clean_fraction=0.25,
275284
max_dropout_fraction=0.1, show=False):
@@ -336,11 +345,11 @@ def clean_windows(X, sfreq, max_bad_chans=0.2, zthresholds=[-3.5, 5],
336345
truncate_quant = [0.0220, 0.6000]
337346
step_sizes = [0.01, 0.01]
338347
shape_range = SHAPE_RANGE
339-
max_bad_chans = np.round(X.shape[0] * max_bad_chans)
348+
max_bad_chans = _round_half_away(X.shape[0] * max_bad_chans)
340349

341350
# set data indices
342351
[nc, ns] = X.shape
343-
N = int(win_len * sfreq)
352+
N = int(_round_half_away(win_len * sfreq))
344353
N_raw = win_len * sfreq # non-truncated N, avoids step-size phase drift
345354
offsets = np.round(np.arange(0, ns - N, N_raw * (1 - win_overlap)))
346355
offsets = offsets.astype(int)
@@ -502,7 +511,7 @@ def asr_calibrate(X, sfreq, cutoff=5, blocksize=100, win_len=0.5,
502511
X, _zf = yulewalk_filter(X, sfreq, ab=None)
503512

504513
# window length for calculating thresholds
505-
N = int(np.round(win_len * sfreq))
514+
N = int(_round_half_away(win_len * sfreq))
506515
if ns < N:
507516
raise ValueError(
508517
f"Calibration data has {ns} samples, shorter than one analysis "

tests/test_asr.py

Lines changed: 29 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@
88
import pytest
99
from scipy import signal
1010

11-
from meegkit.asr import ASR, asr_calibrate, asr_process, clean_windows
11+
from meegkit.asr import ASR, _round_half_away, asr_calibrate, asr_process, clean_windows
1212
from meegkit.utils.asr import SHAPE_RANGE, fit_eeg_distribution, yulewalk, yulewalk_filter
1313
from meegkit.utils.matrix import sliding_window
1414

@@ -141,6 +141,34 @@ def test_yulewalk_filter(n_chans, show=False):
141141
plt.show()
142142

143143

144+
def test_round_half_away():
145+
"""_round_half_away rounds ties away from zero, unlike numpy's banker's rounding."""
146+
# exact-half (tie) behavior: rounds away from zero
147+
assert _round_half_away(2.5) == 3
148+
assert _round_half_away(128.5) == 129
149+
assert _round_half_away(-2.5) == -3
150+
assert _round_half_away(0.5) == 1
151+
assert _round_half_away(1.5) == 2 # np.round gives 2 here too (even)
152+
153+
# divergence from banker's rounding on ties
154+
assert _round_half_away(2.5) != np.round(2.5) # 3 != 2
155+
assert _round_half_away(128.5) != np.round(128.5) # 129 != 128
156+
157+
# non-tie values round normally
158+
assert _round_half_away(2.4) == 2
159+
assert _round_half_away(2.6) == 3
160+
assert _round_half_away(-2.6) == -3
161+
162+
# int() of the result is exact for scalars
163+
assert int(_round_half_away(2.5)) == 3
164+
assert int(_round_half_away(-2.5)) == -3
165+
166+
# array input
167+
ties = np.array([2.5, 128.5, -2.5, 0.5, 1.5])
168+
expected = np.array([3, 129, -3, 1, 2])
169+
np.testing.assert_array_equal(_round_half_away(ties), expected)
170+
171+
144172
def test_asr_functions(show=False, method="riemann"):
145173
"""Test ASR functions (offline use).
146174

0 commit comments

Comments
 (0)