path: root/ep_run/casc_eq_train.py

"""Cascade-EP trainer — EQUILIBRIUM MODE (the true-EP route).
Two-phase (+-beta) relaxation of all layer states to the nudged equilibria via
Gauss-Seidel reverse sweeps (solver choice only; readout is taken AT the relaxed
states with the standard EP formula), weight grad = (1/2beta)[dF/dtheta|+ - dF/dtheta|-].
Inference = plain forward (standard LLM). Twin of casc_bp_train.py (same seed/data)."""
import argparse, math, pickle, time
import numpy as np, torch, torch.nn as nn, torch.nn.functional as F
from pathlib import Path

ap = argparse.ArgumentParser()
ap.add_argument('--tag', default='casc_eq6')
ap.add_argument('--L', type=int, default=6); ap.add_argument('--C', type=int, default=256)
ap.add_argument('--H', type=int, default=8); ap.add_argument('--T', type=int, default=256)
ap.add_argument('--B', type=int, default=24); ap.add_argument('--steps', type=int, default=4000)
ap.add_argument('--lr', type=float, default=3e-4); ap.add_argument('--warmup', type=int, default=200)
ap.add_argument('--beta', type=float, default=0.003); ap.add_argument('--seed', type=int, default=0)
ap.add_argument('--K', type=int, default=3)             # fb (message-passing) rounds
ap.add_argument('--geta', type=float, default=1.0)      # fb mixing (1.0 = undamped)
ap.add_argument('--save_every', type=int, default=1000); ap.add_argument('--log', type=int, default=100)
ap.add_argument('--wandb', default=''); ap.add_argument('--wandb_run', default='')
ap.add_argument('--kmax', type=int, default=8)          # adaptive fb rounds cap
ap.add_argument('--noguard', action='store_true')       # diagnosis: skip only non-finite grads
ap.add_argument('--untie', action='store_true')         # separate readout matrix (untied from tok)
ap.add_argument('--opt', choices=['adamw', 'muon'], default='adamw')
ap.add_argument('--muon_lr', type=float, default=0.02)
ap.add_argument('--tok_init', type=float, default=0.0)  # >0: init tok/pos with this std (GPT-standard 0.02)
ap.add_argument('--compile', action='store_true')       # torch.compile each block (free speed where supported)
ap.add_argument('--sig_every', type=int, default=25)    # tok-sigma refresh interval (amortized)
ap.add_argument('--beta_floor', type=float, default=0.0) # >0: floor beta_t (anti finite-beta SNR collapse at depth)
ap.add_argument('--beta_fixed', action='store_true')     # disable sig^2 schedule, hold beta_t = args.beta constant
ap.add_argument('--cosine', action='store_true')         # warmup then cosine decay to lr_min_ratio*lr over --steps (long runs)
ap.add_argument('--lr_min_ratio', type=float, default=0.1)
ap.add_argument('--qk_norm', action='store_true')        # RMS-norm q,k per head before scores (OLMo2-style; bounds logits, analog-friendly)
ap.add_argument('--final_ln', action='store_true')       # final LayerNorm before readout (standard GPT; bounds sig_tok growth -> keeps beta/estimator healthy on long runs)
ap.add_argument('--resume', default='')                  # path to a ckpt (tok/pos/blocks) to continue from; step taken from ckpt
ap.add_argument('--sig0', type=float, default=-1.0)      # override SIG0 (beta-schedule ref); needed on resume to restore original beta regime
ap.add_argument('--olmo2', action='store_true')          # OLMo2-standard block: norm-AFTER-sublayer RMSNorm, full-width QK-norm, RoPE(500k), SwiGLU, no-bias, untied head, final RMSNorm, 0.02 init
ap.add_argument('--wd', type=float, default=-1.0)        # >=0: grouped weight decay (linear weights+head decay; embeddings/norm-gains none). <0 = legacy uniform 1e-4
ap.add_argument('--zloss', type=float, default=0.0)      # z-loss coefficient on train objective (OLMo2-style logit regularizer); 0 = off
ap.add_argument('--kretry', type=int, default=0)         # >0: on drift-reject, RETRY the batch once with this many fb rounds (diag B: K8 converges the marginal batches) instead of dropping it
ap.add_argument('--bf_late', type=float, default=0.0)    # >0: raise beta_floor to this value from step --bf_late_at (late-training SNR fix; dose-response 2026-07-10)
ap.add_argument('--bf_late_at', type=int, default=25000)
ap.add_argument('--bsign_rand', action='store_true')  # random-sign beta per step (KHS 'random scheme'): averages the O(beta) single-sided bias at single-phase cost
ap.add_argument('--bf16', action='store_true')           # cast model to bf16 (E-accumulation + tok_sigma stay fp32) — the x0.5 cost lever, GATE before production
ap.add_argument('--amp', action='store_true')            # PROPER mixed precision: autocast(bf16) matmuls, fp32 params/states/d/E — amp_gate.py PASSED 2026-07-12 (cos 0.9682 vs fp32 0.9687); --bf16 naive-cast stays DEAD (state quantization, RESULT 11)
ap.add_argument('--dtop_every', type=int, default=1)    # 1 = exact (DEFAULT, BP-parity); 2 = fast mode (~20% cheaper, ~4% CE tax at high lr)
ap.add_argument('--gate_every', type=int, default=200)  # in-training cos(EP,BP) telemetry; <=0 = fully BP-free (no bp_gate at all)
ap.add_argument('--gate_govern', action='store_true')   # let gate cos adjust K/bscale (default: observe-only => training control is BP-free)
args = ap.parse_args()
if args.olmo2:
    args.untie = True
    if args.tok_init <= 0: args.tok_init = 0.02
