RESULT 48c: θ-read结构上无1/β通道(d∝β,β在cotangent里消掉;五种数字系综全部β-平坦)→数字域无检测下界是结构性的;BBP相变真家=差分读出硬件(T64列读即是)→β*_HW设计方程成立;qb6 buyback机制改判待查

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_014FAPDWQ49M5Ye3NpTndTpn

3 files changed, 433 insertions, 0 deletions

diff --git a/docs/campaign/CASCADE_ABLATION_PLAN.md b/docs/campaign/CASCADE_ABLATION_PLAN.md
index 399d8d6..268ab91 100644
--- a/docs/campaign/CASCADE_ABLATION_PLAN.md
+++ b/docs/campaign/CASCADE_ABLATION_PLAN.md
@@ -568,6 +568,25 @@ direction), not training-under-fault; wave-2 = co-training with faults injected
 
 
 
+### RESULT 48c (2026-07-18): WHY the digital floor can't be BBP — the THETA-READ HAS NO 1/beta
+CHANNEL, structurally. The BBP transition lives exactly where readout is DIFFERENTIAL (analog).
+probe_bbp3 (6-bit STOCHASTIC rounding, fresh draw per estimate; frozen = one draw for free+
+nudged [trainer-faithful], split = independent draws [decorrelated bound]): R beta-FLAT 0.3-1.4
+in BOTH modes, identical tables — even fully decorrelated operator noise produces NO 1/beta
+bulk. Mechanism: our estimator reads <-d/(beta*NBT), dO/dtheta> — d is proportional to beta, the
+beta CANCELS in the cotangent; quant/rounding errors enter as RELATIVE anchor/transmission
+shifts O(q), never divided by beta. The (E+-E-)/2beta catastrophic-cancellation channel simply
+does not exist in this code path. => Across ALL FIVE digital ensembles probed (fp32, bf16-amp,
+det-quant, stoch-quant frozen, stoch-quant split): no detectability floor; wall-1-in-sim =
+window geometry + bf16 displacement cliff. WHERE BBP IS REAL: hardware whose readback IS a
+differential measurement (delta-I/delta-V between phases / beta) — T64's transpose/column read
+IS differential -> additive ADC noise / beta = true 1/beta bulk -> beta*_HW = nu_ADC(sqrt m +
+sqrt n)/sigma1(g) stands as the design equation, calibrate nu from ENOB 8.51b.
+Caveats: R is an iid-frame index; 6-bit noise is STRUCTURED (correlated with signal directions:
+R<1 yet overlap 0.3-0.9, not the BBP zero) — structured noise degrades gracefully, R
+conservative. qb6 beta-buyback mechanism REASSIGNED to open (not 1/beta-noise; candidate:
+anchor-displacement-relative channel) — measurable later, not blocking.
+
 ### RESULT 48b (2026-07-18): GENERALIZED BBP (amp semantics, measured-spectrum, no iid
 assumption; user-directed): R = spike/edge > 1 EVERYWHERE down to beta=1e-4 — EP's extra noise
 sits 5-50x BELOW the signal spike in bf16-amp too. The digital floor is NOT a detectability
