From 48ff06dd0f5cfc9c9fa1b371368f39f7e8f14bb2 Mon Sep 17 00:00:00 2001 From: Sergio Date: Sun, 16 Aug 2026 06:21:42 +0200 Subject: [PATCH 1/6] Support explicit autoencoder validation data BlueMath autoencoders derived their validation set by shuffling samples and cutting at validation_split, so the validation partition could not be stated exactly. That is incompatible with chronological validation, where membership must come from the split rather than from a random permutation. fit() now accepts an optional validation_data=(X_validation, y_validation) pair. When 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. Partition resolution moves into one shared _resolve_fit_partitions helper so the base model, OrthogonalAutoencoder, and VariationalAutoencoder cannot drift apart. The sequence models forward the new argument explicitly. Omitting validation_data preserves the historical behaviour exactly, including its use of the global random state. Co-Authored-By: Claude Opus 5 --- bluemath_tk/deeplearning/_base_model.py | 137 ++++++- bluemath_tk/deeplearning/autoencoders.py | 34 +- .../spatiotemporal_autoencoders.py | 2 + .../deeplearning/variational_autoencoders.py | 29 +- .../test_explicit_validation_data.py | 333 ++++++++++++++++++ 5 files changed, 505 insertions(+), 30 deletions(-) create mode 100644 tests/deeplearning/test_explicit_validation_data.py diff --git a/bluemath_tk/deeplearning/_base_model.py b/bluemath_tk/deeplearning/_base_model.py index 921c4dc..c853ebd 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,15 @@ 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: + if len(X) < 2: + raise ValueError( + "Explicit validation_data requires at least two training " + "samples in X." + ) + self._validate_validation_data(X, validation_data) + return + split = int((1 - validation_split) * len(X)) if split < 2: raise ValueError( @@ -283,6 +293,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 +637,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 +665,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 +686,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/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 From bc8b1438a618e97281b81961d1217a3fc0b71dae Mon Sep 17 00:00:00 2001 From: Sergio Date: Sun, 16 Aug 2026 06:21:55 +0200 Subject: [PATCH 2/6] Add common reconstruction benchmark framework Adds bluemath_tk.benchmarking, infrastructure for comparing dimensionality reduction methods on exactly the same held-out samples. Partition membership always comes from a manifest-backed ChronologicalSplit; the runner never creates a split of its own. PCA is fitted on the training partition alone, autoencoders receive the validation partition for early stopping only, and metrics are computed on the test partition only, in the original sample space, using the accepted implementation in bluemath_tk.deeplearning.metrics. Methods reach the runner through one small protocol, so further reconstruction models can be added without touching the scientific core. PCAReconstruction wraps the existing bluemath_tk.datamining.pca.PCA rather than introducing a second PCA; stacking and the inverse reshape both use C order, so sample order and per-sample shape survive the round trip. AutoencoderReconstruction uses only the accepted public model workflow and rejects stochastic variational reconstruction, which would not be comparable with deterministic PCA and autoencoder reconstructions. Methods are supplied as specifications carrying a factory, so a fresh unfitted model is built per run and fitted state cannot leak between runs. Factories are never introspected: the reproducible description of a method comes from explicit, JSON-compatible method_type and configuration fields. The report serializes deterministically. Its identity digest covers what was compared and deliberately excludes measured outcomes, because metric values and wall-clock timings are observational rather than reproducible. Reported latent dimensionality is a dimension ratio, not a bitrate or a storage compression ratio. The metrics import is deferred so that importing bluemath_tk still works without the optional PyTorch dependency. Co-Authored-By: Claude Opus 5 --- bluemath_tk/__init__.py | 2 + bluemath_tk/benchmarking/__init__.py | 25 + bluemath_tk/benchmarking/reconstruction.py | 1352 ++++++++++++++++++++ 3 files changed, 1379 insertions(+) create mode 100644 bluemath_tk/benchmarking/__init__.py create mode 100644 bluemath_tk/benchmarking/reconstruction.py 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..9c79ba7 --- /dev/null +++ b/bluemath_tk/benchmarking/reconstruction.py @@ -0,0 +1,1352 @@ +"""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 math +from collections.abc import Callable, Iterator, Mapping, Sequence +from contextlib import ExitStack, contextmanager +from dataclasses import dataclass, field +from time import perf_counter +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") + + +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 _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)) + available = min(int(X_train.shape[0]), 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) + + +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. + fit_kwargs : dict, optional + Extra keyword arguments forwarded to ``model.fit``. ``validation_data`` + and ``validation_split`` are rejected because the benchmark controls + partition membership. ``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 = ("X", "y", "validation_data", "validation_split") + + def __init__( + self, + model: Any, + latent_dimension: int, + *, + 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." + ) + self.model = model + 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)}." + ) + + self._fit_kwargs = dict(fit_kwargs or {}) + forbidden = sorted( + set(self._FORBIDDEN_FIT_KWARGS).intersection(self._fit_kwargs) + ) + if forbidden: + raise ValueError( + "The benchmark controls partition membership; remove these " + f"fit_kwargs: {forbidden}." + ) + self._fit_kwargs.setdefault("verbose", 0) + + self._predict_kwargs = dict(predict_kwargs or {}) + if self._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." + ) + 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. + X_validation : np.ndarray + Validation samples. Exactly these samples drive the validation loss + and early stopping. + """ + _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") + 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, + 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. + This is recorded verbatim; the factory itself is never introspected. + fit_kwargs : mapping, optional + Extra keyword arguments for ``model.fit``. + predict_kwargs : mapping, optional + Extra keyword arguments for ``model.predict``. + + 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) + frozen_fit_kwargs = dict(fit_kwargs or {}) + frozen_predict_kwargs = dict(predict_kwargs or {}) + + def factory() -> ReconstructionMethod: + return AutoencoderReconstruction( + model_factory(), + latent_dimension=width, + fit_kwargs=frozen_fit_kwargs, + predict_kwargs=frozen_predict_kwargs, + ) + + return BenchmarkMethod( + name=name, + method_type="autoencoder", + latent_dimension=width, + factory=factory, + configuration=configuration, + 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)), + ) + 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}.") + metrics[metric] = _validate_finite_float(f"test_metrics[{metric!r}]", value) + 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. + 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, ...] + 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 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 + 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.") + + object.__setattr__( + self, + "split_identity", + _freeze_json( + _validate_configuration(self.split_identity, name="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.") + 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)}." + ) + 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), + "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), + "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", + "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"]), + 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 + + +@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. + + 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: + devices: list[int] = [] + if torch.cuda.is_available(): # pragma: no cover - needs CUDA + current = torch.cuda.current_device() + devices = [current] + stack.enter_context(torch.random.fork_rng(devices=devices)) + 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 _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, + "dataset_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 seed < 0: + raise ValueError("seed must be non-negative.") + + 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. + 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) + + fit_start = perf_counter() + instance.fit( + X_train, + X_validation if specification.uses_validation_partition else None, + ) + fit_seconds = perf_counter() - fit_start + + reconstruction_start = perf_counter() + reconstruction = instance.reconstruct(X_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, + 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), + ) From 98a6ff6a0ddd0f7acd36787e66d8f7c20d7aec89 Mon Sep 17 00:00:00 2001 From: Sergio Date: Sun, 16 Aug 2026 06:22:06 +0200 Subject: [PATCH 3/6] Add PCA and autoencoder benchmark regression tests Covers the fairness properties the framework claims, not just its outputs. Leakage is attacked structurally rather than through final scores: test values are altered by six orders of magnitude and the fitted PCA components, mean, explained variance, and stacked training matrix must stay bit-identical, and controlled recording methods prove that only training samples reach fitting and only the chronological validation samples reach validation. Also covers exact split membership and partition boundaries, PCA correctness on low-rank data, shape and sample-order round trips for 1D to 4D samples, independence from C or Fortran memory layout, exact agreement with the shared metric implementation, rejection of broadcastable reconstruction shapes, non-finite and non-numeric data, impossible latent dimensions, duplicate method names, input mutation, random-state isolation and reseeding, deterministic report serialization with timing excluded from identity, and a small real StandardAutoencoder plus VariationalAutoencoder integration. Co-Authored-By: Claude Opus 5 --- tests/benchmarking/test_reconstruction.py | 1337 +++++++++++++++++++++ 1 file changed, 1337 insertions(+) create mode 100644 tests/benchmarking/test_reconstruction.py diff --git a/tests/benchmarking/test_reconstruction.py b/tests/benchmarking/test_reconstruction.py new file mode 100644 index 0000000..7f9784f --- /dev/null +++ b/tests/benchmarking/test_reconstruction.py @@ -0,0 +1,1337 @@ +"""Regression tests for the common reconstruction benchmark framework.""" + +from __future__ import annotations + +import copy +import hashlib +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.deeplearning.autoencoders import StandardAutoencoder # noqa: E402 +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") + + +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 + + +# --------------------------------------------------------------------------- +# 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"}, + ) + 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): + 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]) + + +# --------------------------------------------------------------------------- +# 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) + run_reconstruction_benchmark(X, split=split, methods=[specification]) + + assert np.array_equal(X, X_before) + + +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)], + ) + + +def test_pca_rejects_more_components_than_training_samples(): + X = _low_rank_dataset(n_samples=12, sample_shape=(20,), rank=3) + split = _split_for(len(X), fractions=(0.5, 0.25, 0.25)) + + with pytest.raises(ValueError, match="components available"): + run_reconstruction_benchmark( + X, + split=split, + methods=[pca_benchmark_method("too-many", n_components=10)], + ) + + +def test_autoencoder_latent_width_must_match_the_specification(): + class _Model: + k = 3 + is_fitted = False + + def fit(self, X, **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="= 0.0 + assert result.reconstruction_seconds >= 0.0 + + +# --------------------------------------------------------------------------- +# 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_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"]) +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}, + ) + + +def test_benchmark_controlled_fit_kwargs_are_rejected(): + class _Model: + k = 2 + is_fitted = False + + def fit(self, X, **kwargs): + return {} + + def predict(self, X, **kwargs): + return X + + with pytest.raises(ValueError, match="remove these"): + AutoencoderReconstruction( + _Model(), + latent_dimension=2, + fit_kwargs={"validation_split": 0.3}, + ) + + +def test_autoencoder_adapter_requires_fit_and_predict(): + with pytest.raises(TypeError, match="callable fit"): + AutoencoderReconstruction(object(), latent_dimension=2) From 034cff99ad9ae2636f84d10c8bd1b1668dcbdfe8 Mon Sep 17 00:00:00 2001 From: Sergio Date: Sun, 16 Aug 2026 06:22:15 +0200 Subject: [PATCH 4/6] Document PCA and autoencoder benchmarking Explains what the framework measures and, as importantly, what it does not: reconstruction error on held-out samples is not downstream scientific skill, and the reported latent dimensionality ratio is not a bitrate or a storage compression ratio. Documents why chronological partitions matter, how PCA and autoencoders are compared, the one deliberate asymmetry between them, the shared metrics and the BlueMath RMSE convention, how to supply a model through a factory, the explicit validation-data contract, how the test partition stays isolated, deterministic report identity, and current limitations. docs/source/benchmarking.rst is force-added because .gitignore ignores *.rst; docs/source/validation.rst is tracked the same way. Co-Authored-By: Claude Opus 5 --- docs/source/benchmarking.rst | 271 +++++++++++++++++++++++++++++++++++ docs/source/index.rst | 1 + 2 files changed, 272 insertions(+) create mode 100644 docs/source/benchmarking.rst diff --git a/docs/source/benchmarking.rst b/docs/source/benchmarking.rst new file mode 100644 index 0000000..fccc2bf --- /dev/null +++ b/docs/source/benchmarking.rst @@ -0,0 +1,271 @@ +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. + +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 treated asymmetrically in exactly one respect, and the +report records it explicitly through ``uses_validation_partition``: + +* 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. + +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`. + +Common metrics +-------------- + +Metrics are computed with the accepted implementation in +:mod:`bluemath_tk.deeplearning.metrics`, so benchmark numbers agree exactly with +``reconstruction_error`` and with the per-model ``evaluate_reconstruction`` +methods. The available metrics are ``mse``, ``mae``, and ``rmse``, all reported +with ``reduction="mean"``. + +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, 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. + +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. + +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, so peak memory includes one copy of the + train, validation, and test partitions. +* 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 ================== From 8b7ae0a7e19c4db7e05dc6c5a58355f78059ff5d Mon Sep 17 00:00:00 2001 From: Sergio Date: Sun, 16 Aug 2026 06:57:55 +0200 Subject: [PATCH 5/6] Close benchmark fairness gaps found in adversarial review Five defects, each with a regression test that fails without the fix. Autoencoders were trained on contiguous chronological mini-batches, because explicit validation data deliberately preserves the supplied order. PCA has a closed-form solution and pays no such cost, so this was an unrecorded second asymmetry between the two families. Training rows are now permuted once inside the isolated random state, which leaves partition membership untouched, is controlled by the run seed, and is recorded as shuffle_training_data. All methods were handed the same partition arrays, so a method writing to its inputs corrupted every later method and the test partition metrics were scored against it. Each method now receives its own copies, and metrics are always scored against a pristine test partition. fit_kwargs accepted optimizer and criterion. Since a fresh model is built per run, a PyTorch optimizer bound to another model's parameters silently skips them, reporting an untrained network as a legitimate result. Both are now refused. identity_digest() covered only the time axis, so two entirely different datasets with the same shape and coordinates produced the same digest. The report now carries a data_digest over the values, dtype and shape, normalised to C order; the manifest fingerprint is renamed time_axis_fingerprint to say what it actually covers. A model with a fit(self, X, **kwargs) signature absorbed validation_data silently and remained free to build its own random split. An explicit validation_data parameter is now required. Also: PCA rejects component counts above the centered training rank, which is n_train - 1 rather than n_train, so a numerically null component can no longer inflate the reported latent dimensionality; validation data is validated once instead of twice per fit; and two weak tests are replaced, one of which passed while the corruption it was named for was actually occurring. Co-Authored-By: Claude Opus 5 --- bluemath_tk/benchmarking/reconstruction.py | 140 +++++++++- bluemath_tk/deeplearning/_base_model.py | 4 +- docs/source/benchmarking.rst | 91 +++++- tests/benchmarking/test_reconstruction.py | 306 ++++++++++++++++++++- 4 files changed, 502 insertions(+), 39 deletions(-) diff --git a/bluemath_tk/benchmarking/reconstruction.py b/bluemath_tk/benchmarking/reconstruction.py index 9c79ba7..841afc0 100644 --- a/bluemath_tk/benchmarking/reconstruction.py +++ b/bluemath_tk/benchmarking/reconstruction.py @@ -16,7 +16,9 @@ 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 @@ -280,7 +282,10 @@ def fit(self, X_train: np.ndarray, X_validation: np.ndarray | None = None) -> No _validate_sample_data(X_train, name="X_train") sample_shape = tuple(X_train.shape[1:]) n_features = int(np.prod(sample_shape)) - available = min(int(X_train.shape[0]), n_features) + # 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} " @@ -316,6 +321,28 @@ def reconstruct(self, X: np.ndarray) -> np.ndarray: return np.asarray(reconstructed[_PCA_VARIABLE].values, dtype=np.float64) +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. @@ -330,10 +357,19 @@ class AutoencoderReconstruction: 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. ``verbose`` defaults to 0. + 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. @@ -346,13 +382,21 @@ class AutoencoderReconstruction: autoencoder reconstructions is not a like-for-like measurement. """ - _FORBIDDEN_FIT_KWARGS = ("X", "y", "validation_data", "validation_split") + _FORBIDDEN_FIT_KWARGS = ( + "X", + "y", + "criterion", + "optimizer", + "validation_data", + "validation_split", + ) 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, ): @@ -362,7 +406,12 @@ def __init__( 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, @@ -380,8 +429,8 @@ def __init__( ) if forbidden: raise ValueError( - "The benchmark controls partition membership; remove these " - f"fit_kwargs: {forbidden}." + "The benchmark controls partition membership and builds a fresh " + f"model for every run; remove these fit_kwargs: {forbidden}." ) self._fit_kwargs.setdefault("verbose", 0) @@ -416,10 +465,12 @@ def fit(self, X_train: np.ndarray, X_validation: np.ndarray | None = None) -> No Parameters ---------- X_train : np.ndarray - Training samples. Every one of them is used for optimisation. + 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. + and early stopping. Their membership is never changed by shuffling. """ _validate_sample_data(X_train, name="X_train") if X_validation is None: @@ -428,6 +479,8 @@ def fit(self, X_train: np.ndarray, X_validation: np.ndarray | None = None) -> No "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), @@ -557,6 +610,7 @@ def autoencoder_benchmark_method( 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: @@ -576,6 +630,12 @@ def autoencoder_benchmark_method( configuration : mapping, optional JSON-compatible description of the architecture and hyperparameters. This is recorded verbatim; the factory itself is never introspected. + The reserved key ``"shuffle_training_data"`` is added 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``. predict_kwargs : mapping, optional @@ -590,13 +650,23 @@ def autoencoder_benchmark_method( 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) frozen_fit_kwargs = dict(fit_kwargs or {}) frozen_predict_kwargs = dict(predict_kwargs or {}) + recorded = _validate_configuration(configuration, name="configuration") + if "shuffle_training_data" in recorded: + raise ValueError( + "configuration must not set the reserved key " + "'shuffle_training_data'; use the shuffle_training_data argument." + ) + recorded["shuffle_training_data"] = shuffle + def factory() -> ReconstructionMethod: return AutoencoderReconstruction( model_factory(), latent_dimension=width, + shuffle_training_data=shuffle, fit_kwargs=frozen_fit_kwargs, predict_kwargs=frozen_predict_kwargs, ) @@ -606,7 +676,7 @@ def factory() -> ReconstructionMethod: method_type="autoencoder", latent_dimension=width, factory=factory, - configuration=configuration, + configuration=recorded, uses_validation_partition=True, ) @@ -822,6 +892,10 @@ class ReconstructionBenchmarkReport: 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 @@ -844,6 +918,7 @@ class ReconstructionBenchmarkReport: schema_version: int n_samples: int sample_shape: tuple[int, ...] + data_digest: str partition_sizes: Mapping[str, int] metrics: tuple[str, ...] seed: int | None @@ -879,6 +954,13 @@ def __post_init__(self) -> None: 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): @@ -941,6 +1023,7 @@ def to_dict(self) -> dict[str, JsonValue]: "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, @@ -964,6 +1047,7 @@ def identity(self) -> dict[str, JsonValue]: "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, @@ -985,6 +1069,7 @@ def from_dict(cls, payload: Mapping[str, Any]) -> ReconstructionBenchmarkReport: "schema_version", "n_samples", "sample_shape", + "data_digest", "partition_sizes", "metrics", "seed", @@ -1002,6 +1087,7 @@ def from_dict(cls, payload: Mapping[str, Any]) -> ReconstructionBenchmarkReport: 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"], @@ -1050,6 +1136,22 @@ def _isolated_random_state(seed: int | None) -> Iterator[None]: 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, *, @@ -1068,7 +1170,9 @@ def _split_identity( "manifest_schema_version": manifest.schema_version, "method": manifest.method, "n_samples": manifest.n_samples, - "dataset_fingerprint": manifest.dataset_fingerprint, + # 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, @@ -1281,6 +1385,8 @@ def run_reconstruction_benchmark( 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] @@ -1301,15 +1407,20 @@ def run_reconstruction_benchmark( ) built.append(instance) - fit_start = perf_counter() - instance.fit( - X_train, - X_validation if specification.uses_validation_partition else None, + # 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(X_test) + reconstruction = instance.reconstruct(method_test) reconstruction_seconds = perf_counter() - reconstruction_start reconstruction = _validate_reconstruction( @@ -1344,6 +1455,7 @@ def run_reconstruction_benchmark( 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, diff --git a/bluemath_tk/deeplearning/_base_model.py b/bluemath_tk/deeplearning/_base_model.py index c853ebd..65cafed 100644 --- a/bluemath_tk/deeplearning/_base_model.py +++ b/bluemath_tk/deeplearning/_base_model.py @@ -276,12 +276,14 @@ def _validate_fit_inputs( 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." ) - self._validate_validation_data(X, validation_data) return split = int((1 - validation_split) * len(X)) diff --git a/docs/source/benchmarking.rst b/docs/source/benchmarking.rst index fccc2bf..24ea4e1 100644 --- a/docs/source/benchmarking.rst +++ b/docs/source/benchmarking.rst @@ -99,13 +99,27 @@ an explicit message before scikit-learn is reached. accepted public workflow: ``fit`` with explicit chronological validation data and ``predict`` for reconstruction. Model architectures are never modified. -The two families are treated asymmetrically in exactly one respect, and the -report records it explicitly through ``uses_validation_partition``: +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. + 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 ----------------- @@ -125,14 +139,32 @@ A method may also be written from scratch. Anything exposing ``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. + Common metrics -------------- Metrics are computed with the accepted implementation in -:mod:`bluemath_tk.deeplearning.metrics`, so benchmark numbers agree exactly with -``reconstruction_error`` and with the per-model ``evaluate_reconstruction`` -methods. The available metrics are ``mse``, ``mae``, and ``rmse``, all reported -with ``reduction="mean"``. +: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 @@ -228,11 +260,22 @@ Results are returned as a rejected, and no timestamps are written. ``identity()`` and ``identity_digest()`` describe *what was compared*: dataset -shape, 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. +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. Timing and random state ----------------------- @@ -248,6 +291,11 @@ 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. +Two limits of that isolation are worth stating precisely. Python's standard +library ``random`` module is not isolated, so a model that draws from it is +neither seeded nor restored. On a multi-GPU host, only the current CUDA device's +generator is forked and restored. + Current limitations ------------------- @@ -264,8 +312,21 @@ Current limitations * 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, so peak memory includes one copy of the - train, validation, and test partitions. +* 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/tests/benchmarking/test_reconstruction.py b/tests/benchmarking/test_reconstruction.py index 7f9784f..9cfa87b 100644 --- a/tests/benchmarking/test_reconstruction.py +++ b/tests/benchmarking/test_reconstruction.py @@ -4,6 +4,7 @@ import copy import hashlib +import inspect import json import numpy as np @@ -21,7 +22,10 @@ pca_benchmark_method, run_reconstruction_benchmark, ) -from bluemath_tk.deeplearning.autoencoders import StandardAutoencoder # 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, @@ -510,6 +514,28 @@ def test_pca_is_not_given_the_validation_partition(): 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 # --------------------------------------------------------------------------- @@ -540,6 +566,7 @@ def predict(self, X, **kwargs): model_factory=_CapturingModel, latent_dimension=2, configuration={"architecture": "capturing"}, + shuffle_training_data=False, ) run_reconstruction_benchmark(X, split=split, methods=[specification]) @@ -589,6 +616,7 @@ class _Model: is_fitted = False def fit(self, X_train, validation_data=None, **kwargs): + self.is_fitted = True return {} def predict(self, X, **kwargs): @@ -611,6 +639,135 @@ def predict(self, X, **kwargs): 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 # --------------------------------------------------------------------------- @@ -667,9 +824,64 @@ def _destructive(values): return np.array(values, copy=True) specification, _ = _recording_spec(transform=_destructive) - run_reconstruction_benchmark(X, split=split, methods=[specification]) + 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(): @@ -761,24 +973,40 @@ def test_pca_rejects_more_components_than_available(): ) -def test_pca_rejects_more_components_than_training_samples(): +@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=10)], + 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, **kwargs): + def fit(self, X, validation_data=None, **kwargs): return {} def predict(self, X, **kwargs): @@ -1010,6 +1238,11 @@ def test_timings_are_finite_and_non_negative(): 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 @@ -1077,6 +1310,51 @@ def test_identity_digest_ignores_timing_and_metric_values(): 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) @@ -1313,22 +1591,32 @@ def test_stochastic_variational_reconstruction_is_rejected(): ) -def test_benchmark_controlled_fit_kwargs_are_rejected(): +@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, **kwargs): + def fit(self, X, validation_data=None, **kwargs): return {} def predict(self, X, **kwargs): return X - with pytest.raises(ValueError, match="remove these"): + with pytest.raises(ValueError, match="remove these fit_kwargs"): AutoencoderReconstruction( _Model(), latent_dimension=2, - fit_kwargs={"validation_split": 0.3}, + fit_kwargs=fit_kwargs, ) From d3bd8c8b33dddb9d41ecf2c11a10d7de5cea20ee Mon Sep 17 00:00:00 2001 From: Sergio Date: Sun, 16 Aug 2026 16:34:01 +0200 Subject: [PATCH 6/6] Record effective training configuration and harden report validation Three corrections from independent review of the pushed branch. The benchmark identity omitted the effective autoencoder training configuration. fit_kwargs and predict_kwargs were forwarded to the model but never recorded, so changing epochs, learning_rate, batch_size or patience left identity_digest() unchanged and two materially different experiments could claim the same identity. The effective values are now recorded automatically under configuration[fit_kwargs] and configuration[predict_kwargs], with the defaults that materially affect the experiment filled in, so a specification that omits epochs records the same identity as one that passes the default explicitly. Nothing has to be duplicated by hand. Verbosity is deliberately excluded: it controls progress reporting only and cannot change the fitted model, the reconstruction, or any metric. It still reaches the model, it simply does not enter the identity. Values in fit_kwargs and predict_kwargs must now be JSON-compatible. This is a deliberate constraint rather than an implementation limit: a 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. optimizer and criterion remain forbidden, and are now rejected when the specification is built rather than when it first runs. The model factory is still never introspected; the architecture is described by the caller's configuration. Recorded values are deep copies and each run receives freshly rebuilt containers, so mutating a nested caller mapping at any depth after building a specification changes neither later runs nor the recorded identity. Multi-GPU RNG isolation was incomplete. torch.manual_seed reseeds every visible CUDA device, but only the current device was forked, so seeding permanently altered the generators of every other device. Every visible device is now forked and restored. CUDA is still never queried when unavailable, so CPU-only environments perform no CUDA initialization, and the device selection is covered by a monkeypatched regression rather than a GPU-dependent test. Reports rebuilt from JSON now reject internally inconsistent scientific metadata instead of trusting it: partition sizes must sum to n_samples, train, validation and test must each be non-empty, split_identity must record a matching n_samples and well-formed digests, original_scalars_per_sample must equal the product of sample_shape, latent_dimension must equal latent_scalars_per_sample, and no reconstruction metric may be negative. Also validates the benchmark seed against the NumPy range up front, so an out-of-range seed fails before any method is built rather than part way through a run. Co-Authored-By: Claude Opus 5 --- bluemath_tk/benchmarking/reconstruction.py | 258 ++++++++++--- docs/source/benchmarking.rst | 53 ++- tests/benchmarking/test_reconstruction.py | 419 ++++++++++++++++++++- 3 files changed, 677 insertions(+), 53 deletions(-) diff --git a/bluemath_tk/benchmarking/reconstruction.py b/bluemath_tk/benchmarking/reconstruction.py index 841afc0..c7fa5e6 100644 --- a/bluemath_tk/benchmarking/reconstruction.py +++ b/bluemath_tk/benchmarking/reconstruction.py @@ -23,6 +23,7 @@ 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 @@ -57,6 +58,35 @@ _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]: @@ -125,6 +155,23 @@ def _validate_configuration( 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): @@ -321,6 +368,40 @@ def reconstruct(self, X: np.ndarray) -> np.ndarray: 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``. @@ -382,14 +463,7 @@ class AutoencoderReconstruction: autoencoder reconstructions is not a like-for-like measurement. """ - _FORBIDDEN_FIT_KWARGS = ( - "X", - "y", - "criterion", - "optimizer", - "validation_data", - "validation_split", - ) + _FORBIDDEN_FIT_KWARGS = _FORBIDDEN_FIT_KWARGS def __init__( self, @@ -423,24 +497,12 @@ def __init__( f"model latent width k={int(declared_k)}." ) + _reject_forbidden_fit_kwargs(fit_kwargs or {}) self._fit_kwargs = dict(fit_kwargs or {}) - forbidden = sorted( - set(self._FORBIDDEN_FIT_KWARGS).intersection(self._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}." - ) self._fit_kwargs.setdefault("verbose", 0) + _reject_stochastic_prediction(predict_kwargs or {}) self._predict_kwargs = dict(predict_kwargs or {}) - if self._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." - ) self._predict_kwargs.setdefault("verbose", 0) self._history: dict[str, list] | None = None @@ -628,18 +690,27 @@ def autoencoder_benchmark_method( Latent width of the model, cross-checked against ``model.k`` when available. configuration : mapping, optional - JSON-compatible description of the architecture and hyperparameters. - This is recorded verbatim; the factory itself is never introspected. - The reserved key ``"shuffle_training_data"`` is added automatically and - must not be supplied. + 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``. + 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``. + Extra keyword arguments for ``model.predict``, recorded the same way + under ``configuration["predict_kwargs"]``. Returns ------- @@ -651,24 +722,47 @@ def autoencoder_benchmark_method( 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) - frozen_fit_kwargs = dict(fit_kwargs or {}) - frozen_predict_kwargs = dict(predict_kwargs or {}) + + # 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") - if "shuffle_training_data" in recorded: + reserved = sorted(set(_RESERVED_CONFIGURATION_KEYS).intersection(recorded)) + if reserved: raise ValueError( - "configuration must not set the reserved key " - "'shuffle_training_data'; use the shuffle_training_data argument." + 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=frozen_fit_kwargs, - predict_kwargs=frozen_predict_kwargs, + fit_kwargs=_thaw_json(frozen_fit), + predict_kwargs=_thaw_json(frozen_predict), ) return BenchmarkMethod( @@ -744,6 +838,11 @@ def __post_init__(self) -> None: 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", @@ -779,7 +878,13 @@ def __post_init__(self) -> None: 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}.") - metrics[metric] = _validate_finite_float(f"test_metrics[{metric!r}]", value) + 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, @@ -979,6 +1084,19 @@ def __post_init__(self) -> None: 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)) @@ -989,13 +1107,27 @@ def __post_init__(self) -> None: if self.seed < 0: raise ValueError("seed must be non-negative.") - object.__setattr__( - self, - "split_identity", - _freeze_json( - _validate_configuration(self.split_identity, name="split_identity") - ), + 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 @@ -1009,12 +1141,20 @@ def __post_init__(self) -> None: 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]: @@ -1106,6 +1246,19 @@ def _optional_torch() -> Any | 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. @@ -1113,6 +1266,11 @@ def _isolated_random_state(seed: int | None) -> Iterator[None]: 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. @@ -1121,11 +1279,9 @@ def _isolated_random_state(seed: int | None) -> Iterator[None]: with ExitStack() as stack: torch = _optional_torch() if torch is not None: - devices: list[int] = [] - if torch.cuda.is_available(): # pragma: no cover - needs CUDA - current = torch.cuda.current_device() - devices = [current] - stack.enter_context(torch.random.fork_rng(devices=devices)) + stack.enter_context( + torch.random.fork_rng(devices=_forkable_torch_devices(torch)) + ) try: if seed is not None: np.random.seed(seed) @@ -1368,8 +1524,12 @@ def run_reconstruction_benchmark( 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 seed < 0: - raise ValueError("seed must be non-negative.") + 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) diff --git a/docs/source/benchmarking.rst b/docs/source/benchmarking.rst index 24ea4e1..b3d6a39 100644 --- a/docs/source/benchmarking.rst +++ b/docs/source/benchmarking.rst @@ -81,6 +81,21 @@ 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 ------------------------------------- @@ -153,6 +168,23 @@ 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 -------------- @@ -277,6 +309,17 @@ Two fingerprints are recorded, and they cover different things: 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 ----------------------- @@ -291,10 +334,14 @@ 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. -Two limits of that isolation are worth stating precisely. Python's standard +``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. On a multi-GPU host, only the current CUDA device's -generator is forked and restored. +neither seeded nor restored. Current limitations ------------------- diff --git a/tests/benchmarking/test_reconstruction.py b/tests/benchmarking/test_reconstruction.py index 9cfa87b..6ab759c 100644 --- a/tests/benchmarking/test_reconstruction.py +++ b/tests/benchmarking/test_reconstruction.py @@ -22,6 +22,12 @@ 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, @@ -35,6 +41,26 @@ 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), @@ -1514,7 +1540,7 @@ def test_every_method_starts_from_the_same_seeded_state(): assert first_created[0].random_draws == second_created[0].random_draws -@pytest.mark.parametrize("seed", [-1, 1.5, True, "3"]) +@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)) @@ -1623,3 +1649,394 @@ def predict(self, X, **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()