torch.manual_seed(args.seed)
dev = 'cuda' if torch.cuda.is_available() else 'cpu'

DD = Path('/home/yurenh2/ept/ep_run/data/tinystories_bpe')
vocab = pickle.load(open(DD / 'meta.pkl', 'rb'))['vocab_size']

def get_batch(split):
    data = np.memmap(DD / ('train.bin' if split == 'train' else 'val.bin'), dtype=np.uint16, mode='r')
    ix = torch.randint(len(data) - args.T - 1, (args.B,))
    x = torch.stack([torch.from_numpy(data[i:i + args.T].astype(np.int64)) for i in ix])
    y = torch.stack([torch.from_numpy(data[i + 1:i + 1 + args.T].astype(np.int64)) for i in ix])
    return x.to(dev), y.to(dev)

class CausalSelfAttn(nn.Module):
    """explicit MHA (SDPA-backed) so we can QK-norm q,k per head before the scores."""
    def __init__(self, C, H, qk_norm=False):
        super().__init__()
        self.H, self.hd, self.qk_norm = H, C // H, qk_norm
        self.qkv = nn.Linear(C, 3 * C)
        self.proj = nn.Linear(C, C)
        if qk_norm:
            self.q_g = nn.Parameter(torch.ones(self.hd))
            self.k_g = nn.Parameter(torch.ones(self.hd))
    def forward(self, x):
        B, T, C = x.shape
        q, k, v = self.qkv(x).split(C, dim=2)
        q = q.view(B, T, self.H, self.hd).transpose(1, 2)
        k = k.view(B, T, self.H, self.hd).transpose(1, 2)
        v = v.view(B, T, self.H, self.hd).transpose(1, 2)
        if self.qk_norm:  # RMS-norm over head_dim (OLMo2-style), learnable per-dim gain
            q = q * torch.rsqrt(q.pow(2).mean(-1, keepdim=True) + 1e-6) * self.q_g
            k = k * torch.rsqrt(k.pow(2).mean(-1, keepdim=True) + 1e-6) * self.k_g
        y = F.scaled_dot_product_attention(q, k, v, is_causal=True)
        return self.proj(y.transpose(1, 2).contiguous().view(B, T, C))

class Block(nn.Module):
    def __init__(self, C, H, qk_norm=False):
        super().__init__()
        self.ln1, self.ln2 = nn.LayerNorm(C), nn.LayerNorm(C)
        self.attn = CausalSelfAttn(C, H, qk_norm)
        self.ff = nn.Sequential(nn.Linear(C, 4 * C), nn.GELU(), nn.Linear(4 * C, C))
    def forward(self, z, mask=None):
        z = z + self.attn(self.ln1(z))
        return z + self.ff(self.ln2(z))

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   # ~param-match the 4x-GELU MLP (8C^2)
        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):
    """OLMo2 attention: no-bias projs, FULL-WIDTH RMS QK-norm (pre-head-split, HF Olmo2 order), then per-head RoPE."""
    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):
        x1, x2 = x[..., ::2], x[..., 1::2]
        c, s = self.rc[None, None], self.rs[None, None]
        return torch.stack((x1 * c - x2 * s, x1 * s + x2 * c), dim=-1).flatten(-2)
    def forward(self, x):
        B, T, 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(B, T, self.H, self.hd).transpose(1, 2))
        k = self.rope(k.view(B, T, self.H, self.hd).transpose(1, 2))
        v = v.view(B, T, 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(B, T, C))

