diff --git a/CHANGELOG.md b/CHANGELOG.md
index f999b77..031b36f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,21 @@ All notable changes to ChimeraBoost are documented here.
 The format follows [Keep a Changelog](https://keepachangelog.com/).
 
 ## [Unreleased]
+### Added
+- **`store_training_data` and `refresh(X, y)`** on the regressor and the
+  binary classifier (#131, first slice). With `store_training_data=True` a
+  model keeps its training rows in compact form (bins for plain numeric
+  columns, raw values only for the columns that feed cross features,
+  category codes), and `refresh` folds in new rows by replaying every
+  tree's splits on the old and new rows together and refitting only the
+  leaf values. It is the same replay the default full-data refit uses,
+  with the tree count, learning rate and splits pinned; refreshing with no
+  new rows reproduces the model bit for bit. Bagged models, multiclass,
+  `loss="Quantile"` and `random_effects=True` are not supported yet. The
+  default fit's results are unchanged. On a 3.3M-row regression, folding in
+  a day of new rows took 70 s against 232 s for a full refit, at the same
+  accuracy.
+
 ### Changed
 - **The quantile head retrains its winner on all rows by default.** Called
   without an `eval_set`, `ChimeraBoostQuantileRegressor` held back
@@ -49,6 +64,11 @@ The format follows [Keep a Changelog](https://keepachangelog.com/).
   untouched: the exact-output snapshot moves only one quantile
   configuration (183 of 186 pins identical). Record:
   `benchmarks/QUANTILE_PLAN.md`, Phase 2 (Q8, Q9).
+- **Fits with linear leaves are faster, with identical results.** The
+  per-leaf linear fit used to sum each leaf on one thread, and one leaf
+  often holds 40-60% of the rows. Large leaves now spread their sums over
+  several threads, each sum still adding the rows in the same order. A
+  500k-row fit went from 15.9 to 13.4 s, and a refresh from 5.0 to 3.6 s.
 - **The quantile benchmark's NGBoost opponent now uses the RoNGBa
   settings** (Ren, Sun and Wu 2019; #163): trees of up to 31 leaves, a
   learning rate of 0.04 and at most 500 rounds, in place of NGBoost's stock
diff --git a/benchmarks/CAMPAIGN_PLAN.md b/benchmarks/CAMPAIGN_PLAN.md
index 928c49a..d67c77d 100644
--- a/benchmarks/CAMPAIGN_PLAN.md
+++ b/benchmarks/CAMPAIGN_PLAN.md
@@ -1277,6 +1277,11 @@ PR #170 is PARKED OPEN (not merged, not closed); issue #131 stays open.
 The loop moves on to the other issues; if he closes #170, the kernel
 speedup (435e866, bit-identical, 1.19x on linear-leaf fits) is worth
 salvaging as its own PR, since it helps every default fit.
+2026-09-26, the maintainer: "let's move forward with 170, but clean up
+the merge conflicts then i'll merge the PR". Un-parked: main merged into
+the branch; the conflicts were CHANGELOG (entries on both sides, all
+kept), this file and `REFRESH_PLAN.md` (log lines on both sides, all
+kept); no code overlapped. Awaiting his merge.
 
 #### I065 2026-09-24 issue #81 (research cascade: dead self-test anchor, stale `ideas.py` flags; BENCH tooling + test, pre-registered)
 why now: the focus rule's second issue. PR #168 merged (3e02ab1), #84
diff --git a/benchmarks/REFRESH_PLAN.md b/benchmarks/REFRESH_PLAN.md
index 279d4fd..fdf227a 100644
--- a/benchmarks/REFRESH_PLAN.md
+++ b/benchmarks/REFRESH_PLAN.md
@@ -131,3 +131,8 @@ Still open, for their slices:
   implemented"). Slice 1 is PARKED; slices 2-7 are not started. If #170 is
   closed, salvage the bit-identical replay-kernel speedup (435e866) as its
   own PR: it speeds up every linear-leaf fit.
+- 2026-09-26: un-parked. The maintainer: "let's move forward with 170,
+  but clean up the merge conflicts then i'll merge the PR". Main merged
+  into the branch (conflicts only in CHANGELOG and two plan files; no
+  code overlapped), so the salvage note above no longer applies. Slices
+  2-7 are not started.
diff --git a/chimeraboost/booster.py b/chimeraboost/booster.py
index f73f8a2..afa7ac0 100644
--- a/chimeraboost/booster.py
+++ b/chimeraboost/booster.py
@@ -824,6 +824,31 @@ def __init__(self, loss="RMSE", loss_kwargs=None, **kw):
         self.loss_name = loss
         self.loss_kwargs = loss_kwargs or {}
 
+    def replay_kwargs(self):
+        """Constructor kwargs that rebuild an equivalent booster.
+
+        Every ``_BaseBooster.__init__`` parameter (each stored under its own
+        name, so read off ``inspect.signature`` rather than a hand-kept
+        list) plus ``loss`` (from ``loss_name``) and ``loss_kwargs``; lists
+        and dicts are copied.
+        """
+        import inspect
+
+        kw = {}
+        params = inspect.signature(_BaseBooster.__init__).parameters
+        for name in params:
+            if name == "self":
+                continue
+            v = getattr(self, name)
+            if isinstance(v, list):
+                v = list(v)
+            elif isinstance(v, dict):
+                v = dict(v)
+            kw[name] = v
+        kw["loss"] = self.loss_name
+        kw["loss_kwargs"] = dict(self.loss_kwargs)
+        return kw
+
     def _fit_impl(self, X, y, cat_features=None, eval_set=None,
                   sample_weight=None, callbacks=None, prep_cache=None):
         """Fit the additive model.
diff --git a/chimeraboost/sklearn_api.py b/chimeraboost/sklearn_api.py
index c348789..02e74f8 100644
--- a/chimeraboost/sklearn_api.py
+++ b/chimeraboost/sklearn_api.py
@@ -12,6 +12,7 @@
 from .random_effects import (codes_for_labels, estimate_ratio_reml,
                              solve_intercepts)
 from .target_encoding import factorize
+from .training_rows import TrainingRows
 from sklearn.base import BaseEstimator, RegressorMixin, ClassifierMixin
 
 
@@ -60,7 +61,7 @@ def loss(T):
                            "cat_features", "cross_features",
                            "cross_top_columns",
                            "selection_rounds", "refit_full", "refit_members",
-                           "quality", "random_effects"})
+                           "quality", "random_effects", "store_training_data"})
 
 # --- quality: named operating points on the strength/slowdown Pareto --------
 # Evidence: benchmarks/SELECT_PLAN.md. Every recipe only pins parameters that
@@ -235,7 +236,7 @@ def _check_depth(p, name):
 
 
 def _check_flag_params(estimator, p):
-    """The tri-state string/bool flags: refit_full, cross_features, quality."""
+    """The string/bool flags: refit_full, cross_features, quality, store."""
     v = p.get("refit_full")
     if v is not None and v != "replay" and not isinstance(v, (bool, np.bool_)):
         raise ValueError(
@@ -257,6 +258,12 @@ def _check_flag_params(estimator, p):
                 + ", ".join(f"{k} ({n})" for k, n in QUALITY_NAMES.items())
                 + f"; got {v!r}.")
 
+    if "store_training_data" in p:
+        v = p["store_training_data"]
+        if not isinstance(v, (bool, np.bool_)):
+            raise ValueError(
+                f"store_training_data must be True or False; got {v!r}.")
+
 
 def _check_loss_family(p):
     # Regressor-only loss / alpha (the classifier picks its loss automatically).
@@ -1870,6 +1877,174 @@ def _add_callback(callbacks, extra):
     return base + [extra]
 
 
+def _raise_store_unsupported(what):
+    """The slice-1 coverage error for ``store_training_data=True`` fits."""
+    raise NotImplementedError(
+        f"store_training_data=True does not support {what} yet; refresh() "
+        "currently covers single-model regression and binary classification.")
+
+
+def _check_store_fit(est):
+    """Raise unless this ``store_training_data=True`` fit is slice-1 shaped.
+
+    Called inside ``_quality_applied`` so quality 4/5 already resolved to
+    ``n_ensembles``; the multiclass arm runs in the classifier's
+    ``_fit_single`` once the classes are known.
+    """
+    if not est.store_training_data:
+        return
+    if est.n_ensembles and est.n_ensembles > 1:
+        _raise_store_unsupported("n_ensembles > 1 (including quality=4 and 5)")
+    if getattr(est, "loss", None) == "Quantile":
+        _raise_store_unsupported("loss='Quantile'")
+    if getattr(est, "random_effects", False):
+        _raise_store_unsupported("random_effects=True")
+
+
+def _check_store_multiclass(est):
+    """The multiclass arm of the slice-1 gate (classes are known by now)."""
+    if est.store_training_data and est.n_classes_ > 2:
+        _raise_store_unsupported("multiclass classification")
+
+
+def _capture_training_store(est, cap_full, cap_train, refit, *,
+                              classification):
+    """Set ``n_samples_trained_`` and ``_training_data_`` at the end of fit.
+
+    The captured rows are ``cap_full`` (every row) when the full-data refit
+    replaced the model, else ``cap_train`` (the post-split training rows --
+    X as given with an explicit eval_set or ``early_stopping=False``), each
+    an ``(X, y, sample_weight)`` triple in the booster's row order. The
+    captured ``y`` is exactly as the booster saw it (0/1 floats for binary),
+    the weights raw and un-normalized.
+    """
+    cap_X, cap_y, cap_w = cap_full if refit else cap_train
+    est.n_samples_trained_ = int(len(cap_y))
+    if not est.store_training_data:
+        est._training_data_ = None
+        return
+    if classification:
+        # Binary only (multiclass raised before any fitting work).
+        cap_y = (cap_y == est.classes_[1]).astype(np.float64)
+    est._training_data_ = TrainingRows.capture(
+        est.model_.prep_, cap_X, cap_y, cap_w)
+
+
+def _refresh_check_y(y, n_new, classification):
+    """Validate refresh targets; the empty case skips the target-type probe."""
+    if n_new:
+        return _check_y_target(y, n_new, classification)
+    y_arr = np.asarray(y)
+    if y_arr.shape[0] != 0:
+        raise ValueError(
+            "X and y have inconsistent lengths: X has 0 samples, "
+            f"y has {y_arr.shape[0]}.")
+    if y_arr.ndim == 2:
+        if y_arr.shape[1] == 1:
+            return y_arr.ravel()
+        raise ValueError(
+            "Multi-output y is not supported; pass a 1D y of shape "
+            "(n_samples,).")
+    return y_arr
+
+
+def _refresh_check_w(sample_weight, n_new):
+    """Validate refresh weights; the empty case only checks the shape."""
+    if sample_weight is None:
+        return None
+    if n_new:
+        _check_sample_weight_arr(sample_weight, n_new)
+        return np.asarray(sample_weight, dtype=np.float64)
+    sw = np.asarray(sample_weight, dtype=np.float64)
+    if sw.ndim != 1 or sw.shape[0] != 0:
+        raise ValueError(
+            f"sample_weight must be 1D of length 0; got shape {sw.shape}.")
+    return sw
+
+
+def _refresh_check_store_weights(rows0, sample_weight):
+    """The weighted/unweighted refresh presence checks."""
+    had_weights = rows0.sample_weight is not None
+    if had_weights and sample_weight is None:
+        raise ValueError(
+            "This model was fit with sample_weight, so refresh() needs "
+            "sample_weight for the new rows too, in the same units.")
+    if not had_weights and sample_weight is not None:
+        raise ValueError(
+            "This model was fit without sample_weight, so refresh() cannot "
+            "take sample_weight for the new rows; refit with sample_weight "
+            "to use weights.")
+
+
+def _refresh_map_labels(est, y_arr):
+    """Map refresh labels to the fitted 0/1 encoding, rejecting unseen ones."""
+    yv = np.asarray(y_arr)
+    if yv.size == 0:
+        return np.empty(0, dtype=np.float64)
+    try:
+        unseen = np.setdiff1d(np.unique(yv),
+                              np.unique(np.asarray(est.classes_)))
+    except TypeError:   # mixed un-orderable label types
+        seen = set(np.asarray(est.classes_).tolist())
+        unseen = np.array([v for v in dict.fromkeys(yv.tolist())
+                          if v not in seen], dtype=object)
+    if unseen.size:
+        raise ValueError(
+            f"y contains class label(s) {unseen.tolist()} not seen at fit; "
+            "refresh() cannot add classes -- refit instead.")
+    return (yv == est.classes_[1]).astype(np.float64)
+
+
+def _refresh_replay(est, rows, X_all):
+    """Replay the fitted booster's configuration on the combined rows.
+
+    Everything is pinned at its fitted value -- the size-adaptive autos are
+    NOT re-resolved the way ``_refit_on_full`` does -- and the donor's
+    histories are carried over so ``validation_history_`` still reports the
+    curve that chose the budget.
+    """
+    donor = est.model_
+    kw = donor.replay_kwargs()
+    kw["n_estimators"] = len(donor.trees_)
+    kw["learning_rate"] = float(donor.lr_)
+    kw["early_stopping_rounds"] = None
+    kw["replay_donor"] = (donor.trees_, donor.prep_)
+    b = GradientBoosting(**kw)
+    b.fit(X_all, rows.y, cat_features=list(donor.prep_.cat_features_),
+          sample_weight=rows.sample_weight)
+    b.train_history_ = donor.train_history_
+    b.valid_history_ = donor.valid_history_
+    est.model_ = b
+    est._training_data_ = rows
+    est.n_samples_trained_ = rows.n_rows
+    est._shap_importances_cache_ = None
+    if hasattr(est, "expected_value_"):
+        del est.expected_value_
+    return est
+
+
+def _refresh_single(est, X, y, sample_weight, *, classification):
+    """Shared body of the estimators' ``refresh`` (issue #131 slice 1)."""
+    Xv = _check_predict_input(est, X)
+    rows0 = getattr(est, "_training_data_", None)
+    if rows0 is None:
+        raise ValueError(
+            "refresh() needs the rows this model was trained on, but it was "
+            "fit with store_training_data=False. Refit with "
+            "store_training_data=True to enable refresh().")
+    n_new = len(Xv) if Xv is not None else len(X)
+    y_arr = _refresh_check_y(y, n_new, classification)
+    w_new = _refresh_check_w(sample_weight, n_new)
+    _refresh_check_store_weights(rows0, sample_weight)
+    if classification:
+        y_arr = _refresh_map_labels(est, y_arr)
+    else:
+        y_arr = np.asarray(y_arr, dtype=np.float64)
+    prep = est.model_.prep_
+    rows = rows0.append(prep, X, y_arr, w_new)
+    return _refresh_replay(est, rows, rows.rebuild_X(prep))
+
+
 class _RegBoosterFactory:
     """Value-capturing replacement for the regressor `_fit_single`'s old
     `_fit_booster` / `_screen` closures.
@@ -2174,6 +2349,17 @@ class ChimeraBoostRegressor(RegressorMixin, BaseEstimator):
         ``predict`` takes the row groups and adds the fitted intercepts;
         unseen groups get exactly 0. Slice 1: ``loss="RMSE"`` single models
         only (not ``n_ensembles > 1``).
+    store_training_data : bool, default False
+        Keep the rows the final booster's leaves came from, so ``refresh``
+        can fold in new rows later without a from-scratch refit. What is
+        stored is compact, not raw X: plain numeric columns as binner bins,
+        the raw values of only the cross-feature parent columns,
+        categoricals as integer codes plus their per-column categories, the
+        target as the booster saw it, and the raw sample weights. The store
+        rides along in pickles and grows with every ``refresh`` call.
+        Roughly 3x smaller than float64 X for numeric data when cross
+        features engage, about 6.5x when they do not; far smaller for
+        object-dtype categoricals.
 
     Attributes
     ----------
@@ -2212,6 +2398,10 @@ class ChimeraBoostRegressor(RegressorMixin, BaseEstimator):
     group_ratio_ : float or None
         Fitted noise-to-group variance ratio (``inf`` means no group signal
         was found). ``None`` unless fit with ``random_effects=True``.
+    n_samples_trained_ : int
+        Number of rows the final booster's leaves came from: all rows when
+        the full-data refit ran, else the post-split training rows. Grown
+        by ``refresh``.
     """
 
     # Both estimators accept cross_features="always" (the unrefereed forced
@@ -2239,7 +2429,8 @@ def __init__(self, n_estimators=2000, learning_rate=None, depth=None,
                  cat_features=None, quantize_gradients=True,
                  eval_metric=None, delta=1.0, tweedie_variance_power=1.5,
                  refit_full="replay", refit_members=False, quality=None,
-                 adaptive_learning_rate=True, random_effects=False):
+                 adaptive_learning_rate=True, random_effects=False,
+                 store_training_data=False):
         self.n_estimators = n_estimators
         self.learning_rate = learning_rate
         self.depth = depth
@@ -2282,6 +2473,7 @@ def __init__(self, n_estimators=2000, learning_rate=None, depth=None,
         # consulted when learning_rate is None. False == the historical flat 0.1.
         self.adaptive_learning_rate = adaptive_learning_rate
         self.random_effects = random_effects
+        self.store_training_data = store_training_data
 
     def fit(self, X, y, cat_features=None, eval_set=None, groups=None,
             sample_weight=None, callbacks=None):
@@ -2344,6 +2536,7 @@ def fit(self, X, y, cat_features=None, eval_set=None, groups=None,
                 _check_feature_names_match(self, eval_set[0])
 
         with _quality_applied(self):
+            _check_store_fit(self)
             if self.random_effects:
                 if self.loss != "RMSE":
                     raise ValueError(
@@ -2587,6 +2780,7 @@ def _fit_single(self, X, y, cat_features, eval_set, groups, sample_weight,
                 n_re_groups)
             y_full_for_refit = (np.asarray(y_full, dtype=np.float64)
                                 - b_pre[group_codes_full])
+        pre_refit = self.model_
         self._dispatch_reg_refit(kw, loss_kwargs, X_full, y_full_for_refit,
                                  sw_full, cat_features, auto_split,
                                  cat_ctx=full_ctx)
@@ -2600,6 +2794,9 @@ def _fit_single(self, X, y, cat_features, eval_set, groups, sample_weight,
                 self.model_, X_full, y_full, sw_full, group_codes_full,
                 n_re_groups)
 
+        _capture_training_store(
+            self, (X_full, y_full, sw_full), (X, y, sample_weight),
+            self.model_ is not pre_refit, classification=False)
         return self
 
     def _arm_forced_cross(self, fb, select_ll, ll):
@@ -2810,6 +3007,40 @@ def _transform_raw(self, raw):
         tf = getattr(self.model_.loss_, "transform", None)
         return raw if tf is None else tf(raw)
 
+    def refresh(self, X, y, sample_weight=None):
+        """Fold new rows into this fitted model without a from-scratch refit.
+
+        Appends ``(X, y[, sample_weight])`` to the rows stored at fit time
+        (``store_training_data=True``), rebuilds a raw-equivalent training
+        matrix from the store, and replays the FITTED booster's configuration
+        on it -- the same structure-transfer replay the default full-data
+        refit uses. Every tree's structure, the round count and the learning
+        rate stay pinned; leaf values and linear-leaf coefficients are refit
+        against the combined gradients.
+
+        ``refresh`` replays the fitted configuration and ignores any
+        ``set_params`` made after fit. The binner borders, the count and
+        cross selections, the validation history and the fitted selections
+        are unchanged; ``n_samples_trained_`` grows by the number of new
+        rows.
+
+        Parameters
+        ----------
+        X, y : array-like
+            New rows with the same features the model was fit on. Zero rows
+            are allowed, in which case predictions are unchanged.
+        sample_weight : array-like of shape (n_samples,) or None
+            Weights for the new rows, in the same units as the fit weights.
+            Required when the model was fit with weights, rejected when it
+            was not.
+
+        Returns
+        -------
+        self
+        """
+        return _refresh_single(self, X, y, sample_weight,
+                               classification=False)
+
     def predict_raw(self, X):
         """Raw additive score before the loss link and conformal quantile offset.
 
@@ -3235,6 +3466,17 @@ class ChimeraBoostClassifier(ClassifierMixin, BaseEstimator):
         called without its own ``cat_features`` (the fit argument overrides).
         Provided as a constructor argument so ``GridSearchCV``/``Pipeline`` can
         carry it.
+    store_training_data : bool, default False
+        Keep the rows the final booster's leaves came from, so ``refresh``
+        can fold in new rows later without a from-scratch refit. What is
+        stored is compact, not raw X: plain numeric columns as binner bins,
+        the raw values of only the cross-feature parent columns,
+        categoricals as integer codes plus their per-column categories, the
+        target as the booster saw it, and the raw sample weights. The store
+        rides along in pickles and grows with every ``refresh`` call.
+        Roughly 3x smaller than float64 X for numeric data when cross
+        features engage, about 6.5x when they do not; far smaller for
+        object-dtype categoricals.
 
     Attributes
     ----------
@@ -3254,6 +3496,10 @@ class ChimeraBoostClassifier(ClassifierMixin, BaseEstimator):
         Bagged-mode member defaults that were auto-applied (params the user
         left on auto resolve to tuned member values inside a bag; explicit
         values always win). Set only when ``n_ensembles > 1``.
+    n_samples_trained_ : int
+        Number of rows the final booster's leaves came from: all rows when
+        the full-data refit ran, else the post-split training rows. Grown
+        by ``refresh``.
     """
 
     # Not pinned on the classifier: linear_leaves=None is already an auto rule
@@ -3281,7 +3527,8 @@ def __init__(self, n_estimators=2000, learning_rate=None, depth=6,
                  n_ensembles=None, ensemble_n_jobs=-1, max_samples=0.8,
                  cat_features=None, quantize_gradients=True,
                  eval_metric=None, refit_full="replay", refit_members=False,
-                 quality=None, adaptive_learning_rate=True):
+                 quality=None, adaptive_learning_rate=True,
+                 store_training_data=False):
         self.n_estimators = n_estimators
         self.learning_rate = learning_rate
         self.depth = depth
@@ -3319,6 +3566,7 @@ def __init__(self, n_estimators=2000, learning_rate=None, depth=6,
         # Size fade for the auto learning rate, default-on since 0.30.0; only
         # consulted when learning_rate is None. False == the historical flat 0.1.
         self.adaptive_learning_rate = adaptive_learning_rate
+        self.store_training_data = store_training_data
 
     def fit(self, X, y, cat_features=None, eval_set=None, groups=None,
             sample_weight=None, callbacks=None):
@@ -3375,6 +3623,7 @@ def fit(self, X, y, cat_features=None, eval_set=None, groups=None,
                 _check_eval_labels(eval_set, y)
 
         with _quality_applied(self):
+            _check_store_fit(self)
             if self.n_ensembles and self.n_ensembles > 1:
                 if callbacks is not None:
                     raise ValueError(
@@ -3668,6 +3917,7 @@ def _fit_single(self, X, y, cat_features, eval_set, groups, sample_weight,
         (X, y, sample_weight, eval_set, es_active, auto_split,
          X_full, y_full, sw_full, split_idx) = self._resolve_classes_and_split(
             X, y, cat_features, eval_set, groups, sample_weight)
+        _check_store_multiclass(self)
         full_ctx, train_ctx, val_ctx = _shared_cat_ctxs(
             X_full, split_idx, cat_features)
 
@@ -3754,11 +4004,52 @@ def _fit_single(self, X, y, cat_features, eval_set, groups, sample_weight,
             self.temperature_ = _fit_temperature(raw, cal_y, self._multiclass,
                                                  sample_weight=cal_w)
 
+        pre_refit = self.model_
         self._dispatch_cls_refit(kw, X_full, y_full, sw_full, cat_features,
                                  auto_split, cat_ctx=full_ctx)
 
+        _capture_training_store(
+            self, (X_full, y_full, sw_full), (X, y, sample_weight),
+            self.model_ is not pre_refit, classification=True)
         return self
 
+    def refresh(self, X, y, sample_weight=None):
+        """Fold new rows into this fitted binary model without a from-scratch
+        refit.
+
+        Appends ``(X, y[, sample_weight])`` to the rows stored at fit time
+        (``store_training_data=True``), rebuilds a raw-equivalent training
+        matrix from the store, and replays the FITTED booster's configuration
+        on it -- the same structure-transfer replay the default full-data
+        refit uses. Every tree's structure, the round count and the learning
+        rate stay pinned; leaf values and linear-leaf coefficients are refit
+        against the combined gradients.
+
+        ``refresh`` replays the fitted configuration and ignores any
+        ``set_params`` made after fit. The binner borders, the count and
+        cross selections, the validation history, the fitted selections and
+        the frozen ``temperature_`` are unchanged; ``n_samples_trained_``
+        grows by the number of new rows.
+
+        Parameters
+        ----------
+        X, y : array-like
+            New rows with the same features the model was fit on. ``y``
+            carries original class labels, mapped through the fitted
+            ``classes_`` -- labels never seen at fit raise. Zero rows are
+            allowed, in which case predictions are unchanged.
+        sample_weight : array-like of shape (n_samples,) or None
+            Weights for the new rows, in the same units as the fit weights.
+            Required when the model was fit with weights, rejected when it
+            was not.
+
+        Returns
+        -------
+        self
+        """
+        return _refresh_single(self, X, y, sample_weight,
+                               classification=True)
+
     def predict_proba(self, X):
         Xv = _check_predict_input(self, X)
         X = X if Xv is None else Xv
diff --git a/chimeraboost/training_rows.py b/chimeraboost/training_rows.py
new file mode 100644
index 0000000..f8822a9
--- /dev/null
+++ b/chimeraboost/training_rows.py
@@ -0,0 +1,191 @@
+"""Compact storage of the rows a fitted booster's leaves came from.
+
+Plain numeric columns are kept as bins of the fitted binner (feature-major),
+cross parents as raw float64, categoricals as int32 codes plus per-column
+categories in code order, y as float64 and the raw weights. ``rebuild_X``
+turns the store back into a raw-equivalent array that the existing replay
+refit consumes unchanged -- there is deliberately no second preprocessing
+path here.
+"""
+
+import numpy as np
+
+from .binning import BIN_DTYPE
+from .preprocessing import as_model_array
+from .target_encoding import factorize
+
+# What factorize() turns a missing categorical (None, NaN, or refusing
+# self-comparison) into; rebuild_X maps this category back to np.nan.
+_MISSING = "__nan__"
+
+
+def _split_plain_parent(prep):
+    """Split ``prep.num_features_`` into plain (binned) and parent (raw).
+
+    Every numeric column appearing in any of ``prep.cross_pairs`` stays raw,
+    whatever the op: diff/prod crosses are computed from raw floats and gdiff
+    group means are refit on the replay rows. Any future preprocessing input
+    that reads raw numeric values beyond binning (as crosses do) must be
+    added to the raw set here, or refresh loses exactness.
+    """
+    num = list(prep.num_features_)
+    inset = set(num)
+    parents = set()
+    for i, j, _op in prep.cross_pairs:
+        if i in inset:
+            parents.add(i)
+        if j in inset:
+            parents.add(j)
+    return [f for f in num if f not in parents], [f for f in num if f in parents]
+
+
+def _values_for_bins(bins, borders):
+    """A value per bin that bins back identically (binning.py:82-93)."""
+    m = len(borders)
+    out = np.empty(len(bins), dtype=np.float64)
+    if m == 0:
+        out[:] = 0.0
+        out[bins == 1] = np.nan
+        return out
+    miss = bins == m + 1
+    zero = bins == 0
+    out[zero] = np.nextafter(borders[0], -np.inf)
+    mid = ~(miss | zero)
+    out[mid] = borders[bins[mid].astype(np.int64) - 1]
+    out[miss] = np.nan
+    return out
+
+
+def _extend_codes(stored_vals, new_col, stored_codes):
+    """Codes for [stored; new] with unseen categories appended in order."""
+    ext = list(stored_vals)
+    idx = {v: i for i, v in enumerate(ext)}
+    ncodes, ncats = factorize(new_col)
+    loc = np.empty(len(ncats), dtype=np.int64)
+    for li, nc in enumerate(ncats):
+        if nc in idx:
+            loc[li] = idx[nc]
+        else:
+            idx[nc] = len(ext)
+            ext.append(nc)
+            loc[li] = len(ext) - 1
+    new_codes = loc[ncodes].astype(np.int32)
+    full = np.concatenate([stored_codes, new_codes])
+    return full, np.asarray(ext, dtype=object)
+
+
+class TrainingRows:
+    """The rows a fitted booster's leaves came from, in its row order."""
+
+    def __init__(self, n_features, cat_features, num_features, plain_features,
+                 parent_features, plain_bins, parent_raw, cat_codes,
+                 cat_values, y, sample_weight):
+        self.n_features = int(n_features)
+        self.cat_features = list(cat_features)
+        self.num_features = list(num_features)
+        self.plain_features = list(plain_features)
+        self.parent_features = list(parent_features)
+        self.plain_bins = plain_bins
+        self.parent_raw = parent_raw
+        self.cat_codes = cat_codes
+        self.cat_values = cat_values
+        self.y = y
+        self.sample_weight = sample_weight
+
+    @property
+    def n_rows(self):
+        return int(self.y.shape[0])
+
+    @classmethod
+    def capture(cls, prep, X, y, sample_weight=None):
+        """Capture raw rows with a FITTED ``FeaturePreprocessor``."""
+        X = as_model_array(X, bool(prep.cat_features_))
+        y = np.asarray(y, dtype=np.float64)
+        w = (None if sample_weight is None
+             else np.asarray(sample_weight, dtype=np.float64))
+        plain, parents = _split_plain_parent(prep)
+        pos = {f: k for k, f in enumerate(prep.num_features_)}
+        cols = [pos[f] for f in plain]
+        nb = prep.binner_.n_bins_
+        dt = (np.uint8 if all(nb[k] <= 256 for k in cols) else BIN_DTYPE)
+        if cols:
+            bins = np.ascontiguousarray(prep.transform(X)[:, cols].T, dtype=dt)
+        else:
+            bins = np.empty((0, X.shape[0]), dtype=dt)
+        pcols = [pos[f] for f in parents]
+        num = prep._numeric_block(X)
+        if pcols:
+            raw = np.ascontiguousarray(num[:, pcols], dtype=np.float64)
+        else:
+            raw = np.empty((X.shape[0], 0), dtype=np.float64)
+        codes = prep._codes_for_transform(X).astype(np.int32)
+        vals = []
+        for m in prep.cat_maps_:
+            arr = np.empty(len(m), dtype=object)
+            for v, c in m.items():
+                arr[c] = v
+            vals.append(arr)
+        return cls(X.shape[1], prep.cat_features_, prep.num_features_,
+                   plain, parents, bins, raw, codes, vals, y, w)
+
+    def append(self, prep, X_new, y_new, sample_weight_new=None):
+        """A NEW ``TrainingRows`` for [stored; new]; neither mutates."""
+        X_new = as_model_array(X_new, bool(self.cat_features))
+        y_new = np.asarray(y_new, dtype=np.float64)
+        w_new = (None if sample_weight_new is None
+                 else np.asarray(sample_weight_new, dtype=np.float64))
+        if X_new.shape[1] != self.n_features:
+            raise ValueError("X_new has %d columns, stored %d"
+                             % (X_new.shape[1], self.n_features))
+        if (self.sample_weight is None) != (w_new is None):
+            raise ValueError("weighted/unweighted append mismatch")
+        pos = {f: k for k, f in enumerate(prep.num_features_)}
+        cols = [pos[f] for f in self.plain_features]
+        if cols:
+            nb = np.ascontiguousarray(
+                prep.transform(X_new)[:, cols].T,
+                dtype=self.plain_bins.dtype)
+            bins = np.concatenate([self.plain_bins, nb], axis=1)
+        else:
+            bins = np.empty((0, self.n_rows + len(y_new)),
+                            dtype=self.plain_bins.dtype)
+        pcols = [pos[f] for f in self.parent_features]
+        num = prep._numeric_block(X_new)
+        praw = np.concatenate(
+            [self.parent_raw,
+             np.ascontiguousarray(num[:, pcols], dtype=np.float64)
+             if pcols else np.empty((len(y_new), 0))])
+        codes, vals = [], []
+        for j, f in enumerate(self.cat_features):
+            c, v = _extend_codes(self.cat_values[j], X_new[:, f],
+                                 self.cat_codes[:, j])
+            codes.append(c)
+            vals.append(v)
+        stacked = (np.column_stack(codes).astype(np.int32) if codes
+                   else np.empty((self.n_rows + len(y_new), 0),
+                                 dtype=np.int32))
+        y = np.concatenate([self.y, y_new])
+        w = (None if w_new is None
+             else np.concatenate([self.sample_weight, w_new]))
+        return TrainingRows(self.n_features, self.cat_features,
+                            self.num_features, self.plain_features,
+                            self.parent_features, bins, praw, stacked,
+                            vals, y, w)
+
+    def rebuild_X(self, prep):
+        """The raw-equivalent array in the original column layout."""
+        n = self.n_rows
+        X = (np.empty((n, self.n_features), dtype=object)
+             if self.cat_features
+             else np.empty((n, self.n_features), dtype=np.float64))
+        pos = {f: k for k, f in enumerate(prep.num_features_)}
+        for p, f in enumerate(self.plain_features):
+            X[:, f] = _values_for_bins(
+                self.plain_bins[p], prep.binner_.borders_[pos[f]])
+        for j, f in enumerate(self.parent_features):
+            X[:, f] = self.parent_raw[:, j]
+        for j, f in enumerate(self.cat_features):
+            col = self.cat_values[j][self.cat_codes[:, j].astype(np.int64)]
+            col[col == _MISSING] = np.nan
+            X[:, f] = col
+        return X
diff --git a/chimeraboost/tree.py b/chimeraboost/tree.py
index 565e31c..99f6845 100644
--- a/chimeraboost/tree.py
+++ b/chimeraboost/tree.py
@@ -868,95 +868,263 @@ def _linear_leaf_fit(leaf, grad, hess, n_leaves, lin_feats, centers_std, Xb,  #
 
     Notes
     -----
-    Parallel over leaves and bit-identical to the old serial global scan. A
-    stable counting sort groups sample indices by leaf in original order, so each
-    leaf accumulates in exactly the float-add sequence the serial version used --
-    a leaf only ever saw its own samples, in increasing i. Thread-count invariant
-    for the same reason.
-
-    Design values are gathered per sample inside the leaf loop: no (k, n) scratch
-    matrix, and one parallel region holds the JIT compile cost down.
+    Two bit-identical arms on `n <= _SMALL_N`: the small arm uses a serial
+    stable counting sort and parallelizes over leaves; the large arm sorts in
+    parallel, gathers leaf-sorted contiguous buffers once, and parallelizes
+    over (leaf, accumulator-group) tasks so a 60%-of-rows leaf spreads over
+    1+k threads. Every accumulator sums in increasing row index with the same
+    expressions, so both arms match the old serial global scan bit for bit.
     """
     n = leaf.shape[0]
     k = lin_feats.shape[0]
     d = 1 + k
-    coef = np.zeros((n_leaves, d))
+    if n <= _SMALL_N:
+        coef = np.zeros((n_leaves, d))
 
-    # Per-leaf grad/hess totals (for the constant fallback) and counts.
-    counts = np.zeros(n_leaves, dtype=np.int64)
-    Gtot = np.zeros(n_leaves)
-    Htot = np.zeros(n_leaves)
-    for i in range(n):
-        l = leaf[i]
-        counts[l] += 1
-        Gtot[l] += grad[i]
-        Htot[l] += hess[i]
+        # Per-leaf grad/hess totals (for the constant fallback) and counts.
+        counts = np.zeros(n_leaves, dtype=np.int64)
+        Gtot = np.zeros(n_leaves)
+        Htot = np.zeros(n_leaves)
+        for i in range(n):
+            l = leaf[i]
+            counts[l] += 1
+            Gtot[l] += grad[i]
+            Htot[l] += hess[i]
 
-    # Stable counting sort: order[start[l]:start[l+1]] = leaf-l samples in
-    # increasing original index.
-    start = np.zeros(n_leaves + 1, dtype=np.int64)
-    for l in range(n_leaves):
-        start[l + 1] = start[l] + counts[l]
+        # Stable counting sort: order[start[l]:start[l+1]] = leaf-l samples in
+        # increasing original index.
+        start = np.zeros(n_leaves + 1, dtype=np.int64)
+        for l in range(n_leaves):
+            start[l + 1] = start[l] + counts[l]
 
-    pos = start[:n_leaves].copy()
-    order = np.empty(n, dtype=np.int64)
-    for i in range(n):
-        l = leaf[i]
-        order[pos[l]] = i
-        pos[l] += 1
+        pos = start[:n_leaves].copy()
+        order = np.empty(n, dtype=np.int64)
+        for i in range(n):
+            l = leaf[i]
+            order[pos[l]] = i
+            pos[l] += 1
 
-    # Per-leaf normal equations + solve; leaves are independent.
-    for l in prange(n_leaves):
-        if counts[l] == 0:
-            continue
+        # Per-leaf normal equations + solve; leaves are independent.
+        for l in prange(n_leaves):
+            if counts[l] == 0:
+                continue
 
-        if counts[l] < 2 * d or k == 0:
-            if Htot[l] > 0.0:
-                coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
-            continue
+            if counts[l] < 2 * d or k == 0:
+                if Htot[l] > 0.0:
+                    coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
+                continue
 
-        Ml = np.zeros((d, d))
-        rl = np.zeros(d)
-        xrow = np.empty(k)
+            Ml = np.zeros((d, d))
+            rl = np.zeros(d)
+            xrow = np.empty(k)
 
-        for q in range(start[l], start[l + 1]):
-            i = order[q]
-            h = hess[i]
-            g = grad[i]
+            for q in range(start[l], start[l + 1]):
+                i = order[q]
+                h = hess[i]
+                g = grad[i]
+
+                # Standardized design values for this sample; missing bins -> 0.
+                for j in range(k):
+                    f = lin_feats[j]
+                    v = centers_std[f, Xb[f, i]]
+                    xrow[j] = v if np.isfinite(v) else 0.0
+
+                Ml[0, 0] += h
+                rl[0] += -g
+                for j in range(k):
+                    xj = xrow[j]
+                    Ml[0, 1 + j] += h * xj
+                    Ml[1 + j, 0] += h * xj
+                    rl[1 + j] += -g * xj
+                    for jj in range(k):
+                        Ml[1 + j, 1 + jj] += h * xj * xrow[jj]
+
+            Ml[0, 0] += l2_intercept
+            for j in range(1, d):
+                Ml[j, j] += lin_lambda
+            for j in range(d):
+                Ml[j, j] += 1e-9              # jitter: keep the solve well-posed
+
+            beta = _solve_small(Ml, rl)
+            if np.isnan(beta[0]):
+                # Singular pivot (unreachable given the diagonal ridge + jitter):
+                # keep the plain constant Newton value rather than a broken slope.
+                if Htot[l] > 0.0:
+                    coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
+                continue
+
+            for j in range(d):
+                coef[l, j] = lr * beta[j]
+    else:
 
-            # Standardized design values for this sample; missing bins -> 0.
+        # Parallel stable counting sort: per-chunk leaf counts, prefix sums in
+        # (leaf, chunk) order, each chunk scatters its rows in increasing index.
+        n_chunks = (n + 4095) // 4096
+        if n_chunks < 1:
+            n_chunks = 1
+        if n_chunks > 64:
+            n_chunks = 64
+        cnt = np.zeros((n_chunks, n_leaves), dtype=np.int64)
+        for c in prange(n_chunks):
+            lo = c * n // n_chunks
+            hi = (c + 1) * n // n_chunks
+            for i in range(lo, hi):
+                cnt[c, leaf[i]] += 1
+        counts = np.zeros(n_leaves, dtype=np.int64)
+        for l in range(n_leaves):
+            s = 0
+            for c in range(n_chunks):
+                s += cnt[c, l]
+            counts[l] = s
+        start = np.zeros(n_leaves + 1, dtype=np.int64)
+        for l in range(n_leaves):
+            start[l + 1] = start[l] + counts[l]
+        cursor = np.empty((n_chunks, n_leaves), dtype=np.int64)
+        for l in range(n_leaves):
+            s = start[l]
+            for c in range(n_chunks):
+                cursor[c, l] = s
+                s += cnt[c, l]
+        order = np.empty(n, dtype=np.int64)
+        for c in prange(n_chunks):
+            lo = c * n // n_chunks
+            hi = (c + 1) * n // n_chunks
+            for i in range(lo, hi):
+                l = leaf[i]
+                order[cursor[c, l]] = i
+                cursor[c, l] += 1
+
+        # One parallel gather into leaf-sorted contiguous buffers.
+        gs = np.empty(n, dtype=np.float64)
+        hs = np.empty(n, dtype=np.float64)
+        if k > 0:
+            Xd = np.empty((k, n), dtype=np.float64)
+        else:
+            Xd = np.empty((1, n), dtype=np.float64)
+        for q in prange(n):
+            i = order[q]
+            gs[q] = grad[i]
+            hs[q] = hess[i]
             for j in range(k):
                 f = lin_feats[j]
                 v = centers_std[f, Xb[f, i]]
-                xrow[j] = v if np.isfinite(v) else 0.0
+                Xd[j, q] = v if np.isfinite(v) else 0.0
 
-            Ml[0, 0] += h
-            rl[0] += -g
-            for j in range(k):
-                xj = xrow[j]
-                Ml[0, 1 + j] += h * xj
-                Ml[1 + j, 0] += h * xj
-                rl[1 + j] += -g * xj
+        # Task map, accumulator-major so a big leaf's groups spread over threads:
+        # acc 0 = fallback (Gtot, Htot), acc 1 = full intercept group,
+        # acc 2+j = full slope group for feature j.
+        n_full = 0
+        n_fallback = 0
+        for l in range(n_leaves):
+            if counts[l] == 0:
+                continue
+            if counts[l] >= 2 * d and k > 0:
+                n_full += 1
+            else:
+                n_fallback += 1
+        T = n_fallback + n_full * (1 + k)
+        task_leaf = np.empty(T, dtype=np.int64)
+        task_acc = np.empty(T, dtype=np.int64)
+        p = 0
+        for l in range(n_leaves):
+            if counts[l] > 0 and not (counts[l] >= 2 * d and k > 0):
+                task_leaf[p] = l
+                task_acc[p] = 0
+                p += 1
+        for l in range(n_leaves):
+            if counts[l] >= 2 * d and k > 0:
+                task_leaf[p] = l
+                task_acc[p] = 1
+                p += 1
+        for j in range(k):
+            for l in range(n_leaves):
+                if counts[l] >= 2 * d and k > 0:
+                    task_leaf[p] = l
+                    task_acc[p] = 2 + j
+                    p += 1
+
+        Gtot = np.zeros(n_leaves)
+        Htot = np.zeros(n_leaves)
+        Ml_all = np.zeros((n_leaves, d, d))
+        rl_all = np.zeros((n_leaves, d))
+        for t in prange(T):
+            l = task_leaf[t]
+            acc = task_acc[t]
+            lo = start[l]
+            hi = start[l + 1]
+            if acc == 0:
+                sG = 0.0
+                sH = 0.0
+                for q in range(lo, hi):
+                    sG += gs[q]
+                    sH += hs[q]
+                Gtot[l] = sG
+                Htot[l] = sH
+            elif acc == 1:
+                sG = 0.0
+                sH = 0.0
+                sM00 = 0.0
+                sr0 = 0.0
+                for q in range(lo, hi):
+                    g = gs[q]
+                    h = hs[q]
+                    sG += g
+                    sH += h
+                    sM00 += h
+                    sr0 += -g
+                Gtot[l] = sG
+                Htot[l] = sH
+                Ml_all[l, 0, 0] = sM00
+                rl_all[l, 0] = sr0
+            else:
+                j = acc - 2
+                s01 = 0.0
+                s10 = 0.0
+                sr = 0.0
+                sblk = np.zeros(k)
+                for q in range(lo, hi):
+                    h = hs[q]
+                    g = gs[q]
+                    xj = Xd[j, q]
+                    hx = h * xj
+                    s01 += hx
+                    s10 += hx
+                    sr += -g * xj
+                    for jj in range(k):
+                        sblk[jj] += hx * Xd[jj, q]
+                Ml_all[l, 0, 1 + j] = s01
+                Ml_all[l, 1 + j, 0] = s10
+                rl_all[l, 1 + j] = sr
                 for jj in range(k):
-                    Ml[1 + j, 1 + jj] += h * xj * xrow[jj]
-
-        Ml[0, 0] += l2_intercept
-        for j in range(1, d):
-            Ml[j, j] += lin_lambda
-        for j in range(d):
-            Ml[j, j] += 1e-9              # jitter: keep the solve well-posed
-
-        beta = _solve_small(Ml, rl)
-        if np.isnan(beta[0]):
-            # Singular pivot (unreachable given the diagonal ridge + jitter):
-            # keep the plain constant Newton value rather than a broken slope.
-            if Htot[l] > 0.0:
-                coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
-            continue
-
-        for j in range(d):
-            coef[l, j] = lr * beta[j]
+                    Ml_all[l, 1 + j, 1 + jj] = sblk[jj]
 
+        # Same regularization, solve and fallbacks as the small-n arm.
+        coef = np.zeros((n_leaves, d))
+        for l in range(n_leaves):
+            if counts[l] == 0:
+                continue
+            if counts[l] < 2 * d or k == 0:
+                if Htot[l] > 0.0:
+                    coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
+                continue
+            Ml = np.empty((d, d))
+            rl = np.empty(d)
+            for a in range(d):
+                rl[a] = rl_all[l, a]
+                for b in range(d):
+                    Ml[a, b] = Ml_all[l, a, b]
+            Ml[0, 0] += l2_intercept
+            for j in range(1, d):
+                Ml[j, j] += lin_lambda
+            for j in range(d):
+                Ml[j, j] += 1e-9
+            beta = _solve_small(Ml, rl)
+            if np.isnan(beta[0]):
+                if Htot[l] > 0.0:
+                    coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
+                continue
+            for j in range(d):
+                coef[l, j] = lr * beta[j]
     return coef
 
 
@@ -2146,7 +2314,7 @@ def replay_oblivious_tree(donor, Xb, grad, hess, l2, lr, linear_leaves=False,
         return (ObliviousTree(sf, st, np.zeros(1), np.zeros(0)),
                 np.zeros(Xb.shape[1], dtype=np.int64))
 
-    leaf = _assign_leaves(Xb, sf, st)
+    leaf = donor.apply(Xb)
     n_leaves = 1 << len(sf)
     values = _leaf_values(leaf, grad, hess, n_leaves, l2, lr)
 
diff --git a/docs/parameters.md b/docs/parameters.md
index 7ca0a95..1d1384f 100644
--- a/docs/parameters.md
+++ b/docs/parameters.md
@@ -138,6 +138,7 @@ See the User Guide: [early stopping](recipes.md#early-stopping) for `eval_set` a
 | Parameter | Default | Effect |
 |---|---|---|
 | `random_effects` | `False` | Fit one shrunk intercept per `groups=` label on top of the trees (regressor, RMSE single-model only). Small groups pool toward zero; unseen groups predict trees-only. See the User Guide: [grouped data](recipes.md#grouped-data-random-intercepts). |
+| `store_training_data` | `False` | Keep the rows the final model was fit on, in compact form, so `refresh(X, y)` can fold in new rows later by refitting only the leaf values. Regressor and binary classifier, single models only. The stored rows travel with the pickled model and grow with every refresh. See the User Guide: [refreshing with new rows](recipes.md#refreshing-with-new-rows). |
 
 ## System
 
diff --git a/docs/recipes.md b/docs/recipes.md
index 523e875..1157800 100644
--- a/docs/recipes.md
+++ b/docs/recipes.md
@@ -361,6 +361,45 @@ reg = ChimeraBoostRegressor(refit_full=False, random_state=0).fit(X_train, y_tra
 reg = ChimeraBoostRegressor(quality=2, random_state=0).fit(X_train, y_train)
 ```
 
+## Refreshing with new rows
+
+When new labelled rows arrive after a fit, `refresh` folds them in without starting
+over. It keeps every tree's splits and refits only the leaf values on the old rows plus
+the new ones, which is the same replay the default full-data refit uses, so it costs
+about what that refit costs: roughly a third of growing the model again. The model has
+to keep its training rows for this, so ask for it at fit time:
+
+```python
+reg = ChimeraBoostRegressor(store_training_data=True, random_state=0)
+reg.fit(X_old, y_old)
+reg.refresh(X_new, y_new)       # leaves refit on the old and new rows together
+reg.n_samples_trained_          # how many rows the leaves now come from
+```
+
+What a refresh leaves alone:
+
+- The trees' splits, the number of trees and the learning rate. A refresh cannot add a
+  split the original data never called for, so a shift that needs new splits still
+  needs a refit.
+- The binary classifier's probability calibration (`temperature_`).
+- Every setting chosen at fit time. `set_params` after the fit does not change what a
+  refresh does.
+
+Worth knowing:
+
+- The stored rows take much less room than X itself. Plain numeric columns are kept as
+  bin indices, only the few columns that feed cross features keep their raw values, and
+  categories are kept as integer codes. For numeric data expect about a third of the size
+  of X as 64-bit floats. The store travels with the pickled model and grows with every
+  refresh.
+- If the model was fit with `sample_weight`, pass weights for the new rows too, in the
+  same units.
+- The classifier cannot learn new classes this way. Refit instead.
+- Fit with an integer `random_state` if you need refreshes to be exactly reproducible.
+- Bagged models (`n_ensembles > 1`, including `quality=4` and `5`), multiclass
+  classification, `loss="Quantile"` and `random_effects=True` do not support
+  `store_training_data` yet.
+
 ## Early stopping
 
 Early stopping is on by default. With no `eval_set`, the estimator holds out a
diff --git a/tests/test_refresh.py b/tests/test_refresh.py
new file mode 100644
index 0000000..674b0dc
--- /dev/null
+++ b/tests/test_refresh.py
@@ -0,0 +1,778 @@
+"""Stored training rows round-trip (issue #131 slice 1a2).
+
+A fitted booster's rows are kept as binner bins (plain numerics), raw
+floats (cross parents) and int codes (categoricals); ``rebuild_X`` turns
+them back into a raw-equivalent array the existing replay refit consumes
+unchanged. Every comparison here is bit for bit.
+
+Slice 1b (below the line) covers the public ``store_training_data`` /
+``refresh`` API; every comparison there is bit for bit too, never allclose.
+"""
+
+import inspect
+import pickle
+
+import numpy as np
+import pytest
+from sklearn.base import clone
+from sklearn.exceptions import NotFittedError
+
+from chimeraboost import (ChimeraBoostClassifier, ChimeraBoostRegressor,
+                          CustomObjective)
+from chimeraboost.booster import GradientBoosting, _BaseBooster
+from chimeraboost.preprocessing import FeaturePreprocessor
+from chimeraboost.target_encoding import factorize
+from chimeraboost.training_rows import TrainingRows
+
+CAT = [6, 7]
+N_NUM = 6
+
+
+def _cats(rng, n, levels, p_nan):
+    c = np.array(rng.choice(levels, n), dtype=object)
+    c[rng.random(n) < p_nan] = np.nan
+    return c
+
+
+def _target(rng, num, hi, lo):
+    x = np.where(np.isfinite(num), num, 0.0)
+    base = (3.0 * (x[:, 0] > x[:, 1]) + x[:, 2] * x[:, 3]
+            + 0.5 * x[:, 4] - 0.3 * x[:, 5])
+    lo_eff = np.array([{"a": 0.0, "b": 0.5, "c": -0.5, "d": 1.0, "e": -1.0}.get(
+        v, 0.0) for v in lo])
+    hi_eff = np.zeros(len(lo))
+    for k, v in enumerate(hi):
+        if isinstance(v, str) and v.startswith("L"):
+            try:
+                hi_eff[k] = (int(v[1:]) % 10) * 0.2
+            except ValueError:
+                pass
+    return base + lo_eff + hi_eff + 0.1 * rng.standard_normal(len(lo))
+
+
+def _weights(rng, n):
+    w = rng.uniform(0.5, 1.5, n)
+    w[rng.random(n) < 0.01] = 0.0
+    return w
+
+
+def _frame(num, hi, lo):
+    X = np.empty((len(hi), N_NUM + 2), dtype=object)
+    X[:, :N_NUM] = num
+    X[:, 6] = hi
+    X[:, 7] = lo
+    return X
+
+
+def _make_data(seed=0):
+    rng = np.random.default_rng(seed)
+    n_full, n_held, n_new = 2600, 400, 300
+    num = rng.standard_normal((n_full + n_held + n_new, N_NUM))
+    num[rng.random(num.shape) < 0.05] = np.nan
+    hi_levels = [f"L{i}" for i in range(350)]
+    lo_levels = ["a", "b", "c", "d", "e"]
+    hi_full = _cats(rng, n_full, hi_levels, 0.02)
+    lo_full = _cats(rng, n_full, lo_levels, 0.02)
+    hi_held = _cats(rng, n_held, hi_levels, 0.02)
+    lo_held = _cats(rng, n_held, lo_levels, 0.02)
+    hi_new = _cats(rng, n_new, hi_levels, 0.05)
+    lo_new = _cats(rng, n_new, lo_levels, 0.05)
+    unseen = rng.random(n_new) < 0.25
+    hi_new[unseen] = rng.choice([f"NEW{i}" for i in range(10)], unseen.sum())
+    lo_new[rng.random(n_new) < 0.15] = "z"
+    n0, n1 = n_full, n_full + n_held
+    d = {}
+    d["X_full"] = _frame(num[:n_full], hi_full, lo_full)
+    d["X_held"] = _frame(num[n0:n1], hi_held, lo_held)
+    d["X_new"] = _frame(num[n1:], hi_new, lo_new)
+    d["y_full"] = _target(rng, num[:n_full], hi_full, lo_full)
+    d["y_held"] = _target(rng, num[n0:n1], hi_held, lo_held)
+    d["y_new"] = _target(rng, num[n1:], hi_new, lo_new)
+    d["w_full"] = _weights(rng, n_full)
+    d["w_new"] = _weights(rng, n_new)
+    return d
+
+
+@pytest.fixture(scope="module")
+def donor():
+    """A realistic fitted donor: forced cross block, count column, combos."""
+    d = _make_data()
+    est = ChimeraBoostRegressor(random_state=0, n_estimators=100,
+                               cross_features="always", cat_combinations=True,
+                               linear_leaves=True)
+    est.fit(d["X_full"], d["y_full"], cat_features=CAT,
+            sample_weight=d["w_full"])
+    b, prep = est.model_, est.model_.prep_
+    assert len(b.trees_) > 0
+    assert any(t.lin_coef is not None for t in b.trees_)
+    assert prep.count_features_ == [6]
+    assert len(prep.cat_maps_[0]) >= 256
+    ops = {op for _, _, op in b.cross_pairs}
+    assert ops == {"diff", "prod", "gdiff"}, b.cross_pairs
+    # Default refit_full="replay" retrains on every row: the RMSE init is the
+    # FULL-data weighted mean, not the 80% split's.
+    assert b.init_ == pytest.approx(
+        float(np.average(d["y_full"], weights=d["w_full"])))
+    d["est"], d["booster"], d["prep"] = est, b, prep
+    d["rows"] = TrainingRows.capture(prep, d["X_full"], d["y_full"],
+                                     d["w_full"])
+    assert d["rows"].plain_features and d["rows"].parent_features
+    return d
+
+
+def _pinned_prep(donor_booster, donor_prep):
+    """A fresh prep with the replay-refit pinning (_prep_or_replay_matrices)."""
+    fresh = FeaturePreprocessor(
+        donor_booster.max_bins, donor_booster.cat_smoothing,
+        donor_booster.random_state, donor_booster.cat_n_permutations,
+        bool(donor_prep.combo_pairs_), donor_booster.cross_pairs,
+        donor_booster.cat_count_features)
+    fresh._pinned_count_features = list(donor_prep.count_features_)
+    fresh._pinned_cat_counts = (donor_prep.cat_maps_, donor_prep.cat_counts_)
+    return fresh
+
+
+def _fit_pinned(donor, X, y, w):
+    prep = _pinned_prep(donor["booster"], donor["prep"])
+    mat = prep.fit_transform(X, [y], CAT, w,
+                             binner=donor["prep"].binner_)
+    return prep, mat
+
+
+def _assert_prep_equal(p1, p2):
+    assert p1.cat_maps_ == p2.cat_maps_
+    assert p1.combo_pairs_ == p2.combo_pairs_
+    assert p1.combo_maps_ == p2.combo_maps_
+    assert len(p1.gdiff_maps_) == len(p2.gdiff_maps_)
+    for (d1, g1), (d2, g2) in zip(p1.gdiff_maps_, p2.gdiff_maps_):
+        assert d1 == d2
+        assert g1 == g2
+    assert p1.count_features_ == p2.count_features_
+    assert len(p1.cat_counts_) == len(p2.cat_counts_)
+    for c1, c2 in zip(p1.cat_counts_, p2.cat_counts_):
+        np.testing.assert_array_equal(c1, c2)
+    assert len(p1.encoders_) == len(p2.encoders_)
+    for e1, e2 in zip(p1.encoders_, p2.encoders_):
+        assert e1.prior_ == e2.prior_
+        assert e1.n_cat_ == e2.n_cat_
+        for a1, a2 in zip(e1.sums_, e2.sums_):
+            np.testing.assert_array_equal(a1, a2)
+        for a1, a2 in zip(e1.counts_, e2.counts_):
+            np.testing.assert_array_equal(a1, a2)
+    np.testing.assert_array_equal(p1.feature_map_, p2.feature_map_)
+    np.testing.assert_array_equal(p1.is_numeric_binned_, p2.is_numeric_binned_)
+    np.testing.assert_array_equal(p1.n_bins_, p2.n_bins_)
+
+
+def test_round_trip_same_rows(donor):
+    """(a) Rebuilt X refits exactly like the raw rows, same rows."""
+    rows = donor["rows"]
+    assert rows.n_rows == len(donor["y_full"])
+    X_rb = rows.rebuild_X(donor["prep"])
+    assert X_rb.shape == donor["X_full"].shape
+    assert X_rb.dtype == object
+    w = GradientBoosting._normalize_weights(donor["w_full"],
+                                            len(donor["y_full"]))
+    p1, m1 = _fit_pinned(donor, X_rb, donor["y_full"], w)
+    p2, m2 = _fit_pinned(donor, donor["X_full"], donor["y_full"], w)
+    assert m1.dtype == m2.dtype
+    np.testing.assert_array_equal(m1, m2)
+    _assert_prep_equal(p1, p2)
+    np.testing.assert_array_equal(p1.transform(donor["X_held"]),
+                                 p2.transform(donor["X_held"]))
+
+
+def test_round_trip_appended_rows(donor):
+    """(b) The same with new rows appended through ``append``."""
+    prep = donor["prep"]
+    before = ([dict(m) for m in prep.cat_maps_],
+              [(dict(d), g) for d, g in prep.gdiff_maps_],
+              [b.copy() for b in prep.binner_.borders_])
+    rows = donor["rows"]
+    n_old = rows.n_rows
+    rows2 = rows.append(prep, donor["X_new"], donor["y_new"], donor["w_new"])
+    assert rows.n_rows == n_old
+    assert rows2.n_rows == n_old + len(donor["y_new"])
+    assert [dict(m) for m in prep.cat_maps_] == before[0]
+    assert [(dict(d), g) for d, g in prep.gdiff_maps_] == before[1]
+    for b, b0 in zip(prep.binner_.borders_, before[2]):
+        np.testing.assert_array_equal(b, b0)
+    X_cat = np.concatenate([donor["X_full"], donor["X_new"]])
+    y_cat = np.concatenate([donor["y_full"], donor["y_new"]])
+    w_cat = np.concatenate([donor["w_full"], donor["w_new"]])
+    for j, f in enumerate(CAT):
+        codes, cats = factorize(X_cat[:, f])
+        np.testing.assert_array_equal(rows2.cat_codes[:, j], codes)
+        assert list(rows2.cat_values[j]) == list(cats)
+    X_rb = rows2.rebuild_X(prep)
+    assert X_rb.shape == X_cat.shape
+    w = GradientBoosting._normalize_weights(w_cat, len(y_cat))
+    p1, m1 = _fit_pinned(donor, X_rb, y_cat, w)
+    p2, m2 = _fit_pinned(donor, X_cat, y_cat, w)
+    np.testing.assert_array_equal(m1, m2)
+    _assert_prep_equal(p1, p2)
+    np.testing.assert_array_equal(p1.transform(donor["X_held"]),
+                                 p2.transform(donor["X_held"]))
+
+
+def _fit_replay(donor, X, linear):
+    kw = donor["booster"].replay_kwargs()
+    kw["n_estimators"] = len(donor["booster"].trees_)
+    kw["learning_rate"] = float(donor["booster"].lr_)
+    kw["early_stopping_rounds"] = None
+    kw["replay_donor"] = (donor["booster"].trees_, donor["prep"])
+    if linear:
+        kw["linear_leaves"] = True
+    b = GradientBoosting(**kw)
+    b.fit(X, donor["y_full"], cat_features=CAT,
+          sample_weight=donor["w_full"])
+    return b
+
+
+def _assert_booster_equal(b1, b2):
+    assert b1.init_ == b2.init_
+    assert len(b1.trees_) == len(b2.trees_)
+    for t1, t2 in zip(b1.trees_, b2.trees_):
+        np.testing.assert_array_equal(t1.values, t2.values)
+        if t1.lin_coef is None or t2.lin_coef is None:
+            assert t1.lin_coef is None and t2.lin_coef is None
+        else:
+            np.testing.assert_array_equal(t1.lin_coef, t2.lin_coef)
+
+
+@pytest.mark.parametrize("linear", [False, True])
+def test_booster_replay_on_rebuilt(donor, linear):
+    """(c) Replay refit on rebuilt X equals the one on raw rows.
+
+    False is the donor as fitted (linear leaves), True forces
+    ``linear_leaves=True``; replay refits ``lin_coef`` off a linear donor.
+    """
+    X_rb = donor["rows"].rebuild_X(donor["prep"])
+    b1 = _fit_replay(donor, X_rb, linear)
+    b2 = _fit_replay(donor, donor["X_full"], linear)
+    assert any(t.lin_coef is not None for t in b1.trees_)
+    _assert_booster_equal(b1, b2)
+    np.testing.assert_array_equal(b1.predict_raw(donor["X_held"]),
+                                 b2.predict_raw(donor["X_held"]))
+
+
+def test_exhaustive_bin_round_trip(donor):
+    """(d) Every plain-numeric bin maps to a value that bins back to it."""
+    binner = donor["prep"].binner_
+    pos = {f: k for k, f in enumerate(donor["prep"].num_features_)}
+    for f in donor["rows"].plain_features:
+        borders = binner.borders_[pos[f]]
+        m = len(borders)
+        vals = np.empty(m + 2, dtype=np.float64)
+        if m == 0:
+            vals[0], vals[1] = 0.0, np.nan
+        else:
+            vals[0] = np.nextafter(borders[0], -np.inf)
+            vals[1:m + 1] = borders
+            vals[m + 1] = np.nan
+        feat = np.zeros((m + 2, len(binner.borders_)))
+        feat[:, pos[f]] = vals
+        got = binner.transform(feat)[:, pos[f]]
+        np.testing.assert_array_equal(got, np.arange(m + 2))
+
+
+def test_replay_kwargs_covers_init(donor):
+    """(e) ``replay_kwargs`` covers every booster init parameter."""
+    kw = donor["booster"].replay_kwargs()
+    base = set(inspect.signature(_BaseBooster.__init__).parameters) - {"self"}
+    assert set(kw) == base | {"loss", "loss_kwargs"}
+    assert kw["loss"] == donor["booster"].loss_name
+    assert kw["loss_kwargs"] == donor["booster"].loss_kwargs
+    assert kw["cross_pairs"] == donor["booster"].cross_pairs
+    assert kw["cross_pairs"] is not donor["booster"].cross_pairs
+    assert kw["loss_kwargs"] is not donor["booster"].loss_kwargs
+    b2 = GradientBoosting(**kw)
+    assert b2.loss_name == donor["booster"].loss_name
+    assert b2.cross_pairs == donor["booster"].cross_pairs
+
+
+def test_training_rows_pickle_and_dtype(donor):
+    """(f) ``TrainingRows`` pickles exactly; plain bins are uint8."""
+    rows = donor["rows"]
+    assert rows.plain_bins.dtype == np.uint8
+    rt = pickle.loads(pickle.dumps(rows))
+    assert rt.n_features == rows.n_features
+    assert rt.cat_features == rows.cat_features
+    assert rt.num_features == rows.num_features
+    assert rt.plain_features == rows.plain_features
+    assert rt.parent_features == rows.parent_features
+    assert rt.plain_bins.dtype == rows.plain_bins.dtype
+    np.testing.assert_array_equal(rt.plain_bins, rows.plain_bins)
+    np.testing.assert_array_equal(rt.parent_raw, rows.parent_raw)
+    np.testing.assert_array_equal(rt.cat_codes, rows.cat_codes)
+    for v1, v2 in zip(rt.cat_values, rows.cat_values):
+        np.testing.assert_array_equal(v1, v2)
+    np.testing.assert_array_equal(rt.y, rows.y)
+    np.testing.assert_array_equal(rt.sample_weight, rows.sample_weight)
+    assert rt.n_rows == rows.n_rows
+
+
+# ---------------------------------------------------------------------------
+# Slice 1b: the public store_training_data / refresh API (issue #131).
+# ---------------------------------------------------------------------------
+
+
+class _SqErr(CustomObjective):
+    """Squared error through the public custom-objective hook.
+
+    Module level, so instances are picklable.
+    """
+
+    def init(self, y, sample_weight=None):
+        return float(np.average(y, weights=sample_weight))
+
+    def grad_hess(self, y, raw):
+        return raw - y, np.ones_like(raw)
+
+    def eval(self, y, raw, sample_weight=None):
+        return float(np.sqrt(np.average((raw - y) ** 2,
+                                        weights=sample_weight)))
+
+
+def _binary_labels(y):
+    return np.where(np.asarray(y) > np.median(y), "pos", "neg")
+
+
+def _snapshot(est, Xh):
+    """Predictions, init and per-tree leaves to compare a refresh against."""
+    snap = {"pred": est.predict(Xh), "init": est.model_.init_,
+            "n": est.n_samples_trained_,
+            "trees": [(t.values.copy(),
+                       None if t.lin_coef is None else t.lin_coef.copy())
+                      for t in est.model_.trees_]}
+    if hasattr(est, "predict_proba"):
+        snap["proba"] = est.predict_proba(Xh)
+    return snap
+
+
+def _assert_snapshot_equal(est, snap, Xh):
+    np.testing.assert_array_equal(est.predict(Xh), snap["pred"])
+    if "proba" in snap:
+        np.testing.assert_array_equal(est.predict_proba(Xh), snap["proba"])
+    assert est.model_.init_ == snap["init"]
+    assert len(est.model_.trees_) == len(snap["trees"])
+    for t, (v, lc) in zip(est.model_.trees_, snap["trees"]):
+        np.testing.assert_array_equal(t.values, v)
+        if lc is None or t.lin_coef is None:
+            assert lc is None and t.lin_coef is None
+        else:
+            np.testing.assert_array_equal(t.lin_coef, lc)
+    assert est.n_samples_trained_ == snap["n"]
+
+
+# (a) estimator-kwarg overrides per configuration. Every fit also pins an
+# integer random_state, store_training_data=True, and (regressor) the
+# single-fit linear_leaves=False / cross_features=False unless the
+# configuration says otherwise. NaNs ride in every configuration: _make_data
+# salts 5% of the numerics and 2-5% of the categoricals with NaN.
+_A_EST = {
+    "refit-replay": dict(refit_full="replay"),
+    "refit-true": dict(refit_full=True),
+    "refit-false": dict(refit_full=False),
+    "eval-set": {},
+    "no-es": dict(early_stopping=False, n_estimators=20),
+    "count": {},
+    "cross-always": dict(cross_features="always"),
+    "linear": dict(linear_leaves=True),
+    "weighted": {},
+    "ordered": dict(ordered_boosting=True),
+    "mae": dict(loss="MAE"),
+    "poisson": dict(loss="Poisson"),
+    "custom": dict(loss=_SqErr()),
+    "subsample": dict(subsample=0.8, colsample=0.7),
+    "binary": {},
+}
+
+
+def _cfg_parts(cfg, d):
+    """The (a) estimator, fit kwargs and fit target/weights for one config."""
+    base = dict(random_state=0, n_estimators=50, store_training_data=True,
+                linear_leaves=False, cross_features=False)
+    base.update(_A_EST[cfg])
+    fit_kw = dict(cat_features=CAT)
+    y_fit, w_fit = d["y_full"], None
+    cls = ChimeraBoostRegressor
+    if cfg == "eval-set":
+        fit_kw["eval_set"] = (d["X_held"], d["y_held"])
+    if cfg in ("weighted", "binary"):
+        fit_kw["sample_weight"] = d["w_full"]
+        w_fit = d["w_full"]
+    if cfg == "poisson":
+        y_fit = np.exp((y_fit - y_fit.mean()) / y_fit.std())
+    if cfg == "binary":
+        cls = ChimeraBoostClassifier
+        del base["linear_leaves"]
+        y_fit = _binary_labels(d["y_full"])
+    return cls(**base), fit_kw, y_fit, w_fit
+
+
+@pytest.mark.parametrize("cfg", ["refit-replay", "refit-true", "refit-false",
+                                 "eval-set", "no-es", "count", "cross-always",
+                                 "linear", "weighted", "ordered", "mae",
+                                 "poisson", "custom", "subsample", "binary"])
+def test_zero_rows_is_identity(cfg):
+    """(a) refresh() with zero new rows changes nothing, bit for bit."""
+    d = _make_data()
+    est, fit_kw, y_fit, w_fit = _cfg_parts(cfg, d)
+    est.fit(d["X_full"], y_fit, **fit_kw)
+    if cfg == "count":
+        assert est.model_.prep_.count_features_ == [6]
+    if cfg == "cross-always":
+        ops = {op for _, _, op in est.model_.cross_pairs}
+        assert ops == {"diff", "prod", "gdiff"}, est.model_.cross_pairs
+    if cfg == "linear":
+        assert any(t.lin_coef is not None for t in est.model_.trees_)
+    snap = _snapshot(est, d["X_held"])
+    w0 = None if w_fit is None else w_fit[:0]
+    out = est.refresh(d["X_full"][:0], y_fit[:0], sample_weight=w0)
+    assert out is est
+    _assert_snapshot_equal(est, snap, d["X_held"])
+
+
+def test_refresh_equals_manual_replay():
+    """(b) refresh() equals a manual replay refit on the concatenated rows.
+
+    The new rows carry unseen categories, NaN categories and NaN numerics
+    (see _make_data).
+    """
+    d = _make_data()
+    est = ChimeraBoostRegressor(random_state=0, n_estimators=60,
+                                cross_features="always", linear_leaves=True,
+                                store_training_data=True)
+    est.fit(d["X_full"], d["y_full"], cat_features=CAT,
+            sample_weight=d["w_full"])
+    donor = est.model_
+    assert donor.prep_.count_features_ == [6]
+    est.refresh(d["X_new"], d["y_new"], sample_weight=d["w_new"])
+    assert est.n_samples_trained_ == len(d["y_full"]) + len(d["y_new"])
+
+    kw = donor.replay_kwargs()
+    kw["n_estimators"] = len(donor.trees_)
+    kw["learning_rate"] = float(donor.lr_)
+    kw["early_stopping_rounds"] = None
+    kw["replay_donor"] = (donor.trees_, donor.prep_)
+    man = GradientBoosting(**kw)
+    man.fit(np.concatenate([d["X_full"], d["X_new"]]),
+            np.concatenate([d["y_full"], d["y_new"]]),
+            cat_features=CAT,
+            sample_weight=np.concatenate([d["w_full"], d["w_new"]]))
+    _assert_booster_equal(est.model_, man)
+    np.testing.assert_array_equal(est.predict_raw(d["X_held"]),
+                                  man.predict_raw(d["X_held"]))
+
+
+def _pinned_snapshot(est):
+    """Every refresh-pinned quantity the (c) test compares."""
+    m = est.model_
+    return {
+        "n_estimators": est.n_estimators,
+        "best": est.best_iteration_,
+        "lr": m.lr_,
+        "depths": [t.depth for t in m.trees_],
+        "feats": [np.asarray(t.splits_feat).copy() for t in m.trees_],
+        "thrs": [np.asarray(t.splits_thr).copy() for t in m.trees_],
+        "lin_feats": [(None if t.lin_feats is None
+                       else np.asarray(t.lin_feats).copy())
+                      for t in m.trees_],
+        "borders": [b.copy() for b in m.prep_.binner_.borders_],
+        "count": list(m.prep_.count_features_),
+        "cross_pairs": list(m.cross_pairs),
+        "valid": list(m.valid_history_),
+        "n_features": est.n_features_in_,
+        "names": getattr(est, "feature_names_in_", None),
+        "ll_sel": getattr(est, "linear_leaves_selected_", None),
+        "cx_sel": est.cross_features_selected_,
+        "cx_pairs": est.cross_pairs_,
+        "temp": getattr(est, "temperature_", None),
+        "classes": getattr(est, "classes_", None),
+    }
+
+
+def _assert_pinned_equal(est, snap):
+    m = est.model_
+    assert est.n_estimators == snap["n_estimators"]
+    assert est.best_iteration_ == snap["best"]
+    assert m.lr_ == snap["lr"]
+    assert [t.depth for t in m.trees_] == snap["depths"]
+    assert len(m.trees_) == len(snap["feats"])
+    for t, f, th, lf in zip(m.trees_, snap["feats"], snap["thrs"],
+                            snap["lin_feats"]):
+        np.testing.assert_array_equal(np.asarray(t.splits_feat), f)
+        np.testing.assert_array_equal(np.asarray(t.splits_thr), th)
+        if lf is None or t.lin_feats is None:
+            assert lf is None and t.lin_feats is None
+        else:
+            np.testing.assert_array_equal(np.asarray(t.lin_feats), lf)
+    for b, b0 in zip(m.prep_.binner_.borders_, snap["borders"]):
+        np.testing.assert_array_equal(b, b0)
+    assert list(m.prep_.count_features_) == snap["count"]
+    assert list(m.cross_pairs) == snap["cross_pairs"]
+    assert list(m.valid_history_) == snap["valid"]
+    assert est.n_features_in_ == snap["n_features"]
+    if snap["names"] is None:
+        assert getattr(est, "feature_names_in_", None) is None
+    else:
+        np.testing.assert_array_equal(est.feature_names_in_, snap["names"])
+    assert getattr(est, "linear_leaves_selected_", None) == snap["ll_sel"]
+    assert est.cross_features_selected_ == snap["cx_sel"]
+    assert est.cross_pairs_ == snap["cx_pairs"]
+    if snap["temp"] is not None:
+        assert est.temperature_ == snap["temp"]
+    if snap["classes"] is not None:
+        np.testing.assert_array_equal(est.classes_, snap["classes"])
+
+
+@pytest.mark.parametrize("kind", ["regression", "binary"])
+def test_refresh_keeps_pinned_state(kind):
+    """(c) Nothing pinned moves; n_samples_trained_ grows by the new rows."""
+    d = _make_data()
+    if kind == "regression":
+        est = ChimeraBoostRegressor(random_state=0, n_estimators=60,
+                                    cross_features="always",
+                                    linear_leaves=True,
+                                    store_training_data=True)
+        y_fit, y_new = d["y_full"], d["y_new"]
+    else:
+        est = ChimeraBoostClassifier(random_state=0, n_estimators=60,
+                                     cross_features="always",
+                                     store_training_data=True)
+        y_fit = _binary_labels(d["y_full"])
+        y_new = _binary_labels(d["y_new"])
+    est.fit(d["X_full"], y_fit, cat_features=CAT, sample_weight=d["w_full"])
+    assert est.model_.prep_.count_features_ == [6]
+    est.shap_importances(d["X_full"][:20])
+    assert est._shap_importances_cache_ is not None
+    assert hasattr(est, "expected_value_")
+    snap = _pinned_snapshot(est)
+    n0 = est.n_samples_trained_
+    est.refresh(d["X_new"], y_new, sample_weight=d["w_new"])
+    _assert_pinned_equal(est, snap)
+    assert est.n_samples_trained_ == n0 + len(y_new)
+    assert est._shap_importances_cache_ is None
+    assert not hasattr(est, "expected_value_")
+
+
+def test_refresh_preserves_feature_names():
+    """(c) feature_names_in_ survives a refresh (DataFrame fit)."""
+    pd = pytest.importorskip("pandas")
+    d = _make_data()
+    num_cols = [f"n{i}" for i in range(N_NUM)]
+
+    def frame(X):
+        df = pd.DataFrame(np.asarray(X[:, :N_NUM], dtype=np.float64),
+                          columns=num_cols)
+        df["hi"] = list(X[:, 6])
+        df["lo"] = list(X[:, 7])
+        return df
+
+    cols = num_cols + ["hi", "lo"]
+    est = ChimeraBoostRegressor(random_state=0, n_estimators=20,
+                                store_training_data=True)
+    est.fit(frame(d["X_full"]), d["y_full"], cat_features=["hi", "lo"])
+    assert list(est.feature_names_in_) == cols
+    est.refresh(frame(d["X_new"]), d["y_new"])
+    assert list(est.feature_names_in_) == cols
+    assert est.n_samples_trained_ == len(d["y_full"]) + len(d["y_new"])
+
+
+def test_refresh_ignores_post_fit_set_params():
+    """refresh() replays the fitted configuration, not later set_params."""
+    d = _make_data()
+    est = ChimeraBoostRegressor(random_state=0, n_estimators=40,
+                                store_training_data=True)
+    est.fit(d["X_full"][:1500], d["y_full"][:1500], cat_features=CAT)
+    snap = _snapshot(est, d["X_held"])
+    est.set_params(n_estimators=5, learning_rate=0.5, max_bins=16, depth=2,
+                   subsample=0.5, l2_leaf_reg=99.0)
+    est.refresh(d["X_full"][:0], d["y_full"][:0])
+    _assert_snapshot_equal(est, snap, d["X_held"])
+
+
+def test_chained_refresh_equals_single():
+    """(c2) refresh(A) then refresh(B) equals one refresh(A+B), bit for bit."""
+    d = _make_data()
+
+    def fresh():
+        e = ChimeraBoostRegressor(random_state=0, n_estimators=50,
+                                  cross_features="always", linear_leaves=True,
+                                  store_training_data=True)
+        e.fit(d["X_full"], d["y_full"], cat_features=CAT,
+              sample_weight=d["w_full"])
+        return e
+
+    n_a = 120
+    e1 = fresh()
+    e1.refresh(d["X_new"][:n_a], d["y_new"][:n_a],
+               sample_weight=d["w_new"][:n_a])
+    e1.refresh(d["X_new"][n_a:], d["y_new"][n_a:],
+               sample_weight=d["w_new"][n_a:])
+    e2 = fresh()
+    e2.refresh(d["X_new"], d["y_new"], sample_weight=d["w_new"])
+    assert e1.n_samples_trained_ == e2.n_samples_trained_
+    _assert_booster_equal(e1.model_, e2.model_)
+    np.testing.assert_array_equal(e1.predict(d["X_held"]),
+                                  e2.predict(d["X_held"]))
+
+
+def test_pickle_round_trip():
+    """(d) A refreshed-after-round-trip model equals the refreshed original."""
+    d = _make_data()
+    est = ChimeraBoostRegressor(random_state=0, n_estimators=50,
+                                cross_features="always", linear_leaves=True,
+                                store_training_data=True)
+    est.fit(d["X_full"], d["y_full"], cat_features=CAT,
+            sample_weight=d["w_full"])
+    rt = pickle.loads(pickle.dumps(est))
+    assert rt._training_data_.n_rows == est._training_data_.n_rows
+    est.refresh(d["X_new"], d["y_new"], sample_weight=d["w_new"])
+    rt.refresh(d["X_new"], d["y_new"], sample_weight=d["w_new"])
+    _assert_booster_equal(rt.model_, est.model_)
+    np.testing.assert_array_equal(rt.predict(d["X_held"]),
+                                  est.predict(d["X_held"]))
+    assert rt.n_samples_trained_ == est.n_samples_trained_
+    assert est.model_.replay_donor is None
+    assert rt.model_.replay_donor is None
+
+    plain = ChimeraBoostRegressor(random_state=0, n_estimators=10).fit(
+        d["X_full"][:500], d["y_full"][:500], cat_features=CAT)
+    assert plain._training_data_ is None
+
+
+def test_store_training_data_errors():
+    """(e) Every slice-1 error message."""
+    d = _make_data()
+    X, y = d["X_full"][:800], d["y_full"][:800]
+
+    for bad in ("yes", 1, None):
+        with pytest.raises(ValueError, match="store_training_data must be"):
+            ChimeraBoostRegressor(random_state=0,
+                                  store_training_data=bad).fit(
+                X, y, cat_features=CAT)
+
+    with pytest.raises(NotImplementedError, match="n_ensembles > 1"):
+        ChimeraBoostRegressor(random_state=0, n_ensembles=2,
+                              store_training_data=True).fit(
+            X, y, cat_features=CAT)
+
+    with pytest.raises(NotImplementedError, match="n_ensembles > 1"):
+        ChimeraBoostRegressor(random_state=0, quality=4,
+                              store_training_data=True).fit(
+            X, y, cat_features=CAT)
+
+    y3 = np.where(y > np.quantile(y, 2 / 3), "c",
+                  np.where(y > np.quantile(y, 1 / 3), "b", "a"))
+    with pytest.raises(NotImplementedError, match="multiclass classification"):
+        ChimeraBoostClassifier(random_state=0,
+                               store_training_data=True).fit(
+            X, y3, cat_features=CAT)
+
+    with pytest.raises(NotImplementedError, match="loss='Quantile'"):
+        ChimeraBoostRegressor(random_state=0, loss="Quantile",
+                              store_training_data=True).fit(
+            X, y, cat_features=CAT)
+
+    groups = np.arange(len(y)) % 20
+    with pytest.raises(NotImplementedError, match="random_effects=True"):
+        ChimeraBoostRegressor(random_state=0, random_effects=True,
+                              store_training_data=True).fit(
+            X, y, cat_features=CAT, groups=groups)
+
+    with pytest.raises(NotFittedError):
+        ChimeraBoostRegressor().refresh(X[:10], y[:10])
+
+    plain = ChimeraBoostRegressor(random_state=0, n_estimators=10).fit(
+        X, y, cat_features=CAT)
+    with pytest.raises(ValueError, match="store_training_data=False"):
+        plain.refresh(X[:10], y[:10])
+
+
+def test_refresh_weight_and_label_errors():
+    """(e) Weighted/unweighted mismatches and unseen refresh labels."""
+    d = _make_data()
+    X, y = d["X_full"][:800], d["y_full"][:800]
+
+    w = ChimeraBoostRegressor(random_state=0, n_estimators=10,
+                              store_training_data=True).fit(
+        X, y, cat_features=CAT, sample_weight=d["w_full"][:800])
+    with pytest.raises(ValueError, match="needs sample_weight"):
+        w.refresh(X[:5], y[:5])
+
+    u = ChimeraBoostRegressor(random_state=0, n_estimators=10,
+                              store_training_data=True).fit(
+        X, y, cat_features=CAT)
+    with pytest.raises(ValueError, match="cannot take sample_weight"):
+        u.refresh(X[:5], y[:5], sample_weight=np.ones(5))
+
+    yb = _binary_labels(y)
+    clf = ChimeraBoostClassifier(random_state=0, n_estimators=10,
+                                 store_training_data=True).fit(
+        X, yb, cat_features=CAT)
+    bad_y = np.array(["pos", "nope", "neg"], dtype=object)
+    with pytest.raises(ValueError, match="cannot add classes"):
+        clf.refresh(X[:3], bad_y)
+
+
+@pytest.mark.parametrize("cls", [ChimeraBoostRegressor, ChimeraBoostClassifier])
+def test_store_param_round_trip(cls):
+    """(e) clone and get_params round-trip the new parameter."""
+    est = cls(store_training_data=True)
+    assert est.get_params()["store_training_data"] is True
+    assert clone(est).store_training_data is True
+    est.set_params(store_training_data=False)
+    assert est.get_params()["store_training_data"] is False
+    assert cls().get_params()["store_training_data"] is False
+
+
+@pytest.mark.parametrize("kind", ["regression", "binary"])
+def test_default_fit_unaffected(kind):
+    """(e) The default fit is unaffected: flag False vs True agree exactly."""
+    d = _make_data()
+    if kind == "regression":
+        cls, y = ChimeraBoostRegressor, d["y_full"]
+    else:
+        cls, y = ChimeraBoostClassifier, _binary_labels(d["y_full"])
+    kw = dict(random_state=0, n_estimators=50)
+    a = cls(store_training_data=False, **kw).fit(
+        d["X_full"], y, cat_features=CAT, sample_weight=d["w_full"])
+    b = cls(store_training_data=True, **kw).fit(
+        d["X_full"], y, cat_features=CAT, sample_weight=d["w_full"])
+    np.testing.assert_array_equal(a.predict(d["X_held"]), b.predict(d["X_held"]))
+    if kind == "binary":
+        np.testing.assert_array_equal(a.predict_proba(d["X_held"]),
+                                      b.predict_proba(d["X_held"]))
+    assert a._training_data_ is None
+    assert b._training_data_ is not None
+    assert a.n_samples_trained_ == b.n_samples_trained_
+
+
+def test_refresh_usefulness_smoke():
+    """(f) Usefulness smoke (not a gate): refresh lands between 60% and 90%.
+
+    Fit on the first 60% of rows, refresh with the next 30%, and compare
+    test RMSE against the 60% model and a full refit on 90%. Refresh is
+    expected between the two, but that is not asserted.
+    """
+    d = _make_data()
+    X60, y60 = d["X_full"][:1560], d["y_full"][:1560]
+    X30, y30 = d["X_full"][1560:2340], d["y_full"][1560:2340]
+
+    def rmse(est):
+        return float(np.sqrt(np.mean((est.predict(d["X_held"]) - d["y_held"]) ** 2)))
+
+    base = ChimeraBoostRegressor(random_state=0, n_estimators=200,
+                                 store_training_data=True).fit(
+        X60, y60, cat_features=CAT)
+    rmse_base = rmse(base)
+    base.refresh(X30, y30)
+    rmse_ref = rmse(base)
+    full = ChimeraBoostRegressor(random_state=0, n_estimators=200).fit(
+        d["X_full"][:2340], d["y_full"][:2340], cat_features=CAT)
+    rmse_full = rmse(full)
+    print(f"(f) RMSE: 60%={rmse_base:.4f} refreshed={rmse_ref:.4f} "
+          f"90%={rmse_full:.4f}")
+    assert np.isfinite([rmse_base, rmse_ref, rmse_full]).all()
diff --git a/tests/test_replay_kernel.py b/tests/test_replay_kernel.py
new file mode 100644
index 0000000..ade237a
--- /dev/null
+++ b/tests/test_replay_kernel.py
@@ -0,0 +1,210 @@
+"""Bit-identity guards for the fast replay kernel (issue131, I066).
+
+``_linear_leaf_fit_ref`` is a frozen copy of ``tree._linear_leaf_fit`` as of
+2026-09-24, before the parallel rewrite. Every test asserts exact
+``np.array_equal`` against it: any summation-order or expression change fails
+in the last ulp.
+"""
+import numba
+import numpy as np
+import pytest
+from numba import njit, prange
+
+from chimeraboost.tree import (
+    ObliviousTree,
+    _assign_leaves,
+    _leaf_values,
+    _linear_leaf_fit,
+    _solve_small,
+    replay_oblivious_tree,
+)
+
+DEFAULT_THREADS = numba.get_num_threads()
+
+
+@njit(cache=False, parallel=True)
+def _linear_leaf_fit_ref(leaf, grad, hess, n_leaves, lin_feats, centers_std, Xb,  # noqa: C901 -- frozen reference kernel
+                         l2_intercept, lin_lambda, lr):
+    """Frozen copy of tree._linear_leaf_fit (2026-09-24)."""
+    n = leaf.shape[0]
+    k = lin_feats.shape[0]
+    d = 1 + k
+    coef = np.zeros((n_leaves, d))
+
+    counts = np.zeros(n_leaves, dtype=np.int64)
+    Gtot = np.zeros(n_leaves)
+    Htot = np.zeros(n_leaves)
+    for i in range(n):
+        l = leaf[i]
+        counts[l] += 1
+        Gtot[l] += grad[i]
+        Htot[l] += hess[i]
+
+    start = np.zeros(n_leaves + 1, dtype=np.int64)
+    for l in range(n_leaves):
+        start[l + 1] = start[l] + counts[l]
+
+    pos = start[:n_leaves].copy()
+    order = np.empty(n, dtype=np.int64)
+    for i in range(n):
+        l = leaf[i]
+        order[pos[l]] = i
+        pos[l] += 1
+
+    for l in prange(n_leaves):
+        if counts[l] == 0:
+            continue
+
+        if counts[l] < 2 * d or k == 0:
+            if Htot[l] > 0.0:
+                coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
+            continue
+
+        Ml = np.zeros((d, d))
+        rl = np.zeros(d)
+        xrow = np.empty(k)
+
+        for q in range(start[l], start[l + 1]):
+            i = order[q]
+            h = hess[i]
+            g = grad[i]
+
+            for j in range(k):
+                f = lin_feats[j]
+                v = centers_std[f, Xb[f, i]]
+                xrow[j] = v if np.isfinite(v) else 0.0
+
+            Ml[0, 0] += h
+            rl[0] += -g
+            for j in range(k):
+                xj = xrow[j]
+                Ml[0, 1 + j] += h * xj
+                Ml[1 + j, 0] += h * xj
+                rl[1 + j] += -g * xj
+                for jj in range(k):
+                    Ml[1 + j, 1 + jj] += h * xj * xrow[jj]
+
+        Ml[0, 0] += l2_intercept
+        for j in range(1, d):
+            Ml[j, j] += lin_lambda
+        for j in range(d):
+            Ml[j, j] += 1e-9
+
+        beta = _solve_small(Ml, rl)
+        if np.isnan(beta[0]):
+            if Htot[l] > 0.0:
+                coef[l, 0] = -lr * Gtot[l] / (Htot[l] + l2_intercept)
+            continue
+
+        for j in range(d):
+            coef[l, j] = lr * beta[j]
+
+    return coef
+
+
+def _make_inputs(n, n_leaves, k, seed):
+    """Random kernel inputs covering empty, tiny (<2d) and 60%+ leaves."""
+    rng = np.random.default_rng(seed)
+    d = 1 + k
+    n_features, max_bins = 8, 32
+
+    leaf = np.empty(n, dtype=np.int64)
+    n_big = max(int(0.6 * n), 2 * d)
+    if n_big + 5 <= n:
+        n_s1, n_s2 = 2, 3
+    elif n_big + 2 <= n:
+        n_s1, n_s2 = 1, 1
+    else:
+        n_big = n - 2
+        n_s1, n_s2 = 1, 1
+    leaf[:n_big] = 0
+    leaf[n_big:n_big + n_s1] = 1
+    leaf[n_big + n_s1:n_big + n_s1 + n_s2] = 2
+    rest = n - (n_big + n_s1 + n_s2)
+    if rest > 0:
+        leaf[n_big + n_s1 + n_s2:] = rng.integers(3, n_leaves - 2, rest)
+
+    grad = rng.standard_normal(n)
+    hess = rng.random(n) + 0.1
+    hess[rng.random(n) < 0.05] = 0.0
+    lin_feats = np.arange(k, dtype=np.int64)
+    centers_std = rng.standard_normal((n_features, max_bins))
+    centers_std[rng.random((n_features, max_bins)) < 0.05] = np.nan
+    Xb = rng.integers(0, max_bins, size=(n_features, n)).astype(np.uint16)
+    return leaf, grad, hess, lin_feats, centers_std, Xb
+
+
+@pytest.mark.parametrize("n,n_leaves", [(16, 8), (2000, 16), (200_000, 64)])
+@pytest.mark.parametrize("k", [1, 6])
+@pytest.mark.parametrize("n_threads", [1, DEFAULT_THREADS])
+def test_linear_leaf_fit_matches_ref(n, n_leaves, k, n_threads):
+    leaf, grad, hess, lin_feats, cs, Xb = _make_inputs(n, n_leaves, k, 1000 + n + k)
+    prev = numba.get_num_threads()
+    numba.set_num_threads(n_threads)
+    try:
+        ref = _linear_leaf_fit_ref(leaf, grad, hess, n_leaves, lin_feats, cs,
+                                   Xb, 1.0, 1.0, 0.1)
+        got = _linear_leaf_fit(leaf, grad, hess, n_leaves, lin_feats, cs,
+                               Xb, 1.0, 1.0, 0.1)
+    finally:
+        numba.set_num_threads(prev)
+    assert got.shape == ref.shape
+    assert np.array_equal(got, ref)
+
+
+def _make_donor(n_features, max_bins, depth, k, seed):
+    rng = np.random.default_rng(seed)
+    sf = rng.integers(0, n_features, depth).astype(np.int64)
+    st = rng.integers(0, max_bins - 1, depth).astype(np.int64)
+    n_leaves = 1 << depth
+    vals = rng.standard_normal(n_leaves)
+    lin_feats = np.arange(k, dtype=np.int64)
+    lin_coef = np.zeros((n_leaves, 1 + k))
+    return ObliviousTree(sf, st, vals, np.zeros(depth),
+                         lin_feats=lin_feats, lin_coef=lin_coef,
+                         centers_std=None), n_leaves
+
+
+@pytest.mark.parametrize("n,depth", [(2000, 4), (200_000, 6)])
+@pytest.mark.parametrize("k", [1, 6])
+@pytest.mark.parametrize("n_threads", [1, DEFAULT_THREADS])
+def test_replay_matches_ref(n, depth, k, n_threads):
+    rng = np.random.default_rng(77 + n + k)
+    n_features, max_bins = 8, 32
+    donor, n_leaves = _make_donor(n_features, max_bins, depth, k, 5)
+    Xb = rng.integers(0, max_bins, size=(n_features, n)).astype(np.uint16)
+    grad = rng.standard_normal(n)
+    hess = rng.random(n) + 0.1
+    cs = rng.standard_normal((n_features, max_bins))
+    cs[rng.random((n_features, max_bins)) < 0.05] = np.nan
+
+    prev = numba.get_num_threads()
+    numba.set_num_threads(n_threads)
+    try:
+        tree, leaf = replay_oblivious_tree(donor, Xb, grad, hess, 1.0, 0.1,
+                                           linear_leaves=True, centers_std=cs,
+                                           linear_lambda=1.0)
+        leaf_ref = _assign_leaves(Xb, donor.splits_feat, donor.splits_thr)
+        values_ref = _leaf_values(leaf_ref, grad, hess, n_leaves, 1.0, 0.1)
+        coef_ref = _linear_leaf_fit_ref(leaf_ref, grad, hess, n_leaves,
+                                        donor.lin_feats, cs, Xb, 1.0, 1.0, 0.1)
+    finally:
+        numba.set_num_threads(prev)
+    assert np.array_equal(leaf, leaf_ref)
+    assert np.array_equal(tree.values, values_ref)
+    assert np.array_equal(tree.lin_coef, coef_ref)
+
+
+def test_replay_constant_matches_ref():
+    """Linear off: replay still matches (covers the assign swap alone)."""
+    rng = np.random.default_rng(9)
+    n, n_features, max_bins = 200_000, 8, 32
+    donor, n_leaves = _make_donor(n_features, max_bins, 6, 2, 6)
+    Xb = rng.integers(0, max_bins, size=(n_features, n)).astype(np.uint16)
+    grad = rng.standard_normal(n)
+    hess = rng.random(n) + 0.1
+    tree, leaf = replay_oblivious_tree(donor, Xb, grad, hess, 1.0, 0.1)
+    leaf_ref = _assign_leaves(Xb, donor.splits_feat, donor.splits_thr)
+    values_ref = _leaf_values(leaf_ref, grad, hess, n_leaves, 1.0, 0.1)
+    assert np.array_equal(leaf, leaf_ref)
+    assert np.array_equal(tree.values, values_ref)