RESULT 47+48: 涨的量=blocks8-11单模态增益(ρ复刻探针,两血统同构型90%质量,plain 205-210k过1,cent全程平)+b8 logit跑飞70→94; BBP审计=fp32模拟器无加性下界(a=0,overlap 1.000@3e-4)→BBP是硬件设计方程; 三探针入库

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

4 files changed, 514 insertions, 0 deletions

diff --git a/ep_run/probe_bbp.py b/ep_run/probe_bbp.py
new file mode 100644
index 0000000..053f5e5
--- /dev/null
+++ b/ep_run/probe_bbp.py
@@ -0,0 +1,176 @@
+"""BBP floor: is wall-1 a spiked-matrix detectability transition? 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='3e-4,1e-3,3e-3,1e-2')
+ap.add_argument('--nb', 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())
+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):
+    obj = sum((o * c.detach()).sum() for o, c in zip(outs_list, cots))
+    gs = torch.autograd.grad(obj, params, 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, params)]
+
+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 blocks:
+        i = prev.detach().requires_grad_(True)
+        o = b(i)
+        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)
+            o = blocks[l](i)
+            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)
+    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')]
+betas = [float(s) for s in a.betas.split(',')]
+batches = [get_batch() for _ in range(a.nb)]
+# collect per layer: err entry-std per beta (across batches), sigma1(gbp), top-vec overlap per beta
+acc = {i: {b: {'s': [], 'ov': []} for b in betas} for i in SEL}
+sig1 = {i: [] for i in SEL}
+for (x, y) in batches:
+    ref = None
+    for b in betas:
+        g_bp, g_ep = ep_and_bp(x, y, b)
+        for i in SEL:
+            gb, ge = g_bp[i], g_ep[i]
+            if ref is None: pass
+            err = ge - gb
+            acc[i][b]['s'].append(float(err.std()))
+            try:
+                ub = torch.linalg.svd(gb, full_matrices=False).U[:, 0]
+                ue = torch.linalg.svd(ge, full_matrices=False).U[:, 0]
+                acc[i][b]['ov'].append(abs(float(ub @ ue)))
+            except Exception:
+                acc[i][b]['ov'].append(float('nan'))
+        for i in SEL:
+            if b == betas[0]: sig1[i].append(float(torch.linalg.svdvals(g_bp[i])[0]))
+        del g_bp, g_ep
+        torch.cuda.empty_cache()
+print('layer                    m x n      sigma1(gBP)   a(add)      b(mult)   beta*_edge   ov@' +
+      ' ov@'.join(f'{b:g}' for b in betas), flush=True)
+for i in SEL:
+    m, n = params[i].shape
+    s1 = float(np.mean(sig1[i]))
+    ss = np.array([np.mean(acc[i][b]['s']) for b in betas])
+    X = np.vstack([1.0 / np.array(betas), np.ones(len(betas))]).T
+    coef, *_ = np.linalg.lstsq(X, ss, rcond=None)
+    aa, bb = max(coef[0], 0.0), max(coef[1], 0.0)
+    edge = aa * (np.sqrt(m) + np.sqrt(n))
+    bstar = edge / max(s1, 1e-30)
+    ovs = ' '.join(f'{np.nanmean(acc[i][b]["ov"]):.3f}' for b in betas)
+    print(f'{names[i]:22s} {m:5d}x{n:<5d} {s1:11.4g} {aa:10.3g} {bb:10.3g} {bstar:11.3g}   {ovs}', flush=True)
+print('BBP_DONE', flush=True)
diff --git a/ep_run/probe_rhorelax.py b/ep_run/probe_rhorelax.py
new file mode 100644
index 0000000..20acbe5
--- /dev/null
+++ b/ep_run/probe_rhorelax.py
@@ -0,0 +1,144 @@
+"""Replicate the governor's rho meter offline: trainer-faithful NUDGED relax (derive-d via
+top CE grad + down-chain vjp, rebuild via up-chain feedforward + d) at fixed ckpts, K sweeps,
+fixed val batch. Reports per-sweep residual ratio rho (the exact quantity beta_cap gates on)
++ the per-block residual profile of the dominant mode (localization), per lineage x step.
+This is the operator whose contraction collapse killed crown-3; free-state norm audits
+(specaudit, statej) could not see it — the growth may live in curvature/alignment."""
+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('--ckpts', default=('plain:150000,plain:185000,plain:195000,plain:200000,'
+                                    'plain:210000,plain:230000,cent:150000,cent:185000,'
+                                    'cent:195000,cent:200000,cent:210000,cent:230000'))
+ap.add_argument('--K', type=int, default=30)
+ap.add_argument('--beta', type=float, default=3e-3)
+a = ap.parse_args()
+dev = 'cuda'
+torch.manual_seed(7)
+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))
+
+first = True
+for spec in a.ckpts.split(','):
+    lineage, step = spec.split(':'); step = int(step)
+    p = f'runs/fw72m_{lineage}_s{step}.pt'
+    try:
+        ck = torch.load(p, map_location='cpu', weights_only=False)
+    except FileNotFoundError:
+        print(f'{lineage} s{step}: MISSING', flush=True); continue
+    cfg = ck['config']; C, H, L = cfg['C'], cfg['H'], cfg['L']
+    if first:
+        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')
+        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]).to(dev)
+        y = torch.stack([torch.from_numpy(data[i + 1:i + 1 + T].astype(np.int64)) for i in ix]).to(dev)
+        first = False
+    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
+
+    def readout(z): return ln_f(z) @ W_out.t()
+
+    # free feedforward pass (block-wise graphs, EP-style detach between blocks)
+    f_ins, f_outs = [], []
+    prev = tok(x).detach()
+    for b in blocks:
+        i = prev.detach().requires_grad_(True)
+        o = b(i)
+        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
+    rhos, res_list = [], []
+    prof = None
+    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] = (-a.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 = [], []
+        rnum, rden, pblk = 0.0, 0.0, []
+        for l in range(L):
+            i = prev.detach().requires_grad_(True)
+            o = blocks[l](i)
+            n_ins.append(i); n_outs.append(o)
+            znew = o.detach().float() + d[l]
+            dn = float((znew - zs[l]).norm())
+            rnum += dn; rden += float(zs[l].norm()); pblk.append(dn)
+            zs[l] = znew
+            prev = zs[l]
+        ins, outs = n_ins, n_outs
+        res = rnum / max(rden, 1e-9)
+        if res_list: rhos.append(res / max(res_list[-1], 1e-12))
+        res_list.append(res)
+        prof = pblk
+        if not np.isfinite(res) or res > 1e3:
+            print(f'{lineage} s{step//1000}k: DIVERGED at sweep {k} (res {res:.2e})', flush=True)
+            break
+    tail = rhos[-5:] if len(rhos) >= 5 else rhos
+    pn = np.array(prof) / (np.sum(prof) + 1e-30)
+    print(f'{lineage} s{step//1000}k | res0 {res_list[0]:.4f} resK {res_list[-1]:.2e} | '
+          f'rho tail-med {np.median(tail):.4f} max {max(rhos):.4f} | mode blk-profile ' +
+          ' '.join(f'{v:.2f}' for v in pn), flush=True)
+    del tok, blocks, W_out, ln_f, ins, outs, f_ins, f_outs, zs, d
+    torch.cuda.empty_cache()
+print('RHORELAX_DONE', flush=True)
diff --git a/ep_run/probe_specaudit.py b/ep_run/probe_specaudit.py
new file mode 100644
index 0000000..041d068
--- /dev/null
+++ b/ep_run/probe_specaudit.py
@@ -0,0 +1,64 @@
+"""Weight-space spectral audit: which loop-gain component separates the plain vs cent
+lineages before the 196k ceiling collapse? Per-block sigma of Wq/Wk/Wv/proj/w1/w3/w2 and
+rms of qk-norm & post-norm gains, at matched steps. Pure CPU, reads the 5k ckpt ladders."""
+import torch, sys
+
+STEPS_BOTH = [150000, 170000, 185000, 190000, 195000, 200000]
+STEPS_PLAIN_EXTRA = [210000, 220000, 230000]
+L = 12
+
+def sig(w):
+    return float(torch.linalg.svdvals(w.float())[0])
+
+def rms(g):
+    return float(g.float().pow(2).mean().sqrt())
+
+def audit(path):
+    ck = torch.load(path, map_location='cpu', weights_only=False)
+    b = ck['blocks']
+    rows = []
+    for l in range(L):
+        qkv = b[f'{l}.attn.qkv.weight'].float()
+        wq, wk, wv = qkv[0:512], qkv[512:1024], qkv[1024:1536]
+        rows.append(dict(
+            q=sig(wq), k=sig(wk), v=sig(wv), pr=sig(b[f'{l}.attn.proj.weight']),
+            w1=sig(b[f'{l}.ff.w1.weight']), w3=sig(b[f'{l}.ff.w3.weight']),
+            w2=sig(b[f'{l}.ff.w2.weight']),
+            qg=rms(b[f'{l}.attn.qn.g']), kg=rms(b[f'{l}.attn.kn.g']),
+            na=rms(b[f'{l}.na.g']), nf=rms(b[f'{l}.nf.g'])))
+    wout = sig(ck['wout'])
+    return rows, wout
+
+def summarize(tag, step, rows, wout):
+    # aggregates: mean over blocks, top-half mean (6-11), max block, plus the derived products
+    def agg(key, blks):
+        vals = [rows[l][key] for l in blks]
+        return sum(vals) / len(vals)
+    bot, top = range(0, 6), range(6, 12)
+    logit = [rows[l]['qg'] * rows[l]['kg'] for l in range(L)]           # qk-norm logit scale
+    vpath = [rows[l]['v'] * rows[l]['pr'] for l in range(L)]            # attn value-path gain
+    fpath = [max(rows[l]['w1'], rows[l]['w3']) * rows[l]['w2'] for l in range(L)]
+    print(f'{tag} s{step//1000:>3}k | logit bot {agg("qg",bot)*agg("kg",bot):6.3f} top {agg("qg",top)*agg("kg",top):6.3f} max {max(logit):6.3f}(b{logit.index(max(logit))}) '
+          f'| vpath top {sum(vpath[6:])/6:7.2f} max {max(vpath):7.2f}(b{vpath.index(max(vpath))}) '
+          f'| ffn top {sum(fpath[6:])/6:7.2f} max {max(fpath):7.2f}(b{fpath.index(max(fpath))}) '
+          f'| na top {agg("na",top):5.3f} nf top {agg("nf",top):5.3f} | wout {wout:6.1f}', flush=True)
+    return dict(logit=logit, vpath=vpath, fpath=fpath)
+
+res = {}
+for lineage, steps in [('plain', STEPS_BOTH + STEPS_PLAIN_EXTRA), ('cent', STEPS_BOTH + STEPS_PLAIN_EXTRA)]:
+    for s in steps:
+        p = f'runs/fw72m_{lineage}_s{s}.pt'
+        try:
+            rows, wout = audit(p)
+        except FileNotFoundError:
+            print(f'{lineage} s{s}: MISSING', flush=True); continue
+        res[(lineage, s)] = summarize(lineage, s, rows, wout)
+
+# per-block divergence table at 195k (last matched healthy step)
+if ('plain', 195000) in res and ('cent', 195000) in res:
+    print('\nper-block plain/cent ratio at 195k (the pre-collapse fingerprint):', flush=True)
+    p, c = res[('plain', 195000)], res[('cent', 195000)]
+    print('blk  logit-ratio  vpath-ratio  ffn-ratio')
+    for l in range(L):
+        print(f'{l:3d}  {p["logit"][l]/c["logit"][l]:10.3f}  {p["vpath"][l]/c["vpath"][l]:10.3f}  {p["fpath"][l]/c["fpath"][l]:9.3f}')
+print('AUDIT_DONE', flush=True)
diff --git a/ep_run/probe_statej.py b/ep_run/probe_statej.py
new file mode 100644
index 0000000..f25325b
--- /dev/null
+++ b/ep_run/probe_statej.py
@@ -0,0 +1,130 @@
+"""State-side Jacobian audit: per-block sigma(dO_l/dz_in) at the FREE operating state,
+both lineages, matched steps. Weight-space audit (probe_specaudit) showed all weight scales
+FLAT while the beta-normalized loop gain G=0.9/beta_cap grew ~300x -> hypothesis: the growth
+is the CHAIN PRODUCT of per-block state Jacobians (each +tens-of-%, ^12 = hundreds-x).
+Also dumps per-block max|logit| at the state (softcap calibration).
+FD-JVP + autograd-vjp power iteration on J^T J (no forward-mode through SDPA)."""
+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('--ckpts', default=('plain:150000,plain:185000,plain:195000,plain:200000,'
+                                    'plain:210000,plain:230000,cent:150000,cent:185000,'
+                                    'cent:195000,cent:200000,cent:210000,cent:230000'))
+ap.add_argument('--iters', type=int, default=25)
+a = ap.parse_args()
+dev = 'cuda'
+torch.manual_seed(7)
+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 logits_stats(self, x):
+        Bn, Tn, C = x.shape
+        q, k, _ = 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))
+        lg = (q @ k.transpose(-2, -1)) / (self.hd ** 0.5)
+        mask = torch.ones(Tn, Tn, dtype=torch.bool, device=x.device).tril()
+        lg = lg.masked_fill(~mask, 0.0)
+        return float(lg.abs().max()), float(lg.abs().quantile(0.999))
+    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))
+
+def sigma_J(block, z0):
+    """largest singular value of dblock(z)/dz at z0: power iteration on J^T J.
+    Jv by central FD (fp32, scaled eps), J^T u by autograd."""
+    z0 = z0.detach()
+    v = torch.randn_like(z0); v /= v.norm()
+    s = None
+    for _ in range(a.iters):
+        eps = 1e-3 * z0.norm() / (v.norm() * (z0.numel() ** 0.5) + 1e-30) * (z0.numel() ** 0.5)
+        with torch.no_grad():
+            jv = (block(z0 + eps * v) - block(z0 - eps * v)) / (2 * eps)
+        zg = z0.clone().requires_grad_(True)
+        o = block(zg)
+        jtu = torch.autograd.grad((o * jv.detach()).sum(), zg)[0]
+        s = float(jtu.norm().sqrt())          # |J^T J v|^{1/2} -> sigma as v converges
+        v = jtu / (jtu.norm() + 1e-30)
+    return s
+
+first = True
+for spec in a.ckpts.split(','):
+    lineage, step = spec.split(':'); step = int(step)
+    p = f'runs/fw72m_{lineage}_s{step}.pt'
+    try:
+        ck = torch.load(p, map_location='cpu', weights_only=False)
+    except FileNotFoundError:
+        print(f'{lineage} s{step}: MISSING', flush=True); continue
+    cfg = ck['config']; C, H, L = cfg['C'], cfg['H'], cfg['L']
+    if first:
+        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')
+        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]).to(dev)
+        first = False
+    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)
+    with torch.no_grad():
+        zs = [tok(x)]
+        for b in blocks: zs.append(b(zs[-1]))
+    sigs, lgs = [], []
+    for l in range(L):
+        sigs.append(sigma_J(blocks[l], zs[l]))
+        lgs.append(blocks[l].attn.logits_stats(zs[l]))
+    prod = float(np.prod(sigs))
+    top = float(np.prod(sigs[6:]))
+    print(f'{lineage} s{step//1000}k | sig_J per blk ' + ' '.join(f'{s:5.2f}' for s in sigs) +
+          f' | PROD {prod:9.1f} top6 {top:7.1f} | max|logit| ' +
+          ' '.join(f'{m:4.1f}' for m, _ in lgs), flush=True)
+    del tok, blocks, zs
+    torch.cuda.empty_cache()
+print('STATEJ_DONE', flush=True)