class Olmo2Block(nn.Module):
    """OLMo2 reordered norm (norm AFTER each sublayer, inside the residual) — their training-stability change."""
    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, mask=None):
        z = z + self.na(self.attn(z))
        return z + self.nf(self.ff(z))

tok = nn.Embedding(vocab, args.C).to(dev)
pos = nn.Embedding(args.T, args.C).to(dev)
if args.tok_init > 0:
    with torch.no_grad():
        tok.weight.normal_(0, args.tok_init); pos.weight.normal_(0, args.tok_init)
blocks = nn.ModuleList([(Olmo2Block(args.C, args.H, args.T) if args.olmo2 else Block(args.C, args.H, args.qk_norm)) for _ in range(args.L)]).to(dev)
if args.olmo2:
    with torch.no_grad():
        for m in blocks.modules():
            if isinstance(m, nn.Linear): m.weight.normal_(0, 0.02)
if args.compile:
    try:
        for i in range(args.L): blocks[i] = torch.compile(blocks[i], mode='reduce-overhead')
        print('[compile] blocks compiled', flush=True)
    except Exception as e:
        print(f'[compile] disabled ({e})', flush=True)
mask = torch.triu(torch.full((args.T, args.T), float('-inf'), device=dev), 1)
W_out = nn.Parameter(torch.randn(vocab, args.C, device=dev) * 0.02) if args.untie else None
ln_f = (RMSNorm(args.C) if args.olmo2 else (nn.LayerNorm(args.C) if args.final_ln else nn.Identity())).to(dev)
def emb(x):
    return tok(x) if args.olmo2 else tok(x) + pos(torch.arange(args.T, device=dev))[None]
readout = (lambda z: ln_f(z) @ W_out.t()) if args.untie else (lambda z: ln_f(z) @ tok.weight.t())
all_params = list(tok.parameters()) + ([] if args.olmo2 else list(pos.parameters())) + list(blocks.parameters()) + list(ln_f.parameters()) + ([W_out] if args.untie else [])
start_step = 0
if args.resume:
    _ck = torch.load(args.resume, map_location=dev, weights_only=False)
    tok.load_state_dict(_ck['tok']); pos.load_state_dict(_ck['pos']); blocks.load_state_dict(_ck['blocks'])
    if _ck.get('wout') is not None and args.untie:
        with torch.no_grad(): W_out.copy_(_ck['wout'].to(dev))
    if _ck.get('lnf') is not None and not isinstance(ln_f, nn.Identity): ln_f.load_state_dict(_ck['lnf'])
    start_step = int(_ck.get('step', 0))
    print(f'[resume] loaded {args.resume} at step {start_step}', flush=True)
if args.bf16:
    for _m in (tok, pos, blocks):
        _m.to(torch.bfloat16)
    if not isinstance(ln_f, nn.Identity): ln_f.to(torch.bfloat16)
    if args.untie:
        with torch.no_grad(): W_out.data = W_out.data.to(torch.bfloat16)
    print('[bf16] model cast to bfloat16 (E-accum + sigma stay fp32)', flush=True)
if args.opt == 'muon':
    from muon import build_hybrid
    opt, sched = build_hybrid(blocks, all_params, args.lr, args.muon_lr, args.warmup,
                              total_steps=(args.steps if args.cosine else 0), lr_min_ratio=args.lr_min_ratio)
else:
    if args.wd >= 0:   # OLMo2-style grouped decay: linear weights + head decay; embeddings/norm-gains none
        nodecay = {id(p) for p in tok.parameters()} | {id(p) for p in pos.parameters()} | \
                  {id(p) for p in blocks.parameters() if p.ndim < 2} | {id(p) for p in ln_f.parameters()}
        opt = torch.optim.AdamW([
            {'params': [p for p in all_params if id(p) not in nodecay], 'weight_decay': args.wd},
            {'params': [p for p in all_params if id(p) in nodecay], 'weight_decay': 0.0}], lr=args.lr)
    else:
        opt = torch.optim.AdamW(all_params, lr=args.lr, weight_decay=1e-4)
    if args.cosine:
        def _lrlam(s):
            if s < args.warmup: return (s + 1) / max(args.warmup, 1)
            p = min(1.0, (s - args.warmup) / max(1, args.steps - args.warmup))
            return args.lr_min_ratio + 0.5 * (1 - args.lr_min_ratio) * (1 + math.cos(math.pi * p))
        sched = torch.optim.lr_scheduler.LambdaLR(opt, _lrlam)
    else:
        sched = torch.optim.lr_scheduler.LambdaLR(opt, lambda s: min(1.0, (s + 1) / max(args.warmup, 1)))
