diff --git a/bluemath_tk/__init__.py b/bluemath_tk/__init__.py
index 919b62d..4543c06 100644
--- a/bluemath_tk/__init__.py
+++ b/bluemath_tk/__init__.py
@@ -9,6 +9,7 @@
 # Import specific modules instead of using wildcard imports
 from . import (
     additive,
+    benchmarking,
     config,
     core,
     datamining,
@@ -30,6 +31,7 @@
 # Add __all__ variable to control what gets imported when using `from module import *`.
 __all__ = [
     "additive",
+    "benchmarking",
     "config",
     "core",
     "datamining",
diff --git a/bluemath_tk/benchmarking/__init__.py b/bluemath_tk/benchmarking/__init__.py
new file mode 100644
index 0000000..2f834ab
--- /dev/null
+++ b/bluemath_tk/benchmarking/__init__.py
@@ -0,0 +1,25 @@
+"""Reusable benchmarking infrastructure for BlueMath_tk models."""
+
+from .reconstruction import (
+    AutoencoderReconstruction,
+    BenchmarkMethod,
+    MethodBenchmarkResult,
+    PCAReconstruction,
+    ReconstructionBenchmarkReport,
+    ReconstructionMethod,
+    autoencoder_benchmark_method,
+    pca_benchmark_method,
+    run_reconstruction_benchmark,
+)
+
+__all__ = [
+    "AutoencoderReconstruction",
+    "BenchmarkMethod",
+    "MethodBenchmarkResult",
+    "PCAReconstruction",
+    "ReconstructionBenchmarkReport",
+    "ReconstructionMethod",
+    "autoencoder_benchmark_method",
+    "pca_benchmark_method",
+    "run_reconstruction_benchmark",
+]
diff --git a/bluemath_tk/benchmarking/reconstruction.py b/bluemath_tk/benchmarking/reconstruction.py
new file mode 100644
index 0000000..c7fa5e6
--- /dev/null
+++ b/bluemath_tk/benchmarking/reconstruction.py
@@ -0,0 +1,1624 @@
+"""Common reconstruction benchmarking for PCA and autoencoder models.
+
+This module provides the infrastructure required to compare dimensionality
+reduction methods on exactly the same held-out samples. Membership of the
+train, validation, and test partitions always comes from a
+:class:`~bluemath_tk.validation.chronological.ChronologicalSplit`, so the
+benchmark never creates a random split of its own and the test partition never
+reaches any fitting step.
+
+The framework measures *reconstruction* performance only. A low reconstruction
+error does not establish that a method is scientifically better for a
+downstream task, and the reported latent dimensionality is not a storage
+compression ratio.
+"""
+
+from __future__ import annotations
+
+import hashlib
+import inspect
+import math
+import re
+from collections.abc import Callable, Iterator, Mapping, Sequence
+from contextlib import ExitStack, contextmanager
+from dataclasses import dataclass, field
+from time import perf_counter
+from types import MappingProxyType
+from typing import Any, Protocol
+
+import numpy as np
+import xarray as xr
+
+from ..datamining.pca import PCA
+from ..validation.chronological import (
+    ChronologicalSplit,
+    JsonValue,
+    _canonical_json,
+    _freeze_json,
+    _thaw_json,
+    _validate_json_value,
+)
+
+__all__ = [
+    "AutoencoderReconstruction",
+    "BenchmarkMethod",
+    "MethodBenchmarkResult",
+    "PCAReconstruction",
+    "ReconstructionBenchmarkReport",
+    "ReconstructionMethod",
+    "autoencoder_benchmark_method",
+    "pca_benchmark_method",
+    "run_reconstruction_benchmark",
+]
+
+_SCHEMA_VERSION = 1
+_SUPPORTED_METRICS = ("mae", "mse", "rmse")
+_DEFAULT_METRICS = ("mse", "mae", "rmse")
+_METRIC_REDUCTION = "mean"
+_PCA_VARIABLE = "value"
+_PCA_SAMPLE_DIM = "sample"
+_PARTITION_NAMES = ("train", "validation", "test", "excluded")
+_RESERVED_CONFIGURATION_KEYS = (
+    "fit_kwargs",
+    "predict_kwargs",
+    "shuffle_training_data",
+)
+
+# Effective defaults of BaseDeepLearningModel.fit and .predict that materially
+# change the experiment. They are recorded so that two specifications differing
+# only in an omitted argument still describe the same experiment, and so that
+# changing one of them changes the benchmark identity. The literals are pinned
+# here rather than introspected because this module must import without PyTorch;
+# a regression test asserts they still match the model signatures.
+_BENCHMARK_FIT_DEFAULTS: Mapping[str, JsonValue] = MappingProxyType(
+    {
+        "batch_size": 64,
+        "epochs": 500,
+        "learning_rate": 1e-3,
+        "patience": 20,
+    }
+)
+_BENCHMARK_PREDICT_DEFAULTS: Mapping[str, JsonValue] = MappingProxyType(
+    {
+        "batch_size": 64,
+    }
+)
+# Verbosity only controls progress reporting. It cannot change the fitted model,
+# the reconstruction, or any metric, so it is deliberately excluded from the
+# deterministic identity.
+_IDENTITY_EXCLUDED_KWARGS = frozenset({"verbose"})
+
+
+def _load_reconstruction_error() -> Callable[..., Any]:
+    """Import the accepted BlueMath reconstruction metric implementation.
+
+    The import is deferred because ``bluemath_tk.deeplearning.metrics``
+    requires PyTorch, which is an optional dependency. Importing it lazily
+    keeps ``import bluemath_tk`` usable without the deeplearning extra.
+    """
+    try:
+        from ..deeplearning.metrics import reconstruction_error
+    except ImportError as exc:  # pragma: no cover - depends on installation
+        raise ImportError(
+            "Reconstruction benchmarking reuses bluemath_tk.deeplearning.metrics, "
+            "which requires PyTorch. Install the deeplearning extra with "
+            "pip install 'bluemath-tk[deeplearning]'."
+        ) from exc
+    return reconstruction_error
+
+
+def _validate_positive_integer(name: str, value: Any) -> int:
+    if not isinstance(value, int) or isinstance(value, bool):
+        raise TypeError(f"{name} must be an exact non-Boolean integer.")
+    if value < 1:
+        raise ValueError(f"{name} must be a positive integer; got {value}.")
+    return int(value)
+
+
+def _validate_non_empty_string(name: str, value: Any) -> str:
+    if type(value) is not str:
+        raise TypeError(f"{name} must be an exact built-in string.")
+    if not value.strip():
+        raise ValueError(f"{name} must not be empty or blank.")
+    return value
+
+
+def _validate_boolean(name: str, value: Any) -> bool:
+    if type(value) is not bool:
+        raise TypeError(f"{name} must be an exact built-in boolean.")
+    return value
+
+
+def _validate_finite_float(name: str, value: Any) -> float:
+    if isinstance(value, bool) or not isinstance(value, (int, float)):
+        raise TypeError(f"{name} must be a real number.")
+    number = float(value)
+    if not math.isfinite(number):
+        raise ValueError(f"{name} must be finite.")
+    return number
+
+
+def _validate_configuration(
+    configuration: Any,
+    *,
+    name: str,
+) -> dict[str, JsonValue]:
+    """Validate a user-declared, JSON-serializable method configuration."""
+    if configuration is None:
+        return {}
+    if not isinstance(configuration, Mapping):
+        raise TypeError(f"{name} must be a mapping of JSON-compatible values.")
+    payload = {key: value for key, value in configuration.items()}
+    validated = _validate_json_value(payload, path=name)
+    if not isinstance(validated, dict):
+        raise TypeError(f"{name} must be a JSON object.")
+    return validated
+
+
+def _effective_benchmark_kwargs(
+    validated_kwargs: Mapping[str, JsonValue],
+    defaults: Mapping[str, JsonValue],
+) -> dict[str, JsonValue]:
+    """Return the effective keyword arguments that define the experiment.
+
+    Defaults that materially affect the experiment are filled in, so a
+    specification that omits ``epochs`` records the same identity as one that
+    passes the default explicitly. Output-only settings are dropped.
+    """
+    effective: dict[str, JsonValue] = dict(defaults)
+    effective.update(validated_kwargs)
+    for key in _IDENTITY_EXCLUDED_KWARGS:
+        effective.pop(key, None)
+    return effective
+
+
+def _validate_sample_data(X: Any, *, name: str = "X") -> np.ndarray:
+    """Reject datasets the benchmark cannot compare fairly."""
+    if not isinstance(X, np.ndarray):
+        raise TypeError(f"{name} must be a NumPy array.")
+    if X.ndim < 2:
+        raise ValueError(
+            f"{name} must include a leading sample dimension and at least one "
+            "feature dimension. For tabular data use shape "
+            "(n_samples, n_features)."
+        )
+    if X.shape[0] < 1:
+        raise ValueError(f"{name} must contain at least one sample.")
+    if any(dimension < 1 for dimension in X.shape[1:]):
+        raise ValueError(f"Every per-sample dimension of {name} must be positive.")
+    if not np.issubdtype(X.dtype, np.number):
+        raise TypeError(f"{name} must contain numeric values.")
+    if np.issubdtype(X.dtype, np.complexfloating):
+        raise TypeError(
+            f"{name} must contain real-valued data; the shared reconstruction "
+            "metrics do not support complex arrays."
+        )
+    if not np.isfinite(X).all():
+        raise ValueError(f"{name} must not contain NaN or infinite values.")
+    return X
+
+
+def _validate_metrics(metrics: Any) -> tuple[str, ...]:
+    """Validate the requested metric names, preserving the requested order."""
+    if isinstance(metrics, str):
+        raise TypeError("metrics must be a sequence of metric names, not a string.")
+    if not isinstance(metrics, Sequence):
+        raise TypeError("metrics must be a sequence of metric names.")
+    if not metrics:
+        raise ValueError("metrics must request at least one metric.")
+    names: list[str] = []
+    for metric in metrics:
+        name = _validate_non_empty_string("metric name", metric)
+        if name not in _SUPPORTED_METRICS:
+            raise ValueError(
+                f"Unsupported metric {name!r}; supported metrics are "
+                f"{list(_SUPPORTED_METRICS)}."
+            )
+        if name in names:
+            raise ValueError(f"Duplicate metric requested: {name!r}.")
+        names.append(name)
+    return tuple(names)
+
+
+class ReconstructionMethod(Protocol):
+    """Structural interface every benchmarked reconstruction model provides.
+
+    A benchmark method is fitted on the training partition, optionally given
+    the validation partition, and then asked to reconstruct arbitrary samples
+    that share the training per-sample shape.
+
+    Attributes
+    ----------
+    latent_dimension : int
+        Number of latent scalars retained for one sample.
+    is_fitted : bool
+        Whether this instance has already been fitted. The runner rejects
+        instances that are already fitted so that state cannot leak between
+        benchmark runs.
+    """
+
+    latent_dimension: int
+    is_fitted: bool
+
+    def fit(self, X_train: np.ndarray, X_validation: np.ndarray | None) -> None:
+        """Fit the method on the training partition only."""
+        ...
+
+    def reconstruct(self, X: np.ndarray) -> np.ndarray:
+        """Return reconstructions with exactly the shape of ``X``."""
+        ...
+
+
+def _feature_coordinate_names(sample_shape: tuple[int, ...]) -> list[str]:
+    return [f"feature_{index}" for index in range(len(sample_shape))]
+
+
+def _to_pca_dataset(X: np.ndarray) -> xr.Dataset:
+    """Wrap ``(n_samples, ...)`` data in the Dataset layout the PCA class needs."""
+    sample_shape = tuple(X.shape[1:])
+    coordinate_names = _feature_coordinate_names(sample_shape)
+    dimensions = [_PCA_SAMPLE_DIM, *coordinate_names]
+    coordinates: dict[str, np.ndarray] = {
+        _PCA_SAMPLE_DIM: np.arange(X.shape[0]),
+    }
+    for name, size in zip(coordinate_names, sample_shape):
+        coordinates[name] = np.arange(size)
+    return xr.Dataset({_PCA_VARIABLE: (dimensions, X)}, coords=coordinates)
+
+
+class PCAReconstruction:
+    """Benchmark adapter around :class:`bluemath_tk.datamining.pca.PCA`.
+
+    Samples of shape ``(n_samples, d1, ..., dm)`` are presented to the existing
+    BlueMath PCA implementation as a single stacked variable. Stacking and the
+    inverse reshape both use C order, so sample order and the per-sample shape
+    survive the round trip unchanged.
+
+    Parameters
+    ----------
+    n_components : int
+        Number of principal components to retain. This is the latent
+        dimensionality used for comparison against autoencoders.
+    scale_data : bool, optional
+        When True, the BlueMath PCA standardizes features with a
+        ``StandardScaler`` fitted on the training partition only. Default is
+        False, which leaves only the centering that scikit-learn's PCA performs
+        intrinsically. The conservative default keeps the comparison against
+        autoencoders free of preprocessing that they do not receive.
+
+    Notes
+    -----
+    PCA has no early stopping and no validation-driven model selection, so this
+    adapter ignores the validation partition entirely and fits on the training
+    partition alone.
+    """
+
+    def __init__(self, n_components: int, *, scale_data: bool = False):
+        self.n_components = _validate_positive_integer("n_components", n_components)
+        self.scale_data = _validate_boolean("scale_data", scale_data)
+        self._pca: PCA | None = None
+        self._sample_shape: tuple[int, ...] | None = None
+
+    @property
+    def latent_dimension(self) -> int:
+        """Return the number of retained principal components."""
+        return self.n_components
+
+    @property
+    def is_fitted(self) -> bool:
+        """Return whether the underlying PCA model has been fitted."""
+        return self._pca is not None and bool(self._pca.is_fitted)
+
+    @property
+    def pca(self) -> PCA:
+        """Return the fitted BlueMath PCA instance."""
+        if self._pca is None:
+            raise ValueError("The PCA benchmark method has not been fitted yet.")
+        return self._pca
+
+    def fit(self, X_train: np.ndarray, X_validation: np.ndarray | None = None) -> None:
+        """Fit PCA on the training partition only.
+
+        Parameters
+        ----------
+        X_train : np.ndarray
+            Training samples with shape ``(n_train, ...)``.
+        X_validation : np.ndarray, optional
+            Ignored. PCA performs no validation-driven model selection.
+        """
+        _validate_sample_data(X_train, name="X_train")
+        sample_shape = tuple(X_train.shape[1:])
+        n_features = int(np.prod(sample_shape))
+        # PCA centers the training matrix, so its rank is at most n_train - 1.
+        # Allowing n_train components would retain a numerically null component
+        # and inflate the reported latent dimensionality.
+        available = min(int(X_train.shape[0]) - 1, n_features)
+        if self.n_components > available:
+            raise ValueError(
+                f"n_components={self.n_components} exceeds the {available} "
+                "components available from a training partition of "
+                f"{X_train.shape[0]} samples with {n_features} scalars per "
+                "sample."
+            )
+
+        self._sample_shape = sample_shape
+        self._pca = PCA(n_components=self.n_components)
+        self._pca.fit(
+            data=_to_pca_dataset(X_train),
+            vars_to_stack=[_PCA_VARIABLE],
+            coords_to_stack=_feature_coordinate_names(sample_shape),
+            pca_dim_for_rows=_PCA_SAMPLE_DIM,
+            scale_data=self.scale_data,
+        )
+
+    def reconstruct(self, X: np.ndarray) -> np.ndarray:
+        """Project ``X`` onto the fitted components and invert the projection."""
+        if self._pca is None or self._sample_shape is None:
+            raise ValueError(
+                "The PCA benchmark method must be fitted before reconstructing."
+            )
+        _validate_sample_data(X, name="X")
+        if tuple(X.shape[1:]) != self._sample_shape:
+            raise ValueError(
+                f"Expected per-sample shape {self._sample_shape}, got "
+                f"{tuple(X.shape[1:])}."
+            )
+        principal_components = self._pca.transform(data=_to_pca_dataset(X))
+        reconstructed = self._pca.inverse_transform(PCs=principal_components)
+        return np.asarray(reconstructed[_PCA_VARIABLE].values, dtype=np.float64)
+
+
+_FORBIDDEN_FIT_KWARGS = (
+    "X",
+    "y",
+    "criterion",
+    "optimizer",
+    "validation_data",
+    "validation_split",
+)
+
+
+def _reject_forbidden_fit_kwargs(fit_kwargs: Mapping[str, Any]) -> None:
+    """Reject fit arguments the benchmark must control itself."""
+    if not isinstance(fit_kwargs, Mapping):
+        raise TypeError("fit_kwargs must be a mapping.")
+    forbidden = sorted(set(_FORBIDDEN_FIT_KWARGS).intersection(fit_kwargs))
+    if forbidden:
+        raise ValueError(
+            "The benchmark controls partition membership and builds a fresh "
+            f"model for every run; remove these fit_kwargs: {forbidden}."
+        )
+
+
+def _reject_stochastic_prediction(predict_kwargs: Mapping[str, Any]) -> None:
+    """Reject stochastic reconstruction, which is not comparable."""
+    if not isinstance(predict_kwargs, Mapping):
+        raise TypeError("predict_kwargs must be a mapping.")
+    if predict_kwargs.get("stochastic"):
+        raise ValueError(
+            "Stochastic reconstruction is not comparable with the "
+            "deterministic PCA and autoencoder reconstructions used by this "
+            "benchmark. Remove stochastic=True from predict_kwargs."
+        )
+
+
+def _require_explicit_validation_data_support(model: Any) -> None:
+    """Reject models whose ``fit`` would swallow ``validation_data``.
+
+    A ``fit(self, X, **kwargs)`` signature accepts ``validation_data`` silently
+    and is free to build its own random validation split, which is exactly the
+    leakage this framework exists to prevent. Requiring the parameter to be
+    declared makes that impossible to do by accident.
+    """
+    try:
+        signature = inspect.signature(model.fit)
+    except (TypeError, ValueError):  # pragma: no cover - exotic callables
+        return
+    parameter = signature.parameters.get("validation_data")
+    if parameter is None or parameter.kind is inspect.Parameter.VAR_KEYWORD:
+        raise TypeError(
+            f"{type(model).__name__}.fit() must declare an explicit "
+            "validation_data parameter. Without it the benchmark cannot prove "
+            "that the chronological validation partition is the one actually "
+            "used, because **kwargs would silently absorb it."
+        )
+
+
+class AutoencoderReconstruction:
+    """Benchmark adapter around a BlueMath autoencoder.
+
+    The adapter uses the accepted public workflow only: ``model.fit(...)`` with
+    explicit chronological validation data, and ``model.predict(...)`` for
+    deterministic reconstruction. Model architectures are never modified.
+
+    Parameters
+    ----------
+    model : object
+        An unfitted BlueMath autoencoder exposing ``fit`` and ``predict``.
+    latent_dimension : int
+        The latent width declared for this model. When the model exposes ``k``
+        the two values must agree.
+    shuffle_training_data : bool, optional
+        When True (the default), the training samples are permuted once before
+        fitting so that mini-batches are not contiguous blocks of adjacent
+        timestamps. The permutation is drawn inside the benchmark's isolated
+        random state, so it is controlled by the run ``seed`` and never touches
+        the caller's random state. Validation membership is unaffected. Set to
+        False to train on the chronological order exactly as supplied.
+    fit_kwargs : dict, optional
+        Extra keyword arguments forwarded to ``model.fit``. ``validation_data``
+        and ``validation_split`` are rejected because the benchmark controls
+        partition membership, and ``optimizer`` and ``criterion`` are rejected
+        because a stateful object built for one model instance must not be
+        reused by the fresh instance of another run. ``verbose`` defaults to 0.
+    predict_kwargs : dict, optional
+        Extra keyword arguments forwarded to ``model.predict``. ``verbose``
+        defaults to 0.
+
+    Notes
+    -----
+    For variational autoencoders, ``predict`` defaults to the deterministic
+    posterior-mean reconstruction. ``stochastic=True`` is rejected here because
+    comparing a single stochastic draw against deterministic PCA and
+    autoencoder reconstructions is not a like-for-like measurement.
+    """
+
+    _FORBIDDEN_FIT_KWARGS = _FORBIDDEN_FIT_KWARGS
+
+    def __init__(
+        self,
+        model: Any,
+        latent_dimension: int,
+        *,
+        shuffle_training_data: bool = True,
+        fit_kwargs: Mapping[str, Any] | None = None,
+        predict_kwargs: Mapping[str, Any] | None = None,
+    ):
+        for attribute in ("fit", "predict"):
+            if not callable(getattr(model, attribute, None)):
+                raise TypeError(
+                    f"model must expose a callable {attribute}() method to be "
+                    "benchmarked as an autoencoder."
+                )
+        _require_explicit_validation_data_support(model)
+        self.model = model
+        self.shuffle_training_data = _validate_boolean(
+            "shuffle_training_data",
+            shuffle_training_data,
+        )
+        self._latent_dimension = _validate_positive_integer(
+            "latent_dimension",
+            latent_dimension,
+        )
+        declared_k = getattr(model, "k", None)
+        if declared_k is not None and int(declared_k) != self._latent_dimension:
+            raise ValueError(
+                f"latent_dimension={self._latent_dimension} contradicts the "
+                f"model latent width k={int(declared_k)}."
+            )
+
+        _reject_forbidden_fit_kwargs(fit_kwargs or {})
+        self._fit_kwargs = dict(fit_kwargs or {})
+        self._fit_kwargs.setdefault("verbose", 0)
+
+        _reject_stochastic_prediction(predict_kwargs or {})
+        self._predict_kwargs = dict(predict_kwargs or {})
+        self._predict_kwargs.setdefault("verbose", 0)
+        self._history: dict[str, list] | None = None
+
+    @property
+    def latent_dimension(self) -> int:
+        """Return the declared latent width of the wrapped model."""
+        return self._latent_dimension
+
+    @property
+    def is_fitted(self) -> bool:
+        """Return whether the wrapped model reports itself as fitted."""
+        return bool(getattr(self.model, "is_fitted", False))
+
+    @property
+    def history(self) -> dict[str, list] | None:
+        """Return the training history returned by the last fit call."""
+        return self._history
+
+    def fit(self, X_train: np.ndarray, X_validation: np.ndarray | None = None) -> None:
+        """Fit on the training partition, validating on the given samples only.
+
+        Parameters
+        ----------
+        X_train : np.ndarray
+            Training samples. Every one of them is used for optimisation. When
+            ``shuffle_training_data`` is True they are permuted first, so that
+            mini-batches are not contiguous blocks of adjacent timestamps.
+        X_validation : np.ndarray
+            Validation samples. Exactly these samples drive the validation loss
+            and early stopping. Their membership is never changed by shuffling.
+        """
+        _validate_sample_data(X_train, name="X_train")
+        if X_validation is None:
+            raise ValueError(
+                "Autoencoder benchmarking requires the validation partition for "
+                "early stopping. Declare uses_validation_partition=True."
+            )
+        _validate_sample_data(X_validation, name="X_validation")
+        if self.shuffle_training_data:
+            X_train = X_train[np.random.permutation(len(X_train))]
+        self._history = self.model.fit(
+            X_train,
+            validation_data=(X_validation, None),
+            **self._fit_kwargs,
+        )
+
+    def reconstruct(self, X: np.ndarray) -> np.ndarray:
+        """Return the model's deterministic reconstruction of ``X``."""
+        _validate_sample_data(X, name="X")
+        return np.asarray(self.model.predict(X, **self._predict_kwargs))
+
+
+@dataclass(frozen=True)
+class BenchmarkMethod:
+    """Reproducible specification of one benchmarked reconstruction method.
+
+    Attributes
+    ----------
+    name : str
+        Unique, human-readable identifier used in the benchmark report.
+    method_type : str
+        Explicit, user-supplied family label such as ``"pca"`` or
+        ``"autoencoder"``. The benchmark never infers this by introspecting the
+        factory, because callables cannot be serialized reproducibly.
+    latent_dimension : int
+        Latent scalars retained per sample. The runner cross-checks this
+        against the value reported by the constructed method.
+    factory : callable
+        Zero-argument callable returning a fresh, unfitted
+        :class:`ReconstructionMethod`. A new instance is built for every run so
+        that fitted state cannot leak between runs.
+    configuration : mapping, optional
+        JSON-compatible description of the method configuration, recorded
+        verbatim in the report.
+    uses_validation_partition : bool, optional
+        Whether the method consumes the validation partition, for example for
+        early stopping. Default is True. PCA declares False because it performs
+        no validation-driven model selection.
+    """
+
+    name: str
+    method_type: str
+    latent_dimension: int
+    factory: Callable[[], ReconstructionMethod]
+    configuration: Mapping[str, JsonValue] = field(default_factory=dict)
+    uses_validation_partition: bool = True
+
+    def __post_init__(self) -> None:
+        """Validate and freeze the specification after construction."""
+        object.__setattr__(self, "name", _validate_non_empty_string("name", self.name))
+        object.__setattr__(
+            self,
+            "method_type",
+            _validate_non_empty_string("method_type", self.method_type),
+        )
+        object.__setattr__(
+            self,
+            "latent_dimension",
+            _validate_positive_integer("latent_dimension", self.latent_dimension),
+        )
+        if not callable(self.factory):
+            raise TypeError("factory must be a zero-argument callable.")
+        object.__setattr__(
+            self,
+            "uses_validation_partition",
+            _validate_boolean(
+                "uses_validation_partition",
+                self.uses_validation_partition,
+            ),
+        )
+        object.__setattr__(
+            self,
+            "configuration",
+            _freeze_json(
+                _validate_configuration(self.configuration, name="configuration")
+            ),
+        )
+
+
+def pca_benchmark_method(
+    name: str,
+    *,
+    n_components: int,
+    scale_data: bool = False,
+) -> BenchmarkMethod:
+    """Build a PCA benchmark specification using the existing BlueMath PCA.
+
+    Parameters
+    ----------
+    name : str
+        Unique identifier for this method in the report.
+    n_components : int
+        Number of principal components, used as the latent dimensionality.
+    scale_data : bool, optional
+        Standardize features using a scaler fitted on the training partition
+        only. Default is False.
+
+    Returns
+    -------
+    BenchmarkMethod
+        A specification whose factory builds a fresh
+        :class:`PCAReconstruction`.
+    """
+    components = _validate_positive_integer("n_components", n_components)
+    scale = _validate_boolean("scale_data", scale_data)
+
+    def factory() -> ReconstructionMethod:
+        return PCAReconstruction(n_components=components, scale_data=scale)
+
+    return BenchmarkMethod(
+        name=name,
+        method_type="pca",
+        latent_dimension=components,
+        factory=factory,
+        configuration={
+            "implementation": "bluemath_tk.datamining.pca.PCA",
+            "n_components": components,
+            "scale_data": scale,
+        },
+        uses_validation_partition=False,
+    )
+
+
+def autoencoder_benchmark_method(
+    name: str,
+    *,
+    model_factory: Callable[[], Any],
+    latent_dimension: int,
+    configuration: Mapping[str, JsonValue] | None = None,
+    shuffle_training_data: bool = True,
+    fit_kwargs: Mapping[str, Any] | None = None,
+    predict_kwargs: Mapping[str, Any] | None = None,
+) -> BenchmarkMethod:
+    """Build an autoencoder benchmark specification.
+
+    Parameters
+    ----------
+    name : str
+        Unique identifier for this method in the report.
+    model_factory : callable
+        Zero-argument callable returning a fresh, unfitted BlueMath
+        autoencoder. Supplying a factory rather than an instance guarantees
+        that no fitted state is shared between runs.
+    latent_dimension : int
+        Latent width of the model, cross-checked against ``model.k`` when
+        available.
+    configuration : mapping, optional
+        JSON-compatible description of the architecture and hyperparameters,
+        recorded verbatim. The factory itself is never introspected, so this is
+        where the architecture must be described. The reserved keys
+        ``"fit_kwargs"``, ``"predict_kwargs"``, and ``"shuffle_training_data"``
+        are filled in automatically and must not be supplied.
+    shuffle_training_data : bool, optional
+        Permute the training samples once before fitting, so that mini-batches
+        are not contiguous blocks of adjacent timestamps. Default is True. The
+        choice is recorded in the method configuration.
+    fit_kwargs : mapping, optional
+        Extra keyword arguments for ``model.fit``. Values must be
+        JSON-compatible so that the effective training configuration can be
+        recorded, because a configuration that cannot be recorded cannot be
+        part of a reproducible identity. The effective values, including
+        defaults such as ``epochs`` and ``learning_rate``, are recorded under
+        ``configuration["fit_kwargs"]`` and therefore change the benchmark
+        identity. ``verbose`` is excluded because it only controls progress
+        reporting.
+    predict_kwargs : mapping, optional
+        Extra keyword arguments for ``model.predict``, recorded the same way
+        under ``configuration["predict_kwargs"]``.
+
+    Returns
+    -------
+    BenchmarkMethod
+        A specification whose factory builds a fresh
+        :class:`AutoencoderReconstruction`.
+    """
+    if not callable(model_factory):
+        raise TypeError("model_factory must be a zero-argument callable.")
+    width = _validate_positive_integer("latent_dimension", latent_dimension)
+    shuffle = _validate_boolean("shuffle_training_data", shuffle_training_data)
+
+    # Reject the benchmark-controlled arguments before JSON validation, so the
+    # error names the real problem instead of complaining about serialization.
+    _reject_forbidden_fit_kwargs(fit_kwargs or {})
+    _reject_stochastic_prediction(predict_kwargs or {})
+
+    # _validate_configuration rebuilds every container, so the recorded values
+    # are a deep copy. Later mutation of the caller's mappings, at any nesting
+    # depth, cannot change what this specification runs or records.
+    validated_fit = _validate_configuration(fit_kwargs, name="fit_kwargs")
+    validated_predict = _validate_configuration(predict_kwargs, name="predict_kwargs")
+    frozen_fit = _freeze_json(validated_fit)
+    frozen_predict = _freeze_json(validated_predict)
+
+    recorded = _validate_configuration(configuration, name="configuration")
+    reserved = sorted(set(_RESERVED_CONFIGURATION_KEYS).intersection(recorded))
+    if reserved:
+        raise ValueError(
+            f"configuration must not set the reserved keys {reserved}; they are "
+            "recorded automatically from the shuffle_training_data, fit_kwargs, "
+            "and predict_kwargs arguments."
+        )
+    recorded["shuffle_training_data"] = shuffle
+    recorded["fit_kwargs"] = _effective_benchmark_kwargs(
+        validated_fit,
+        _BENCHMARK_FIT_DEFAULTS,
+    )
+    recorded["predict_kwargs"] = _effective_benchmark_kwargs(
+        validated_predict,
+        _BENCHMARK_PREDICT_DEFAULTS,
+    )
+
+    def factory() -> ReconstructionMethod:
+        # Thawing rebuilds plain containers on every call, so no run can observe
+        # or mutate the keyword arguments used by another run.
+        return AutoencoderReconstruction(
+            model_factory(),
+            latent_dimension=width,
+            shuffle_training_data=shuffle,
+            fit_kwargs=_thaw_json(frozen_fit),
+            predict_kwargs=_thaw_json(frozen_predict),
+        )
+
+    return BenchmarkMethod(
+        name=name,
+        method_type="autoencoder",
+        latent_dimension=width,
+        factory=factory,
+        configuration=recorded,
+        uses_validation_partition=True,
+    )
+
+
+@dataclass(frozen=True)
+class MethodBenchmarkResult:
+    """Reconstruction result for one method on the test partition.
+
+    Attributes
+    ----------
+    name : str
+        Method identifier.
+    method_type : str
+        Declared method family.
+    latent_dimension : int
+        Latent scalars retained per sample.
+    uses_validation_partition : bool
+        Whether the method consumed the validation partition.
+    original_scalars_per_sample : int
+        Number of scalars in one input sample.
+    latent_scalars_per_sample : int
+        Number of scalars in one latent representation.
+    latent_dimensionality_ratio : float
+        ``latent_scalars_per_sample / original_scalars_per_sample``. This is a
+        dimensionality ratio only. It is not a bitrate, a storage compression
+        ratio, or a compressed file size, because it ignores latent precision,
+        quantisation, entropy coding, and model parameter storage.
+    test_metrics : mapping
+        Reconstruction metrics computed on the test partition only.
+    fit_seconds : float
+        Observed wall-clock fitting time from a monotonic clock.
+    reconstruction_seconds : float
+        Observed wall-clock reconstruction time from a monotonic clock.
+    configuration : mapping
+        The specification configuration recorded verbatim.
+    """
+
+    name: str
+    method_type: str
+    latent_dimension: int
+    uses_validation_partition: bool
+    original_scalars_per_sample: int
+    latent_scalars_per_sample: int
+    latent_dimensionality_ratio: float
+    test_metrics: Mapping[str, float]
+    fit_seconds: float
+    reconstruction_seconds: float
+    configuration: Mapping[str, JsonValue] = field(default_factory=dict)
+
+    def __post_init__(self) -> None:
+        """Validate and freeze the result after construction."""
+        object.__setattr__(self, "name", _validate_non_empty_string("name", self.name))
+        object.__setattr__(
+            self,
+            "method_type",
+            _validate_non_empty_string("method_type", self.method_type),
+        )
+        for attribute in (
+            "latent_dimension",
+            "original_scalars_per_sample",
+            "latent_scalars_per_sample",
+        ):
+            object.__setattr__(
+                self,
+                attribute,
+                _validate_positive_integer(attribute, getattr(self, attribute)),
+            )
+        if self.latent_dimension != self.latent_scalars_per_sample:
+            raise ValueError(
+                f"latent_dimension {self.latent_dimension} must equal "
+                f"latent_scalars_per_sample {self.latent_scalars_per_sample}."
+            )
+        object.__setattr__(
+            self,
+            "uses_validation_partition",
+            _validate_boolean(
+                "uses_validation_partition",
+                self.uses_validation_partition,
+            ),
+        )
+        ratio = _validate_finite_float(
+            "latent_dimensionality_ratio",
+            self.latent_dimensionality_ratio,
+        )
+        expected_ratio = (
+            self.latent_scalars_per_sample / self.original_scalars_per_sample
+        )
+        if ratio != expected_ratio:
+            raise ValueError(
+                "latent_dimensionality_ratio must equal "
+                "latent_scalars_per_sample / original_scalars_per_sample "
+                f"({expected_ratio!r}); got {ratio!r}."
+            )
+        object.__setattr__(self, "latent_dimensionality_ratio", ratio)
+        for attribute in ("fit_seconds", "reconstruction_seconds"):
+            seconds = _validate_finite_float(attribute, getattr(self, attribute))
+            if seconds < 0:
+                raise ValueError(f"{attribute} must not be negative.")
+            object.__setattr__(self, attribute, seconds)
+
+        if not isinstance(self.test_metrics, Mapping) or not self.test_metrics:
+            raise TypeError("test_metrics must be a non-empty mapping.")
+        metrics: dict[str, float] = {}
+        for key, value in self.test_metrics.items():
+            metric = _validate_non_empty_string("test_metrics key", key)
+            if metric not in _SUPPORTED_METRICS:
+                raise ValueError(f"Unsupported metric in test_metrics: {metric!r}.")
+            score = _validate_finite_float(f"test_metrics[{metric!r}]", value)
+            if score < 0.0:
+                raise ValueError(
+                    f"test_metrics[{metric!r}] must not be negative; MSE, MAE, "
+                    f"and RMSE are non-negative by construction. Got {score!r}."
+                )
+            metrics[metric] = score
+        object.__setattr__(self, "test_metrics", _freeze_json(metrics))
+        object.__setattr__(
+            self,
+            "configuration",
+            _freeze_json(
+                _validate_configuration(self.configuration, name="configuration")
+            ),
+        )
+
+    def to_dict(self) -> dict[str, JsonValue]:
+        """Return a JSON-compatible dictionary describing this result."""
+        return {
+            "name": self.name,
+            "method_type": self.method_type,
+            "latent_dimension": self.latent_dimension,
+            "uses_validation_partition": self.uses_validation_partition,
+            "original_scalars_per_sample": self.original_scalars_per_sample,
+            "latent_scalars_per_sample": self.latent_scalars_per_sample,
+            "latent_dimensionality_ratio": self.latent_dimensionality_ratio,
+            "test_metrics": _thaw_json(self.test_metrics),
+            "timing": {
+                "fit_seconds": self.fit_seconds,
+                "reconstruction_seconds": self.reconstruction_seconds,
+            },
+            "configuration": _thaw_json(self.configuration),
+        }
+
+    def identity(self) -> dict[str, JsonValue]:
+        """Return the configuration identity, excluding measured outcomes."""
+        return {
+            "name": self.name,
+            "method_type": self.method_type,
+            "latent_dimension": self.latent_dimension,
+            "uses_validation_partition": self.uses_validation_partition,
+            "original_scalars_per_sample": self.original_scalars_per_sample,
+            "latent_scalars_per_sample": self.latent_scalars_per_sample,
+            "configuration": _thaw_json(self.configuration),
+        }
+
+    @classmethod
+    def from_dict(cls, payload: Mapping[str, Any]) -> MethodBenchmarkResult:
+        """Construct a validated result from a dictionary."""
+        if not isinstance(payload, Mapping):
+            raise TypeError("Result payload must be a mapping.")
+        required = {
+            "name",
+            "method_type",
+            "latent_dimension",
+            "uses_validation_partition",
+            "original_scalars_per_sample",
+            "latent_scalars_per_sample",
+            "latent_dimensionality_ratio",
+            "test_metrics",
+            "timing",
+            "configuration",
+        }
+        _require_exact_fields(payload, required=required, label="Result")
+        timing = payload["timing"]
+        if type(timing) is not dict:
+            raise TypeError("Result timing must be an exact JSON object.")
+        _require_exact_fields(
+            timing,
+            required={"fit_seconds", "reconstruction_seconds"},
+            label="Result timing",
+        )
+        if type(payload["test_metrics"]) is not dict:
+            raise TypeError("Result test_metrics must be an exact JSON object.")
+        if type(payload["configuration"]) is not dict:
+            raise TypeError("Result configuration must be an exact JSON object.")
+        return cls(
+            name=payload["name"],
+            method_type=payload["method_type"],
+            latent_dimension=payload["latent_dimension"],
+            uses_validation_partition=payload["uses_validation_partition"],
+            original_scalars_per_sample=payload["original_scalars_per_sample"],
+            latent_scalars_per_sample=payload["latent_scalars_per_sample"],
+            latent_dimensionality_ratio=payload["latent_dimensionality_ratio"],
+            test_metrics=payload["test_metrics"],
+            fit_seconds=timing["fit_seconds"],
+            reconstruction_seconds=timing["reconstruction_seconds"],
+            configuration=payload["configuration"],
+        )
+
+
+def _require_exact_fields(
+    payload: Mapping[str, Any],
+    *,
+    required: set[str],
+    label: str,
+) -> None:
+    if any(type(key) is not str for key in payload):
+        raise TypeError(f"{label} field names must be exact built-in strings.")
+    missing = sorted(required.difference(payload))
+    extra = sorted(set(payload).difference(required))
+    if missing:
+        raise ValueError(f"{label} is missing required fields: {missing}.")
+    if extra:
+        raise ValueError(f"{label} contains unsupported fields: {extra}.")
+
+
+@dataclass(frozen=True)
+class ReconstructionBenchmarkReport:
+    """Complete record of one reconstruction benchmark run.
+
+    Attributes
+    ----------
+    schema_version : int
+        Report schema version.
+    n_samples : int
+        Total samples in the benchmarked dataset.
+    sample_shape : tuple of int
+        Per-sample shape, excluding the leading sample dimension.
+    data_digest : str
+        SHA-256 digest of the benchmarked values, their dtype, and their shape.
+        The split manifest fingerprints the time axis only, so this is what
+        establishes that two runs compared the same data.
+    partition_sizes : mapping
+        Sample counts for the train, validation, test, and excluded partitions.
+    metrics : tuple of str
+        Metric names in the order they were requested.
+    seed : int or None
+        Seed applied inside an isolated random state, if any.
+    split_identity : mapping
+        Stable identity of the chronological split that produced the
+        partitions.
+    results : tuple of MethodBenchmarkResult
+        One result per benchmarked method, in specification order.
+
+    Notes
+    -----
+    The report deliberately provides no ranking or "best method" field. It
+    measures reconstruction error on held-out samples, which is not the same as
+    downstream scientific skill.
+    """
+
+    schema_version: int
+    n_samples: int
+    sample_shape: tuple[int, ...]
+    data_digest: str
+    partition_sizes: Mapping[str, int]
+    metrics: tuple[str, ...]
+    seed: int | None
+    split_identity: Mapping[str, JsonValue]
+    results: tuple[MethodBenchmarkResult, ...]
+
+    def __post_init__(self) -> None:
+        """Validate and freeze the report after construction."""
+        if type(self.schema_version) is not int:
+            raise TypeError("schema_version must be an exact non-Boolean integer.")
+        if self.schema_version != _SCHEMA_VERSION:
+            raise ValueError(
+                f"Unsupported benchmark report schema version "
+                f"{self.schema_version}; expected {_SCHEMA_VERSION}."
+            )
+        object.__setattr__(
+            self,
+            "n_samples",
+            _validate_positive_integer("n_samples", self.n_samples),
+        )
+        if isinstance(self.sample_shape, (str, bytes)) or not isinstance(
+            self.sample_shape, Sequence
+        ):
+            raise TypeError("sample_shape must be a sequence of positive integers.")
+        object.__setattr__(
+            self,
+            "sample_shape",
+            tuple(
+                _validate_positive_integer("sample_shape entry", dimension)
+                for dimension in self.sample_shape
+            ),
+        )
+        if not self.sample_shape:
+            raise ValueError("sample_shape must contain at least one dimension.")
+
+        if type(self.data_digest) is not str:
+            raise TypeError("data_digest must be an exact built-in string.")
+        if re.fullmatch(r"[0-9a-f]{64}", self.data_digest) is None:
+            raise ValueError(
+                "data_digest must be a lowercase 64-character SHA-256 hex digest."
+            )
+
+        if not isinstance(self.partition_sizes, Mapping):
+            raise TypeError("partition_sizes must be a mapping.")
+        if set(self.partition_sizes) != set(_PARTITION_NAMES):
+            raise ValueError(
+                f"partition_sizes must define exactly {list(_PARTITION_NAMES)}."
+            )
+        sizes: dict[str, int] = {}
+        for partition in _PARTITION_NAMES:
+            value = self.partition_sizes[partition]
+            if type(value) is not int or isinstance(value, bool):
+                raise TypeError(
+                    f"partition_sizes[{partition!r}] must be an exact integer."
+                )
+            if value < 0:
+                raise ValueError(
+                    f"partition_sizes[{partition!r}] must not be negative."
+                )
+            sizes[partition] = value
+        for partition in ("train", "validation", "test"):
+            if sizes[partition] < 1:
+                raise ValueError(
+                    f"partition_sizes[{partition!r}] must contain at least one "
+                    "sample; a benchmark cannot fit or score an empty partition."
+                )
+        total = sum(sizes.values())
+        if total != self.n_samples:
+            raise ValueError(
+                f"partition_sizes sum to {total}, but the report describes "
+                f"{self.n_samples} samples. Every sample must be classified as "
+                "train, validation, test, or excluded."
+            )
+        object.__setattr__(self, "partition_sizes", _freeze_json(sizes))
+
+        object.__setattr__(self, "metrics", _validate_metrics(self.metrics))
+
+        if self.seed is not None:
+            if type(self.seed) is not int or isinstance(self.seed, bool):
+                raise TypeError("seed must be an exact non-Boolean integer or None.")
+            if self.seed < 0:
+                raise ValueError("seed must be non-negative.")
+
+        split_identity = _validate_configuration(
+            self.split_identity,
+            name="split_identity",
+        )
+        if "n_samples" not in split_identity:
+            raise ValueError("split_identity must record n_samples.")
+        if split_identity["n_samples"] != self.n_samples:
+            raise ValueError(
+                f"split_identity records {split_identity['n_samples']!r} samples, "
+                f"but the report describes {self.n_samples}."
+            )
+        for field_name in ("time_axis_fingerprint", "partition_digest"):
+            value = split_identity.get(field_name)
+            if value is not None and (
+                type(value) is not str or re.fullmatch(r"[0-9a-f]{64}", value) is None
+            ):
+                raise ValueError(
+                    f"split_identity[{field_name!r}] must be a lowercase "
+                    "64-character SHA-256 hex digest."
+                )
+        object.__setattr__(self, "split_identity", _freeze_json(split_identity))
+
+        if isinstance(self.results, (str, bytes)) or not isinstance(
+            self.results, Sequence
+        ):
+            raise TypeError("results must be a sequence of MethodBenchmarkResult.")
+        results = tuple(self.results)
+        if not results:
+            raise ValueError("results must contain at least one method result.")
+        if any(not isinstance(result, MethodBenchmarkResult) for result in results):
+            raise TypeError("Every result must be a MethodBenchmarkResult.")
+        names = [result.name for result in results]
+        if len(set(names)) != len(names):
+            raise ValueError("Benchmark method names must be unique within a report.")
+        expected_scalars = int(np.prod(self.sample_shape))
+        for result in results:
+            if set(result.test_metrics) != set(self.metrics):
+                raise ValueError(
+                    f"Result {result.name!r} does not report exactly the "
+                    f"requested metrics {list(self.metrics)}."
+                )
+            if result.original_scalars_per_sample != expected_scalars:
+                raise ValueError(
+                    f"Result {result.name!r} records "
+                    f"{result.original_scalars_per_sample} scalars per sample, "
+                    f"but sample_shape {list(self.sample_shape)} contains "
+                    f"{expected_scalars}."
+                )
+        object.__setattr__(self, "results", results)
+
+    def to_dict(self) -> dict[str, JsonValue]:
+        """Return a JSON-compatible dictionary describing the complete run."""
+        return {
+            "schema_version": self.schema_version,
+            "n_samples": self.n_samples,
+            "sample_shape": list(self.sample_shape),
+            "data_digest": self.data_digest,
+            "partition_sizes": _thaw_json(self.partition_sizes),
+            "metrics": list(self.metrics),
+            "seed": self.seed,
+            "split_identity": _thaw_json(self.split_identity),
+            "results": [result.to_dict() for result in self.results],
+        }
+
+    def to_json(self, *, indent: int | None = 2) -> str:
+        """Serialize the complete run deterministically as strict JSON."""
+        return _canonical_json(self.to_dict(), indent=indent) + "\n"
+
+    def identity(self) -> dict[str, JsonValue]:
+        """Return the deterministic identity of the benchmark configuration.
+
+        The identity answers "what was compared, on which samples". It
+        deliberately excludes measured outcomes: metric values and wall-clock
+        timings are observational and are not reproducible bit for bit across
+        machines, library versions, or devices.
+        """
+        return {
+            "schema_version": self.schema_version,
+            "n_samples": self.n_samples,
+            "sample_shape": list(self.sample_shape),
+            "data_digest": self.data_digest,
+            "partition_sizes": _thaw_json(self.partition_sizes),
+            "metrics": list(self.metrics),
+            "seed": self.seed,
+            "split_identity": _thaw_json(self.split_identity),
+            "methods": [result.identity() for result in self.results],
+        }
+
+    def identity_digest(self) -> str:
+        """Return a SHA-256 digest of the deterministic benchmark identity."""
+        payload = _canonical_json(self.identity()).encode("utf-8")
+        return hashlib.sha256(payload).hexdigest()
+
+    @classmethod
+    def from_dict(cls, payload: Mapping[str, Any]) -> ReconstructionBenchmarkReport:
+        """Construct a validated report from a dictionary."""
+        if not isinstance(payload, Mapping):
+            raise TypeError("Report payload must be a mapping.")
+        required = {
+            "schema_version",
+            "n_samples",
+            "sample_shape",
+            "data_digest",
+            "partition_sizes",
+            "metrics",
+            "seed",
+            "split_identity",
+            "results",
+        }
+        _require_exact_fields(payload, required=required, label="Report")
+        for name in ("sample_shape", "metrics", "results"):
+            if type(payload[name]) is not list:
+                raise TypeError(f"Report {name} must be an exact JSON list.")
+        for name in ("partition_sizes", "split_identity"):
+            if type(payload[name]) is not dict:
+                raise TypeError(f"Report {name} must be an exact JSON object.")
+        return cls(
+            schema_version=payload["schema_version"],
+            n_samples=payload["n_samples"],
+            sample_shape=tuple(payload["sample_shape"]),
+            data_digest=payload["data_digest"],
+            partition_sizes=payload["partition_sizes"],
+            metrics=tuple(payload["metrics"]),
+            seed=payload["seed"],
+            split_identity=payload["split_identity"],
+            results=tuple(
+                MethodBenchmarkResult.from_dict(result) for result in payload["results"]
+            ),
+        )
+
+
+def _optional_torch() -> Any | None:
+    try:
+        import torch
+    except ImportError:  # pragma: no cover - depends on installation
+        return None
+    return torch
+
+
+def _forkable_torch_devices(torch: Any) -> list[int]:
+    """Return every CUDA device whose generator benchmark seeding would change.
+
+    ``torch.manual_seed`` seeds all visible CUDA devices, so forking only the
+    current device would leave the other devices permanently reseeded.
+    ``torch.cuda.is_available`` does not initialize CUDA, so CPU-only
+    environments return an empty list without touching the driver.
+    """
+    if not torch.cuda.is_available():
+        return []
+    return list(range(torch.cuda.device_count()))
+
+
+@contextmanager
+def _isolated_random_state(seed: int | None) -> Iterator[None]:
+    """Run a block with an isolated, optionally seeded random state.
+
+    The caller's global NumPy random state and PyTorch generator states are
+    restored on exit, so benchmarking never perturbs surrounding code.
+
+    ``torch.manual_seed`` seeds every visible CUDA device, not only the current
+    one, so every visible device is forked and restored. When CUDA is
+    unavailable no device is touched, which keeps CPU-only environments free of
+    any CUDA initialization.
+
+    Seeding makes a run repeatable on the same machine, device, and library
+    versions. It does not guarantee bitwise-identical PyTorch results across
+    devices, because algorithm selection and reduction order may differ.
+    """
+    numpy_state = np.random.get_state()
+    with ExitStack() as stack:
+        torch = _optional_torch()
+        if torch is not None:
+            stack.enter_context(
+                torch.random.fork_rng(devices=_forkable_torch_devices(torch))
+            )
+        try:
+            if seed is not None:
+                np.random.seed(seed)
+                if torch is not None:
+                    torch.manual_seed(seed)
+            yield
+        finally:
+            np.random.set_state(numpy_state)
+
+
+def _data_digest(X: np.ndarray) -> str:
+    """Return a layout-independent SHA-256 digest of the benchmarked values.
+
+    The split manifest fingerprints the time axis, never the data values, so a
+    separate digest is needed before two runs can be claimed to have compared
+    the same samples. Normalising to C order first makes the digest independent
+    of whether the caller supplied C-ordered or Fortran-ordered data.
+    """
+    contiguous = np.ascontiguousarray(X)
+    digest = hashlib.sha256()
+    digest.update(str(contiguous.dtype.str).encode("utf-8"))
+    digest.update(str(contiguous.shape).encode("utf-8"))
+    digest.update(memoryview(contiguous).cast("B"))
+    return digest.hexdigest()
+
+
+def _split_identity(
+    split: ChronologicalSplit,
+    *,
+    time_axis_verified: bool,
+) -> dict[str, JsonValue]:
+    """Return a stable identity for the partitions produced by ``split``."""
+    manifest = split.manifest
+    partitions = {
+        "train": list(manifest.train_indices),
+        "validation": list(manifest.validation_indices),
+        "test": list(manifest.test_indices),
+        "excluded": list(manifest.excluded_indices),
+    }
+    digest = hashlib.sha256(_canonical_json(partitions).encode("utf-8")).hexdigest()
+    return {
+        "manifest_schema_version": manifest.schema_version,
+        "method": manifest.method,
+        "n_samples": manifest.n_samples,
+        # Named as the manifest names it. This fingerprints the time axis only;
+        # the report's separate data_digest covers the benchmarked values.
+        "time_axis_fingerprint": manifest.dataset_fingerprint,
+        "time_kind": manifest.time_kind,
+        "axis_mode": manifest.axis_mode,
+        "partition_digest": digest,
+        "time_axis_verified": time_axis_verified,
+    }
+
+
+def _validate_methods(methods: Any) -> tuple[BenchmarkMethod, ...]:
+    if isinstance(methods, (str, bytes)) or not isinstance(methods, Sequence):
+        raise TypeError("methods must be a sequence of BenchmarkMethod objects.")
+    specifications = tuple(methods)
+    if not specifications:
+        raise ValueError("methods must contain at least one BenchmarkMethod.")
+    if any(
+        not isinstance(specification, BenchmarkMethod)
+        for specification in specifications
+    ):
+        raise TypeError("Every entry in methods must be a BenchmarkMethod.")
+    names = [specification.name for specification in specifications]
+    duplicates = sorted({name for name in names if names.count(name) > 1})
+    if duplicates:
+        raise ValueError(
+            f"Benchmark method names must be unique; duplicates: {duplicates}."
+        )
+    return specifications
+
+
+def _build_method(specification: BenchmarkMethod) -> ReconstructionMethod:
+    """Instantiate one method and verify it satisfies the benchmark contract."""
+    instance = specification.factory()
+    if instance is None:
+        raise TypeError(f"The factory for method {specification.name!r} returned None.")
+    for attribute in ("fit", "reconstruct"):
+        if not callable(getattr(instance, attribute, None)):
+            raise TypeError(
+                f"Method {specification.name!r} must expose a callable "
+                f"{attribute}() method."
+            )
+    latent = getattr(instance, "latent_dimension", None)
+    if latent is None:
+        raise TypeError(
+            f"Method {specification.name!r} must expose a latent_dimension."
+        )
+    if _validate_positive_integer("latent_dimension", latent) != (
+        specification.latent_dimension
+    ):
+        raise ValueError(
+            f"Method {specification.name!r} reports latent dimension "
+            f"{int(latent)}, but the specification declares "
+            f"{specification.latent_dimension}."
+        )
+    if getattr(instance, "is_fitted", False):
+        raise ValueError(
+            f"The factory for method {specification.name!r} returned an already "
+            "fitted instance. Factories must return a fresh, unfitted model so "
+            "that state cannot leak between benchmark runs."
+        )
+    return instance
+
+
+def _validate_reconstruction(
+    reconstruction: Any,
+    reference: np.ndarray,
+    *,
+    name: str,
+) -> np.ndarray:
+    """Reject reconstructions NumPy would otherwise broadcast into shape."""
+    if not isinstance(reconstruction, np.ndarray):
+        raise TypeError(
+            f"Method {name!r} must return a NumPy array from reconstruct()."
+        )
+    if tuple(reconstruction.shape) != tuple(reference.shape):
+        raise ValueError(
+            f"Method {name!r} returned reconstruction shape "
+            f"{tuple(reconstruction.shape)}, but the test partition has shape "
+            f"{tuple(reference.shape)}. Broadcasting is never applied."
+        )
+    if not np.issubdtype(reconstruction.dtype, np.number):
+        raise TypeError(f"Method {name!r} must return numeric reconstructions.")
+    if np.issubdtype(reconstruction.dtype, np.complexfloating):
+        raise TypeError(f"Method {name!r} must return real-valued reconstructions.")
+    if not np.isfinite(reconstruction).all():
+        raise ValueError(f"Method {name!r} returned non-finite reconstruction values.")
+    return reconstruction
+
+
+def _test_metrics(
+    y_true: np.ndarray,
+    y_pred: np.ndarray,
+    metrics: tuple[str, ...],
+) -> dict[str, float]:
+    reconstruction_error = _load_reconstruction_error()
+    values: dict[str, float] = {}
+    for metric in metrics:
+        values[metric] = float(
+            reconstruction_error(
+                y_true,
+                y_pred,
+                metric=metric,
+                reduction=_METRIC_REDUCTION,
+            )
+        )
+    return values
+
+
+def run_reconstruction_benchmark(
+    X: np.ndarray,
+    *,
+    split: ChronologicalSplit,
+    methods: Sequence[BenchmarkMethod],
+    metrics: Sequence[str] = _DEFAULT_METRICS,
+    seed: int | None = None,
+    sample_times: Any | None = None,
+    sample_start_times: Any | None = None,
+    sample_end_times: Any | None = None,
+) -> ReconstructionBenchmarkReport:
+    """Compare reconstruction methods on identical chronological partitions.
+
+    Every method is fitted on ``X[split.train_indices]``. Methods that declare
+    ``uses_validation_partition`` additionally receive
+    ``X[split.validation_indices]`` for early stopping. Metrics are computed
+    only on ``X[split.test_indices]``, in the original sample space, using the
+    shared implementation in :mod:`bluemath_tk.deeplearning.metrics`.
+
+    Parameters
+    ----------
+    X : np.ndarray
+        Dataset of shape ``(n_samples, ...)``. It is never modified, and every
+        partition handed to a method is an independent copy.
+    split : ChronologicalSplit
+        Manifest-backed chronological partitions. The benchmark never creates a
+        split of its own.
+    methods : sequence of BenchmarkMethod
+        Method specifications with unique names.
+    metrics : sequence of str, optional
+        Metric names from ``("mse", "mae", "rmse")``. Default is
+        ``("mse", "mae", "rmse")``.
+    seed : int, optional
+        Seed applied inside an isolated random state before each method is
+        built and fitted. The caller's random state is restored afterwards.
+    sample_times, sample_start_times, sample_end_times : array-like, optional
+        The time coordinates that order ``X``. When supplied, the split
+        manifest is validated against them, which proves the manifest is being
+        replayed against the dataset it was created from rather than a
+        reordered or different dataset.
+
+    Returns
+    -------
+    ReconstructionBenchmarkReport
+        The complete record of the run.
+
+    Raises
+    ------
+    ValueError
+        If the data, split, methods, metrics, or any reconstruction violates
+        the benchmark contract.
+
+    Notes
+    -----
+    Domain-specific normalisation is the caller's responsibility in this first
+    framework release. Any preprocessing must be applied identically to every
+    compared method and fitted on the training partition alone.
+
+    The reported ``latent_dimensionality_ratio`` compares latent scalars with
+    input scalars. It is not a bitrate and not a storage compression ratio.
+    """
+    _validate_sample_data(X, name="X")
+    if not isinstance(split, ChronologicalSplit):
+        raise TypeError(
+            "split must be a bluemath_tk.validation.ChronologicalSplit, so that "
+            "partition membership is always backed by a validated manifest."
+        )
+    manifest = split.manifest
+    if manifest.n_samples != int(X.shape[0]):
+        raise ValueError(
+            f"The split describes {manifest.n_samples} samples, but X contains "
+            f"{int(X.shape[0])}."
+        )
+
+    time_axis_verified = any(
+        coordinates is not None
+        for coordinates in (sample_times, sample_start_times, sample_end_times)
+    )
+    if time_axis_verified:
+        manifest.validate_against(
+            sample_times=sample_times,
+            sample_start_times=sample_start_times,
+            sample_end_times=sample_end_times,
+        )
+
+    requested_metrics = _validate_metrics(metrics)
+    specifications = _validate_methods(methods)
+    if seed is not None:
+        if type(seed) is not int or isinstance(seed, bool):
+            raise TypeError("seed must be an exact non-Boolean integer or None.")
+        if not 0 <= seed < 2**32:
+            # Checked here rather than left to NumPy, so an out-of-range seed
+            # fails before any method is built instead of part way through a run.
+            raise ValueError(
+                f"seed must lie in [0, 2**32) to be usable as a NumPy seed; got {seed}."
+            )
+
+    train_indices = np.asarray(split.train_indices)
+    validation_indices = np.asarray(split.validation_indices)
+    test_indices = np.asarray(split.test_indices)
+    for name, indices in (
+        ("train_indices", train_indices),
+        ("validation_indices", validation_indices),
+        ("test_indices", test_indices),
+    ):
+        if indices.size == 0:
+            raise ValueError(f"split.{name} must not be empty.")
+        if int(indices.max()) >= int(X.shape[0]) or int(indices.min()) < 0:
+            raise ValueError(f"split.{name} contains an index outside X.")
+
+    # Fancy indexing copies, so no method can reach or mutate the caller's X.
+    # These references stay pristine: each method receives its own copies, and
+    # metrics are always computed against this untouched test partition.
+    X_train = X[train_indices]
+    X_validation = X[validation_indices]
+    X_test = X[test_indices]
+
+    sample_shape = tuple(int(dimension) for dimension in X.shape[1:])
+    original_scalars = int(np.prod(sample_shape))
+
+    built: list[ReconstructionMethod] = []
+    results: list[MethodBenchmarkResult] = []
+    for specification in specifications:
+        with _isolated_random_state(seed):
+            instance = _build_method(specification)
+            if any(instance is other for other in built):
+                raise ValueError(
+                    f"The factory for method {specification.name!r} returned an "
+                    "instance already used by another method. Every method must "
+                    "get its own model."
+                )
+            built.append(instance)
+
+            # Every method gets its own copies, so a method that writes to its
+            # inputs cannot corrupt the partitions seen by later methods.
+            method_train = X_train.copy()
+            method_validation = (
+                X_validation.copy() if specification.uses_validation_partition else None
+            )
+            method_test = X_test.copy()
+
+            fit_start = perf_counter()
+            instance.fit(method_train, method_validation)
+            fit_seconds = perf_counter() - fit_start
+
+            reconstruction_start = perf_counter()
+            reconstruction = instance.reconstruct(method_test)
+            reconstruction_seconds = perf_counter() - reconstruction_start
+
+        reconstruction = _validate_reconstruction(
+            reconstruction,
+            X_test,
+            name=specification.name,
+        )
+        results.append(
+            MethodBenchmarkResult(
+                name=specification.name,
+                method_type=specification.method_type,
+                latent_dimension=specification.latent_dimension,
+                uses_validation_partition=specification.uses_validation_partition,
+                original_scalars_per_sample=original_scalars,
+                latent_scalars_per_sample=specification.latent_dimension,
+                latent_dimensionality_ratio=(
+                    specification.latent_dimension / original_scalars
+                ),
+                test_metrics=_test_metrics(
+                    X_test,
+                    reconstruction,
+                    requested_metrics,
+                ),
+                fit_seconds=max(fit_seconds, 0.0),
+                reconstruction_seconds=max(reconstruction_seconds, 0.0),
+                configuration=_thaw_json(specification.configuration),
+            )
+        )
+
+    counts = split.counts
+    return ReconstructionBenchmarkReport(
+        schema_version=_SCHEMA_VERSION,
+        n_samples=int(X.shape[0]),
+        sample_shape=sample_shape,
+        data_digest=_data_digest(X),
+        partition_sizes={name: int(counts[name]) for name in _PARTITION_NAMES},
+        metrics=requested_metrics,
+        seed=seed,
+        split_identity=_split_identity(split, time_axis_verified=time_axis_verified),
+        results=tuple(results),
+    )
diff --git a/bluemath_tk/deeplearning/_base_model.py b/bluemath_tk/deeplearning/_base_model.py
index 921c4dc..65cafed 100644
--- a/bluemath_tk/deeplearning/_base_model.py
+++ b/bluemath_tk/deeplearning/_base_model.py
@@ -235,6 +235,7 @@ def _validate_fit_inputs(
         batch_size: int,
         epochs: int,
         patience: int,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None = None,
     ) -> None:
         """Validate common training inputs before building the model."""
         if not isinstance(X, np.ndarray):
@@ -257,7 +258,7 @@ def _validate_fit_inputs(
         self._validate_finite_array(X, "X")
         self._validate_finite_array(y, "y")
 
-        if (
+        if validation_data is None and (
             not isinstance(validation_split, (int, float))
             or isinstance(validation_split, bool)
             or not np.isfinite(float(validation_split))
@@ -274,6 +275,17 @@ def _validate_fit_inputs(
             if not isinstance(value, int) or isinstance(value, bool) or value < 1:
                 raise ValueError(f"{name} must be a positive integer.")
 
+        if validation_data is not None:
+            # The pair itself is validated once, in _resolve_fit_partitions,
+            # which runs immediately after this method and before the model is
+            # built.
+            if len(X) < 2:
+                raise ValueError(
+                    "Explicit validation_data requires at least two training "
+                    "samples in X."
+                )
+            return
+
         split = int((1 - validation_split) * len(X))
         if split < 2:
             raise ValueError(
@@ -283,6 +295,97 @@ def _validate_fit_inputs(
         if len(X) - split < 1:
             raise ValueError("The validation split must contain at least one sample.")
 
+    def _validate_validation_data(
+        self,
+        X: np.ndarray,
+        validation_data: tuple[np.ndarray, np.ndarray | None],
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """Validate explicit validation data and return the resolved arrays.
+
+        Parameters
+        ----------
+        X : np.ndarray
+            The training inputs, used to check the per-sample contract.
+        validation_data : tuple
+            An ``(X_validation, y_validation)`` pair. ``y_validation`` may be
+            ``None``, in which case the model's default reconstruction target
+            is derived from ``X_validation``.
+
+        Returns
+        -------
+        tuple of np.ndarray
+            The validated ``(X_validation, y_validation)`` arrays.
+        """
+        if not isinstance(validation_data, tuple) or len(validation_data) != 2:
+            raise TypeError(
+                "validation_data must be an (X_validation, y_validation) tuple; "
+                "pass y_validation=None to reconstruct X_validation itself."
+            )
+        X_validation, y_validation = validation_data
+        if not isinstance(X_validation, np.ndarray):
+            raise TypeError("validation_data[0] must be a NumPy array.")
+        if X_validation.ndim != X.ndim:
+            raise ValueError(
+                f"validation_data[0] must have {X.ndim} dimensions to match X; "
+                f"got {X_validation.ndim}."
+            )
+        if len(X_validation) < 1:
+            raise ValueError("validation_data[0] must contain at least one sample.")
+        if tuple(X_validation.shape[1:]) != tuple(X.shape[1:]):
+            raise ValueError(
+                "validation_data[0] per-sample shape "
+                f"{tuple(X_validation.shape[1:])} does not match the training "
+                f"per-sample shape {tuple(X.shape[1:])}."
+            )
+        self._validate_finite_array(X_validation, "validation_data[0]")
+
+        if y_validation is None:
+            y_validation = self._get_reconstruction_target(X_validation)
+        if not isinstance(y_validation, np.ndarray):
+            raise TypeError("validation_data[1] must be a NumPy array or None.")
+        if len(y_validation) != len(X_validation):
+            raise ValueError(
+                "validation_data arrays must contain the same number of samples; "
+                f"got {len(X_validation)} and {len(y_validation)}."
+            )
+        self._validate_target_shape(X_validation, y_validation)
+        self._validate_finite_array(y_validation, "validation_data[1]")
+        return X_validation, y_validation
+
+    def _resolve_fit_partitions(
+        self,
+        X: np.ndarray,
+        y: np.ndarray,
+        validation_split: float,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None,
+    ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+        """Return the training and validation arrays used by one fit call.
+
+        With ``validation_data=None`` the historical behaviour is preserved
+        exactly: one random permutation of ``X`` is cut at ``validation_split``
+        using the global NumPy random state.
+
+        When explicit ``validation_data`` is supplied, ``validation_split`` is
+        ignored, every sample of ``X`` is used for optimisation in the given
+        order, and the global NumPy random state is left untouched. This makes
+        the validation membership exactly reproducible, which is required for
+        chronological validation.
+        """
+        if validation_data is None:
+            indices = np.arange(len(X))
+            np.random.shuffle(indices)
+            split = int((1 - validation_split) * len(X))
+            train_indices, validation_indices = indices[:split], indices[split:]
+            return (
+                X[train_indices],
+                y[train_indices],
+                X[validation_indices],
+                y[validation_indices],
+            )
+
+        X_validation, y_validation = self._validate_validation_data(X, validation_data)
+        return X, y, X_validation, y_validation
+
     def _get_init_config(self) -> dict:
         """Collect constructor parameters needed to recreate this model."""
         config = {}
@@ -536,9 +639,21 @@ def fit(
         criterion: nn.Module | None = None,
         patience: int = 20,
         verbose: int = 1,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None = None,
         **kwargs,
     ) -> dict[str, list]:
-        """Fit a reconstruction model with finite, sample-weighted losses."""
+        """Fit a reconstruction model with finite, sample-weighted losses.
+
+        Parameters
+        ----------
+        validation_data : tuple, optional
+            An explicit ``(X_validation, y_validation)`` pair. When supplied,
+            ``validation_split`` is ignored, all of ``X`` is used for
+            optimisation, and exactly these samples drive the validation loss
+            and early stopping. ``y_validation`` may be ``None`` to reconstruct
+            ``X_validation`` itself. Default is None, which keeps the historical
+            random ``validation_split`` behaviour.
+        """
         learning_rate = self._validate_learning_rate(learning_rate)
         if not isinstance(X, np.ndarray):
             raise TypeError("X must be a NumPy array.")
@@ -552,7 +667,14 @@ def fit(
             batch_size,
             epochs,
             patience,
+            validation_data=validation_data,
         )
+        (
+            X_train_array,
+            y_train_array,
+            X_validation_array,
+            y_validation_array,
+        ) = self._resolve_fit_partitions(X, y, validation_split, validation_data)
         self._validate_or_set_build_input_shape(tuple(X.shape))
         self.is_fitted = False
 
@@ -566,22 +688,17 @@ def fit(
         if criterion is None:
             criterion = nn.MSELoss()
 
-        indices = np.arange(len(X))
-        np.random.shuffle(indices)
-        split = int((1 - validation_split) * len(X))
-        train_indices, validation_indices = indices[:split], indices[split:]
-
         X_train = torch.as_tensor(
-            X[train_indices], dtype=torch.float32, device=self.device
+            X_train_array, dtype=torch.float32, device=self.device
         )
         y_train = torch.as_tensor(
-            y[train_indices], dtype=torch.float32, device=self.device
+            y_train_array, dtype=torch.float32, device=self.device
         )
         X_validation = torch.as_tensor(
-            X[validation_indices], dtype=torch.float32, device=self.device
+            X_validation_array, dtype=torch.float32, device=self.device
         )
         y_validation = torch.as_tensor(
-            y[validation_indices], dtype=torch.float32, device=self.device
+            y_validation_array, dtype=torch.float32, device=self.device
         )
 
         history = {"train_loss": [], "val_loss": []}
diff --git a/bluemath_tk/deeplearning/autoencoders.py b/bluemath_tk/deeplearning/autoencoders.py
index 2705a15..4640c3e 100644
--- a/bluemath_tk/deeplearning/autoencoders.py
+++ b/bluemath_tk/deeplearning/autoencoders.py
@@ -398,9 +398,18 @@ def fit(
         criterion: Optional[nn.Module] = None,
         patience: int = 20,
         verbose: int = 1,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None = None,
         **kwargs,
     ) -> Dict[str, list]:
-        """Fit with orthogonality and latent-decorrelation penalties."""
+        """Fit with orthogonality and latent-decorrelation penalties.
+
+        Parameters
+        ----------
+        validation_data : tuple, optional
+            An explicit ``(X_validation, y_validation)`` pair. When supplied,
+            ``validation_split`` is ignored and exactly these samples drive the
+            validation objective and early stopping. Default is None.
+        """
         learning_rate = self._validate_learning_rate(learning_rate)
         if not isinstance(X, np.ndarray):
             raise TypeError("X must be a NumPy array.")
@@ -413,7 +422,14 @@ def fit(
             batch_size,
             epochs,
             patience,
+            validation_data=validation_data,
         )
+        (
+            X_train_array,
+            y_train_array,
+            X_validation_array,
+            y_validation_array,
+        ) = self._resolve_fit_partitions(X, y, validation_split, validation_data)
         self._validate_or_set_build_input_shape(tuple(X.shape))
         self.is_fitted = False
 
@@ -426,21 +442,17 @@ def fit(
         if criterion is None:
             criterion = nn.MSELoss()
 
-        indices = np.arange(len(X))
-        np.random.shuffle(indices)
-        split = int((1 - validation_split) * len(X))
-        train_indices, validation_indices = indices[:split], indices[split:]
         X_train = torch.as_tensor(
-            X[train_indices], dtype=torch.float32, device=self.device
+            X_train_array, dtype=torch.float32, device=self.device
         )
         y_train = torch.as_tensor(
-            y[train_indices], dtype=torch.float32, device=self.device
+            y_train_array, dtype=torch.float32, device=self.device
         )
         X_validation = torch.as_tensor(
-            X[validation_indices], dtype=torch.float32, device=self.device
+            X_validation_array, dtype=torch.float32, device=self.device
         )
         y_validation = torch.as_tensor(
-            y[validation_indices], dtype=torch.float32, device=self.device
+            y_validation_array, dtype=torch.float32, device=self.device
         )
 
         history = {"train_loss": [], "val_loss": []}
@@ -1210,6 +1222,7 @@ def fit(
         criterion: nn.Module | None = None,
         patience: int = 20,
         verbose: int = 1,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None = None,
         **kwargs,
     ) -> dict[str, list]:
         """Fit the model to reconstruct the complete input sequence."""
@@ -1228,6 +1241,7 @@ def fit(
             criterion=criterion,
             patience=patience,
             verbose=verbose,
+            validation_data=validation_data,
             **kwargs,
         )
 
@@ -1491,6 +1505,7 @@ def fit(
         criterion: nn.Module | None = None,
         patience: int = 20,
         verbose: int = 1,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None = None,
         **kwargs,
     ) -> dict[str, list]:
         """Fit the model to reconstruct the complete input sequence."""
@@ -1509,6 +1524,7 @@ def fit(
             criterion=criterion,
             patience=patience,
             verbose=verbose,
+            validation_data=validation_data,
             **kwargs,
         )
 
diff --git a/bluemath_tk/deeplearning/spatiotemporal_autoencoders.py b/bluemath_tk/deeplearning/spatiotemporal_autoencoders.py
index 2e3492b..a0a9e0e 100644
--- a/bluemath_tk/deeplearning/spatiotemporal_autoencoders.py
+++ b/bluemath_tk/deeplearning/spatiotemporal_autoencoders.py
@@ -154,6 +154,7 @@ def fit(
         criterion: nn.Module | None = None,
         patience: int = 20,
         verbose: int = 1,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None = None,
         **kwargs,
     ) -> dict[str, list]:
         """Fit the model to reconstruct the complete input sequence."""
@@ -171,6 +172,7 @@ def fit(
             criterion=criterion,
             patience=patience,
             verbose=verbose,
+            validation_data=validation_data,
             **kwargs,
         )
 
diff --git a/bluemath_tk/deeplearning/variational_autoencoders.py b/bluemath_tk/deeplearning/variational_autoencoders.py
index 964467f..b58bfa4 100644
--- a/bluemath_tk/deeplearning/variational_autoencoders.py
+++ b/bluemath_tk/deeplearning/variational_autoencoders.py
@@ -223,6 +223,7 @@ def fit(
         criterion: nn.Module | None = None,
         patience: int = 20,
         verbose: int = 1,
+        validation_data: tuple[np.ndarray, np.ndarray | None] | None = None,
         **kwargs,
     ) -> dict[str, list]:
         """Fit the VAE with stochastic train and validation objectives.
@@ -231,6 +232,13 @@ def fit(
         used for training and controls early stopping. The separate
         ``val_deterministic_reconstruction_loss`` reports posterior-mean
         reconstruction for stable scientific comparison.
+
+        Parameters
+        ----------
+        validation_data : tuple, optional
+            An explicit ``(X_validation, y_validation)`` pair. When supplied,
+            ``validation_split`` is ignored and exactly these samples drive the
+            validation objective and early stopping. Default is None.
         """
         learning_rate = self._validate_learning_rate(learning_rate)
         if not isinstance(X, np.ndarray):
@@ -245,7 +253,14 @@ def fit(
             batch_size,
             epochs,
             patience,
+            validation_data=validation_data,
         )
+        (
+            X_train_array,
+            y_train_array,
+            X_validation_array,
+            y_validation_array,
+        ) = self._resolve_fit_partitions(X, y, validation_split, validation_data)
         self._validate_or_set_build_input_shape(tuple(X.shape))
         self.is_fitted = False
 
@@ -266,29 +281,23 @@ def fit(
                 "reconstruction criterion."
             )
 
-        indices = np.arange(len(X))
-        np.random.shuffle(indices)
-        split = int((1 - validation_split) * len(X))
-        train_indices = indices[:split]
-        validation_indices = indices[split:]
-
         X_train = torch.as_tensor(
-            X[train_indices],
+            X_train_array,
             dtype=torch.float32,
             device=self.device,
         )
         y_train = torch.as_tensor(
-            y[train_indices],
+            y_train_array,
             dtype=torch.float32,
             device=self.device,
         )
         X_validation = torch.as_tensor(
-            X[validation_indices],
+            X_validation_array,
             dtype=torch.float32,
             device=self.device,
         )
         y_validation = torch.as_tensor(
-            y[validation_indices],
+            y_validation_array,
             dtype=torch.float32,
             device=self.device,
         )
diff --git a/docs/source/benchmarking.rst b/docs/source/benchmarking.rst
new file mode 100644
index 0000000..b3d6a39
--- /dev/null
+++ b/docs/source/benchmarking.rst
@@ -0,0 +1,379 @@
+PCA and autoencoder reconstruction benchmarking
+===============================================
+
+Comparing a PCA baseline against an autoencoder is only meaningful when both
+methods are fitted on the same training samples and scored on the same held-out
+samples. ``bluemath_tk.benchmarking`` provides that shared infrastructure so
+that reconstruction results are comparable rather than accidental.
+
+What the framework measures
+---------------------------
+
+The framework measures **reconstruction error on a held-out test partition**,
+in the original sample space, for any number of dimensionality-reduction
+methods that share one latent budget.
+
+It deliberately does not measure downstream scientific skill. A method with the
+lowest reconstruction error is not automatically the best choice for a physical
+diagnostic, an extreme-value analysis, or a forecast. The report therefore
+contains no ranking and no "best method" field.
+
+Why chronological partitions matter
+-----------------------------------
+
+Random splits of a time series place neighbouring, highly correlated samples on
+both sides of the split, so a model can reach a low error by memorising almost
+identical neighbours. Partition membership therefore always comes from
+:func:`~bluemath_tk.validation.chronological.split_chronologically` or from a
+replayed :class:`~bluemath_tk.validation.chronological.ValidationSplitManifest`.
+
+The benchmark never creates a split of its own, and it refuses to run without a
+manifest-backed :class:`~bluemath_tk.validation.chronological.ChronologicalSplit`.
+
+Basic use
+---------
+
+.. code-block:: python
+
+   import numpy as np
+   from bluemath_tk.benchmarking import (
+       autoencoder_benchmark_method,
+       pca_benchmark_method,
+       run_reconstruction_benchmark,
+   )
+   from bluemath_tk.deeplearning.autoencoders import StandardAutoencoder
+   from bluemath_tk.validation import split_chronologically
+
+   rng = np.random.default_rng(0)
+   times = np.arange("2000-01", "2005-01", dtype="datetime64[M]")
+   latent = rng.normal(size=(times.size, 3))
+   mixing = rng.normal(size=(3, 12))
+   X = (latent @ mixing).reshape(times.size, 3, 4)
+
+   split = split_chronologically(sample_times=times, fractions=(0.6, 0.2, 0.2))
+
+   report = run_reconstruction_benchmark(
+       X,
+       split=split,
+       methods=[
+           pca_benchmark_method("pca-k3", n_components=3),
+           autoencoder_benchmark_method(
+               "standard-ae-k3",
+               model_factory=lambda: StandardAutoencoder(k=3, hidden_dims=[32]),
+               latent_dimension=3,
+               configuration={
+                   "architecture": "StandardAutoencoder",
+                   "hidden_dims": [32],
+               },
+               fit_kwargs={"epochs": 50, "batch_size": 16},
+           ),
+       ],
+       metrics=("mse", "mae", "rmse"),
+       seed=0,
+       sample_times=times,
+   )
+
+   for result in report.results:
+       print(result.name, result.test_metrics)
+
+Passing ``sample_times`` is optional but recommended. When supplied, the split
+manifest is validated against those coordinates, which proves the manifest is
+being replayed against the dataset it was created from rather than against a
+reordered or different dataset.
+
+The recorded configuration of the autoencoder above describes the **effective**
+training run, not only what the caller happened to spell out. Defaults that were
+never passed are filled in, so the serialized identity is complete:
+
+.. code-block:: python
+
+   >>> report.results[1].configuration["fit_kwargs"]
+   {'batch_size': 16, 'epochs': 50, 'learning_rate': 0.001, 'patience': 20}
+   >>> report.results[1].configuration["predict_kwargs"]
+   {'batch_size': 64}
+
+Changing ``epochs``, ``learning_rate``, ``batch_size``, or ``patience``
+therefore changes ``identity_digest()``. Nothing has to be duplicated by hand
+inside ``configuration``.
+
+How PCA and autoencoders are compared
+-------------------------------------
+
+Every method is presented to the runner through the same small interface: it is
+fitted, then asked to reconstruct samples. Adding a new reconstruction model
+later therefore does not change the scientific core of the benchmark.
+
+``PCAReconstruction`` wraps the existing
+:class:`bluemath_tk.datamining.pca.PCA` implementation. Samples of shape
+``(n_samples, d1, ..., dm)`` are presented as one stacked variable. Stacking and
+the inverse reshape both use C order, so sample order and the per-sample shape
+survive the round trip unchanged. Impossible component counts are rejected with
+an explicit message before scikit-learn is reached.
+
+``AutoencoderReconstruction`` wraps a BlueMath autoencoder and uses only its
+accepted public workflow: ``fit`` with explicit chronological validation data
+and ``predict`` for reconstruction. Model architectures are never modified.
+
+The two families are inherently different procedures, so a few asymmetries are
+unavoidable. Each one is recorded rather than hidden:
+
+* PCA has no early stopping and no validation-driven model selection, so it is
+  fitted on the training partition alone and never sees the validation
+  partition. Autoencoders receive the validation partition for early stopping
+  only. This is recorded per method as ``uses_validation_partition``.
+* PCA has a closed-form solution, while autoencoders are fitted by mini-batch
+  gradient descent. Training samples are therefore permuted once before an
+  autoencoder is fitted, so that mini-batches are not contiguous blocks of
+  adjacent timestamps, which would otherwise make batch statistics and gradient
+  estimates reflect temporally correlated neighbours. The permutation is drawn
+  inside the benchmark's isolated random state, so it is controlled by the run
+  ``seed`` and never touches the caller's random state. It reorders training
+  rows only: partition membership, and in particular validation membership, is
+  unchanged. This is recorded per method as ``shuffle_training_data`` and can be
+  disabled with ``shuffle_training_data=False``.
+
+No other asymmetry is introduced. Both families see the same training samples,
+neither sees the test partition, and both are scored by the same metric code on
+the same test samples.
+
+Supplying a model
+-----------------
+
+Methods are supplied as specifications, not as instances. Each specification
+carries a **factory** that returns a fresh, unfitted model. A new instance is
+built for every run, so fitted state cannot leak between runs, and the runner
+rejects a factory that returns an already fitted or already used instance.
+
+Factories are never introspected. Because an arbitrary Python callable cannot be
+serialized reproducibly, the reproducible description of a method comes from the
+explicit, user-supplied ``method_type`` and ``configuration`` fields, which must
+be JSON-compatible and are recorded verbatim.
+
+A method may also be written from scratch. Anything exposing
+``latent_dimension``, ``is_fitted``, ``fit(X_train, X_validation)``, and
+``reconstruct(X)`` satisfies the ``ReconstructionMethod`` protocol and can be
+wrapped in a :class:`~bluemath_tk.benchmarking.BenchmarkMethod`.
+
+Two constructor arguments are refused for autoencoders. ``validation_data`` and
+``validation_split`` are refused because the benchmark controls partition
+membership. ``optimizer`` and ``criterion`` are refused because a stateful
+object built for one model instance would be silently reused by the fresh
+instance of the next run; a PyTorch optimizer bound to another model's
+parameters skips them without raising, which would report an untrained network
+as a legitimate result.
+
+A wrapped autoencoder must also declare ``validation_data`` explicitly in its
+``fit`` signature. A ``fit(self, X, **kwargs)`` signature would absorb the
+argument silently and remain free to build its own random validation split,
+which is exactly the leakage this framework exists to prevent, so it is
+rejected.
+
+Everything else in ``fit_kwargs`` and ``predict_kwargs`` must be
+JSON-compatible, and is rejected with a clear error otherwise. This is a
+deliberate constraint rather than an implementation limit: a training setting
+that cannot be recorded cannot form part of a reproducible identity, and
+guessing at a ``repr`` of an arbitrary object would produce an identity that
+looks precise while meaning nothing.
+
+The recorded values are deep copies. Mutating a nested mapping you passed as
+``configuration``, ``fit_kwargs``, or ``predict_kwargs`` after building the
+specification changes neither what later runs execute nor what the report
+records, and each run receives its own containers.
+
+This automatic recording belongs to ``autoencoder_benchmark_method``. If you
+assemble a :class:`~bluemath_tk.benchmarking.BenchmarkMethod` by hand around a
+custom method, the ``configuration`` you supply is the whole of what gets
+recorded, so it must describe everything that defines the experiment.
+
+Common metrics
+--------------
+
+Metrics are computed with the accepted implementation in
+:mod:`bluemath_tk.deeplearning.metrics`. The available metrics are ``mse``,
+``mae``, and ``rmse``, all reported with ``reduction="mean"``.
+
+Benchmark numbers are therefore bit-identical to a direct call to
+``reconstruction_error(y_true, y_pred, metric=..., reduction="mean")``. They
+agree with the per-model ``evaluate_reconstruction`` summaries to floating-point
+rounding rather than bit for bit, because those summarise per-sample errors
+through a different reduction path.
+
+Note that ``rmse`` follows the BlueMath convention: it is the mean over samples
+of the per-sample root-mean-square error, which is not the same quantity as the
+square root of the reported ``mse``.
+
+Latent dimensionality is not storage compression
+------------------------------------------------
+
+Each result reports:
+
+.. code-block:: text
+
+   original_scalars_per_sample
+   latent_scalars_per_sample
+   latent_dimensionality_ratio
+
+``latent_dimensionality_ratio`` is ``latent_scalars_per_sample`` divided by
+``original_scalars_per_sample``. It is a **dimensionality** ratio only.
+
+It is explicitly **not** a bitrate, a storage compression ratio, an entropy
+coding result, or a compressed file size, because it ignores latent numeric
+precision, quantisation, entropy coding, and the storage cost of the model
+parameters themselves. Quantisation and real storage accounting belong to a
+later contribution.
+
+How the test partition is kept isolated
+---------------------------------------
+
+The test partition never reaches any fitting step:
+
+* PCA is fitted on ``X[split.train_indices]`` only, including the optional
+  ``StandardScaler`` statistics when ``scale_data=True``.
+* Autoencoder optimisation receives ``X[split.train_indices]`` only.
+* Autoencoder validation loss and early stopping receive
+  ``X[split.validation_indices]`` only.
+* Metrics are computed on ``X[split.test_indices]`` only.
+
+Explicit validation data
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+BlueMath autoencoders historically derived their validation set by shuffling
+the samples and cutting at ``validation_split``, which cannot express a
+chronological validation partition. ``fit`` therefore accepts an explicit
+``validation_data`` pair:
+
+.. code-block:: python
+
+   model.fit(X_train, validation_data=(X_validation, None))
+
+When ``validation_data`` is supplied, ``validation_split`` is ignored, all of
+``X`` is used for optimisation in the order given, exactly the supplied samples
+drive the validation loss and early stopping, and the global NumPy random state
+is left untouched. Passing ``y_validation=None`` reconstructs ``X_validation``
+itself. The historical ``validation_split`` behaviour is unchanged when
+``validation_data`` is omitted.
+
+Preprocessing
+~~~~~~~~~~~~~
+
+Domain-specific normalisation is the caller's responsibility in this first
+release. The framework applies no shared preprocessing of its own, which
+guarantees that no method receives a transformation the others do not. If you
+normalise, fit the transformation on the training partition alone and apply the
+identical transformation to every compared method.
+
+The one exception is intrinsic to PCA itself: scikit-learn's PCA always centers
+the data internally. The additional ``StandardScaler`` step of the BlueMath PCA
+is exposed as ``scale_data`` and defaults to ``False``, so it is never applied
+silently.
+
+Variational autoencoders
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+For a fair comparison, reconstruction must be deterministic. ``predict`` on
+:class:`~bluemath_tk.deeplearning.variational_autoencoders.VariationalAutoencoder`
+uses the posterior mean by default, and the benchmark adapter rejects
+``stochastic=True`` rather than silently comparing one stochastic draw against
+deterministic PCA and autoencoder reconstructions.
+
+Reproducibility
+---------------
+
+Results are returned as a
+:class:`~bluemath_tk.benchmarking.ReconstructionBenchmarkReport`:
+
+.. code-block:: python
+
+   payload = report.to_dict()
+   text = report.to_json()
+   digest = report.identity_digest()
+
+``to_json`` is strict and deterministic: keys are sorted, NaN and infinity are
+rejected, and no timestamps are written.
+
+``identity()`` and ``identity_digest()`` describe *what was compared*: dataset
+shape, a digest of the data values, partition sizes, split identity, requested
+metrics, seed, and every method specification. They deliberately exclude
+measured outcomes. Metric values and wall-clock timings are observational and
+are not reproducible bit for bit across machines, library versions, or devices,
+so including them in a reproducibility identity would make that identity
+meaningless.
+
+Two fingerprints are recorded, and they cover different things:
+
+* ``split_identity["time_axis_fingerprint"]`` comes from the split manifest and
+  covers the **time coordinates** only.
+* ``data_digest`` covers the **benchmarked values**, their dtype, and their
+  shape. It is normalised to C order, so C-ordered and Fortran-ordered copies of
+  the same data produce the same digest.
+
+Both are needed: identical time coordinates do not imply identical data.
+
+A report rebuilt with ``from_dict`` is re-validated across fields, not only
+field by field, so internally inconsistent scientific metadata is rejected
+rather than silently trusted. The partition sizes must sum to ``n_samples``;
+train, validation, and test must each be non-empty; ``split_identity`` must
+record the same ``n_samples`` as the report and well-formed digests; each
+result's ``original_scalars_per_sample`` must equal the product of
+``sample_shape``; ``latent_dimension`` must equal ``latent_scalars_per_sample``;
+``latent_dimensionality_ratio`` must equal the ratio it claims to be; and no
+reconstruction metric may be negative, since MSE, MAE, and RMSE are
+non-negative by construction.
+
+Timing and random state
+-----------------------
+
+``fit_seconds`` and ``reconstruction_seconds`` come from
+:func:`time.perf_counter`. They are observational measurements of one run on one
+machine and should not be treated as deterministic benchmark outputs.
+
+Each method is built and fitted inside an isolated random state. The caller's
+global NumPy state and PyTorch generator states are restored afterwards. When
+``seed`` is supplied, every method starts from the same seeded state. Seeding
+makes a run repeatable on the same machine, device, and library versions; it
+does not guarantee bitwise-identical PyTorch results across devices, because
+algorithm selection and reduction order may differ.
+
+``torch.manual_seed`` reseeds every visible CUDA device, not only the current
+one, so every visible CUDA device is forked and restored. When CUDA is
+unavailable no device is queried at all, which keeps CPU-only environments free
+of any CUDA initialization.
+
+One limit of that isolation is worth stating precisely: Python's standard
+library ``random`` module is not isolated, so a model that draws from it is
+neither seeded nor restored.
+
+Current limitations
+-------------------
+
+* Reconstruction error only. Downstream scientific evaluation is not included.
+* One fixed chronological split per run. Rolling-origin, expanding-window, and
+  walk-forward validation are not included.
+* No hyperparameter optimisation. Architectures and hyperparameters are supplied
+  by the caller.
+* PCA is benchmarked with an explicit integer component count. The
+  explained-variance-ratio mode of :class:`bluemath_tk.datamining.pca.PCA` is not
+  exposed here, because a fixed latent budget is what makes the comparison
+  against an autoencoder latent dimension meaningful.
+* No latent quantisation, entropy coding, bitrate, or storage-size accounting.
+* No shared preprocessing, and no dataset-specific loaders or downloads.
+* Metrics require the ``deeplearning`` extra, because they are reused from
+  :mod:`bluemath_tk.deeplearning.metrics`.
+* Partitions are materialised as copies, and each method receives its own copies
+  so that a method writing to its inputs cannot affect later methods. Peak
+  memory therefore includes two copies of the train, validation, and test
+  partitions.
+* :class:`bluemath_tk.datamining.pca.PCA` logs at INFO and WARNING level on every
+  fit and every transform, and creates a ``logs/`` directory in the working
+  directory. A benchmark run surfaces that existing behaviour and offers no way
+  to quiet it; ``verbose=0`` applies to autoencoders only.
+* ``PCAReconstruction`` does not pass ``random_state`` to scikit-learn. For large
+  problems ``svd_solver="auto"`` may select randomized SVD, which draws from the
+  global NumPy random state; supply a benchmark ``seed`` if you need that to be
+  reproducible.
+* ``uses_validation_partition`` is a declaration by the specification, not a
+  measurement. It states whether a method is handed the validation partition; it
+  cannot verify what the method then does with it.
+* The runner cannot verify that ``X`` is ordered by the manifest's time axis
+  unless the time coordinates are passed through ``sample_times``,
+  ``sample_start_times``, or ``sample_end_times``.
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 7017bd9..2fe10f2 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -15,6 +15,7 @@ BlueMath-tk: A Python Library for Coastal Climate Hazards
    contribute
    modules
    validation
+   benchmarking
 
 Indices and tables
 ==================
diff --git a/tests/benchmarking/test_reconstruction.py b/tests/benchmarking/test_reconstruction.py
new file mode 100644
index 0000000..6ab759c
--- /dev/null
+++ b/tests/benchmarking/test_reconstruction.py
@@ -0,0 +1,2042 @@
+"""Regression tests for the common reconstruction benchmark framework."""
+
+from __future__ import annotations
+
+import copy
+import hashlib
+import inspect
+import json
+
+import numpy as np
+import pytest
+
+torch = pytest.importorskip("torch")
+
+from bluemath_tk.benchmarking import (  # noqa: E402
+    AutoencoderReconstruction,
+    BenchmarkMethod,
+    MethodBenchmarkResult,
+    PCAReconstruction,
+    ReconstructionBenchmarkReport,
+    autoencoder_benchmark_method,
+    pca_benchmark_method,
+    run_reconstruction_benchmark,
+)
+from bluemath_tk.benchmarking.reconstruction import (  # noqa: E402
+    _BENCHMARK_FIT_DEFAULTS,
+    _BENCHMARK_PREDICT_DEFAULTS,
+    _isolated_random_state,
+)
+from bluemath_tk.deeplearning._base_model import BaseDeepLearningModel  # noqa: E402
+from bluemath_tk.deeplearning.autoencoders import (  # noqa: E402
+    OrthogonalAutoencoder,
+    StandardAutoencoder,
+)
+from bluemath_tk.deeplearning.metrics import reconstruction_error  # noqa: E402
+from bluemath_tk.deeplearning.variational_autoencoders import (  # noqa: E402
+    VariationalAutoencoder,
+)
+from bluemath_tk.validation import split_chronologically  # noqa: E402
+
+METRICS = ("mse", "mae", "rmse")
+
+
+class _JsonModel:
+    """Minimal duck-typed autoencoder used for identity and mutation tests."""
+
+    k = 2
+    is_fitted = False
+
+    def __init__(self):
+        self.fit_calls = []
+        self.predict_calls = []
+
+    def fit(self, X_train, validation_data=None, **kwargs):
+        self.fit_calls.append(copy.deepcopy(kwargs))
+        self.is_fitted = True
+        return {}
+
+    def predict(self, X, **kwargs):
+        self.predict_calls.append(copy.deepcopy(kwargs))
+        return np.array(X, copy=True)
+
+
+def _low_rank_dataset(
+    n_samples: int = 60,
+    sample_shape: tuple[int, ...] = (3, 4),
+    rank: int = 3,
+    seed: int = 0,
+) -> np.ndarray:
+    """Build deterministic synthetic data lying exactly on a low-rank subspace."""
+    rng = np.random.default_rng(seed)
+    n_features = int(np.prod(sample_shape))
+    latent = rng.normal(size=(n_samples, rank))
+    mixing = rng.normal(size=(rank, n_features))
+    return (latent @ mixing).reshape(n_samples, *sample_shape)
+
+
+def _split_for(n_samples: int, fractions=(0.6, 0.2, 0.2)):
+    return split_chronologically(
+        sample_times=np.arange(n_samples),
+        fractions=fractions,
+    )
+
+
+def _digest(array: np.ndarray) -> str:
+    contiguous = np.ascontiguousarray(array, dtype=np.float64)
+    return hashlib.sha256(contiguous.tobytes()).hexdigest()
+
+
+class _RecordingMethod:
+    """Controlled benchmark method that records every array it is handed."""
+
+    def __init__(self, latent_dimension=2, transform=None, draw_random=False):
+        self.latent_dimension = latent_dimension
+        self.is_fitted = False
+        self._transform = transform
+        self._draw_random = draw_random
+        self.fit_train = None
+        self.fit_train_object = None
+        self.fit_validation = None
+        self.reconstruct_inputs = []
+        self.reconstruct_input_objects = []
+        self.random_draws = []
+
+    def fit(self, X_train, X_validation):
+        self.fit_train_object = X_train
+        self.fit_train = np.array(X_train, copy=True)
+        self.fit_validation = (
+            None if X_validation is None else np.array(X_validation, copy=True)
+        )
+        if self._draw_random:
+            self.random_draws.append(float(np.random.rand()))
+            self.random_draws.append(float(torch.rand(1).item()))
+        self.is_fitted = True
+
+    def reconstruct(self, X):
+        self.reconstruct_input_objects.append(X)
+        self.reconstruct_inputs.append(np.array(X, copy=True))
+        if self._transform is None:
+            return np.array(X, copy=True)
+        return self._transform(X)
+
+
+def _recording_spec(
+    name: str = "recorder",
+    *,
+    latent_dimension: int = 2,
+    transform=None,
+    uses_validation_partition: bool = True,
+    draw_random: bool = False,
+    declared_latent_dimension: int | None = None,
+):
+    """Return a specification plus the list receiving every built instance."""
+    created: list[_RecordingMethod] = []
+
+    def factory():
+        instance = _RecordingMethod(
+            latent_dimension=latent_dimension,
+            transform=transform,
+            draw_random=draw_random,
+        )
+        created.append(instance)
+        return instance
+
+    specification = BenchmarkMethod(
+        name=name,
+        method_type="controlled-fake",
+        latent_dimension=(
+            latent_dimension
+            if declared_latent_dimension is None
+            else declared_latent_dimension
+        ),
+        factory=factory,
+        configuration={"kind": "controlled-fake"},
+        uses_validation_partition=uses_validation_partition,
+    )
+    return specification, created
+
+
+def _pca_spec(name: str = "pca", *, n_components: int = 3, scale_data: bool = False):
+    """Return a PCA specification plus the list receiving every built adapter."""
+    created: list[PCAReconstruction] = []
+
+    def factory():
+        instance = PCAReconstruction(
+            n_components=n_components,
+            scale_data=scale_data,
+        )
+        created.append(instance)
+        return instance
+
+    specification = BenchmarkMethod(
+        name=name,
+        method_type="pca",
+        latent_dimension=n_components,
+        factory=factory,
+        configuration={"n_components": n_components, "scale_data": scale_data},
+        uses_validation_partition=False,
+    )
+    return specification, created
+
+
+# ---------------------------------------------------------------------------
+# A. PCA correctness baseline
+# ---------------------------------------------------------------------------
+
+
+def test_pca_reconstructs_low_rank_data_with_sufficient_components():
+    X = _low_rank_dataset(rank=3)
+    split = _split_for(len(X))
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[pca_benchmark_method("pca-k3", n_components=3)],
+    )
+
+    result = report.results[0]
+    assert result.test_metrics["mse"] < 1e-20
+    assert result.test_metrics["mae"] < 1e-10
+    assert result.test_metrics["rmse"] < 1e-10
+
+
+def test_pca_adapter_matches_the_existing_pca_api_semantics():
+    X = _low_rank_dataset(rank=3)
+    split = _split_for(len(X))
+    X_train = X[split.train_indices]
+    X_test = X[split.test_indices]
+
+    adapter = PCAReconstruction(n_components=3)
+    adapter.fit(X_train)
+    reconstruction = adapter.reconstruct(X_test)
+
+    # The adapter must not invent a second PCA: it delegates to the fitted
+    # scikit-learn estimator owned by bluemath_tk.datamining.pca.PCA.
+    estimator = adapter.pca.pca
+    expected = estimator.inverse_transform(
+        estimator.transform(X_test.reshape(len(X_test), -1))
+    ).reshape(X_test.shape)
+    assert np.array_equal(reconstruction, expected)
+    assert adapter.pca.is_fitted is True
+    assert adapter.latent_dimension == 3
+
+
+def test_fewer_components_than_rank_degrades_but_stays_finite():
+    X = _low_rank_dataset(rank=3)
+    split = _split_for(len(X))
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[
+            pca_benchmark_method("pca-k1", n_components=1),
+            pca_benchmark_method("pca-k3", n_components=3),
+        ],
+    )
+
+    poor, good = report.results
+    assert poor.test_metrics["mse"] > good.test_metrics["mse"]
+    assert np.isfinite(poor.test_metrics["mse"])
+
+
+# ---------------------------------------------------------------------------
+# B. Shape round trip and sample ordering
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("sample_shape", [(7,), (3, 4), (2, 3, 4), (2, 1, 3, 2)])
+def test_pca_preserves_sample_shape_for_n_dimensional_samples(sample_shape):
+    X = _low_rank_dataset(n_samples=40, sample_shape=sample_shape, rank=2)
+    split = _split_for(len(X))
+    X_train = X[split.train_indices]
+    X_test = X[split.test_indices]
+
+    adapter = PCAReconstruction(n_components=2)
+    adapter.fit(X_train)
+    reconstruction = adapter.reconstruct(X_test)
+
+    assert reconstruction.shape == X_test.shape
+    assert np.allclose(reconstruction, X_test, atol=1e-8)
+
+
+def test_pca_preserves_sample_order():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(3, 4), rank=3)
+    split = _split_for(len(X))
+    X_train = X[split.train_indices]
+    X_test = X[split.test_indices]
+
+    adapter = PCAReconstruction(n_components=3)
+    adapter.fit(X_train)
+    reconstruction = adapter.reconstruct(X_test)
+
+    aligned = np.max(np.abs(reconstruction - X_test), axis=tuple(range(1, X.ndim)))
+    rolled = np.max(
+        np.abs(reconstruction - np.roll(X_test, 1, axis=0)),
+        axis=tuple(range(1, X.ndim)),
+    )
+    assert np.all(aligned < 1e-8)
+    assert np.all(rolled > 1e-6)
+
+
+def test_results_are_independent_of_the_input_memory_layout():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(3, 4), rank=3)
+    fortran = np.asfortranarray(X)
+    assert fortran.flags.f_contiguous
+    assert np.array_equal(fortran, X)
+    split = _split_for(len(X))
+    method = pca_benchmark_method("pca-k3", n_components=3)
+
+    c_order = run_reconstruction_benchmark(X, split=split, methods=[method])
+    f_order = run_reconstruction_benchmark(fortran, split=split, methods=[method])
+
+    assert dict(c_order.results[0].test_metrics) == dict(
+        f_order.results[0].test_metrics
+    )
+    assert c_order.identity_digest() == f_order.identity_digest()
+
+
+def test_pca_reconstruction_is_identical_for_fortran_ordered_samples():
+    X = _low_rank_dataset(n_samples=30, sample_shape=(3, 5), rank=4)
+    train, test = X[:20], X[20:]
+
+    c_adapter = PCAReconstruction(n_components=4)
+    c_adapter.fit(train)
+    f_adapter = PCAReconstruction(n_components=4)
+    f_adapter.fit(np.asfortranarray(train))
+
+    assert np.array_equal(
+        c_adapter.reconstruct(test),
+        f_adapter.reconstruct(np.asfortranarray(test)),
+    )
+
+
+def test_pca_reconstruction_uses_c_order_flattening():
+    X = _low_rank_dataset(n_samples=30, sample_shape=(3, 5), rank=4)
+    adapter = PCAReconstruction(n_components=4)
+    adapter.fit(X[:20])
+
+    stacked = adapter.pca.stacked_data_matrix
+    assert np.array_equal(stacked, X[:20].reshape(20, -1))
+
+
+# ---------------------------------------------------------------------------
+# C. Shared metric agreement
+# ---------------------------------------------------------------------------
+
+
+def test_benchmark_metrics_agree_with_shared_metric_implementation():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec(transform=lambda values: values + 0.25)
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[specification],
+        metrics=METRICS,
+    )
+
+    X_test = X[split.test_indices]
+    result = report.results[0]
+    for metric in METRICS:
+        expected = float(
+            reconstruction_error(
+                X_test,
+                X_test + 0.25,
+                metric=metric,
+                reduction="mean",
+            )
+        )
+        assert result.test_metrics[metric] == expected
+
+
+def test_reported_metrics_preserve_the_requested_order():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec()
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[specification],
+        metrics=("rmse", "mse"),
+    )
+    assert report.metrics == ("rmse", "mse")
+    assert set(report.results[0].test_metrics) == {"rmse", "mse"}
+
+
+def test_metrics_are_computed_only_on_the_test_partition():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, created = _recording_spec()
+
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    reconstruct_inputs = created[0].reconstruct_inputs
+    assert len(reconstruct_inputs) == 1
+    assert np.array_equal(reconstruct_inputs[0], X[split.test_indices])
+
+
+# ---------------------------------------------------------------------------
+# D. Chronological split membership
+# ---------------------------------------------------------------------------
+
+
+def test_runner_uses_exactly_the_supplied_chronological_partitions():
+    X = _low_rank_dataset(n_samples=53)
+    split = _split_for(len(X), fractions=(0.5, 0.25, 0.25))
+    specification, created = _recording_spec()
+
+    report = run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    instance = created[0]
+    assert np.array_equal(instance.fit_train, X[split.train_indices])
+    assert np.array_equal(instance.fit_validation, X[split.validation_indices])
+    assert np.array_equal(instance.reconstruct_inputs[0], X[split.test_indices])
+    assert report.partition_sizes["train"] == int(split.train_indices.size)
+    assert report.partition_sizes["validation"] == int(split.validation_indices.size)
+    assert report.partition_sizes["test"] == int(split.test_indices.size)
+
+
+def test_partitions_have_no_off_by_one_boundaries():
+    X = _low_rank_dataset(n_samples=41)
+    split = _split_for(len(X), fractions=(0.6, 0.2, 0.2))
+    specification, created = _recording_spec()
+
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    instance = created[0]
+    assert len(instance.fit_train) == int(split.train_indices.size)
+    assert len(instance.fit_validation) == int(split.validation_indices.size)
+    # The first validation sample must not appear in training and the last
+    # training sample must not appear in validation.
+    assert not np.array_equal(instance.fit_train[-1], X[split.validation_indices[0]])
+    assert np.array_equal(instance.fit_train[-1], X[split.train_indices[-1]])
+    assert np.array_equal(instance.fit_validation[0], X[split.validation_indices[0]])
+
+
+def test_runner_rejects_a_split_describing_a_different_sample_count():
+    X = _low_rank_dataset(n_samples=40)
+    split = _split_for(30)
+    specification, _ = _recording_spec()
+
+    with pytest.raises(ValueError, match="describes 30 samples"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+def test_runner_requires_a_manifest_backed_chronological_split():
+    X = _low_rank_dataset()
+    specification, _ = _recording_spec()
+
+    with pytest.raises(TypeError, match="ChronologicalSplit"):
+        run_reconstruction_benchmark(
+            X,
+            split={"train": [0], "validation": [1], "test": [2]},
+            methods=[specification],
+        )
+
+
+def test_optional_time_axis_verification_detects_a_replayed_wrong_dataset():
+    X = _low_rank_dataset(n_samples=40)
+    times = np.arange(40)
+    split = split_chronologically(sample_times=times, fractions=(0.6, 0.2, 0.2))
+    specification, _ = _recording_spec()
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[specification],
+        sample_times=times,
+    )
+    assert report.split_identity["time_axis_verified"] is True
+
+    specification, _ = _recording_spec()
+    with pytest.raises(ValueError, match="fingerprint"):
+        run_reconstruction_benchmark(
+            X,
+            split=split,
+            methods=[specification],
+            sample_times=times + 1,
+        )
+
+
+def test_time_axis_verification_is_recorded_as_false_when_not_requested():
+    X = _low_rank_dataset(n_samples=40)
+    split = _split_for(40)
+    specification, _ = _recording_spec()
+
+    report = run_reconstruction_benchmark(X, split=split, methods=[specification])
+    assert report.split_identity["time_axis_verified"] is False
+
+
+# ---------------------------------------------------------------------------
+# E. Leakage adversary
+# ---------------------------------------------------------------------------
+
+
+def test_altering_test_values_does_not_change_fitted_pca_state():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+
+    first_spec, first_created = _pca_spec()
+    run_reconstruction_benchmark(X, split=split, methods=[first_spec])
+
+    attacked = X.copy()
+    attacked[split.test_indices] = attacked[split.test_indices] * 1e6 + 12345.0
+    second_spec, second_created = _pca_spec()
+    run_reconstruction_benchmark(attacked, split=split, methods=[second_spec])
+
+    original = first_created[0].pca.pca
+    adversarial = second_created[0].pca.pca
+    assert np.array_equal(original.components_, adversarial.components_)
+    assert np.array_equal(original.mean_, adversarial.mean_)
+    assert np.array_equal(original.explained_variance_, adversarial.explained_variance_)
+    assert np.array_equal(
+        first_created[0].pca.stacked_data_matrix,
+        second_created[0].pca.stacked_data_matrix,
+    )
+
+
+def test_altering_test_values_does_not_change_the_data_reaching_fitting():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+
+    first_spec, first_created = _recording_spec()
+    run_reconstruction_benchmark(X, split=split, methods=[first_spec])
+
+    attacked = X.copy()
+    attacked[split.test_indices] = 9.9e5
+    second_spec, second_created = _recording_spec()
+    run_reconstruction_benchmark(attacked, split=split, methods=[second_spec])
+
+    assert _digest(first_created[0].fit_train) == _digest(second_created[0].fit_train)
+    assert _digest(first_created[0].fit_validation) == _digest(
+        second_created[0].fit_validation
+    )
+
+
+def test_pca_scaler_statistics_come_from_the_training_partition_only():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, created = _pca_spec(n_components=3, scale_data=True)
+
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    scaler = created[0].pca.scaler
+    train_flat = X[split.train_indices].reshape(int(split.train_indices.size), -1)
+    assert np.allclose(scaler.mean_, train_flat.mean(axis=0))
+    assert not np.allclose(scaler.mean_, X.reshape(len(X), -1).mean(axis=0))
+
+
+def test_pca_is_not_given_the_validation_partition():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, created = _recording_spec(uses_validation_partition=False)
+
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    assert created[0].fit_validation is None
+    assert pca_benchmark_method("p", n_components=2).uses_validation_partition is False
+
+
+def test_pca_adapter_ignores_any_validation_data_it_is_handed():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    X_train = X[split.train_indices]
+
+    without = PCAReconstruction(n_components=3)
+    without.fit(X_train, None)
+    with_validation = PCAReconstruction(n_components=3)
+    with_validation.fit(X_train, X[split.validation_indices])
+
+    assert np.array_equal(
+        without.pca.pca.components_,
+        with_validation.pca.pca.components_,
+    )
+    assert np.array_equal(without.pca.pca.mean_, with_validation.pca.pca.mean_)
+    assert np.array_equal(
+        without.pca.stacked_data_matrix,
+        with_validation.pca.stacked_data_matrix,
+    )
+    assert without.pca.stacked_data_matrix.shape[0] == int(split.train_indices.size)
+
+
+# ---------------------------------------------------------------------------
+# F. Validation-set fidelity
+# ---------------------------------------------------------------------------
+
+
+def test_autoencoder_adapter_forwards_the_exact_validation_partition():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=3)
+    split = _split_for(len(X))
+    captured = {}
+
+    class _CapturingModel:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X_train, validation_data=None, **kwargs):
+            captured["train"] = np.array(X_train, copy=True)
+            captured["validation"] = np.array(validation_data[0], copy=True)
+            captured["validation_target"] = validation_data[1]
+            captured["kwargs"] = dict(kwargs)
+            self.is_fitted = True
+            return {"train_loss": [0.0], "val_loss": [0.0]}
+
+        def predict(self, X, **kwargs):
+            return np.array(X, copy=True)
+
+    specification = autoencoder_benchmark_method(
+        "capturing",
+        model_factory=_CapturingModel,
+        latent_dimension=2,
+        configuration={"architecture": "capturing"},
+        shuffle_training_data=False,
+    )
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    assert np.array_equal(captured["train"], X[split.train_indices])
+    assert np.array_equal(captured["validation"], X[split.validation_indices])
+    assert captured["validation_target"] is None
+    assert captured["kwargs"]["verbose"] == 0
+
+
+def test_test_samples_never_reach_autoencoder_fitting():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=3)
+    split = _split_for(len(X))
+    seen = []
+
+    class _WatchingModel:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X_train, validation_data=None, **kwargs):
+            seen.append(np.array(X_train, copy=True))
+            seen.append(np.array(validation_data[0], copy=True))
+            self.is_fitted = True
+            return {}
+
+        def predict(self, X, **kwargs):
+            return np.array(X, copy=True)
+
+    specification = autoencoder_benchmark_method(
+        "watching",
+        model_factory=_WatchingModel,
+        latent_dimension=2,
+    )
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    test_rows = {row.tobytes() for row in X[split.test_indices]}
+    for array in seen:
+        for row in array:
+            assert row.tobytes() not in test_rows
+
+
+def test_autoencoder_adapter_requires_a_validation_partition():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=3)
+    split = _split_for(len(X))
+
+    class _Model:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X_train, validation_data=None, **kwargs):
+            self.is_fitted = True
+            return {}
+
+        def predict(self, X, **kwargs):
+            return np.array(X, copy=True)
+
+    base = autoencoder_benchmark_method(
+        "no-validation",
+        model_factory=_Model,
+        latent_dimension=2,
+    )
+    specification = BenchmarkMethod(
+        name=base.name,
+        method_type=base.method_type,
+        latent_dimension=base.latent_dimension,
+        factory=base.factory,
+        uses_validation_partition=False,
+    )
+
+    with pytest.raises(ValueError, match="requires the validation partition"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+def test_models_that_would_swallow_validation_data_are_rejected():
+    class _Swallowing:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X, **kwargs):
+            return {}
+
+        def predict(self, X, **kwargs):
+            return X
+
+    with pytest.raises(TypeError, match="explicit\\s+validation_data parameter"):
+        AutoencoderReconstruction(_Swallowing(), latent_dimension=2)
+
+
+def test_every_shipped_autoencoder_declares_validation_data():
+    for model_class in (
+        StandardAutoencoder,
+        OrthogonalAutoencoder,
+        VariationalAutoencoder,
+    ):
+        parameter = inspect.signature(model_class.fit).parameters.get("validation_data")
+        assert parameter is not None
+        assert parameter.kind is not inspect.Parameter.VAR_KEYWORD
+
+
+def test_training_data_is_shuffled_by_default_without_changing_membership():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=3)
+    split = _split_for(len(X))
+    captured = {}
+
+    class _CapturingModel:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X_train, validation_data=None, **kwargs):
+            captured["train"] = np.array(X_train, copy=True)
+            captured["validation"] = np.array(validation_data[0], copy=True)
+            self.is_fitted = True
+            return {}
+
+        def predict(self, X, **kwargs):
+            return np.array(X, copy=True)
+
+    shuffled = autoencoder_benchmark_method(
+        "shuffled",
+        model_factory=_CapturingModel,
+        latent_dimension=2,
+    )
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[shuffled],
+        seed=0,
+    )
+    expected_train = X[split.train_indices]
+    assert not np.array_equal(captured["train"], expected_train)
+    # Shuffling reorders the training rows but never changes membership.
+    assert sorted(row.tobytes() for row in captured["train"]) == sorted(
+        row.tobytes() for row in expected_train
+    )
+    assert np.array_equal(captured["validation"], X[split.validation_indices])
+    assert report.results[0].configuration["shuffle_training_data"] is True
+
+    ordered = autoencoder_benchmark_method(
+        "ordered",
+        model_factory=_CapturingModel,
+        latent_dimension=2,
+        shuffle_training_data=False,
+    )
+    report = run_reconstruction_benchmark(X, split=split, methods=[ordered], seed=0)
+    assert np.array_equal(captured["train"], expected_train)
+    assert report.results[0].configuration["shuffle_training_data"] is False
+
+
+def test_training_shuffle_is_reproducible_and_isolated():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=3)
+    split = _split_for(len(X))
+    seen = []
+
+    class _CapturingModel:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X_train, validation_data=None, **kwargs):
+            seen.append(np.array(X_train, copy=True))
+            self.is_fitted = True
+            return {}
+
+        def predict(self, X, **kwargs):
+            return np.array(X, copy=True)
+
+    def _method(name):
+        return autoencoder_benchmark_method(
+            name,
+            model_factory=_CapturingModel,
+            latent_dimension=2,
+        )
+
+    np.random.seed(2024)
+    state_before = np.random.get_state()
+    run_reconstruction_benchmark(X, split=split, methods=[_method("a")], seed=4)
+    run_reconstruction_benchmark(X, split=split, methods=[_method("b")], seed=4)
+    state_after = np.random.get_state()
+
+    assert np.array_equal(seen[0], seen[1])
+    assert np.array_equal(state_before[1], state_after[1])
+
+
+def test_reserved_shuffle_configuration_key_is_rejected():
+    class _Model:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X, validation_data=None, **kwargs):
+            return {}
+
+        def predict(self, X, **kwargs):
+            return X
+
+    with pytest.raises(ValueError, match="reserved key"):
+        autoencoder_benchmark_method(
+            "clash",
+            model_factory=_Model,
+            latent_dimension=2,
+            configuration={"shuffle_training_data": False},
+        )
+
+
+# ---------------------------------------------------------------------------
+# G. No mutation of caller state
+# ---------------------------------------------------------------------------
+
+
+def test_runner_does_not_mutate_data_split_or_configuration():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    X_before = X.copy()
+    train_before = split.train_indices.copy()
+    validation_before = split.validation_indices.copy()
+    test_before = split.test_indices.copy()
+    configuration = {"nested": {"values": [1, 2, 3]}}
+    configuration_before = copy.deepcopy(configuration)
+
+    specification = BenchmarkMethod(
+        name="pca",
+        method_type="pca",
+        latent_dimension=3,
+        factory=lambda: PCAReconstruction(n_components=3),
+        configuration=configuration,
+        uses_validation_partition=False,
+    )
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    assert np.array_equal(X, X_before)
+    assert np.array_equal(split.train_indices, train_before)
+    assert np.array_equal(split.validation_indices, validation_before)
+    assert np.array_equal(split.test_indices, test_before)
+    assert configuration == configuration_before
+    assert split.train_indices.flags.writeable is False
+    assert split.test_indices.flags.writeable is False
+
+
+def test_methods_receive_copies_that_do_not_share_memory_with_the_input():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, created = _recording_spec()
+
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    instance = created[0]
+    assert not np.shares_memory(instance.fit_train_object, X)
+    assert not np.shares_memory(instance.reconstruct_input_objects[0], X)
+
+
+def test_a_method_that_mutates_its_partition_cannot_corrupt_the_input():
+    X = _low_rank_dataset()
+    X_before = X.copy()
+    split = _split_for(len(X))
+
+    def _destructive(values):
+        values[...] = 0.0
+        return np.array(values, copy=True)
+
+    specification, _ = _recording_spec(transform=_destructive)
+    report = run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    assert np.array_equal(X, X_before)
+    # Metrics must be scored against the pristine test partition, not against
+    # the copy the method just zeroed.
+    expected = float(
+        reconstruction_error(
+            X[split.test_indices],
+            np.zeros_like(X[split.test_indices]),
+            metric="mse",
+            reduction="mean",
+        )
+    )
+    assert report.results[0].test_metrics["mse"] == expected
+    assert report.results[0].test_metrics["mse"] > 0.0
+
+
+def test_a_mutating_method_cannot_corrupt_later_methods():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+
+    class _Vandal:
+        latent_dimension = 2
+        is_fitted = False
+
+        def fit(self, X_train, X_validation):
+            X_train[...] = 0.0
+            if X_validation is not None:
+                X_validation[...] = 0.0
+            self.is_fitted = True
+
+        def reconstruct(self, X_predict):
+            X_predict[...] = 0.0
+            return np.zeros_like(X_predict)
+
+    vandal = BenchmarkMethod(
+        name="vandal",
+        method_type="controlled-fake",
+        latent_dimension=2,
+        factory=_Vandal,
+        uses_validation_partition=True,
+    )
+    victim, victim_created = _recording_spec(name="victim")
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[vandal, victim],
+    )
+
+    instance = victim_created[0]
+    assert np.array_equal(instance.fit_train, X[split.train_indices])
+    assert np.array_equal(instance.fit_validation, X[split.validation_indices])
+    assert np.array_equal(instance.reconstruct_inputs[0], X[split.test_indices])
+    # The victim reconstructs its input perfectly, so a corrupted test
+    # partition would have shown up as a spuriously perfect score too.
+    assert report.results[1].test_metrics["mse"] == 0.0
+    assert report.results[0].test_metrics["mse"] > 0.0
+
+
+def test_each_run_builds_a_fresh_method_instance():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, created = _recording_spec()
+
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    assert len(created) == 2
+    assert created[0] is not created[1]
+
+
+def test_runner_rejects_an_already_fitted_instance():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    shared = PCAReconstruction(n_components=3)
+    shared.fit(X[split.train_indices])
+
+    specification = BenchmarkMethod(
+        name="reused",
+        method_type="pca",
+        latent_dimension=3,
+        factory=lambda: shared,
+        uses_validation_partition=False,
+    )
+    with pytest.raises(ValueError, match="already fitted instance"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+def test_runner_rejects_one_instance_shared_by_two_methods():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+
+    class _StatelessMethod:
+        """A method that never reports being fitted, defeating the fit guard."""
+
+        latent_dimension = 2
+        is_fitted = False
+
+        def fit(self, X_train, X_validation):
+            return None
+
+        def reconstruct(self, X):
+            return np.array(X, copy=True)
+
+    shared = _StatelessMethod()
+    methods = [
+        BenchmarkMethod(
+            name=name,
+            method_type="controlled-fake",
+            latent_dimension=2,
+            factory=lambda: shared,
+            uses_validation_partition=False,
+        )
+        for name in ("first", "second")
+    ]
+    with pytest.raises(ValueError, match="already used by another method"):
+        run_reconstruction_benchmark(X, split=split, methods=methods)
+
+
+# ---------------------------------------------------------------------------
+# H. Invalid latent dimensions
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("n_components", [0, -1, -10])
+def test_non_positive_latent_dimension_is_rejected(n_components):
+    with pytest.raises(ValueError, match="positive integer"):
+        pca_benchmark_method("bad", n_components=n_components)
+
+
+@pytest.mark.parametrize("n_components", [1.5, True, "3", None])
+def test_non_integer_latent_dimension_is_rejected(n_components):
+    with pytest.raises(TypeError, match="exact non-Boolean integer"):
+        pca_benchmark_method("bad", n_components=n_components)
+
+
+def test_pca_rejects_more_components_than_available():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(5,), rank=3)
+    split = _split_for(len(X))
+
+    with pytest.raises(ValueError, match="components available"):
+        run_reconstruction_benchmark(
+            X,
+            split=split,
+            methods=[pca_benchmark_method("too-many", n_components=6)],
+        )
+
+
+@pytest.mark.parametrize("n_components", [6, 7, 10])
+def test_pca_rejects_more_components_than_the_centered_training_rank(n_components):
+    X = _low_rank_dataset(n_samples=12, sample_shape=(20,), rank=3)
+    split = _split_for(len(X), fractions=(0.5, 0.25, 0.25))
+    assert int(split.train_indices.size) == 6
+
+    # Centering costs one degree of freedom, so 6 training samples support at
+    # most 5 components.
+    with pytest.raises(ValueError, match="components available"):
+        run_reconstruction_benchmark(
+            X,
+            split=split,
+            methods=[pca_benchmark_method("too-many", n_components=n_components)],
+        )
+
+
+def test_pca_accepts_the_largest_meaningful_component_count():
+    X = _low_rank_dataset(n_samples=12, sample_shape=(20,), rank=3)
+    split = _split_for(len(X), fractions=(0.5, 0.25, 0.25))
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[pca_benchmark_method("max-k", n_components=5)],
+    )
+    assert report.results[0].latent_dimension == 5
+
+
+def test_autoencoder_latent_width_must_match_the_specification():
+    class _Model:
+        k = 3
+        is_fitted = False
+
+        def fit(self, X, validation_data=None, **kwargs):
+            return {}
+
+        def predict(self, X, **kwargs):
+            return X
+
+    with pytest.raises(ValueError, match="contradicts the model latent width"):
+        AutoencoderReconstruction(_Model(), latent_dimension=5)
+
+
+def test_method_reporting_a_different_latent_dimension_is_rejected():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec(
+        latent_dimension=2,
+        declared_latent_dimension=4,
+    )
+
+    with pytest.raises(ValueError, match="reports latent dimension"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+# ---------------------------------------------------------------------------
+# I. Reconstruction shape mismatch
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize(
+    "transform",
+    [
+        lambda values: values[:, :1, :],
+        lambda values: values.mean(axis=1, keepdims=True),
+        lambda values: values[:1],
+        lambda values: values.reshape(len(values), -1),
+        lambda values: np.zeros((len(values), 1, 1)),
+    ],
+)
+def test_broadcastable_but_wrong_reconstruction_shape_is_rejected(transform):
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec(transform=transform)
+
+    with pytest.raises(ValueError, match="Broadcasting is never applied"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+def test_non_array_reconstruction_is_rejected():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec(transform=lambda values: values.tolist())
+
+    with pytest.raises(TypeError, match="NumPy array"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+def test_complex_reconstruction_is_rejected():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec(
+        transform=lambda values: values.astype(np.complex128),
+    )
+
+    with pytest.raises(TypeError, match="real-valued"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+# ---------------------------------------------------------------------------
+# J. Non-finite and invalid input data
+# ---------------------------------------------------------------------------
+
+
+@pytest.mark.parametrize("bad_value", [np.nan, np.inf, -np.inf])
+def test_non_finite_datasets_are_rejected(bad_value):
+    X = _low_rank_dataset()
+    X = X.copy()
+    X[3, 1, 2] = bad_value
+    split = _split_for(len(X))
+    specification, _ = _recording_spec()
+
+    with pytest.raises(ValueError, match="NaN or infinite"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+def test_non_finite_reconstruction_is_rejected():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec(
+        transform=lambda values: np.full_like(values, np.nan),
+    )
+
+    with pytest.raises(ValueError, match="non-finite reconstruction"):
+        run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+
+def test_complex_and_non_numeric_datasets_are_rejected():
+    split = _split_for(30)
+    specification, _ = _recording_spec()
+
+    complex_data = np.ones((30, 4), dtype=np.complex128)
+    with pytest.raises(TypeError, match="real-valued"):
+        run_reconstruction_benchmark(
+            complex_data,
+            split=split,
+            methods=[specification],
+        )
+
+    text_data = np.full((30, 4), "a", dtype="<U1")
+    with pytest.raises(TypeError, match="numeric values"):
+        run_reconstruction_benchmark(text_data, split=split, methods=[specification])
+
+
+def test_empty_and_one_dimensional_datasets_are_rejected():
+    split = _split_for(30)
+    specification, _ = _recording_spec()
+
+    with pytest.raises(ValueError, match="leading sample dimension"):
+        run_reconstruction_benchmark(
+            np.arange(30.0),
+            split=split,
+            methods=[specification],
+        )
+
+    with pytest.raises(ValueError, match="at least one sample"):
+        run_reconstruction_benchmark(
+            np.empty((0, 4)),
+            split=split,
+            methods=[specification],
+        )
+
+    with pytest.raises(ValueError, match="per-sample dimension"):
+        run_reconstruction_benchmark(
+            np.empty((30, 0)),
+            split=split,
+            methods=[specification],
+        )
+
+
+# ---------------------------------------------------------------------------
+# K. Duplicate method names and invalid specifications
+# ---------------------------------------------------------------------------
+
+
+def test_duplicate_method_names_are_rejected():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    first, _ = _recording_spec(name="same")
+    second, _ = _recording_spec(name="same")
+
+    with pytest.raises(ValueError, match="unique"):
+        run_reconstruction_benchmark(X, split=split, methods=[first, second])
+
+
+def test_empty_or_invalid_method_sequences_are_rejected():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+
+    with pytest.raises(ValueError, match="at least one BenchmarkMethod"):
+        run_reconstruction_benchmark(X, split=split, methods=[])
+
+    with pytest.raises(TypeError, match="must be a BenchmarkMethod"):
+        run_reconstruction_benchmark(X, split=split, methods=[object()])
+
+
+def test_invalid_metric_requests_are_rejected():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec()
+
+    with pytest.raises(ValueError, match="Unsupported metric"):
+        run_reconstruction_benchmark(
+            X,
+            split=split,
+            methods=[specification],
+            metrics=("r2",),
+        )
+
+    with pytest.raises(ValueError, match="Duplicate metric"):
+        run_reconstruction_benchmark(
+            X,
+            split=split,
+            methods=[specification],
+            metrics=("mse", "mse"),
+        )
+
+    with pytest.raises(TypeError, match="not a string"):
+        run_reconstruction_benchmark(
+            X,
+            split=split,
+            methods=[specification],
+            metrics="mse",
+        )
+
+
+def test_non_serializable_configuration_is_rejected():
+    with pytest.raises(TypeError, match="non-JSON value"):
+        BenchmarkMethod(
+            name="bad",
+            method_type="fake",
+            latent_dimension=1,
+            factory=lambda: None,
+            configuration={"callable": len},
+        )
+
+
+def test_factory_must_be_callable():
+    with pytest.raises(TypeError, match="zero-argument callable"):
+        BenchmarkMethod(
+            name="bad",
+            method_type="fake",
+            latent_dimension=1,
+            factory="not-callable",
+        )
+
+
+# ---------------------------------------------------------------------------
+# L. Timing
+# ---------------------------------------------------------------------------
+
+
+def test_timings_are_finite_and_non_negative():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec()
+
+    report = run_reconstruction_benchmark(X, split=split, methods=[specification])
+    result = report.results[0]
+
+    assert np.isfinite(result.fit_seconds)
+    assert np.isfinite(result.reconstruction_seconds)
+    assert result.fit_seconds >= 0.0
+    assert result.reconstruction_seconds >= 0.0
+
+    timing = report.to_dict()["results"][0]["timing"]
+    assert set(timing) == {"fit_seconds", "reconstruction_seconds"}
+    assert all(type(value) is float for value in timing.values())
+    assert all(np.isfinite(value) and value >= 0.0 for value in timing.values())
+
+
+# ---------------------------------------------------------------------------
+# M. Reproducible metadata
+# ---------------------------------------------------------------------------
+
+
+def _example_report() -> ReconstructionBenchmarkReport:
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    return run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[
+            pca_benchmark_method("pca-k2", n_components=2),
+            pca_benchmark_method("pca-k3", n_components=3),
+        ],
+        seed=5,
+    )
+
+
+def test_report_dictionary_and_json_round_trip():
+    report = _example_report()
+
+    assert ReconstructionBenchmarkReport.from_dict(report.to_dict()).to_dict() == (
+        report.to_dict()
+    )
+    parsed = json.loads(report.to_json())
+    assert ReconstructionBenchmarkReport.from_dict(parsed).to_dict() == report.to_dict()
+
+
+def test_report_json_is_deterministic_and_sorted():
+    report = _example_report()
+
+    first = report.to_json()
+    assert first == report.to_json()
+    assert first.endswith("\n")
+
+    payload = json.loads(first)
+    assert list(payload) == sorted(payload)
+    assert list(payload["results"][0]) == sorted(payload["results"][0])
+
+
+def test_identity_excludes_timing_and_measured_metrics():
+    report = _example_report()
+    identity = report.identity()
+
+    serialized = json.dumps(identity, sort_keys=True)
+    assert "fit_seconds" not in serialized
+    assert "reconstruction_seconds" not in serialized
+    assert "test_metrics" not in serialized
+    assert "timing" not in serialized
+
+
+def test_identity_digest_ignores_timing_and_metric_values():
+    report = _example_report()
+    payload = report.to_dict()
+
+    tampered = copy.deepcopy(payload)
+    tampered["results"][0]["timing"]["fit_seconds"] = 999.0
+    tampered["results"][0]["timing"]["reconstruction_seconds"] = 888.0
+    tampered["results"][0]["test_metrics"]["mse"] = 42.0
+
+    rebuilt = ReconstructionBenchmarkReport.from_dict(tampered)
+    assert rebuilt.identity_digest() == report.identity_digest()
+    assert rebuilt.to_dict() != payload
+
+
+def test_identity_digest_distinguishes_different_datasets():
+    split = _split_for(60)
+    method = pca_benchmark_method("pca-k2", n_components=2)
+
+    first = run_reconstruction_benchmark(
+        _low_rank_dataset(seed=0),
+        split=split,
+        methods=[method],
+    )
+    second = run_reconstruction_benchmark(
+        _low_rank_dataset(seed=1),
+        split=split,
+        methods=[method],
+    )
+
+    assert first.data_digest != second.data_digest
+    assert first.identity_digest() != second.identity_digest()
+
+
+def test_data_digest_is_stable_and_layout_independent():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    method = pca_benchmark_method("pca-k2", n_components=2)
+
+    first = run_reconstruction_benchmark(X, split=split, methods=[method])
+    second = run_reconstruction_benchmark(X.copy(), split=split, methods=[method])
+    fortran = run_reconstruction_benchmark(
+        np.asfortranarray(X),
+        split=split,
+        methods=[method],
+    )
+
+    assert first.data_digest == second.data_digest == fortran.data_digest
+    assert len(first.data_digest) == 64
+
+
+def test_report_rejects_a_malformed_data_digest():
+    report = _example_report()
+    payload = report.to_dict()
+    payload["data_digest"] = "not-a-digest"
+
+    with pytest.raises(ValueError, match="SHA-256 hex digest"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+def test_identity_digest_changes_when_the_split_changes():
+    X = _low_rank_dataset()
+    method = pca_benchmark_method("pca-k2", n_components=2)
+
+    first = run_reconstruction_benchmark(
+        X,
+        split=_split_for(len(X), fractions=(0.6, 0.2, 0.2)),
+        methods=[method],
+    )
+    second = run_reconstruction_benchmark(
+        X,
+        split=_split_for(len(X), fractions=(0.5, 0.25, 0.25)),
+        methods=[method],
+    )
+    assert first.identity_digest() != second.identity_digest()
+
+
+def test_split_identity_is_stable_for_the_same_split():
+    X = _low_rank_dataset()
+    method = pca_benchmark_method("pca-k2", n_components=2)
+
+    first = run_reconstruction_benchmark(X, split=_split_for(len(X)), methods=[method])
+    second = run_reconstruction_benchmark(X, split=_split_for(len(X)), methods=[method])
+    assert dict(first.split_identity) == dict(second.split_identity)
+
+
+def test_report_rejects_malformed_payloads():
+    report = _example_report()
+    payload = report.to_dict()
+
+    missing = {key: value for key, value in payload.items() if key != "metrics"}
+    with pytest.raises(ValueError, match="missing required fields"):
+        ReconstructionBenchmarkReport.from_dict(missing)
+
+    extra = dict(payload)
+    extra["unexpected"] = 1
+    with pytest.raises(ValueError, match="unsupported fields"):
+        ReconstructionBenchmarkReport.from_dict(extra)
+
+    wrong_version = dict(payload)
+    wrong_version["schema_version"] = 99
+    with pytest.raises(ValueError, match="Unsupported benchmark report schema"):
+        ReconstructionBenchmarkReport.from_dict(wrong_version)
+
+    duplicated = copy.deepcopy(payload)
+    duplicated["results"].append(copy.deepcopy(duplicated["results"][0]))
+    with pytest.raises(ValueError, match="unique"):
+        ReconstructionBenchmarkReport.from_dict(duplicated)
+
+
+def test_result_rejects_non_finite_and_negative_values():
+    with pytest.raises(ValueError, match="must be finite"):
+        MethodBenchmarkResult(
+            name="bad",
+            method_type="pca",
+            latent_dimension=1,
+            uses_validation_partition=False,
+            original_scalars_per_sample=4,
+            latent_scalars_per_sample=1,
+            latent_dimensionality_ratio=0.25,
+            test_metrics={"mse": float("nan")},
+            fit_seconds=0.1,
+            reconstruction_seconds=0.1,
+        )
+
+    with pytest.raises(ValueError, match="must not be negative"):
+        MethodBenchmarkResult(
+            name="bad",
+            method_type="pca",
+            latent_dimension=1,
+            uses_validation_partition=False,
+            original_scalars_per_sample=4,
+            latent_scalars_per_sample=1,
+            latent_dimensionality_ratio=0.25,
+            test_metrics={"mse": 1.0},
+            fit_seconds=-1.0,
+            reconstruction_seconds=0.1,
+        )
+
+
+def test_latent_dimensionality_ratio_is_a_dimension_ratio():
+    X = _low_rank_dataset(sample_shape=(3, 4))
+    split = _split_for(len(X))
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[pca_benchmark_method("pca-k3", n_components=3)],
+    )
+    result = report.results[0]
+    assert result.original_scalars_per_sample == 12
+    assert result.latent_scalars_per_sample == 3
+    assert result.latent_dimensionality_ratio == pytest.approx(0.25)
+
+
+def test_result_rejects_an_inconsistent_latent_dimensionality_ratio():
+    with pytest.raises(ValueError, match="must equal"):
+        MethodBenchmarkResult(
+            name="misleading",
+            method_type="pca",
+            latent_dimension=3,
+            uses_validation_partition=False,
+            original_scalars_per_sample=12,
+            latent_scalars_per_sample=3,
+            latent_dimensionality_ratio=0.01,
+            test_metrics={"mse": 1.0},
+            fit_seconds=0.0,
+            reconstruction_seconds=0.0,
+        )
+
+
+# ---------------------------------------------------------------------------
+# N. Random-state isolation
+# ---------------------------------------------------------------------------
+
+
+def test_benchmark_does_not_leak_random_state_to_the_caller():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec(draw_random=True)
+
+    np.random.seed(1234)
+    torch.manual_seed(4321)
+    numpy_before = np.random.get_state()
+    torch_before = torch.random.get_rng_state().clone()
+
+    run_reconstruction_benchmark(X, split=split, methods=[specification], seed=99)
+
+    numpy_after = np.random.get_state()
+    torch_after = torch.random.get_rng_state()
+    assert numpy_before[0] == numpy_after[0]
+    assert np.array_equal(numpy_before[1], numpy_after[1])
+    assert numpy_before[2:] == numpy_after[2:]
+    assert torch.equal(torch_before, torch_after)
+
+
+def test_seeding_makes_the_random_draws_reproducible():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+
+    first_spec, first_created = _recording_spec(draw_random=True)
+    run_reconstruction_benchmark(X, split=split, methods=[first_spec], seed=7)
+    second_spec, second_created = _recording_spec(draw_random=True)
+    run_reconstruction_benchmark(X, split=split, methods=[second_spec], seed=7)
+
+    assert first_created[0].random_draws == second_created[0].random_draws
+
+
+def test_every_method_starts_from_the_same_seeded_state():
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    first, first_created = _recording_spec(name="first", draw_random=True)
+    second, second_created = _recording_spec(name="second", draw_random=True)
+
+    run_reconstruction_benchmark(X, split=split, methods=[first, second], seed=3)
+
+    assert first_created[0].random_draws == second_created[0].random_draws
+
+
+@pytest.mark.parametrize("seed", [-1, 1.5, True, "3", 2**32, 2**40])
+def test_invalid_seeds_are_rejected(seed):
+    X = _low_rank_dataset()
+    split = _split_for(len(X))
+    specification, _ = _recording_spec()
+
+    with pytest.raises((TypeError, ValueError)):
+        run_reconstruction_benchmark(
+            X,
+            split=split,
+            methods=[specification],
+            seed=seed,
+        )
+
+
+# ---------------------------------------------------------------------------
+# O. Tiny real autoencoder integration and VAE determinism
+# ---------------------------------------------------------------------------
+
+
+def test_benchmark_runs_against_a_real_small_autoencoder():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=2)
+    split = _split_for(len(X))
+
+    report = run_reconstruction_benchmark(
+        X,
+        split=split,
+        methods=[
+            pca_benchmark_method("pca-k2", n_components=2),
+            autoencoder_benchmark_method(
+                "standard-ae-k2",
+                model_factory=lambda: StandardAutoencoder(k=2, hidden_dims=[8]),
+                latent_dimension=2,
+                configuration={
+                    "architecture": "StandardAutoencoder",
+                    "hidden_dims": [8],
+                },
+                fit_kwargs={"epochs": 3, "batch_size": 8, "patience": 2},
+            ),
+        ],
+        seed=0,
+    )
+
+    assert [result.name for result in report.results] == ["pca-k2", "standard-ae-k2"]
+    for result in report.results:
+        assert set(result.test_metrics) == set(METRICS)
+        assert all(np.isfinite(value) for value in result.test_metrics.values())
+        assert result.latent_dimension == 2
+    assert report.results[0].uses_validation_partition is False
+    assert report.results[1].uses_validation_partition is True
+
+
+def test_variational_autoencoder_reconstruction_is_deterministic_by_default():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=2)
+    split = _split_for(len(X))
+    model = VariationalAutoencoder(k=2, hidden_dims=[8], beta=0.0)
+    adapter = AutoencoderReconstruction(model, latent_dimension=2)
+
+    adapter.fit(X[split.train_indices], X[split.validation_indices])
+    X_test = X[split.test_indices]
+    first = adapter.reconstruct(X_test)
+    second = adapter.reconstruct(X_test)
+
+    assert np.array_equal(first, second)
+
+
+def test_stochastic_variational_reconstruction_is_rejected():
+    model = VariationalAutoencoder(k=2, hidden_dims=[8])
+
+    with pytest.raises(ValueError, match="Stochastic reconstruction"):
+        AutoencoderReconstruction(
+            model,
+            latent_dimension=2,
+            predict_kwargs={"stochastic": True},
+        )
+
+
+@pytest.mark.parametrize(
+    "fit_kwargs",
+    [
+        {"validation_split": 0.3},
+        {"validation_data": (np.zeros((2, 6)), None)},
+        {"optimizer": object()},
+        {"criterion": object()},
+        {"y": np.zeros((2, 6))},
+    ],
+)
+def test_benchmark_controlled_fit_kwargs_are_rejected(fit_kwargs):
+    class _Model:
+        k = 2
+        is_fitted = False
+
+        def fit(self, X, validation_data=None, **kwargs):
+            return {}
+
+        def predict(self, X, **kwargs):
+            return X
+
+    with pytest.raises(ValueError, match="remove these fit_kwargs"):
+        AutoencoderReconstruction(
+            _Model(),
+            latent_dimension=2,
+            fit_kwargs=fit_kwargs,
+        )
+
+
+def test_autoencoder_adapter_requires_fit_and_predict():
+    with pytest.raises(TypeError, match="callable fit"):
+        AutoencoderReconstruction(object(), latent_dimension=2)
+
+
+# ---------------------------------------------------------------------------
+# P. Effective training configuration contributes to the identity
+# ---------------------------------------------------------------------------
+
+
+def _autoencoder_identity(**kwargs) -> str:
+    """Return the identity digest of a run differing only in the given kwargs."""
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=2)
+    split = _split_for(len(X))
+    specification = autoencoder_benchmark_method(
+        "ae",
+        model_factory=_JsonModel,
+        latent_dimension=2,
+        configuration={"architecture": "JsonModel"},
+        **kwargs,
+    )
+    report = run_reconstruction_benchmark(X, split=split, methods=[specification])
+    return report.identity_digest()
+
+
+@pytest.mark.parametrize(
+    "changed",
+    [
+        {"epochs": 7},
+        {"batch_size": 4},
+        {"learning_rate": 0.05},
+        {"patience": 3},
+        {"epochs": 7, "learning_rate": 0.05},
+    ],
+)
+def test_changing_effective_fit_kwargs_changes_the_identity(changed):
+    baseline = _autoencoder_identity(fit_kwargs={"epochs": 5})
+    altered = _autoencoder_identity(fit_kwargs={"epochs": 5, **changed})
+    assert altered != baseline
+
+
+def test_changing_predict_kwargs_changes_the_identity():
+    baseline = _autoencoder_identity(predict_kwargs={"batch_size": 64})
+    altered = _autoencoder_identity(predict_kwargs={"batch_size": 8})
+    assert altered != baseline
+
+
+def test_omitting_a_default_records_the_same_identity_as_passing_it():
+    explicit = _autoencoder_identity(
+        fit_kwargs=dict(_BENCHMARK_FIT_DEFAULTS),
+        predict_kwargs=dict(_BENCHMARK_PREDICT_DEFAULTS),
+    )
+    omitted = _autoencoder_identity()
+    assert explicit == omitted
+
+
+def test_verbosity_does_not_change_the_identity():
+    quiet = _autoencoder_identity(fit_kwargs={"epochs": 5, "verbose": 0})
+    loud = _autoencoder_identity(
+        fit_kwargs={"epochs": 5, "verbose": 2},
+        predict_kwargs={"verbose": 1},
+    )
+    assert quiet == loud
+
+
+def test_timing_does_not_change_the_identity():
+    first = _autoencoder_identity(fit_kwargs={"epochs": 5})
+    second = _autoencoder_identity(fit_kwargs={"epochs": 5})
+    assert first == second
+
+
+def test_effective_fit_configuration_is_recorded_in_the_report():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=2)
+    split = _split_for(len(X))
+    specification = autoencoder_benchmark_method(
+        "ae",
+        model_factory=_JsonModel,
+        latent_dimension=2,
+        configuration={"architecture": "JsonModel"},
+        fit_kwargs={"epochs": 5, "verbose": 0},
+        predict_kwargs={"batch_size": 8},
+    )
+    report = run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    recorded = report.results[0].configuration
+    assert recorded["architecture"] == "JsonModel"
+    assert dict(recorded["fit_kwargs"]) == {
+        "batch_size": 64,
+        "epochs": 5,
+        "learning_rate": 1e-3,
+        "patience": 20,
+    }
+    assert dict(recorded["predict_kwargs"]) == {"batch_size": 8}
+    # Verbosity is excluded from the identity but still reaches the model.
+    assert "verbose" not in recorded["fit_kwargs"]
+    payload = json.loads(report.to_json())
+    assert payload["results"][0]["configuration"]["fit_kwargs"]["epochs"] == 5
+
+
+def test_recorded_defaults_still_match_the_model_signatures():
+    fit_parameters = inspect.signature(BaseDeepLearningModel.fit).parameters
+    for name, value in _BENCHMARK_FIT_DEFAULTS.items():
+        assert fit_parameters[name].default == value, name
+
+    predict_parameters = inspect.signature(BaseDeepLearningModel.predict).parameters
+    for name, value in _BENCHMARK_PREDICT_DEFAULTS.items():
+        assert predict_parameters[name].default == value, name
+
+
+def test_verbosity_still_reaches_the_model_even_though_it_is_not_recorded():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=2)
+    split = _split_for(len(X))
+    built = []
+
+    def factory():
+        model = _JsonModel()
+        built.append(model)
+        return model
+
+    specification = autoencoder_benchmark_method(
+        "ae",
+        model_factory=factory,
+        latent_dimension=2,
+        fit_kwargs={"epochs": 5, "verbose": 3},
+    )
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    assert built[0].fit_calls[0]["verbose"] == 3
+    assert built[0].fit_calls[0]["epochs"] == 5
+    assert built[0].predict_calls[0]["verbose"] == 0
+
+
+@pytest.mark.parametrize(
+    "key", ["fit_kwargs", "predict_kwargs", "shuffle_training_data"]
+)
+def test_reserved_configuration_keys_are_rejected(key):
+    with pytest.raises(ValueError, match="reserved key"):
+        autoencoder_benchmark_method(
+            "clash",
+            model_factory=_JsonModel,
+            latent_dimension=2,
+            configuration={key: 1},
+        )
+
+
+@pytest.mark.parametrize(
+    "kwargs",
+    [
+        {"fit_kwargs": {"callback": len}},
+        {"fit_kwargs": {"nested": {"bad": object()}}},
+        {"predict_kwargs": {"hook": len}},
+        {"fit_kwargs": {"epochs": float("nan")}},
+    ],
+)
+def test_non_serializable_fit_or_predict_kwargs_are_rejected(kwargs):
+    with pytest.raises((TypeError, ValueError)):
+        autoencoder_benchmark_method(
+            "bad",
+            model_factory=_JsonModel,
+            latent_dimension=2,
+            **kwargs,
+        )
+
+
+@pytest.mark.parametrize("key", ["optimizer", "criterion", "validation_split", "y"])
+def test_forbidden_fit_kwargs_are_rejected_at_specification_time(key):
+    with pytest.raises(ValueError, match="remove these fit_kwargs"):
+        autoencoder_benchmark_method(
+            "bad",
+            model_factory=_JsonModel,
+            latent_dimension=2,
+            fit_kwargs={key: object()},
+        )
+
+
+def test_stochastic_prediction_is_rejected_at_specification_time():
+    with pytest.raises(ValueError, match="Stochastic reconstruction"):
+        autoencoder_benchmark_method(
+            "bad",
+            model_factory=_JsonModel,
+            latent_dimension=2,
+            predict_kwargs={"stochastic": True},
+        )
+
+
+def test_nested_fit_kwargs_cannot_be_mutated_after_specification():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=2)
+    split = _split_for(len(X))
+    nested = {"schedule": {"warmup": 5}}
+    caller_fit_kwargs = {"epochs": 5, "extra": nested}
+    built = []
+
+    def factory():
+        model = _JsonModel()
+        built.append(model)
+        return model
+
+    specification = autoencoder_benchmark_method(
+        "ae",
+        model_factory=factory,
+        latent_dimension=2,
+        fit_kwargs=caller_fit_kwargs,
+    )
+    before = run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    # Mutate the caller's mappings at every nesting depth after construction.
+    nested["schedule"]["warmup"] = 9999
+    nested["injected"] = True
+    caller_fit_kwargs["epochs"] = 1234
+
+    after = run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    assert before.identity_digest() == after.identity_digest()
+    assert built[0].fit_calls[0] == built[1].fit_calls[0]
+    assert built[1].fit_calls[0]["extra"] == {"schedule": {"warmup": 5}}
+    assert built[1].fit_calls[0]["epochs"] == 5
+
+
+def test_each_run_receives_independent_keyword_argument_containers():
+    X = _low_rank_dataset(n_samples=40, sample_shape=(6,), rank=2)
+    split = _split_for(len(X))
+    seen = []
+
+    class _MutatingModel(_JsonModel):
+        def fit(self, X_train, validation_data=None, **kwargs):
+            # Snapshot what this run was handed before vandalising it, so the
+            # assertion cannot be satisfied by this run's own mutation.
+            seen.append((kwargs, copy.deepcopy(kwargs)))
+            kwargs["extra"]["schedule"]["warmup"] = -1
+            self.is_fitted = True
+            return {}
+
+    specification = autoencoder_benchmark_method(
+        "ae",
+        model_factory=_MutatingModel,
+        latent_dimension=2,
+        fit_kwargs={"extra": {"schedule": {"warmup": 5}}},
+    )
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+    run_reconstruction_benchmark(X, split=split, methods=[specification])
+
+    first_kwargs, first_snapshot = seen[0]
+    second_kwargs, second_snapshot = seen[1]
+    assert first_kwargs is not second_kwargs
+    assert first_kwargs["extra"] is not second_kwargs["extra"]
+    assert first_kwargs["extra"]["schedule"] is not second_kwargs["extra"]["schedule"]
+    # Run 1 zeroed its own container; run 2 must still start from the spec.
+    assert first_snapshot["extra"]["schedule"]["warmup"] == 5
+    assert second_snapshot["extra"]["schedule"]["warmup"] == 5
+    assert first_kwargs["extra"]["schedule"]["warmup"] == -1
+
+
+# ---------------------------------------------------------------------------
+# Q. Multi-device random-state isolation
+# ---------------------------------------------------------------------------
+
+
+def _fork_rng_spy(monkeypatch):
+    """Record the devices requested from fork_rng, forking only the CPU."""
+    recorded = {}
+    real_fork_rng = torch.random.fork_rng
+
+    def spy(*args, **kwargs):
+        devices = kwargs.get("devices", args[0] if args else None)
+        recorded["devices"] = None if devices is None else list(devices)
+        return real_fork_rng(devices=[])
+
+    monkeypatch.setattr(torch.random, "fork_rng", spy)
+    return recorded
+
+
+def test_every_visible_cuda_device_is_forked(monkeypatch):
+    recorded = _fork_rng_spy(monkeypatch)
+    monkeypatch.setattr(torch.cuda, "is_available", lambda: True)
+    monkeypatch.setattr(torch.cuda, "device_count", lambda: 3)
+
+    with _isolated_random_state(11):
+        pass
+
+    # torch.manual_seed reseeds every visible device, so every visible device
+    # must be forked, not only the current one.
+    assert recorded["devices"] == [0, 1, 2]
+
+
+def test_no_cuda_device_is_touched_when_cuda_is_unavailable(monkeypatch):
+    recorded = _fork_rng_spy(monkeypatch)
+    monkeypatch.setattr(torch.cuda, "is_available", lambda: False)
+
+    def _forbidden():
+        raise AssertionError("device_count must not be called without CUDA.")
+
+    monkeypatch.setattr(torch.cuda, "device_count", _forbidden)
+
+    with _isolated_random_state(11):
+        pass
+
+    assert recorded["devices"] == []
+
+
+def test_isolated_random_state_restores_cpu_generators():
+    np.random.seed(5)
+    torch.manual_seed(6)
+    numpy_before = np.random.get_state()
+    torch_before = torch.random.get_rng_state().clone()
+
+    with _isolated_random_state(77):
+        np.random.rand(10)
+        torch.rand(10)
+
+    assert np.array_equal(numpy_before[1], np.random.get_state()[1])
+    assert torch.equal(torch_before, torch.random.get_rng_state())
+
+
+# ---------------------------------------------------------------------------
+# R. Serialized cross-field validation
+# ---------------------------------------------------------------------------
+
+
+def test_partition_sizes_must_sum_to_the_sample_count():
+    payload = _example_report().to_dict()
+    payload["partition_sizes"]["excluded"] = 7
+
+    with pytest.raises(ValueError, match="partition_sizes sum to"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+@pytest.mark.parametrize("partition", ["train", "validation", "test"])
+def test_scientific_partitions_must_not_be_empty(partition):
+    payload = _example_report().to_dict()
+    payload["partition_sizes"]["excluded"] += payload["partition_sizes"][partition]
+    payload["partition_sizes"][partition] = 0
+
+    with pytest.raises(ValueError, match="at least one\\s+sample"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+def test_split_identity_sample_count_must_agree_with_the_report():
+    payload = _example_report().to_dict()
+    payload["split_identity"]["n_samples"] = 999
+
+    with pytest.raises(ValueError, match="split_identity records"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+def test_split_identity_must_record_a_sample_count():
+    payload = _example_report().to_dict()
+    del payload["split_identity"]["n_samples"]
+
+    with pytest.raises(ValueError, match="must record n_samples"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+@pytest.mark.parametrize("field", ["time_axis_fingerprint", "partition_digest"])
+def test_split_identity_digests_must_be_well_formed(field):
+    payload = _example_report().to_dict()
+    payload["split_identity"][field] = "nope"
+
+    with pytest.raises(ValueError, match="SHA-256 hex digest"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+def test_original_scalars_must_agree_with_the_sample_shape():
+    payload = _example_report().to_dict()
+    payload["results"][0]["original_scalars_per_sample"] = 11
+    payload["results"][0]["latent_dimensionality_ratio"] = (
+        payload["results"][0]["latent_scalars_per_sample"] / 11
+    )
+
+    with pytest.raises(ValueError, match="scalars per sample"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+def test_latent_dimension_must_agree_with_latent_scalars():
+    payload = _example_report().to_dict()
+    payload["results"][0]["latent_dimension"] = 5
+
+    with pytest.raises(ValueError, match="must equal\\s+latent_scalars_per_sample"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+@pytest.mark.parametrize("metric", ["mse", "mae", "rmse"])
+def test_negative_reconstruction_metrics_are_rejected(metric):
+    payload = _example_report().to_dict()
+    payload["results"][0]["test_metrics"][metric] = -0.5
+
+    with pytest.raises(ValueError, match="must not be negative"):
+        ReconstructionBenchmarkReport.from_dict(payload)
+
+
+def test_a_consistent_report_still_round_trips_after_the_new_checks():
+    report = _example_report()
+    rebuilt = ReconstructionBenchmarkReport.from_dict(report.to_dict())
+    assert rebuilt.to_dict() == report.to_dict()
+    assert rebuilt.identity_digest() == report.identity_digest()
diff --git a/tests/deeplearning/test_explicit_validation_data.py b/tests/deeplearning/test_explicit_validation_data.py
new file mode 100644
index 0000000..dbc4515
--- /dev/null
+++ b/tests/deeplearning/test_explicit_validation_data.py
@@ -0,0 +1,333 @@
+"""Regression tests for explicit autoencoder validation data.
+
+``validation_data`` was added so that chronological validation membership can
+be supplied exactly. These tests pin both the new semantics and the unchanged
+behaviour of the historical ``validation_split`` path.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+torch = pytest.importorskip("torch")
+
+from bluemath_tk.deeplearning.autoencoders import (  # noqa: E402
+    ConvLSTMAutoencoder,
+    OrthogonalAutoencoder,
+    StandardAutoencoder,
+)
+from bluemath_tk.deeplearning.spatiotemporal_autoencoders import (  # noqa: E402
+    SpatialTokenConvLSTMTransformerAutoencoder,
+)
+from bluemath_tk.deeplearning.variational_autoencoders import (  # noqa: E402
+    VariationalAutoencoder,
+)
+
+
+def _dataset(n_samples: int = 40, n_features: int = 6, seed: int = 0) -> np.ndarray:
+    rng = np.random.default_rng(seed)
+    latent = rng.normal(size=(n_samples, 2))
+    mixing = rng.normal(size=(2, n_features))
+    return latent @ mixing
+
+
+def _sequence_dataset(n_samples: int = 12, seed: int = 0) -> np.ndarray:
+    rng = np.random.default_rng(seed)
+    return rng.normal(size=(n_samples, 2, 1, 4, 4))
+
+
+class _PartitionRecorder:
+    """Mixin capturing exactly what the resolved fit partitions contain."""
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.resolved = None
+
+    def _resolve_fit_partitions(self, X, y, validation_split, validation_data):
+        resolved = super()._resolve_fit_partitions(
+            X,
+            y,
+            validation_split,
+            validation_data,
+        )
+        self.resolved = tuple(np.array(part, copy=True) for part in resolved)
+        return resolved
+
+
+class _RecordingStandardAutoencoder(_PartitionRecorder, StandardAutoencoder):
+    pass
+
+
+class _RecordingOrthogonalAutoencoder(_PartitionRecorder, OrthogonalAutoencoder):
+    pass
+
+
+class _RecordingVariationalAutoencoder(_PartitionRecorder, VariationalAutoencoder):
+    pass
+
+
+RECORDING_MODELS = [
+    _RecordingStandardAutoencoder,
+    _RecordingOrthogonalAutoencoder,
+    _RecordingVariationalAutoencoder,
+]
+
+
+@pytest.mark.parametrize("model_class", RECORDING_MODELS)
+def test_explicit_validation_data_uses_exactly_the_supplied_samples(model_class):
+    X = _dataset()
+    X_train = X[:30]
+    X_validation = X[30:36]
+    X_test = X[36:]
+
+    model = model_class(k=2, hidden_dims=[8])
+    model.fit(
+        X_train,
+        validation_data=(X_validation, None),
+        epochs=2,
+        batch_size=8,
+        patience=2,
+        verbose=0,
+    )
+
+    resolved_train, resolved_train_y, resolved_val, resolved_val_y = model.resolved
+    assert np.array_equal(resolved_train, X_train)
+    assert np.array_equal(resolved_train_y, X_train)
+    assert np.array_equal(resolved_val, X_validation)
+    assert np.array_equal(resolved_val_y, X_validation)
+
+    # No test sample reaches optimisation or validation.
+    test_rows = {row.tobytes() for row in X_test}
+    for array in (resolved_train, resolved_val):
+        for row in array:
+            assert row.tobytes() not in test_rows
+
+
+@pytest.mark.parametrize("model_class", RECORDING_MODELS)
+def test_explicit_validation_data_preserves_training_order(model_class):
+    X = _dataset()
+    model = model_class(k=2, hidden_dims=[8])
+    model.fit(
+        X[:30],
+        validation_data=(X[30:36], None),
+        epochs=2,
+        batch_size=8,
+        patience=2,
+        verbose=0,
+    )
+    assert np.array_equal(model.resolved[0], X[:30])
+
+
+@pytest.mark.parametrize("model_class", RECORDING_MODELS)
+def test_explicit_validation_data_leaves_the_global_numpy_state_untouched(model_class):
+    X = _dataset()
+    model = model_class(k=2, hidden_dims=[8])
+
+    np.random.seed(17)
+    state_before = np.random.get_state()
+    model.fit(
+        X[:30],
+        validation_data=(X[30:36], None),
+        epochs=2,
+        batch_size=8,
+        patience=2,
+        verbose=0,
+    )
+    state_after = np.random.get_state()
+
+    assert state_before[0] == state_after[0]
+    assert np.array_equal(state_before[1], state_after[1])
+    assert state_before[2:] == state_after[2:]
+
+
+@pytest.mark.parametrize("model_class", RECORDING_MODELS)
+def test_validation_split_path_is_unchanged_and_still_shuffles(model_class):
+    X = _dataset()
+    model = model_class(k=2, hidden_dims=[8])
+
+    np.random.seed(3)
+    state_before = np.random.get_state()
+    model.fit(X, validation_split=0.25, epochs=2, batch_size=8, patience=2, verbose=0)
+    state_after = np.random.get_state()
+
+    resolved_train, _, resolved_val, _ = model.resolved
+    assert len(resolved_train) == int(0.75 * len(X))
+    assert len(resolved_val) == len(X) - int(0.75 * len(X))
+    # The historical path consumes the global random state.
+    assert not np.array_equal(state_before[1], state_after[1])
+    # The random split is not the chronological tail.
+    assert not np.array_equal(resolved_val, X[-len(resolved_val) :])
+
+
+def test_validation_split_selection_is_reproducible_for_a_fixed_seed():
+    X = _dataset()
+
+    def _resolved_validation():
+        model = _RecordingStandardAutoencoder(k=2, hidden_dims=[8])
+        np.random.seed(101)
+        model.fit(
+            X,
+            validation_split=0.2,
+            epochs=1,
+            batch_size=8,
+            patience=1,
+            verbose=0,
+        )
+        return model.resolved[2]
+
+    assert np.array_equal(_resolved_validation(), _resolved_validation())
+
+
+def test_validation_data_takes_precedence_over_validation_split():
+    X = _dataset()
+    model = _RecordingStandardAutoencoder(k=2, hidden_dims=[8])
+    model.fit(
+        X[:30],
+        validation_split=0.9,
+        validation_data=(X[30:34], None),
+        epochs=1,
+        batch_size=8,
+        patience=1,
+        verbose=0,
+    )
+    assert np.array_equal(model.resolved[0], X[:30])
+    assert np.array_equal(model.resolved[2], X[30:34])
+
+
+def test_explicit_validation_targets_are_honoured():
+    X = _dataset()
+    targets = X[30:34] * 2.0
+    model = _RecordingStandardAutoencoder(k=2, hidden_dims=[8])
+    model.fit(
+        X[:30],
+        validation_data=(X[30:34], targets),
+        epochs=1,
+        batch_size=8,
+        patience=1,
+        verbose=0,
+    )
+    assert np.array_equal(model.resolved[3], targets)
+
+
+@pytest.mark.parametrize(
+    ("validation_data", "error", "message"),
+    [
+        ([1, 2], TypeError, "must be an"),
+        ((1, 2, 3), TypeError, "must be an"),
+        (("not-an-array", None), TypeError, "must be a NumPy array"),
+        ((np.empty((0, 6)), None), ValueError, "at least one sample"),
+        ((np.zeros((4, 5)), None), ValueError, "per-sample shape"),
+        ((np.zeros((4, 6, 1)), None), ValueError, "dimensions to match X"),
+        ((np.full((4, 6), np.nan), None), ValueError, "only finite values"),
+        ((np.full((4, 6), np.inf), None), ValueError, "only finite values"),
+        ((np.zeros((4, 6)), np.zeros((3, 6))), ValueError, "same number of samples"),
+        ((np.zeros((4, 6)), np.zeros((4, 5))), ValueError, "incompatible"),
+        ((np.zeros((4, 6)), "not-an-array"), TypeError, "NumPy array or None"),
+    ],
+)
+def test_invalid_validation_data_is_rejected(validation_data, error, message):
+    X = _dataset()
+    model = StandardAutoencoder(k=2, hidden_dims=[8])
+
+    with pytest.raises(error, match=message):
+        model.fit(
+            X[:30],
+            validation_data=validation_data,
+            epochs=1,
+            batch_size=8,
+            patience=1,
+            verbose=0,
+        )
+
+
+def test_explicit_validation_data_requires_at_least_two_training_samples():
+    X = _dataset()
+    model = StandardAutoencoder(k=2, hidden_dims=[8])
+
+    with pytest.raises(ValueError, match="at least two training"):
+        model.fit(
+            X[:1],
+            validation_data=(X[30:34], None),
+            epochs=1,
+            batch_size=8,
+            patience=1,
+            verbose=0,
+        )
+
+
+def test_validation_split_is_still_validated_when_no_validation_data_is_given():
+    X = _dataset()
+    model = StandardAutoencoder(k=2, hidden_dims=[8])
+
+    with pytest.raises(ValueError, match="strictly between 0 and 1"):
+        model.fit(X, validation_split=1.5, epochs=1, batch_size=8, verbose=0)
+
+
+def test_validation_split_bounds_are_ignored_with_explicit_validation_data():
+    X = _dataset()
+    model = _RecordingStandardAutoencoder(k=2, hidden_dims=[8])
+
+    model.fit(
+        X[:30],
+        validation_split=0.0,
+        validation_data=(X[30:34], None),
+        epochs=1,
+        batch_size=8,
+        patience=1,
+        verbose=0,
+    )
+    assert np.array_equal(model.resolved[2], X[30:34])
+
+
+def test_sequence_models_forward_validation_data_through_their_wrappers():
+    X = _sequence_dataset()
+
+    builders = (
+        lambda: ConvLSTMAutoencoder(k=2),
+        lambda: SpatialTokenConvLSTMTransformerAutoencoder(
+            k=2,
+            spatial_pool_size=(1, 1),
+            d_model=8,
+            n_heads=2,
+            n_layers=1,
+        ),
+    )
+    for builder in builders:
+        captured = {}
+        model = builder()
+        original = model._resolve_fit_partitions
+
+        def _spy(X_fit, y, validation_split, validation_data, _original=original):
+            resolved = _original(X_fit, y, validation_split, validation_data)
+            captured["resolved"] = tuple(np.array(part, copy=True) for part in resolved)
+            return resolved
+
+        model._resolve_fit_partitions = _spy
+        model.fit(
+            X[:8],
+            validation_data=(X[8:10], None),
+            epochs=1,
+            batch_size=4,
+            patience=1,
+            verbose=0,
+        )
+        assert np.array_equal(captured["resolved"][0], X[:8])
+        assert np.array_equal(captured["resolved"][2], X[8:10])
+
+
+def test_fitting_twice_with_explicit_validation_data_stays_consistent():
+    X = _dataset()
+    model = StandardAutoencoder(k=2, hidden_dims=[8])
+    for _ in range(2):
+        history = model.fit(
+            X[:30],
+            validation_data=(X[30:36], None),
+            epochs=2,
+            batch_size=8,
+            patience=2,
+            verbose=0,
+        )
+        assert len(history["val_loss"]) >= 1
+        assert all(np.isfinite(value) for value in history["val_loss"])
+    assert model.is_fitted is True