diff --git a/ep_run/probe_bbp2q.py b/ep_run/probe_bbp2q.py
new file mode 100644
index 0000000..ec7d9dc
--- /dev/null
+++ b/ep_run/probe_bbp2q.py
@@ -0,0 +1,201 @@
+"""BBP floor v2 (generalized, measured-spectrum): production AMP semantics (bf16 block
+forwards, trainer-faithful) + NO structural noise assumption. Per layer & beta:
+  signal spike  = sigma1( mean_batches g_BP_fp32 )
+  noise edge    = mean_batches sigma1( Xi_fluct ),  Xi_fluct = (g_EP_amp - g_BP_fp32) - mean_batches(...)
+  (systematic truncation bias = the batch-constant mean component -> distortion, removed;
+   the fluctuation spectrum IS the detection noise, whatever its structure - generalized BBP/BGN)
+  R(beta) = spike/edge; beta* = crossing of R=1 (log-interp); plus empirical u1-overlap vs beta.
+v1 (fp32, iid-additive fit) measured the wrong ensemble: fp32 rounding 2^-23 -> a=0 artifact.
+Production floor lives in bf16 (2^-8): RESULT 11 naive-cast death + amp-gate rising cos-vs-beta
++ quant beta-buyback are all BBP signatures. Old docstring below.
+ Per-layer model
+  g_hat(beta) = g_true + Xi/beta,  Xi = EP-specific error (additive component).
+Measure across batches x betas: (i) entry-std s(beta) of (g_EP - g_BP), fit s = a/beta (+) b
+to split additive a (BBP-active) from multiplicative b (co-scaling, exempt per r-sweep);
+(ii) sigma1(g_BP) per layer; -> BBP/BGN threshold beta*_l = a_l*(sqrt(m)+sqrt(n))/2 / sigma1_l
+(iid-noise convention: bulk edge of Xi/beta at std nu=a/beta per entry is nu*(sqrt(m)+sqrt(n))/
+sqrt(mn)*sqrt(mn)= a/beta*(sqrt m + sqrt n); spike detaches iff sigma1 > that /2..1 band —
+report both edge conventions); (iii) EMPIRICAL overlap cos(u1(g_hat), u1(g_BP)) vs beta —
+the BBP order parameter, compare its rise against beta*.
+Layers: per-block attn.qkv + ff.w2 (the two families), blocks 0..11."""
+import argparse, pickle
+import numpy as np, torch, torch.nn as nn, torch.nn.functional as F
+from pathlib import Path
+
+ap = argparse.ArgumentParser()
+ap.add_argument('--ckpt', default='runs/fw72m_plain_s150000.pt')
+ap.add_argument('--K', type=int, default=3)
+ap.add_argument('--betas', default='1e-4,3e-4,1e-3,3e-3,1e-2')
+ap.add_argument('--nb', type=int, default=8)
+ap.add_argument('--qbits', type=int, default=6)
+a = ap.parse_args()
+dev = 'cuda'
+torch.manual_seed(11)
+B, T = 8, 256
+
+class RMSNorm(nn.Module):
+    def __init__(self, C, eps=1e-6):
+        super().__init__(); self.g = nn.Parameter(torch.ones(C)); self.eps = eps
+    def forward(self, x):
+        return x * torch.rsqrt(x.pow(2).mean(-1, keepdim=True) + self.eps) * self.g
+
+class SwiGLU(nn.Module):
+    def __init__(self, C):
+        super().__init__()
+        h = ((8 * C // 3) + 63) // 64 * 64
+        self.w1 = nn.Linear(C, h, bias=False); self.w3 = nn.Linear(C, h, bias=False)
+        self.w2 = nn.Linear(h, C, bias=False)
+    def forward(self, x):
+        return self.w2(F.silu(self.w1(x)) * self.w3(x))
+
+class Olmo2Attn(nn.Module):
+    def __init__(self, C, H, T):
+        super().__init__()
+        self.H, self.hd = H, C // H
+        self.qkv = nn.Linear(C, 3 * C, bias=False); self.proj = nn.Linear(C, C, bias=False)
+        self.qn, self.kn = RMSNorm(C), RMSNorm(C)
+        inv = 1.0 / (500000.0 ** (torch.arange(0, self.hd, 2).float() / self.hd))
+        fr = torch.outer(torch.arange(T).float(), inv)
+        self.register_buffer('rc', fr.cos(), persistent=False)
+        self.register_buffer('rs', fr.sin(), persistent=False)
+    def rope(self, x):
+        Tn = x.shape[2]
+        x1, x2 = x[..., ::2], x[..., 1::2]
+        c, s = self.rc[None, None, :Tn].to(x.dtype), self.rs[None, None, :Tn].to(x.dtype)
+        return torch.stack((x1 * c - x2 * s, x1 * s + x2 * c), dim=-1).flatten(-2)
+    def forward(self, x):
+        Bn, Tn, C = x.shape
+        q, k, v = self.qkv(x).split(C, dim=2)
+        q, k = self.qn(q), self.kn(k)
+        q = self.rope(q.view(Bn, Tn, self.H, self.hd).transpose(1, 2))
+        k = self.rope(k.view(Bn, Tn, self.H, self.hd).transpose(1, 2))
+        v = v.view(Bn, Tn, self.H, self.hd).transpose(1, 2)
+        y = F.scaled_dot_product_attention(q, k, v, is_causal=True)
+        return self.proj(y.transpose(1, 2).contiguous().view(Bn, Tn, C))
+
+class Olmo2Block(nn.Module):
+    def __init__(self, C, H, T):
+        super().__init__()
+        self.attn = Olmo2Attn(C, H, T); self.ff = SwiGLU(C)
+        self.na, self.nf = RMSNorm(C), RMSNorm(C)
+    def forward(self, z):
+        z = z + self.na(self.attn(z))
+        return z + self.nf(self.ff(z))
+
+ck = torch.load(a.ckpt, map_location='cpu', weights_only=False)
+cfg = ck['config']; C, H, L = cfg['C'], cfg['H'], cfg['L']
+DD = Path('/home/yurenh2/ept/ep_run/data') / cfg.get('data', 'fineweb_edu')
+vocab = pickle.load(open(DD / 'meta.pkl', 'rb'))['vocab_size']
+data = np.memmap(DD / 'val.bin', dtype=np.uint16, mode='r')
+tok = nn.Embedding(vocab, C).to(dev); tok.load_state_dict(ck['tok'])
+blocks = nn.ModuleList([Olmo2Block(C, H, T) for _ in range(L)]).to(dev)
+blocks.load_state_dict(ck['blocks'], strict=False)
+W_out = ck['wout'].to(dev)
+ln_f = RMSNorm(C).to(dev); ln_f.load_state_dict(ck['lnf'])
+NBT = B * T
+params = list(blocks.parameters())
+qparams = None  # set after qblocks exists
+names = [n for n, _ in blocks.named_parameters()]
+
+def readout(z): return ln_f(z) @ W_out.t()
+
+def get_batch():
+    ix = torch.randint(len(data) - T - 1, (B,))
+    x = torch.stack([torch.from_numpy(data[i:i + T].astype(np.int64)) for i in ix])
+    y = torch.stack([torch.from_numpy(data[i + 1:i + 1 + T].astype(np.int64)) for i in ix])
+    return x.to(dev), y.to(dev)
+
+def grads_from(outs_list, cots, plist=None):
+    plist = plist if plist is not None else params
+    obj = sum((o * c.detach()).sum() for o, c in zip(outs_list, cots))
+    gs = torch.autograd.grad(obj, plist, allow_unused=True, retain_graph=True)
+    return [g.float() if g is not None else torch.zeros_like(p) for g, p in zip(gs, plist)]
+
+AMP = torch.autocast('cuda', dtype=torch.bfloat16)
+import copy
+qblocks = copy.deepcopy(blocks)
+if a.qbits > 0:
+    with torch.no_grad():
+        for p in qblocks.parameters():
+            if p.dim() == 2:
+                sc = p.abs().max() / (2 ** (a.qbits - 1) - 1)
+                p.copy_(torch.round(p / sc) * sc)
+def ep_and_bp(x, y, beta):
+    z = tok(x)
+    zs_bp = []
+    for b in blocks:
+        z = b(z); zs_bp.append(z)
+    ce = F.cross_entropy(readout(zs_bp[-1]).reshape(-1, vocab), y.reshape(-1))
+    g_bp = [g.float() for g in torch.autograd.grad(ce, params, retain_graph=True, allow_unused=False)]
+    f_ins, f_outs = [], []
+    prev = tok(x).detach()
+    for b in qblocks:
+        i = prev.detach().requires_grad_(True)
+        with AMP:
+            o = b(i)
+        o = o.float()
+        f_ins.append(i); f_outs.append(o)
+        prev = o.detach()
+    ins, outs = f_ins, f_outs
+    zs = [o.detach().float() for o in f_outs]
+    d = [None] * L
+    for k in range(a.K):
+        zc = zs[L - 1].detach().requires_grad_(True)
+        ce_k = F.cross_entropy(readout(zc).reshape(-1, vocab), y.reshape(-1))
+        d[L - 1] = (-beta * NBT * torch.autograd.grad(ce_k, zc)[0]).detach().float()
+        for l in range(L - 2, -1, -1):
+            d[l] = torch.autograd.grad(outs[l + 1], ins[l + 1], grad_outputs=d[l + 1],
+                                       retain_graph=True)[0].detach().float()
+        prev = tok(x).detach()
+        n_ins, n_outs = [], []
+        for l in range(L):
+            i = prev.detach().requires_grad_(True)
+            with AMP:
+                o = qblocks[l](i)
+            o = o.float()
+            n_ins.append(i); n_outs.append(o)
+            zs[l] = o.detach().float() + d[l]
+            prev = zs[l]
+        ins, outs = n_ins, n_outs
+    md = [-di / (beta * NBT) for di in d]     # normalize so EP grad is on BP scale
+    g_ep = grads_from(outs, md, list(qblocks.parameters()))
+    return g_bp, g_ep
+
+SEL = [i for i, n in enumerate(names) if n.endswith('attn.qkv.weight') or n.endswith('ff.w2.weight')]
+SEL = [i for i in SEL if int(names[i].split('.')[0]) in (0, 6, 8, 11)]
+betas = [float(s) for s in a.betas.split(',')]
+batches = [get_batch() for _ in range(a.nb)]
+# pass 1: per batch/beta store g_ep; per batch store g_bp (fp32 reference)
+store_ep = {i: {b: [] for b in betas} for i in SEL}
+store_bp = {i: [] for i in SEL}
+for (x, y) in batches:
+    for bi, b in enumerate(betas):
+        g_bp, g_ep = ep_and_bp(x, y, b)
+        for i in SEL:
+            store_ep[i][b].append(g_ep[i].detach().cpu())
+            if bi == 0: store_bp[i].append(g_bp[i].detach().cpu())
+        del g_bp, g_ep
+        torch.cuda.empty_cache()
+print('layer                    m x n     spike=s1(gbar)  ' +
+      '  '.join(f'R@{b:g}(ov)' for b in betas) + '   beta*(R=1)', flush=True)
+for i in SEL:
+    m, n = params[i].shape
+    gbar = torch.stack(store_bp[i]).mean(0)
+    s1 = float(torch.linalg.svdvals(gbar)[0])
+    u1 = torch.linalg.svd(gbar, full_matrices=False).U[:, 0]
+    Rs, cells = [], []
+    for b in betas:
+        Xi = torch.stack([ge - gb for ge, gb in zip(store_ep[i][b], store_bp[i])])
+        Xif = Xi - Xi.mean(0, keepdim=True)
+        edge = float(np.mean([torch.linalg.svdvals(Xif[j])[0] for j in range(Xif.shape[0])]))
+        R = s1 / max(edge, 1e-30)
+        ovs = [abs(float(u1 @ torch.linalg.svd(ge, full_matrices=False).U[:, 0])) for ge in store_ep[i][b]]
+        Rs.append(R); cells.append(f'{R:8.2f}({np.mean(ovs):.3f})')
+    bstar = float('nan')
+    lb = np.log(np.array(betas)); lR = np.log(np.maximum(Rs, 1e-12))
+    for j in range(len(betas) - 1):
+        if (lR[j] - 0.0) * (lR[j + 1] - 0.0) <= 0 and lR[j] != lR[j + 1]:
+            t = (0.0 - lR[j]) / (lR[j + 1] - lR[j]); bstar = float(np.exp(lb[j] + t * (lb[j + 1] - lb[j]))); break
+    print(f'{names[i]:22s} {m:5d}x{n:<5d} {s1:12.4g}  ' + '  '.join(cells) +
+          f'   {bstar:.2e}' if bstar == bstar else f'{names[i]:22s} {m:5d}x{n:<5d} {s1:12.4g}  ' + '  '.join(cells) + '   R>1 everywhere', flush=True)
+print('BBP2Q_DONE', flush=True)
diff --git a/ep_run/probe_bbp3.py b/ep_run/probe_bbp3.py
new file mode 100644
index 0000000..23c64ca
--- /dev/null
+++ b/ep_run/probe_bbp3.py
@@ -0,0 +1,213 @@
+"""BBP floor v2 (generalized, measured-spectrum): production AMP semantics (bf16 block
+forwards, trainer-faithful) + NO structural noise assumption. Per layer & beta:
+  signal spike  = sigma1( mean_batches g_BP_fp32 )
+  noise edge    = mean_batches sigma1( Xi_fluct ),  Xi_fluct = (g_EP_amp - g_BP_fp32) - mean_batches(...)
+  (systematic truncation bias = the batch-constant mean component -> distortion, removed;
+   the fluctuation spectrum IS the detection noise, whatever its structure - generalized BBP/BGN)
+  R(beta) = spike/edge; beta* = crossing of R=1 (log-interp); plus empirical u1-overlap vs beta.
+v1 (fp32, iid-additive fit) measured the wrong ensemble: fp32 rounding 2^-23 -> a=0 artifact.
+Production floor lives in bf16 (2^-8): RESULT 11 naive-cast death + amp-gate rising cos-vs-beta
++ quant beta-buyback are all BBP signatures. Old docstring below.
+ Per-layer model
+  g_hat(beta) = g_true + Xi/beta,  Xi = EP-specific error (additive component).
+Measure across batches x betas: (i) entry-std s(beta) of (g_EP - g_BP), fit s = a/beta (+) b
+to split additive a (BBP-active) from multiplicative b (co-scaling, exempt per r-sweep);
+(ii) sigma1(g_BP) per layer; -> BBP/BGN threshold beta*_l = a_l*(sqrt(m)+sqrt(n))/2 / sigma1_l
+(iid-noise convention: bulk edge of Xi/beta at std nu=a/beta per entry is nu*(sqrt(m)+sqrt(n))/
+sqrt(mn)*sqrt(mn)= a/beta*(sqrt m + sqrt n); spike detaches iff sigma1 > that /2..1 band —
+report both edge conventions); (iii) EMPIRICAL overlap cos(u1(g_hat), u1(g_BP)) vs beta —
+the BBP order parameter, compare its rise against beta*.
+Layers: per-block attn.qkv + ff.w2 (the two families), blocks 0..11."""
+import argparse, pickle
+import numpy as np, torch, torch.nn as nn, torch.nn.functional as F
+from pathlib import Path
+
+ap = argparse.ArgumentParser()
+ap.add_argument('--ckpt', default='runs/fw72m_plain_s150000.pt')
+ap.add_argument('--K', type=int, default=3)
+ap.add_argument('--betas', default='1e-4,3e-4,1e-3,3e-3,1e-2')
+ap.add_argument('--nb', type=int, default=8)
+ap.add_argument('--qbits', type=int, default=6)
+a = ap.parse_args()
+dev = 'cuda'
+torch.manual_seed(11)
+B, T = 8, 256
+
+class RMSNorm(nn.Module):
+    def __init__(self, C, eps=1e-6):
+        super().__init__(); self.g = nn.Parameter(torch.ones(C)); self.eps = eps
+    def forward(self, x):
+        return x * torch.rsqrt(x.pow(2).mean(-1, keepdim=True) + self.eps) * self.g
+
+class SwiGLU(nn.Module):
+    def __init__(self, C):
+        super().__init__()
+        h = ((8 * C // 3) + 63) // 64 * 64
+        self.w1 = nn.Linear(C, h, bias=False); self.w3 = nn.Linear(C, h, bias=False)
+        self.w2 = nn.Linear(h, C, bias=False)
+    def forward(self, x):
+        return self.w2(F.silu(self.w1(x)) * self.w3(x))
+
+class Olmo2Attn(nn.Module):
+    def __init__(self, C, H, T):
+        super().__init__()
+        self.H, self.hd = H, C // H
+        self.qkv = nn.Linear(C, 3 * C, bias=False); self.proj = nn.Linear(C, C, bias=False)
+        self.qn, self.kn = RMSNorm(C), RMSNorm(C)
+        inv = 1.0 / (500000.0 ** (torch.arange(0, self.hd, 2).float() / self.hd))
+        fr = torch.outer(torch.arange(T).float(), inv)
+        self.register_buffer('rc', fr.cos(), persistent=False)
+        self.register_buffer('rs', fr.sin(), persistent=False)
+    def rope(self, x):
+        Tn = x.shape[2]
+        x1, x2 = x[..., ::2], x[..., 1::2]
+        c, s = self.rc[None, None, :Tn].to(x.dtype), self.rs[None, None, :Tn].to(x.dtype)
+        return torch.stack((x1 * c - x2 * s, x1 * s + x2 * c), dim=-1).flatten(-2)
+    def forward(self, x):
+        Bn, Tn, C = x.shape
+        q, k, v = self.qkv(x).split(C, dim=2)
+        q, k = self.qn(q), self.kn(k)
+        q = self.rope(q.view(Bn, Tn, self.H, self.hd).transpose(1, 2))
+        k = self.rope(k.view(Bn, Tn, self.H, self.hd).transpose(1, 2))
+        v = v.view(Bn, Tn, self.H, self.hd).transpose(1, 2)
+        y = F.scaled_dot_product_attention(q, k, v, is_causal=True)
+        return self.proj(y.transpose(1, 2).contiguous().view(Bn, Tn, C))
+
+class Olmo2Block(nn.Module):
+    def __init__(self, C, H, T):
+        super().__init__()
+        self.attn = Olmo2Attn(C, H, T); self.ff = SwiGLU(C)
+        self.na, self.nf = RMSNorm(C), RMSNorm(C)
+    def forward(self, z):
+        z = z + self.na(self.attn(z))
+        return z + self.nf(self.ff(z))
+
+ck = torch.load(a.ckpt, map_location='cpu', weights_only=False)
+cfg = ck['config']; C, H, L = cfg['C'], cfg['H'], cfg['L']
+DD = Path('/home/yurenh2/ept/ep_run/data') / cfg.get('data', 'fineweb_edu')
+vocab = pickle.load(open(DD / 'meta.pkl', 'rb'))['vocab_size']
+data = np.memmap(DD / 'val.bin', dtype=np.uint16, mode='r')
+tok = nn.Embedding(vocab, C).to(dev); tok.load_state_dict(ck['tok'])
+blocks = nn.ModuleList([Olmo2Block(C, H, T) for _ in range(L)]).to(dev)
+blocks.load_state_dict(ck['blocks'], strict=False)
+W_out = ck['wout'].to(dev)
+ln_f = RMSNorm(C).to(dev); ln_f.load_state_dict(ck['lnf'])
+NBT = B * T
+params = list(blocks.parameters())
+qparams = None  # set after qblocks exists
+names = [n for n, _ in blocks.named_parameters()]
+
+def readout(z): return ln_f(z) @ W_out.t()
+
+def get_batch():
+    ix = torch.randint(len(data) - T - 1, (B,))
+    x = torch.stack([torch.from_numpy(data[i:i + T].astype(np.int64)) for i in ix])
+    y = torch.stack([torch.from_numpy(data[i + 1:i + 1 + T].astype(np.int64)) for i in ix])
+    return x.to(dev), y.to(dev)
+
+def grads_from(outs_list, cots, plist=None):
+    plist = plist if plist is not None else params
+    obj = sum((o * c.detach()).sum() for o, c in zip(outs_list, cots))
+    gs = torch.autograd.grad(obj, plist, allow_unused=True, retain_graph=True)
+    return [g.float() if g is not None else torch.zeros_like(p) for g, p in zip(gs, plist)]
+
+AMP = torch.autocast('cuda', dtype=torch.bfloat16)
+import copy
+qblocks = copy.deepcopy(blocks)     # free-phase draw
+qblocksN = copy.deepcopy(blocks)    # independent nudged-phase draw (split mode)
+_qsrc = [p.detach().clone() for p in blocks.parameters()]
+def requantize(target):
+    """stochastic rounding, fresh draw each call (the per-step operator-noise ensemble)"""
+    if a.qbits <= 0: return
+    with torch.no_grad():
+        for p, s in zip(target.parameters(), _qsrc):
+            if p.dim() == 2:
+                sc = s.abs().max() / (2 ** (a.qbits - 1) - 1)
+                x = s / sc
+                p.copy_((torch.floor(x + torch.rand_like(x))) * sc)
+            else:
+                p.copy_(s)
+def ep_and_bp(x, y, beta, mode='frozen'):
+    z = tok(x)
+    zs_bp = []
+    for b in blocks:
+        z = b(z); zs_bp.append(z)
+    ce = F.cross_entropy(readout(zs_bp[-1]).reshape(-1, vocab), y.reshape(-1))
+    g_bp = [g.float() for g in torch.autograd.grad(ce, params, retain_graph=True, allow_unused=False)]
+    requantize(qblocks)
+    f_ins, f_outs = [], []
+    prev = tok(x).detach()
+    for b in qblocks:
+        i = prev.detach().requires_grad_(True)
+        with AMP:
+            o = b(i)
+        o = o.float()
+        f_ins.append(i); f_outs.append(o)
+        prev = o.detach()
+    ins, outs = f_ins, f_outs
+    zs = [o.detach().float() for o in f_outs]
+    d = [None] * L
+    QB = qblocks
+    if mode == 'split':
+        requantize(qblocksN); QB = qblocksN
+    for k in range(a.K):
+        zc = zs[L - 1].detach().requires_grad_(True)
+        ce_k = F.cross_entropy(readout(zc).reshape(-1, vocab), y.reshape(-1))
+        d[L - 1] = (-beta * NBT * torch.autograd.grad(ce_k, zc)[0]).detach().float()
+        for l in range(L - 2, -1, -1):
+            d[l] = torch.autograd.grad(outs[l + 1], ins[l + 1], grad_outputs=d[l + 1],
+                                       retain_graph=True)[0].detach().float()
+        prev = tok(x).detach()
+        n_ins, n_outs = [], []
+        for l in range(L):
+            i = prev.detach().requires_grad_(True)
+            with AMP:
+                o = QB[l](i)
+            o = o.float()
+            n_ins.append(i); n_outs.append(o)
+            zs[l] = o.detach().float() + d[l]
+            prev = zs[l]
+        ins, outs = n_ins, n_outs
+    md = [-di / (beta * NBT) for di in d]     # normalize so EP grad is on BP scale
+    g_ep = grads_from(outs, md, list(QB.parameters()))
+    return g_bp, g_ep
+
+SEL = [i for i, n in enumerate(names) if n.endswith('attn.qkv.weight') or n.endswith('ff.w2.weight')]
+SEL = [i for i in SEL if (names[i].split('.')[0], names[i].split('.')[1]) in
+       (('0','attn'), ('8','attn'), ('6','ff'), ('11','ff'))]
+betas = [float(s) for s in a.betas.split(',')]
+REP = 6
+x, y = get_batch()
+g_bp_ref = None
+print(f'qbits={a.qbits} stochastic rounding | one batch, {REP} rounding draws per estimate', flush=True)
+for mode in ('frozen', 'split'):
+    store = {i: {b: [] for b in betas} for i in SEL}
+    for b in betas:
+        for r in range(REP):
+            g_bp, g_ep = ep_and_bp(x, y, b, mode)
+            if g_bp_ref is None: g_bp_ref = [g.detach().cpu() for g in g_bp]
+            for i in SEL: store[i][b].append(g_ep[i].detach().cpu())
+            del g_bp, g_ep
+            torch.cuda.empty_cache()
+    print(f'--- mode={mode} (frozen=trainer-faithful one draw/estimate; split=decorrelated free vs nudged) ---', flush=True)
+    for i in SEL:
+        m, n = params[i].shape
+        gb = g_bp_ref[i]
+        s1 = float(torch.linalg.svdvals(gb)[0])
+        u1 = torch.linalg.svd(gb, full_matrices=False).U[:, 0]
+        cells, Rs = [], []
+        for b in betas:
+            Xi = torch.stack(store[i][b])
+            Xif = Xi - Xi.mean(0, keepdim=True)
+            edge = float(np.mean([torch.linalg.svdvals(Xif[j])[0] for j in range(Xif.shape[0])]))
+            R = s1 / max(edge, 1e-30)
+            ovs = [abs(float(u1 @ torch.linalg.svd(ge, full_matrices=False).U[:, 0])) for ge in store[i][b]]
+            Rs.append(R); cells.append(f'{R:8.2f}({np.mean(ovs):.3f})')
+        lb = np.log(np.array(betas)); lR = np.log(np.maximum(Rs, 1e-12))
+        tag = 'no-crossing(R<1 all)' if max(Rs) < 1 else ('no-crossing(R>1 all)' if min(Rs) > 1 else '')
+        if not tag:
+            for j in range(len(betas) - 1):
+                if (lR[j]) * (lR[j + 1]) <= 0 and lR[j] != lR[j + 1]:
+                    t = (0.0 - lR[j]) / (lR[j + 1] - lR[j]); tag = f'beta*={float(np.exp(lb[j] + t * (lb[j + 1] - lb[j]))):.2e}'; break
+        print(f'{names[i]:22s} {m:5d}x{n:<5d} s1 {s1:9.3g}  ' + '  '.join(cells) + f'   {tag}', flush=True)
+print('BBP3_DONE', flush=True)