NBT = args.B * args.T

def obj_loss(logits2d, y1d):
    """train objective: CE (+ optional z-loss). Used in the nudge force, theta-readout and bp_gate
    so EP tracks BP on the SAME objective; evaluate() stays pure CE for comparability."""
    l = F.cross_entropy(logits2d, y1d)
    if args.zloss > 0:
        l = l + args.zloss * (torch.logsumexp(logits2d.float(), -1) ** 2).mean()
    return l

def free_states_graphed(x):
    """free forward, keeping per-layer graphs (in_l, out_l) so round-1 backward vjps reuse them."""
    with torch.no_grad():
        z0 = emb(x)
    ins, outs, zs = [], [], []
    prev = z0
    with torch.autocast('cuda', dtype=torch.bfloat16, enabled=args.amp):
        for b in blocks:
            i = prev.detach().requires_grad_(True)
            o = b(i, mask)
            ins.append(i); outs.append(o); zs.append(o.detach().float())
            prev = zs[-1]
    return z0, zs, ins, outs

@torch.no_grad()
def tok_sigma(iters=8):
    """top singular value of tok.weight (power iteration on the raw matrix)."""
    W = (W_out if args.untie else tok.weight).float()
    v = torch.randn(W.shape[1], device=dev); v /= v.norm()
    sig = 1.0
    for _ in range(iters):
        u = W @ v; u /= max(u.norm(), 1e-12)
        v = W.t() @ u; sig = v.norm(); v /= max(sig, 1e-12)
    return float(sig)

def relax(z0, zs, ins, outs, y, beta, K, x):
    """K fb rounds with GRAPH REUSE + two dedups: (a) the top CE force d_top is refreshed on
    even rounds only (states move O(beta) per round -> O(beta^2) error); (b) the LAST rebuild
    keeps graphs (layer-0 fed a graphed emb) and returns (ins, outs) so the theta-readout
    reuses them instead of re-running a full graphed chain."""
    d = [None] * args.L
    for k in range(K):
        if k % args.dtop_every == 0 or d[args.L - 1] is None:
            zc = zs[args.L - 1].detach().requires_grad_(True)
            ce = obj_loss(readout(zc).reshape(-1, vocab), y.reshape(-1))
            d[args.L - 1] = (-beta * NBT * torch.autograd.grad(ce, zc)[0]).detach()
        for l in range(args.L - 2, -1, -1):
            d[l] = torch.autograd.grad(outs[l + 1], ins[l + 1], grad_outputs=d[l + 1].to(outs[l + 1].dtype))[0].detach().float()
        last = (k + 1 == K)
        prev = z0
        n_ins, n_outs = [], []
        with torch.autocast('cuda', dtype=torch.bfloat16, enabled=args.amp):
            for l in range(args.L):
                if last and l == 0:
                    i = emb(x)   # graphed emb for the readout's E-path
                else:
                    i = prev.detach().requires_grad_(True)
                o = blocks[l](i, mask)
                zs[l] = (o.detach().float() + d[l])
                n_ins.append(i); n_outs.append(o)
                prev = zs[l]
        ins, outs = n_ins, n_outs
    return zs, outs

def dFdtheta(zs, x, y, beta):
    """theta-readout at FIXED states. Not used by the training loop (relax reuses its own
    graphs); kept as the INVARIANT-TEST surface for test_bp_free.py. Self-sealing: inputs
    are detached here so the local-graph property holds for any caller."""
    zs = [z.detach() for z in zs]
    prev = emb(x)
    E = 0.0
    for z, b in zip(zs, blocks):
        E = E + 0.5 * ((z - b(prev, mask)) ** 2).sum()
        prev = z   # zs detached at entry => blocks l>0 get detached inputs; block 0 gets the graphed emb
    obj = E / NBT + beta * F.cross_entropy(readout(zs[-1]).reshape(-1, vocab), y.reshape(-1))
    gs = torch.autograd.grad(obj, all_params, allow_unused=True)
    return [g if g is not None else None for g in gs]


