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..2d4a75f 100644
--- a/README.md
+++ b/README.md
@@ -58,15 +58,24 @@ $ python evaluation.py config/deeprecon_pyfastl2lir_alpha100_vgg19_allunits.yaml
 
 ```shell
 # Training of decoding models
-$ python train_decoder_sklearn_ridge.py config/deeprecon_pyfastl2lir_alpha100_vgg19_allunits.yaml
+$ python train_decoder_sklearn_ridge.py config/deeprecon_sklearn_ridge_alpha100_vgg19_allunits.yaml
 
 # Prediction of DNN features
-$ python predict_feature.py config/deeprecon_pyfastl2lir_alpha100_vgg19_allunits.yaml
+$ python predict_feature.py config/deeprecon_sklearn_ridge_alpha100_vgg19_allunits.yaml
 
 # Evaluation
-$ python evaluation.py config/deeprecon_pyfastl2lir_alpha100_vgg19_allunits.yaml
+$ python evaluation.py config/deeprecon_sklearn_ridge_alpha100_vgg19_allunits.yaml
 ```
 
+The scikit-learn Ridge decoder is trained and stored in a factorized form: the
+model maps brain activity onto the training stimulus basis, and prediction
+combines its coefficients with the training features. This is mathematically
+identical to regressing the features directly, but the stored decoder is much
+smaller and training does not scale with the feature dimension.
+`predict_feature.py` therefore reads the training features
+(`decoder.features.paths`, already set in the example config). See
+`ridge_factorization.py` for the details.
+
 ### Cross-validation feature decoding
 
 - Training: `cv_train_decoder_fastl2lir.py` (example for scikit-learn Ridge regression)
