diff --git a/examples/aeris/__init__.py b/examples/aeris/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/examples/aeris/aeris.toml b/examples/aeris/aeris.toml
new file mode 100644
index 0000000..e3489d8
--- /dev/null
+++ b/examples/aeris/aeris.toml
@@ -0,0 +1,69 @@
+# aeris.toml  —  AERIS continuous-learning
+seed = 42
+device = "auto"
+multi_gpu = false
+verbosity = "INFO:1"
+
+[model]
+name = "aeris_init.pt"
+pretrained_path = "examples/aeris/model"
+max_ckpts = 100
+ckpts_path = "output/aeris/"
+
+[data]
+name = "aeris"
+path = "examples/aeris/data/aeris_dataset.csv"
+batch_size = 16
+
+[train]
+batch_size = 300
+num_workers = 4
+init_lr = 1e-3
+max_iter = 100
+grad_accumulation_steps = 1
+
+[continual_learning]
+update_mode = "ewc_online"
+mix_historic_data = true
+
+# JVP regularization (used when update_mode = "jvp_reg")
+jvp_rho_theta = 0.05
+jvp_rho_x = 1.0
+jvp_data_sign = 1.0
+
+# EWC (used when update_mode = "ewc_online")
+ewc_lambda = 1000.0
+ewc_ema_decay = 0.95
+
+# KFAC (used when update_mode = "kfac_online")
+kfac_lambda = 1e-2
+kfac_ema_decay = 0.95
+
+[drift_detection]
+detector_name = "KSWINDetector"
+detection_interval = 2
+aggregation = "mean"
+metric_index = 1
+reset_after_learning = true
+max_stream_updates = 84
+
+# --- ADWINDetector ---
+adwin_delta = 0.9
+adwin_minor_threshold = 0.1
+adwin_moderate_threshold = 0.6
+
+# --- KSWINDetector ---
+kswin_alpha = 0.05
+kswin_window_size = 100
+kswin_stat_size = 30
+
+# --- PageHinkleyDetector ---
+ph_min_instances = 30
+ph_delta = 0.005
+ph_threshold = 0.6
+ph_alpha = 0.9999
+
+[logging]
+backend = "wandb"
+experiment_name = "aeris-cl"
+# mlflow_tracking_uri = "http://127.0.0.1:5000"
diff --git a/examples/aeris/model.py b/examples/aeris/model.py
new file mode 100644
index 0000000..361953f
--- /dev/null
+++ b/examples/aeris/model.py
@@ -0,0 +1,349 @@
+# examples/aeris/model.py
+"""AERIS model harness for the BaseSim continuous-learning framework.
+
+This harness wraps a 8-layer neural network trained to predict enthalpy per atom from a given fuel material."""
+
+import gc
+import math
+import torch
+import numpy as np
+from typing import Tuple, Optional, List, Any
+from torch import nn, Tensor
+from torch.optim import Optimizer
+from torch.utils.data import DataLoader, ConcatDataset, TensorDataset
+
+from apeiron.model.torch_model_harness import BaseModelHarness
+from apeiron.config.configuration import Config
+
+from examples.aeris.utils import (
+    load_datasets,
+    make_loader,
+    load_pretrained_model,
+    split_into_windows,
+)
+
+# Aeris model architecture used for prediction
+class AerisFullStructure(nn.Module):
+    def __init__(self, input_dim, dropout=0.1):
+        super().__init__()
+        first_layer = min(1024, max(512, input_dim * 2))
+        self.layers = nn.Sequential(
+            nn.Linear(input_dim, first_layer), nn.ReLU(), nn.BatchNorm1d(first_layer),
+            nn.Linear(first_layer, first_layer), nn.ReLU(), nn.Dropout(dropout),
+            nn.Linear(first_layer, 512), nn.ReLU(), nn.BatchNorm1d(512),
+            nn.Linear(512, 512), nn.ReLU(), nn.Dropout(dropout),
+            nn.Linear(512, 256), nn.ReLU(), nn.BatchNorm1d(256),
+            nn.Linear(256, 256), nn.ReLU(), nn.Dropout(dropout),
+            nn.Linear(256, 128), nn.ReLU(), nn.BatchNorm1d(128),
+            nn.Linear(128, 64), nn.ReLU(), nn.Dropout(dropout),
+            nn.Linear(64, 32), nn.ReLU(), nn.Linear(32, 1)
+        )
+
+    def forward(self, x):
+        return self.layers(x)
+
+    def train(self, mode: bool = True):
+        """Set training mode, but always keep BatchNorm layers in eval mode.
+
+        The ``jvp_reg`` continual-learning updater runs this model through
+        ``torch.func.jvp``/``grad`` via ``functional_call``. In training mode
+        ``BatchNorm1d`` performs an in-place ``num_batches_tracked.add_(1)`` on a
+        captured buffer, which functorch transforms forbid. Keeping BatchNorm in
+        eval mode (frozen running stats) both avoids that crash and prevents the
+        small, drifted CL batches from corrupting the normalization statistics.
+        Dropout still follows ``mode`` normally.
+        """
+        super().train(mode)
+        for m in self.modules():
+            if isinstance(m, nn.BatchNorm1d):
+                m.eval()
+        return self
+
+
+# Fraction of each time window reserved for validation
+_VAL_FRACTION: float = 0.2
+
+
+class AERIS(BaseModelHarness):
+    """
+    Continuous-learning harness for the AERIS prediction model.
+    """
+
+    def __init__(self, cfg: Config):
+        # ----- build model ---------------------------------------------------
+        ckpt = load_pretrained_model(
+            cfg.model.pretrained_path, cfg.model.name, device=cfg.device
+        )
+
+        # Expect full harness checkpoint format
+        if not isinstance(ckpt, dict) or "model_state_dict" not in ckpt:
+            raise ValueError(
+                f"Checkpoint '{cfg.model.name}' must be a full harness checkpoint "
+                "with 'model_state_dict', 'feature_names', 'scaler', and 'input_dim'. "
+                "Legacy weight-only checkpoints are no longer supported."
+            )
+
+        state_dict = ckpt["model_state_dict"]
+        feature_names = ckpt["feature_names"]
+        scaler = ckpt["scaler"]
+        input_dim = int(ckpt["input_dim"])
+
+        model = AerisFullStructure(input_dim=input_dim)
+        model.load_state_dict(state_dict)
+        model.to(cfg.device)
+        model.eval()
+
+        super().__init__(cfg=cfg, model=model)
+
+        self._feature_names = feature_names
+        self._scaler = scaler
+        self._input_dim = input_dim
+
+        # ----- data loaders  -------------------------------------
+        X, y = load_datasets(cfg.data.path, cfg.data.name, feature_names, input_dim)
+        # X shape: (n_samples, 245) y shape: (n_samples, 1)
+
+        # scale (must match training)
+        X_scaled = scaler.transform(X).astype(np.float32)
+        X_tensor = torch.tensor(X_scaled, dtype=torch.float32)
+        y_tensor = torch.tensor(y, dtype=torch.float32).view(-1, 1)
+
+        self.windows = split_into_windows(X_tensor, y_tensor)
+        print(f"Prepared {len(self.windows)} time windows for streaming. Each window has ~{self.windows[0][0].shape[0]} samples.")
+
+        # ----- optional base/initial training data ---------------------------
+        # The pre-drift training split the model was originally fit on. When
+        # cfg.data.base_train_path is set it is blended into full retrains via
+        # get_base_train_dataloaders() so the from-scratch model does not forget
+        # the base distribution. Featurized + scaled exactly like the stream.
+        self._base_train_ds: Optional[TensorDataset] = None
+        self._base_val_ds: Optional[TensorDataset] = None
+        base_path = getattr(cfg.data, "base_train_path", "")
+        if base_path:
+            Xb, yb = load_datasets(base_path, cfg.data.name, feature_names, input_dim)
+            Xb_scaled = scaler.transform(Xb).astype(np.float32)
+            Xb_t = torch.tensor(Xb_scaled, dtype=torch.float32)
+            yb_t = torch.tensor(yb, dtype=torch.float32).view(-1, 1)
+            # Shuffle once (seeded) so the val slice isn't a biased tail.
+            gen = torch.Generator().manual_seed(cfg.seed)
+            perm = torch.randperm(Xb_t.shape[0], generator=gen)
+            Xb_t, yb_t = Xb_t[perm], yb_t[perm]
+            n = Xb_t.shape[0]
+            n_val = max(1, int(n * _VAL_FRACTION))
+            n_train = n - n_val
+            self._base_train_ds = TensorDataset(Xb_t[:n_train], yb_t[:n_train])
+            self._base_val_ds = TensorDataset(Xb_t[n_train:], yb_t[n_train:])
+            print(
+                f"Loaded {n} base training samples from {base_path} "
+                "(blended into full retrains)."
+            )
+
+        # ----- eval metrics (prediction) -------------------------------------
+        self._y_var_ref = self._reference_variance(y_tensor)
+        self.eval_metrics = {
+            "mse": self.mse_metric(),
+            "mae": self.get_criterion(),
+            "r2": self.r2_metric(),
+            "nrmse": self.nrmse_metric(),
+        }
+        self.higher_is_better = {
+            "mse": False, "mae": False, "r2": True, "nrmse": False,
+        }
+
+        # ----- streaming state -----------------------------------------------
+        self.window_idx: int = 0
+        self.history_windows: List[Tuple[Tensor, Tensor]] = []
+
+        self._cur_train_loader: Optional[DataLoader] = None
+        self._cur_val_loader: Optional[DataLoader] = None
+        self._cur_stream_loader: Optional[DataLoader] = None
+
+    def _reference_variance(self, y_stream: Tensor) -> float:
+        """Fixed denominator for R^2 / NRMSE metrics."""
+        y_ref = (
+            self._base_train_ds.tensors[1]
+            if self._base_train_ds is not None
+            else y_stream
+        )
+        var = float(y_ref.float().var(unbiased=False).item())
+        if var <= 0.0:
+            raise ValueError("Reference target variance is 0; R^2/NRMSE undefined.")
+        return var
+
+    def r2_metric(self):
+        """1 - MSE/var_ref against a constant denominator.
+
+        Being affine in MSE means the sample-weighted mean that
+        BaseModelHarness.eval() computes is exactly the pooled R^2 -- which
+        would NOT hold if each batch normalized by its own variance.
+        """
+        var_ref = self._y_var_ref
+
+        def _r2(y_hat: Tensor, y: Tensor) -> Tensor:
+            return 1.0 - torch.mean((y_hat - y) ** 2) / var_ref
+
+        return _r2
+
+    def nrmse_metric(self):
+        """RMSE / std_ref, same reference as r2_metric."""
+        std_ref = math.sqrt(self._y_var_ref)
+
+        def _nrmse(y_hat: Tensor, y: Tensor) -> Tensor:
+            return torch.sqrt(torch.mean((y_hat - y) ** 2)) / std_ref
+
+        return _nrmse
+
+    def get_optmizer(self) -> Optimizer:  # noqa: D102  (spelling kept for ABC)
+        weight_decay = 1e-7
+        return torch.optim.AdamW(self.model.parameters(), lr=self.cfg.train.init_lr, weight_decay=weight_decay)
+
+    def mse_metric(self):  # noqa: D102
+        return nn.MSELoss()
+
+    def get_criterion(self):
+        return nn.L1Loss()
+
+    def get_stream_dataloader(self):
+        assert self._cur_stream_loader is not None
+        return self._cur_stream_loader
+
+    def get_train_dataloaders(self) -> Tuple[DataLoader, DataLoader]:  # noqa: D102
+        assert self._cur_train_loader is not None and self._cur_val_loader is not None
+        return self._cur_train_loader, self._cur_val_loader
+
+    def get_hist_dataloaders(
+        self,
+    ) -> Tuple[Optional[DataLoader], Optional[DataLoader]]:
+        """Return loaders over all previously-seen time windows.
+
+        Returns ``(None, None)`` until at least two windows have been served.
+        """
+        if self.window_idx <= 1:
+            return None, None
+
+        # Concatenate all history windows
+        hist_train_views: List[TensorDataset] = []
+        hist_val_views: List[TensorDataset] = []
+
+        for X_w, y_w in self.history_windows:
+            n = X_w.shape[0]
+            n_val = max(1, int(n * _VAL_FRACTION))
+            n_train = n - n_val
+            hist_train_views.append(TensorDataset(X_w[:n_train], y_w[:n_train]))
+            hist_val_views.append(TensorDataset(X_w[n_train:], y_w[n_train:]))
+
+        ds_hist_train: ConcatDataset[Any] = ConcatDataset(hist_train_views)
+        ds_hist_val: ConcatDataset[Any] = ConcatDataset(hist_val_views)
+
+        bs = self.cfg.train.batch_size
+        nw = self.cfg.train.num_workers
+        pin = torch.cuda.is_available()
+        return (
+            make_loader(
+                ds_hist_train, bs, shuffle=True, num_workers=nw, pin_memory=pin
+            ),
+            make_loader(ds_hist_val, bs, shuffle=False, num_workers=nw, pin_memory=pin),
+        )
+
+    def get_base_train_dataloaders(
+        self,
+    ) -> Tuple[Optional[DataLoader], Optional[DataLoader]]:
+        """Return (train, val) loaders over the initial training split.
+
+        Returns ``(None, None)`` unless ``cfg.data.base_train_path`` was set.
+        Loaders are built on demand (only during a full retrain) so no worker
+        processes are held open for the rest of the run.
+        """
+        if self._base_train_ds is None:
+            return None, None
+
+        bs = self.cfg.train.batch_size
+        nw = self.cfg.train.num_workers
+        pin = torch.cuda.is_available()
+        train_loader = make_loader(
+            self._base_train_ds, bs, shuffle=True, num_workers=nw, pin_memory=pin
+        )
+        val_loader = (
+            make_loader(
+                self._base_val_ds, bs, shuffle=False, num_workers=nw, pin_memory=pin
+            )
+            if self._base_val_ds is not None
+            else None
+        )
+        return train_loader, val_loader
+
+    def update_data_stream(self) -> None:
+        """Advance to the next chronological time window.
+
+        The current window is added to the history, and new train/val loaders
+        are built from the upcoming window.
+        """
+        self._dispose_current_loaders()
+
+        if self.window_idx >= len(self.windows):
+            print(
+                f"Warning: All {len(self.windows)} time windows exhausted; "
+                "wrapping around to the first window."
+            )
+            self.window_idx = 0
+
+        X_w, y_w = self.windows[self.window_idx]
+
+        # Archive previous window in history (skip the very first call)
+        if self.window_idx > 0:
+            prev_X, prev_y = self.windows[self.window_idx - 1]
+            # Only add if not already stored (idempotency guard)
+            if len(self.history_windows) < self.window_idx:
+                self.history_windows.append((prev_X, prev_y))
+        # Train / val split (last _VAL_FRACTION chronologically)
+        n = X_w.shape[0]
+        n_val = max(1, int(n * _VAL_FRACTION))
+        n_train = n - n_val
+
+        ds_train = TensorDataset(X_w[:n_train], y_w[:n_train])
+        ds_val = TensorDataset(X_w[n_train:], y_w[n_train:])
+
+        bs = self.cfg.train.batch_size
+        nw = self.cfg.train.num_workers
+        pin = torch.cuda.is_available()
+
+        self._cur_train_loader = make_loader(
+            ds_train, bs, shuffle=True, num_workers=nw, pin_memory=pin
+        )
+        self._cur_val_loader = make_loader(
+            ds_val, bs, shuffle=False, num_workers=nw, pin_memory=pin
+        )
+
+        bs = self.cfg.data.batch_size
+        self._cur_stream_loader = make_loader(
+            ds_train, bs, shuffle=True, num_workers=nw, pin_memory=pin
+        )
+        self.window_idx += 1
+
+    def build_checkpoint_payload(self) -> dict[str, Any]:
+        """Save full harness checkpoint matching load_pretrained_model format.
+
+        Overrides BaseModelHarness to include feature names, scaler, and input_dim
+        so that saved checkpoints can be loaded directly without reference checkpoints.
+        """
+        return {
+            "model_state_dict": self.model.state_dict(),
+            "feature_names": self._feature_names,
+            "scaler": self._scaler,
+            "input_dim": self._input_dim,
+        }
+
+    # --------------------------------------------------------------------- #
+    # Helpers
+    # --------------------------------------------------------------------- #
+    def _dispose_current_loaders(self) -> None:
+        if self._cur_train_loader is not None:
+            del self._cur_train_loader
+            self._cur_train_loader = None
+        if self._cur_val_loader is not None:
+            del self._cur_val_loader
+            self._cur_val_loader = None
+        if self._cur_stream_loader is not None:
+            del self._cur_stream_loader
+            self._cur_stream_loader = None
+        gc.collect()
diff --git a/examples/aeris/utils.py b/examples/aeris/utils.py
new file mode 100644
index 0000000..46c3221
--- /dev/null
+++ b/examples/aeris/utils.py
@@ -0,0 +1,234 @@
+# examples/aeris/utils.py
+"""Utility functions for the AERIS continuous-learning example.
+
+Expected directory layout (pointed to by ``cfg.data.path``)::
+
+    <data_dir>/
+        dataset.csv            # data that will be parsed by the SIM framework
+        aeris_model.pt        # AERIS pre-trained model
+
+Featurization uses the *fast path* and must stay byte-for-byte consistent with
+``examples/aeris/scripts/make_drift_split.py:build_features`` -- the same code
+that produced the training features. 227 of the 233 features are pre-computed
+columns pulled directly from the CSV; the remaining 6 lattice params are parsed
+from the ``structure`` string. No matminer recompute (that path could disagree
+with the pre-computed columns the model was trained on).
+"""
+
+import os
+import glob
+import re
+from typing import Dict, List, Tuple, Any, Optional
+
+import numpy as np
+import pandas as pd
+import torch
+from torch import Tensor
+from torch.utils.data import DataLoader, Dataset
+
+
+def load_pretrained_model(
+    data_path: str, model_name: str, device: str = "cpu"
+) -> dict[str, Any]:
+    """Load the pretrained AERIS model.
+
+    Parameters
+    ----------
+    data_path:
+        Directory containing the model.
+    model_name:
+        The name of the pretrained model.
+    device:
+        Device to map the scalers to.
+
+    Returns
+    -------
+    model_info = {
+        'model_state_dict': model.state_dict(),
+        'input_dim': input_dim,
+        'feature_names': feature_names,
+        'scaler': scaler,
+        'metrics': {'mae': mae, 'rmse': rmse, 'r2': r2},
+        'history': history
+    }
+    """
+    ckpt = None
+    if os.path.exists(data_path):
+        ckpt = torch.load(
+            os.path.join(data_path, model_name), map_location=device, weights_only=False
+        )
+    if ckpt is None:
+        raise FileNotFoundError("No model found at path: " + data_path)
+    return ckpt
+
+
+# -----------------------------
+# Fast-path featurization
+# (mirrors scripts/make_drift_split.py so the harness featurizes inputs exactly
+#  the way the model was trained)
+# -----------------------------
+LATTICE_KEYS = [
+    "lattice_a", "lattice_b", "lattice_c",
+    "lattice_alpha", "lattice_beta", "lattice_gamma",
+]
+
+
+def _parse_lattice(struct_str: Any) -> Dict[str, float]:
+    """Extract the 6 lattice params from a pymatgen structure string."""
+    r = {k: 0.0 for k in LATTICE_KEYS}
+    s = str(struct_str)
+    abc = re.search(r"abc\s*:\s*([\d.]+)\s+([\d.]+)\s+([\d.]+)", s)
+    ang = re.search(r"angles\s*:\s*([\d.]+)\s+([\d.]+)\s+([\d.]+)", s)
+    if abc:
+        r["lattice_a"], r["lattice_b"], r["lattice_c"] = map(float, abc.groups())
+    if ang:
+        r["lattice_alpha"], r["lattice_beta"], r["lattice_gamma"] = map(float, ang.groups())
+    return r
+
+
+def _build_X_fast(df: pd.DataFrame, feature_names: List[str]) -> np.ndarray:
+    """Assemble the ``(N, len(feature_names))`` matrix via the fast path.
+
+    227 features are pulled directly from pre-computed CSV columns; the 6 lattice
+    params are parsed from the ``structure`` string. Any feature not found is
+    left at 0. Identical to ``make_drift_split.build_features``.
+    """
+    n = len(df)
+    X = np.zeros((n, len(feature_names)), dtype=np.float32)
+    col_idx = {f: j for j, f in enumerate(feature_names)}
+
+    present = [f for f in feature_names if f in df.columns]
+    sub = df[present].apply(pd.to_numeric, errors="coerce").to_numpy(np.float32)
+    for k, f in enumerate(present):
+        X[:, col_idx[f]] = sub[:, k]
+
+    if "structure" in df.columns:
+        lat = np.array(
+            [list(_parse_lattice(s).values()) for s in df["structure"].tolist()],
+            dtype=np.float32,
+        )
+        for k, f in enumerate(LATTICE_KEYS):
+            if f in col_idx:
+                X[:, col_idx[f]] = lat[:, k]
+
+    return np.nan_to_num(X, nan=0.0, posinf=1e6, neginf=-1e6)
+
+
+def load_datasets(
+    data_path: str, dataset_name: str, feature_names: List[str], input_dim: int
+) -> Tuple[np.ndarray, Optional[np.ndarray]]:
+    """Load the dataset used by the model.
+
+    Features are assembled in the exact ``feature_names`` order via the fast
+    path (pre-computed columns + parsed lattice params), matching how the model
+    was trained. Rows with a missing target are dropped.
+
+    Returns
+    -------
+    X: numpy.ndarray of shape (n_samples, n_features) dtype float32 (unscaled)
+    y: numpy.ndarray of shape (n_samples, 1) dtype float32
+
+    Note: scaling is intentionally NOT applied here. The caller (model harness)
+    applies the saved scaler from the checkpoint via scaler.transform().
+    """
+    dataset_pattern = os.path.join(data_path)
+    dataset_files: List[str] = glob.glob(dataset_pattern)
+    if not dataset_files:
+        raise FileNotFoundError(f"No dataset files matched pattern: {dataset_pattern}")
+
+    dfs = [pd.read_csv(fp, low_memory=False) for fp in dataset_files]
+    dataset: pd.DataFrame = pd.concat(dfs, ignore_index=True)
+
+    target_col = "formation_energy_per_atom"
+    if target_col not in dataset.columns:
+        raise KeyError(f"Missing target column '{target_col}'")
+    dataset = dataset[dataset[target_col].notna()].reset_index(drop=True)
+
+    X = _build_X_fast(dataset, feature_names)
+    y = dataset[target_col].to_numpy(np.float32).reshape(-1, 1)
+
+    if X.shape[1] != input_dim:
+        raise ValueError(
+            f"Checkpoint input_dim={input_dim} but built X has {X.shape[1]} features."
+        )
+
+    return X, y
+
+
+# Default number of samples per time window.  Can be overridden by the caller.
+DEFAULT_WINDOW_SIZE: int = 500
+
+def split_into_windows(
+    X: Tensor,
+    y: Tensor,
+    window_size: int = DEFAULT_WINDOW_SIZE,
+) -> List[Tuple[Tensor, Tensor]]:
+    """Split chronologically-ordered tensors into non-overlapping windows.
+
+    Any leftover samples that don't fill a complete window are appended as
+    a final (smaller) window so no data is discarded.
+
+    Parameters
+    ----------
+    X:
+        Input features ``[N, D]``.
+    y:
+        Targets ``[N, T]``.
+    window_size:
+        Number of samples per window.
+
+    Returns
+    -------
+    List of ``(X_chunk, y_chunk)`` tuples.
+    """
+    n = X.shape[0]
+    windows: List[Tuple[Tensor, Tensor]] = []
+    for start in range(0, n, window_size):
+        end = min(start + window_size, n)
+        windows.append((X[start:end], y[start:end]))
+    return windows
+
+
+def make_loader(
+    ds: Dataset,
+    batch_size: int,
+    shuffle: bool,
+    num_workers: int = 4,
+    pin_memory: bool = True,
+    persistent_workers: bool = True,
+    prefetch_factor: int = 2,
+) -> DataLoader:
+    """Build a ``DataLoader`` from a ``Dataset``.
+
+    Parameters
+    ----------
+    ds:
+        The base dataset.
+    batch_size:
+        Batch size.
+    shuffle:
+        Whether to shuffle.
+    num_workers:
+        Number of data-loading workers.
+    pin_memory:
+        Pin CUDA memory for faster transfers.
+    persistent_workers:
+        Keep worker processes alive between iterations.
+    prefetch_factor:
+        Samples to prefetch per worker.
+
+    Returns
+    -------
+    DataLoader
+    """
+    kwargs: dict = dict(batch_size=batch_size, shuffle=shuffle, drop_last=False)
+    if num_workers > 0:
+        kwargs.update(
+            dict(
+                num_workers=num_workers,
+                pin_memory=pin_memory,
+                persistent_workers=persistent_workers,
+                prefetch_factor=prefetch_factor,
+            )
+        )
+    return DataLoader(ds, **kwargs)  # type: ignore[arg-type]
diff --git a/examples/mnist/mnist.toml b/examples/mnist/mnist.toml
index 12449b1..d7b288d 100644
--- a/examples/mnist/mnist.toml
+++ b/examples/mnist/mnist.toml
@@ -26,7 +26,7 @@ grad_accumulation_steps = 2
 
 [continual_learning]
 update_mode = "kfac_online"