SIG0 = None
BGEN = torch.Generator().manual_seed(args.seed + 990)   # separate RNG: sign flips must not shift the data stream
GOV = {'K': None, 'bscale': 1.0, 'gema': None, 'drift': 0.0, 'gn': 0.0, 'sig': 0.0}
def ep_step(x, y):
    """single-sided EP with a QUALITY-GOVERNED estimator: beta_t = beta0*bscale*sig0^2/sig^2,
    K = GOV['K'] fb rounds; guard = finiteness + drift + grad-norm sanity only."""
    global SIG0
    if GOV['K'] is None: GOV['K'] = args.K
    if GOV.get('step', 0) % args.sig_every == 0 or GOV.get('sig', 0) == 0:
        GOV['sig'] = tok_sigma()
    GOV['step'] = GOV.get('step', 0) + 1
    sig = GOV['sig']
    if SIG0 is None: SIG0 = args.sig0 if args.sig0 > 0 else sig
    beta_t = args.beta * GOV['bscale'] * (SIG0 * SIG0) / max(sig * sig, 1e-9)
    if args.beta_fixed: beta_t = args.beta * GOV['bscale']
    fl = args.beta_floor
    if args.bf_late > 0.0 and GOV.get('step', 0) >= args.bf_late_at: fl = args.bf_late
    if fl > 0.0: beta_t = max(beta_t, fl)
    if args.bsign_rand and torch.rand((), generator=BGEN).item() < 0.5: beta_t = -beta_t
    z0, zs, ins, outs = free_states_graphed(x)
    zs_free = [z.clone() for z in zs]
    free_ce = F.cross_entropy(readout(zs_free[-1]).reshape(-1, vocab), y.reshape(-1)).item()
    zp, last_outs = relax(z0, zs, ins, outs, y, +beta_t, GOV['K'], x)
    with torch.no_grad():
        drift = sum(float((a - b).norm()) for a, b in zip(zp, zs_free)) / max(
            sum(float(b.norm()) for b in zs_free), 1e-9)
    if (not math.isfinite(drift)) or (drift > 0.5 and not args.noguard):
        ok_retry = False
        if args.kretry > 0 and math.isfinite(drift) and not args.noguard:
            GOV['skr'] = GOV.get('skr', 0) + 1   # marginal batch: retry once with deeper relaxation
            z0, zs, ins, outs = free_states_graphed(x)
            zs_free = [z.clone() for z in zs]
            zp, last_outs = relax(z0, zs, ins, outs, y, +beta_t, args.kretry, x)
            with torch.no_grad():
                drift = sum(float((a - b).norm()) for a, b in zip(zp, zs_free)) / max(
                    sum(float(b.norm()) for b in zs_free), 1e-9)
            ok_retry = math.isfinite(drift) and drift <= 0.5
        if not ok_retry:
            GOV['skd'] = GOV.get('skd', 0) + 1   # drift-guard reject (relaxation non-convergence)
            for p in all_params: p.grad = None
            return free_ce, beta_t, GOV['K'], False
    GOV['drift'] = drift
    E = 0.0
    for z, o in zip(zp, last_outs): E = E + 0.5 * ((z.detach().float() - o.float()) ** 2).sum()   # fp32 accumulation (bf16-safe; no-op in fp32)
    obj = E / (NBT * beta_t) + obj_loss(readout(zp[-1].detach()).reshape(-1, vocab), y.reshape(-1))
    gs = torch.autograd.grad(obj, all_params, allow_unused=True)
    gn = 0.0
    for g in gs:
        if g is not None: gn += float((g ** 2).sum())
    gn = gn ** 0.5
    if GOV['gema'] is None: GOV['gema'] = gn
    GOV['gema'] = 0.99 * GOV['gema'] + 0.01 * gn        # EMA always updates (frozen-ref bugfix)
    GOV['gn'] = gn
    if not math.isfinite(gn) or (gn > 8 * GOV['gema'] and not args.noguard):
        GOV['skg'] = GOV.get('skg', 0) + 1   # gn-EMA-guard reject (gradient-magnitude spike)
        for p in all_params: p.grad = None
        return free_ce, beta_t, GOV['K'], False
    for p, g in zip(all_params, gs):
        p.grad = g
    return free_ce, beta_t, GOV['K'], True

