Reproducible Pipelines

A reproducible pipeline preserves the information needed to rerun data preparation, fitting, and evaluation. Exact output bytes, agreement within a numerical tolerance, stable decisions, and replication of a statistical conclusion are different targets. Choose the target before deciding which differences count as a failure.

Floating point is not associative

import numpy as np

g = np.random.default_rng(0)
x = g.normal(size=1_000_000).astype(np.float32)
orders = {"as generated": x, "reversed": x[::-1],
          "sorted": np.sort(x), "shuffled": g.permutation(x)}
sums = {k: float(np.add.reduce(v, dtype=np.float32)) for k, v in orders.items()}
for k, v in sums.items():
    print(f"{k:14s} {v!r}")
print(f"spread between orderings: {max(sums.values()) - min(sums.values()):.6e}")
print(f"float64 reference:        {np.add.reduce(x, dtype=np.float64):.10f}")
# as generated   998.5707397460938
# reversed       998.5708618164062
# sorted         998.5911865234375
# shuffled       998.5706176757812
# spread between orderings: 2.056885e-02
# float64 reference:        998.5706627697

The four float32 reductions differ in this run. Float32 represents a limited set of numbers, so intermediate additions round; changing their grouping can change the result. For a smaller float32 example, \((10^8-10^8)+1=1\) versus \(10^8+(-10^8+1)=0\). NumPy can use partial pairwise summation, so the sorted-array result should not be explained as if it necessarily used one left-to-right running total. The float64 reduction is a higher-precision reference for these stored float32 inputs, not an exact sum of the original unrounded draws.

Changing reduction order can change the result, but need not do so. Threading, chunking, row order, library kernels, and hardware can affect numerical operations. A random seed controls a random-number stream; it does not fix all of these choices.

import numpy as np
from threadpoolctl import threadpool_info, threadpool_limits

g = np.random.default_rng(0)
A = g.normal(size=(1000, 1000)); B = g.normal(size=(1000, 1000))
reference = None
for requested in (1, 2, 4):
    with threadpool_limits(limits=requested, user_api="blas"):
        pools = [p for p in threadpool_info() if p["user_api"] == "blas"]
        if not pools:
            raise RuntimeError("No supported BLAS pool found; thread control was not verified")
        result = A @ B
        if reference is None:
            reference = result.copy()
        print("requested", requested, "reported", [p["num_threads"] for p in pools])
        print("same values", np.array_equal(reference, result),
              "max difference", f"{np.max(np.abs(reference-result)):.3e}")
# requested 1 reported [1]
# same values True max difference 0.000e+00
# requested 2 reported [2]
# same values True max difference 0.000e+00
# requested 4 reported [4]
# same values False max difference 1.066e-13

This block actually changes the supported BLAS pool limit and reports it; the earlier summation block does not use BLAS matrix multiplication. A reported limit is not proof that every operation used that many workers. Results can be identical across limits on some builds. OMP_NUM_THREADS is not a universal BLAS control: OpenBLAS, MKL, and other backends have their own controls, often read at initialization. Record the active runtime as well as requested settings. In this run, one and two threads agree, while four differ by at most about 1.07e-13. This experiment’s numerical outputs are specific to the installed backend; another build can produce a different comparison.

Choose tolerances in terms of the quantity and decision that matter. A tiny change can flip a threshold decision when a score is sufficiently close to its cutoff; bitwise agreement can also be useful for debugging. Parallel execution is not inherently nondeterministic, and serial execution alone does not ensure agreement across platforms. The input-perturbation exercise is a local sensitivity check, not a test of all effects of changing BLAS threads during training.

Reproducible rankings need an explicit tie-break

Even unchanged scores can produce a different selected set if ties have no stable secondary key. This example scores the training rows deliberately to create many ties; it is not a generalization evaluation.

from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier

X, y = make_classification(n_samples=40000, n_features=20, n_informative=8,
                           class_sep=0.6, random_state=0)
p = RandomForestClassifier(200, random_state=0).fit(X, y).predict_proba(X)[:, 1]
print(f"distinct probability values observed: {len(np.unique(p))}")
print(f"rows sharing the maximum value {p.max():.4f}: {np.sum(p == p.max())}")

print(f"{'k':>7} {'top-k overlap after reordering the rows':>40}")
for k in (10, 100, 1000, 5000):
    base = set(np.argsort(-p, kind="stable")[:k])
    ov = []
    for rep in range(20):
        perm = np.random.default_rng(rep).permutation(len(p))
        top = perm[np.argsort(-p[perm], kind="stable")[:k]]         # same rows, different arrival order
        ov.append(len(base & set(top)) / k)
    print(f"{k:7d} {np.mean(ov):40.4f}")
# distinct probability values observed: 172
# rows sharing the maximum value 1.0000: 1363
#       k  top-k overlap after reordering the rows
#      10                                   0.0050
#     100                                   0.0725
#    1000                                   0.7353
#    5000                                   0.9341