-mix_historic_data = false
+mix_historic_data = true
 
 # jvp_reg
 jvp_rho_theta = 0.05
diff --git a/examples/utils.py b/examples/utils.py
index 0cde7b7..915d970 100644
--- a/examples/utils.py
+++ b/examples/utils.py
@@ -15,6 +15,10 @@ def get_example(cfg: Config) -> BaseModelHarness:
         from examples.imagenet.model import IMAGENET_VISION
 
         return IMAGENET_VISION(cfg=cfg)
+    elif cfg.data.name == "aeris":
+        from examples.aeris.model import AERIS
+
+        return AERIS(cfg=cfg)
     else:
         raise NotImplementedError(
             f"Example for dataset {cfg.data.name} is not implemented."
diff --git a/src/apeiron/model/torch_model_harness.py b/src/apeiron/model/torch_model_harness.py
index 07c2c29..1768ee0 100644
--- a/src/apeiron/model/torch_model_harness.py
+++ b/src/apeiron/model/torch_model_harness.py
@@ -222,13 +222,27 @@ def task_diagonals(self) -> List[List[float]]:
     def ckpts_enabled(self) -> bool:
         return self.cfg.model.max_ckpts > 0 and bool(self.cfg.model.ckpts_path)
 
+    def build_checkpoint_payload(self) -> Any:
+        """Build the checkpoint object to save.
+
+        Subclasses can override this to include additional metadata beyond weights
+        (e.g., preprocessing scalers, feature names, architecture parameters)
+        so that saved checkpoints match the format expected by the loader.
+
+        Returns
+        -------
+        By default, returns ``model.state_dict()`` (weights only).
+        """
+        return self.model.state_dict()
+
     def save_ckpt(self, event: int) -> str:
         """Persist model state, evict oldest when over budget."""
         d = Path(self.cfg.model.ckpts_path)
         d.mkdir(parents=True, exist_ok=True)
 
         fname = f"drift_adaptation_{event}.pt"
-        torch.save(self.model.state_dict(), d / fname)
+        payload = self.build_checkpoint_payload()
+        torch.save(payload, d / fname)
         (d / "latest").write_text(fname)
 
         # Guillotine the oldest survivors