def bp_gate(x, y):
    """true BP grads for telemetry cos (called before opt.step; reads p.grad separately)."""
    z = emb(x)
    for b in blocks: z = b(z, mask)
    ce = obj_loss(readout(z).reshape(-1, vocab), y.reshape(-1))
    return list(torch.autograd.grad(ce, all_params, allow_unused=True))

@torch.no_grad()
def evaluate(nb=6):
    tot = 0.0
    for _ in range(nb):
        x, y = get_batch('val')
        z = emb(x)
        for b in blocks: z = b(z, mask)
        tot += F.cross_entropy(readout(z).reshape(-1, vocab), y.reshape(-1)).item()
    return tot / nb

wb = None
if args.wandb:
    try:
        import wandb as _w
        wb = _w.init(project=args.wandb, name=args.wandb_run or args.tag, id=args.wandb_run or args.tag,
                     resume='allow', config=vars(args))
    except Exception as e:
        print(f'[wandb] disabled ({e})', flush=True)

n = sum(p.numel() for p in all_params)
print(f'[{args.tag}] cascade-EP(EQUILIBRIUM/fb) L{args.L} C{args.C} T{args.T} beta={args.beta} '
      f'K={args.K} geta={args.geta} | {n/1e6:.2f}M | {dev}', flush=True)
best, t0 = 1e9, time.time()
skips = 0
for _ in range(start_step): sched.step()   # advance LR schedule to the resumed step
for step in range(start_step, args.steps + 1):
    x, y = get_batch('train')
    ce, beta_t, rounds, ok = ep_step(x, y)
    if not ok: skips += 1
    gcos = float('nan')
    if args.gate_every > 0 and step % args.gate_every == 0 and ok:
        gbp = bp_gate(x, y)
        num = den1 = den2 = 0.0
        for p, g in zip(all_params, gbp):
            if p.grad is None or g is None: continue
            num += float((p.grad * g).sum()); den1 += float((p.grad ** 2).sum()); den2 += float((g ** 2).sum())
        gcos = num / max((den1 ** 0.5) * (den2 ** 0.5), 1e-12)
        if args.gate_govern:                              # opt-in: BP-informed control flow
            if gcos < 0.97:
                GOV['K'] = min(GOV['K'] + 2, args.kmax); GOV['bscale'] = max(GOV['bscale'] * 0.7, 0.05)
            elif gcos > 0.995 and GOV['K'] > args.K:
                GOV['K'] -= 1; GOV['bscale'] = min(GOV['bscale'] * 1.05, 1.0)
    torch.nn.utils.clip_grad_norm_(all_params, 1.0)
    opt.step(); sched.step(); opt.zero_grad(set_to_none=True)
    if step % args.log == 0:
        val = evaluate(); best = min(best, val)
        gtag = '' if math.isnan(gcos) else f' cos={gcos:.4f}'
        print(f'step {step:5d}/{args.steps} | train {ce:.4f} val {val:.4f} (best {best:.4f}) '
              f'| beta={beta_t:.2e} K={rounds} skips={skips}(d{GOV.get("skd",0)}/g{GOV.get("skg",0)}/r{GOV.get("skr",0)}){gtag} '
              f'drift={GOV["drift"]:.3f} gn={GOV["gn"]:.2e} sig={GOV["sig"]:.1f} | {step/max(time.time()-t0,1e-9):.3f} it/s', flush=True)
        if wb is not None:
            try: wb.log({'train_ce': ce, 'val_ce': val, 'best': best, 'beta_t': beta_t,
                         'rounds': rounds, 'skips': skips, 'gate_cos': (None if math.isnan(gcos) else gcos)}, step=step)
            except Exception: pass
    if step % args.save_every == 0 and step > 0:
        torch.save({'tok': tok.state_dict(), 'pos': pos.state_dict(), 'blocks': blocks.state_dict(),
                    'wout': (W_out.detach().cpu() if args.untie else None),
                    'lnf': (ln_f.state_dict() if not isinstance(ln_f, nn.Identity) else None),
                    'step': step, 'val': best, 'config': vars(args)}, Path('runs') / f'{args.tag}_s{step}.pt')
print(f'[{args.tag}] DONE best val CE {best:.4f} (BP twin 2.9746; zil-diagnostic 3.3236)', flush=True)
if wb is not None:
    try: wb.summary['best_val_ce'] = best; wb.finish()
    except Exception: pass