@@ -85,6 +94,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/config/README.md b/config/README.md
index 311ed2c..b9c2588 100644
--- a/config/README.md
+++ b/config/README.md
@@ -17,7 +17,7 @@ These settings are used by the decoder training script (e.g., `train_decoder_fas
   - `decoder.path`: Defines the path where the decoder data is located. The path can include placeholders like `${decoder.name}` and `${decoder.features.name}` to dynamically set paths based on the feature name.
   - `decoder.parameters`:
     - `alpha`: A regularization parameter set to 100.
-    - `chunk_axis`: Indicates that chunking is done along axis 1.
+    - `chunk_axis`: Indicates that chunking is done along axis 1. For the sklearn Ridge decoder it applies to the prediction step only.
 2. fMRI data
   - This specifies fMRI data used for the decoder training.
   - `decoder.fmri.name`: Refers to the dataset name, here "ImageNetTraining_fmriprep_volume_native".
@@ -27,7 +27,7 @@ These settings are used by the decoder training script (e.g., `train_decoder_fas
 3. Features
   - This specifies features used for the decoder training.
   - `decoder.features.name`: Refers to the pre-trained DNN model used.
-  - `decoder.features.paths`: Provides the path to the feature data.
+  - `decoder.features.paths`: Provides the path to the feature data. The factorized sklearn Ridge decoder reads these at prediction time as well, as the features of the training stimuli.
   - `decoder.features.layers`: Specifies the layers of the DNN from which features will be extracted.
 
 ## Decoded features
diff --git a/config/deeprecon_sklearn_ridge_alpha100_vgg19_allunits.yaml b/config/deeprecon_sklearn_ridge_alpha100_vgg19_allunits.yaml
index cf3babd..20f6390 100644
--- a/config/deeprecon_sklearn_ridge_alpha100_vgg19_allunits.yaml
+++ b/config/deeprecon_sklearn_ridge_alpha100_vgg19_allunits.yaml
@@ -5,6 +5,7 @@ decoder:
   path: ./data/feature_decoders/ImageNetTraining/${decoder.name}/${decoder.features.name}
   parameters:
     alpha: 100
+    # Axis along which the feature combination is blocked at prediction time.
     chunk_axis: 1
 
   fmri:
@@ -34,6 +35,7 @@ decoder:
 
   features:
     name: caffe/VGG19
+    # Features of the training stimuli; also read by `predict_feature.py`.
     paths:
       - ./data/features/ImageNetTraining/${decoder.features.name}
     layers:
diff --git a/predict_feature.py b/predict_feature.py
index 098a1b9..1a4fefd 100644
--- a/predict_feature.py
+++ b/predict_feature.py
@@ -1,12 +1,20 @@
-'''DNN Feature decoding - feature prediction script.'''
+'''DNN Feature decoding - feature prediction script.
 
+Decoders trained by ``train_decoder_sklearn_ridge.py`` are factorized (see
+``ridge_factorization``), so prediction runs in two steps: predict one
+coefficient per training stimulus, then combine those coefficients with the
+training features of the layer being decoded.  Decoders saved by the earlier,
+direct implementation are detected automatically and predicted as before.
+'''
 
-from typing import Dict, List, Optional
+
+from typing import Optional, Sequence
 
 from itertools import product
 import os
 import shutil
 from time import time
+import warnings
 
 import bdpy
 from bdpy.dataform import load_array, save_array
@@ -14,9 +22,147 @@
 from bdpy.ml import ModelTest
 from bdpy.pipeline.config import init_hydra_cfg
 from bdpy.util import makedir_ifnot
-from fastl2lir import FastL2LiR
 import numpy as np
 
+from ridge_factorization import (
+    FEATURE_INDEX_FILE,
+    TrainingFeatureLoader,
+    brain_model_dir,
+    check_column_correspondence,
+    combine_features,
+    feature_statistics,
+    is_factorized_model_dir,
+    load_factorized_metadata,
+    layer_model_dir,
+    load_train_labels,
+    model_file_path,
+    resolve_feature_index_path,
+    save_array_atomic,
+    verify_train_labels,
+)
+
+
+# Prediction back ends #######################################################
+
+def _training_features(model_dir: str, layer: str,
+                       feature_loader: TrainingFeatureLoader):
+    """Load the training features the decoder's coefficients refer to.
+
+    Returns ``(train_labels, features)``, one feature row per training stimulus
+    in the order the coefficient columns were fitted in.
+    """
+    metadata = load_factorized_metadata(model_dir)
+    train_labels = load_train_labels(model_dir)
+    verify_train_labels(metadata, train_labels, model_dir)
+    return train_labels, feature_loader.get(layer, train_labels)
+
+
+class FeatureStatistics(object):
+    """``y_mean``/``y_norm`` of the training features, cached per label set.
+
+    They depend on ``(layer, ordered training labels)`` only -- never on the
+    ROI -- and at DeepRecon scale each computation is two passes over a 15 GB
+    tensor, so they are computed once and reused for every decoder directory.
+    """
+
+    def __init__(self):
+        self._cache = {}
+        self.computations = 0  # observable in tests
+
+    def get(self, layer: str, train_labels: Sequence[str],
+            features: np.ndarray):
+        key = (layer, tuple(train_labels))
+        if key not in self._cache:
+            self._cache[key] = feature_statistics(features)
+            self.computations += 1
+        return self._cache[key]
+
+
+def _materialize_feature_statistics(decoder_path: str, layer: str,
+                                    subject: str, roi: str, model_dir: str,
+                                    feature_loader: TrainingFeatureLoader,
+                                    statistics: FeatureStatistics) -> None:
+    """Fill in the ``y_mean``/``y_norm`` that ``evaluation.py`` reads.
+
+    Compatibility sidecars, not parameters of the model: prediction never reads
+    them, and training cannot produce them because it never opens a feature
+    file, so the step that has the features writes them.  Only when missing --
+    one decoder directory belongs to one training feature configuration -- and
+    atomically, since parallel prediction workers can reach the same path.
+    """
+    directory = layer_model_dir(decoder_path, layer, subject, roi)
+    targets = {key: os.path.join(directory, '%s.mat' % key)
+               for key in ('y_mean', 'y_norm')}
+    if all(os.path.exists(path) for path in targets.values()):
+        return
+
+    train_labels, features = _training_features(model_dir, layer,
+                                                feature_loader)
+    values = dict(zip(('y_mean', 'y_norm'),
+                      statistics.get(layer, train_labels, features)))
+
+    for key, path in sorted(targets.items()):
+        if os.path.exists(path):
+            continue
+        try:
+            save_array_atomic(path, values[key], key=key, dtype=np.float32)
+            print('Saved %s' % path)
+        except Exception:
+            warnings.warn('Failed to save %s. The decoder directory may be '
+                          'read-only; evaluation.py will not find the training '
+                          'feature statistics.' % path)
+
+
+def _predict_factorized(model_dir: str, brain: np.ndarray,
+                        layer: str, chunk_axis: Optional[int],
+                        feature_loader: TrainingFeatureLoader) -> np.ndarray:
+    '''Predict features with a factorized decoder.
+
+    ``brain`` must already be normalized with the decoder's ``x_mean`` /
+    ``x_norm``.  The feature normalization parameters are not needed: they
+    cancel algebraically (``ridge_factorization``), so the prediction is the
+    linear combination of the raw training features.
+    '''
+    train_labels, features = _training_features(model_dir, layer,
+                                                feature_loader)
+
+    test = ModelTest(None, brain)
+    test.model_format = 'pickle'
+    test.model_path = model_file_path(model_dir)
+    test.dtype = np.float32
+    test.chunk_axis = None  # the target is the stimulus basis; never chunked
+
+    coefficients = test.run()
+    check_column_correspondence(np.asarray(coefficients).shape[1],
+                                train_labels, features, model_dir)
+
+    # np.asarray rather than .astype: the training features can be very large
+    # and are already float32, so this must not copy them.
+    return combine_features(np.asarray(coefficients, dtype=np.float32),
+                            np.asarray(features, dtype=np.float32),
+                            chunk_axis=chunk_axis)
+
+
+def _predict_legacy(model_dir: str, brain: np.ndarray,
+                    chunk_axis: Optional[int]) -> np.ndarray:
+    '''Predict features with a decoder saved by the direct implementation.
+
+    One Ridge model per chunk of the feature dimensions, and the prediction is
+    in normalized feature space, so it has to be un-normalized here.
+    '''
+    feat_mean = load_array(os.path.join(model_dir, 'y_mean.mat'), key='y_mean')  # shape = (1, shape_features)
+    feat_norm = load_array(os.path.join(model_dir, 'y_norm.mat'), key='y_norm')  # shape = (1, shape_features)
+
+    test = ModelTest(None, brain)
+    test.model_format = 'pickle'
+    test.model_path = model_dir
+    test.dtype = np.float32
+    test.chunk_axis = chunk_axis
+
+    feat_pred = test.run()
+
+    return feat_pred * feat_norm + feat_mean
+
 
 # Main #######################################################################
 
@@ -31,6 +177,7 @@ def featdec_predict(
         excluded_labels=[],
         average_sample=True,
         chunk_axis=1,
+        training_features_paths=None,
         analysis_name="feature_prediction"
 ):
     '''Feature prediction.
@@ -39,6 +186,8 @@ def featdec_predict(
 
     - fmri_data
     - feature_decoder_dir
+    - training_features_paths: directories holding the features of the training
+      stimuli.  Required for decoders in the factorized format.
 
     Output:
 
@@ -64,21 +213,43 @@ def featdec_predict(
     data_brain = {sbj: bdpy.BData(dat_file[0])
                   for sbj, dat_file in fmri_data.items()}
 
+    # Training features of the factorized decoders.  Only the layer currently
+    # being decoded is held in memory (see TrainingFeatureLoader).
+    feature_loader = TrainingFeatureLoader(
+        training_features_paths or [], feature_index_file=feature_index_file)
+    statistics = FeatureStatistics()
+
     # Initialize directories -------------------------------------------
     makedir_ifnot(output_dir)
     makedir_ifnot('tmp')
 
     # Save feature index -----------------------------------------------------
     if feature_index_file is not None:
-        feature_index_save_file = os.path.join(output_dir, 'feature_index.mat')
-        shutil.copy(feature_index_file, feature_index_save_file)
+        # The path as given wins; only when it does not resolve is the feature
+        # store consulted, which is where a factorized decoder's index lives
+        # (`TrainingFeatureLoader` hands the same resolved path to `Features`).
+        feature_index_source = feature_index_file
+        if not os.path.exists(feature_index_source) and training_features_paths:
+            feature_index_source = resolve_feature_index_path(
+                training_features_paths[0], feature_index_file)
+        feature_index_save_file = os.path.join(output_dir, FEATURE_INDEX_FILE)
+        shutil.copy(feature_index_source, feature_index_save_file)
         print('Saved %s' % feature_index_save_file)
 
     # Analysis loop ----------------------------------------------------
     print('----------------------------------------')
     print('Analysis loop')
 
+    # `layer` is the outer loop so that the training features of one layer are
+    # read once and reused for every subject and ROI.
+    current_layer = None
     for layer, sbj, roi in product(layers, fmri_data, rois):
+        if layer != current_layer:
+            # Release the previous layer's training features before loading the
+            # next ones: at most one layer is ever resident.
+            feature_loader.release()
+            current_layer = layer
+
         print('--------------------')
         print('Feature:    %s' % layer)
         print('Subject:    %s' % sbj)
@@ -89,10 +260,36 @@ def featdec_predict(
         analysis_id = analysis_name + '-' + sbj + '-' + roi + '-' + layer
         results_dir_prediction = os.path.join(output_dir, layer, sbj, roi)
 
+        # The brain-side directory is the only discriminator: once prediction
+        # has run once, a factorized decoder's per-layer directory looks like a
+        # legacy decoder directory from the outside.
+        shared_dir = brain_model_dir(decoder_path, sbj, roi)
+        factorized = is_factorized_model_dir(shared_dir)
+        model_dir = (shared_dir if factorized
+                     else layer_model_dir(decoder_path, layer, sbj, roi))
+
+        # Before the "already done" skip below: a run whose features were
+        # already decoded would otherwise skip forever and the statistics would
+        # never appear. Only for factorized decoders -- a legacy decoder wrote
+        # its own at training time -- and only when the features are available.
+        if factorized and training_features_paths:
+            _materialize_feature_statistics(
+                decoder_path, layer, sbj, roi, model_dir, feature_loader,
+                statistics)
+
         if os.path.exists(results_dir_prediction):
             print('%s is already done. Skipped.' % analysis_id)
             continue
 
+        # Before `makedir_ifnot`: a run that cannot proceed must not leave an
+        # output directory behind, or the next, correct run would skip it.
+        if factorized and not training_features_paths:
+            raise ValueError(
+                '%s is a factorized decoder, which needs the features of the '
+                'training stimuli. Pass training_features_paths (the '
+                'decoder.features.paths of the config used for training).'
+                % model_dir)
+
         makedir_ifnot(results_dir_prediction)
 
         distcomp_db = os.path.join('./tmp', analysis_name + '.db')
@@ -127,16 +324,10 @@ def featdec_predict(
 
         print('Elapsed time (data preparation): %f' % (time() - start_time))
 
-        # Model directory
-        # ---------------
-        model_dir = os.path.join(decoder_path, layer, sbj, roi, 'model')
-
         # Preprocessing
         # -------------
         brain_mean = load_array(os.path.join(model_dir, 'x_mean.mat'), key='x_mean')  # shape = (1, n_voxels)
         brain_norm = load_array(os.path.join(model_dir, 'x_norm.mat'), key='x_norm')  # shape = (1, n_voxels)
-        feat_mean = load_array(os.path.join(model_dir, 'y_mean.mat'), key='y_mean')  # shape = (1, shape_features)
-        feat_norm = load_array(os.path.join(model_dir, 'y_norm.mat'), key='y_norm')  # shape = (1, shape_features)
 
         brain = (brain - brain_mean) / brain_norm
 
@@ -146,20 +337,15 @@ def featdec_predict(
 
         start_time = time()
 
-        test = ModelTest(None, brain)
-        test.model_format = 'pickle'
-        test.model_path = model_dir
-        test.dtype = np.float32
-        test.chunk_axis = chunk_axis
-
-        feat_pred = test.run()
+        if factorized:
+            feat_pred = _predict_factorized(
+                model_dir, brain, layer, chunk_axis, feature_loader)
+        else:
+            print('Legacy decoder format detected in %s' % model_dir)
+            feat_pred = _predict_legacy(model_dir, brain, chunk_axis)
 
         print('Total elapsed time (prediction): %f' % (time() - start_time))
 
-        # Postprocessing
-        # --------------
-        feat_pred = feat_pred * feat_norm + feat_mean
-
         # Save results
         # ------------
         print('Saving results')
@@ -183,6 +369,8 @@ def featdec_predict(
 
         distcomp.unlock(analysis_id)
 
+    feature_loader.release()
+
     print('%s finished.' % analysis_name)
 
     return output_dir
@@ -216,6 +404,13 @@ def featdec_predict(
     average_sample = cfg["decoded_feature"]["parameters"]["average_sample"]
     excluded_labels = cfg.decoded_feature.fmri.get("exclude_labels", [])
 
+    # Features of the training stimuli. The factorized decoder references them
+    # instead of storing a copy, so prediction needs them too.
+    training_features_paths = cfg.decoded_feature.decoder.get(
+        "training_features_paths", None)
+    if training_features_paths is None:
+        training_features_paths = cfg["decoder"]["features"]["paths"]
+
     featdec_predict(
         test_fmri_data,
         decoder_path,
@@ -227,5 +422,6 @@ def featdec_predict(
         excluded_labels=excluded_labels,
         average_sample=average_sample,
         chunk_axis=cfg["decoder"]["parameters"]["chunk_axis"],
+        training_features_paths=training_features_paths,
         analysis_name=analysis_name
     )
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/ridge_factorization.py b/ridge_factorization.py
new file mode 100644
index 0000000..70ea42e
--- /dev/null
+++ b/ridge_factorization.py
@@ -0,0 +1,532 @@
+'''Factorized formulation of the sklearn Ridge feature decoder.
+
+Motivation
+----------
+
+The direct decoder regresses brain activity on DNN features and stores a
+coefficient matrix of shape ``(d_out, d_in)``.  For DeepRecon-scale features
+that matrix is enormous (``d_out`` reaches ~10^7 per VGG-19 layer against
+``d_in`` ~10^4 voxels), even though the training targets live in a space of at
+most ``M`` dimensions: every stimulus is presented several times, so the target
+matrix is ``M`` unique feature vectors repeated over ``N`` trials.
+
+Notation
+--------
+
+``X`` (N x p)
+    Raw brain data, one row per fMRI trial, with labels ``l_1 ... l_N``.
+``L_1 ... L_M``
+    The unique stimulus labels, ``np.unique(brain_labels)``.
+``R`` (N x M)
+    One-hot assignment, ``R[i, j] = 1`` iff ``l_i == L_j``.  Every row sums to
+    one, i.e. ``R @ 1 == 1``.
+``F`` (M x d)
+    Features of the unique stimuli, rows ordered like ``L``.
+``m``, ``s`` (1 x d)
+    Feature mean / SD over the ``M`` unique rows (``ddof=1``).
+
+The direct decoder fits ``Ridge(alpha)`` on ``X_n = (X - x_mean) / x_norm`` and
+``Y_n = (R F - m) / s``, then un-normalizes the prediction as ``* s + m``.  Note
+``Y_n = R F_n`` with ``F_n = (F - m) / s``, exactly, because ``R @ 1 == 1``.
+
+Equivalence
+-----------
+
+With ``fit_intercept=True`` the sklearn Ridge prediction is affine in the
+target.  Writing ``X_c = X_n - 1 x_bar_n`` and
+
+    S = (X_t - 1 x_bar_n) (X_c^T X_c + alpha I)^-1 X_c^T        (n_test x N)
+
+which depends on the brain data only,
+
+    predict(X_t) = S Y + (1 - S 1) y_bar,    y_bar = (1/N) 1^T Y
+
+Substituting ``Y = R F_n`` and ``y_bar = r_bar F_n`` with ``r_bar = (1/N) 1^T R``,
+
+    predict = [ S R + (1 - S 1) r_bar ] F_n = C F_n
+
+and ``C = S R + (1 - S 1) r_bar`` is precisely what the *same* estimator returns
+when it is fitted on ``R`` instead:
+
+    C = Ridge(alpha).fit(X_n, R).predict(X_t)                   (n_test x M)
+
+Row sums: ``C 1 = S R 1 + (1 - S 1)(r_bar 1) = S 1 + (1 - S 1) = 1``.  The final
+un-normalization therefore collapses:
+
+    (C F_n) * s + m = C (F - 1 m) + m = C F - (C 1) m + m = C F
+
+So the whole decoder is ``Ridge(alpha).fit(X_n, R).predict(X_t) @ F``: the
+feature mean/SD normalization cancels exactly, the meaning of ``alpha`` is
+unchanged (same ``X``, same penalty), and the intercept is still fitted and
+applied -- just on the ``M``-dimensional target.
+
+What is stored shrinks from ``d_out x d_in`` to ``M x d_in`` and no longer
+depends on the feature layer, so a single fit per (subject, ROI) serves every
+layer.  ``F`` itself is not copied into the decoder: it already exists on disk
+as the training feature directory, and prediction loads it from there.
+
+Note also that the fit is on ``R``, which is built from the trial labels alone,
+so *training never reads* ``F``: what it learns is the map from brain activity
+to coefficients over the training stimuli, plus the identity and order of those
+stimuli.  So the artifact records those labels and nothing about the features,
+and the one thing it has to guarantee is
+
+    coefficient column i is multiplied by the feature row of training stimulus i
+
+Three cheap checks cover it: ``verify_train_labels`` (the stored label list is
+the one the coefficients were fitted with, in order),
+``TrainingFeatureLoader.get`` (exactly one feature row from exactly one store
+per requested label, assembled in the requested order), and
+``check_column_correspondence`` at the contraction.
+
+This is a statement about the *fit*, not a licence to point one decoder
+directory at several feature sets: the surrounding pipeline stores one set of
+feature statistics per decoder directory for ``evaluation.py``, so one training
+feature configuration means one decoder directory.  The layout follows from the
+same two facts -- the model does not depend on the layer, the statistics do:
+
+    <decoder>/<subject>/<roi>/model/          the brain-side model, once
+    <decoder>/<layer>/<subject>/<roi>/model/  y_mean.mat, y_norm.mat
+
+The statistics are compatibility artifacts for ``evaluation.py``, not
+parameters of the model -- prediction never reads them, because the feature
+normalization cancels -- and they are written by ``predict_feature.py``, the
+step that has the features.
+'''
+
+
+from typing import Any, Dict, List, Optional, Sequence, Tuple
+
+import hashlib
+import json
+import os
+import pickle
+import uuid
+
+import numpy as np
+import yaml
+from bdpy.dataform import Features, save_array
+from bdpy.util import makedir_ifnot
+from sklearn.linear_model import Ridge
+
+
+# Model artifact ##############################################################
+
+FORMAT_NAME = 'sklearn_ridge_factorized'
+FORMAT_VERSION = 1
+
+MODEL_FILE = 'model.pkl.gz'          # bdpy ModelTraining pickle layout
+LABELS_FILE = 'train_labels.json'
+META_FILE = 'factorized.yaml'
+FEATURE_INDEX_FILE = 'feature_index.mat'
+
+# Identifies how `labels_digest` canonicalizes its input, so the recipe can
+# change later without silently comparing incomparable digests.
+LABELS_DIGEST_SPEC = 'sha256/labels-v1'
+
+PICKLE_PROTOCOL = 4                  # matches bdpy.ml.ModelTraining
+
+# bdpy ModelTraining/ModelTest chunk the target only when Y.ndim >= chunk_ndim + 1.
+CHUNK_NDIM = 2
+
+
+# Brain-side model ############################################################
+
+def one_hot_assignment(brain_labels: Sequence[Any],
+                       unique_labels: Sequence[Any],
+                       dtype: Any = np.float64) -> np.ndarray:
+    '''Return the assignment matrix ``R`` (n_trials x n_unique_stimuli).
+
+    ``R[i, j] == 1`` iff trial ``i`` presented ``unique_labels[j]``.  This is
+    the same alignment bdpy's ``Y_sort`` index expresses; expressed as a matrix
+    it makes ``Y == R @ F`` explicit.
+    '''
+    unique_labels = list(unique_labels)
+    position = {label: j for j, label in enumerate(unique_labels)}
+    if len(position) != len(unique_labels):
+        raise ValueError('unique_labels contains duplicates')
+
+    assignment = np.zeros((len(brain_labels), len(unique_labels)), dtype=dtype)
+    for i, label in enumerate(brain_labels):
+        try:
+            assignment[i, position[label]] = 1
+        except KeyError:
+            raise ValueError('Label %r of trial %d is not in unique_labels'
+                             % (label, i))
+    return assignment
+
+
+def normalize_brain_for_training(brain: np.ndarray, mean: np.ndarray,
+                                 std: np.ndarray) -> np.ndarray:
+    '''z-score the brain data the way ``bdpy.ml.ModelTraining`` does.
+
+    Infinities produced by zero-variance voxels are replaced by zero, matching
+    ``ModelTraining.run()``.  NaNs are deliberately left alone: the original
+    implementation does not clean them either.
+    '''
+    with np.errstate(invalid='ignore', divide='ignore'):
+        normalized = (brain - mean) / std
+    normalized[np.isinf(normalized)] = 0
+    return normalized
+
+
+def fit_stimulus_ridge(brain_normalized: np.ndarray,
+                       brain_labels: Sequence[Any],
+                       unique_labels: Sequence[Any],
+                       alpha: float = 100,
+                       dtype: Any = np.float32) -> Ridge:
+    '''Fit ``Ridge(alpha)`` from normalized brain data onto the stimulus basis.
+
+    The target is the one-hot assignment matrix ``R``, so the fitted model maps
+    a brain pattern to ``M`` coefficients, one per unique training stimulus.
+    It does not depend on the feature layer.
+    '''
+    assignment = one_hot_assignment(brain_labels, unique_labels)
+    model = Ridge(alpha=alpha)
+    model.fit(brain_normalized.astype(dtype), assignment.astype(dtype))
+    return model
+
+
+# Feature-side combination ####################################################
+
+def combine_features(coefficients: np.ndarray, features: np.ndarray,
+                     chunk_axis: Optional[int] = None) -> np.ndarray:
+    '''Compute ``coefficients @ features`` over the stimulus axis.
+
+    ``features`` has the training features of the unique stimuli along axis 0
+    and may be multidimensional (e.g. ``(M, channels, height, width)``).  When
+    ``chunk_axis`` is given and the features are at least 3-D, the product is
+    evaluated one index at a time along that axis and written straight into the
+    output, so no intermediate copy of the whole product is created.  The
+    contraction runs over axis 0 only, so every output element is
+    mathematically independent of the blocking; it is not bit-identical, since
+    BLAS may accumulate a width-1 slice in a different order than the full
+    array (a float32 rounding difference, ~1e-7 relative).
+    '''
+    coefficients = np.asarray(coefficients)
+    features = np.asarray(features)
+
+    if coefficients.ndim != 2:
+        raise ValueError('coefficients must be 2-D, got %d dimensions'
+                         % coefficients.ndim)
+    if features.ndim < 2:
+        raise ValueError('features must be at least 2-D, got %d dimensions'
+                         % features.ndim)
+    if coefficients.shape[1] != features.shape[0]:
+        raise ValueError(
+            'coefficients has %d stimulus columns but features has %d rows'
+            % (coefficients.shape[1], features.shape[0]))
+
+    if chunk_axis is None or features.ndim < CHUNK_NDIM + 1:
+        return np.tensordot(coefficients, features, axes=([1], [0]))
+
+    out_shape = (coefficients.shape[0],) + features.shape[1:]
+    out = np.empty(out_shape, dtype=np.result_type(coefficients, features))
+    selector: List[Any] = [slice(None)] * out.ndim
+    for i_chunk in range(features.shape[chunk_axis]):
+        chunk = np.take(features, [i_chunk], axis=chunk_axis)
+        selector[chunk_axis] = slice(i_chunk, i_chunk + 1)
+        out[tuple(selector)] = np.tensordot(coefficients, chunk,
+                                            axes=([1], [0]))
+    return out
+
+
+# Feature selection ###########################################################
+
+def resolve_feature_index_path(features_path: str,
+                               feature_index_file: str) -> str:
+    '''Where the feature index of ``features_path`` lives.
+
+    The single resolution rule for the whole pipeline, so that the loader and
+    the copy of the index cannot disagree: the index belongs to the feature
+    store and is resolved inside ``features_path`` unless it is absolute.
+    '''
+    if os.path.isabs(feature_index_file):
+        return feature_index_file
+    return os.path.join(features_path, feature_index_file)
+
+
+# Feature statistics ##########################################################
+
+def feature_statistics(features: np.ndarray):
+    '''``(mean, std)`` of the training features, shaped ``(1, *shape)``.
+
+    The same definition the direct implementation used: per-unit statistics over
+    the **unique** stimuli, with ``ddof=1``.
+    '''
+    mean = np.mean(features, axis=0)[np.newaxis, :]
+    std = np.std(features, axis=0, ddof=1)[np.newaxis, :]
+    return mean, std
+
+
+def save_array_atomic(path: str, array: np.ndarray, key: str,
+                      dtype: Any = np.float32) -> None:
+    '''Write a ``.mat`` array so that concurrent writers cannot tear it.
+
+    The feature statistics are materialized during prediction, and parallel
+    prediction workers against one decoder are a supported mode, so two
+    processes can reach the same target path.  Writing to a unique temporary
+    name in the same directory and then ``os.replace``-ing it is atomic within a
+    filesystem: a reader sees either the old file or the complete new one, never
+    a partial write, and a failed write leaves any existing file untouched.
+    '''
+    directory = os.path.dirname(os.path.abspath(path))
+    makedir_ifnot(directory)
+    # The name must keep the '.mat' suffix: hdf5storage appends one otherwise,
+    # and the file would not be where os.replace looks for it.
+    temporary = os.path.join(
+        directory,
+        '.tmp-%d-%s-%s' % (os.getpid(), uuid.uuid4().hex[:8],
+                           os.path.basename(path)))
+    try:
+        save_array(temporary, array, key=key, dtype=dtype, sparse=False)
+        os.replace(temporary, path)
+    except BaseException:
+        if os.path.exists(temporary):
+            os.unlink(temporary)
+        raise
+
+
+# Stimulus correspondence #####################################################
+
+def labels_digest(labels: Sequence[Any]) -> str:
+    '''Digest of the **ordered** training labels.
+
+    The stored coefficients are indexed by training stimulus, so the label order
+    is what makes ``coeff @ F`` meaningful.  The labels are already part of the
+    artifact, so pinning them costs nothing and does not read any feature file.
+    A reordered or edited ``train_labels.json`` no longer matches.
+    '''
+    digest = hashlib.sha256()
+    digest.update(b'featdec-labels-v1|')
+    for label in labels:
+        digest.update(str(label).encode('utf-8'))
+        digest.update(b'\0')
+    return digest.hexdigest()
+
+
+def verify_train_labels(metadata: Dict[str, Any], train_labels: Sequence[Any],
+                        model_dir: str) -> None:
+    '''Check the stored labels are the ones the coefficients were fitted with.'''
+    spec = metadata.get('labels_digest_spec')
+    if spec != LABELS_DIGEST_SPEC:
+        raise ValueError(
+            '%s records its label digest as %r, but this code computes %r.'
+            % (model_dir, spec, LABELS_DIGEST_SPEC))
+
+    actual = labels_digest(train_labels)
+    if actual != metadata['labels_digest']:
+        raise ValueError(
+            'The training labels of %s do not match the order its coefficients '
+            'were fitted with (digest %s, expected %s). Each coefficient column '
+            'belongs to one training stimulus, so a reordered or edited %s '
+            'would silently combine the wrong features.'
+            % (model_dir, actual[:16], metadata['labels_digest'][:16],
+               LABELS_FILE))
+
+
+def check_column_correspondence(n_columns: int, train_labels: Sequence[Any],
+                                features: np.ndarray,
+                                model_dir: str) -> None:
+    '''Assert the decoder's only requirement, at the point it is relied on.
+
+    Coefficient column ``i`` belongs to training stimulus ``i``, so it has to
+    meet the feature row of that stimulus.  The label order is checked by
+    ``verify_train_labels`` and the row assembly by
+    ``TrainingFeatureLoader.get``; this catches the remaining way the three can
+    disagree -- a model whose width is not the number of stored labels -- with a
+    readable message instead of a numpy shape error inside the contraction.
+    '''
+    n_labels = len(train_labels)
+    n_rows = int(np.asarray(features).shape[0])
+    if n_columns == n_labels == n_rows:
+        return
+    raise ValueError(
+        'Coefficient column i must be combined with the feature row of '
+        'training stimulus i, but %s predicts %d coefficient(s) while %s lists '
+        '%d training stimulus/stimuli and %d feature row(s) were loaded.'
+        % (model_dir, n_columns, LABELS_FILE, n_labels, n_rows))
+
+
+# Training feature access #####################################################
+
+class TrainingFeatureLoader(object):
+    '''Load unique-stimulus training features, one layer at a time.
+
+    Prediction needs the training features ``F`` of the layer being decoded.
+    Holding one layer costs the same memory the training step already needs, but
+    holding several would not scale, so this cache keeps **exactly one** layer
+    resident: requesting a different layer releases the previous array first.
+    '''
+
+    def __init__(self, features_paths: Sequence[str],
+                 feature_index_file: Optional[str] = None):
+        self._paths = list(features_paths)
+        self._feature_index_file = feature_index_file
+        self._stores: Optional[List[Features]] = None
+        self._cache_key: Optional[Tuple[str, Tuple[Any, ...]]] = None
+        self._cache: Optional[np.ndarray] = None
+        self.load_count = 0  # number of times features were read from disk
+
+    @property
+    def cached_layers(self) -> List[str]:
+        '''Layers currently held in memory (never more than one).'''
+        return [] if self._cache_key is None else [self._cache_key[0]]
+
+    def features(self) -> List[Features]:
+        if self._stores is None:
+            if self._feature_index_file is None:
+                self._stores = [Features(path) for path in self._paths]
+            else:
+                self._stores = [
+                    Features(path, feature_index=resolve_feature_index_path(
+                        path, self._feature_index_file))
+                    for path in self._paths]
+        return self._stores
+
+    def get(self, layer: str, labels: Sequence[Any]) -> np.ndarray:
+        '''Features of ``labels`` for ``layer``, shape ``(len(labels), ...)``.
+
+        Row ``i`` is the features of ``labels[i]`` -- which is what makes the
+        coefficients meaningful, so it is enforced per label rather than checked
+        in aggregate.  ``bdpy.dataform.utils.get_multi_features`` cannot be used
+        for this: it appends a row for *every* store holding a label and skips
+        labels no store holds, so one duplicated and one missing label cancel
+        out in the total row count while every row after them is shifted against
+        its coefficient column.
+        '''
+        key = (layer, tuple(labels))
+        if self._cache_key == key:
+            assert self._cache is not None
+            return self._cache
+
+        # Release the previous layer before allocating the next one.
+        self.release()
+
+        rows = []
+        for label in labels:
+            sources = [store for store in self.features()
+                       if label in store.labels]
+            if len(sources) != 1:
+                raise ValueError(
+                    'Training stimulus %r must come from exactly one feature '
+                    'store, but %d of %d store(s) provide it: %s'
+                    % (label, len(sources), len(self._paths),
+                       ', '.join(self._paths)))
+            row = sources[0].get(layer=layer, label=label)
+            if row.shape[0] != 1:
+                raise ValueError(
+                    'Training stimulus %r contributes %d feature rows for '
+                    'layer %r; one row per stimulus is required.'
+                    % (label, row.shape[0], layer))
+            rows.append(row)
+
+        features = np.vstack(rows)
+        self._cache_key = key
+        self._cache = features
+        self.load_count += 1
+        return features
+
+    def release(self) -> None:
+        '''Drop the cached layer.'''
+        self._cache = None
+        self._cache_key = None
+
+
+# Serialization ###############################################################
+
+def brain_model_dir(decoder_path: str, subject: str, roi: str) -> str:
+    """Where the brain-side model of one (subject, ROI) lives."""
+    return os.path.join(decoder_path, subject, roi, 'model')
+
+
+def layer_model_dir(decoder_path: str, layer: str, subject: str,
+                    roi: str) -> str:
+    """The per-layer directory: a legacy decoder, or the statistics sidecars.
+
+    ``evaluation.py`` reads ``y_mean``/``y_norm`` from exactly this path.
+    """
+    return os.path.join(decoder_path, layer, subject, roi, 'model')
+
+
+def is_factorized_model_dir(model_dir: str) -> bool:
+    '''True if ``model_dir`` holds a factorized decoder.
+
+    The presence of ``factorized.yaml`` is the only discriminator, and it is
+    checked on the **brain-side** directory. A factorized decoder's per-layer
+    directories contain nothing but the statistics sidecars, which from the
+    outside look exactly like a legacy decoder directory, so deciding by the
+    per-layer path would misread a decoder as legacy as soon as prediction has
+    run against it once.
+    '''
+    return os.path.isfile(os.path.join(model_dir, META_FILE))
+
+
+def save_factorized_model(model_dir: str, model: Ridge,
+                          train_labels: Sequence[Any],
+                          metadata: Optional[Dict[str, Any]] = None) -> List[str]:
+    '''Write the factorized decoder into ``model_dir``.
+
+    The estimator itself is stored in bdpy's ``ModelTraining`` pickle layout
+    (``{'model': ..., 'y_shape': ...}``) so that ``bdpy.ml.ModelTest`` reads it
+    without modification.  ``train_labels`` records the stimulus order, i.e. the
+    row order of ``F`` the coefficients refer to.  ``factorized.yaml`` marks the
+    format and is what distinguishes a factorized decoder from a legacy one.
+    '''
+    makedir_ifnot(model_dir)
+    train_labels = [str(label) for label in train_labels]
+
+    model_path = os.path.join(model_dir, MODEL_FILE)
+    with open(model_path, 'wb') as f:
+        pickle.dump({'model': model, 'y_shape': (len(train_labels),)}, f,
+                    protocol=PICKLE_PROTOCOL)
+
+    labels_path = os.path.join(model_dir, LABELS_FILE)
+    with open(labels_path, 'w') as f:
+        json.dump(train_labels, f, indent=1)
+
+    meta: Dict[str, Any] = {
+        'format': FORMAT_NAME,
+        'version': FORMAT_VERSION,
+        'model_file': MODEL_FILE,
+        'train_labels_file': LABELS_FILE,
+        'n_train_stimuli': len(train_labels),
+    }
+    meta.update(metadata or {})
+    meta_path = os.path.join(model_dir, META_FILE)
+    with open(meta_path, 'w') as f:
+        f.write(yaml.dump(meta, default_flow_style=False, sort_keys=True))
+
+    return [model_path, labels_path, meta_path]
+
+
+def load_factorized_metadata(model_dir: str) -> Dict[str, Any]:
+    '''Read ``factorized.yaml``.'''
+    with open(os.path.join(model_dir, META_FILE), 'r') as f:
+        meta = yaml.safe_load(f)
+    if not isinstance(meta, dict) or meta.get('format') != FORMAT_NAME:
+        raise ValueError('%s is not a %s decoder' % (model_dir, FORMAT_NAME))
+    if meta.get('version') != FORMAT_VERSION:
+        raise ValueError(
+            'Unsupported %s version %r in %s (this code writes version %d)'
+            % (FORMAT_NAME, meta.get('version'), model_dir, FORMAT_VERSION))
+    return meta
+
+
+def load_train_labels(model_dir: str) -> List[str]:
+    '''Read the stimulus order the stored coefficients refer to.'''
+    meta = load_factorized_metadata(model_dir)
+    labels_file = meta.get('train_labels_file', LABELS_FILE)
+    with open(os.path.join(model_dir, labels_file), 'r') as f:
+        labels = json.load(f)
+    if len(labels) != meta['n_train_stimuli']:
+        raise ValueError(
+            '%s lists %d labels but the metadata says %d'
+            % (labels_file, len(labels), meta['n_train_stimuli']))
+    return [str(label) for label in labels]
+
+
+def model_file_path(model_dir: str) -> str:
+    '''Path of the pickled brain-side Ridge model.'''
+    meta = load_factorized_metadata(model_dir)
+    return os.path.join(model_dir, meta.get('model_file', MODEL_FILE))
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/<analysis_name>.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 0000000..fa374e9
Binary files /dev/null and b/tests/data/golden/sklearn_ridge_pipeline.npz differ
diff --git a/tests/generate_golden.py b/tests/generate_golden.py
new file mode 100644
index 0000000..35abb50
--- /dev/null
+++ b/tests/generate_golden.py
@@ -0,0 +1,93 @@
+"""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)
+        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)
+                    values = pipeline.read_norm_params(decoder_dir, subject,
+                                                       roi)
+                    values.update(pipeline.read_feature_statistics(
+                        decoder_dir, layer, subject, roi))
+                    for key, value in values.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..0f14687
--- /dev/null
+++ b/tests/helpers/legacy_ridge.py
@@ -0,0 +1,203 @@
+"""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
+
+# What the *direct* implementation stored per decoder. The factorized decoder
+# keeps only the brain half (`pipeline.BRAIN_NORM_KEYS`); the feature half is a
+# property of the features a prediction combines and lives next to the decoded
+# features.
+NORM_KEYS = ('x_mean', 'x_norm', 'y_mean', 'y_norm')
+
+
+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)
+
+
+def write_legacy_decoder(model_dir: str, trained: Dict[str, Any]) -> None:
+    """Write ``trained`` in the pre-factorization on-disk decoder format.
+
+    One pickle per chunk (``%08d.pkl.gz``, or ``model.pkl.gz`` when the target
+    was not chunked) holding ``{'model': ..., 'y_shape': ...}``, the four
+    normalization ``.mat`` files, and an ``info.yaml`` completion marker.  Used
+    to test that the current prediction script can still read decoders produced
+    by the direct implementation.
+    """
+    import os
+    import pickle
+
+    import yaml
+    from bdpy.dataform import save_array
+
+    os.makedirs(model_dir, exist_ok=True)
+
+    for key in ('x_mean', 'x_norm', 'y_mean', 'y_norm'):
+        save_array(os.path.join(model_dir, key + '.mat'), trained[key],
+                   key=key, dtype=np.float32, sparse=False)
+
+    for i, entry in enumerate(trained['models']):
+        name = ('%08d.pkl.gz' % i) if trained['chunking'] else 'model.pkl.gz'
+        with open(os.path.join(model_dir, name), 'wb') as f:
+            pickle.dump({'model': entry['model'], 'y_shape': entry['y_shape']},
+                        f, protocol=4)
+
+    with open(os.path.join(model_dir, 'info.yaml'), 'w') as f:
+        f.write(yaml.dump({'_status': {'computation_id': 'legacy',
+                                       'computation_status': 'done'}},
+                          default_flow_style=False))
diff --git a/tests/helpers/pipeline.py b/tests/helpers/pipeline.py
new file mode 100644
index 0000000..2e1705f
--- /dev/null
+++ b/tests/helpers/pipeline.py
@@ -0,0 +1,156 @@
+"""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 evaluation
+import predict_feature
+import ridge_factorization
+import train_decoder_sklearn_ridge
+
+BRAIN_NORM_KEYS = ('x_mean', 'x_norm')
+FEATURE_NORM_KEYS = ('y_mean', 'y_norm')
+
+
+def run_training(dataset, output_dir: str, alpha: float = 100,
+                 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``.
+
+    No feature-side argument at all -- not even the layers: the factorized fit
+    is on the stimulus basis and is shared by every layer.
+    """
+    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,
+        output_dir=output_dir,
+        rois=dict(dataset.rois) if rois is None else dict(rois),
+        label_key=dataset.label_key,
+        alpha=alpha,
+        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] = (),
+                   features_paths: Optional[Sequence[str]] = None,
+                   feature_index_file: Optional[str] = None,
+                   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.  Pass
+    ``features_paths=[]`` to run without them.
+    """
+    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,
+        feature_index_file=feature_index_file,
+        analysis_name=analysis_name,
+    )
+    if 'training_features_paths' in inspect.signature(
+            predict_feature.featdec_predict).parameters:
+        # An explicit empty list means "run without the training features",
+        # which is a supported case for an already-finished prediction.
+        if features_paths is None:
+            features_paths = [dataset.train_features_dir]
+        kwargs['training_features_paths'] = list(features_paths)
+    return predict_feature.featdec_predict(fmri, decoder_dir, **kwargs)
+
+
+def brain_model_dir(decoder_dir: str, subject: str, roi: str) -> str:
+    """The shared brain-side model of a factorized decoder."""
+    return ridge_factorization.brain_model_dir(decoder_dir, subject, roi)
+
+
+def layer_model_dir(decoder_dir: str, layer: str, subject: str,
+                    roi: str) -> str:
+    """A legacy decoder directory, or a factorized decoder's sidecars."""
+    return ridge_factorization.layer_model_dir(decoder_dir, layer, subject, roi)
+
+
+def read_norm_params(decoder_dir: str, subject: str, roi: str,
+                     keys: Sequence[str] = BRAIN_NORM_KEYS
+                     ) -> Dict[str, np.ndarray]:
+    """Load the brain-side normalization parameters of a factorized decoder."""
+    directory = brain_model_dir(decoder_dir, subject, roi)
+    return {key: load_array(os.path.join(directory, key + '.mat'), key=key)
+            for key in keys}
+
+
+def read_feature_statistics(decoder_dir: str, layer: str, subject: str,
+                            roi: str) -> Dict[str, np.ndarray]:
+    """Load the per-layer statistics sidecars ``evaluation.py`` reads."""
+    directory = layer_model_dir(decoder_dir, layer, subject, roi)
+    return {key: load_array(os.path.join(directory, key + '.mat'), key=key)
+            for key in FEATURE_NORM_KEYS}
+
+
+def read_legacy_norm_params(decoder_dir: str, layer: str, subject: str,
+                            roi: str) -> Dict[str, np.ndarray]:
+    """All four parameters of a decoder in the direct (per-layer) format."""
+    directory = layer_model_dir(decoder_dir, layer, subject, roi)
+    return {key: load_array(os.path.join(directory, key + '.mat'), key=key)
+            for key in BRAIN_NORM_KEYS + FEATURE_NORM_KEYS}
+
+
+def run_evaluation(dataset, decoder_dir: str, decoded_dir: str,
+                   true_features_dir: str,
+                   layers: Optional[Sequence[str]] = None,
+                   subjects: Optional[Sequence[str]] = None,
+                   rois: Optional[Dict[str, str]] = None,
+                   output_file: Optional[str] = None,
+                   feature_index_file: Optional[str] = None,
+                   average_sample: bool = True) -> str:
+    """Evaluate decoded features with ``evaluation.py``."""
+    return evaluation.featdec_eval(
+        decoded_dir,
+        true_features_dir,
+        output_file=output_file or os.path.join(decoded_dir, 'evaluation.db'),
+        subjects=list(subjects or dataset.test_fmri),
+        rois=list(rois or dataset.rois),
+        layers=list(dataset.layers if layers is None else layers),
+        feature_index_file=feature_index_file,
+        feature_decoder_path=decoder_dir,
+        average_sample=average_sample,
+    )
+
+
+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..ea7e28e
--- /dev/null
+++ b/tests/helpers/synthetic.py
@@ -0,0 +1,274 @@
+"""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 ``<root>/<layer>/<label>.mat`` (key ``feat``)."""
+    for layer in layer_shapes:
+        layer_dir = os.path.join(root, layer)
+        os.makedirs(layer_dir, exist_ok=True)
+        for i, label in enumerate(labels):
+            save_array(os.path.join(layer_dir, label + '.mat'),
+                       arrays[layer][i][np.newaxis],
+                       key='feat', dtype=np.float32, sparse=False)
+
+
+def make_dataset(root: str, seed: int = 0,
+                 constant_feature_unit: bool = False) -> SyntheticDataset:
+    """Create a complete synthetic dataset under ``root``.
+
+    Parameters
+    ----------
+    root
+        Directory to create the dataset in.
+    seed
+        Seed for the random voxel/feature values.
+    constant_feature_unit
+        If True, one feature unit of every layer is constant across the training
+        stimuli, i.e. it has zero training variance.  Used to pin the behavior of
+        the decoder for such units.
+    """
+    rng = np.random.RandomState(seed)
+
+    n_train = len(TRAIN_REPEATS)
+    n_test = len(TEST_REPEATS)
+    unique_train_labels = ['train%04d' % (i + 1) for i in range(n_train)]
+    unique_test_labels = ['test%04d' % (i + 1) for i in range(n_test)]
+    label_to_number = {lb: i + 1 for i, lb in
+                       enumerate(unique_train_labels + unique_test_labels)}
+
+    train_trial_labels = _expand(unique_train_labels, TRAIN_REPEATS)
+    test_trial_labels = _expand(unique_test_labels, TEST_REPEATS)
+
+    # Features of the unique training stimuli.
+    features: Dict[str, np.ndarray] = {}
+    for layer, shape in LAYER_SHAPES.items():
+        feat = rng.randn(n_train, *shape).astype(np.float32) * 2.0 + 1.0
+        if constant_feature_unit:
+            flat = feat.reshape(n_train, -1)
+            flat[:, 0] = 7.0
+            feat = flat.reshape((n_train,) + shape)
+        features[layer] = feat
+
+    # Brain data, weakly coupled to the features so the decoder has signal.
+    train_brain: Dict[str, np.ndarray] = {}
+    train_labels: Dict[str, List[str]] = {}
+    test_brain: Dict[str, np.ndarray] = {}
+    test_labels: Dict[str, List[str]] = {}
+
+    driver = features['fc_like'].reshape(n_train, -1)
+    mixing = rng.randn(driver.shape[1], N_VOXELS)
+
+    for subject in ('sub-01', 'sub-02'):
+        idx = [unique_train_labels.index(lb) for lb in train_trial_labels]
+        signal = driver[idx] @ mixing
+        train_brain[subject] = (
+            signal + rng.randn(len(train_trial_labels), N_VOXELS) * 0.5)
+        train_labels[subject] = list(train_trial_labels)
+        test_brain[subject] = rng.randn(len(test_trial_labels), N_VOXELS) * 3.0
+        test_labels[subject] = list(test_trial_labels)
+
+    # sub-01's training trials are split over two files, sub-02's over one.
+    train_fmri: Dict[str, List[str]] = {}
+    split = 7
+    sub01 = [os.path.join(root, 'fmri', 'sub-01_train_run1.h5'),
+             os.path.join(root, 'fmri', 'sub-01_train_run2.h5')]
+    write_bdata(sub01[0], train_brain['sub-01'][:split],
+                train_labels['sub-01'][:split], label_to_number)
+    write_bdata(sub01[1], train_brain['sub-01'][split:],
+                train_labels['sub-01'][split:], label_to_number)
+    train_fmri['sub-01'] = sub01
+
+    sub02 = [os.path.join(root, 'fmri', 'sub-02_train.h5')]
+    write_bdata(sub02[0], train_brain['sub-02'], train_labels['sub-02'],
+                label_to_number)
+    train_fmri['sub-02'] = sub02
+
+    # ``predict_feature.py`` reads only the first file of each subject.
+    test_fmri: Dict[str, List[str]] = {}
+    for subject in ('sub-01', 'sub-02'):
+        path = os.path.join(root, 'fmri', '%s_test.h5' % subject)
+        write_bdata(path, test_brain[subject], test_labels[subject],
+                    label_to_number)
+        test_fmri[subject] = [path]
+
+    train_features_dir = os.path.join(root, 'features', 'train')
+    write_features(train_features_dir, LAYER_SHAPES, unique_train_labels,
+                   features)
+
+    return SyntheticDataset(
+        root=root,
+        train_fmri=train_fmri,
+        test_fmri=test_fmri,
+        train_features_dir=train_features_dir,
+        train_brain=train_brain,
+        train_labels=train_labels,
+        test_brain=test_brain,
+        test_labels=test_labels,
+        features=features,
+        unique_train_labels=unique_train_labels,
+        unique_test_labels=unique_test_labels,
+    )
+
+
+def add_layer(dataset: SyntheticDataset, name: str, shape: Tuple[int, ...],
+              seed: int = 1) -> np.ndarray:
+    """Add one more feature layer to an existing dataset's feature tree."""
+    rng = np.random.RandomState(seed)
+    n_stimuli = len(dataset.unique_train_labels)
+    array = rng.randn(n_stimuli, *shape).astype(np.float32)
+    write_features(dataset.train_features_dir, {name: shape},
+                   dataset.unique_train_labels, {name: array})
+    dataset.layer_shapes[name] = shape
+    dataset.features[name] = array
+    return array
+
+
+def write_test_features(dataset: SyntheticDataset, seed: int = 21) -> str:
+    """Write a ``Features`` tree for the *test* stimuli and return its path.
+
+    ``evaluation.py`` compares the decoded features against the true features
+    of the test stimuli, which the training dataset does not contain.
+    """
+    root = os.path.join(dataset.root, 'features', 'test')
+    write_feature_tree(root, dataset.layer_shapes, dataset.unique_test_labels,
+                       seed=seed)
+    return root
+
+
+def write_feature_index(path: str,
+                        index_by_layer: Dict[str, Sequence[int]]) -> None:
+    """Write a feature-index ``.mat`` as ``Features(feature_index=...)`` expects.
+
+    bdpy reads it as ``hdf5storage.loadmat(path)['index']`` and then indexes the
+    result by layer name, so ``index`` is a mapping from layer to a 1-D array of
+    zero-based column indices into the C-order-flattened features.
+    """
+    import hdf5storage
+
+    payload = {layer: np.asarray(index).ravel()
+               for layer, index in index_by_layer.items()}
+    os.makedirs(os.path.dirname(path) or '.', exist_ok=True)
+    hdf5storage.savemat(path, {'index': payload}, format='7.3',
+                        oned_as='column', store_python_metadata=True)
+
+
+def write_feature_tree(root: str, layer_shapes: Dict[str, Tuple[int, ...]],
+                       labels: Sequence[str], seed: int = 0
+                       ) -> Dict[str, np.ndarray]:
+    """Write an independent ``Features`` tree and return its arrays."""
+    rng = np.random.RandomState(seed)
+    arrays = {layer: rng.randn(len(labels), *shape).astype(np.float32)
+              for layer, shape in layer_shapes.items()}
+    write_features(root, layer_shapes, labels, arrays)
+    return arrays
diff --git a/tests/test_chunking.py b/tests/test_chunking.py
new file mode 100644
index 0000000..40de131
--- /dev/null
+++ b/tests/test_chunking.py
@@ -0,0 +1,195 @@
+"""Chunking of the feature dimension.
+
+Chunking used to be a *training* concern: bdpy fitted one Ridge per index along
+``chunk_axis`` and flattened each chunk with ``order='F'``.  In the factorized
+decoder the fitted model no longer touches the feature dimensions at all, so
+chunking is purely a prediction-time memory device: the linear combination with
+the training features is evaluated one index at a time along ``chunk_axis``.
+
+Either way ``chunk_axis`` must not change the result, and each feature position
+must depend only on its own training features.  The shipped config only enables
+``fc6``/``fc7``/``fc8``, which are 2-D and never reach this path, so it is
+pinned explicitly here.
+"""
+
+from __future__ import annotations
+
+import os
+import pickle
+
+import numpy as np
+import pytest
+
+from ridge_factorization import combine_features
+from tests.helpers import pipeline
+
+ALPHA = 100
+
+
+def _model_files(decoder_dir, layer, subject='sub-01', roi='VC'):
+    directory = pipeline.brain_model_dir(decoder_dir, subject, roi)
+    return sorted(f for f in os.listdir(directory) if f.endswith('.pkl.gz'))
+
+
+def test_one_model_per_layer_regardless_of_feature_dimensionality(dataset,
+                                                                  decoder_dir):
+    """No per-chunk model files: the stored model is the stimulus-basis Ridge."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    for layer in dataset.layers:
+        assert _model_files(decoder_dir, layer) == ['model.pkl.gz']
+
+
+def test_stored_model_targets_the_stimulus_basis(dataset, decoder_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    n_stimuli = len(np.unique(dataset.labels('sub-01', 'train')))
+    for layer in dataset.layers:
+        directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+        with open(os.path.join(directory, 'model.pkl.gz'), 'rb') as f:
+            payload = pickle.load(f)
+        assert payload['y_shape'] == (n_stimuli,)
+        n_voxels = dataset.brain('sub-01', 'VC').shape[1]
+        assert payload['model'].coef_.shape == (n_stimuli, n_voxels)
+
+
+@pytest.mark.parametrize('chunk_axis', [1, 2, None])
+def test_chunk_axis_does_not_change_the_prediction(dataset, tmp_path,
+                                                   chunk_axis):
+    """``chunk_axis`` blocks the computation; it must not move the numbers."""
+    layer, subject, roi = 'conv_like', 'sub-01', 'VC'
+
+    reference_decoder = str(tmp_path / 'decoders_ref')
+    reference_decoded = str(tmp_path / 'decoded_ref')
+    pipeline.run_training(dataset, reference_decoder, alpha=ALPHA,
+                          analysis_name='train_ref')
+    pipeline.run_prediction(dataset, reference_decoder, reference_decoded,
+                            chunk_axis=None, layers=[layer],
+                            analysis_name='predict_ref')
+    reference = pipeline.read_decoded_features(
+        reference_decoded, layer, subject, roi, dataset.unique_test_labels)
+
+    decoder_dir = str(tmp_path / ('decoders_%s' % chunk_axis))
+    decoded_dir = str(tmp_path / ('decoded_%s' % chunk_axis))
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA,
+                          analysis_name='train_%s' % chunk_axis)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=chunk_axis, layers=[layer],
+                            analysis_name='predict_%s' % chunk_axis)
+    actual = pipeline.read_decoded_features(
+        decoded_dir, layer, subject, roi, dataset.unique_test_labels)
+
+    np.testing.assert_allclose(actual, reference, rtol=1e-5, atol=1e-6)
+
+
+@pytest.mark.parametrize('chunk_axis', [1, 2, 3])
+def test_combine_features_matches_the_unblocked_contraction(chunk_axis):
+    """Blocking contracts over the same axis, to within rounding.
+
+    The contraction runs over the stimulus axis only, so every output element is
+    mathematically independent of the blocking.  It is not guaranteed to be
+    bit-identical: BLAS may accumulate a width-1 slice in a different order than
+    the full array, which shows up in float32 (see the float32 case below).
+    """
+    rng = np.random.RandomState(0)
+    coefficients = rng.randn(5, 7)
+    features = rng.randn(7, 3, 4, 2)
+
+    whole = combine_features(coefficients, features, chunk_axis=None)
+    blocked = combine_features(coefficients, features, chunk_axis=chunk_axis)
+
+    assert blocked.shape == whole.shape
+    np.testing.assert_allclose(blocked, whole, rtol=1e-12, atol=1e-12)
+
+
+def test_combine_features_blocking_error_stays_at_rounding_level():
+    """At realistic float32 sizes the two orders differ only by rounding."""
+    rng = np.random.RandomState(0)
+    coefficients = rng.randn(50, 600).astype(np.float32)
+    features = rng.randn(600, 512, 14, 14).astype(np.float32)
+
+    whole = combine_features(coefficients, features, chunk_axis=None)
+    blocked = combine_features(coefficients, features, chunk_axis=1)
+
+    np.testing.assert_allclose(blocked, whole, rtol=1e-5,
+                               atol=1e-5 * np.max(np.abs(whole)))
+
+
+def test_combine_features_ignores_chunk_axis_for_2d_features():
+    """2-D features are below bdpy's chunking threshold, as before."""
+    rng = np.random.RandomState(0)
+    coefficients = rng.randn(4, 6)
+    features = rng.randn(6, 9)
+
+    np.testing.assert_array_equal(
+        combine_features(coefficients, features, chunk_axis=1),
+        combine_features(coefficients, features, chunk_axis=None))
+
+
+def test_each_feature_channel_depends_only_on_its_own_features(dataset,
+                                                               tmp_path):
+    """A per-position shift of one channel reappears at exactly those positions.
+
+    The perturbation has to differ *within* the channel. A uniform shift of the
+    whole channel would still pass if a flatten/unflatten pair disagreed about C
+    versus Fortran order, because a transposition of identical values is
+    invisible; with distinct values per position, a mismatch would move the
+    response to the wrong position inside the channel.
+
+    A constant per-position shift passes through the decoder exactly: it moves
+    ``y_mean`` by ``delta`` and leaves ``y_norm`` untouched, so the normalized
+    target is unchanged and the un-normalization adds ``delta`` back. In the
+    factorized decoder the same follows from the stimulus coefficients summing
+    to one.
+    """
+    import shutil
+
+    from bdpy.dataform import load_array, save_array
+
+    layer, subject, roi = 'conv_like', 'sub-01', 'VC'
+
+    baseline_decoder = str(tmp_path / 'decoders_base')
+    baseline_decoded = str(tmp_path / 'decoded_base')
+    pipeline.run_training(dataset, baseline_decoder, alpha=ALPHA, analysis_name='train_base')
+    pipeline.run_prediction(dataset, baseline_decoder, baseline_decoded,
+                            chunk_axis=1, layers=[layer],
+                            analysis_name='predict_base')
+    baseline = pipeline.read_decoded_features(
+        baseline_decoded, layer, subject, roi, dataset.unique_test_labels)
+
+    # Copy the feature tree and perturb channel 1 only, by a different amount at
+    # each position within the channel.
+    delta = np.array([[1.0, 2.0], [4.0, 8.0]], dtype=np.float32)
+    assert delta.shape == dataset.layer_shapes[layer][1:]
+    assert not np.allclose(delta, delta.T), \
+        'delta must be asymmetric, or an order mix-up would be invisible'
+
+    perturbed_features = str(tmp_path / 'features_perturbed')
+    shutil.copytree(dataset.train_features_dir, perturbed_features)
+    for label in dataset.unique_train_labels:
+        path = os.path.join(perturbed_features, layer, '%s.mat' % label)
+        array = load_array(path, key='feat')
+        array[:, 1] = array[:, 1] + delta
+        save_array(path, array, key='feat', dtype=np.float32, sparse=False)
+
+    dataset.train_features_dir = perturbed_features
+    perturbed_decoder = str(tmp_path / 'decoders_perturbed')
+    perturbed_decoded = str(tmp_path / 'decoded_perturbed')
+    pipeline.run_training(dataset, perturbed_decoder, alpha=ALPHA,
+                          analysis_name='train_perturbed')
+    pipeline.run_prediction(dataset, perturbed_decoder, perturbed_decoded,
+                            chunk_axis=1, layers=[layer],
+                            analysis_name='predict_perturbed')
+    perturbed = pipeline.read_decoded_features(
+        perturbed_decoded, layer, subject, roi, dataset.unique_test_labels)
+
+    shift = perturbed - baseline
+
+    # Channel 1 moves by delta, position for position.
+    np.testing.assert_allclose(
+        shift[:, 1], np.broadcast_to(delta, shift[:, 1].shape),
+        rtol=0, atol=1e-4)
+
+    # Nothing else moves.
+    other = np.delete(shift, 1, axis=1)
+    np.testing.assert_allclose(other, np.zeros_like(other), rtol=0, atol=1e-4)
diff --git a/tests/test_evaluation.py b/tests/test_evaluation.py
new file mode 100644
index 0000000..67b5401
--- /dev/null
+++ b/tests/test_evaluation.py
@@ -0,0 +1,54 @@
+"""The factorized pipeline still satisfies the evaluation file contract.
+
+``evaluation.py`` is not changed by the factorization, and it reads the
+training-feature statistics from
+``<decoder>/<layer>/<subject>/<roi>/model/y_{mean,norm}.mat``.  The factorized
+decoder cannot write them at training time -- it never opens a feature file --
+so ``predict_feature.py`` fills them in.  This test exists to catch that
+contract breaking: nothing else in the suite runs ``evaluation.py``.
+"""
+
+from __future__ import annotations
+
+import os
+
+import numpy as np
+
+from tests.helpers import pipeline, synthetic
+
+ALPHA = 100
+CHUNK_AXIS = 1
+LAYER = 'fc_like'
+ONE_SUBJECT = ['sub-01']
+ONE_ROI = {'VC': 'ROI_VC = 1'}
+
+
+def test_train_predict_evaluate_end_to_end(dataset, tmp_path):
+    """Train, predict, then run the unchanged ``evaluation.py`` over it."""
+    decoder_dir = str(tmp_path / 'decoders')
+    decoded_dir = str(tmp_path / 'decoded')
+
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA,
+                          subjects=ONE_SUBJECT, rois=ONE_ROI)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS, layers=[LAYER],
+                            subjects=ONE_SUBJECT, rois=ONE_ROI)
+
+    # The statistics are where evaluation.py looks them up, unchanged:
+    # <decoder>/<layer>/<subject>/<roi>/model/.
+    for key in ('y_mean', 'y_norm'):
+        assert os.path.isfile(os.path.join(
+            pipeline.layer_model_dir(decoder_dir, LAYER, 'sub-01', 'VC'),
+            '%s.mat' % key))
+
+    output = pipeline.run_evaluation(
+        dataset, decoder_dir, decoded_dir, synthetic.write_test_features(dataset),
+        layers=[LAYER], subjects=ONE_SUBJECT, rois=ONE_ROI)
+
+    from evaluation import ResultsStore
+    results = ResultsStore(output)
+    for metric in ('profile_correlation', 'pattern_correlation',
+                   'identification_accuracy'):
+        value = results.get(layer=LAYER, subject='sub-01', roi='VC',
+                            metric=metric)
+        assert value is not None and np.asarray(value).size > 0, metric
diff --git a/tests/test_factorized_artifacts.py b/tests/test_factorized_artifacts.py
new file mode 100644
index 0000000..f0e0020
--- /dev/null
+++ b/tests/test_factorized_artifacts.py
@@ -0,0 +1,129 @@
+"""The factorized decoder artifact: metadata, layer independence, robustness."""
+
+from __future__ import annotations
+
+import json
+import os
+
+import numpy as np
+import pytest
+import yaml
+
+import ridge_factorization
+from tests.helpers import pipeline, synthetic
+
+ALPHA = 100
+ONE_SUBJECT = ['sub-01']
+ONE_ROI = {'VC': 'ROI_VC = 1'}
+
+
+def test_metadata_describes_the_format(dataset, decoder_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    with open(os.path.join(directory, 'factorized.yaml')) as f:
+        meta = yaml.safe_load(f)
+
+    assert meta['format'] == ridge_factorization.FORMAT_NAME
+    assert meta['version'] == ridge_factorization.FORMAT_VERSION
+    assert meta['alpha'] == float(ALPHA)
+    assert meta['subject'] == 'sub-01'
+    assert meta['roi'] == 'VC'
+    assert meta['dtype'] == 'float32'
+    assert meta['n_trials'] == len(dataset.labels('sub-01', 'train'))
+    assert meta['n_voxels'] == dataset.brain('sub-01', 'VC').shape[1]
+    assert meta['n_train_stimuli'] == len(
+        np.unique(dataset.labels('sub-01', 'train')))
+
+    # The decoder holds coefficients over the training stimuli, so the artifact
+    # pins their identity and order -- and nothing else about the features,
+    # which it never reads. Not even a layer: the model is shared by all of
+    # them.
+    labels = [str(lb) for lb in np.unique(dataset.labels('sub-01', 'train'))]
+    assert meta['labels_digest_spec'] == \
+        ridge_factorization.LABELS_DIGEST_SPEC
+    assert meta['labels_digest'] == ridge_factorization.labels_digest(labels)
+
+    for absent in ('layer', 'chunk_axis', 'feature_source', 'feature_shape',
+                   'feature_digest', 'feature_digest_spec',
+                   'feature_index_file'):
+        assert absent not in meta, absent
+
+
+def test_train_labels_record_the_feature_row_order(dataset, decoder_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    expected = [str(lb) for lb in np.unique(dataset.labels('sub-01', 'train'))]
+    directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    with open(os.path.join(directory, 'train_labels.json')) as f:
+        labels = json.load(f)
+    assert labels == expected
+    assert ridge_factorization.load_train_labels(directory) == expected
+
+
+def test_the_brain_side_model_is_stored_once(dataset, decoder_dir,
+                                             decoded_dir):
+    """The stimulus-basis Ridge does not depend on the layer, so it is not
+    copied under every layer -- the per-layer directories that prediction
+    creates hold only the statistics ``evaluation.py`` reads.
+    """
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir)
+
+    for subject in dataset.train_fmri:
+        for roi in dataset.rois:
+            assert os.path.isfile(os.path.join(
+                pipeline.brain_model_dir(decoder_dir, subject, roi),
+                'model.pkl.gz'))
+            for layer in dataset.layers:
+                sidecars = pipeline.layer_model_dir(decoder_dir, layer,
+                                                    subject, roi)
+                assert sorted(os.listdir(sidecars)) == ['y_mean.mat',
+                                                        'y_norm.mat']
+
+
+def test_model_size_does_not_grow_with_the_feature_dimension(dataset,
+                                                             decoder_dir,
+                                                             decoded_dir):
+    """The stored model scales with n_train_stimuli x n_voxels only.
+
+    This is the whole point of the factorization: the direct implementation
+    stored n_features x n_voxels coefficients, so a 50x larger layer meant a
+    50x larger model file. Here the layer does not enter the model at all --
+    only the sidecars grow with it.
+    """
+    small_units = int(np.prod(dataset.layer_shapes['fc_like']))
+    synthetic.add_layer(dataset, 'big_like', (small_units * 50,))
+
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            layers=['fc_like', 'big_like'])
+
+    directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    assert sorted(os.listdir(directory)) == [
+        'factorized.yaml', 'info.yaml', 'model.pkl.gz', 'train_labels.json',
+        'x_mean.mat', 'x_norm.mat',
+    ]
+
+    statistics = {
+        layer: pipeline.read_feature_statistics(decoder_dir, layer, 'sub-01',
+                                                'VC')['y_mean'].size
+        for layer in ('fc_like', 'big_like')}
+    assert statistics['big_like'] == 50 * statistics['fc_like']
+
+
+def test_metadata_version_is_checked(dataset, decoder_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA,
+                          subjects=ONE_SUBJECT, rois=ONE_ROI)
+
+    directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    meta_path = os.path.join(directory, 'factorized.yaml')
+    with open(meta_path) as f:
+        meta = yaml.safe_load(f)
+    meta['version'] = ridge_factorization.FORMAT_VERSION + 1
+    with open(meta_path, 'w') as f:
+        f.write(yaml.dump(meta))
+
+    with pytest.raises(ValueError, match='Unsupported'):
+        ridge_factorization.load_factorized_metadata(directory)
+
diff --git a/tests/test_label_alignment.py b/tests/test_label_alignment.py
new file mode 100644
index 0000000..d5e0ca7
--- /dev/null
+++ b/tests/test_label_alignment.py
@@ -0,0 +1,80 @@
+"""Alignment between fMRI trials and unique-stimulus features.
+
+The training script derives ``feat_labels = np.unique(brain_labels)``, loads one
+feature row per unique stimulus, and expands those rows back to trial order with
+``feat_index`` (bdpy's ``Y_sort``).  That expansion is exactly the one-hot
+assignment ``Y = R F`` the factorized decoder relies on, so it is pinned here.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from tests.helpers import legacy_ridge, pipeline
+
+
+def test_assignment_index_maps_every_trial_to_its_stimulus(dataset):
+    brain_labels = dataset.labels('sub-01', 'train')
+    feat_labels = list(np.unique(brain_labels))
+
+    index = legacy_ridge.assignment_index(brain_labels, feat_labels)
+
+    assert index.shape == (len(brain_labels),)
+    assert [feat_labels[i] for i in index] == list(brain_labels)
+
+
+def test_unique_labels_are_sorted_and_deduplicated(dataset):
+    brain_labels = dataset.labels('sub-01', 'train')
+    feat_labels = list(np.unique(brain_labels))
+
+    assert feat_labels == sorted(set(brain_labels))
+    assert len(feat_labels) < len(brain_labels)
+
+
+def test_repetitions_are_unbalanced_in_the_fixture(dataset):
+    """The fixture must keep unbalanced repeats, or several tests go blind."""
+    brain_labels = dataset.labels('sub-01', 'train')
+    counts = [brain_labels.count(lb) for lb in np.unique(brain_labels)]
+
+    assert len(set(counts)) > 1, counts
+
+
+def test_expansion_repeats_feature_rows_per_trial(dataset):
+    brain_labels = dataset.labels('sub-01', 'train')
+    feat_labels = list(np.unique(brain_labels))
+    feat = dataset.feature_matrix('conv_like', feat_labels)
+
+    index = legacy_ridge.assignment_index(brain_labels, feat_labels)
+    expanded = feat[index]
+
+    assert expanded.shape == (len(brain_labels),) + feat.shape[1:]
+    for trial, label in enumerate(brain_labels):
+        np.testing.assert_array_equal(
+            expanded[trial], feat[feat_labels.index(label)])
+
+
+def test_trial_order_does_not_matter_for_the_fitted_decoder(dataset,
+                                                            tmp_path):
+    """Permuting (brain trial, label) pairs together leaves predictions intact."""
+    from tests.helpers import synthetic
+
+    layer, subject, roi = 'fc_like', 'sub-01', 'VC'
+    brain = dataset.brain(subject, roi, 'train')
+    brain_labels = dataset.labels(subject, 'train')
+    feat_labels = list(np.unique(brain_labels))
+    feat = dataset.feature_matrix(layer, feat_labels)
+
+    trained = legacy_ridge.legacy_train(brain, brain_labels, feat, feat_labels,
+                                        alpha=100, chunk_axis=1)
+
+    order = np.random.RandomState(1).permutation(len(brain_labels))
+    shuffled = legacy_ridge.legacy_train(
+        brain[order], [brain_labels[i] for i in order], feat, feat_labels,
+        alpha=100, chunk_axis=1)
+
+    test_brain, _ = legacy_ridge.average_test_brain(
+        dataset.brain(subject, roi, 'test'), dataset.labels(subject, 'test'))
+    np.testing.assert_allclose(
+        legacy_ridge.legacy_predict(trained, test_brain),
+        legacy_ridge.legacy_predict(shuffled, test_brain),
+        rtol=1e-4, atol=1e-5)
diff --git a/tests/test_lazy_feature_statistics.py b/tests/test_lazy_feature_statistics.py
new file mode 100644
index 0000000..ac1a4a1
--- /dev/null
+++ b/tests/test_lazy_feature_statistics.py
@@ -0,0 +1,312 @@
+"""The training-feature statistics: compatibility sidecars, written lazily.
+
+``evaluation.py`` normalizes the pattern correlation with the per-unit mean and
+SD of the training features, read from
+``<decoder>/<layer>/<subject>/<roi>/model/``.  The factorized decoder cannot
+produce them: training never opens a feature file, and prediction does not need
+them either (the feature normalization cancels).  So ``predict_feature.py`` --
+the step that has the features -- writes them there, purely for the unchanged
+evaluation pipeline.
+
+Three things make that work: the write happens *before* the "results directory
+already exists" skip, so a run whose features were already decoded still gets
+them; it is atomic, because parallel prediction workers are a supported mode and
+two runs with different analysis names are not excluded by the DistComp lock;
+and the statistics are computed once per ``(layer, ordered training labels)``
+and then written into every decoder directory that needs them.  Existing files
+are never overwritten: one decoder directory belongs to one training feature
+configuration.
+"""
+
+from __future__ import annotations
+
+import os
+import threading
+
+import numpy as np
+import pytest
+from bdpy.dataform import load_array, save_array
+
+import predict_feature
+import ridge_factorization
+from ridge_factorization import save_array_atomic
+from tests.helpers import pipeline
+
+ALPHA = 100
+CHUNK_AXIS = 1
+LAYER = 'fc_like'
+ONE_SUBJECT = ['sub-01']
+ONE_ROI = {'VC': 'ROI_VC = 1'}
+KEYS = ('y_mean', 'y_norm')
+
+
+def _train(dataset, decoder_dir, **kwargs):
+    return pipeline.run_training(
+        dataset, decoder_dir, alpha=ALPHA, subjects=ONE_SUBJECT, rois=ONE_ROI, **kwargs)
+
+
+def _predict(dataset, decoder_dir, decoded_dir, **kwargs):
+    return pipeline.run_prediction(
+        dataset, decoder_dir, decoded_dir, chunk_axis=CHUNK_AXIS,
+        layers=[LAYER], subjects=ONE_SUBJECT, rois=ONE_ROI, **kwargs)
+
+
+def _statistics_paths(decoder_dir, roi='VC'):
+    directory = pipeline.layer_model_dir(decoder_dir, LAYER, 'sub-01', roi)
+    return {key: os.path.join(directory, '%s.mat' % key) for key in KEYS}
+
+
+# The materialization ##########################################################
+
+def test_a_second_prediction_run_does_not_rewrite_the_statistics(dataset,
+                                                                 tmp_path):
+    """The cheap existence check comes first, so nothing is recomputed."""
+    decoder_dir = str(tmp_path / 'decoders')
+    decoded_dir = str(tmp_path / 'decoded')
+    _train(dataset, decoder_dir)
+    _predict(dataset, decoder_dir, decoded_dir)
+
+    paths = _statistics_paths(decoder_dir)
+    before = {key: os.stat(path).st_mtime_ns for key, path in paths.items()}
+
+    _predict(dataset, decoder_dir, decoded_dir, analysis_name='again')
+
+    after = {key: os.stat(path).st_mtime_ns for key, path in paths.items()}
+    assert after == before
+
+
+def test_the_statistics_are_materialized_past_the_already_done_skip(
+        dataset, tmp_path, capsys):
+    """The reason the write runs before the skip.
+
+    A prediction whose features were already decoded would otherwise skip
+    forever, leaving ``evaluation.py`` permanently unable to run against it.
+    """
+    decoder_dir = str(tmp_path / 'decoders')
+    decoded_dir = str(tmp_path / 'decoded')
+    _train(dataset, decoder_dir)
+    _predict(dataset, decoder_dir, decoded_dir)
+
+    paths = _statistics_paths(decoder_dir)
+    for path in paths.values():
+        os.unlink(path)
+
+    decoded_before = sorted(
+        (name, os.stat(os.path.join(decoded_dir, LAYER, 'sub-01', 'VC',
+                                    name)).st_mtime_ns)
+        for name in os.listdir(os.path.join(decoded_dir, LAYER, 'sub-01',
+                                            'VC')))
+
+    capsys.readouterr()
+    _predict(dataset, decoder_dir, decoded_dir, analysis_name='rerun')
+
+    # Every (layer, subject, ROI) hit the skip ...
+    assert 'is already done. Skipped.' in capsys.readouterr().out
+    decoded_after = sorted(
+        (name, os.stat(os.path.join(decoded_dir, LAYER, 'sub-01', 'VC',
+                                    name)).st_mtime_ns)
+        for name in os.listdir(os.path.join(decoded_dir, LAYER, 'sub-01',
+                                            'VC')))
+    assert decoded_after == decoded_before
+
+    # ... and the statistics came back anyway, with the right values.
+    feat = dataset.feature_matrix(LAYER, dataset.unique_train_labels)
+    np.testing.assert_allclose(
+        load_array(paths['y_mean'], key='y_mean'),
+        np.mean(feat, axis=0)[np.newaxis].astype(np.float32),
+        rtol=1e-5, atol=1e-6)
+    np.testing.assert_allclose(
+        load_array(paths['y_norm'], key='y_norm'),
+        np.std(feat, axis=0, ddof=1)[np.newaxis].astype(np.float32),
+        rtol=1e-5, atol=1e-6)
+
+
+def test_a_finished_prediction_reruns_without_the_training_features(dataset,
+                                                                    tmp_path):
+    """Nothing left to decode, so nothing to verify or materialize either."""
+    decoder_dir = str(tmp_path / 'decoders')
+    decoded_dir = str(tmp_path / 'decoded')
+    _train(dataset, decoder_dir)
+    _predict(dataset, decoder_dir, decoded_dir)
+
+    _predict(dataset, decoder_dir, decoded_dir, features_paths=[],
+             analysis_name='no_features')
+
+
+def test_an_unwritable_output_directory_warns_instead_of_failing(dataset,
+                                                                 tmp_path,
+                                                                 monkeypatch):
+    """Decoding must still work when the statistics cannot be written.
+
+    The write is made to fail here rather than by permission bits, since the
+    test suite may run as a user for whom no permission bit denies anything.
+    """
+    decoder_dir = str(tmp_path / 'decoders')
+    decoded_dir = str(tmp_path / 'decoded')
+    _train(dataset, decoder_dir)
+
+    def read_only(*args, **kwargs):
+        raise IOError('Read-only file system')
+
+    monkeypatch.setattr(predict_feature, 'save_array_atomic', read_only)
+    with pytest.warns(UserWarning, match='read-only'):
+        _predict(dataset, decoder_dir, decoded_dir)
+
+    assert not os.path.exists(_statistics_paths(decoder_dir)['y_mean'])
+    assert pipeline.read_decoded_features(
+        decoded_dir, LAYER, 'sub-01', 'VC',
+        dataset.unique_test_labels).shape[0] == len(dataset.unique_test_labels)
+
+
+def test_the_statistics_are_computed_once_per_layer_and_label_set(
+        dataset, decoder_dir, decoded_dir, monkeypatch):
+    """Not once per decoder: they do not depend on the ROI.
+
+    At DeepRecon scale each computation is two passes over a 15 GB tensor, and
+    there are 9 ROIs per subject.
+    """
+    computed = []
+    original = ridge_factorization.feature_statistics
+
+    def spy(features):
+        computed.append(features.shape)
+        return original(features)
+
+    monkeypatch.setattr(predict_feature, 'feature_statistics', spy)
+
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    # Both synthetic subjects see the same stimuli, so one computation per
+    # layer covers every (subject, ROI) of the run ...
+    assert len(dataset.rois) > 1 and len(dataset.test_fmri) > 1
+    assert len(computed) == len(dataset.layers)
+
+    # ... while every decoder directory still gets its own copy, because that
+    # is where evaluation.py looks.
+    for layer in dataset.layers:
+        for subject in dataset.test_fmri:
+            for roi in dataset.rois:
+                assert pipeline.read_feature_statistics(decoder_dir, layer,
+                                                        subject, roi)
+
+
+def test_existing_statistics_are_never_overwritten(dataset, tmp_path):
+    """One decoder directory belongs to one training feature configuration.
+
+    Re-pointing it at another ``Y_tr`` or feature index is not a supported
+    workflow, so prediction fills the sidecars in when they are missing and
+    leaves them alone otherwise -- it never decides that what is there is
+    wrong.
+    """
+    decoder_dir = str(tmp_path / 'decoders')
+    _train(dataset, decoder_dir)
+    _predict(dataset, decoder_dir, str(tmp_path / 'decoded'))
+
+    paths = _statistics_paths(decoder_dir)
+    planted = np.full((1, 3), 7.0, dtype=np.float32)
+    for path, key in ((paths[key], key) for key in KEYS):
+        save_array(path, planted, key=key, dtype=np.float32, sparse=False)
+
+    _predict(dataset, decoder_dir, str(tmp_path / 'decoded_again'),
+             analysis_name='again')
+
+    for key in KEYS:
+        np.testing.assert_array_equal(load_array(paths[key], key=key), planted)
+
+
+# save_array_atomic ############################################################
+
+def _temp_files(directory):
+    return [name for name in os.listdir(directory) if name.startswith('.tmp')]
+
+
+def test_atomic_save_writes_a_loadable_file(tmp_path):
+    """The temp name has to keep the '.mat' suffix that hdf5storage appends."""
+    path = str(tmp_path / 'y_mean.mat')
+    array = np.arange(6, dtype=np.float32)[np.newaxis]
+
+    save_array_atomic(path, array, key='y_mean', dtype=np.float32)
+
+    np.testing.assert_array_equal(load_array(path, key='y_mean'), array)
+    assert _temp_files(str(tmp_path)) == []
+
+
+def test_atomic_save_replaces_an_existing_file(tmp_path):
+    path = str(tmp_path / 'y_mean.mat')
+    save_array_atomic(path, np.zeros((1, 6), dtype=np.float32), key='y_mean')
+
+    replacement = np.arange(6, dtype=np.float32)[np.newaxis]
+    save_array_atomic(path, replacement, key='y_mean')
+
+    np.testing.assert_array_equal(load_array(path, key='y_mean'), replacement)
+    assert _temp_files(str(tmp_path)) == []
+
+
+def test_a_failed_atomic_save_leaves_the_original_intact(tmp_path,
+                                                         monkeypatch):
+    path = str(tmp_path / 'y_mean.mat')
+    original = np.arange(6, dtype=np.float32)[np.newaxis]
+    save_array_atomic(path, original, key='y_mean')
+
+    def broken(*args, **kwargs):
+        raise IOError('disk full')
+
+    monkeypatch.setattr(ridge_factorization, 'save_array', broken)
+    with pytest.raises(IOError):
+        save_array_atomic(path, np.zeros((1, 6), dtype=np.float32),
+                          key='y_mean')
+
+    np.testing.assert_array_equal(load_array(path, key='y_mean'), original)
+    assert _temp_files(str(tmp_path)) == []
+
+
+def test_a_partially_written_temp_file_is_cleaned_up(tmp_path, monkeypatch):
+    """The failure that would otherwise leave rubbish next to the target."""
+    path = str(tmp_path / 'y_mean.mat')
+    real_save = ridge_factorization.save_array
+
+    def write_then_fail(save_file, array, **kwargs):
+        real_save(save_file, array, **kwargs)
+        raise IOError('interrupted after the write')
+
+    monkeypatch.setattr(ridge_factorization, 'save_array', write_then_fail)
+    with pytest.raises(IOError):
+        save_array_atomic(path, np.zeros((1, 6), dtype=np.float32),
+                          key='y_mean')
+
+    assert not os.path.exists(path)
+    assert _temp_files(str(tmp_path)) == []
+
+
+def test_concurrent_atomic_saves_leave_one_valid_file(tmp_path):
+    """Two writers past the existence check, as parallel workers can be.
+
+    Each writer builds its own temp name, so the loser's rename is the only
+    thing that races, and a reader sees one complete file either way.
+    """
+    path = str(tmp_path / 'y_norm.mat')
+    candidates = [np.full((1, 6), value, dtype=np.float32)
+                  for value in (1.0, 2.0)]
+    start = threading.Barrier(len(candidates))
+    errors = []
+
+    def write(array):
+        try:
+            start.wait()
+            save_array_atomic(path, array, key='y_norm')
+        except BaseException as error:  # pragma: no cover - reported below
+            errors.append(error)
+
+    threads = [threading.Thread(target=write, args=(array,))
+               for array in candidates]
+    for thread in threads:
+        thread.start()
+    for thread in threads:
+        thread.join()
+
+    assert errors == []
+    loaded = load_array(path, key='y_norm')
+    assert any(np.array_equal(loaded, array) for array in candidates)
+    assert _temp_files(str(tmp_path)) == []
diff --git a/tests/test_legacy_model_compat.py b/tests/test_legacy_model_compat.py
new file mode 100644
index 0000000..a3c3faf
--- /dev/null
+++ b/tests/test_legacy_model_compat.py
@@ -0,0 +1,196 @@
+"""Reading decoders saved by the direct (pre-factorization) implementation.
+
+New training runs write the factorized format, but large decoders produced by
+the previous implementation may already exist, so ``predict_feature.py`` keeps a
+compatibility path.  The format is detected from the presence of
+``factorized.yaml``.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from ridge_factorization import is_factorized_model_dir
+from tests.helpers import legacy_ridge, pipeline
+
+ALPHA = 100
+CHUNK_AXIS = 1
+
+
+def build_legacy_decoder(dataset, decoder_dir, layers=None,
+                         subjects=('sub-01',), rois=('VC',),
+                         chunk_axis=CHUNK_AXIS):
+    """Train with the reference implementation and save the legacy layout."""
+    trained = {}
+    for layer in (layers or dataset.layers):
+        for subject in subjects:
+            for roi in rois:
+                brain_labels = dataset.labels(subject, 'train')
+                feat_labels = list(np.unique(brain_labels))
+                entry = legacy_ridge.legacy_train(
+                    dataset.brain(subject, roi, 'train'), brain_labels,
+                    dataset.feature_matrix(layer, feat_labels), feat_labels,
+                    alpha=ALPHA, chunk_axis=chunk_axis)
+                legacy_ridge.write_legacy_decoder(
+                    pipeline.layer_model_dir(decoder_dir, layer, subject, roi),
+                    entry)
+                trained[(layer, subject, roi)] = entry
+    return trained
+
+
+def test_format_detection(dataset, tmp_path):
+    factorized_dir = str(tmp_path / 'factorized')
+    pipeline.run_training(dataset, factorized_dir, alpha=ALPHA)
+    assert is_factorized_model_dir(
+        pipeline.brain_model_dir(factorized_dir, 'sub-01', 'VC'))
+
+    legacy_dir = str(tmp_path / 'legacy')
+    build_legacy_decoder(dataset, legacy_dir, layers=['fc_like'])
+    assert not is_factorized_model_dir(
+        pipeline.brain_model_dir(legacy_dir, 'sub-01', 'VC'))
+    assert not is_factorized_model_dir(
+        pipeline.layer_model_dir(legacy_dir, 'fc_like', 'sub-01', 'VC'))
+
+
+@pytest.mark.parametrize('layer', ['fc_like', 'conv_like'])
+def test_legacy_decoder_predictions_are_unchanged(dataset, tmp_path, layer):
+    """A legacy decoder still decodes exactly as it used to, chunked or not."""
+    decoder_dir = str(tmp_path / 'legacy_decoders')
+    decoded_dir = str(tmp_path / 'legacy_decoded')
+
+    trained = build_legacy_decoder(dataset, decoder_dir, layers=[layer])
+
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS, layers=[layer],
+                            subjects=['sub-01'], rois={'VC': 'ROI_VC = 1'})
+
+    test_brain, test_labels = legacy_ridge.average_test_brain(
+        dataset.brain('sub-01', 'VC', 'test'), dataset.labels('sub-01', 'test'))
+    expected = legacy_ridge.legacy_predict(trained[(layer, 'sub-01', 'VC')],
+                                           test_brain)
+    actual = pipeline.read_decoded_features(decoded_dir, layer, 'sub-01', 'VC',
+                                            test_labels)
+
+    np.testing.assert_allclose(actual, expected.astype(np.float32),
+                               rtol=1e-5, atol=1e-6)
+
+
+def test_legacy_decoder_needs_no_training_features(dataset, tmp_path):
+    """The legacy path is self-contained: it must not require ``F`` on disk."""
+    import predict_feature
+
+    decoder_dir = str(tmp_path / 'legacy_decoders')
+    decoded_dir = str(tmp_path / 'legacy_decoded')
+    build_legacy_decoder(dataset, decoder_dir, layers=['fc_like'])
+
+    predict_feature.featdec_predict(
+        {'sub-01': dataset.test_fmri['sub-01']}, decoder_dir,
+        output_dir=decoded_dir,
+        rois={'VC': 'ROI_VC = 1'},
+        label_key=dataset.label_key,
+        layers=['fc_like'],
+        average_sample=True,
+        chunk_axis=CHUNK_AXIS,
+        training_features_paths=None,
+        analysis_name='legacy_without_features')
+
+    names = pipeline.decoded_feature_labels(decoded_dir, 'fc_like', 'sub-01',
+                                            'VC')
+    assert names == dataset.unique_test_labels
+
+
+def test_factorized_decoder_requires_training_features(dataset, tmp_path):
+    """A clear error beats a confusing one -- and it must leave no output.
+
+    The check runs before the output directory is created, so the failed run
+    cannot make the next, correct one print "already done" and skip.
+    """
+    import os
+
+    import predict_feature
+
+    decoder_dir = str(tmp_path / 'decoders')
+    decoded_dir = str(tmp_path / 'decoded')
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    def predict(training_features_paths, analysis_name):
+        return predict_feature.featdec_predict(
+            {'sub-01': dataset.test_fmri['sub-01']}, decoder_dir,
+            output_dir=decoded_dir,
+            rois={'VC': 'ROI_VC = 1'},
+            label_key=dataset.label_key,
+            layers=['fc_like'],
+            average_sample=True,
+            chunk_axis=CHUNK_AXIS,
+            training_features_paths=training_features_paths,
+            analysis_name=analysis_name)
+
+    with pytest.raises(ValueError, match='training_features_paths'):
+        predict(None, 'factorized_without_features')
+
+    assert not os.path.exists(os.path.join(decoded_dir, 'fc_like', 'sub-01',
+                                           'VC'))
+
+    predict([dataset.train_features_dir], 'factorized_with_features')
+
+    assert pipeline.decoded_feature_labels(
+        decoded_dir, 'fc_like', 'sub-01', 'VC') == dataset.unique_test_labels
+
+
+def test_legacy_prediction_still_copies_the_feature_index(dataset, tmp_path):
+    """The pre-factorization contract: given an index, it lands in the output.
+
+    The direct implementation copied ``feature_index_file`` as given, and a
+    legacy decoder is predicted without training features, so the copy cannot
+    depend on them.
+    """
+    import os
+
+    from tests.helpers import synthetic
+
+    decoder_dir = str(tmp_path / 'legacy_decoders')
+    decoded_dir = str(tmp_path / 'legacy_decoded')
+    build_legacy_decoder(dataset, decoder_dir, layers=['fc_like'])
+
+    # Written into the working directory, which is what the old code resolved
+    # against (conftest chdirs into the test's own tmp_path).
+    synthetic.write_feature_index('index_local.mat', {'fc_like': [0, 2]})
+
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS, layers=['fc_like'],
+                            subjects=['sub-01'], rois={'VC': 'ROI_VC = 1'},
+                            features_paths=[],
+                            feature_index_file='index_local.mat',
+                            analysis_name='legacy_with_index')
+
+    assert os.path.isfile(os.path.join(decoded_dir, 'feature_index.mat'))
+
+
+def test_legacy_and_factorized_decoders_agree(dataset, tmp_path):
+    """The compatibility path and the new path give the same numbers."""
+    layer = 'conv_like'
+
+    legacy_dir = str(tmp_path / 'legacy_decoders')
+    legacy_decoded = str(tmp_path / 'legacy_decoded')
+    build_legacy_decoder(dataset, legacy_dir, layers=[layer])
+    pipeline.run_prediction(dataset, legacy_dir, legacy_decoded,
+                            chunk_axis=CHUNK_AXIS, layers=[layer],
+                            subjects=['sub-01'], rois={'VC': 'ROI_VC = 1'},
+                            analysis_name='predict_legacy')
+    legacy_pred = pipeline.read_decoded_features(
+        legacy_decoded, layer, 'sub-01', 'VC', dataset.unique_test_labels)
+
+    new_dir = str(tmp_path / 'new_decoders')
+    new_decoded = str(tmp_path / 'new_decoded')
+    pipeline.run_training(dataset, new_dir, alpha=ALPHA,
+                          subjects=['sub-01'], rois={'VC': 'ROI_VC = 1'},
+                          analysis_name='train_new')
+    pipeline.run_prediction(dataset, new_dir, new_decoded,
+                            chunk_axis=CHUNK_AXIS, layers=[layer],
+                            subjects=['sub-01'], rois={'VC': 'ROI_VC = 1'},
+                            analysis_name='predict_new')
+    new_pred = pipeline.read_decoded_features(
+        new_decoded, layer, 'sub-01', 'VC', dataset.unique_test_labels)
+
+    np.testing.assert_allclose(new_pred, legacy_pred, rtol=1e-4, atol=1e-5)
diff --git a/tests/test_legacy_reference.py b/tests/test_legacy_reference.py
new file mode 100644
index 0000000..7ad10ba
--- /dev/null
+++ b/tests/test_legacy_reference.py
@@ -0,0 +1,97 @@
+"""Pin the reference implementation against the real scripts and the goldens.
+
+``tests/helpers/legacy_ridge.py`` re-implements the direct brain -> feature
+Ridge decoder in plain numpy + scikit-learn.  These tests are what make it
+trustworthy: they show it reproduces what the actual scripts produce, so it can
+later be used as the yardstick for the factorized implementation.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from tests.conftest import GOLDEN_DIR
+from tests.helpers import legacy_ridge, pipeline
+
+import os
+
+ALPHA = 100
+CHUNK_AXIS = 1
+
+
+def reference_prediction(dataset, layer, subject, roi, alpha=ALPHA,
+                         chunk_axis=CHUNK_AXIS):
+    """Train and predict with the reference implementation."""
+    brain = dataset.brain(subject, roi, 'train')
+    brain_labels = dataset.labels(subject, 'train')
+    feat_labels = list(np.unique(brain_labels))
+    feat = dataset.feature_matrix(layer, feat_labels)
+
+    trained = legacy_ridge.legacy_train(
+        brain, brain_labels, feat, feat_labels, alpha=alpha,
+        chunk_axis=chunk_axis)
+
+    test_brain, test_labels = legacy_ridge.average_test_brain(
+        dataset.brain(subject, roi, 'test'), dataset.labels(subject, 'test'))
+    return trained, legacy_ridge.legacy_predict(trained, test_brain), test_labels
+
+
+@pytest.fixture
+def pipeline_outputs(dataset, decoder_dir, decoded_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+    return decoder_dir, decoded_dir
+
+
+def test_reference_matches_scripts(dataset, pipeline_outputs):
+    """The reference reproduces the shipped scripts, chunked and unchunked."""
+    decoder_dir, decoded_dir = pipeline_outputs
+
+    for layer in dataset.layers:
+        for subject in sorted(dataset.train_fmri):
+            for roi in sorted(dataset.rois):
+                trained, expected, test_labels = reference_prediction(
+                    dataset, layer, subject, roi)
+
+                actual = pipeline.read_decoded_features(
+                    decoded_dir, layer, subject, roi, test_labels)
+                np.testing.assert_allclose(
+                    actual, expected.astype(np.float32),
+                    rtol=1e-5, atol=1e-6,
+                    err_msg='prediction mismatch for %s/%s/%s'
+                            % (layer, subject, roi))
+
+                saved = pipeline.read_norm_params(decoder_dir, subject, roi)
+                saved.update(pipeline.read_feature_statistics(
+                    decoder_dir, layer, subject, roi))
+                for key in legacy_ridge.NORM_KEYS:
+                    np.testing.assert_allclose(
+                        saved[key], trained[key].astype(np.float32),
+                        rtol=1e-6, atol=1e-7,
+                        err_msg='%s mismatch for %s/%s/%s'
+                                % (key, layer, subject, roi))
+
+
+def test_reference_matches_golden(dataset):
+    """The reference reproduces the committed pre-factorization fixtures."""
+    golden = np.load(os.path.join(GOLDEN_DIR, 'sklearn_ridge_pipeline.npz'))
+    assert int(golden['_alpha']) == ALPHA
+    assert int(golden['_chunk_axis']) == CHUNK_AXIS
+
+    for layer in [str(x) for x in golden['_layers']]:
+        for subject in [str(x) for x in golden['_subjects']]:
+            for roi in [str(x) for x in golden['_rois']]:
+                trained, expected, test_labels = reference_prediction(
+                    dataset, layer, subject, roi)
+                assert test_labels == [str(x) for x in golden['_test_labels']]
+
+                prefix = '%s|%s|%s|' % (layer, subject, roi)
+                np.testing.assert_allclose(
+                    golden[prefix + 'pred'], expected.astype(np.float32),
+                    rtol=1e-5, atol=1e-6)
+                for key in legacy_ridge.NORM_KEYS:
+                    np.testing.assert_allclose(
+                        golden[prefix + key], trained[key].astype(np.float32),
+                        rtol=1e-6, atol=1e-7)
diff --git a/tests/test_normalization.py b/tests/test_normalization.py
new file mode 100644
index 0000000..c73229e
--- /dev/null
+++ b/tests/test_normalization.py
@@ -0,0 +1,122 @@
+"""Normalization semantics of the factorized decoder.
+
+Brain normalization is unchanged.  Feature normalization is where the
+factorization removes work: the round trip ``(Y - y_mean) / y_norm`` at training
+followed by ``* y_norm + y_mean`` at prediction cancels exactly (see
+``ridge_factorization``), so the factorized prediction combines the raw training
+features.  ``y_mean.mat`` / ``y_norm.mat`` still exist in the decoder's
+per-layer directory because ``evaluation.py`` normalizes the pattern
+correlation with them, but nothing in the decoding needs them, so
+``predict_feature.py`` writes them as compatibility sidecars (see
+``tests/test_lazy_feature_statistics.py``).
+"""
+
+from __future__ import annotations
+
+import os
+
+import numpy as np
+
+from bdpy.dataform import load_array, save_array
+from ridge_factorization import normalize_brain_for_training
+from tests.helpers import pipeline
+
+ALPHA = 100
+CHUNK_AXIS = 1
+
+
+def test_brain_normalization_zeroes_infinities():
+    """Matches ``ModelTraining.run()``: infinities become zero, NaNs do not."""
+    brain = np.array([[1.0, 2.0, 3.0],
+                      [1.0, 4.0, 3.0]])
+    mean = np.array([[0.0, 3.0, 3.0]])
+    std = np.array([[0.0, 1.0, 0.0]])
+
+    normalized = normalize_brain_for_training(brain, mean, std)
+
+    # Column 0: nonzero deviation over zero SD -> inf -> replaced by 0.
+    np.testing.assert_array_equal(normalized[:, 0], [0.0, 0.0])
+    # Column 1: ordinary z-score.
+    np.testing.assert_allclose(normalized[:, 1], [-1.0, 1.0])
+    # Column 2: 0/0 -> NaN, deliberately left as is.
+    assert np.all(np.isnan(normalized[:, 2]))
+
+
+def test_feature_statistics_are_materialized_by_prediction(
+        dataset, decoder_dir, decoded_dir):
+    """Training cannot compute them -- it never reads the features.
+
+    ``evaluation.py`` reads them from ``<layer>/<subject>/<roi>/model/``, so
+    the definition has to stay exactly what the direct implementation wrote:
+    per-unit statistics over the unique stimuli, with ``ddof=1``.
+    """
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    for layer in dataset.layers:
+        for subject in dataset.train_fmri:
+            for roi in dataset.rois:
+                assert not os.path.exists(os.path.join(
+                    pipeline.layer_model_dir(decoder_dir, layer, subject, roi),
+                    'y_mean.mat'))
+
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    for layer in dataset.layers:
+        for subject in dataset.test_fmri:
+            feat = dataset.feature_matrix(
+                layer, list(np.unique(dataset.labels(subject, 'train'))))
+            for roi in dataset.rois:
+                params = pipeline.read_feature_statistics(decoder_dir, layer,
+                                                          subject, roi)
+                np.testing.assert_allclose(
+                    params['y_mean'],
+                    np.mean(feat, axis=0)[np.newaxis].astype(np.float32),
+                    rtol=1e-5, atol=1e-6)
+                np.testing.assert_allclose(
+                    params['y_norm'],
+                    np.std(feat, axis=0, ddof=1)[np.newaxis].astype(np.float32),
+                    rtol=1e-5, atol=1e-6)
+
+
+def test_prediction_does_not_depend_on_the_feature_statistics(dataset,
+                                                              tmp_path):
+    """Corrupting ``y_mean``/``y_norm`` must not change what is decoded.
+
+    A direct demonstration that the feature normalization cancels: the
+    factorized prediction never reads these files.  They are corrupted in the
+    decoder and the features are decoded again into a fresh directory, so the
+    second run really does see them.
+    """
+    decoder_dir = str(tmp_path / 'decoders')
+
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, str(tmp_path / 'intact'),
+                            chunk_axis=CHUNK_AXIS,
+                            analysis_name='predict_intact')
+
+    intact = {layer: pipeline.read_decoded_features(
+        str(tmp_path / 'intact'), layer, 'sub-01', 'VC',
+        dataset.unique_test_labels) for layer in dataset.layers}
+
+    for layer in dataset.layers:
+        for subject in dataset.test_fmri:
+            for roi in dataset.rois:
+                directory = pipeline.layer_model_dir(decoder_dir, layer,
+                                                     subject, roi)
+                for key in ('y_mean', 'y_norm'):
+                    path = os.path.join(directory, '%s.mat' % key)
+                    array = load_array(path, key=key)
+                    save_array(path, array * 3.0 + 17.0, key=key,
+                               dtype=np.float32, sparse=False)
+
+    pipeline.run_prediction(dataset, decoder_dir, str(tmp_path / 'corrupted'),
+                            chunk_axis=CHUNK_AXIS,
+                            analysis_name='predict_corrupted')
+
+    for layer in dataset.layers:
+        np.testing.assert_array_equal(
+            pipeline.read_decoded_features(str(tmp_path / 'corrupted'), layer,
+                                           'sub-01', 'VC',
+                                           dataset.unique_test_labels),
+            intact[layer])
diff --git a/tests/test_pipeline_golden.py b/tests/test_pipeline_golden.py
new file mode 100644
index 0000000..1b94e09
--- /dev/null
+++ b/tests/test_pipeline_golden.py
@@ -0,0 +1,64 @@
+"""End-to-end regression test against the committed golden fixtures.
+
+The fixtures in ``tests/data/golden/`` were produced by the direct
+(pre-factorization) sklearn Ridge implementation; regenerate them only with
+``tests/generate_golden.py`` and only when the numerical output is *meant* to
+change.
+"""
+
+from __future__ import annotations
+
+import os
+
+import numpy as np
+import pytest
+
+from tests.conftest import GOLDEN_DIR
+from tests.helpers import legacy_ridge, pipeline
+
+GOLDEN_FILE = os.path.join(GOLDEN_DIR, 'sklearn_ridge_pipeline.npz')
+
+
+@pytest.fixture(scope='module')
+def golden():
+    return np.load(GOLDEN_FILE)
+
+
+def test_golden_fixture_is_present(golden):
+    assert len(golden.files) > 0
+
+
+def test_pipeline_reproduces_golden(dataset, decoder_dir, decoded_dir, golden):
+    alpha = int(golden['_alpha'])
+    chunk_axis = int(golden['_chunk_axis'])
+
+    pipeline.run_training(dataset, decoder_dir, alpha=alpha)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=chunk_axis)
+
+    test_labels = [str(x) for x in golden['_test_labels']]
+
+    for layer in [str(x) for x in golden['_layers']]:
+        for subject in [str(x) for x in golden['_subjects']]:
+            for roi in [str(x) for x in golden['_rois']]:
+                prefix = '%s|%s|%s|' % (layer, subject, roi)
+
+                actual = pipeline.read_decoded_features(
+                    decoded_dir, layer, subject, roi, test_labels)
+                np.testing.assert_allclose(
+                    actual, golden[prefix + 'pred'], rtol=1e-4, atol=1e-5,
+                    err_msg='decoded features drifted for %s/%s/%s'
+                            % (layer, subject, roi))
+
+                saved = pipeline.read_norm_params(decoder_dir, subject, roi)
+                # The feature statistics are written by prediction now, but
+                # into the same place and with the same values the direct
+                # implementation produced.
+                saved.update(pipeline.read_feature_statistics(
+                    decoder_dir, layer, subject, roi))
+                for key in legacy_ridge.NORM_KEYS:
+                    np.testing.assert_allclose(
+                        saved[key], golden[prefix + key],
+                        rtol=1e-6, atol=1e-7,
+                        err_msg='%s drifted for %s/%s/%s'
+                                % (key, layer, subject, roi))
diff --git a/tests/test_predict_feature.py b/tests/test_predict_feature.py
new file mode 100644
index 0000000..40094f9
--- /dev/null
+++ b/tests/test_predict_feature.py
@@ -0,0 +1,129 @@
+"""Behavior of ``predict_feature.featdec_predict``.
+
+Pins the decoded-feature output contract: one ``<label>.mat`` per predicted
+sample under ``<output>/<layer>/<subject>/<roi>/``, key ``feat``, shape
+``(1, *feature_shape)``, ``float32`` -- the layout
+``bdpy.dataform.DecodedFeatures`` (and therefore ``evaluation.py``) expects.
+"""
+
+from __future__ import annotations
+
+import os
+
+import numpy as np
+
+from tests.helpers import legacy_ridge, pipeline
+
+ALPHA = 100
+CHUNK_AXIS = 1
+
+
+def test_output_layout_and_dtype(dataset, decoder_dir, decoded_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    from bdpy.dataform import load_array
+
+    for layer in dataset.layers:
+        feature_shape = dataset.layer_shapes[layer]
+        for subject in dataset.test_fmri:
+            for roi in dataset.rois:
+                directory = os.path.join(decoded_dir, layer, subject, roi)
+                names = pipeline.decoded_feature_labels(decoded_dir, layer,
+                                                        subject, roi)
+                assert names == dataset.unique_test_labels
+                for label in names:
+                    array = load_array(
+                        os.path.join(directory, '%s.mat' % label), key='feat')
+                    assert array.shape == (1,) + feature_shape
+                    assert array.dtype == np.float32
+
+
+def test_average_sample_averages_repetitions(dataset, decoder_dir,
+                                             decoded_dir):
+    """With ``average_sample`` on, one file per unique stimulus is written."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS, average_sample=True)
+
+    names = pipeline.decoded_feature_labels(decoded_dir, dataset.layers[0],
+                                            'sub-01', 'VC')
+    assert names == dataset.unique_test_labels
+    assert len(names) < len(dataset.labels('sub-01', 'test'))
+
+
+def test_single_trial_output_names(dataset, decoder_dir, decoded_dir):
+    """With ``average_sample`` off, files are ``sample%06d-<label>.mat``."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS, average_sample=False)
+
+    trial_labels = dataset.labels('sub-01', 'test')
+    expected = sorted(legacy_ridge.single_trial_labels(trial_labels))
+    names = pipeline.decoded_feature_labels(decoded_dir, dataset.layers[0],
+                                            'sub-01', 'VC')
+    assert names == expected
+
+
+def test_excluded_labels_are_dropped_when_averaging(dataset, decoder_dir,
+                                                    decoded_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    excluded = [dataset.unique_test_labels[0]]
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS, average_sample=True,
+                            excluded_labels=excluded)
+
+    names = pipeline.decoded_feature_labels(decoded_dir, dataset.layers[0],
+                                            'sub-01', 'VC')
+    assert names == dataset.unique_test_labels[1:]
+
+
+def test_existing_output_directory_is_skipped(dataset, decoder_dir,
+                                              decoded_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    layer, subject, roi = dataset.layers[0], 'sub-01', 'VC'
+    target = os.path.join(decoded_dir, layer, subject, roi,
+                          '%s.mat' % dataset.unique_test_labels[0])
+    mtime = os.stat(target).st_mtime_ns
+
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS,
+                            analysis_name='test_prediction_rerun')
+
+    assert os.stat(target).st_mtime_ns == mtime
+
+
+def test_prediction_uses_the_decoders_brain_normalization(dataset, decoder_dir,
+                                                          decoded_dir):
+    """Test brain data is z-scored with the *training* mean/SD, not its own."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    layer, subject, roi = 'fc_like', 'sub-01', 'VC'
+    brain_labels = dataset.labels(subject, 'train')
+    feat_labels = list(np.unique(brain_labels))
+    trained = legacy_ridge.legacy_train(
+        dataset.brain(subject, roi, 'train'), brain_labels,
+        dataset.feature_matrix(layer, feat_labels), feat_labels,
+        alpha=ALPHA, chunk_axis=CHUNK_AXIS)
+
+    test_brain, test_labels = legacy_ridge.average_test_brain(
+        dataset.brain(subject, roi, 'test'), dataset.labels(subject, 'test'))
+    expected = legacy_ridge.legacy_predict(trained, test_brain)
+    actual = pipeline.read_decoded_features(decoded_dir, layer, subject, roi,
+                                            test_labels)
+    np.testing.assert_allclose(actual, expected.astype(np.float32),
+                               rtol=1e-5, atol=1e-6)
+
+    # Z-scoring the test data with its own statistics would give something else.
+    self_scaled = ((test_brain - np.mean(test_brain, axis=0)[np.newaxis, :])
+                   / np.std(test_brain, axis=0, ddof=1)[np.newaxis, :])
+    other = legacy_ridge.legacy_predict(
+        {**trained, 'x_mean': np.zeros_like(trained['x_mean']),
+         'x_norm': np.ones_like(trained['x_norm'])}, self_scaled)
+    assert not np.allclose(actual, other, rtol=1e-3)
diff --git a/tests/test_ridge_factorization.py b/tests/test_ridge_factorization.py
new file mode 100644
index 0000000..8c44616
--- /dev/null
+++ b/tests/test_ridge_factorization.py
@@ -0,0 +1,235 @@
+"""Exact equivalence of the direct and factorized Ridge formulations.
+
+These tests work in ``float64`` on small arrays, so the tolerances can be tight
+enough to be a real statement about the algebra rather than about float32
+rounding.  See ``ridge_factorization`` for the derivation; the identities
+checked here are the three steps it rests on:
+
+1. ``Ridge(alpha).fit(X, R @ F).predict(Xt) == Ridge(alpha).fit(X, R).predict(Xt) @ F``
+2. the rows of the predicted stimulus coefficients sum to one, and
+3. therefore the feature mean/SD normalization cancels.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+from sklearn.linear_model import Ridge
+
+from ridge_factorization import combine_features, one_hot_assignment
+
+RTOL = 1e-9
+ATOL = 1e-9
+
+
+def make_problem(n_stimuli=6, repeats=(3, 3, 2, 2, 1, 1), n_voxels=5,
+                 n_features=7, n_test=4, seed=0):
+    """Random brain/feature problem with repeated stimulus presentations."""
+    assert len(repeats) == n_stimuli
+    rng = np.random.RandomState(seed)
+
+    labels = []
+    for stimulus, count in enumerate(repeats):
+        labels.extend([stimulus] * count)
+    unique_labels = list(range(n_stimuli))
+
+    assignment = one_hot_assignment(labels, unique_labels)
+    features = rng.randn(n_stimuli, n_features)
+    brain = rng.randn(len(labels), n_voxels)
+    brain_test = rng.randn(n_test, n_voxels)
+    return brain, brain_test, assignment, features
+
+
+def direct_prediction(brain, brain_test, assignment, features, alpha):
+    """Fit on the expanded target ``Y = R F`` -- the original formulation."""
+    model = Ridge(alpha=alpha)
+    model.fit(brain, assignment @ features)
+    return model.predict(brain_test)
+
+
+def factorized_prediction(brain, brain_test, assignment, features, alpha,
+                          chunk_axis=None):
+    """Fit on the stimulus basis ``R``, then combine with ``F``."""
+    model = Ridge(alpha=alpha)
+    model.fit(brain, assignment)
+    coefficients = model.predict(brain_test)
+    return combine_features(coefficients, features, chunk_axis=chunk_axis)
+
+
+@pytest.mark.parametrize('alpha', [0.0, 1e-3, 1.0, 100.0, 1e6])
+def test_equivalent_for_a_range_of_alpha(alpha):
+    brain, brain_test, assignment, features = make_problem()
+
+    np.testing.assert_allclose(
+        factorized_prediction(brain, brain_test, assignment, features, alpha),
+        direct_prediction(brain, brain_test, assignment, features, alpha),
+        rtol=RTOL, atol=ATOL)
+
+
+@pytest.mark.parametrize('repeats', [
+    (5, 5, 5, 5, 5, 5),      # balanced repetitions
+    (3, 3, 2, 2, 1, 1),      # unbalanced repetitions
+    (1, 1, 1, 1, 1, 1),      # no repetition at all (M == N)
+    (10, 1, 1, 1, 1, 1),     # strongly unbalanced
+])
+def test_equivalent_for_any_repetition_structure(repeats):
+    brain, brain_test, assignment, features = make_problem(repeats=repeats)
+
+    np.testing.assert_allclose(
+        factorized_prediction(brain, brain_test, assignment, features, 100.0),
+        direct_prediction(brain, brain_test, assignment, features, 100.0),
+        rtol=RTOL, atol=ATOL)
+
+
+@pytest.mark.parametrize('n_features', [1, 2, 50])
+def test_equivalent_for_single_and_multi_output(n_features):
+    brain, brain_test, assignment, features = make_problem(
+        n_features=n_features)
+
+    actual = factorized_prediction(brain, brain_test, assignment, features,
+                                   100.0)
+    expected = direct_prediction(brain, brain_test, assignment, features,
+                                 100.0)
+
+    # scikit-learn ravels a single-column target, so for n_features == 1 the
+    # direct formulation returns shape (n_test,) where the factorized one keeps
+    # (n_test, 1).  The values are what the equivalence is about; the real
+    # pipeline always has n_features >= 2 anyway.
+    np.testing.assert_allclose(actual, expected.reshape(actual.shape),
+                               rtol=RTOL, atol=ATOL)
+
+
+@pytest.mark.parametrize('n_voxels,n_stimuli,repeats', [
+    (3, 6, (3, 3, 2, 2, 1, 1)),     # fewer voxels than stimuli
+    (40, 6, (3, 3, 2, 2, 1, 1)),    # more voxels than trials (dual solver)
+])
+def test_equivalent_regardless_of_problem_shape(n_voxels, n_stimuli, repeats):
+    brain, brain_test, assignment, features = make_problem(
+        n_stimuli=n_stimuli, repeats=repeats, n_voxels=n_voxels)
+
+    np.testing.assert_allclose(
+        factorized_prediction(brain, brain_test, assignment, features, 100.0),
+        direct_prediction(brain, brain_test, assignment, features, 100.0),
+        rtol=RTOL, atol=ATOL)
+
+
+def test_equivalent_for_multidimensional_features():
+    brain, brain_test, assignment, _ = make_problem()
+    rng = np.random.RandomState(3)
+    features = rng.randn(assignment.shape[1], 3, 2, 2)
+
+    flat = features.reshape(features.shape[0], -1)
+    expected = direct_prediction(brain, brain_test, assignment, flat, 100.0)
+    expected = expected.reshape((brain_test.shape[0],) + features.shape[1:])
+
+    for chunk_axis in (None, 1, 2, 3):
+        actual = factorized_prediction(brain, brain_test, assignment, features,
+                                       100.0, chunk_axis=chunk_axis)
+        assert actual.shape == expected.shape
+        np.testing.assert_allclose(actual, expected, rtol=RTOL, atol=ATOL)
+
+
+def test_predicted_stimulus_coefficients_sum_to_one():
+    """``C @ 1 == 1``: the identity that makes the feature mean cancel."""
+    brain, brain_test, assignment, _ = make_problem()
+
+    for alpha in (1e-3, 1.0, 100.0, 1e6):
+        model = Ridge(alpha=alpha).fit(brain, assignment)
+        coefficients = model.predict(brain_test)
+        np.testing.assert_allclose(coefficients.sum(axis=1),
+                                   np.ones(brain_test.shape[0]),
+                                   rtol=RTOL, atol=ATOL)
+
+
+def test_feature_normalization_cancels():
+    """Normalizing the target and un-normalizing the prediction is a no-op.
+
+    The direct pipeline fitted on ``(R F - m) / s`` and un-normalized with
+    ``* s + m``.  Because the stimulus coefficients sum to one, that round trip
+    returns exactly the un-normalized linear combination.
+    """
+    brain, brain_test, assignment, features = make_problem()
+
+    mean = np.mean(features, axis=0)[np.newaxis, :]
+    std = np.std(features, axis=0, ddof=1)[np.newaxis, :]
+    assert np.all(std > 0)
+
+    normalized = direct_prediction(brain, brain_test, assignment,
+                                   (features - mean) / std, 100.0)
+    round_tripped = normalized * std + mean
+
+    np.testing.assert_allclose(
+        factorized_prediction(brain, brain_test, assignment, features, 100.0),
+        round_tripped, rtol=RTOL, atol=ATOL)
+
+
+def test_intercept_is_still_fitted_and_applied():
+    """The factorization keeps sklearn's intercept, it does not drop it."""
+    brain, brain_test, assignment, features = make_problem()
+
+    model = Ridge(alpha=100.0).fit(brain, assignment)
+    assert np.any(np.abs(model.intercept_) > 1e-6)
+
+    with_intercept = model.predict(brain_test)
+    without_intercept = brain_test @ model.coef_.T
+    assert not np.allclose(with_intercept, without_intercept)
+
+
+def test_intercept_is_what_makes_the_feature_normalization_cancel():
+    """Locate precisely what ``fit_intercept=True`` (the default) buys.
+
+    The factorization itself does not need the intercept: Ridge prediction is
+    linear in the target either way, so ``S(R F) == (S R) F`` holds with
+    ``fit_intercept=False`` too.  What the intercept buys is ``C @ 1 == 1``, and
+    *that* is what makes the feature mean/SD round trip cancel.  Since the
+    implementation drops that round trip, it does depend on the default.
+    """
+    brain, brain_test, assignment, features = make_problem()
+
+    direct = Ridge(alpha=100.0, fit_intercept=False)
+    direct.fit(brain, assignment @ features)
+    basis = Ridge(alpha=100.0, fit_intercept=False)
+    basis.fit(brain, assignment)
+    coefficients = basis.predict(brain_test)
+
+    # The two formulations still agree without an intercept ...
+    np.testing.assert_allclose(combine_features(coefficients, features),
+                               direct.predict(brain_test),
+                               rtol=RTOL, atol=ATOL)
+
+    # ... but the rows no longer sum to one, so the un-normalization would not
+    # cancel any more.
+    assert not np.allclose(coefficients.sum(axis=1), 1.0)
+
+    mean = np.mean(features, axis=0)[np.newaxis, :]
+    std = np.std(features, axis=0, ddof=1)[np.newaxis, :]
+    normalized = Ridge(alpha=100.0, fit_intercept=False).fit(
+        brain, assignment @ ((features - mean) / std)).predict(brain_test)
+    assert not np.allclose(normalized * std + mean,
+                           combine_features(coefficients, features),
+                           rtol=1e-6)
+
+
+def test_one_hot_assignment_matches_the_expansion():
+    labels = ['b', 'a', 'b', 'c', 'a', 'b']
+    unique_labels = sorted(set(labels))
+
+    assignment = one_hot_assignment(labels, unique_labels)
+
+    assert assignment.shape == (len(labels), len(unique_labels))
+    np.testing.assert_array_equal(assignment.sum(axis=1), np.ones(len(labels)))
+
+    rng = np.random.RandomState(0)
+    features = rng.randn(len(unique_labels), 4)
+    expected = features[[unique_labels.index(lb) for lb in labels]]
+    np.testing.assert_array_equal(assignment @ features, expected)
+
+
+def test_one_hot_assignment_rejects_unknown_labels():
+    with pytest.raises(ValueError, match='not in unique_labels'):
+        one_hot_assignment(['a', 'z'], ['a', 'b'])
+
+
+def test_one_hot_assignment_rejects_duplicate_basis_labels():
+    with pytest.raises(ValueError, match='duplicates'):
+        one_hot_assignment(['a'], ['a', 'a'])
diff --git a/tests/test_serialization.py b/tests/test_serialization.py
new file mode 100644
index 0000000..d2aa613
--- /dev/null
+++ b/tests/test_serialization.py
@@ -0,0 +1,78 @@
+"""Round trip: train -> serialize -> load -> predict."""
+
+from __future__ import annotations
+
+import os
+import pickle
+
+import numpy as np
+
+import ridge_factorization
+from ridge_factorization import combine_features, load_train_labels
+from tests.helpers import legacy_ridge, pipeline
+
+ALPHA = 100
+CHUNK_AXIS = 1
+
+
+def test_reloaded_model_reproduces_the_saved_prediction(dataset, decoder_dir,
+                                                        decoded_dir):
+    """Recomputing the decode by hand from the artifact gives the same output."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    for layer in dataset.layers:
+        subject, roi = 'sub-01', 'VC'
+        directory = pipeline.brain_model_dir(decoder_dir, subject, roi)
+
+        with open(ridge_factorization.model_file_path(directory), 'rb') as f:
+            model = pickle.load(f)['model']
+        train_labels = load_train_labels(directory)
+        params = pipeline.read_norm_params(decoder_dir, subject, roi,
+                                           keys=pipeline.BRAIN_NORM_KEYS)
+
+        test_brain, test_labels = legacy_ridge.average_test_brain(
+            dataset.brain(subject, roi, 'test'),
+            dataset.labels(subject, 'test'))
+        normalized = ((test_brain - params['x_mean']) / params['x_norm'])
+        coefficients = model.predict(normalized.astype(np.float32))
+        features = dataset.feature_matrix(layer, train_labels)
+        expected = combine_features(coefficients, features.astype(np.float32),
+                                    chunk_axis=CHUNK_AXIS)
+
+        actual = pipeline.read_decoded_features(decoded_dir, layer, subject,
+                                                roi, test_labels)
+        np.testing.assert_allclose(actual, expected.astype(np.float32),
+                                   rtol=1e-6, atol=1e-7)
+
+
+def test_prediction_is_reproducible_across_runs(dataset, decoder_dir,
+                                                tmp_path):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    outputs = []
+    for i in range(2):
+        decoded_dir = str(tmp_path / ('decoded_%d' % i))
+        pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                                chunk_axis=CHUNK_AXIS, layers=['fc_like'],
+                                analysis_name='predict_%d' % i)
+        outputs.append(pipeline.read_decoded_features(
+            decoded_dir, 'fc_like', 'sub-01', 'VC',
+            dataset.unique_test_labels))
+
+    np.testing.assert_array_equal(outputs[0], outputs[1])
+
+
+def test_artifact_file_set(dataset, decoder_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    assert sorted(os.listdir(directory)) == [
+        'factorized.yaml',
+        'info.yaml',
+        'model.pkl.gz',
+        'train_labels.json',
+        'x_mean.mat',
+        'x_norm.mat',
+    ]
diff --git a/tests/test_stimulus_correspondence.py b/tests/test_stimulus_correspondence.py
new file mode 100644
index 0000000..f873e50
--- /dev/null
+++ b/tests/test_stimulus_correspondence.py
@@ -0,0 +1,306 @@
+"""The decoder's only requirement about the features it is combined with.
+
+The factorized fit is on the *stimulus basis*: it maps brain activity to one
+coefficient per training stimulus and never sees a feature value.  What the
+stored coefficients therefore require of the features they are combined with is
+exactly one thing:
+
+    coefficient column i is multiplied by the feature row of training stimulus i
+
+These tests cover the three places that can break it -- the stored label order,
+the row assembly, and the width of the model -- and that the decoder's own
+training features go through them untouched.
+
+(The algebra also holds for any other feature matrix with the same row order;
+that is a property of `C @ F`, tested in ``tests/test_ridge_factorization.py``,
+not a supported way to use a decoder *directory*, which carries one set of
+feature statistics for ``evaluation.py``.)
+"""
+
+from __future__ import annotations
+
+import json
+import os
+import pickle
+
+import numpy as np
+import pytest
+import yaml
+
+import ridge_factorization
+from ridge_factorization import (
+    TrainingFeatureLoader,
+    labels_digest,
+    resolve_feature_index_path,
+)
+from tests.helpers import legacy_ridge, pipeline, synthetic
+
+ALPHA = 100
+CHUNK_AXIS = 1
+ONE_SUBJECT = ['sub-01']
+ONE_ROI = {'VC': 'ROI_VC = 1'}
+LAYER = 'fc_like'
+
+
+def _train(dataset, decoder_dir, **kwargs):
+    return pipeline.run_training(
+        dataset, decoder_dir, alpha=ALPHA, subjects=ONE_SUBJECT, rois=ONE_ROI, **kwargs)
+
+
+def _predict(dataset, decoder_dir, decoded_dir, **kwargs):
+    return pipeline.run_prediction(
+        dataset, decoder_dir, decoded_dir, chunk_axis=CHUNK_AXIS,
+        layers=[LAYER], subjects=ONE_SUBJECT, rois=ONE_ROI, **kwargs)
+
+
+def _stimulus_coefficients(dataset, decoder_dir):
+    """Run the decoder's brain-side model by hand, on averaged test trials."""
+    directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    with open(ridge_factorization.model_file_path(directory), 'rb') as f:
+        model = pickle.load(f)['model']
+    params = pipeline.read_norm_params(decoder_dir, 'sub-01', 'VC')
+    test_brain, _ = legacy_ridge.average_test_brain(
+        dataset.brain('sub-01', 'VC', 'test'), dataset.labels('sub-01', 'test'))
+    normalized = (test_brain - params['x_mean']) / params['x_norm']
+    return model.predict(normalized.astype(np.float32))
+
+
+# The invariant, positively #####################################################
+
+def test_matching_features_are_accepted(dataset, decoder_dir, decoded_dir):
+    """The checks must not reject the decoder's own training features."""
+    _train(dataset, decoder_dir)
+    _predict(dataset, decoder_dir, decoded_dir)
+
+    brain_labels = dataset.labels('sub-01', 'train')
+    feat_labels = list(np.unique(brain_labels))
+    trained = legacy_ridge.legacy_train(
+        dataset.brain('sub-01', 'VC', 'train'), brain_labels,
+        dataset.feature_matrix(LAYER, feat_labels), feat_labels,
+        alpha=ALPHA, chunk_axis=CHUNK_AXIS)
+    test_brain, test_labels = legacy_ridge.average_test_brain(
+        dataset.brain('sub-01', 'VC', 'test'), dataset.labels('sub-01', 'test'))
+
+    np.testing.assert_allclose(
+        pipeline.read_decoded_features(decoded_dir, LAYER, 'sub-01', 'VC',
+                                       test_labels),
+        legacy_ridge.legacy_predict(trained, test_brain).astype(np.float32),
+        rtol=1e-5, atol=1e-6)
+
+
+# The stored label order ########################################################
+
+def test_permuted_train_labels_are_rejected(dataset, tmp_path):
+    """The case that rules out per-column statistics as the guard.
+
+    Reversing ``train_labels.json`` leaves every feature mean and SD untouched
+    while pairing each coefficient column with the wrong stimulus, so only a
+    digest of the *ordered* labels can see it -- and it needs no feature read.
+    """
+    decoder_dir = str(tmp_path / 'decoders')
+    _train(dataset, decoder_dir)
+    model_dir = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    labels_path = os.path.join(model_dir, 'train_labels.json')
+    with open(labels_path) as f:
+        labels = json.load(f)
+    with open(labels_path, 'w') as f:
+        json.dump(list(reversed(labels)), f)
+
+    ordered = dataset.feature_matrix(LAYER, labels)
+    permuted = dataset.feature_matrix(LAYER, list(reversed(labels)))
+    # (Only the float32 summation order differs, hence a tolerance rather than
+    # equality.)
+    np.testing.assert_allclose(np.mean(permuted, axis=0),
+                               np.mean(ordered, axis=0), rtol=1e-6)
+    np.testing.assert_allclose(np.std(permuted, axis=0, ddof=1),
+                               np.std(ordered, axis=0, ddof=1), rtol=1e-6)
+
+    with pytest.raises(ValueError, match='do not match the order'):
+        _predict(dataset, decoder_dir, str(tmp_path / 'decoded'),
+                 analysis_name='permuted')
+
+
+def test_a_renamed_train_label_is_rejected(dataset, tmp_path):
+    """Editing one label keeps the count but changes which stimulus is meant."""
+    decoder_dir = str(tmp_path / 'decoders')
+    _train(dataset, decoder_dir)
+    labels_path = os.path.join(
+        pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC'),
+        'train_labels.json')
+    with open(labels_path) as f:
+        labels = json.load(f)
+    labels[0] = labels[-1]
+    with open(labels_path, 'w') as f:
+        json.dump(labels, f)
+
+    with pytest.raises(ValueError, match='do not match the order'):
+        _predict(dataset, decoder_dir, str(tmp_path / 'decoded'),
+                 analysis_name='renamed')
+
+
+def test_a_model_of_the_wrong_width_is_rejected(dataset, tmp_path):
+    """The last way the three can disagree, once the digests agree.
+
+    The label list, its digest and ``n_train_stimuli`` are all made consistent
+    here, so what fires is the check at the contraction rather than an earlier
+    one -- and its message names the invariant, not a digest.
+    """
+    decoder_dir = str(tmp_path / 'decoders')
+    _train(dataset, decoder_dir)
+    model_dir = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+
+    labels_path = os.path.join(model_dir, 'train_labels.json')
+    with open(labels_path) as f:
+        labels = json.load(f)
+    truncated = labels[:-1]
+    with open(labels_path, 'w') as f:
+        json.dump(truncated, f)
+
+    meta_path = os.path.join(model_dir, 'factorized.yaml')
+    with open(meta_path) as f:
+        meta = yaml.safe_load(f)
+    meta['labels_digest'] = labels_digest(truncated)
+    meta['n_train_stimuli'] = len(truncated)
+    with open(meta_path, 'w') as f:
+        yaml.dump(meta, f)
+
+    with pytest.raises(ValueError, match='column i must be combined'):
+        _predict(dataset, decoder_dir, str(tmp_path / 'decoded'),
+                 analysis_name='narrow')
+
+
+# The row assembly ##############################################################
+
+def test_a_missing_training_stimulus_is_rejected(dataset, tmp_path):
+    decoder_dir = str(tmp_path / 'decoders')
+    _train(dataset, decoder_dir)
+
+    partial = str(tmp_path / 'partial_features')
+    synthetic.write_feature_tree(partial, {LAYER: dataset.layer_shapes[LAYER]},
+                                 dataset.unique_train_labels[:-1], seed=3)
+
+    with pytest.raises(ValueError, match='exactly one feature store'):
+        _predict(dataset, decoder_dir, str(tmp_path / 'decoded'),
+                 features_paths=[partial], analysis_name='missing')
+
+
+def test_a_duplicated_stimulus_cannot_hide_behind_a_missing_one(dataset,
+                                                                tmp_path):
+    """Why the row count alone is not a check.
+
+    Two stores: the first holds every stimulus but the last, the second holds
+    only the first stimulus.  ``bdpy``'s ``get_multi_features`` emits one row
+    for the duplicated label from each store and none for the missing one,
+    landing on exactly ``len(labels)`` rows with every row after the duplicate
+    shifted against its coefficient column.
+    """
+    from bdpy.dataform import Features
+    from bdpy.dataform.utils import get_multi_features
+
+    decoder_dir = str(tmp_path / 'decoders')
+    _train(dataset, decoder_dir)
+
+    labels = dataset.unique_train_labels
+    shapes = {LAYER: dataset.layer_shapes[LAYER]}
+    first = str(tmp_path / 'store_a')
+    second = str(tmp_path / 'store_b')
+    synthetic.write_feature_tree(first, shapes, labels[:-1], seed=4)
+    synthetic.write_feature_tree(second, shapes, labels[:1], seed=5)
+
+    stores = [Features(first), Features(second)]
+    assert get_multi_features(stores, LAYER, labels=labels).shape[0] == \
+        len(labels), 'the bdpy behavior this test is about has changed'
+
+    loader = TrainingFeatureLoader([first, second])
+    with pytest.raises(ValueError, match='exactly one feature store'):
+        loader.get(LAYER, labels)
+
+    with pytest.raises(ValueError, match='exactly one feature store'):
+        _predict(dataset, decoder_dir, str(tmp_path / 'decoded'),
+                 features_paths=[first, second], analysis_name='duplicated')
+
+
+def test_a_multi_row_feature_file_is_rejected(dataset, tmp_path):
+    """One row per stimulus, or the rows after it are shifted."""
+    from bdpy.dataform import save_array
+
+    labels = dataset.unique_train_labels
+    store = str(tmp_path / 'store')
+    synthetic.write_feature_tree(store, {LAYER: (4,)}, labels, seed=6)
+    save_array(os.path.join(store, LAYER, '%s.mat' % labels[0]),
+               np.zeros((2, 4), dtype=np.float32), key='feat',
+               dtype=np.float32, sparse=False)
+
+    loader = TrainingFeatureLoader([store])
+    with pytest.raises(ValueError, match='contributes 2 feature rows'):
+        loader.get(LAYER, labels)
+
+
+def test_the_loader_returns_one_row_per_label_in_order(dataset, tmp_path):
+    """Two stores that partition the stimuli: still one row each, in order."""
+    labels = dataset.unique_train_labels
+    shapes = {LAYER: dataset.layer_shapes[LAYER]}
+    first = str(tmp_path / 'store_a')
+    second = str(tmp_path / 'store_b')
+    left = synthetic.write_feature_tree(first, shapes, labels[:2], seed=8)
+    right = synthetic.write_feature_tree(second, shapes, labels[2:], seed=9)
+
+    loader = TrainingFeatureLoader([first, second])
+    requested = list(reversed(labels))
+    features = loader.get(LAYER, requested)
+
+    by_label = dict(zip(labels[:2], left[LAYER]))
+    by_label.update(zip(labels[2:], right[LAYER]))
+    np.testing.assert_array_equal(
+        features, np.vstack([by_label[lb][np.newaxis] for lb in requested]))
+
+
+# The feature index is prediction-side selection ################################
+
+def test_index_path_resolves_inside_the_feature_directory():
+    """Not against the current directory, which is where the old copy looked."""
+    assert resolve_feature_index_path(
+        os.path.join('data', 'features'), 'index_random1000.mat') == os.path.join(
+            'data', 'features', 'index_random1000.mat')
+
+
+def test_absolute_index_path_is_left_alone():
+    absolute = os.path.join(os.sep, 'somewhere', 'index.mat')
+    assert resolve_feature_index_path('data/features', absolute) == absolute
+
+
+def test_the_feature_index_selects_units_at_prediction_time(dataset, tmp_path):
+    """The decoder knows nothing about it, and does not need to.
+
+    Training never reads a feature file, so an index is not part of the
+    decoder: it selects which units of ``F`` prediction combines, and the
+    decoded features come out with the selected width.
+    """
+    decoder_dir = str(tmp_path / 'decoders')
+    decoded_dir = str(tmp_path / 'decoded')
+    index = [3, 0]
+    synthetic.write_feature_index(
+        os.path.join(dataset.train_features_dir, 'index_two.mat'),
+        {LAYER: index})
+
+    _train(dataset, decoder_dir)
+    assert not os.path.exists(os.path.join(decoder_dir, 'feature_index.mat'))
+
+    _predict(dataset, decoder_dir, decoded_dir,
+             feature_index_file='index_two.mat', analysis_name='indexed')
+
+    # The index travels with the decoded features, as in the fastl2lir path.
+    assert os.path.isfile(os.path.join(decoded_dir, 'feature_index.mat'))
+
+    _, test_labels = legacy_ridge.average_test_brain(
+        dataset.brain('sub-01', 'VC', 'test'), dataset.labels('sub-01', 'test'))
+    predicted = pipeline.read_decoded_features(decoded_dir, LAYER, 'sub-01',
+                                               'VC', test_labels)
+    assert predicted.shape == (len(test_labels), len(index))
+
+    train_labels = ridge_factorization.load_train_labels(
+        pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC'))
+    selected = dataset.feature_matrix(LAYER, train_labels)[:, index]
+    expected = _stimulus_coefficients(dataset, decoder_dir) @ selected
+    np.testing.assert_allclose(predicted, expected.astype(np.float32),
+                               rtol=1e-5, atol=1e-6)
diff --git a/tests/test_train_decoder_sklearn_ridge.py b/tests/test_train_decoder_sklearn_ridge.py
new file mode 100644
index 0000000..3ca6282
--- /dev/null
+++ b/tests/test_train_decoder_sklearn_ridge.py
@@ -0,0 +1,184 @@
+"""Behavior of ``train_decoder_sklearn_ridge.featdec_sklearn_ridge_train``.
+
+These tests pin the decoder *artifact contract* -- directory layout, file names,
+array keys, dtypes and the exact definition of the saved normalization
+parameters -- independently of how the regression itself is computed.
+``evaluation.py`` reads ``y_mean.mat``/``y_norm.mat`` straight out of
+``<decoder>/<layer>/<subject>/<roi>/model/``, so this layout is a public
+interface, not an implementation detail.
+"""
+
+from __future__ import annotations
+
+import os
+
+import numpy as np
+import yaml
+
+from tests.helpers import pipeline
+
+ALPHA = 100
+CHUNK_AXIS = 1
+
+
+def test_creates_one_model_dir_per_subject_and_roi(dataset, decoder_dir):
+    """One brain-side model per (subject, ROI), and no per-layer copy of it.
+
+    The fit is on the training stimulus basis, which does not depend on the
+    layer, so there is nothing layer-specific to store.
+    """
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    for subject in dataset.train_fmri:
+        for roi in dataset.rois:
+            directory = pipeline.brain_model_dir(decoder_dir, subject, roi)
+            assert os.path.isdir(directory), directory
+            for key in pipeline.BRAIN_NORM_KEYS:
+                assert os.path.isfile(os.path.join(directory, key + '.mat'))
+            for key in pipeline.FEATURE_NORM_KEYS:
+                assert not os.path.exists(os.path.join(directory,
+                                                       key + '.mat'))
+            assert os.path.isfile(os.path.join(directory, 'info.yaml'))
+
+    # Training does not create the per-layer directories at all: they hold the
+    # statistics sidecars, which prediction writes.
+    assert sorted(os.listdir(decoder_dir)) == sorted(dataset.train_fmri)
+
+
+def test_normalization_parameter_shapes_and_dtype(dataset, decoder_dir):
+    """Only the brain half: the decoder stores no feature statistics."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    for subject in dataset.train_fmri:
+        for roi in dataset.rois:
+            params = pipeline.read_norm_params(decoder_dir, subject, roi)
+            n_voxels = dataset.brain(subject, roi).shape[1]
+            assert params['x_mean'].shape == (1, n_voxels)
+            assert params['x_norm'].shape == (1, n_voxels)
+            for key in pipeline.BRAIN_NORM_KEYS:
+                assert params[key].dtype == np.float32
+
+
+def test_brain_normalization_uses_ddof_one_over_all_trials(dataset,
+                                                           decoder_dir):
+    """Brain mean/SD are per-voxel over *every* trial, with ``ddof=1``."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    subject, roi = 'sub-01', 'VC'
+    brain = dataset.brain(subject, roi, 'train')
+    params = pipeline.read_norm_params(decoder_dir, subject, roi,
+                                       keys=pipeline.BRAIN_NORM_KEYS)
+
+    np.testing.assert_allclose(
+        params['x_mean'], np.mean(brain, axis=0)[np.newaxis, :].astype(
+            np.float32), rtol=1e-6)
+    np.testing.assert_allclose(
+        params['x_norm'],
+        np.std(brain, axis=0, ddof=1)[np.newaxis, :].astype(np.float32),
+        rtol=1e-6)
+    # ddof=0 would be a different number for this data.
+    assert not np.allclose(params['x_norm'],
+                           np.std(brain, axis=0, ddof=0)[np.newaxis, :],
+                           rtol=1e-4)
+
+
+def test_feature_statistics_are_over_unique_stimuli(dataset, decoder_dir,
+                                                    decoded_dir):
+    """``y_mean``/``y_norm`` weight every stimulus once, not once per trial.
+
+    The repetitions in this dataset are unbalanced, so the two definitions
+    genuinely differ and the distinction is observable.  Prediction writes
+    them, into the per-layer directory ``evaluation.py`` reads, because
+    training never opens a feature file.
+    """
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir)
+
+    subject = 'sub-01'
+    brain_labels = dataset.labels(subject, 'train')
+    unique_labels = list(np.unique(brain_labels))
+
+    for layer in dataset.layers:
+        params = pipeline.read_feature_statistics(decoder_dir, layer, subject,
+                                                  'VC')
+        unique_feat = dataset.feature_matrix(layer, unique_labels)
+        trial_feat = dataset.feature_matrix(layer, brain_labels)
+
+        np.testing.assert_allclose(
+            params['y_mean'],
+            np.mean(unique_feat, axis=0)[np.newaxis].astype(np.float32),
+            rtol=1e-5, atol=1e-6)
+        np.testing.assert_allclose(
+            params['y_norm'],
+            np.std(unique_feat, axis=0, ddof=1)[np.newaxis].astype(np.float32),
+            rtol=1e-5, atol=1e-6)
+        assert not np.allclose(
+            params['y_mean'], np.mean(trial_feat, axis=0)[np.newaxis],
+            rtol=1e-4)
+
+
+def test_multi_file_subject_uses_every_file(dataset, decoder_dir):
+    """sub-01's trials are split over two ``.h5`` files; both must be used."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    roi = 'VC'
+    params = pipeline.read_norm_params(decoder_dir, 'sub-01', roi,
+                                       keys=pipeline.BRAIN_NORM_KEYS)
+    all_trials = dataset.brain('sub-01', roi, 'train')
+
+    np.testing.assert_allclose(
+        params['x_mean'],
+        np.mean(all_trials, axis=0)[np.newaxis, :].astype(np.float32),
+        rtol=1e-6)
+    # Only the first file would give a different mean.
+    assert not np.allclose(params['x_mean'],
+                           np.mean(all_trials[:7], axis=0)[np.newaxis, :],
+                           rtol=1e-4)
+
+
+def test_info_yaml_marks_completion(dataset, decoder_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    for subject in dataset.train_fmri:
+        for roi in dataset.rois:
+            info_file = os.path.join(
+                pipeline.brain_model_dir(decoder_dir, subject, roi),
+                'info.yaml')
+            with open(info_file) as f:
+                info = yaml.safe_load(f)
+            assert info['_status']['computation_status'] == 'done'
+
+
+def test_completed_analyses_are_skipped_on_rerun(dataset, decoder_dir):
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    watched = []
+    for subject in dataset.train_fmri:
+        for roi in dataset.rois:
+            directory = pipeline.brain_model_dir(decoder_dir, subject, roi)
+            for name in sorted(os.listdir(directory)):
+                path = os.path.join(directory, name)
+                watched.append((path, os.stat(path).st_mtime_ns))
+
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA,
+                          analysis_name='test_training_rerun')
+
+    for path, mtime in watched:
+        assert os.stat(path).st_mtime_ns == mtime, \
+            '%s was rewritten although the analysis was already done' % path
+
+
+def test_model_artifacts_are_readable(dataset, decoder_dir):
+    """The brain-side ``model/`` holds a loadable model file."""
+    import pickle
+
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    directory = pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC')
+    model_files = sorted(f for f in os.listdir(directory)
+                         if f.endswith('.pkl.gz'))
+    assert model_files
+    for name in model_files:
+        with open(os.path.join(directory, name), 'rb') as f:
+            payload = pickle.load(f)
+        assert 'model' in payload and 'y_shape' in payload
diff --git a/tests/test_training_feature_loading.py b/tests/test_training_feature_loading.py
new file mode 100644
index 0000000..19207b8
--- /dev/null
+++ b/tests/test_training_feature_loading.py
@@ -0,0 +1,170 @@
+"""What the training features cost, and who reads them.
+
+The factorized fit maps brain activity onto the training *stimulus basis*, so it
+needs the trial labels and no feature values: **training never opens a feature
+file**.  Prediction does, and there the invariant is that **no more than one
+feature array is ever alive at a time**: for VGG19 ``conv1_1`` with 1200
+training stimuli a single ``F`` is already ~15 GB.  Asserting on
+``TrainingFeatureLoader.cached_layers`` alone cannot check that -- it is a
+one-element list by construction, and it says nothing about references the
+*caller* still holds -- so these tests track the arrays with ``weakref`` and
+count the ones that are actually still alive.
+"""
+
+from __future__ import annotations
+
+import gc
+import os
+import weakref
+
+import numpy as np
+import pytest
+
+from ridge_factorization import TrainingFeatureLoader
+from tests.helpers import pipeline, synthetic
+
+ALPHA = 100
+CHUNK_AXIS = 1
+
+
+class FeatureArrayTracker(object):
+    """Track how many arrays returned by the loader are still reachable."""
+
+    def __init__(self):
+        self.refs = []
+        self.alive_on_entry = []
+        self.max_alive = 0
+        self.layers = []         # every request, cache hits included
+        self.loaded_layers = []  # only the requests that hit the disk
+
+    def live_count(self) -> int:
+        gc.collect()
+        return len({id(obj) for obj in (ref() for ref in self.refs)
+                    if obj is not None})
+
+    def observe(self, layer, features):
+        self.layers.append(layer)
+        self.refs.append(weakref.ref(features))
+        self.max_alive = max(self.max_alive, self.live_count())
+
+
+@pytest.fixture
+def tracker(monkeypatch):
+    tracker = FeatureArrayTracker()
+    original = TrainingFeatureLoader.get
+
+    def spy(self, layer, labels):
+        tracker.alive_on_entry.append(tracker.live_count())
+        loads_before = self.load_count
+        features = original(self, layer, labels)
+        if self.load_count != loads_before:
+            tracker.loaded_layers.append(layer)
+        tracker.observe(layer, features)
+        return features
+
+    monkeypatch.setattr(TrainingFeatureLoader, 'get', spy)
+    return tracker
+
+
+@pytest.fixture
+def no_feature_reads(monkeypatch):
+    """Make any attempt to read a feature file fail the test."""
+    def forbidden(self, layer, labels):
+        raise AssertionError('the training features were read for layer %r'
+                             % layer)
+
+    monkeypatch.setattr(TrainingFeatureLoader, 'get', forbidden)
+
+
+def test_loader_holds_a_single_layer(dataset):
+    loader = TrainingFeatureLoader([dataset.train_features_dir])
+    labels = dataset.unique_train_labels
+
+    first = loader.get('fc_like', labels)
+    assert loader.cached_layers == ['fc_like']
+
+    # Same request: served from the cache, no extra read.
+    assert loader.get('fc_like', labels) is first
+    assert loader.load_count == 1
+
+    loader.get('conv_like', labels)
+    assert loader.cached_layers == ['conv_like']
+    assert loader.load_count == 2
+
+    loader.release()
+    assert loader.cached_layers == []
+
+
+def test_loader_drops_its_reference_before_loading_the_next_layer(dataset):
+    loader = TrainingFeatureLoader([dataset.train_features_dir])
+    labels = dataset.unique_train_labels
+
+    ref = weakref.ref(loader.get('fc_like', labels))
+    loader.get('conv_like', labels)
+    gc.collect()
+
+    assert ref() is None
+
+
+def test_loader_returns_rows_in_the_requested_order(dataset):
+    loader = TrainingFeatureLoader([dataset.train_features_dir])
+    labels = list(reversed(dataset.unique_train_labels))
+
+    features = loader.get('conv_like', labels)
+
+    np.testing.assert_array_equal(features,
+                                  dataset.feature_matrix('conv_like', labels))
+
+
+def test_training_does_not_open_a_feature_file(dataset, decoder_dir,
+                                               no_feature_reads):
+    """The point of the factorized fit: it needs the labels, not the features.
+
+    Reading them was worth ~3 TB over a DeepRecon run (5 subjects x 9 ROIs x
+    9 layers) purely to recompute identical per-column statistics.
+    """
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    for layer in dataset.layers:
+        assert os.path.isfile(os.path.join(
+            pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC'),
+            'model.pkl.gz'))
+
+
+def test_training_reads_nothing_by_default_and_one_layer_at_a_time(
+        dataset, decoder_dir, tracker):
+    """Same statement measured through the loader rather than an exception."""
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+
+    assert tracker.layers == []
+    assert tracker.loaded_layers == []
+
+
+def test_prediction_never_holds_two_layers_of_features(dataset, decoder_dir,
+                                                       decoded_dir, tracker):
+    synthetic.add_layer(dataset, 'extra_like', (5,))
+
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    assert tracker.max_alive == 1
+    assert tracker.live_count() == 0
+
+
+def test_prediction_reads_each_layer_once(dataset, decoder_dir, decoded_dir,
+                                          tracker):
+    """One disk read per layer, not one per (layer, subject, ROI).
+
+    Every (layer, subject, ROI) asks for the features twice on a first run --
+    once to write the statistics sidecars its decoder directory needs, once to
+    decode -- and the loader serves all of that from the one resident layer.
+    """
+    pipeline.run_training(dataset, decoder_dir, alpha=ALPHA)
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir,
+                            chunk_axis=CHUNK_AXIS)
+
+    combinations = (len(dataset.layers) * len(dataset.test_fmri)
+                    * len(dataset.rois))
+    assert len(tracker.layers) == 2 * combinations
+    assert tracker.loaded_layers == list(reversed(dataset.layers))
diff --git a/tests/test_zero_variance_features.py b/tests/test_zero_variance_features.py
new file mode 100644
index 0000000..c4b7d1e
--- /dev/null
+++ b/tests/test_zero_variance_features.py
@@ -0,0 +1,71 @@
+"""Feature units with zero variance across the training stimuli.
+
+The *direct* decoder z-scored the target as ``(Y - y_mean) / y_norm`` and then
+cleaned up only infinities (``Y[np.isinf(Y)] = 0``).  A unit that is constant
+across the M unique training stimuli has ``y_norm == 0`` and ``Y - y_mean == 0``,
+so the division yielded ``NaN``, which survived the clean-up and scikit-learn
+rejected it.  That behavior is recorded here against the reference
+implementation, which is what the golden fixtures were produced with.
+
+The factorized decoder does not reproduce it. It never normalizes the features
+and decodes a constant unit correctly as ``coeff @ F``, and since training no
+longer reads the features it cannot detect the condition at all. The divergence
+is deliberate and is called out in the PR.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+import ridge_factorization
+
+from tests.helpers import legacy_ridge, pipeline
+
+
+def test_fixture_really_has_a_constant_unit(dataset_with_constant_unit):
+    dataset = dataset_with_constant_unit
+    for layer in dataset.layers:
+        feat = dataset.features[layer].reshape(
+            dataset.features[layer].shape[0], -1)
+        assert np.std(feat[:, 0], ddof=1) == 0.0
+
+
+def test_factorized_training_accepts_constant_feature_units(
+        dataset_with_constant_unit, decoder_dir):
+    """The deliberate divergence: the factorized path no longer rejects them."""
+    pipeline.run_training(dataset_with_constant_unit, decoder_dir,
+                          alpha=100)
+
+    assert ridge_factorization.is_factorized_model_dir(
+        pipeline.brain_model_dir(decoder_dir, 'sub-01', 'VC'))
+
+
+def test_factorized_prediction_decodes_constant_units_exactly(
+        dataset_with_constant_unit, decoder_dir, decoded_dir):
+    """A constant unit decodes to its constant value rather than failing."""
+    dataset = dataset_with_constant_unit
+    pipeline.run_training(dataset, decoder_dir, alpha=100, subjects=['sub-01'],
+                          rois={'VC': 'ROI_VC = 1'})
+    pipeline.run_prediction(dataset, decoder_dir, decoded_dir, chunk_axis=1,
+                            layers=['fc_like'], subjects=['sub-01'],
+                            rois={'VC': 'ROI_VC = 1'})
+
+    predicted = pipeline.read_decoded_features(
+        decoded_dir, 'fc_like', 'sub-01', 'VC', dataset.unique_test_labels)
+    constant = dataset.features['fc_like'].reshape(
+        dataset.features['fc_like'].shape[0], -1)[0, 0]
+    np.testing.assert_allclose(predicted[:, 0], constant, rtol=1e-5)
+
+
+def test_reference_implementation_rejects_them_too(dataset_with_constant_unit):
+    dataset = dataset_with_constant_unit
+    subject, roi, layer = 'sub-01', 'VC', 'fc_like'
+    brain_labels = dataset.labels(subject, 'train')
+    feat_labels = list(np.unique(brain_labels))
+
+    with pytest.raises(ValueError, match='NaN'):
+        legacy_ridge.legacy_train(
+            dataset.brain(subject, roi, 'train'), brain_labels,
+            dataset.feature_matrix(layer, feat_labels), feat_labels,
+            alpha=100, chunk_axis=1)
diff --git a/train_decoder_sklearn_ridge.py b/train_decoder_sklearn_ridge.py
index 0f2e623..9f3d60a 100644
--- a/train_decoder_sklearn_ridge.py
+++ b/train_decoder_sklearn_ridge.py
@@ -1,39 +1,138 @@
-'''DNN Feature decoding - decoders training script.'''
+'''DNN Feature decoding - decoders training script.
 
+Trains a scikit-learn ``Ridge`` decoder in the factorized form described in
+``ridge_factorization``: the model maps brain activity onto the training
+stimulus basis, so the fit needs the trial labels only.  This script reads no
+feature data and takes no feature-side parameter.
+'''
 
-from typing import Dict, List, Optional
+
+from typing import Any, Dict, List, Optional, Sequence
 
 from itertools import product
 import os
-import shutil
 from time import time
 import warnings
 
 import bdpy
 from bdpy.bdata.utils import select_data_multi_bdatas, get_labels_multi_bdatas
-from bdpy.dataform import Features, save_array
-from bdpy.dataform.utils import get_multi_features
+from bdpy.dataform import save_array
 from bdpy.distcomp import DistComp
-from bdpy.ml import ModelTraining
 from bdpy.pipeline.config import init_hydra_cfg
 from bdpy.util import makedir_ifnot
-from sklearn.linear_model import Ridge
 import numpy as np
 import yaml
 
+from ridge_factorization import (
+    LABELS_DIGEST_SPEC,
+    brain_model_dir,
+    fit_stimulus_ridge,
+    labels_digest,
+    normalize_brain_for_training,
+    save_factorized_model,
+)
+
+
+DTYPE = np.float32
+
+
+# Helpers ####################################################################
+
+def _load_info(info_file: str) -> Dict[str, Any]:
+    '''Read a model ``info.yaml``, retrying while it reads back as empty.'''
+    with open(info_file, 'r') as f:
+        info = yaml.safe_load(f)
+    while info is None:
+        warnings.warn('Failed to load info from %s. Retrying...'
+                      % info_file)
+        with open(info_file, 'r') as f:
+            info = yaml.safe_load(f)
+    return info
+
+
+def _is_done(model_dir: str) -> bool:
+    '''True if ``model_dir`` holds a completed decoder.'''
+    info_file = os.path.join(model_dir, 'info.yaml')
+    if not os.path.exists(info_file):
+        return False
+    info = _load_info(info_file)
+    status = info.get('_status', {}) if isinstance(info, dict) else {}
+    return status.get('computation_status') == 'done'
+
+
+def _mark_done(model_dir: str, computation_id: str) -> None:
+    '''Record the completion status in ``info.yaml``.
+
+    ``bdpy.ml.ModelTraining`` used to write this; the factorized training does
+    not go through it, so the same marker is written here instead.
+    '''
+    info_file = os.path.join(model_dir, 'info.yaml')
+    info = _load_info(info_file) if os.path.exists(info_file) else {}
+    if not isinstance(info, dict):
+        info = {}
+    info.setdefault('_status', {})
+    info['_status'].update({
+        'computation_id': computation_id,
+        'computation_status': 'done',
+    })
+    with open(info_file, 'w') as f:
+        f.write(yaml.dump(info, default_flow_style=False))
+
+
+def _save_norm_params(model_dir: str,
+                      norm_param: Dict[str, np.ndarray]) -> None:
+    '''Save the brain normalization parameters ``x_mean``/``x_norm``.
+
+    The feature statistics are written by ``predict_feature.py``: this script
+    never opens a feature file.
+    '''
+    for key in sorted(norm_param):
+        save_file = os.path.join(model_dir, key + '.mat')
+        if not os.path.exists(save_file):
+            try:
+                save_array(save_file, norm_param[key], key=key, dtype=DTYPE,
+                           sparse=False)
+                print('Saved %s' % save_file)
+            except Exception:
+                warnings.warn('Failed to save %s. Possibly double running.'
+                              % save_file)
+
+
+def _write_brain_decoder(model_dir: str, model, feat_labels: Sequence[Any],
+                         brain_mean: np.ndarray, brain_norm: np.ndarray,
+                         alpha: float, subject: str, roi: str,
+                         n_trials: int, n_voxels: int) -> None:
+    """Write the brain-side decoder of one (subject, ROI)."""
+    metadata = {
+        'alpha': float(alpha),
+        'dtype': np.dtype(DTYPE).name,
+        'labels_digest': labels_digest(feat_labels),
+        'labels_digest_spec': LABELS_DIGEST_SPEC,
+        'n_trials': int(n_trials),
+        'n_voxels': int(n_voxels),
+        'roi': str(roi),
+        'subject': str(subject),
+    }
+
+    # Save normalization parameters
+    # -----------------------------
+    print('Saving normalization parameters.')
+    _save_norm_params(model_dir, {'x_mean': brain_mean, 'x_norm': brain_norm})
+
+    # Save the model
+    # --------------
+    save_factorized_model(model_dir, model, feat_labels, metadata=metadata)
+    print('Saved %s' % model_dir)
+
 
 # Main #######################################################################
 
 def featdec_sklearn_ridge_train(
         fmri_data: Dict[str, List[str]],
-        features_paths: List[str],
         output_dir: str = './feature_decoders',
         rois: Optional[Dict[str, str]] = None,
         label_key: Optional[str] = None,
-        layers: Optional[List[str]] = None,
-        feature_index_file: Optional[str] = None,
         alpha: int = 100,
-        chunk_axis: int = 1,
         analysis_name: str = "feature_decoder_training"
 ):
     '''Feature decoder training.
@@ -41,33 +140,22 @@ def featdec_sklearn_ridge_train(
     Input:
 
     - fmri_data
-    - features_paths
 
     Output:
 
     - output_dir
 
-    Parameters:
-
-    TBA
-
     Note:
 
-    If Y.ndim >= 3, Y is divided into chunks along `chunk_axis`.
-    Note that Y[0] should be sample dimension.
+    One model per (subject, ROI), shared by every feature layer, written to
+    `<output_dir>/<subject>/<roi>/model/`.
     '''
     if rois is None:
         rois = {}
-    if layers is None:
-        layers = []
-
-    layers = layers[::-1]  # Start training from deep layers
 
     # Print info -------------------------------------------------------------
     print('Subjects:        %s' % list(fmri_data.keys()))
     print('ROIs:            %s' % list(rois.keys()))
-    print('Target features: %s' % features_paths)
-    print('Layers:          %s' % layers)
     print('')
 
     # Load data --------------------------------------------------------------
@@ -76,51 +164,37 @@ def featdec_sklearn_ridge_train(
 
     data_brain = {sbj: [bdpy.BData(f) for f in data_files] for sbj, data_files in fmri_data.items()}
 
-    if feature_index_file is not None:
-        data_features = [Features(f, feature_index=os.path.join(f, feature_index_file)) for f in features_paths]
-    else:
-        data_features = [Features(f) for f in features_paths]
-
     # Initialize directories -------------------------------------------------
     makedir_ifnot(output_dir)
     makedir_ifnot('tmp')
 
-    # Save feature index -----------------------------------------------------
-    if feature_index_file is not None:
-        feature_index_save_file = os.path.join(output_dir, 'feature_index.mat')
-        shutil.copy(feature_index_file, feature_index_save_file)
-        print('Saved %s' % feature_index_save_file)
+    # Distributed computation setup ------------------------------------------
+    distcomp_db = os.path.join('./tmp', analysis_name + '.db')
+    distcomp = DistComp(backend='sqlite3', db_path=distcomp_db)
 
     # Analysis loop ----------------------------------------------------------
     print('----------------------------------------')
     print('Analysis loop')
 
-    for layer, sbj, roi in np.random.permutation(list(product(layers, fmri_data.keys(), rois.keys()))):
+    for sbj, roi in np.random.permutation(list(product(fmri_data.keys(), rois.keys()))):
         print('--------------------')
-        print('Feature:    %s' % layer)
         print('Subject:    %s' % sbj)
         print('ROI:        %s' % roi)
 
-        # Setup
-        # -----
-        analysis_id = analysis_name + '-' + sbj + '-' + roi + '-' + layer
-        model_dir = os.path.join(output_dir, layer, sbj, roi, 'model')
+        analysis_id = analysis_name + '-' + sbj + '-' + roi
+        model_dir = brain_model_dir(output_dir, sbj, roi)
+
+        if _is_done(model_dir):
+            print('%s is already done and skipped' % analysis_id)
+            continue
+
         makedir_ifnot(model_dir)
 
-        # Check whether the analysis has been done or not.
-        info_file = os.path.join(model_dir, 'info.yaml')
-        if os.path.exists(info_file):
-            with open(info_file, 'r') as f:
-                info = yaml.safe_load(f)
-            while info is None:
-                warnings.warn('Failed to load info from %s. Retrying...'
-                              % info_file)
-                with open(info_file, 'r') as f:
-                    info = yaml.safe_load(f)
-            if '_status' in info and 'computation_status' in info['_status']:
-                if info['_status']['computation_status'] == 'done':
-                    print('%s is already done and skipped' % analysis_id)
-                    continue
+        # DistComp.lock() returns True if the computation is not locked and
+        # successfully locked.
+        if not distcomp.lock(analysis_id):
+            print('%s is already running. Skipped.' % analysis_id)
+            continue
 
         # Preparing data
         # --------------
@@ -132,9 +206,8 @@ def featdec_sklearn_ridge_train(
         brain = select_data_multi_bdatas(data_brain[sbj], rois[roi])
         brain_labels = get_labels_multi_bdatas(data_brain[sbj], label_key)
 
-        # Features
+        # Unique training stimuli, in the order the training features are read.
         feat_labels = np.unique(brain_labels)
-        feat = get_multi_features(data_features, layer, labels=feat_labels)
 
         # Use brain data that has a label included in feature data
         brain = np.vstack([_b for _b, bl in zip(brain, brain_labels) if bl in feat_labels])
@@ -149,62 +222,25 @@ def featdec_sklearn_ridge_train(
         brain_mean = np.mean(brain, axis=0)[np.newaxis, :]  # np.newaxis was added to match Matlab outputs
         brain_norm = np.std(brain, axis=0, ddof=1)[np.newaxis, :]
 
-        # Normalize features
-        feat_mean = np.mean(feat, axis=0)[np.newaxis, :]
-        feat_norm = np.std(feat, axis=0, ddof=1)[np.newaxis, :]
-
-        # Index to sort features by brain data (matching samples)
-        # -----------------------------------------
-        feat_index = np.array([np.where(np.array(feat_labels) == bl) for bl in brain_labels]).flatten()
-
-        # Save normalization parameters
-        # -----------------------------
-        print('Saving normalization parameters.')
-        norm_param = {
-            'x_mean': brain_mean, 'y_mean': feat_mean,
-            'x_norm': brain_norm, 'y_norm': feat_norm
-        }
-        save_targets = [u'x_mean', u'y_mean', u'x_norm', u'y_norm']
-        for sv in save_targets:
-            save_file = os.path.join(model_dir, sv + '.mat')
-            if not os.path.exists(save_file):
-                try:
-                    save_array(save_file, norm_param[sv], key=sv, dtype=np.float32, sparse=False)
-                    print('Saved %s' % save_file)
-                except Exception:
-                    warnings.warn('Failed to save %s. Possibly double running.' % save_file)
-
-        # Preparing learning
-        # ------------------
-        model = Ridge(alpha=alpha)
-
-        # Distributed computation setup
-        # -----------------------------
-        makedir_ifnot('./tmp')
-        distcomp_db = os.path.join('./tmp', analysis_name + '.db')
-        distcomp = DistComp(backend='sqlite3', db_path=distcomp_db)
-
         # Model training
         # --------------
-        print('Model training')
+        print('Model training (brain -> training stimulus basis)')
         start_time = time()
 
-        train = ModelTraining(model, brain, feat)
-        train.id = analysis_id
+        brain_normalized = normalize_brain_for_training(brain, brain_mean,
+                                                        brain_norm)
+        model = fit_stimulus_ridge(brain_normalized, brain_labels, feat_labels,
+                                   alpha=alpha, dtype=DTYPE)
 
-        train.X_normalize = {'mean': brain_mean, 'std': brain_norm}
-        train.Y_normalize = {'mean': feat_mean, 'std': feat_norm}
-        train.Y_sort = {'index': feat_index}
+        print('Elapsed time (model training): %f' % (time() - start_time))
 
-        train.dtype = np.float32
-        train.chunk_axis = chunk_axis
-        train.save_format = 'pickle'
-        train.save_path = model_dir
-        train.distcomp = distcomp
+        _write_brain_decoder(
+            model_dir, model, feat_labels, brain_mean, brain_norm, alpha,
+            sbj, roi, n_trials=brain.shape[0], n_voxels=brain.shape[1])
 
-        train.run()
+        _mark_done(model_dir, analysis_id)
 
-        print('Total elapsed time (model training): %f' % (time() - start_time))
+        distcomp.unlock(analysis_id)
 
     print('%s finished.' % analysis_name)
 
@@ -229,21 +265,15 @@ def featdec_sklearn_ridge_train(
     }
     label_key = cfg["decoder"]["fmri"]["label_key"]
 
-    training_features = cfg["decoder"]["features"]["paths"]
-    layers = cfg["decoder"]["features"]["layers"]
-    feature_index_file = cfg.decoder.features.get("index_file", None)
-
+    # Nothing under `decoder.features` is read here; `predict_feature.py`
+    # consumes those keys.
     decoder_dir = cfg["decoder"]["path"]
 
     featdec_sklearn_ridge_train(
         training_fmri,
-        training_features,
         output_dir=decoder_dir,
         rois=rois,
         label_key=label_key,
-        layers=layers,
-        feature_index_file=feature_index_file,
         alpha=cfg["decoder"]["parameters"]["alpha"],
-        chunk_axis=cfg["decoder"]["parameters"]["chunk_axis"],
         analysis_name=analysis_name
     )
diff --git a/uv.lock b/uv.lock
index c9ac265..26c6c4e 100644
--- a/uv.lock
+++ b/uv.lock
@@ -218,6 +218,18 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/e7/05/c19819d5e3d95294a6f5947fb9b9629efb316b96de511b418c53d245aae6/cycler-0.12.1-py3-none-any.whl", hash = "sha256:85cef7cff222d8644161529808465972e51340599459b8ac3ccbac5a854e0d30", size = 8321, upload-time = "2023-10-07T05:32:16.783Z" },
 ]
 
+[[package]]
+name = "exceptiongroup"
+version = "1.3.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "typing-extensions", marker = "python_full_version < '3.11'" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/50/79/66800aadf48771f6b62f7eb014e352e5d06856655206165d775e675a02c9/exceptiongroup-1.3.1.tar.gz", hash = "sha256:8b412432c6055b0b7d14c310000ae93352ed6754f70fa8f7c34141f91c4e3219", size = 30371, upload-time = "2025-11-21T23:01:54.787Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/8a/0e/97c33bf5009bdbac74fd2beace167cab3f978feb69cc36f1ef79360d6c4e/exceptiongroup-1.3.1-py3-none-any.whl", hash = "sha256:a7a39a3bd276781e98394987d3a5701d0c4edffb633bb7a5144577f82c773598", size = 16740, upload-time = "2025-11-21T23:01:53.443Z" },
+]
+
 [[package]]
 name = "fastl2lir"
 version = "0.11"
@@ -243,9 +255,21 @@ dependencies = [
     { name = "matplotlib" },
     { name = "numpy" },
     { name = "pyyaml" },
+    { name = "scikit-learn", version = "1.7.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" },
+    { name = "scikit-learn", version = "1.8.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" },
     { name = "tqdm" },
 ]
 
+[package.optional-dependencies]
+dev = [
+    { name = "pytest" },
+]
+
+[package.dev-dependencies]
+dev = [
+    { name = "pytest" },
+]
+
 [package.metadata]
 requires-dist = [
     { name = "bdpy", specifier = ">=0.24" },
@@ -253,9 +277,15 @@ requires-dist = [
     { name = "hydra-core" },
     { name = "matplotlib" },
     { name = "numpy", specifier = "<2.0.0" },
+    { name = "pytest", marker = "extra == 'dev'", specifier = ">=7.0" },
     { name = "pyyaml" },
+    { name = "scikit-learn" },
     { name = "tqdm" },
 ]
+provides-extras = ["dev"]
+
+[package.metadata.requires-dev]
+dev = [{ name = "pytest", specifier = ">=7.0" }]
 
 [[package]]
 name = "fonttools"
@@ -399,6 +429,15 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/c6/50/e0edd38dcd63fb26a8547f13d28f7a008bc4a3fd4eb4ff030673f22ad41a/hydra_core-1.3.2-py3-none-any.whl", hash = "sha256:fa0238a9e31df3373b35b0bfb672c34cc92718d21f81311d8996a16de1141d8b", size = 154547, upload-time = "2023-02-23T18:33:40.801Z" },
 ]
 
+[[package]]
+name = "iniconfig"
+version = "2.3.0"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/72/34/14ca021ce8e5dfedc35312d08ba8bf51fdd999c576889fc2c24cb97f4f10/iniconfig-2.3.0.tar.gz", hash = "sha256:c76315c77db068650d49c5b56314774a7804df16fee4402c1f19d6d15d8c4730", size = 20503, upload-time = "2025-10-18T21:55:43.219Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/cb/b1/3846dd7f199d53cb17f49cba7e651e9ce294d8497c8c150530ed11865bb8/iniconfig-2.3.0-py3-none-any.whl", hash = "sha256:f631c04d2c48c52b84d0d0549c99ff3859c98df65b3101406327ecc7d53fbf12", size = 7484, upload-time = "2025-10-18T21:55:41.639Z" },
+]
+
 [[package]]
 name = "joblib"
 version = "1.5.3"
@@ -890,6 +929,24 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/bc/60/5382c03e1970de634027cee8e1b7d39776b778b81812aaf45b694dfe9e28/pillow-12.2.0-pp311-pypy311_pp73-win_amd64.whl", hash = "sha256:bfa9c230d2fe991bed5318a5f119bd6780cda2915cca595393649fc118ab895e", size = 7080946, upload-time = "2026-04-01T14:46:11.734Z" },
 ]
 
+[[package]]
+name = "pluggy"
+version = "1.6.0"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/f9/e2/3e91f31a7d2b083fe6ef3fa267035b518369d9511ffab804f839851d2779/pluggy-1.6.0.tar.gz", hash = "sha256:7dcc130b76258d33b90f61b658791dede3486c3e6bfb003ee5c9bfb396dd22f3", size = 69412, upload-time = "2025-05-15T12:30:07.975Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/54/20/4d324d65cc6d9205fabedc306948156824eb9f0ee1633355a8f7ec5c66bf/pluggy-1.6.0-py3-none-any.whl", hash = "sha256:e920276dd6813095e9377c0bc5566d94c932c33b27a3e3945d8389c374dd4746", size = 20538, upload-time = "2025-05-15T12:30:06.134Z" },
+]
+
+[[package]]
+name = "pygments"
+version = "2.21.0"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/49/2e/ced460408999b33da6b31b0021b0f37d329e202d4169aeb164493778f25b/pygments-2.21.0.tar.gz", hash = "sha256:610ca751c9bc2492b38eb9a38a7fbc93edbbb2d7182edaf34e66ae493dee5c8c", size = 5005329, upload-time = "2026-08-17T08:02:48.824Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/71/46/17f022dd3e953bf20a04a028a21ec746d942f8d2af30fa0f124fa0e6a684/pygments-2.21.0-py3-none-any.whl", hash = "sha256:2363c69b61c4a97c838da3b130dcd6468f4848992b21a82f2a63ec34377137d9", size = 1250147, upload-time = "2026-08-17T08:02:44.912Z" },
+]
+
 [[package]]
 name = "pyparsing"
 version = "3.3.2"
@@ -899,6 +956,24 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/10/bd/c038d7cc38edc1aa5bf91ab8068b63d4308c66c4c8bb3cbba7dfbc049f9c/pyparsing-3.3.2-py3-none-any.whl", hash = "sha256:850ba148bd908d7e2411587e247a1e4f0327839c40e2e5e6d05a007ecc69911d", size = 122781, upload-time = "2026-01-21T03:57:55.912Z" },
 ]
 
+[[package]]
+name = "pytest"
+version = "9.1.1"
+source = { registry = "https://pypi.org/simple" }
+dependencies = [
+    { name = "colorama", marker = "sys_platform == 'win32'" },
+    { name = "exceptiongroup", marker = "python_full_version < '3.11'" },
+    { name = "iniconfig" },
+    { name = "packaging" },
+    { name = "pluggy" },
+    { name = "pygments" },
+    { name = "tomli", marker = "python_full_version < '3.11'" },
+]
+sdist = { url = "https://files.pythonhosted.org/packages/e4/47/b9efed96c114afcfa3c9d3fe98a76a1d14c74a9e266d397cf6eb64be5e01/pytest-9.1.1.tar.gz", hash = "sha256:1088fbde8f2b49d95a549a195707afa7a76a3ce9bcadc26b6d71f0ffda5fe313", size = 1636369, upload-time = "2026-06-19T10:58:32.857Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/24/25/1de2678b631f5a49215c6c96fff41ba892b0a34df68d6d80292b1b48aa7f/pytest-9.1.1-py3-none-any.whl", hash = "sha256:37a86b45efb9a47a61a36449063e8e18d0cab3161329fc099eb21783169c4f0c", size = 386536, upload-time = "2026-06-19T10:58:31.347Z" },
+]
+
 [[package]]
 name = "python-dateutil"
 version = "2.9.0.post0"
@@ -1245,6 +1320,60 @@ wheels = [
     { url = "https://files.pythonhosted.org/packages/32/d5/f9a850d79b0851d1d4ef6456097579a9005b31fea68726a4ae5f2d82ddd9/threadpoolctl-3.6.0-py3-none-any.whl", hash = "sha256:43a0b8fd5a2928500110039e43a5eed8480b918967083ea48dc3ab9f13c4a7fb", size = 18638, upload-time = "2025-03-13T13:49:21.846Z" },
 ]
 
+[[package]]
+name = "tomli"
+version = "2.4.1"
+source = { registry = "https://pypi.org/simple" }
+sdist = { url = "https://files.pythonhosted.org/packages/22/de/48c59722572767841493b26183a0d1cc411d54fd759c5607c4590b6563a6/tomli-2.4.1.tar.gz", hash = "sha256:7c7e1a961a0b2f2472c1ac5b69affa0ae1132c39adcb67aba98568702b9cc23f", size = 17543, upload-time = "2026-03-25T20:22:03.828Z" }
+wheels = [
+    { url = "https://files.pythonhosted.org/packages/f4/11/db3d5885d8528263d8adc260bb2d28ebf1270b96e98f0e0268d32b8d9900/tomli-2.4.1-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:f8f0fc26ec2cc2b965b7a3b87cd19c5c6b8c5e5f436b984e85f486d652285c30", size = 154704, upload-time = "2026-03-25T20:21:10.473Z" },
+    { url = "https://files.pythonhosted.org/packages/6d/f7/675db52c7e46064a9aa928885a9b20f4124ecb9bc2e1ce74c9106648d202/tomli-2.4.1-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:4ab97e64ccda8756376892c53a72bd1f964e519c77236368527f758fbc36a53a", size = 149454, upload-time = "2026-03-25T20:21:12.036Z" },
+    { url = "https://files.pythonhosted.org/packages/61/71/81c50943cf953efa35bce7646caab3cf457a7d8c030b27cfb40d7235f9ee/tomli-2.4.1-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:96481a5786729fd470164b47cdb3e0e58062a496f455ee41b4403be77cb5a076", size = 237561, upload-time = "2026-03-25T20:21:13.098Z" },
+    { url = "https://files.pythonhosted.org/packages/48/c1/f41d9cb618acccca7df82aaf682f9b49013c9397212cb9f53219e3abac37/tomli-2.4.1-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:5a881ab208c0baf688221f8cecc5401bd291d67e38a1ac884d6736cbcd8247e9", size = 243824, upload-time = "2026-03-25T20:21:14.569Z" },
+    { url = "https://files.pythonhosted.org/packages/22/e4/5a816ecdd1f8ca51fb756ef684b90f2780afc52fc67f987e3c61d800a46d/tomli-2.4.1-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:47149d5bd38761ac8be13a84864bf0b7b70bc051806bc3669ab1cbc56216b23c", size = 242227, upload-time = "2026-03-25T20:21:15.712Z" },
+    { url = "https://files.pythonhosted.org/packages/6b/49/2b2a0ef529aa6eec245d25f0c703e020a73955ad7edf73e7f54ddc608aa5/tomli-2.4.1-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:ec9bfaf3ad2df51ace80688143a6a4ebc09a248f6ff781a9945e51937008fcbc", size = 247859, upload-time = "2026-03-25T20:21:17.001Z" },
+    { url = "https://files.pythonhosted.org/packages/83/bd/6c1a630eaca337e1e78c5903104f831bda934c426f9231429396ce3c3467/tomli-2.4.1-cp311-cp311-win32.whl", hash = "sha256:ff2983983d34813c1aeb0fa89091e76c3a22889ee83ab27c5eeb45100560c049", size = 97204, upload-time = "2026-03-25T20:21:18.079Z" },
+    { url = "https://files.pythonhosted.org/packages/42/59/71461df1a885647e10b6bb7802d0b8e66480c61f3f43079e0dcd315b3954/tomli-2.4.1-cp311-cp311-win_amd64.whl", hash = "sha256:5ee18d9ebdb417e384b58fe414e8d6af9f4e7a0ae761519fb50f721de398dd4e", size = 108084, upload-time = "2026-03-25T20:21:18.978Z" },
+    { url = "https://files.pythonhosted.org/packages/b8/83/dceca96142499c069475b790e7913b1044c1a4337e700751f48ed723f883/tomli-2.4.1-cp311-cp311-win_arm64.whl", hash = "sha256:c2541745709bad0264b7d4705ad453b76ccd191e64aa6f0fc66b69a293a45ece", size = 95285, upload-time = "2026-03-25T20:21:20.309Z" },
+    { url = "https://files.pythonhosted.org/packages/c1/ba/42f134a3fe2b370f555f44b1d72feebb94debcab01676bf918d0cb70e9aa/tomli-2.4.1-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:c742f741d58a28940ce01d58f0ab2ea3ced8b12402f162f4d534dfe18ba1cd6a", size = 155924, upload-time = "2026-03-25T20:21:21.626Z" },
+    { url = "https://files.pythonhosted.org/packages/dc/c7/62d7a17c26487ade21c5422b646110f2162f1fcc95980ef7f63e73c68f14/tomli-2.4.1-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:7f86fd587c4ed9dd76f318225e7d9b29cfc5a9d43de44e5754db8d1128487085", size = 150018, upload-time = "2026-03-25T20:21:23.002Z" },
+    { url = "https://files.pythonhosted.org/packages/5c/05/79d13d7c15f13bdef410bdd49a6485b1c37d28968314eabee452c22a7fda/tomli-2.4.1-cp312-cp312-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:ff18e6a727ee0ab0388507b89d1bc6a22b138d1e2fa56d1ad494586d61d2eae9", size = 244948, upload-time = "2026-03-25T20:21:24.04Z" },
+    { url = "https://files.pythonhosted.org/packages/10/90/d62ce007a1c80d0b2c93e02cab211224756240884751b94ca72df8a875ca/tomli-2.4.1-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:136443dbd7e1dee43c68ac2694fde36b2849865fa258d39bf822c10e8068eac5", size = 253341, upload-time = "2026-03-25T20:21:25.177Z" },
+    { url = "https://files.pythonhosted.org/packages/1a/7e/caf6496d60152ad4ed09282c1885cca4eea150bfd007da84aea07bcc0a3e/tomli-2.4.1-cp312-cp312-musllinux_1_2_aarch64.whl", hash = "sha256:5e262d41726bc187e69af7825504c933b6794dc3fbd5945e41a79bb14c31f585", size = 248159, upload-time = "2026-03-25T20:21:26.364Z" },
+    { url = "https://files.pythonhosted.org/packages/99/e7/c6f69c3120de34bbd882c6fba7975f3d7a746e9218e56ab46a1bc4b42552/tomli-2.4.1-cp312-cp312-musllinux_1_2_x86_64.whl", hash = "sha256:5cb41aa38891e073ee49d55fbc7839cfdb2bc0e600add13874d048c94aadddd1", size = 253290, upload-time = "2026-03-25T20:21:27.46Z" },
+    { url = "https://files.pythonhosted.org/packages/d6/2f/4a3c322f22c5c66c4b836ec58211641a4067364f5dcdd7b974b4c5da300c/tomli-2.4.1-cp312-cp312-win32.whl", hash = "sha256:da25dc3563bff5965356133435b757a795a17b17d01dbc0f42fb32447ddfd917", size = 98141, upload-time = "2026-03-25T20:21:28.492Z" },
+    { url = "https://files.pythonhosted.org/packages/24/22/4daacd05391b92c55759d55eaee21e1dfaea86ce5c571f10083360adf534/tomli-2.4.1-cp312-cp312-win_amd64.whl", hash = "sha256:52c8ef851d9a240f11a88c003eacb03c31fc1c9c4ec64a99a0f922b93874fda9", size = 108847, upload-time = "2026-03-25T20:21:29.386Z" },
+    { url = "https://files.pythonhosted.org/packages/68/fd/70e768887666ddd9e9f5d85129e84910f2db2796f9096aa02b721a53098d/tomli-2.4.1-cp312-cp312-win_arm64.whl", hash = "sha256:f758f1b9299d059cc3f6546ae2af89670cb1c4d48ea29c3cacc4fe7de3058257", size = 95088, upload-time = "2026-03-25T20:21:30.677Z" },
+    { url = "https://files.pythonhosted.org/packages/07/06/b823a7e818c756d9a7123ba2cda7d07bc2dd32835648d1a7b7b7a05d848d/tomli-2.4.1-cp313-cp313-macosx_10_13_x86_64.whl", hash = "sha256:36d2bd2ad5fb9eaddba5226aa02c8ec3fa4f192631e347b3ed28186d43be6b54", size = 155866, upload-time = "2026-03-25T20:21:31.65Z" },
+    { url = "https://files.pythonhosted.org/packages/14/6f/12645cf7f08e1a20c7eb8c297c6f11d31c1b50f316a7e7e1e1de6e2e7b7e/tomli-2.4.1-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:eb0dc4e38e6a1fd579e5d50369aa2e10acfc9cace504579b2faabb478e76941a", size = 149887, upload-time = "2026-03-25T20:21:33.028Z" },
+    { url = "https://files.pythonhosted.org/packages/5c/e0/90637574e5e7212c09099c67ad349b04ec4d6020324539297b634a0192b0/tomli-2.4.1-cp313-cp313-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:c7f2c7f2b9ca6bdeef8f0fa897f8e05085923eb091721675170254cbc5b02897", size = 243704, upload-time = "2026-03-25T20:21:34.51Z" },
+    { url = "https://files.pythonhosted.org/packages/10/8f/d3ddb16c5a4befdf31a23307f72828686ab2096f068eaf56631e136c1fdd/tomli-2.4.1-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:f3c6818a1a86dd6dca7ddcaaf76947d5ba31aecc28cb1b67009a5877c9a64f3f", size = 251628, upload-time = "2026-03-25T20:21:36.012Z" },
+    { url = "https://files.pythonhosted.org/packages/e3/f1/dbeeb9116715abee2485bf0a12d07a8f31af94d71608c171c45f64c0469d/tomli-2.4.1-cp313-cp313-musllinux_1_2_aarch64.whl", hash = "sha256:d312ef37c91508b0ab2cee7da26ec0b3ed2f03ce12bd87a588d771ae15dcf82d", size = 247180, upload-time = "2026-03-25T20:21:37.136Z" },
+    { url = "https://files.pythonhosted.org/packages/d3/74/16336ffd19ed4da28a70959f92f506233bd7cfc2332b20bdb01591e8b1d1/tomli-2.4.1-cp313-cp313-musllinux_1_2_x86_64.whl", hash = "sha256:51529d40e3ca50046d7606fa99ce3956a617f9b36380da3b7f0dd3dd28e68cb5", size = 251674, upload-time = "2026-03-25T20:21:38.298Z" },
+    { url = "https://files.pythonhosted.org/packages/16/f9/229fa3434c590ddf6c0aa9af64d3af4b752540686cace29e6281e3458469/tomli-2.4.1-cp313-cp313-win32.whl", hash = "sha256:2190f2e9dd7508d2a90ded5ed369255980a1bcdd58e52f7fe24b8162bf9fedbd", size = 97976, upload-time = "2026-03-25T20:21:39.316Z" },
+    { url = "https://files.pythonhosted.org/packages/6a/1e/71dfd96bcc1c775420cb8befe7a9d35f2e5b1309798f009dca17b7708c1e/tomli-2.4.1-cp313-cp313-win_amd64.whl", hash = "sha256:8d65a2fbf9d2f8352685bc1364177ee3923d6baf5e7f43ea4959d7d8bc326a36", size = 108755, upload-time = "2026-03-25T20:21:40.248Z" },
+    { url = "https://files.pythonhosted.org/packages/83/7a/d34f422a021d62420b78f5c538e5b102f62bea616d1d75a13f0a88acb04a/tomli-2.4.1-cp313-cp313-win_arm64.whl", hash = "sha256:4b605484e43cdc43f0954ddae319fb75f04cc10dd80d830540060ee7cd0243cd", size = 95265, upload-time = "2026-03-25T20:21:41.219Z" },
+    { url = "https://files.pythonhosted.org/packages/3c/fb/9a5c8d27dbab540869f7c1f8eb0abb3244189ce780ba9cd73f3770662072/tomli-2.4.1-cp314-cp314-macosx_10_15_x86_64.whl", hash = "sha256:fd0409a3653af6c147209d267a0e4243f0ae46b011aa978b1080359fddc9b6cf", size = 155726, upload-time = "2026-03-25T20:21:42.23Z" },
+    { url = "https://files.pythonhosted.org/packages/62/05/d2f816630cc771ad836af54f5001f47a6f611d2d39535364f148b6a92d6b/tomli-2.4.1-cp314-cp314-macosx_11_0_arm64.whl", hash = "sha256:a120733b01c45e9a0c34aeef92bf0cf1d56cfe81ed9d47d562f9ed591a9828ac", size = 149859, upload-time = "2026-03-25T20:21:43.386Z" },
+    { url = "https://files.pythonhosted.org/packages/ce/48/66341bdb858ad9bd0ceab5a86f90eddab127cf8b046418009f2125630ecb/tomli-2.4.1-cp314-cp314-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:559db847dc486944896521f68d8190be1c9e719fced785720d2216fe7022b662", size = 244713, upload-time = "2026-03-25T20:21:44.474Z" },
+    { url = "https://files.pythonhosted.org/packages/df/6d/c5fad00d82b3c7a3ab6189bd4b10e60466f22cfe8a08a9394185c8a8111c/tomli-2.4.1-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:01f520d4f53ef97964a240a035ec2a869fe1a37dde002b57ebc4417a27ccd853", size = 252084, upload-time = "2026-03-25T20:21:45.62Z" },
+    { url = "https://files.pythonhosted.org/packages/00/71/3a69e86f3eafe8c7a59d008d245888051005bd657760e96d5fbfb0b740c2/tomli-2.4.1-cp314-cp314-musllinux_1_2_aarch64.whl", hash = "sha256:7f94b27a62cfad8496c8d2513e1a222dd446f095fca8987fceef261225538a15", size = 247973, upload-time = "2026-03-25T20:21:46.937Z" },
+    { url = "https://files.pythonhosted.org/packages/67/50/361e986652847fec4bd5e4a0208752fbe64689c603c7ae5ea7cb16b1c0ca/tomli-2.4.1-cp314-cp314-musllinux_1_2_x86_64.whl", hash = "sha256:ede3e6487c5ef5d28634ba3f31f989030ad6af71edfb0055cbbd14189ff240ba", size = 256223, upload-time = "2026-03-25T20:21:48.467Z" },
+    { url = "https://files.pythonhosted.org/packages/8c/9a/b4173689a9203472e5467217e0154b00e260621caa227b6fa01feab16998/tomli-2.4.1-cp314-cp314-win32.whl", hash = "sha256:3d48a93ee1c9b79c04bb38772ee1b64dcf18ff43085896ea460ca8dec96f35f6", size = 98973, upload-time = "2026-03-25T20:21:49.526Z" },
+    { url = "https://files.pythonhosted.org/packages/14/58/640ac93bf230cd27d002462c9af0d837779f8773bc03dee06b5835208214/tomli-2.4.1-cp314-cp314-win_amd64.whl", hash = "sha256:88dceee75c2c63af144e456745e10101eb67361050196b0b6af5d717254dddf7", size = 109082, upload-time = "2026-03-25T20:21:50.506Z" },
+    { url = "https://files.pythonhosted.org/packages/d5/2f/702d5e05b227401c1068f0d386d79a589bb12bf64c3d2c72ce0631e3bc49/tomli-2.4.1-cp314-cp314-win_arm64.whl", hash = "sha256:b8c198f8c1805dc42708689ed6864951fd2494f924149d3e4bce7710f8eb5232", size = 96490, upload-time = "2026-03-25T20:21:51.474Z" },
+    { url = "https://files.pythonhosted.org/packages/45/4b/b877b05c8ba62927d9865dd980e34a755de541eb65fffba52b4cc495d4d2/tomli-2.4.1-cp314-cp314t-macosx_10_15_x86_64.whl", hash = "sha256:d4d8fe59808a54658fcc0160ecfb1b30f9089906c50b23bcb4c69eddc19ec2b4", size = 164263, upload-time = "2026-03-25T20:21:52.543Z" },
+    { url = "https://files.pythonhosted.org/packages/24/79/6ab420d37a270b89f7195dec5448f79400d9e9c1826df982f3f8e97b24fd/tomli-2.4.1-cp314-cp314t-macosx_11_0_arm64.whl", hash = "sha256:7008df2e7655c495dd12d2a4ad038ff878d4ca4b81fccaf82b714e07eae4402c", size = 160736, upload-time = "2026-03-25T20:21:53.674Z" },
+    { url = "https://files.pythonhosted.org/packages/02/e0/3630057d8eb170310785723ed5adcdfb7d50cb7e6455f85ba8a3deed642b/tomli-2.4.1-cp314-cp314t-manylinux2014_aarch64.manylinux_2_17_aarch64.manylinux_2_28_aarch64.whl", hash = "sha256:1d8591993e228b0c930c4bb0db464bdad97b3289fb981255d6c9a41aedc84b2d", size = 270717, upload-time = "2026-03-25T20:21:55.129Z" },
+    { url = "https://files.pythonhosted.org/packages/7a/b4/1613716072e544d1a7891f548d8f9ec6ce2faf42ca65acae01d76ea06bb0/tomli-2.4.1-cp314-cp314t-manylinux2014_x86_64.manylinux_2_17_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:734e20b57ba95624ecf1841e72b53f6e186355e216e5412de414e3c51e5e3c41", size = 278461, upload-time = "2026-03-25T20:21:56.228Z" },
+    { url = "https://files.pythonhosted.org/packages/05/38/30f541baf6a3f6df77b3df16b01ba319221389e2da59427e221ef417ac0c/tomli-2.4.1-cp314-cp314t-musllinux_1_2_aarch64.whl", hash = "sha256:8a650c2dbafa08d42e51ba0b62740dae4ecb9338eefa093aa5c78ceb546fcd5c", size = 274855, upload-time = "2026-03-25T20:21:57.653Z" },
+    { url = "https://files.pythonhosted.org/packages/77/a3/ec9dd4fd2c38e98de34223b995a3b34813e6bdadf86c75314c928350ed14/tomli-2.4.1-cp314-cp314t-musllinux_1_2_x86_64.whl", hash = "sha256:504aa796fe0569bb43171066009ead363de03675276d2d121ac1a4572397870f", size = 283144, upload-time = "2026-03-25T20:21:59.089Z" },
+    { url = "https://files.pythonhosted.org/packages/ef/be/605a6261cac79fba2ec0c9827e986e00323a1945700969b8ee0b30d85453/tomli-2.4.1-cp314-cp314t-win32.whl", hash = "sha256:b1d22e6e9387bf4739fbe23bfa80e93f6b0373a7f1b96c6227c32bef95a4d7a8", size = 108683, upload-time = "2026-03-25T20:22:00.214Z" },
+    { url = "https://files.pythonhosted.org/packages/12/64/da524626d3b9cc40c168a13da8335fe1c51be12c0a63685cc6db7308daae/tomli-2.4.1-cp314-cp314t-win_amd64.whl", hash = "sha256:2c1c351919aca02858f740c6d33adea0c5deea37f9ecca1cc1ef9e884a619d26", size = 121196, upload-time = "2026-03-25T20:22:01.169Z" },
+    { url = "https://files.pythonhosted.org/packages/5a/cd/e80b62269fc78fc36c9af5a6b89c835baa8af28ff5ad28c7028d60860320/tomli-2.4.1-cp314-cp314t-win_arm64.whl", hash = "sha256:eab21f45c7f66c13f2a9e0e1535309cee140182a9cdae1e041d02e47291e8396", size = 100393, upload-time = "2026-03-25T20:22:02.137Z" },
+    { url = "https://files.pythonhosted.org/packages/7b/61/cceae43728b7de99d9b847560c262873a1f6c98202171fd5ed62640b494b/tomli-2.4.1-py3-none-any.whl", hash = "sha256:0d85819802132122da43cb86656f8d1f8c6587d54ae7dcaf30e90533028b49fe", size = 14583, upload-time = "2026-03-25T20:22:03.012Z" },
+]
+
 [[package]]
 name = "tqdm"
 version = "4.67.3"