The distinct-value count describes these predictions, not the full output range of every 200-tree forest. A forest averages tree leaf class proportions; impure leaves can produce more than a simple grid of vote fractions. Here 1,363 rows share the maximum score, and the mean top-10 overlap after reordering is only 0.0050. The stable sort used here deliberately preserves arrival order within each tie. NumPy’s default argsort is not guaranteed to do that. Overlap is the mean fraction of selected IDs retained over 20 permutations, not a probability estimate for every possible deployment ordering.

Sort by descending score and then an immutable, unique record ID. The IDs below represent identities assigned before reordering; regenerating them from each incoming row position would defeat the fix. Define policies for duplicate IDs and nonfinite scores. A deterministic tie-break makes a decision repeatable, but does not by itself make the tie policy appropriate.

ids = np.arange(len(p))
assert np.unique(ids).size == ids.size and np.isfinite(p).all()
for k in (10, 100, 1000, 5000):
    baseline = ids[np.lexsort((ids, -p))[:k]]
    overlaps = []
    for rep in range(20):
        perm = np.random.default_rng(rep).permutation(len(p))
        selected = ids[perm][np.lexsort((ids[perm], -p[perm]))[:k]]
        assert np.array_equal(baseline, selected)
        overlaps.append(len(set(baseline) & set(selected)) / k)
    print(k, f"stable-ID overlap {np.mean(overlaps):.4f}")
# 10 stable-ID overlap 1.0000
# 100 stable-ID overlap 1.0000
# 1000 stable-ID overlap 1.0000
# 5000 stable-ID overlap 1.0000

Save the preprocessing and the prediction rule together

A pipeline chains operations in their execution order. Here missing values are filled with training medians, features are standardized using training statistics, and a logistic model is fitted. Calling predict_proba reuses those fitted transformations. Split the raw rows first; during cross-validation, fit the entire pipeline within each training fold. A pipeline cannot repair a split that already leaks people, time, or future information.

import json, hashlib, tempfile
from pathlib import Path
import joblib
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification

Xrun, yrun = make_classification(n_samples=1000, n_features=10, random_state=17)
Xrun[::17, 0] = np.nan
train_ids, test_ids = train_test_split(np.arange(len(yrun)), test_size=.3,
                                      stratify=yrun, random_state=23)
config = {"data_seed": 17, "split_seed": 23, "C": 1.0, "solver": "lbfgs", "max_iter": 1000}
pipe = make_pipeline(SimpleImputer(strategy="median"), StandardScaler(),
                     LogisticRegression(C=config["C"], solver=config["solver"], max_iter=config["max_iter"]))
pipe.fit(Xrun[train_ids], yrun[train_ids])
expected = pipe.predict_proba(Xrun[test_ids])
with tempfile.TemporaryDirectory() as directory:
    root = Path(directory)
    np.savez(root / "data.npz", X=Xrun, y=yrun, train_ids=train_ids, test_ids=test_ids)
    joblib.dump(pipe, root / "pipeline.joblib")
    (root / "config.json").write_text(json.dumps(config, sort_keys=True))
    names = ["data.npz", "pipeline.joblib", "config.json"]
    hashes = {name: hashlib.sha256((root/name).read_bytes()).hexdigest() for name in names}
    assert all(hashlib.sha256((root/name).read_bytes()).hexdigest() == h for name,h in hashes.items())
    with np.load(root / "data.npz") as saved:
        restored = joblib.load(root / "pipeline.joblib")
        actual = restored.predict_proba(saved["X"][saved["test_ids"]])
    print("restored probabilities match", np.array_equal(expected, actual))
    print("saved artifacts", sorted(hashes))
print("training-only median", np.isclose(pipe[0].statistics_[0], np.nanmedian(Xrun[train_ids, 0])))
# restored probabilities match True
# saved artifacts ['config.json', 'data.npz', 'pipeline.joblib']
# training-only median True

The snapshot stores the row order and exact split IDs, while the saved pipeline includes the fitted imputer, scaler, and classifier. This example restores predictions in the same process and environment; it does not prove independent retraining or cross-version compatibility. The temporary directory is deleted when the block finishes. Use a retained, versioned run directory in a real project. Load joblib artifacts only from trusted sources, and retain the compatible environment needed to read them.

A handoff also needs the runnable entry point, source revision and any uncommitted changes, feature names and order, resolved training configuration, metric and threshold definitions, data provenance and snapshot location, and an environment specification. A hash can verify an available artifact but cannot recover a missing one. Seeds alone do not preserve a split if the input rows change order; stable record IDs or the saved split are needed. For a fresh-process check, load the bundle without relying on notebook state, then compare agreed outputs.

Exercises

1. How far does a 1e-12 difference travel? Hold fitted models fixed, perturb their inputs, and compare probabilities and decisions.

You should get: small probability changes with no flips on this sample, plus a boundary example that does flip.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression

X, y = make_classification(n_samples=40000, n_features=20, n_informative=8,
                           class_sep=0.6, random_state=0)
