From 6938115f9423be246355677a79e3f7446fe437ad Mon Sep 17 00:00:00 2001 From: Kenya Otsuka Date: Mon, 7 Sep 2026 09:48:49 +0000 Subject: [PATCH 1/2] test: Add pytest infrastructure and regression tests for the sklearn Ridge decoder The repository had no automated tests, which made it unsafe to change the decoding implementation. This adds the first test suite, pinning the current behavior of `train_decoder_sklearn_ridge.py` and `predict_feature.py`. - Add `scikit-learn` to the project dependencies. It is imported by `train_decoder_sklearn_ridge.py` but was never declared, so `uv sync` produced an environment that could not run the script. - Add a `dev` dependency group with pytest and `[tool.pytest.ini_options]`, including `pythonpath = ["."]` so the top-level analysis modules are importable from the tests without any `sys.path` manipulation. - `tests/helpers/synthetic.py`: build tiny `BData` / `Features` datasets on disk, with repeated stimuli (unbalanced repeats), two ROIs, a subject split over two `.h5` files, and both a 2-D `fc`-like and a 3-D `conv`-like feature layer so the chunking / `order='F'` path is exercised. - `tests/helpers/legacy_ridge.py`: a plain numpy + scikit-learn reference implementation of the direct brain -> feature Ridge decoder, mirroring bdpy's `ModelTraining` / `ModelTest` semantics exactly. - `tests/generate_golden.py` and `tests/data/golden/`: regression fixtures recording the pipeline output on a fixed-seed synthetic dataset. Run it as `python -m tests.generate_golden`. - Tests covering the artifact contract, decoded-feature output format, normalization parameters, label alignment, chunking, and the behavior of feature units with zero training variance. - Add a GitHub Actions workflow running the suite on push and pull requests. No change to the decoding implementation. --- .github/workflows/test.yml | 27 +++ .gitignore | 3 + README.md | 16 ++ pyproject.toml | 21 ++ tests/__init__.py | 0 tests/conftest.py | 45 ++++ tests/data/golden/sklearn_ridge_pipeline.npz | Bin 0 -> 12230 bytes tests/generate_golden.py | 91 ++++++++ tests/helpers/__init__.py | 0 tests/helpers/legacy_ridge.py | 164 ++++++++++++++ tests/helpers/pipeline.py | 99 +++++++++ tests/helpers/synthetic.py | 221 +++++++++++++++++++ tests/test_chunking.py | 155 +++++++++++++ tests/test_label_alignment.py | 80 +++++++ tests/test_legacy_reference.py | 97 ++++++++ tests/test_pipeline_golden.py | 61 +++++ tests/test_predict_feature.py | 135 +++++++++++ tests/test_train_decoder_sklearn_ridge.py | 185 ++++++++++++++++ tests/test_zero_variance_features.py | 47 ++++ uv.lock | 129 +++++++++++ 20 files changed, 1576 insertions(+) create mode 100644 .github/workflows/test.yml create mode 100644 tests/__init__.py create mode 100644 tests/conftest.py create mode 100644 tests/data/golden/sklearn_ridge_pipeline.npz create mode 100644 tests/generate_golden.py create mode 100644 tests/helpers/__init__.py create mode 100644 tests/helpers/legacy_ridge.py create mode 100644 tests/helpers/pipeline.py create mode 100644 tests/helpers/synthetic.py create mode 100644 tests/test_chunking.py create mode 100644 tests/test_label_alignment.py create mode 100644 tests/test_legacy_reference.py create mode 100644 tests/test_pipeline_golden.py create mode 100644 tests/test_predict_feature.py create mode 100644 tests/test_train_decoder_sklearn_ridge.py create mode 100644 tests/test_zero_variance_features.py diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml new file mode 100644 index 0000000..0c43d7e --- /dev/null +++ b/.github/workflows/test.yml @@ -0,0 +1,27 @@ +name: tests + +on: + push: + branches: [main, dev] + pull_request: + +jobs: + pytest: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: Install uv + uses: astral-sh/setup-uv@v5 + with: + enable-cache: true + + # bdpy declares Requires-Python <3.12. + - name: Set up Python + run: uv python install 3.11 + + - name: Install dependencies + run: uv sync --group dev --python 3.11 + + - name: Run tests + run: uv run --python 3.11 pytest diff --git a/.gitignore b/.gitignore index 2d4689a..6c8452f 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,6 @@ .venv .ipynb_checkpoints +__pycache__/ +*.py[cod] +.pytest_cache/ diff --git a/README.md b/README.md index d50e4b7..df7d46b 100644 --- a/README.md +++ b/README.md @@ -85,6 +85,22 @@ $ python cv_predict_feature_fastl2lir.py config/deeprecon_cv_pyfastl2lir_alpha10 $ python cv_evaluation.py config/deeprecon_cv_pyfastl2lir_alpha100_vgg19_allunits.yaml ``` +### Tests + +The test suite runs on small synthetic data generated on the fly; no downloaded +dataset is required. + +```shell +# Install the test dependencies and run the suite +$ uv sync --group dev +$ uv run pytest +``` + +`tests/data/golden/` holds regression fixtures recording the numerical output of +the decoding pipeline. Regenerate them with +`uv run python -m tests.generate_golden` only when the expected output is meant +to change, and say so explicitly in the commit message. + ## References - Horikawa and Kamitani (2017) Generic decoding of seen and imagined objects using hierarchical visual features. *Nature Communications* 8:15037. https://www.nature.com/articles/ncomms15037 diff --git a/pyproject.toml b/pyproject.toml index b879cb1..68d0a59 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -12,8 +12,29 @@ dependencies = [ "matplotlib", "numpy<2.0.0", "pyyaml", + "scikit-learn", "tqdm", ] +[dependency-groups] +dev = [ + "pytest>=7.0", +] + [tool.uv] package = false + +[tool.pytest.ini_options] +testpaths = ["tests"] +addopts = "-ra" +# The analysis scripts are top-level modules rather than an installed package, +# so the repository root has to be importable. `pythonpath` states that +# explicitly instead of relying on pytest's implicit rootdir insertion. +pythonpath = ["."] +markers = [ + "slow: tests that take more than a few seconds", +] +filterwarnings = [ + "ignore::DeprecationWarning", + "ignore::PendingDeprecationWarning", +] diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..ae95565 --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,45 @@ +"""Shared fixtures for the feature-decoding tests.""" + +from __future__ import annotations + +import os + +import pytest + +from tests.helpers import synthetic + +GOLDEN_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'data', 'golden') + + +@pytest.fixture(autouse=True) +def _run_in_tmp_path(tmp_path, monkeypatch): + """Run every test in its own directory. + + Both scripts create ``./tmp/.db`` (the DistComp sqlite lock + database) relative to the current working directory. + """ + monkeypatch.chdir(tmp_path) + + +@pytest.fixture +def dataset(tmp_path): + """A small synthetic dataset with repeated stimuli.""" + return synthetic.make_dataset(str(tmp_path / 'data'), seed=0) + + +@pytest.fixture +def dataset_with_constant_unit(tmp_path): + """Synthetic dataset where one feature unit has zero training variance.""" + return synthetic.make_dataset(str(tmp_path / 'data'), seed=0, + constant_feature_unit=True) + + +@pytest.fixture +def decoder_dir(tmp_path): + return str(tmp_path / 'feature_decoders') + + +@pytest.fixture +def decoded_dir(tmp_path): + return str(tmp_path / 'decoded_features') diff --git a/tests/data/golden/sklearn_ridge_pipeline.npz b/tests/data/golden/sklearn_ridge_pipeline.npz new file mode 100644 index 0000000000000000000000000000000000000000..fa374e9e57cd263509ef06bfd3fb73c821f8a05a GIT binary patch literal 12230 zcmds730PCd7ET0$qAXPe6_D7f)T*esy(R@zd|Ings<@!gG=Uh|BC-Vy6l$L!?i(&p zgsRp06cA)7YM^xi7lc|WvM52QTTK)n^`cD!Zw9$=@8n*BPoelO-*7*^kNNLE|D1E? z%*=6_tD`%X!C>?TZ5Tspz1yz2^)X$o zUN45ka}2T!Oqr>sEbv(u5`)9(N&m(f7()i0D?`M-JbsvWP*>5goux^GlRkL`PNNT` z5hN6Pp#=W|&i-JN8!nT1$4ZZml2NS+8~y+qqzvz78@?9ML=_&d;z#pAdnM&@;G=;mKBJ zZdw=~={!RLV7P$o+)LoWp283~u?t3ubuS>iq9ucia9%C*7PzY2aFScOEY`y5V5d3R z8nale*=&p{n94TN5r|m2=tN#AUbhJ+WCDc>76t`_*%S!<+O(=^<>Rh}msY}4 zpEoN!jOQ=pi*f7U-pX=GhirGP(5hlI$6#cO^?B)}2929!Y-t`|J3eInk6%8TJ6gdc zhlPV>vn8D&Dqa=n-=#!P=@AFFy6pK2tfs zr};YjTHbZVcK^BsW2MEL8k9~xPbwS($~}z{`$)gQ{9e24Yj{N2yJ|E$_xs&H`QZCZ zr^tu;DjTs>ao#zAJAA%urDH>blXSJklo*>``4tj-*VW^te-t=8zA;CNjK){l%)m!8 zl-Ak<0^1|K3G?NY$X|J8kQEQi%QP?DnJ>+qk$Ul2wvp6qQB}#l^y&lLnjG2eXxaI7 z(PrGz0oSly=Z;nec|53!tP9?#-f6s0A-bu^o?7+jAv<`)O^1=vuhhDADs0tzhE|i# zNUG-lelCsoQX_Y-m zoDi$!0hjJ};LPA4wt4u3jS1w@5X0)N< z_`gCqhx5Ho9p4vtv*~H|93y0!V+Z}>r=2&^WV?hn8Yx^N?r!>8{cEmK5U7Y8SM6Pr zv#(;h)75MI=*Rps$%Tp76|Pyi^+$_xn$Jy-P(>DU6#4)0$TYvLd|Y$ZzOX6gxmcC4 zWQ>u=z0Zwo)=sIqZJH?Hoyw@Hj+dTm`q1}17M>Kn4{$RRBlNAozh$kAv!IE6WpN=!gSId2Dk3KhznP_XdHTT`C z3aehfSs^^NWE8D}zKwx_-aOuxjAF;unP#{RQG~hHx#^DJ5Nl&)WfwZ4Jnc@EFn>Wz z#ho9r)6S=Enhn-`iq#pgMWI7<9$W_2>VM8g@MOqC3oeR z{D-qvn^WZF-j|)zo5R1o;jmd$)^yt0TIFVbztYjeVeNg=T2Ti%?k( zT5M4+ldUeRe(|TPKH`J*rtC z?R9z*_kK)*tOz^eZK23KollN$JY)E0@$q;Yi(jgpT~mnXTO{MjMHD@*-Vm`zLq=R>Tp>2_A*b=T16 z=e8?dnuAX#E-f1<-R0__ve&!X7-KpvBXZKy-Ri_fE*dn5hZN9?S9tdLBMN$&L<7aM zi8p$Zf^Iv|UB(AC+3C3L(a|I)6J0{{=%;pyknq8lVY$xbJipj;0L2g578Obr(iqk|WL8D9y{^=jw5q(37cbljzBKObUD5WdZCB+78d+N9@21BqOWL z*Zb!*Rz1JH^^xlib%uPlVo~0KQ%feg3HC-5${k!KdG;kqw=s`%A|lI$Fhex7d=fXG2dvSpv}F;Ck&s9yO%}`pxD`V z#nNUwa@0a;g3b*ro_DkI2`}LMV>j7pbwRnUw79?!8?t4#0xP<7^2u5YVg9_ZEfIP- z8|)DWBknjUM$s)x+w1K&=~xpD1Yg?z*4-x2{isRzFz2HqNmTZ4S=lDY*ZShFr?IkW z><+ua@lv@q7ZiIB(PLu184@I;ZIb{rCNJw3EB^qfZ+X%@lN0!#ast|BLyL7XIQua z>oVX+`+$Z&QOrKzvu$M@Xg64I0gEXvfz}uA1Kv^0SlF$-`U33;3n5^Q1DpgyQpZuf z;P(h+5olLfs{ohBWYEY|e}Y}{$_2C=tSEq$0yUtSfQCO&%mq528&Y3@_J*5jxQsx% zX5(GKJBnEkdTYD?N73N5JT=;WPzgbYhWi(a29K-Y9XQmI4I2JLF=yz}p|yt;6b;^8 zQ=^%%=+Wr7!MkI43k+|BhvLM*JBnFNi-xD4C^NX-f}4K}Olit#kU5(F8{HXms4* zJx2}TqG|o(L4>B}MmGmABcrz1FF-ic!0OQY;ki7_d6GdRTiKv@=& z2%}eu@?@+*V6dbE-!?3|@zGm)9t~FKR44NiT2b?&w{*DI$8{&qU_Ss&HAXt+iP4>m z^ou4^hn7x^?PR14)*77tM@zS7Qz3Oub26t|8YDVPfLlCVak8by45Q-IuD2e@HyMm^ O;717Zg1PU)zy1qbM_KLw literal 0 HcmV?d00001 diff --git a/tests/generate_golden.py b/tests/generate_golden.py new file mode 100644 index 0000000..92cecb9 --- /dev/null +++ b/tests/generate_golden.py @@ -0,0 +1,91 @@ +"""Regenerate the golden regression fixtures for the sklearn Ridge decoder. + +The fixtures record the normalization parameters and decoded features produced +by ``train_decoder_sklearn_ridge.py`` + ``predict_feature.py`` on a fixed-seed +synthetic dataset. They were generated from the *direct* (pre-factorization) +implementation and are what the factorized implementation must reproduce. + +Only run this when the expected numerical output is intended to change, and say +so explicitly in the commit message. + +Note that regeneration is not bit-for-bit stable across environments: the +pipeline runs in ``float32`` and the Ridge solve's summation order depends on +how BLAS partitions the work, which moves the last bits (observed ~5e-8 +relative). That is why the tests compare against these fixtures with a +tolerance rather than exactly, and why regenerating them without an intended +behavior change only adds noise. + + uv run python -m tests.generate_golden +""" + +from __future__ import annotations + +import argparse +import os +import tempfile + +import numpy as np + +from tests.helpers import pipeline, synthetic + +ALPHA = 100 +CHUNK_AXIS = 1 +SEED = 0 +GOLDEN_FILE = os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'data', 'golden', 'sklearn_ridge_pipeline.npz') + + +def build(work_dir: str) -> dict: + """Run the full pipeline in ``work_dir`` and collect everything of interest.""" + cwd = os.getcwd() + os.chdir(work_dir) + try: + dataset = synthetic.make_dataset(os.path.join(work_dir, 'data'), + seed=SEED) + decoder_dir = os.path.join(work_dir, 'feature_decoders') + decoded_dir = os.path.join(work_dir, 'decoded_features') + + pipeline.run_training(dataset, decoder_dir, alpha=ALPHA, + chunk_axis=CHUNK_AXIS) + pipeline.run_prediction(dataset, decoder_dir, decoded_dir, + chunk_axis=CHUNK_AXIS) + + golden = { + '_alpha': np.array(ALPHA), + '_chunk_axis': np.array(CHUNK_AXIS), + '_seed': np.array(SEED), + '_layers': np.array(dataset.layers), + '_subjects': np.array(sorted(dataset.train_fmri)), + '_rois': np.array(sorted(dataset.rois)), + '_test_labels': np.array(dataset.unique_test_labels), + } + for layer in dataset.layers: + for subject in sorted(dataset.train_fmri): + for roi in sorted(dataset.rois): + prefix = '%s|%s|%s|' % (layer, subject, roi) + for key, value in pipeline.read_norm_params( + decoder_dir, layer, subject, roi).items(): + golden[prefix + key] = value + golden[prefix + 'pred'] = pipeline.read_decoded_features( + decoded_dir, layer, subject, roi, + dataset.unique_test_labels) + return golden + finally: + os.chdir(cwd) + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('-o', '--output', default=GOLDEN_FILE) + args = parser.parse_args() + + with tempfile.TemporaryDirectory() as work_dir: + golden = build(work_dir) + + os.makedirs(os.path.dirname(args.output), exist_ok=True) + np.savez_compressed(args.output, **golden) + print('Saved %s (%d arrays)' % (args.output, len(golden))) + + +if __name__ == '__main__': + main() diff --git a/tests/helpers/__init__.py b/tests/helpers/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/helpers/legacy_ridge.py b/tests/helpers/legacy_ridge.py new file mode 100644 index 0000000..3ffe4c7 --- /dev/null +++ b/tests/helpers/legacy_ridge.py @@ -0,0 +1,164 @@ +"""Reference implementation of the *direct* sklearn Ridge feature decoder. + +This mirrors, in plain numpy + scikit-learn, exactly what +``train_decoder_sklearn_ridge.py`` + ``predict_feature.py`` did before the +decoder was factorized, including every detail bdpy's +``ModelTraining``/``ModelTest`` contribute: + +* brain z-scoring with ``ddof=1`` followed by ``X[np.isinf(X)] = 0``; +* feature z-scoring computed over the **unique** stimuli; +* one Ridge per index along ``chunk_axis`` when the feature array is at least + 3-D (``Y.ndim >= chunk_ndim + 1`` with bdpy's default ``chunk_ndim=2``); +* the chunk taken with ``np.take(..., [i], axis=chunk_axis)`` so the axis is + kept with size 1, and the normalization parameters sliced the same way; +* ``Y = Y[feat_index]`` applied *after* normalization -- this is the repetition + expansion that turns M unique feature rows into N trial-aligned rows; +* ``reshape(n, -1, order='F')`` before fitting and the inverse reshape after + predicting; +* ``float32`` casts for both X and Y; +* ``np.concatenate(..., axis=chunk_axis)`` to reassemble the chunks; +* un-normalization as ``y_pred * y_norm + y_mean``. + +It exists so that the factorized implementation can be compared against the old +math directly, and so that the old math itself is pinned by tests against the +real scripts (see ``tests/test_legacy_reference.py``). +""" + +from __future__ import annotations + +from typing import Any, Dict, List, Optional, Sequence +import copy + +import numpy as np +from sklearn.linear_model import Ridge + +CHUNK_NDIM = 2 # bdpy ModelTraining default + + +def brain_norm_params(brain: np.ndarray): + """``(mean, std)`` of the brain data, shaped ``(1, n_voxels)``.""" + mean = np.mean(brain, axis=0)[np.newaxis, :] + std = np.std(brain, axis=0, ddof=1)[np.newaxis, :] + return mean, std + + +def feature_norm_params(feat: np.ndarray): + """``(mean, std)`` of the features, shaped ``(1, *feature_shape)``.""" + mean = np.mean(feat, axis=0)[np.newaxis, :] + std = np.std(feat, axis=0, ddof=1)[np.newaxis, :] + return mean, std + + +def assignment_index(brain_labels: Sequence[str], + feat_labels: Sequence[str]) -> np.ndarray: + """Row of ``feat`` matching each brain trial (the ``Y_sort`` index).""" + return np.array([np.where(np.array(feat_labels) == bl) + for bl in brain_labels]).flatten() + + +def normalize_brain_for_training(brain: np.ndarray, mean: np.ndarray, + std: np.ndarray) -> np.ndarray: + """z-score as ``ModelTraining`` does it, including the ``inf`` clean-up.""" + with np.errstate(invalid='ignore', divide='ignore'): + x = (brain - mean) / std + x[np.isinf(x)] = 0 + return x + + +def use_chunking(feat: np.ndarray, chunk_axis: Optional[int]) -> bool: + if chunk_axis is None: + return False + return feat.ndim >= CHUNK_NDIM + 1 + + +def legacy_train(brain: np.ndarray, brain_labels: Sequence[str], + feat: np.ndarray, feat_labels: Sequence[str], + alpha: float = 100, chunk_axis: Optional[int] = 1, + dtype: Any = np.float32) -> Dict[str, Any]: + """Fit the direct brain -> feature Ridge decoder, chunk by chunk.""" + x_mean, x_norm = brain_norm_params(brain) + y_mean, y_norm = feature_norm_params(feat) + feat_index = assignment_index(brain_labels, feat_labels) + + x = normalize_brain_for_training(brain, x_mean, x_norm) + + chunking = use_chunking(feat, chunk_axis) + chunk_index: Sequence[Optional[int]] + chunk_index = range(feat.shape[chunk_axis]) if chunking else [None] + + models: List[Dict[str, Any]] = [] + for i_chunk in chunk_index: + if chunking: + y = np.take(feat, [i_chunk], axis=chunk_axis) + chunk_mean = np.take(y_mean, [i_chunk], axis=chunk_axis) + chunk_norm = np.take(y_norm, [i_chunk], axis=chunk_axis) + else: + y = feat + chunk_mean, chunk_norm = y_mean, y_norm + + with np.errstate(invalid='ignore', divide='ignore'): + y = (y - chunk_mean) / chunk_norm + y[np.isinf(y)] = 0 + y = y[feat_index] + + y_shape = y.shape[1:] + if y.ndim > 2: + y = y.reshape(y.shape[0], -1, order='F') + + model = Ridge(alpha=alpha) + model.fit(x.astype(dtype), y.astype(dtype)) + models.append({'model': model, 'y_shape': y_shape}) + + return { + 'models': models, + 'x_mean': x_mean, 'x_norm': x_norm, + 'y_mean': y_mean, 'y_norm': y_norm, + 'feat_index': feat_index, + 'chunk_axis': chunk_axis, + 'chunking': chunking, + } + + +def legacy_predict(trained: Dict[str, Any], brain_test: np.ndarray, + dtype: Any = np.float32) -> np.ndarray: + """Predict features from *raw* test brain data, as ``predict_feature`` did. + + Note the deliberate asymmetry with training: the prediction script does not + apply the ``inf`` clean-up after z-scoring. + """ + x = (brain_test - trained['x_mean']) / trained['x_norm'] + x = x.astype(dtype) + + preds = [] + for entry in trained['models']: + y_pred = entry['model'].predict(x) + if y_pred.shape[1:] != entry['y_shape']: + y_pred = y_pred.reshape((y_pred.shape[0],) + entry['y_shape'], + order='F') + preds.append(y_pred) + + if trained['chunk_axis'] is None: + pred = preds[0] + else: + pred = np.concatenate(preds, axis=trained['chunk_axis']) + + return pred * trained['y_norm'] + trained['y_mean'] + + +def average_test_brain(brain: np.ndarray, labels: Sequence[str], + excluded_labels: Sequence[str] = ()): + """Trial-average the test brain data, as ``predict_feature`` does.""" + unique = [lb for lb in np.unique(labels) if lb not in excluded_labels] + averaged = np.vstack([ + np.mean(brain[(np.array(labels) == lb).flatten(), :], axis=0) + for lb in unique]) + return averaged, list(unique) + + +def single_trial_labels(labels: Sequence[str]) -> List[str]: + """Output labels used by ``predict_feature`` when ``average_sample`` is off.""" + return ['sample{:06}-{}'.format(i + 1, lb) for i, lb in enumerate(labels)] + + +def clone_trained(trained: Dict[str, Any]) -> Dict[str, Any]: + return copy.deepcopy(trained) diff --git a/tests/helpers/pipeline.py b/tests/helpers/pipeline.py new file mode 100644 index 0000000..1b8a5b3 --- /dev/null +++ b/tests/helpers/pipeline.py @@ -0,0 +1,99 @@ +"""Run the real training / prediction scripts on a synthetic dataset. + +Shared by the tests and by ``tests/generate_golden.py`` so that the golden +fixtures and the tests exercise exactly the same code path. +""" + +from __future__ import annotations + +from typing import Any, Dict, List, Optional, Sequence +import inspect +import os + +import numpy as np +from bdpy.dataform import load_array + +import predict_feature +import train_decoder_sklearn_ridge + +NORM_KEYS = ('x_mean', 'x_norm', 'y_mean', 'y_norm') + + +def run_training(dataset, output_dir: str, alpha: float = 100, + chunk_axis: Optional[int] = 1, + layers: Optional[Sequence[str]] = None, + subjects: Optional[Sequence[str]] = None, + rois: Optional[Dict[str, str]] = None, + analysis_name: str = 'test_training') -> str: + """Train decoders with ``train_decoder_sklearn_ridge.py``.""" + fmri = {sbj: paths for sbj, paths in dataset.train_fmri.items() + if subjects is None or sbj in subjects} + return train_decoder_sklearn_ridge.featdec_sklearn_ridge_train( + fmri, + [dataset.train_features_dir], + output_dir=output_dir, + rois=dict(dataset.rois) if rois is None else dict(rois), + label_key=dataset.label_key, + layers=list(dataset.layers if layers is None else layers), + alpha=alpha, + chunk_axis=chunk_axis, + analysis_name=analysis_name, + ) + + +def run_prediction(dataset, decoder_dir: str, output_dir: str, + chunk_axis: Optional[int] = 1, + layers: Optional[Sequence[str]] = None, + subjects: Optional[Sequence[str]] = None, + rois: Optional[Dict[str, str]] = None, + average_sample: bool = True, + excluded_labels: Sequence[str] = (), + analysis_name: str = 'test_prediction') -> str: + """Predict features with ``predict_feature.py``. + + ``training_features_paths`` is passed only when the script accepts it, so + this helper works both before and after the decoder is factorized. + """ + fmri = {sbj: paths for sbj, paths in dataset.test_fmri.items() + if subjects is None or sbj in subjects} + kwargs: Dict[str, Any] = dict( + output_dir=output_dir, + rois=dict(dataset.rois) if rois is None else dict(rois), + label_key=dataset.label_key, + layers=list(dataset.layers if layers is None else layers), + excluded_labels=list(excluded_labels), + average_sample=average_sample, + chunk_axis=chunk_axis, + analysis_name=analysis_name, + ) + if 'training_features_paths' in inspect.signature( + predict_feature.featdec_predict).parameters: + kwargs['training_features_paths'] = [dataset.train_features_dir] + return predict_feature.featdec_predict(fmri, decoder_dir, **kwargs) + + +def model_dir(decoder_dir: str, layer: str, subject: str, roi: str) -> str: + return os.path.join(decoder_dir, layer, subject, roi, 'model') + + +def read_norm_params(decoder_dir: str, layer: str, subject: str, + roi: str) -> Dict[str, np.ndarray]: + directory = model_dir(decoder_dir, layer, subject, roi) + return {key: load_array(os.path.join(directory, key + '.mat'), key=key) + for key in NORM_KEYS} + + +def read_decoded_features(output_dir: str, layer: str, subject: str, roi: str, + labels: Sequence[str]) -> np.ndarray: + """Stack the saved per-label ``.mat`` files into one array.""" + directory = os.path.join(output_dir, layer, subject, roi) + return np.vstack([ + load_array(os.path.join(directory, '%s.mat' % label), key='feat') + for label in labels]) + + +def decoded_feature_labels(output_dir: str, layer: str, subject: str, + roi: str) -> List[str]: + directory = os.path.join(output_dir, layer, subject, roi) + return sorted(os.path.splitext(f)[0] for f in os.listdir(directory) + if f.endswith('.mat')) diff --git a/tests/helpers/synthetic.py b/tests/helpers/synthetic.py new file mode 100644 index 0000000..60fb3a9 --- /dev/null +++ b/tests/helpers/synthetic.py @@ -0,0 +1,221 @@ +"""Build tiny synthetic on-disk datasets for the feature-decoding tests. + +Everything is small enough to live in a temporary directory: a handful of +``bdpy.BData`` HDF5 files and a ``bdpy.dataform.Features`` tree. No part of the +real DeepRecon dataset is needed. + +The dataset intentionally reproduces the properties the decoders depend on: + +* stimuli are **repeated** across fMRI trials, with an *unbalanced* number of + repetitions per stimulus (3, 3, 2, 2, 1, 1); +* one subject's training trials are split over **two** ``.h5`` files, so + ``select_data_multi_bdatas`` / ``get_labels_multi_bdatas`` are exercised; +* two ROIs of different sizes; +* a 1-D ``fc``-like feature layer (no chunking) *and* a 3-D ``conv``-like layer, + so the ``chunk_axis`` / ``order='F'`` reshape path is exercised. The shipped + config only enables ``fc6``/``fc7``/``fc8``, which never reaches that path. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import Dict, List, Sequence, Tuple +import os + +import numpy as np +from bdpy import BData +from bdpy.dataform import save_array + +# Feature layers: name -> feature shape of a single stimulus. +LAYER_SHAPES: Dict[str, Tuple[int, ...]] = { + 'conv_like': (3, 2, 2), + 'fc_like': (4,), +} + +N_VOXELS = 8 +ROIS: Dict[str, str] = {'VC': 'ROI_VC = 1', 'HALF': 'ROI_HALF = 1'} +ROI_SIZES: Dict[str, int] = {'VC': N_VOXELS, 'HALF': N_VOXELS // 2} +LABEL_KEY = 'stimulus_name' + +# Unbalanced repetitions: 6 unique training stimuli, 12 trials. +TRAIN_REPEATS: Tuple[int, ...] = (3, 3, 2, 2, 1, 1) +# 3 unique test stimuli, 5 trials. +TEST_REPEATS: Tuple[int, ...] = (2, 2, 1) + + +@dataclass +class SyntheticDataset: + """Paths and in-memory copies of a synthetic decoding dataset.""" + + root: str + train_fmri: Dict[str, List[str]] + test_fmri: Dict[str, List[str]] + train_features_dir: str + train_brain: Dict[str, np.ndarray] = field(default_factory=dict) + train_labels: Dict[str, List[str]] = field(default_factory=dict) + test_brain: Dict[str, np.ndarray] = field(default_factory=dict) + test_labels: Dict[str, List[str]] = field(default_factory=dict) + features: Dict[str, np.ndarray] = field(default_factory=dict) + unique_train_labels: List[str] = field(default_factory=list) + unique_test_labels: List[str] = field(default_factory=list) + rois: Dict[str, str] = field(default_factory=lambda: dict(ROIS)) + label_key: str = LABEL_KEY + layer_shapes: Dict[str, Tuple[int, ...]] = field( + default_factory=lambda: dict(LAYER_SHAPES)) + + @property + def layers(self) -> List[str]: + return sorted(self.layer_shapes) + + def brain(self, subject: str, roi: str, split: str = 'train') -> np.ndarray: + """Voxel data for one subject/ROI, as the scripts would select it.""" + data = self.train_brain if split == 'train' else self.test_brain + return data[subject][:, :ROI_SIZES[roi]] + + def labels(self, subject: str, split: str = 'train') -> List[str]: + return (self.train_labels if split == 'train' else self.test_labels)[subject] + + def feature_matrix(self, layer: str, labels: Sequence[str]) -> np.ndarray: + """Rows of ``features[layer]`` in the order of ``labels``.""" + index = [self.unique_train_labels.index(lb) for lb in labels] + return self.features[layer][index] + + +def _expand(labels: Sequence[str], repeats: Sequence[int]) -> List[str]: + """Repeat each label ``repeats[i]`` times, interleaved as in a real run.""" + trials: List[str] = [] + remaining = list(repeats) + while any(r > 0 for r in remaining): + for i, label in enumerate(labels): + if remaining[i] > 0: + trials.append(label) + remaining[i] -= 1 + return trials + + +def write_bdata(path: str, brain: np.ndarray, labels: Sequence[str], + label_to_number: Dict[str, int]) -> None: + """Write one ``BData`` HDF5 file with ROI metadata and a label vmap.""" + n_voxels = brain.shape[1] + bdata = BData() + bdata.add(np.asarray(brain, dtype=float), 'VoxelData') + bdata.add(np.array([[label_to_number[lb]] for lb in labels], dtype=float), + LABEL_KEY) + bdata.add_metadata('ROI_VC', np.ones(n_voxels), where='VoxelData') + half = np.zeros(n_voxels) + half[:ROI_SIZES['HALF']] = 1 + bdata.add_metadata('ROI_HALF', half, where='VoxelData') + bdata.add_vmap(LABEL_KEY, {v: k for k, v in label_to_number.items()}) + os.makedirs(os.path.dirname(path), exist_ok=True) + bdata.save(path) + + +def write_features(root: str, layer_shapes: Dict[str, Tuple[int, ...]], + labels: Sequence[str], + arrays: Dict[str, np.ndarray]) -> None: + """Write a ``Features`` tree ``//