Upload provenance_reconstruct.py

provenance_reconstruct.py

@@ -0,0 +1,846 @@
+#!/usr/bin/env python3
+"""
+provenance_reconstruct.py
+=========================
+
+Forensic provenance reconstruction for bare numbers.
+
+Given a value that arrived without provenance — a float from a solver,
+a token count from an API — infer the most likely provenance that
+produced it, with calibrated confidence.
+
+This is the consumer-side complement to numerical_provenance.py and
+token_provenance.py.  Those tools produce values that carry their
+provenance.  This tool reconstructs provenance for values that lost
+it at an API boundary, a log line, or a file copy.
+
+Pure stdlib + numpy.  No downloads.  No tokenizer files.  No model hub.
+The "tokenizer" is a 12-line local function.  The "model" is a small
+MLP trained in seconds on data generated by the tools themselves.
+
+The output is deliberately marked:
+
+    observed provenance   [ok: jacobi_converged -> participation_ratio]
+    inferred provenance   [CONVERGED? p=0.87]     <-- this tool
+    OOD reconstruction    [CONVERGED?? OOD p=1.00]
+
+A downstream consumer can require observed provenance and refuse
+inferred provenance.  The '?' is the honesty marker.  The '??' plus
+'OOD' means the input is outside the training distribution and the
+confidence is not meaningful.
+
+Findings this file demonstrates
+-------------------------------
+
+The token task is a clean positive result: the residual against a
+free local estimator equals the template offset exactly, so a tiny
+model learns the three-way classification at ~100% accuracy.  This is
+arithmetic dressed up as learning; the model is packaging, not magic.
+
+The numeric task is a clean negative result: once the leaking solver-
+config features are removed, a bare value + matrix size carries very
+little recoverable convergence information.  Accuracy lands near the
+majority-class baseline.  This is a finding, not a failure.  The value
+of the tool is its *calibration* and its *OOD marker* — it says
+"[? p=0.61]" on the inputs where reconstruction is genuinely
+uncertain, and "[?? OOD p=1.00]" on the inputs where the confidence
+is not meaningful at all.
+
+Runnable.  Requires only numpy.
+"""
+
+from __future__ import annotations
+
+import math
+import random
+import re
+from dataclasses import dataclass
+from typing import List, Sequence, Tuple
+
+import numpy as np
+
+
+# =====================================================================
+# §1  Utilities
+# =====================================================================
+
+def decimal_digits(x: float) -> int:
+    """Significant decimal digits visible in repr(x)."""
+    s = f"{x:.15g}"
+    if "e" in s or "E" in s:
+        mantissa = s.split("e")[0].split("E")[0]
+    else:
+        mantissa = s
+    if "." in mantissa:
+        frac = mantissa.split(".")[1].rstrip("0")
+        return len(frac)
+    return 0
+
+
+def punct_ratio(text: str) -> float:
+    if not text:
+        return 0.0
+    n = sum(1 for ch in text if not ch.isalnum() and not ch.isspace())
+    return n / len(text)
+
+
+def digit_ratio(text: str) -> float:
+    if not text:
+        return 0.0
+    return sum(1 for ch in text if ch.isdigit()) / len(text)
+
+
+# =====================================================================
+# §2  Inlined numerical primitives
+#     (in the real repo these are imports from numerical_provenance)
+# =====================================================================
+
+def jacobi_eigen(A: List[List[float]], max_iter: int,
+                 tol: float = 1e-10) -> Tuple[List[float], dict]:
+    n = len(A)
+    B = [row[:] for row in A]
+    converged = False
+    iters = 0
+
+    for it in range(max_iter):
+        iters = it + 1
+        off = 0.0
+        p, q = 0, 1
+        for i in range(n):
+            for j in range(i + 1, n):
+                if abs(B[i][j]) > off:
+                    off = abs(B[i][j])
+                    p, q = i, j
+        if off < tol:
+            converged = True
+            break
+        theta = 0.5 * math.atan2(2 * B[p][q], B[p][p] - B[q][q])
+        c, s = math.cos(theta), math.sin(theta)
+        for i in range(n):
+            if i == p or i == q:
+                continue
+            bip, biq = B[i][p], B[i][q]
+            B[i][p] = c * bip - s * biq
+            B[p][i] = B[i][p]
+            B[i][q] = s * bip + c * biq
+            B[q][i] = B[i][q]
+        bpp, bqq, bpq = B[p][p], B[q][q], B[p][q]
+        B[p][p] = c * c * bpp - 2 * c * s * bpq + s * s * bqq
+        B[q][q] = s * s * bpp + 2 * c * s * bpq + c * c * bqq
+        B[p][q] = B[q][p] = 0.0
+
+    off_sq = 0.0
+    for i in range(n):
+        for j in range(i + 1, n):
+            off_sq += 2 * B[i][j] ** 2
+    residual = math.sqrt(off_sq)
+
+    eigenvalues = sorted([B[i][i] for i in range(n)], reverse=True)
+    return eigenvalues, {
+        "converged": converged,
+        "iterations": iters,
+        "residual": residual,
+    }
+
+
+def participation_ratio(eigs: Sequence[float]) -> float:
+    pos = [e for e in eigs if e > 1e-12]
+    if not pos:
+        return 0.0
+    s = sum(pos)
+    ss = sum(e * e for e in pos)
+    return (s * s) / ss if ss > 0 else 0.0
+
+
+def build_symmetric(n: int, rng: random.Random) -> List[List[float]]:
+    A = [[0.0] * n for _ in range(n)]
+    for i in range(n):
+        for j in range(i, n):
+            v = rng.gauss(0, 1)
+            A[i][j] = v
+            A[j][i] = v
+    return A
+
+
+# =====================================================================
+# §3  Inlined token primitives
+#     (in the real repo these are imports from token_provenance)
+# =====================================================================
+
+_TOKEN_RE = re.compile(r"[A-Za-z]+|\d+|[^\sA-Za-z\d]")
+
+
+def bpe_like(text: str) -> int:
+    """Local, free, no-download token estimator.
+
+    This is the "tokenizer" the reconstructor uses.  Alpha runs split
+    into ~4-char chunks, digits into ~2-char chunks, punctuation
+    counts as one.  It is not exact — but it is the *reference* the
+    residual features subtract from.
+    """
+    n = 0
+    for tok in _TOKEN_RE.findall(text):
+        if tok.isalpha():
+            n += max(1, (len(tok) + 3) // 4)
+        elif tok.isdigit():
+            n += max(1, (len(tok) + 1) // 2)
+        else:
+            n += 1
+    return n
+
+
+# =====================================================================
+# §4  Feature extraction
+# =====================================================================
+
+NUM_FEAT_NAMES = (
+    "log10|value|",
+    "decimal_digits",
+    "n",
+    "value / (n+1)",
+    "value_round",
+    "value_small",
+    "value_in_range_n",
+    "value_cluster",
+)
+
+TOK_FEAT_NAMES = (
+    "count",
+    "count / chars",
+    "count / words",
+    "count - bpe_like",
+    "count - bpe_like - 11",
+    "count - bpe_like - 7",
+    "punct_ratio",
+    "digit_ratio",
+    "chars",
+    "words",
+    "bpe_like",
+    "count % 8",
+)
+
+
+def numeric_features(value: float, n: int) -> np.ndarray:
+    """Observables about a bare numerical output.
+
+    Deliberately excludes max_iter and tol: those are solver config,
+    not observable at reconstruction time in the honest setting.  The
+    convergence state must be inferred from the *shape* of the value
+    and the matrix size, not from the request that produced it.
+
+    If a feature that encodes max_iter or tol is added here, accuracy
+    jumps to ~99% and the model is cheating.  The self-test checks for
+    that.
+    """
+    log_v = math.log10(abs(value) + 1e-30)
+    digits = decimal_digits(value)
+
+    # PR <= n always.  A value outside [1, n] is suspicious.
+    in_range = 1.0 if 1.0 <= abs(value) <= float(n) + 1e-9 else 0.0
+
+    # Converged participation ratios cluster near divisors of n
+    # (eigenvalue spectra of random symmetric matrices have effective
+    # rank roughly proportional to n).  Distance to the nearest
+    # expected cluster, normalized by n.
+    clusters = (n, n / 2, n / 4, n / 8)
+    cluster_dist = min(abs(value - c) for c in clusters) / max(1.0, n)
+
+    return np.array([
+        log_v,
+        digits,
+        float(n),
+        value / (n + 1.0),
+        1.0 if abs(value - round(value)) < 1e-9 else 0.0,
+        1.0 if abs(value) < 1.0 else 0.0,
+        in_range,
+        cluster_dist,
+    ], dtype=np.float64)
+
+
+def token_features(count: int, text: str) -> np.ndarray:
+    """Observables about a bare token count.
+
+    The residual features are exact arithmetic.  The template offset
+    is a fixed integer, so the model learns a threshold lookup and
+    reports 100% accuracy.  This is not a bug — it is a finding.  The
+    token side is arithmetic dressed as learning.
+    """
+    chars = len(text)
+    words = len(text.split())
+    local = bpe_like(text)
+    return np.array([
+        float(count),
+        count / max(1, chars),
+        count / max(1, words),
+        float(count - local),
+        float(count - local - 11),
+        float(count - local - 7),
+        punct_ratio(text),
+        digit_ratio(text),
+        float(chars),
+        float(words),
+        float(local),
+        float(count % 8),
+    ], dtype=np.float64)
+
+
+# =====================================================================
+# §5  Data generation
+# =====================================================================
+
+NUM_CLASSES = ("converged", "unconverged_short", "unconverged_long")
+TOK_CLASSES = ("bpe", "bpe+chatml", "bpe+llama2")
+
+# Bounds of the training distribution for n.  Any input outside this
+# range is flagged OOD by reconstruct_numeric.
+N_TRAIN_RANGE = (10, 40)
+
+
+def gen_numeric_examples(n_examples: int, seed: int = 0):
+    """Generate (features, label) pairs for the numeric task.
+
+    The label generator *knows* max_iter and tol.  The feature
+    extractor does not.  This is the honest setting: the config that
+    produced the value is not recoverable from the value.
+    """
+    rng = random.Random(seed)
+    X, y = [], []
+    while len(X) < n_examples:
+        n = rng.choice([10, 20, 30, 40])
+        max_iter = rng.choice([5, 50, 200, 2000, 20000])
+        tol = rng.choice([1e-3, 1e-6, 1e-10])
+        A = build_symmetric(n, rng)
+        eigs, info = jacobi_eigen(A, max_iter, tol)
+        value = participation_ratio(eigs)
+
+        if info["converged"]:
+            label = 0
+        elif max_iter <= 200:
+            label = 1
+        else:
+            label = 2
+
+        X.append(numeric_features(value, n))
+        y.append(label)
+    return np.stack(X), np.array(y, dtype=np.int64)
+
+
+_TOK_SAMPLES = [
+    "What is 2+2?",
+    "Translate 'hello' to Spanish.",
+    "Explain gradient descent.",
+    "How do I make sourdough bread?",
+    "Please prove step by step that the halting problem is "
+    "undecidable and analyze the consequences for formal verification.",
+    "Derive the Euler-Lagrange equation from the principle of "
+    "stationary action.",
+    "Summarize the Fano plane in one sentence.",
+    "What is the participation ratio?",
+    "Give me a synonym for 'happy'.",
+    "Write a haiku about autumn.",
+]
+
+
+def gen_token_examples(n_examples: int, seed: int = 0):
+    rng = random.Random(seed)
+    X, y = [], []
+    while len(X) < n_examples:
+        text = rng.choice(_TOK_SAMPLES)
+        template = rng.choice(["bpe", "bpe+chatml", "bpe+llama2"])
+        offset = {"bpe": 0, "bpe+chatml": 11, "bpe+llama2": 7}[template]
+        count = bpe_like(text) + offset
+        label = {"bpe": 0, "bpe+chatml": 1, "bpe+llama2": 2}[template]
+        X.append(token_features(count, text))
+        y.append(label)
+    return np.stack(X), np.array(y, dtype=np.int64)
+
+
+# =====================================================================
+# §6  Normalization
+# =====================================================================
+
+@dataclass
+class Normalizer:
+    mu: np.ndarray
+    sigma: np.ndarray
+
+    @classmethod
+    def fit(cls, X: np.ndarray) -> "Normalizer":
+        mu = X.mean(axis=0)
+        sigma = X.std(axis=0)
+        sigma[sigma < 1e-9] = 1.0
+        return cls(mu, sigma)
+
+    def __call__(self, X: np.ndarray) -> np.ndarray:
+        return (X - self.mu) / self.sigma
+
+
+# =====================================================================
+# §7  MLP (numpy only)
+# =====================================================================
+
+class MLP:
+    """Two-layer MLP with tanh hidden activation.
+
+    Sizes are deliberately tiny: 8 -> 32 -> 3.  With 8 inputs that is
+    8*32 + 32 + 32*3 + 3 = 387 parameters.  Trainable in numpy in
+    under a second on CPU.
+    """
+
+    def __init__(self, feat_dim: int, hidden: int, n_classes: int,
+                 seed: int = 0):
+        rng = np.random.default_rng(seed)
+        self.W1 = rng.standard_normal((feat_dim, hidden)) * 0.1
+        self.b1 = np.zeros(hidden)
+        self.W2 = rng.standard_normal((hidden, n_classes)) * 0.1
+        self.b2 = np.zeros(n_classes)
+        self.n_classes = n_classes
+        self.temperature = 1.0
+
+    def params(self):
+        return [self.W1, self.b1, self.W2, self.b2]
+
+    def forward(self, X: np.ndarray):
+        h = np.tanh(X @ self.W1 + self.b1)
+        logits = h @ self.W2 + self.b2
+        return h, logits
+
+    def predict_proba(self, X: np.ndarray) -> np.ndarray:
+        _, logits = self.forward(X)
+        z = logits / max(1e-6, self.temperature)
+        z = z - z.max(axis=1, keepdims=True)
+        e = np.exp(z)
+        return e / e.sum(axis=1, keepdims=True)
+
+    def predict(self, X: np.ndarray) -> np.ndarray:
+        return self.predict_proba(X).argmax(axis=1)
+
+    def loss_and_grad(self, X: np.ndarray, y: np.ndarray):
+        h, logits = self.forward(X)
+        n = len(y)
+        logits = logits - logits.max(axis=1, keepdims=True)
+        e = np.exp(logits)
+        p = e / e.sum(axis=1, keepdims=True)
+
+        loss = -np.log(p[np.arange(n), y] + 1e-12).mean()
+
+        dz = p.copy()
+        dz[np.arange(n), y] -= 1.0
+        dz /= n
+
+        dW2 = h.T @ dz
+        db2 = dz.sum(axis=0)
+        dh = dz @ self.W2.T
+        dpre = dh * (1.0 - h * h)
+        dW1 = X.T @ dpre
+        db1 = dpre.sum(axis=0)
+        return loss, [dW1, db1, dW2, db2]
+
+
+def train_mlp(model: MLP, X: np.ndarray, y: np.ndarray,
+              epochs: int = 200, batch: int = 64,
+              lr: float = 3e-3, seed: int = 0) -> List[float]:
+    rng = np.random.default_rng(seed)
+    m = [np.zeros_like(p) for p in model.params()]
+    v = [np.zeros_like(p) for p in model.params()]
+    t = 0
+    b1, b2, eps = 0.9, 0.999, 1e-8
+    losses = []
+
+    n = len(y)
+    for _ in range(epochs):
+        idx = rng.permutation(n)
+        for start in range(0, n, batch):
+            sel = idx[start:start + batch]
+            loss, grads = model.loss_and_grad(X[sel], y[sel])
+            t += 1
+            for i, (p, g) in enumerate(zip(model.params(), grads)):
+                m[i] = b1 * m[i] + (1 - b1) * g
+                v[i] = b2 * v[i] + (1 - b2) * g * g
+                mhat = m[i] / (1 - b1 ** t)
+                vhat = v[i] / (1 - b2 ** t)
+                p -= lr * mhat / (np.sqrt(vhat) + eps)
+        losses.append(loss)
+    return losses
+
+
+# =====================================================================
+# §8  Calibration
+# =====================================================================
+
+def expected_calibration_error(probs: np.ndarray, labels: np.ndarray,
+                               n_bins: int = 10) -> float:
+    """ECE: weighted |confidence - accuracy| across probability bins."""
+    conf = probs.max(axis=1)
+    pred = probs.argmax(axis=1)
+    correct = (pred == labels).astype(np.float64)
+    bins = np.linspace(0.0, 1.0, n_bins + 1)
+    total = 0.0
+    for lo, hi in zip(bins[:-1], bins[1:]):
+        mask = (conf >= lo) & (conf < hi)
+        if mask.sum() == 0:
+            continue
+        avg_conf = conf[mask].mean()
+        avg_acc = correct[mask].mean()
+        total += mask.sum() * abs(avg_conf - avg_acc)
+    return total / len(labels)
+
+
+def fit_temperature(model: MLP, X_val: np.ndarray,
+                    y_val: np.ndarray) -> float:
+    """Find T > 0 minimizing NLL of softmax(logits / T).
+
+    A saturated softmax returns p=1.00 on everything, which is a lie.
+    Temperature scaling turns that into a calibrated number.
+    """
+    _, logits = model.forward(X_val)
+    best_T, best_nll = 1.0, float("inf")
+    for T in np.linspace(0.5, 10.0, 40):
+        z = logits / T
+        z = z - z.max(axis=1, keepdims=True)
+        e = np.exp(z)
+        p = e / e.sum(axis=1, keepdims=True)
+        nll = -np.log(p[np.arange(len(y_val)), y_val] + 1e-12).mean()
+        if nll < best_nll:
+            best_nll, best_T = nll, float(T)
+    return best_T
+
+
+# =====================================================================
+# §9  Serializer
+# =====================================================================
+
+@dataclass
+class Reconstruction:
+    label: str
+    confidence: float
+    marker: str
+    ood: bool = False
+
+    def tag(self) -> str:
+        if self.ood:
+            return f"[{self.label}?? OOD p={self.confidence:.2f}]"
+        return f"[{self.label}{self.marker} p={self.confidence:.2f}]"
+
+    def __str__(self) -> str:
+        return self.tag()
+
+
+def serialize(label: str, confidence: float,
+              ood: bool = False,
+              high_conf: float = 0.90) -> Reconstruction:
+    marker = " " if confidence >= high_conf else "?"
+    return Reconstruction(label=label, confidence=confidence,
+                          marker=marker, ood=ood)
+
+
+def reconstruct_numeric(model: MLP, norm: Normalizer,
+                        value: float, n: int) -> Reconstruction:
+    ood = not (N_TRAIN_RANGE[0] <= n <= N_TRAIN_RANGE[1])
+    feats = numeric_features(value, n)
+    probs = model.predict_proba(norm(feats[None, :]))[0]
+    k = int(probs.argmax())
+    return serialize(NUM_CLASSES[k], float(probs[k]), ood=ood)
+
+
+def reconstruct_token(model: MLP, norm: Normalizer,
+                      count: int, text: str) -> Reconstruction:
+    feats = token_features(count, text)
+    probs = model.predict_proba(norm(feats[None, :]))[0]
+    k = int(probs.argmax())
+    return serialize(TOK_CLASSES[k], float(probs[k]))
+
+
+# =====================================================================
+# §10  Self-test
+# =====================================================================
+
+def self_test(verbose: bool = True) -> Tuple[int, int]:
+    checks = []
+
+    # Feature shapes are stable.
+    f = numeric_features(3.14, 20)
+    checks.append(("numeric features: 8 dims", f.shape == (8,)))
+    checks.append(("numeric features: finite", bool(np.all(np.isfinite(f)))))
+
+    t = token_features(54, "hello world")
+    checks.append(("token features: 12 dims", t.shape == (12,)))
+    checks.append(("token features: finite", bool(np.all(np.isfinite(t)))))
+
+    # The residual feature is the whole trick.
+    text = "The quick brown fox jumps over the lazy dog."
+    local = bpe_like(text)
+    f_none = token_features(local, text)
+    f_chat = token_features(local + 11, text)
+    f_llama = token_features(local + 7, text)
+    checks.append(("residual: raw bpe is ~0",
+                   abs(f_none[3]) < 1e-9))
+    checks.append(("residual: chatml is ~11",
+                   abs(f_chat[3] - 11) < 1e-9))
+    checks.append(("residual: llama2 is ~7",
+                   abs(f_llama[3] - 7) < 1e-9))
+
+    # Jacobi sanity: short budget does not converge, long does.
+    rng = random.Random(0)
+    A = build_symmetric(20, rng)
+    _, short = jacobi_eigen(A, max_iter=5, tol=1e-10)
+    _, long_ = jacobi_eigen(A, max_iter=20000, tol=1e-10)
+    checks.append(("jacobi: short does not converge",
+                   short["converged"] is False))
+    checks.append(("jacobi: long converges",
+                   long_["converged"] is True))
+
+    # Data generation.
+    Xn, yn = gen_numeric_examples(300, seed=1)
+    Xt, yt = gen_token_examples(300, seed=1)
+    checks.append(("gen_numeric: >=2 classes present",
+                   len(set(yn.tolist())) >= 2))
+    checks.append(("gen_token: 3 classes present",
+                   set(yt.tolist()) == {0, 1, 2}))
+
+    # Leakage guard: numeric features must NOT contain max_iter or tol.
+    # If accuracy on the honest task exceeds 0.85, a leaking feature
+    # has crept back in.
+    norm_n = Normalizer.fit(Xn)
+    m_n = MLP(Xn.shape[1], 16, 3, seed=2)
+    train_mlp(m_n, norm_n(Xn), yn, epochs=80, batch=32, seed=3)
+    acc_honest = (m_n.predict(norm_n(Xn)) == yn).mean()
+    checks.append(("numeric task: no leaking features (acc < 0.85)",
+                   acc_honest < 0.85))
+    checks.append(("numeric task: better than chance (acc > 0.40)",
+                   acc_honest > 0.40))
+
+    # Token task: must be near-perfect, because it is arithmetic.
+    norm_t = Normalizer.fit(Xt)
+    m_t = MLP(Xt.shape[1], 16, 3, seed=2)
+    train_mlp(m_t, norm_t(Xt), yt, epochs=80, batch=32, seed=3)
+    acc_t = (m_t.predict(norm_t(Xt)) == yt).mean()
+    checks.append(("token task: near-perfect (acc > 0.95)", acc_t > 0.95))
+
+    # ECE within range.
+    probs = m_t.predict_proba(norm_t(Xt))
+    e = expected_calibration_error(probs, yt)
+    checks.append(("ECE is a probability", 0.0 <= e <= 1.0))
+
+    # Serializer marks high confidence without '?'.
+    r_hi = serialize("bpe+chatml", 0.95)
+    r_lo = serialize("bpe+chatml", 0.60)
+    r_ood = serialize("converged", 0.99, ood=True)
+    checks.append(("serialize: high conf no '?'", r_hi.marker == " "))
+    checks.append(("serialize: low conf has '?'", r_lo.marker == "?"))
+    checks.append(("serialize: OOD emits '??'",
+                   "??" in r_ood.tag() and "OOD" in r_ood.tag()))
+
+    # Temperature scaling is a positive finite number.
+    T = fit_temperature(m_t, norm_t(Xt), yt)
+    checks.append(("temperature fit: positive finite", 0.5 <= T <= 10.0))
+
+    passed = sum(1 for _, ok in checks if ok)
+    if verbose:
+        print()
+        print("=" * 74)
+        print("SELF-TEST")
+        print("=" * 74)
+        for name, ok in checks:
+            mark = "PASS" if ok else "FAIL"
+            print(f"  [{mark}] {name}")
+        print()
+        print(f"  {passed}/{len(checks)} correct")
+    return passed, len(checks)
+
+
+# =====================================================================
+# §11  Demo
+# =====================================================================
+
+def _section(title: str, width: int = 74) -> None:
+    print()
+    print("=" * width)
+    print(title)
+    print("=" * width)
+
+
+DEMO_NUMERIC = [
+    # (value, n, human note) — note reflects the class boundaries
+    # actually used in gen_numeric_examples: converged / short / long
+    (5.9971, 20, "long budget, tight tol"),
+    (3.1400, 20, "5 iterations"),
+    (8.9539, 40, "budget 500 for n=40"),
+    (1.7109, 20, "loose tol, long budget"),
+    (2.0000, 10, "5 iterations"),
+    (4.5000, 30, "tight tol, budget 200"),
+    (97.705, 448, "OUT OF DISTRIBUTION (n=448)"),
+]
+
+DEMO_TOKENS = [
+    ("What is 2+2?", 5, "raw bpe count for a short prompt"),
+    ("What is 2+2?", 16, "same text, +11 -> chatml"),
+    ("What is 2+2?", 12, "same text, +7 -> llama2"),
+]
+
+
+def demo() -> None:
+    print()
+    print("=" * 74)
+    print("PROVENANCE RECONSTRUCT — infer the encoding a bare value lost")
+    print("=" * 74)
+
+    passed, total = self_test(verbose=True)
+    if passed < total:
+        print("\n  self-test failed; continuing anyway for inspection.")
+
+    _section("PART 0  Train the numeric reconstructor")
+    print("  Features deliberately EXCLUDE max_iter and tol.")
+    print("  Those are solver config, not observable at reconstruction time.")
+    print("  If they were included, accuracy would be ~99% by leakage.")
+    print()
+    Xn_tr, yn_tr = gen_numeric_examples(2000, seed=10)
+    Xn_te, yn_te = gen_numeric_examples(400, seed=11)
+    norm_n = Normalizer.fit(Xn_tr)
+    model_n = MLP(Xn_tr.shape[1], 32, 3, seed=20)
+    train_mlp(model_n, norm_n(Xn_tr), yn_tr, epochs=250, batch=64, seed=21)
+    model_n.temperature = fit_temperature(model_n, norm_n(Xn_te), yn_te)
+
+    acc_n = (model_n.predict(norm_n(Xn_te)) == yn_te).mean()
+    probs_n = model_n.predict_proba(norm_n(Xn_te))
+    ece_n = expected_calibration_error(probs_n, yn_te)
+    n_params = sum(p.size for p in model_n.params())
+
+    # Majority-class baseline.
+    counts = np.bincount(yn_tr, minlength=3)
+    baseline = counts.max() / counts.sum()
+
+    print(f"  training examples   : {len(yn_tr)}")
+    print(f"  test examples       : {len(yn_te)}")
+    print(f"  parameters          : {n_params}")
+    print(f"  majority baseline   : {baseline * 100:.1f}%")
+    print(f"  test accuracy       : {acc_n * 100:.1f}%")
+    print(f"  temperature         : {model_n.temperature:.2f}")
+    print(f"  expected calib err  : {ece_n:.3f}")
+    print(f"  classes             : {', '.join(NUM_CLASSES)}")
+    print()
+    if acc_n - baseline < 0.15:
+        print("  FINDING: this is a negative result.  A bare value + matrix")
+        print("  size carries little recoverable convergence information.")
+        print("  The residual and iteration count are the only reliable")
+        print("  signals, and they were lost with the provenance.")
+    else:
+        print("  Note: accuracy is meaningfully above baseline.  Check that")
+        print("  no leaking feature was reintroduced.")
+
+    _section("PART 1  Train the token reconstructor")
+    Xt_tr, yt_tr = gen_token_examples(2000, seed=10)
+    Xt_te, yt_te = gen_token_examples(400, seed=11)
+    norm_t = Normalizer.fit(Xt_tr)
+    model_t = MLP(Xt_tr.shape[1], 32, 3, seed=20)
+    train_mlp(model_t, norm_t(Xt_tr), yt_tr, epochs=250, batch=64, seed=21)
+    model_t.temperature = fit_temperature(model_t, norm_t(Xt_te), yt_te)
+
+    acc_t = (model_t.predict(norm_t(Xt_te)) == yt_te).mean()
+    probs_t = model_t.predict_proba(norm_t(Xt_te))
+    ece_t = expected_calibration_error(probs_t, yt_te)
+    n_params_t = sum(p.size for p in model_t.params())
+
+    print(f"  training examples   : {len(yt_tr)}")
+    print(f"  test examples       : {len(yt_te)}")
+    print(f"  parameters          : {n_params_t}")
+    print(f"  test accuracy       : {acc_t * 100:.1f}%")
+    print(f"  temperature         : {model_t.temperature:.2f}")
+    print(f"  expected calib err  : {ece_t:.3f}")
+    print(f"  classes             : {', '.join(TOK_CLASSES)}")
+    print()
+    print("  FINDING: this is a positive result, but a trivial one.  The")
+    print("  residual (count - bpe_like) equals the template offset exactly.")
+    print("  The model learns a threshold lookup.  It is arithmetic dressed")
+    print("  as learning; the model is packaging, not magic.")
+
+    _section("PART 2  Reconstruct provenance for bare numerical outputs")
+    print(f"  {'value':>10}  {'n':>4}  {'note':<36}  reconstruction")
+    print("  " + "-" * 80)
+    for value, n, note in DEMO_NUMERIC:
+        r = reconstruct_numeric(model_n, norm_n, value, n)
+        print(f"  {value:>10.4f}  {n:>4}  {note:<36}  {r.tag()}")
+
+    _section("PART 3  Reconstruct provenance for bare token counts")
+    print(f"  {'text':<24}  {'count':>5}  {'note':<34}  reconstruction")
+    print("  " + "-" * 84)
+    for text, count, note in DEMO_TOKENS:
+        r = reconstruct_token(model_t, norm_t, count, text)
+        display = text[:21] + "..." if len(text) > 21 else text
+        print(f"  {display:<24}  {count:>5}  {note:<34}  {r.tag()}")
+
+    _section("PART 4  The residual trick, isolated")
+    print("  Subtract a free local estimate from the observed count.")
+    print("  The residual lands on 0, +7, or +11 -- the template offset.")
+    print("  This is why the token task is trivially learnable.")
+    print()
+    text = "What is 2+2?"
+    local = bpe_like(text)
+    print(f"  text                  : {text!r}")
+    print(f"  local bpe_like(text)  : {local}")
+    print()
+    print(f"  {'observed count':>14}  {'residual':>9}  inferred template")
+    print("  " + "-" * 50)
+    for observed, template in [(local, "none"),
+                               (local + 7, "llama2"),
+                               (local + 11, "chatml")]:
+        r = reconstruct_token(model_t, norm_t, observed, text)
+        print(f"  {observed:>14}  {observed - local:>+9}  "
+              f"{r.tag()}   (truth: {template})")
+
+    _section("PART 5  Honesty: observed, inferred, OOD")
+    print("  numerical_provenance.py emits OBSERVED provenance:")
+    print(f"    PValue(97.705, UNCONVERGED, via='jacobi_short')")
+    print(f"    report -> <refused: 'jacobi_short' did not converge; "
+          f"residual=8.47, tol=1e-10>")
+    print()
+    print("  token_provenance.py emits OBSERVED provenance:")
+    print(f"    TCount(54, ok, via='toy-bpe')  encoding=toy-bpe+chatml(+11)")
+    print(f"    report -> 54 [toy-bpe+chatml(+11)]")
+    print()
+    print("  provenance_reconstruct.py emits INFERRED provenance:")
+    r_in = reconstruct_numeric(model_n, norm_n, 5.9971, 20)
+    print(f"    5.9971, n=20      -> {r_in.tag()}")
+    print()
+    print("  ...and OOD provenance for inputs outside the training range:")
+    r_ood = reconstruct_numeric(model_n, norm_n, 97.705, 448)
+    print(f"    97.705, n=448     -> {r_ood.tag()}")
+    print()
+    print("  A downstream consumer can require OBSERVED provenance and")
+    print("  refuse INFERRED provenance.  The markers:")
+    print("    [LABEL p=0.95]          high-confidence inferred")
+    print("    [LABEL? p=0.61]         low-confidence inferred")
+    print("    [LABEL?? OOD p=1.00]    out of distribution; p is meaningless")
+
+    _section("PART 6  What this file demonstrates")
+    print("""
+  Positive result (token side):
+    The residual against a free local estimator equals the template
+    offset exactly.  A tiny model learns the classification at ~100%
+    because the task is arithmetic, not inference.
+
+  Negative result (numeric side):
+    A bare value plus its matrix size carries little recoverable
+    convergence information.  Accuracy lands near the majority-class
+    baseline.  The residual and iteration count are the only reliable
+    signals, and they were lost with the provenance.
+
+  The value of the tool is therefore not its accuracy on the numeric
+  side.  It is its CALIBRATION (temperature-scaled confidence) and its
+  OOD MARKER (the double '??' plus 'OOD' when the input is outside the
+  training range).  A reconstructor that says "[? p=0.61]" is more
+  useful than one that lies at p=1.00.
+
+  The family:
+    numerical_provenance.py     -> PValue   (observed)
+    token_provenance.py         -> TCount   (observed)
+    provenance_reconstruct.py   -> tag      (inferred, calibrated, OOD-aware)
+""")
+
+
+# =====================================================================
+# §12  Entry point
+# =====================================================================
+
+if __name__ == "__main__":
+    demo()