g = np.random.default_rng(0)
Xp = X * (1 + 1e-12 * g.normal(size=X.shape))        # a chosen relative perturbation

for name, m in (("RandomForest", RandomForestClassifier(200, random_state=0)),
                ("LogisticRegression", LogisticRegression(max_iter=3000))):
    m.fit(X, y)
    p1, p2 = m.predict_proba(X)[:, 1], m.predict_proba(Xp)[:, 1]
    print(f"{name:20s} max |dp| {np.max(np.abs(p1 - p2)):.3e}"
          f"   label flips {np.sum((p1 >= 0.5) != (p2 >= 0.5))}/{len(y)}")
from sklearn.tree import DecisionTreeClassifier
stump = DecisionTreeClassifier(max_depth=1, random_state=0).fit([[0.0], [2.0]], [0, 1])
boundary = np.float32(stump.tree_.threshold[0])
nearby = np.array([boundary, np.nextafter(boundary, np.float32(np.inf))], dtype=np.float32)
print("adjacent float32 boundary predictions", stump.predict(nearby[:, None]))
# RandomForest         max |dp| 0.000e+00   label flips 0/40000
# LogisticRegression   max |dp| 3.471e-12   label flips 0/40000
# adjacent float32 boundary predictions [0 1]

2. Which seed is actually moving the number? Compare model-seed variability within each split with variation across split averages.

You should get: a crossed table that shows whether seed sensitivity differs across splits.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score

Xseed, yseed = make_classification(n_samples=4000, n_features=20, n_informative=8,
                                  class_sep=.5, flip_y=.05, random_state=0)
scores = np.empty((4, 4))
for split_seed in range(4):
    Xtr, Xte, ytr, yte = train_test_split(Xseed, yseed, test_size=.3,
                                       random_state=split_seed, stratify=yseed)
    for model_seed in range(4):
        model = RandomForestClassifier(80, random_state=model_seed).fit(Xtr, ytr)
        scores[split_seed, model_seed] = roc_auc_score(yte, model.predict_proba(Xte)[:, 1])
print("rows=split seed, columns=model seed")
print(np.round(scores, 4))
print("within-split model SD", np.round(scores.std(axis=1, ddof=1), 4))
print("split means", np.round(scores.mean(axis=1), 4))
print("SD of split means", f"{scores.mean(axis=1).std(ddof=1):.4f}")
# rows=split seed, columns=model seed
# [[0.9161 0.9132 0.9125 0.9178]
#  [0.9161 0.921  0.9168 0.9193]
#  [0.9246 0.926  0.9255 0.9241]
#  [0.9029 0.9023 0.9029 0.9039]]
# within-split model SD [0.0025 0.0023 0.0008 0.0007]
# split means [0.9149 0.9183 0.9251 0.903 ]
# SD of split means 0.0092

3. Write down what an environment is. Capture runtime details and array fingerprints, then identify what the record still needs.

You should get: a manifest, and a clear view of which entries earn their place.

import sys, os, platform, json, hashlib, io, contextlib
import numpy as np, sklearn, scipy
from threadpoolctl import threadpool_info

def array_record(value):
    arr = np.asarray(value)
    if arr.dtype.kind not in "biufc":
        raise TypeError("This example hashes numeric arrays only")
    header = json.dumps({"shape": list(arr.shape), "dtype": arr.dtype.str}, sort_keys=True).encode()
    digest = hashlib.sha256(len(header).to_bytes(8, "big") + header + np.ascontiguousarray(arr).tobytes()).hexdigest()
    return {"shape": list(arr.shape), "dtype": arr.dtype.str, "sha256": digest}

def manifest(X, y):
    output = io.StringIO()
    with contextlib.redirect_stdout(output):
        np.show_config()
    return {"python": sys.version, "platform": platform.platform(),
            "machine": platform.machine(), "numpy": np.__version__,
            "scipy": scipy.__version__, "sklearn": sklearn.__version__,
            "threads_requested": {k: os.environ.get(k) for k in
                ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS")},
            "threadpools_reported": threadpool_info(), "numpy_build": output.getvalue(),
            "X": array_record(X), "y": array_record(y)}

m = manifest(Xrun, yrun)
serialized = json.dumps(m, indent=2)
print("manifest JSON roundtrip", json.loads(serialized) == m)
x0 = np.arange(6, dtype=np.float64)
print("shape changes hash", array_record(x0)["sha256"] != array_record(x0.reshape(2,3))["sha256"])
print("dtype changes hash", array_record(x0)["sha256"] != array_record(x0.view(np.int64))["sha256"])
print("copy preserves hash", array_record(x0)["sha256"] == array_record(x0.copy())["sha256"])
# manifest JSON roundtrip True
# shape changes hash True
# dtype changes hash True
# copy preserves hash True

References

NumPy summation; Scikit-learn: Common pitfalls and recommended practices; Scikit-learn model persistence.