From bce15c973799117afc1ab924c26bacad37096108 Mon Sep 17 00:00:00 2001 From: Constantinos Eleftheriou Date: Mon, 12 Oct 2026 15:37:57 +0100 Subject: [PATCH 1/5] vectorise population and rolling LDA decoder --- src/mesoscopy/process/decoding.py | 135 +++++++++++++++++++++++++++++- tests/test_process.py | 94 +++++++++++++++++++++ 2 files changed, 227 insertions(+), 2 deletions(-) diff --git a/src/mesoscopy/process/decoding.py b/src/mesoscopy/process/decoding.py index 8d7144a..9560af1 100644 --- a/src/mesoscopy/process/decoding.py +++ b/src/mesoscopy/process/decoding.py @@ -413,11 +413,21 @@ def _standardised(x: npt.NDArray[np.float64], train: npt.NDArray[np.float64]) -> mean = np.einsum("pfn,nr->pfr", train, x) / n centred = x[None, None] - mean[:, :, None, :] var = np.einsum("pfn,pfnr->pfr", train, centred**2) / n - # StandardScaler's rule for a feature that is constant up to rounding. + return centred / _scale(var, mean, n)[:, :, None, :] + + +def _scale( + var: npt.NDArray[np.float64], mean: npt.NDArray[np.float64], n: npt.NDArray[np.float64] +) -> npt.NDArray[np.float64]: + """Standard deviation as `StandardScaler` takes it, 1 for a feature that is constant up to rounding. + + Returns: + npt.NDArray[np.float64]: Same shape as `var`. + """ eps = np.finfo(np.float64).eps sd = np.sqrt(var) sd[var <= n * eps * var + (n * mean * eps) ** 2] = 1.0 - return centred / sd[:, :, None, :] + return sd def _test_fold(values: npt.NDArray[np.float64], fold_of: npt.NDArray[np.int64]) -> npt.NDArray[np.float64]: @@ -535,6 +545,75 @@ def objective(w: npt.NDArray[np.float64], b: npt.NDArray[np.float64]) -> npt.NDA return xs * w[:, :, None, :] + b[:, :, None, :] +def population_lda_decisions( + xs: npt.NDArray[np.float64], positive: npt.NDArray[np.float64], train: npt.NDArray[np.float64] +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Decision values and weights of linear discriminants with equal priors, one per problem, fold and feature set. + + Matches scikit-learn's `lsqr` solver with `shrinkage="auto"`. Each class's covariance is the Ledoit-Wolf + estimate of its standardised training trials, scaled back, and the pooled covariance is the mean of the two. + + Args: + xs (npt.NDArray[np.float64]): Standardised features from `_standardised`, shape `(n_problems, n_folds, + n_trials, n_sets, n_features)`. + positive (npt.NDArray[np.float64]): 1 for positive trials, shape `(n_problems, n_trials)`. + train (npt.NDArray[np.float64]): 1 for training trials, shape `(n_problems, n_folds, n_trials)`. + + Returns: + tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: Decision values of every trial under every fit, + shape `(n_problems, n_folds, n_trials, n_sets)`, positive favouring the positive class, and the weights, + shape `(n_problems, n_folds, n_sets, n_features)`. + """ + means, covariances = [], [] + for weight in (train * positive[:, None, :], train * (1 - positive)[:, None, :]): + n = weight.sum(axis=-1)[:, :, None, None] + mean = np.einsum("pfn,pfnsk->pfsk", weight, xs) / n + centred = (xs - mean[:, :, None]) * weight[..., None, None] + scale = _scale((centred**2).sum(axis=2) / n, mean, n) + shrunk = _ledoit_wolf(centred / scale[:, :, None], n) + covariances.append(scale[..., :, None] * shrunk * scale[..., None, :]) + means.append(mean) + coef = _solve(0.5 * (covariances[0] + covariances[1]), means[0] - means[1]) + decisions = np.einsum("pfnsk,pfsk->pfns", xs - 0.5 * (means[0] + means[1])[:, :, None], coef) + return decisions, coef + + +def _ledoit_wolf(z: npt.NDArray[np.float64], n: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + """Ledoit-Wolf covariance of each class's standardised trials, as `sklearn.covariance.ledoit_wolf` takes it. + + `z` holds the trials, zero outside the class, with shape `(n_problems, n_folds, n_trials, n_sets, + n_features)`, and `n` the class sizes. + + Returns: + npt.NDArray[np.float64]: Shape `(n_problems, n_folds, n_sets, n_features, n_features)`. + """ + n_features = z.shape[-1] + samples = np.moveaxis(z, 2, -1) + gram = samples @ np.swapaxes(samples, -1, -2) + n = n[..., 0] + covariance = gram / n[..., None, None] + variance = np.diagonal(covariance, axis1=-2, axis2=-1).sum(axis=-1) + mu = variance / n_features + gram_squares = (gram**2).sum(axis=(-2, -1)) / n**2 + beta = (((z**2).sum(axis=-1) ** 2).sum(axis=2) / n - gram_squares) / (n_features * n) + delta = (gram_squares - 2 * mu * variance + n_features * mu**2) / n_features + beta = np.minimum(beta, delta) + shrinkage = np.divide(beta, delta, out=np.zeros_like(beta), where=beta != 0) + return (1 - shrinkage)[..., None, None] * covariance + (shrinkage * mu)[..., None, None] * np.eye(n_features) + + +def _solve(a: npt.NDArray[np.float64], b: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + """Solutions of `a @ x = b` for stacked symmetric matrices, through the pseudo-inverse if any is singular. + + Returns: + npt.NDArray[np.float64]: Same shape as `b`. + """ + try: + return np.linalg.solve(a, b[..., None])[..., 0] + except np.linalg.LinAlgError: + return (np.linalg.pinv(a, rcond=np.finfo(np.float64).eps, hermitian=True) @ b[..., None])[..., 0] + + def region_decisions( features: npt.NDArray[np.float64], y: npt.NDArray, fold_of: npt.NDArray[np.int64], decoder: str ) -> npt.NDArray[np.float64]: @@ -573,6 +652,58 @@ def region_decisions( return decisions +def population_decisions( + features: npt.NDArray[np.float64], y: npt.NDArray, fold_of: npt.NDArray[np.int64], decoder: str +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Held-out decision values and weights of multi-feature decoders, for every problem and feature set at once. + + Problems and folds are as for `region_decisions`. Each feature set, such as the regions of one window, gets + its own decoder. + + Args: + features (npt.NDArray[np.float64]): Features of shape `(n_trials, n_sets, n_features)`. + y (npt.NDArray): Class per trial per problem, shape `(n_problems, n_trials)`. + fold_of (npt.NDArray[np.int64]): Test fold per trial per problem, shape `(n_problems, n_trials)`. + decoder (str): One of `DECODERS`. + + Returns: + tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: Decision values of shape `(n_problems, n_trials, + n_sets)`, and the weights of the standardised features in each fit, shape `(n_problems, n_folds, n_sets, + n_features)`. + + Raises: + ValueError: If `decoder` is not one of `DECODERS`. + NotImplementedError: For the logistic regression, which is not vectorised yet. + """ + if decoder not in DECODERS: + msg = f"Unknown decoder {decoder!r}; expected one of {', '.join(DECODERS)}." + raise ValueError(msg) + if decoder == "logistic": + msg = "The population logistic regression is not vectorised yet." + raise NotImplementedError(msg) + folds = int(fold_of.max()) + 1 + n_trials, n_sets, n_features = features.shape + # Elements per set of the standardised features or of the covariances, whichever is larger. + per_set = folds * n_features * max(n_trials, n_features) + set_chunk = max(1, min(n_sets, _CHUNK_ELEMENTS // per_set)) + chunk = max(1, _CHUNK_ELEMENTS // (per_set * set_chunk)) + decisions = np.empty((len(y), n_trials, n_sets)) + weights = np.empty((len(y), folds, n_sets, n_features)) + for first in range(0, len(y), chunk): + last = min(first + chunk, len(y)) + fold = fold_of[first:last] + positive = np.asarray(y[first:last], dtype=np.float64) + train = (fold[:, None, :] != np.arange(folds)[None, :, None]).astype(np.float64) + for start in range(0, n_sets, set_chunk): + sets = slice(start, min(start + set_chunk, n_sets)) + x = features[:, sets] + xs = _standardised(x.reshape(n_trials, -1), train).reshape(*train.shape, *x.shape[1:]) + fitted, coef = population_lda_decisions(xs, positive, train) + decisions[first:last, :, sets] = _test_fold(fitted, fold) + weights[first:last, :, sets] = coef + return decisions, weights + + def problem_scores( decisions: npt.NDArray[np.float64], y: npt.NDArray, fold_of: npt.NDArray[np.int64] ) -> dict[str, npt.NDArray[np.float64]]: diff --git a/tests/test_process.py b/tests/test_process.py index 3e162db..cc3cea9 100644 --- a/tests/test_process.py +++ b/tests/test_process.py @@ -5260,6 +5260,100 @@ def test_unknown_decoder_raises(self): pdc.region_decisions(np.ones((36, 1)), y[None], fold_of[:1], "svm") +class TestPopulationDecisions: + @staticmethod + def _problems(): + noise = np.random.default_rng(0).standard_normal((36, 3)) + y, _ = pdc.trial_labels(_gonogo_info(), "stim") + splits = pdc.cross_validation_splits(y, 5, 2, 42) + return np.column_stack([_decoding_features(), noise]), y, splits, pdc.fold_ids(splits, 5, len(y)) + + @staticmethod + def _assert_matches_sklearn(features, y, splits, decisions, weights): + for i, (train, test) in enumerate(splits): + model = pdc.make_decoder("lda").fit(features[train], y[train]) + np.testing.assert_allclose(decisions[i // 5, test, 0], model.decision_function(features[test]), atol=1e-10) + np.testing.assert_allclose(weights[i // 5, i % 5, 0], model[-1].coef_[0], atol=1e-10) + + def test_lda_matches_sklearn(self): + features, y, splits, fold_of = self._problems() + decisions, weights = pdc.population_decisions(features[:, None], np.broadcast_to(y, (2, 36)), fold_of, "lda") + assert decisions.shape == (2, 36, 1) + assert weights.shape == (2, 5, 1, 6) + self._assert_matches_sklearn(features, y, splits, decisions, weights) + + def test_lda_more_features_than_trials(self): + # Shrinkage keeps the covariance invertible with fewer trials per class than features. + _, y, splits, fold_of = self._problems() + features = 0.5 * y[:, None] + np.random.default_rng(1).standard_normal((36, 40)) + decisions, weights = pdc.population_decisions(features[:, None], y[None], fold_of[:1], "lda") + self._assert_matches_sklearn(features, y, splits[:5], decisions, weights) + + def test_lda_constant_feature(self): + features, y, splits, fold_of = self._problems() + features[:, 1] = 2.0 + decisions, weights = pdc.population_decisions(features[:, None], y[None], fold_of[:1], "lda") + self._assert_matches_sklearn(features, y, splits[:5], decisions, weights) + assert (weights[..., 1] == 0).all() + + def test_lda_one_feature_matches_region_decisions(self): + features, y, _, fold_of = self._problems() + labels = np.broadcast_to(y, (2, 36)) + decisions, _ = pdc.population_decisions(features[:, :, None], labels, fold_of, "lda") + np.testing.assert_allclose(decisions, pdc.region_decisions(features, labels, fold_of, "lda"), atol=1e-10) + + def test_sets_are_independent(self): + features, y, _, fold_of = self._problems() + sets = np.stack([features, features**2, -features], axis=1) + decisions, weights = pdc.population_decisions(sets, y[None], fold_of[:1], "lda") + for s in range(3): + alone, alone_weights = pdc.population_decisions(sets[:, [s]], y[None], fold_of[:1], "lda") + np.testing.assert_allclose(decisions[..., [s]], alone) + np.testing.assert_allclose(weights[:, :, [s]], alone_weights) + + def test_problems_are_independent(self): + features, y, _, fold_of = self._problems() + shuffled = np.random.default_rng(0).permutation(y) + together, _ = pdc.population_decisions(features[:, None], np.stack([y, shuffled]), fold_of, "lda") + alone, _ = pdc.population_decisions(features[:, None], shuffled[None], fold_of[1:], "lda") + np.testing.assert_allclose(together[1], alone[0]) + assert not np.allclose(together[0], together[1]) + + def test_chunks_match(self, monkeypatch): + features, y, _, fold_of = self._problems() + sets = np.stack([features, features**2, -features], axis=1) + labels = np.stack([y, np.random.default_rng(0).permutation(y)]) + whole = pdc.population_decisions(sets, labels, fold_of, "lda") + monkeypatch.setattr(pdc, "_CHUNK_ELEMENTS", 1) + chunked = pdc.population_decisions(sets, labels, fold_of, "lda") + np.testing.assert_allclose(chunked[0], whole[0]) + np.testing.assert_allclose(chunked[1], whole[1]) + + def test_constant_features_predict_nothing(self): + # Every covariance is zero, so the solve falls back to the pseudo-inverse. + _, y, _, fold_of = self._problems() + decisions, weights = pdc.population_decisions(np.ones((36, 1, 2)), y[None], fold_of[:1], "lda") + np.testing.assert_allclose(decisions, 0.0, atol=1e-12) + np.testing.assert_allclose(weights, 0.0, atol=1e-12) + + def test_separable_features(self): + _, y, _, fold_of = self._problems() + features = y[:, None] + 0.01 * np.random.default_rng(0).standard_normal((36, 3)) + decisions, _ = pdc.population_decisions(features[:, None], y[None], fold_of[:1], "lda") + assert np.isfinite(decisions).all() + assert ((decisions[0, :, 0] > 0) == (y == 1)).all() + + def test_logistic_not_vectorised_yet(self): + _, y, _, fold_of = self._problems() + with pytest.raises(NotImplementedError, match="not vectorised yet"): + pdc.population_decisions(np.ones((36, 1, 2)), y[None], fold_of[:1], "logistic") + + def test_unknown_decoder_raises(self): + _, y, _, fold_of = self._problems() + with pytest.raises(ValueError, match="Unknown decoder 'svm'"): + pdc.population_decisions(np.ones((36, 1, 2)), y[None], fold_of[:1], "svm") + + class TestProblemScores: def test_matches_sklearn_metrics(self): from sklearn.metrics import balanced_accuracy_score From 9bf6ba26f91ad69e840fd8682276263576efed63 Mon Sep 17 00:00:00 2001 From: Constantinos Eleftheriou Date: Mon, 12 Oct 2026 15:57:08 +0100 Subject: [PATCH 2/5] vectorise population and rolling logistic decoder --- src/mesoscopy/process/decoding.py | 171 ++++++++++++++++-------------- tests/test_process.py | 116 +++++++++++++------- 2 files changed, 168 insertions(+), 119 deletions(-) diff --git a/src/mesoscopy/process/decoding.py b/src/mesoscopy/process/decoding.py index 9560af1..f4a2912 100644 --- a/src/mesoscopy/process/decoding.py +++ b/src/mesoscopy/process/decoding.py @@ -472,79 +472,6 @@ def lda_decisions( return xs * coef[:, :, None, :] + intercept[:, :, None, :] -def logistic_decisions( - xs: npt.NDArray[np.float64], - positive: npt.NDArray[np.float64], - train: npt.NDArray[np.float64], - c: float = 1.0, - tolerance: float = 1e-8, - max_iter: int = 100, -) -> npt.NDArray[np.float64]: - """Decision values of one-feature balanced L2 logistic regressions, one per problem, fold and feature. - - Newton's method with backtracking on `0.5 * w**2 / c + sum(s_i * log(1 + exp(-y_i * (w * x_i + b))))`, the - objective scikit-learn's `lbfgs` solver minimises, with `s_i` the balanced class weights. - - Args: - xs (npt.NDArray[np.float64]): Standardised features from `_standardised`, shape `(n_problems, n_folds, - n_trials, n_features)`. - positive (npt.NDArray[np.float64]): 1 for positive trials, shape `(n_problems, n_trials)`. - train (npt.NDArray[np.float64]): 1 for training trials, shape `(n_problems, n_folds, n_trials)`. - c (float, optional): Inverse L2 strength. Defaults to 1. - tolerance (float, optional): Largest parameter step at convergence. Defaults to 1e-8. - max_iter (int, optional): Newton iterations. Defaults to 100. - - Returns: - npt.NDArray[np.float64]: Decision values of every trial under every fit, same shape as `xs`. - """ - from scipy.special import expit - - n_train = train.sum(axis=-1, keepdims=True) - n_positive = (train * positive[:, None, :]).sum(axis=-1, keepdims=True) - class_weight = np.where( - positive[:, None, :] > 0, n_train / (2 * n_positive), n_train / (2 * (n_train - n_positive)) - ) - weight = (train * class_weight)[..., None] - sign = (2 * positive - 1)[:, None, :, None] - target = positive[:, None, :, None] - shape = (*xs.shape[:2], xs.shape[3]) - w = np.zeros(shape) - b = np.zeros(shape) - - def objective(w: npt.NDArray[np.float64], b: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: - z = xs * w[:, :, None, :] + b[:, :, None, :] - return 0.5 * w**2 / c + (weight * np.logaddexp(0, -sign * z)).sum(axis=2) - - value = objective(w, b) - for _ in range(max_iter): - z = xs * w[:, :, None, :] + b[:, :, None, :] - p = expit(z) - residual = weight * (p - target) - curvature = weight * p * (1 - p) - g_w = w / c + (residual * xs).sum(axis=2) - g_b = residual.sum(axis=2) - h_ww = 1 / c + (curvature * xs**2).sum(axis=2) - h_wb = (curvature * xs).sum(axis=2) - h_bb = curvature.sum(axis=2) - det = np.maximum(h_ww * h_bb - h_wb**2, np.finfo(float).tiny) - step_w = (h_bb * g_w - h_wb * g_b) / det - step_b = (h_ww * g_b - h_wb * g_w) / det - scale = np.ones(shape) - for _ in range(20): - candidate = objective(w - scale * step_w, b - scale * step_b) - # Allow for rounding in the objective, or the search stalls at the optimum. - worse = candidate > value * (1 + 1e-12) + 1e-12 - if not worse.any(): - break - scale[worse] *= 0.5 - w -= scale * step_w - b -= scale * step_b - value = objective(w, b) if worse.any() else candidate - if max(np.abs(scale * step_w).max(), np.abs(scale * step_b).max()) < tolerance: - break - return xs * w[:, :, None, :] + b[:, :, None, :] - - def population_lda_decisions( xs: npt.NDArray[np.float64], positive: npt.NDArray[np.float64], train: npt.NDArray[np.float64] ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: @@ -614,13 +541,97 @@ def _solve(a: npt.NDArray[np.float64], b: npt.NDArray[np.float64]) -> npt.NDArra return (np.linalg.pinv(a, rcond=np.finfo(np.float64).eps, hermitian=True) @ b[..., None])[..., 0] +def population_logistic_decisions( + xs: npt.NDArray[np.float64], + positive: npt.NDArray[np.float64], + train: npt.NDArray[np.float64], + c: float = 1.0, + tolerance: float = 1e-8, + max_iter: int = 100, +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Decision values and weights of balanced L2 logistic regressions, one per problem, fold and feature set. + + Newton's method with backtracking on `0.5 * |w|**2 / c + sum(s_i * log(1 + exp(-y_i * (w @ x_i + b))))`, the + objective scikit-learn's `lbfgs` solver minimises, with `s_i` the balanced class weights. + + Args: + xs (npt.NDArray[np.float64]): Standardised features from `_standardised`, shape `(n_problems, n_folds, + n_trials, n_sets, n_features)`. + positive (npt.NDArray[np.float64]): 1 for positive trials, shape `(n_problems, n_trials)`. + train (npt.NDArray[np.float64]): 1 for training trials, shape `(n_problems, n_folds, n_trials)`. + c (float, optional): Inverse L2 strength. Defaults to 1. + tolerance (float, optional): Largest parameter step at convergence. Defaults to 1e-8. + max_iter (int, optional): Newton iterations. Defaults to 100. + + Returns: + tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: Decision values of every trial under every fit, + shape `(n_problems, n_folds, n_trials, n_sets)`, and the weights, shape `(n_problems, n_folds, n_sets, + n_features)`. + """ + from scipy.special import expit + + x = np.ascontiguousarray(np.moveaxis(xs, 2, -2)) + n_train = train.sum(axis=-1, keepdims=True) + n_positive = (train * positive[:, None, :]).sum(axis=-1, keepdims=True) + class_weight = np.where( + positive[:, None, :] > 0, n_train / (2 * n_positive), n_train / (2 * (n_train - n_positive)) + ) + weight = (train * class_weight)[:, :, None, :] + sign = (2 * positive - 1)[:, None, None, :] + target = positive[:, None, None, :] + n_features = x.shape[-1] + shape = x.shape[:3] + w = np.zeros((*shape, n_features)) + b = np.zeros(shape) + hessian = np.zeros((*shape, n_features + 1, n_features + 1)) + penalty = np.eye(n_features) / c + + def decision(w: npt.NDArray[np.float64], b: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + return (x @ w[..., None])[..., 0] + b[..., None] + + def objective(w: npt.NDArray[np.float64], b: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: + return 0.5 * (w**2).sum(axis=-1) / c + (weight * np.logaddexp(0, -sign * decision(w, b))).sum(axis=-1) + + value = objective(w, b) + for _ in range(max_iter): + p = expit(decision(w, b)) + residual = weight * (p - target) + curvature = weight * p * (1 - p) + weighted = curvature[..., None] * x + border = weighted.sum(axis=-2) + hessian[..., :-1, :-1] = np.swapaxes(x, -1, -2) @ weighted + penalty + hessian[..., :-1, -1] = border + hessian[..., -1, :-1] = border + hessian[..., -1, -1] = curvature.sum(axis=-1) + gradient = np.concatenate( + [w / c + (residual[..., None, :] @ x)[..., 0, :], residual.sum(axis=-1)[..., None]], axis=-1 + ) + step = _solve(hessian, gradient) + step_w, step_b = step[..., :-1], step[..., -1] + scale = np.ones(shape) + for _ in range(20): + candidate = objective(w - scale[..., None] * step_w, b - scale * step_b) + # Allow for rounding in the objective, or the search stalls at the optimum. + worse = candidate > value * (1 + 1e-12) + 1e-12 + if not worse.any(): + break + scale[worse] *= 0.5 + w -= scale[..., None] * step_w + b -= scale * step_b + value = objective(w, b) if worse.any() else candidate + if np.abs(scale[..., None] * step).max() < tolerance: + break + return np.moveaxis(decision(w, b), -1, 2), w + + def region_decisions( features: npt.NDArray[np.float64], y: npt.NDArray, fold_of: npt.NDArray[np.int64], decoder: str ) -> npt.NDArray[np.float64]: """Held-out decision values of one-feature decoders, for every problem and feature at once. A problem is one labelling of the trials with one k-fold assignment; the fits of a fold train on the other - folds and score the held-out trials. + folds and score the held-out trials. The logistic regressions go through `population_decisions` with one + feature per set. Args: features (npt.NDArray[np.float64]): Features of shape `(n_trials, n_features)`, one decoder per column. @@ -637,6 +648,8 @@ def region_decisions( if decoder not in DECODERS: msg = f"Unknown decoder {decoder!r}; expected one of {', '.join(DECODERS)}." raise ValueError(msg) + if decoder == "logistic": + return population_decisions(features[:, :, None], y, fold_of, decoder)[0] folds = int(fold_of.max()) + 1 n_trials, n_features = features.shape chunk = max(1, _CHUNK_ELEMENTS // (folds * n_trials * n_features)) @@ -647,8 +660,7 @@ def region_decisions( positive = np.asarray(y[first:last], dtype=np.float64) train = (fold[:, None, :] != np.arange(folds)[None, :, None]).astype(np.float64) xs = _standardised(features, train) - fitted = lda_decisions(xs, positive, train) if decoder == "lda" else logistic_decisions(xs, positive, train) - decisions[first:last] = _test_fold(fitted, fold) + decisions[first:last] = _test_fold(lda_decisions(xs, positive, train), fold) return decisions @@ -673,14 +685,11 @@ def population_decisions( Raises: ValueError: If `decoder` is not one of `DECODERS`. - NotImplementedError: For the logistic regression, which is not vectorised yet. """ if decoder not in DECODERS: msg = f"Unknown decoder {decoder!r}; expected one of {', '.join(DECODERS)}." raise ValueError(msg) - if decoder == "logistic": - msg = "The population logistic regression is not vectorised yet." - raise NotImplementedError(msg) + fit = population_lda_decisions if decoder == "lda" else population_logistic_decisions folds = int(fold_of.max()) + 1 n_trials, n_sets, n_features = features.shape # Elements per set of the standardised features or of the covariances, whichever is larger. @@ -698,7 +707,7 @@ def population_decisions( sets = slice(start, min(start + set_chunk, n_sets)) x = features[:, sets] xs = _standardised(x.reshape(n_trials, -1), train).reshape(*train.shape, *x.shape[1:]) - fitted, coef = population_lda_decisions(xs, positive, train) + fitted, coef = fit(xs, positive, train) decisions[first:last, :, sets] = _test_fold(fitted, fold) weights[first:last, :, sets] = coef return decisions, weights diff --git a/tests/test_process.py b/tests/test_process.py index cc3cea9..1a12aa6 100644 --- a/tests/test_process.py +++ b/tests/test_process.py @@ -5269,85 +5269,125 @@ def _problems(): return np.column_stack([_decoding_features(), noise]), y, splits, pdc.fold_ids(splits, 5, len(y)) @staticmethod - def _assert_matches_sklearn(features, y, splits, decisions, weights): + def _assert_matches_sklearn(decoder, features, y, splits, decisions, weights): + from sklearn.linear_model import LogisticRegression + from sklearn.pipeline import make_pipeline + from sklearn.preprocessing import StandardScaler + for i, (train, test) in enumerate(splits): - model = pdc.make_decoder("lda").fit(features[train], y[train]) - np.testing.assert_allclose(decisions[i // 5, test, 0], model.decision_function(features[test]), atol=1e-10) - np.testing.assert_allclose(weights[i // 5, i % 5, 0], model[-1].coef_[0], atol=1e-10) + if decoder == "lda": + model, atol = pdc.make_decoder("lda"), 1e-10 + else: + # lbfgs stops about 1e-6 short of the optimum even at this tolerance; see test_logistic_reaches_optimum. + logistic = LogisticRegression(C=1.0, class_weight="balanced", tol=1e-12, max_iter=100000) + model, atol = make_pipeline(StandardScaler(), logistic), 1e-5 + model.fit(features[train], y[train]) + np.testing.assert_allclose(decisions[i // 5, test, 0], model.decision_function(features[test]), atol=atol) + np.testing.assert_allclose(weights[i // 5, i % 5, 0], model[-1].coef_[0], atol=atol) - def test_lda_matches_sklearn(self): + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_matches_sklearn(self, decoder): features, y, splits, fold_of = self._problems() - decisions, weights = pdc.population_decisions(features[:, None], np.broadcast_to(y, (2, 36)), fold_of, "lda") + decisions, weights = pdc.population_decisions(features[:, None], np.broadcast_to(y, (2, 36)), fold_of, decoder) assert decisions.shape == (2, 36, 1) assert weights.shape == (2, 5, 1, 6) - self._assert_matches_sklearn(features, y, splits, decisions, weights) + self._assert_matches_sklearn(decoder, features, y, splits, decisions, weights) - def test_lda_more_features_than_trials(self): - # Shrinkage keeps the covariance invertible with fewer trials per class than features. + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_more_features_than_trials(self, decoder): + # Shrinkage or the L2 penalty keeps each fit well posed with fewer trials per class than features. _, y, splits, fold_of = self._problems() features = 0.5 * y[:, None] + np.random.default_rng(1).standard_normal((36, 40)) - decisions, weights = pdc.population_decisions(features[:, None], y[None], fold_of[:1], "lda") - self._assert_matches_sklearn(features, y, splits[:5], decisions, weights) + decisions, weights = pdc.population_decisions(features[:, None], y[None], fold_of[:1], decoder) + self._assert_matches_sklearn(decoder, features, y, splits[:5], decisions, weights) - def test_lda_constant_feature(self): + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_constant_feature(self, decoder): features, y, splits, fold_of = self._problems() features[:, 1] = 2.0 - decisions, weights = pdc.population_decisions(features[:, None], y[None], fold_of[:1], "lda") - self._assert_matches_sklearn(features, y, splits[:5], decisions, weights) + decisions, weights = pdc.population_decisions(features[:, None], y[None], fold_of[:1], decoder) + self._assert_matches_sklearn(decoder, features, y, splits[:5], decisions, weights) assert (weights[..., 1] == 0).all() + @staticmethod + def _assert_logistic_optimum(fold_of, y, atol=1e-9): + # The gradient of the objective vanishes in every fit. + features = 0.5 * y[:, None] + np.random.default_rng(1).standard_normal((36, 40)) + positive = y[None].astype(np.float64) + train = (fold_of[:1, None, :] != np.arange(5)[None, :, None]).astype(np.float64) + xs = pdc._standardised(features, train) + decisions, weights = pdc.population_logistic_decisions(xs[:, :, :, None], positive, train) + n_train = train.sum(axis=-1, keepdims=True) + n_positive = (train * positive[:, None]).sum(axis=-1, keepdims=True) + balanced = np.where(positive[:, None] > 0, n_train / (2 * n_positive), n_train / (2 * (n_train - n_positive))) + residual = train * balanced * (sst.logistic.cdf(decisions[..., 0]) - positive[:, None]) + np.testing.assert_allclose(weights[:, :, 0] + np.einsum("pfn,pfnk->pfk", residual, xs), 0.0, atol=atol) + np.testing.assert_allclose(residual.sum(axis=-1), 0.0, atol=atol) + + def test_logistic_reaches_optimum(self): + _, y, _, fold_of = self._problems() + self._assert_logistic_optimum(fold_of, y) + + def test_logistic_backtracks_overshooting_steps(self, monkeypatch): + # Newton steps three times too long overshoot, and the line search halves them back. Convergence is then + # linear, so the step tolerance leaves a larger gradient. + _, y, _, fold_of = self._problems() + solve = pdc._solve + monkeypatch.setattr(pdc, "_solve", lambda a, b: 3 * solve(a, b)) + self._assert_logistic_optimum(fold_of, y, atol=1e-5) + def test_lda_one_feature_matches_region_decisions(self): features, y, _, fold_of = self._problems() labels = np.broadcast_to(y, (2, 36)) decisions, _ = pdc.population_decisions(features[:, :, None], labels, fold_of, "lda") np.testing.assert_allclose(decisions, pdc.region_decisions(features, labels, fold_of, "lda"), atol=1e-10) - def test_sets_are_independent(self): + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_sets_are_independent(self, decoder): features, y, _, fold_of = self._problems() sets = np.stack([features, features**2, -features], axis=1) - decisions, weights = pdc.population_decisions(sets, y[None], fold_of[:1], "lda") + decisions, weights = pdc.population_decisions(sets, y[None], fold_of[:1], decoder) for s in range(3): - alone, alone_weights = pdc.population_decisions(sets[:, [s]], y[None], fold_of[:1], "lda") - np.testing.assert_allclose(decisions[..., [s]], alone) - np.testing.assert_allclose(weights[:, :, [s]], alone_weights) + alone, alone_weights = pdc.population_decisions(sets[:, [s]], y[None], fold_of[:1], decoder) + np.testing.assert_allclose(decisions[..., [s]], alone, atol=1e-12) + np.testing.assert_allclose(weights[:, :, [s]], alone_weights, atol=1e-12) - def test_problems_are_independent(self): + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_problems_are_independent(self, decoder): features, y, _, fold_of = self._problems() shuffled = np.random.default_rng(0).permutation(y) - together, _ = pdc.population_decisions(features[:, None], np.stack([y, shuffled]), fold_of, "lda") - alone, _ = pdc.population_decisions(features[:, None], shuffled[None], fold_of[1:], "lda") - np.testing.assert_allclose(together[1], alone[0]) + together, _ = pdc.population_decisions(features[:, None], np.stack([y, shuffled]), fold_of, decoder) + alone, _ = pdc.population_decisions(features[:, None], shuffled[None], fold_of[1:], decoder) + np.testing.assert_allclose(together[1], alone[0], atol=1e-12) assert not np.allclose(together[0], together[1]) - def test_chunks_match(self, monkeypatch): + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_chunks_match(self, monkeypatch, decoder): features, y, _, fold_of = self._problems() sets = np.stack([features, features**2, -features], axis=1) labels = np.stack([y, np.random.default_rng(0).permutation(y)]) - whole = pdc.population_decisions(sets, labels, fold_of, "lda") + whole = pdc.population_decisions(sets, labels, fold_of, decoder) monkeypatch.setattr(pdc, "_CHUNK_ELEMENTS", 1) - chunked = pdc.population_decisions(sets, labels, fold_of, "lda") - np.testing.assert_allclose(chunked[0], whole[0]) - np.testing.assert_allclose(chunked[1], whole[1]) + chunked = pdc.population_decisions(sets, labels, fold_of, decoder) + np.testing.assert_allclose(chunked[0], whole[0], atol=1e-12) + np.testing.assert_allclose(chunked[1], whole[1], atol=1e-12) - def test_constant_features_predict_nothing(self): - # Every covariance is zero, so the solve falls back to the pseudo-inverse. + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_constant_features_predict_nothing(self, decoder): + # For the discriminant every covariance is zero, so the solve falls back to the pseudo-inverse. _, y, _, fold_of = self._problems() - decisions, weights = pdc.population_decisions(np.ones((36, 1, 2)), y[None], fold_of[:1], "lda") + decisions, weights = pdc.population_decisions(np.ones((36, 1, 2)), y[None], fold_of[:1], decoder) np.testing.assert_allclose(decisions, 0.0, atol=1e-12) np.testing.assert_allclose(weights, 0.0, atol=1e-12) - def test_separable_features(self): + @pytest.mark.parametrize("decoder", pdc.DECODERS) + def test_separable_features(self, decoder): _, y, _, fold_of = self._problems() features = y[:, None] + 0.01 * np.random.default_rng(0).standard_normal((36, 3)) - decisions, _ = pdc.population_decisions(features[:, None], y[None], fold_of[:1], "lda") + decisions, _ = pdc.population_decisions(features[:, None], y[None], fold_of[:1], decoder) assert np.isfinite(decisions).all() assert ((decisions[0, :, 0] > 0) == (y == 1)).all() - def test_logistic_not_vectorised_yet(self): - _, y, _, fold_of = self._problems() - with pytest.raises(NotImplementedError, match="not vectorised yet"): - pdc.population_decisions(np.ones((36, 1, 2)), y[None], fold_of[:1], "logistic") - def test_unknown_decoder_raises(self): _, y, _, fold_of = self._problems() with pytest.raises(ValueError, match="Unknown decoder 'svm'"): From 1020855e92e1fa3d962294bdd09bcbc675714a68 Mon Sep 17 00:00:00 2001 From: Constantinos Eleftheriou Date: Mon, 12 Oct 2026 16:48:27 +0100 Subject: [PATCH 3/5] update decoding docs --- docs/how-to/decoding.md | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/how-to/decoding.md b/docs/how-to/decoding.md index 4ca7b94..fe19460 100644 --- a/docs/how-to/decoding.md +++ b/docs/how-to/decoding.md @@ -100,7 +100,7 @@ takes in activity after it. The chance level is taken once per decoder, region and group on the `--baseline` window (default every sample before the event), since shuffled labels give the same chance level in every window. The shuffle columns therefore repeat along a time course and `p_value` is read against them in every window. The rolling table is the slowest -part of the command, since the population decoder is fitted in every window. +part of the command, since every decoder is fitted in every window. ## Reading the table @@ -128,5 +128,7 @@ tables.table.query("decoder == 'lda' and group == 'all'").nlargest(5, "balanced_ `decoding.decoding_table` returns a `DecodingTables` with the scores in `table` and the population weights in `weights`. `decoding.rolling_table` returns the rolling table. `decoding.region_decisions` fits the one-feature -decoders of every region, fold and shuffle at once from a `(n_trials, n_features)` array and the fold of each trial, -and `decoding.fit_scores` fits one scikit-learn decoder per fold. +decoders of every region, fold and shuffle at once from a `(n_trials, n_features)` array and the fold of each trial. +`decoding.population_decisions` does the same for the population decoders from a `(n_trials, n_sets, n_features)` +array, with one decoder per set of features such as the regions of one window. `decoding.fit_scores` fits one +scikit-learn decoder and acts as a reference to compare against. From dd6e9fae5aef33765b462ce590b572ed1213c794 Mon Sep 17 00:00:00 2001 From: Constantinos Eleftheriou Date: Mon, 12 Oct 2026 16:49:58 +0100 Subject: [PATCH 4/5] use numpy vectorised decoders instead of scikit --- src/mesoscopy/process/decoding.py | 78 +++++++++++++++++-------------- tests/test_process.py | 55 ++++++++++++++++++++++ 2 files changed, 97 insertions(+), 36 deletions(-) diff --git a/src/mesoscopy/process/decoding.py b/src/mesoscopy/process/decoding.py index f4a2912..01194f0 100644 --- a/src/mesoscopy/process/decoding.py +++ b/src/mesoscopy/process/decoding.py @@ -23,9 +23,9 @@ Each region, and all regions together, is used to classify the trials of a go/no-go session with a logistic regression and a linear discriminant decoder under repeated stratified cross-validation, against a label-shuffle -null, over the cue-to-response epoch or in windows sliding over the peri-event window. The one-feature decoders are -fitted in numpy over every region, window, fold and shuffle at once; the population decoder goes through -scikit-learn. +null, over the cue-to-response epoch or in windows sliding over the peri-event window. The decoders are fitted in +numpy over every region, window, fold and shuffle at once, and `fit_scores` keeps the scikit-learn fits they are +tested against. """ from __future__ import annotations @@ -223,7 +223,7 @@ def fit_scores( decoder: str, folds: int, ) -> dict[str, typing.Any]: - """Cross-validated scores and weights of one decoder. + """Cross-validated scores and weights of one scikit-learn decoder, the reference for `population_decisions`. Each repeat scores every trial once, so the confusion counts are summed over the folds of a repeat and averaged over repeats, and the AUROC is taken over each repeat's pooled decision values. @@ -835,7 +835,7 @@ def group_decoding( All fits share the same splits, and all shuffles the same permutations, so the regions are compared on the same test trials. The one-feature decoders go through `region_decisions`, the population decoder through - `fit_scores`. + `population_decisions`. Args: features (npt.NDArray[np.float64]): Features of the group's trials, shape `(n_trials, n_regions)`. @@ -864,13 +864,13 @@ def group_decoding( result: dict[str, dict[str, typing.Any]] = {} for decoder in DECODERS: scores: dict[str, typing.Any] = {} - observed = problem_scores( - region_decisions(x, observed_labels, observed_folds, decoder), observed_labels, observed_folds - ) + decisions, weights = _all_decisions(x, observed_labels, observed_folds, decoder) + observed = problem_scores(decisions, observed_labels, observed_folds) null = None if shuffles: - null = problem_scores(region_decisions(x, null_labels, null_folds, decoder), null_labels, null_folds) - for i, column in enumerate(columns): + null_decisions, _ = _all_decisions(x, null_labels, null_folds, decoder) + null = problem_scores(null_decisions, null_labels, null_folds) + for i, name in enumerate([*(regions[column] for column in columns), POPULATION]): summary = _column_summary(observed, i) if null is not None: summary |= null_summary( @@ -878,29 +878,31 @@ def group_decoding( null["accuracy"][:, :, i].mean(axis=1), null["balanced_accuracy"][:, :, i].mean(axis=1), ) - scores[regions[column]] = summary - - population = fit_scores(x, y, splits, decoder, folds) - summary = score_summary(population) - if shuffles: - null_scores = [ - fit_scores(x, row, s, decoder, folds) for row, s in zip(null_labels, null_splits, strict=True) - ] - summary |= null_summary( - summary["balanced_accuracy"], - np.array([n["accuracy"].mean() for n in null_scores]), - np.array([n["balanced_accuracy"].mean() for n in null_scores]), - ) - scores[POPULATION] = summary + scores[name] = summary + weights = weights.reshape(len(splits), len(columns)) mean = np.full(len(regions), np.nan) sd = np.full(len(regions), np.nan) - mean[columns] = population["weights"].mean(axis=0) - sd[columns] = population["weights"].std(axis=0, ddof=1) if len(splits) > 1 else np.nan + mean[columns] = weights.mean(axis=0) + sd[columns] = weights.std(axis=0, ddof=1) if len(splits) > 1 else np.nan scores["weights"] = (mean, sd) result[decoder] = scores return result +def _all_decisions( + x: npt.NDArray[np.float64], y: npt.NDArray, fold_of: npt.NDArray[np.int64], decoder: str +) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Held-out decision values of each region's decoder and, last, of the population decoder. + + Returns: + tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: Decision values of shape `(n_problems, n_trials, + n_regions + 1)` from `region_decisions` and `population_decisions`, and the population weights of each + fit, shape `(n_problems, n_folds, n_regions)`. + """ + population, weights = population_decisions(x[:, None], y, fold_of, decoder) + return np.concatenate([region_decisions(x, y, fold_of, decoder), population], axis=-1), weights[:, :, 0] + + def labelled_trials( perievent: pd.DataFrame, trials: pd.DataFrame | None = None ) -> tuple[npt.NDArray[np.float64], pd.DataFrame]: @@ -1082,11 +1084,8 @@ def _null_scores( """ labels, splits = shuffled_labels(y, folds, shuffles, rng) fold_of = np.array([fold_ids(s, folds, len(y))[0] for s in splits], dtype=np.int64) - scores = problem_scores(region_decisions(baseline, labels, fold_of, decoder), labels, fold_of) - population = [fit_scores(baseline, row, s, decoder, folds) for row, s in zip(labels, splits, strict=True)] - accuracy = np.column_stack([scores["accuracy"].mean(axis=1), [p["accuracy"].mean() for p in population]]) - balanced = scores["balanced_accuracy"].mean(axis=1) - return accuracy, np.column_stack([balanced, [p["balanced_accuracy"].mean() for p in population]]) + scores = problem_scores(_all_decisions(baseline, labels, fold_of, decoder)[0], labels, fold_of) + return scores["accuracy"].mean(axis=1), scores["balanced_accuracy"].mean(axis=1) def rolling_decoding( @@ -1101,9 +1100,9 @@ def rolling_decoding( ) -> dict[str, list[list[dict[str, float]]]]: """Scores of every decoder on each region and on all regions in every window, for one trial group. - The one-feature decoders of every window and region are fitted at once through `region_decisions` and the - population decoder is fitted per window through `fit_scores`. The null is taken once on `baseline` and applied - to every window. + The one-feature decoders of every window and region are fitted at once through `region_decisions`, and the + population decoders of every window through `population_decisions`. The null is taken once on `baseline` and + applied to every window. Args: features (npt.NDArray[np.float64]): Window features of the group's trials, shape `(n_trials, n_windows, @@ -1130,7 +1129,14 @@ def rolling_decoding( x = features[:, :, columns].reshape(n_trials, n_windows * len(columns)) result: dict[str, list[list[dict[str, float]]]] = {} for decoder in DECODERS: - observed = problem_scores(region_decisions(x, labels, fold_of, decoder), labels, fold_of) + decisions = np.concatenate( + [ + region_decisions(x, labels, fold_of, decoder), + population_decisions(features[:, :, columns], labels, fold_of, decoder)[0], + ], + axis=-1, + ) + observed = problem_scores(decisions, labels, fold_of) null = _null_scores(baseline[:, columns], y, folds, shuffles, rng, decoder) if shuffles else None per_region: list[list[dict[str, float]]] = [[{} for _ in range(n_windows)] for _ in range(n_regions + 1)] for i, column in enumerate(columns): @@ -1140,7 +1146,7 @@ def rolling_decoding( summary |= null_summary(summary["balanced_accuracy"], null[0][:, i], null[1][:, i]) per_region[column][w] = summary for w in range(n_windows): - summary = score_summary(fit_scores(features[:, w, columns], y, splits, decoder, folds)) + summary = _column_summary(observed, n_windows * len(columns) + w) if null is not None: summary |= null_summary(summary["balanced_accuracy"], null[0][:, -1], null[1][:, -1]) per_region[n_regions][w] = summary diff --git a/tests/test_process.py b/tests/test_process.py index 1a12aa6..4316dd7 100644 --- a/tests/test_process.py +++ b/tests/test_process.py @@ -5737,6 +5737,37 @@ def test_null_seed(self): other = pdc.decoding_table(_decoding_perievent(), repeats=2, shuffles=10, seed=1).table assert not first["shuffle_balanced_accuracy"].equals(other["shuffle_balanced_accuracy"]) + def test_null_matches_fit_scores(self): + # The population null, refitted with scikit-learn on the same shuffles. + perievent = _decoding_perievent() + table = pdc.decoding_table(perievent, repeats=2, shuffles=10).table + features = _decoding_features(perievent) + y, _ = pdc.trial_labels(_gonogo_info(), "stim") + labels, splits = pdc.shuffled_labels(y, 5, 10, np.random.default_rng(42)) + for decoder in pdc.DECODERS: + null = [pdc.fit_scores(features, row, s, decoder, 5) for row, s in zip(labels, splits, strict=True)] + row = self._row(table, decoder, pdc.POPULATION, "all") + expected = pdc.null_summary( + row["balanced_accuracy"], + np.array([n["accuracy"].mean() for n in null]), + np.array([n["balanced_accuracy"].mean() for n in null]), + ) + for score in pdc.NULL_SCORES: + assert row[score] == pytest.approx(expected[score]) + + def test_weights_match_fit_scores(self): + perievent = _decoding_perievent() + weights = pdc.decoding_table(perievent, repeats=2, shuffles=0).weights + features = _decoding_features(perievent) + y, _ = pdc.trial_labels(_gonogo_info(), "stim") + splits = pdc.cross_validation_splits(y, 5, 2, 42) + # scikit-learn's lbfgs stops short of the logistic optimum at its default tolerance. + for decoder, atol in (("lda", 1e-10), ("logistic", 1e-3)): + fitted = pdc.fit_scores(features, y, splits, decoder, 5)["weights"] + rows = weights[(weights["decoder"] == decoder) & (weights["group"] == "all")] + np.testing.assert_allclose(rows["weight"], fitted.mean(axis=0), atol=atol) + np.testing.assert_allclose(rows["weight_sd"], fitted.std(axis=0, ddof=1), atol=atol) + class TestWindowFeatures: def test_means(self): @@ -5859,6 +5890,30 @@ def test_null(self): assert fitted["p_value"].between(1 / 21, 1).all() assert table[table["group"] == "resp-nopush"][list(pdc.NULL_SCORES)].isna().all().all() + def test_null_matches_fit_scores(self): + # The population null, refitted with scikit-learn on the same shuffles of the pre-event means. + perievent = _decoding_perievent() + table = pdc.rolling_table(perievent, repeats=2, shuffles=10) + info = _gonogo_info() + unbounded = np.full(len(info), np.inf) + cube = pc.epoch_cube( + perievent, info["trial_index"].to_numpy(), METRICS_GRID, CONNECTIVITY_REGIONS, -unbounded, unbounded + ) + baseline = np.nanmean(cube[:, METRICS_GRID < 0], axis=1) + y, _ = pdc.trial_labels(info, "stim") + rng = np.random.default_rng(42) + for decoder in pdc.DECODERS: + labels, splits = pdc.shuffled_labels(y, 5, 10, rng) + null = [pdc.fit_scores(baseline, row, s, decoder, 5) for row, s in zip(labels, splits, strict=True)] + row = self._rows(table, decoder, pdc.POPULATION, "all").iloc[0] + expected = pdc.null_summary( + row["balanced_accuracy"], + np.array([n["accuracy"].mean() for n in null]), + np.array([n["balanced_accuracy"].mean() for n in null]), + ) + for score in pdc.NULL_SCORES: + assert row[score] == pytest.approx(expected[score]) + def test_baseline(self): perievent = _decoding_perievent() default = pdc.rolling_table(perievent, repeats=2, shuffles=10) From 7b483f3584425aaeb9f4a1690cc7e6a82241bec1 Mon Sep 17 00:00:00 2001 From: Constantinos Eleftheriou Date: Mon, 12 Oct 2026 16:50:16 +0100 Subject: [PATCH 5/5] update changelog --- CHANGELOG.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6ae08be..1883773 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,10 @@ Versions follow [Semantic Versioning](https://semver.org) (`..