path: root/experiments/minimal_aux_compression.py

"""
Phase 10A.7: Minimal Auxiliary Compression Ablation.

Core question: Can gain come from compressed/frozen/per-layer representations
instead of a full input-conditioned network?

8 branches from the same DFA checkpoint at t0=5:
1. continue_DFA                  — pure DFA baseline
2. blend_random_trainable        — standard Vec (10A.5/10A.6 reference)
3. blend_zero_target_trainable   — Vec trained with loss=||a_aux||^2 (from 10A.6)
4. blend_zero_target_normmatched — zero_target + blockwise norm matching to random_trainable
5. blend_perlayer_vector         — no network; per-block nn.Parameter v_l, broadcast over batch
6. blend_random_freeze_after_1   — random Vec, train 1 epoch then freeze
7. blend_random_freeze_after_5   — random Vec, train 5 epochs then freeze
8. blend_random_freeze_after_10  — random Vec, train 10 epochs then freeze
"""
import os
import sys
import json
import argparse
import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.utils.data import DataLoader
import torchvision
import torchvision.transforms as transforms
import copy

sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))

from models.residual_mlp import ResidualMLP
from models.value_net import SinusoidalTimeEmbed
from metrics.credit_metrics import cosine_similarity_batch, perturbation_correlation


# ---------------------------------------------------------------------------
# Auxiliary network architectures
# ---------------------------------------------------------------------------

class VectorCreditNet(nn.Module):
    """Standard Vec: takes (h, t, s) -> d_hidden credit vector."""
    def __init__(self, d_hidden, s_dim, time_embed_dim=32, hidden_dim=256, num_layers=3):
        super().__init__()
        self.ln = nn.LayerNorm(d_hidden)
        self.time_embed = SinusoidalTimeEmbed(time_embed_dim)
        input_dim = d_hidden + time_embed_dim + s_dim
        layers = []
        for i in range(num_layers):
            in_d = input_dim if i == 0 else hidden_dim
            layers.append(nn.Linear(in_d, hidden_dim))
            layers.append(nn.GELU())
        layers.append(nn.Linear(hidden_dim, d_hidden))
        self.net = nn.Sequential(*layers)

    def forward(self, h, t, s):
        return self.net(torch.cat([self.ln(h), self.time_embed(t), s], dim=-1))


class PerLayerVector(nn.Module):
    """No network: each block l has a trainable nn.Parameter v_l of shape (d_hidden,).
    All samples in a batch receive the same v_l (broadcast).
    forward(h, t, s, block_idx) returns v_l expanded to (batch, d_hidden).
    """
    def __init__(self, d_hidden, num_blocks):
        super().__init__()
        # Initialize with small random values
        self.vectors = nn.ParameterList(
            [nn.Parameter(torch.randn(d_hidden) * 0.01) for _ in range(num_blocks)]
        )
        self._block_idx = 0

    def set_block(self, l):
        self._block_idx = l

    def forward(self, h, t, s):
        batch = h.size(0)
        return self.vectors[self._block_idx].unsqueeze(0).expand(batch, -1)


# ---------------------------------------------------------------------------
# Data
# ---------------------------------------------------------------------------

def get_cifar10(batch_size=128):
    transform_train = transforms.Compose([
        transforms.RandomCrop(32, padding=4), transforms.RandomHorizontalFlip(),
        transforms.ToTensor(),
        transforms.Normalize((0.4914, 0.4822, 0.4465), (0.2470, 0.2435, 0.2616))])
    transform_test = transforms.Compose([
        transforms.ToTensor(),
        transforms.Normalize((0.4914, 0.4822, 0.4465), (0.2470, 0.2435, 0.2616))])
    trainset = torchvision.datasets.CIFAR10(
        root='./data', train=True, download=True, transform=transform_train)
    testset = torchvision.datasets.CIFAR10(
        root='./data', train=False, download=True, transform=transform_test)
    return (DataLoader(trainset, batch_size=batch_size, shuffle=True,
                       num_workers=4, pin_memory=True),
            DataLoader(testset, batch_size=batch_size, shuffle=False,
                       num_workers=4, pin_memory=True))


# ---------------------------------------------------------------------------
# Evaluation helpers
# ---------------------------------------------------------------------------

def evaluate(model, test_loader, device):
    model.eval(); c, t = 0, 0
    with torch.no_grad():
        for x, y in test_loader:
            x = x.view(x.size(0), -1).to(device); y = y.to(device)
            c += (model(x).argmax(1) == y).sum().item(); t += x.size(0)
    return c / t


def compute_diagnostics(model, aux_net, Bs, test_loader, device, credit_mode, alpha=0.75):
    """Compute mean Gamma (BP cosine) and mean rho (perturbation correlation)."""
    model.eval()
    if aux_net is not None:
        aux_net.eval()
    L = model.num_blocks

    for x, y in test_loader:
        x = x.view(x.size(0), -1).to(device); y = y.to(device); break
    batch = x.size(0)

    # BP pass for hidden gradients (offline eval only, not used for training)
    was_frozen = not next(model.parameters()).requires_grad
    if was_frozen:
        for p in model.parameters(): p.requires_grad_(True)
    model.zero_grad()
    lo, hbp = model(x, return_hidden=True)
    for l in range(L + 1): hbp[l].retain_grad()
    F.cross_entropy(lo, y).backward()
    bp = {l: hbp[l].grad.detach().clone() for l in range(L + 1)}
    if was_frozen:
        for p in model.parameters(): p.requires_grad_(False)

    with torch.no_grad():
        lo2, hi = model(x, return_hidden=True)
        eT = lo2.softmax(-1); eT[torch.arange(batch), y] -= 1; s = eT.detach()

    gammas, rhos = [], []
    for l in range(L):
        h_l = hi[l].detach()
        t_l = torch.full((batch,), l / L, device=device)

        if credit_mode == 'dfa':
            a_l = (s @ Bs[l].T).detach()
        elif credit_mode == 'blend' and aux_net is not None:
            a_dfa = (s @ Bs[l].T).detach()
            if isinstance(aux_net, PerLayerVector):
                aux_net.set_block(l)
            a_aux = aux_net(h_l, t_l, s).detach()
            rd = (a_dfa ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
            rv = (a_aux ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
            a_l = alpha * a_aux / rv + (1 - alpha) * a_dfa / rd
        else:
            a_l = (s @ Bs[l].T).detach()

        gammas.append(cosine_similarity_batch(a_l, bp[l]))

        def make_fwd(sl):
            def f(h):
                with torch.no_grad():
                    c = h
                    for i in range(sl, L):
                        c = c + model.blocks[i](c)
                    return F.cross_entropy(
                        model.out_head(model.out_ln(c)), y, reduction='none')
            return f

        rhos.append(perturbation_correlation(h_l, a_l, make_fwd(l), epsilon=1e-3, M=16))

    return float(np.mean(gammas)), float(np.mean(rhos))


# ---------------------------------------------------------------------------
# DFA training + checkpoint
# ---------------------------------------------------------------------------

def train_dfa_get_checkpoint(model, train_loader, test_loader, device,
                              total_epochs, t0, lr, wd):
    d = model.d_hidden; L = model.num_blocks
    Bs = [torch.randn(d, 10, device=device) / np.sqrt(10) for _ in range(L)]
    block_opts = [optim.AdamW(b.parameters(), lr=lr, weight_decay=wd)
                  for b in model.blocks]
    embed_opt = optim.AdamW(model.embed.parameters(), lr=lr, weight_decay=wd)
    head_opt = optim.AdamW(
        list(model.out_head.parameters()) + list(model.out_ln.parameters()),
        lr=lr, weight_decay=wd)
    scheds = ([optim.lr_scheduler.CosineAnnealingLR(o, T_max=total_epochs)
               for o in block_opts] +
              [optim.lr_scheduler.CosineAnnealingLR(embed_opt, T_max=total_epochs),
               optim.lr_scheduler.CosineAnnealingLR(head_opt, T_max=total_epochs)])
    ckpt = None
    for epoch in range(1, total_epochs + 1):
        model.train(); tl, c, t = 0, 0, 0
        for x, y in train_loader:
            x = x.view(x.size(0), -1).to(device); y = y.to(device); b = x.size(0)
            with torch.no_grad():
                lo, hi = model(x, return_hidden=True); lv = F.cross_entropy(lo, y)
                eT = lo.softmax(-1); eT[torch.arange(b), y] -= 1
            hL = hi[-1].detach()
            lo2 = F.cross_entropy(model.out_head(model.out_ln(hL)), y)
            head_opt.zero_grad(); lo2.backward(); head_opt.step()
            for l in range(L):
                a = (eT @ Bs[l].T).detach()
                rm = (a ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                f = model.blocks[l](hi[l].detach())
                ll = (f * (a / rm)).sum(-1).mean()
                block_opts[l].zero_grad(); ll.backward()
                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
                block_opts[l].step()
            a0 = (eT @ Bs[0].T).detach()
            r0 = (a0 ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
            el = (model.embed(x) * (a0 / r0)).sum(-1).mean()
            embed_opt.zero_grad(); el.backward(); embed_opt.step()
            tl += lv.item() * b; c += (lo.argmax(1) == y).sum().item(); t += b
        for s in scheds: s.step()
        if epoch == t0:
            acc = evaluate(model, test_loader, device)
            ckpt = {'model': copy.deepcopy(model.state_dict()),
                    'Bs': [B.clone() for B in Bs], 'acc': acc}
            print(f"  [DFA] Checkpoint at epoch {t0}: acc={acc:.4f}")
        if epoch % 10 == 0:
            print(f"  [DFA] Epoch {epoch}: acc={evaluate(model, test_loader, device):.4f}")
    return Bs, ckpt


# ---------------------------------------------------------------------------
# Norm matching: estimate per-block blend RMS from random_trainable after k epochs
# ---------------------------------------------------------------------------

def estimate_normmatched_gammas(model_init_state, Bs, train_loader, test_loader, device,
                                 t0, total_epochs, alpha, lr, lr_fb, wd, M,
                                 collect_epochs=10, input_dim=3072, d=256, L=4):
    """Run random_trainable for collect_epochs after handoff, collect per-block RMS
    of blended credits. Returns per-block gamma (scalar) for norm matching."""
    torch.manual_seed(42 + 9999)
    model_tmp = ResidualMLP(input_dim, d, 10, L).to(device)
    model_tmp.load_state_dict(model_init_state)

    torch.manual_seed(42 + 7777)
    vec_tmp = VectorCreditNet(d_hidden=d, s_dim=10).to(device)

    block_opts = [optim.AdamW(b.parameters(), lr=lr, weight_decay=wd)
                  for b in model_tmp.blocks]
    embed_opt = optim.AdamW(model_tmp.embed.parameters(), lr=lr, weight_decay=wd)
    head_opt = optim.AdamW(
        list(model_tmp.out_head.parameters()) + list(model_tmp.out_ln.parameters()),
        lr=lr, weight_decay=wd)
    vec_opt = optim.Adam(vec_tmp.parameters(), lr=lr_fb)
    scheds = ([optim.lr_scheduler.CosineAnnealingLR(o, T_max=total_epochs)
               for o in block_opts] +
              [optim.lr_scheduler.CosineAnnealingLR(embed_opt, T_max=total_epochs),
               optim.lr_scheduler.CosineAnnealingLR(head_opt, T_max=total_epochs)])
    for _ in range(t0):
        for s in scheds: s.step()

    eps_pert = 1e-3
    # Accumulate per-block blend RMS over all collect_epochs
    block_rms_accum = [[] for _ in range(L)]

    for epoch in range(t0 + 1, t0 + collect_epochs + 1):
        model_tmp.train(); vec_tmp.train()
        for x, y in train_loader:
            x = x.view(x.size(0), -1).to(device); y = y.to(device); batch = x.size(0)
            with torch.no_grad():
                lo, hi = model_tmp(x, return_hidden=True); lv = F.cross_entropy(lo, y)
                eT = lo.softmax(-1); eT[torch.arange(batch), y] -= 1; s = eT.detach()
            hL = hi[-1].detach()

            # Train Vec with standard perturbation loss
            t_L = torch.ones(batch, device=device)
            a_term = vec_tmp(hL, t_L, s)
            hL_req = hL.clone().requires_grad_(True)
            ce = F.cross_entropy(model_tmp.out_head(model_tmp.out_ln(hL_req)), y,
                                  reduction='sum')
            dL = torch.autograd.grad(ce, hL_req)[0].detach()
            loss_term = ((a_term - dL) ** 2).sum(-1).mean()
            lt = np.random.randint(0, L)
            h_l = hi[lt].detach(); t_l = torch.full((batch,), lt / L, device=device)
            a_l = vec_tmp(h_l, t_l, s)
            lp2 = torch.tensor(0.0, device=device)
            for _ in range(M):
                v = torch.randn_like(h_l); v = v / (v.norm(-1, keepdim=True) + 1e-8)
                with torch.no_grad():
                    lp = F.cross_entropy(
                        model_tmp.forward_from_layer(h_l + eps_pert * v, lt), y, reduction='none')
                    lm = F.cross_entropy(
                        model_tmp.forward_from_layer(h_l - eps_pert * v, lt), y, reduction='none')
                    gj = (lp - lm) / (2 * eps_pert)
                lp2 = lp2 + (((a_l * v).sum(-1) - gj.detach()) ** 2).mean()
            lp2 /= M
            vl = loss_term + lp2
            vec_opt.zero_grad(); vl.backward()
            torch.nn.utils.clip_grad_norm_(vec_tmp.parameters(), 1.0); vec_opt.step()

            # Compute and record blend RMS per block
            with torch.no_grad():
                for l in range(L):
                    a_dfa = (eT @ Bs[l].T).detach()
                    h_bl = hi[l].detach(); t_bl = torch.full((batch,), l / L, device=device)
                    a_vec = vec_tmp(h_bl, t_bl, s).detach()
                    rms_d = (a_dfa ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                    rms_v = (a_vec ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                    a_blend = alpha * a_vec / rms_v + (1 - alpha) * a_dfa / rms_d
                    block_rms_accum[l].append((a_blend ** 2).mean().sqrt().item())

            # Update head and blocks (needed to keep training realistic)
            lo2 = F.cross_entropy(model_tmp.out_head(model_tmp.out_ln(hL)), y)
            head_opt.zero_grad(); lo2.backward(); head_opt.step()
            for l in range(L):
                a_dfa = (eT @ Bs[l].T).detach()
                with torch.no_grad():
                    h_bl = hi[l].detach(); t_bl = torch.full((batch,), l / L, device=device)
                    a_vec = vec_tmp(h_bl, t_bl, s)
                    rms_d = (a_dfa ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                    rms_v = (a_vec ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                    a_blend = alpha * a_vec / rms_v + (1 - alpha) * a_dfa / rms_d
                rm = (a_blend ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                f = model_tmp.blocks[l](hi[l].detach())
                ll = (f * (a_blend / rm)).sum(-1).mean()
                block_opts[l].zero_grad(); ll.backward()
                torch.nn.utils.clip_grad_norm_(model_tmp.blocks[l].parameters(), 1.0)
                block_opts[l].step()
            a0_dfa = (eT @ Bs[0].T).detach()
            with torch.no_grad():
                h0 = hi[0].detach(); t0_t = torch.full((batch,), 0.0, device=device)
                av0 = vec_tmp(h0, t0_t, s)
                rd0 = (a0_dfa ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                rv0 = (av0 ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                a0b = alpha * av0 / rv0 + (1 - alpha) * a0_dfa / rd0
            r0 = (a0b ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
            el = (model_tmp.embed(x) * (a0b / r0)).sum(-1).mean()
            embed_opt.zero_grad(); el.backward(); embed_opt.step()
        for sch in scheds: sch.step()
        print(f"  [norm_match_calibration] Epoch {epoch}: "
              f"block_rms={[round(np.mean(block_rms_accum[l]), 5) for l in range(L)]}")

    # Per-block target RMS
    target_rms = [float(np.mean(block_rms_accum[l])) for l in range(L)]
    print(f"  [norm_match_calibration] Target RMS per block: {[round(r, 5) for r in target_rms]}")
    del model_tmp, vec_tmp
    return target_rms


# ---------------------------------------------------------------------------
# Branch runner
# ---------------------------------------------------------------------------

def run_branch(model, aux_net, Bs, train_loader, test_loader, device,
               t0, total_epochs, branch_type, alpha, lr, lr_fb, wd, M,
               branch_name='', freeze_after=None, normmatched_target_rms=None):
    """
    Run a training branch from a loaded checkpoint.

    branch_type options:
      'dfa'                       — pure DFA
      'blend_trainable'           — blend with Vec trained online (perturbation targets)
      'blend_zero_target'         — blend with Vec trained with ||a_aux||^2
      'blend_zero_normmatched'    — zero_target + blockwise norm matching
      'blend_perlayer_vector'     — no network, per-block nn.Parameter v_l (broadcast over batch)
      'blend_freeze_after_k'      — random Vec, train for freeze_after epochs then freeze

    freeze_after: int or None — for 'blend_freeze_after_k', number of epochs to train before freeze
    normmatched_target_rms: list of float, one per block — for 'blend_zero_normmatched'
    """
    d = model.d_hidden; L = model.num_blocks; eps_pert = 1e-3

    # Determine which aux nets need training
    trainable_types = {'blend_trainable', 'blend_zero_target', 'blend_zero_normmatched',
                       'blend_perlayer_vector', 'blend_freeze_after_k'}
    aux_trained = (branch_type in trainable_types) and (aux_net is not None)

    block_opts = [optim.AdamW(b.parameters(), lr=lr, weight_decay=wd)
                  for b in model.blocks]
    embed_opt = optim.AdamW(model.embed.parameters(), lr=lr, weight_decay=wd)
    head_opt = optim.AdamW(
        list(model.out_head.parameters()) + list(model.out_ln.parameters()),
        lr=lr, weight_decay=wd)
    # For perlayer_vector, only optimize the vectors (not a full network)
    if branch_type == 'blend_perlayer_vector' and isinstance(aux_net, PerLayerVector):
        aux_opt = optim.Adam(aux_net.parameters(), lr=lr_fb)
    elif aux_trained:
        aux_opt = optim.Adam(aux_net.parameters(), lr=lr_fb)
    else:
        aux_opt = None

    scheds = ([optim.lr_scheduler.CosineAnnealingLR(o, T_max=total_epochs)
               for o in block_opts] +
              [optim.lr_scheduler.CosineAnnealingLR(embed_opt, T_max=total_epochs),
               optim.lr_scheduler.CosineAnnealingLR(head_opt, T_max=total_epochs)])
    # Advance schedulers to match checkpoint epoch
    for _ in range(t0):
        for s in scheds: s.step()

    log = {'test_acc': [], 'train_loss': [], 'gamma': [], 'rho': [], 'alpha_eff': []}
    diag_epochs = set(
        list(range(t0 + 1, min(t0 + 6, total_epochs + 1))) +
        [t0 + 8, t0 + 10, t0 + 15, t0 + 20] +
        list(range(t0 + 10, total_epochs + 1, 10)) +
        [total_epochs])

    frozen = False  # Track whether Vec has been frozen (for freeze_after_k)

    for epoch in range(t0 + 1, total_epochs + 1):
        # Handle freeze_after_k: freeze after freeze_after epochs of training
        if branch_type == 'blend_freeze_after_k' and freeze_after is not None:
            epochs_trained = epoch - t0 - 1  # epochs completed since handoff
            if epochs_trained >= freeze_after and not frozen:
                if aux_net is not None:
                    aux_net.requires_grad_(False)
                    aux_net.eval()
                aux_opt = None
                frozen = True
                print(f"    [{branch_name}] Freezing Vec at epoch {epoch} "
                      f"(after {freeze_after} training epochs)")

        model.train()
        if aux_net is not None and (aux_opt is not None or not frozen):
            if aux_opt is not None:
                aux_net.train()
            else:
                aux_net.eval()
        elif aux_net is not None:
            aux_net.eval()

        tl, c, t = 0, 0, 0
        epoch_aux_norms, epoch_dfa_norms = [], []

        for x, y in train_loader:
            x = x.view(x.size(0), -1).to(device); y = y.to(device); batch = x.size(0)
            with torch.no_grad():
                lo, hi = model(x, return_hidden=True); lv = F.cross_entropy(lo, y)
                eT = lo.softmax(-1); eT[torch.arange(batch), y] -= 1; s = eT.detach()
            hL = hi[-1].detach()

            # ----------------------------------------------------------------
            # Train auxiliary network (if applicable)
            # ----------------------------------------------------------------
            if aux_opt is not None:
                if branch_type in ('blend_trainable', 'blend_freeze_after_k'):
                    # Standard perturbation targets
                    t_L = torch.ones(batch, device=device)
                    a_term = aux_net(hL, t_L, s)
                    hL_req = hL.clone().requires_grad_(True)
                    ce = F.cross_entropy(
                        model.out_head(model.out_ln(hL_req)), y, reduction='sum')
                    dL = torch.autograd.grad(ce, hL_req)[0].detach()
                    loss_term = ((a_term - dL) ** 2).sum(-1).mean()
                    lt = np.random.randint(0, L)
                    h_l = hi[lt].detach()
                    t_l = torch.full((batch,), lt / L, device=device)
                    a_l = aux_net(h_l, t_l, s)
                    lp2 = torch.tensor(0.0, device=device)
                    for _ in range(M):
                        v = torch.randn_like(h_l)
                        v = v / (v.norm(-1, keepdim=True) + 1e-8)
                        with torch.no_grad():
                            lp = F.cross_entropy(
                                model.forward_from_layer(h_l + eps_pert * v, lt),
                                y, reduction='none')
                            lm = F.cross_entropy(
                                model.forward_from_layer(h_l - eps_pert * v, lt),
                                y, reduction='none')
                            gj = (lp - lm) / (2 * eps_pert)
                        lp2 = lp2 + (((a_l * v).sum(-1) - gj.detach()) ** 2).mean()
                    lp2 /= M
                    vl = loss_term + lp2

                elif branch_type in ('blend_zero_target', 'blend_zero_normmatched'):
                    # Minimize ||a_aux||^2 — teaches the network to output zero
                    lt = np.random.randint(0, L)
                    h_l = hi[lt].detach()
                    t_l = torch.full((batch,), lt / L, device=device)
                    a_l = aux_net(h_l, t_l, s)
                    vl = (a_l ** 2).sum(-1).mean()

                elif branch_type == 'blend_perlayer_vector':
                    # Per-layer vector: train v_l with perturbation loss
                    # v_l is used as gradient surrogate for a random layer
                    lt = np.random.randint(0, L)
                    h_l = hi[lt].detach()
                    t_l = torch.full((batch,), lt / L, device=device)
                    aux_net.set_block(lt)
                    # v_l is the per-layer vector (same for all samples in batch)
                    a_l = aux_net(h_l, t_l, s)  # (batch, d) — same v_lt broadcast
                    lp2 = torch.tensor(0.0, device=device)
                    for _ in range(M):
                        v = torch.randn_like(h_l)
                        v = v / (v.norm(-1, keepdim=True) + 1e-8)
                        with torch.no_grad():
                            lp = F.cross_entropy(
                                model.forward_from_layer(h_l + eps_pert * v, lt),
                                y, reduction='none')
                            lm = F.cross_entropy(
                                model.forward_from_layer(h_l - eps_pert * v, lt),
                                y, reduction='none')
                            gj = (lp - lm) / (2 * eps_pert)
                        # <v_l, v_dir> — v_l is shared across batch, v is per-sample
                        # (a_l * v).sum(-1) computes dot product per sample
                        lp2 = lp2 + (((a_l * v).sum(-1) - gj.detach()) ** 2).mean()
                    lp2 /= M
                    vl = lp2

                else:
                    vl = None

                if vl is not None:
                    aux_opt.zero_grad(); vl.backward()
                    torch.nn.utils.clip_grad_norm_(aux_net.parameters(), 1.0)
                    aux_opt.step()

            # ----------------------------------------------------------------
            # Compute credits for each block
            # ----------------------------------------------------------------
            dfa_credits = [(eT @ Bs[l].T).detach() for l in range(L)]
            credits = []
            for l in range(L):
                a_dfa = dfa_credits[l]
                rms_d = (a_dfa ** 2).mean(-1, keepdim=True).sqrt() + 1e-6

                if branch_type == 'dfa':
                    credits.append(a_dfa / rms_d)
                else:
                    # All blend branches
                    h_l = hi[l].detach()
                    t_l = torch.full((batch,), l / L, device=device)
                    with torch.no_grad():
                        if isinstance(aux_net, PerLayerVector):
                            aux_net.set_block(l)
                        a_aux = aux_net(h_l, t_l, s).detach()
                    rms_v = (a_aux ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                    a_blend = alpha * a_aux / rms_v + (1 - alpha) * a_dfa / rms_d

                    # For norm-matched zero_target: scale blended credit by per-block gamma
                    if branch_type == 'blend_zero_normmatched' and normmatched_target_rms is not None:
                        current_rms = (a_blend ** 2).mean().sqrt().item()
                        gamma_nm = normmatched_target_rms[l] / (current_rms + 1e-8)
                        a_blend = a_blend * gamma_nm

                    credits.append(a_blend)

                # Track norms for alpha_eff
                a_c = credits[-1]
                if branch_type == 'dfa':
                    epoch_aux_norms.append(0.0)
                    epoch_dfa_norms.append(a_c.norm().item())
                else:
                    a_dfa_n = a_dfa / rms_d
                    rms_v2 = (a_aux ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                    epoch_aux_norms.append((alpha * a_aux / rms_v2).norm().item())
                    epoch_dfa_norms.append(((1 - alpha) * a_dfa_n).norm().item())

            # ----------------------------------------------------------------
            # Update output head (local exact gradient — allowed)
            # ----------------------------------------------------------------
            lo2 = F.cross_entropy(model.out_head(model.out_ln(hL)), y)
            head_opt.zero_grad(); lo2.backward(); head_opt.step()

            # ----------------------------------------------------------------
            # Update blocks with local surrogate
            # ----------------------------------------------------------------
            for l in range(L):
                a = credits[l]
                rm = (a ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
                f = model.blocks[l](hi[l].detach())
                ll = (f * (a / rm)).sum(-1).mean()
                block_opts[l].zero_grad(); ll.backward()
                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
                block_opts[l].step()

            # Update embedding with block-0 credit
            a0 = credits[0]
            r0 = (a0 ** 2).mean(-1, keepdim=True).sqrt() + 1e-6
            el = (model.embed(x) * (a0 / r0)).sum(-1).mean()
            embed_opt.zero_grad(); el.backward(); embed_opt.step()

            tl += lv.item() * batch; c += (lo.argmax(1) == y).sum().item(); t += batch

        for sch in scheds: sch.step()
        ta = evaluate(model, test_loader, device)
        log['test_acc'].append(ta); log['train_loss'].append(tl / t)

        mean_aux = np.mean(epoch_aux_norms) if epoch_aux_norms else 0.0
        mean_dfa = np.mean(epoch_dfa_norms) if epoch_dfa_norms else 1.0
        aeff = mean_aux / (mean_aux + mean_dfa + 1e-12)
        log['alpha_eff'].append((epoch, aeff))

        if epoch in diag_epochs:
            cm = 'blend' if branch_type != 'dfa' else 'dfa'
            gamma, rho = compute_diagnostics(
                model, aux_net if branch_type != 'dfa' else None,
                Bs, test_loader, device, cm, alpha)
            log['gamma'].append((epoch, gamma)); log['rho'].append((epoch, rho))
            if epoch <= t0 + 15 or epoch % 20 == 0 or epoch == total_epochs:
                frozen_str = ' [FROZEN]' if frozen else ''
                print(f"    [{branch_name}]{frozen_str} Ep {epoch}: acc={ta:.4f}, "
                      f"G={gamma:.4f}, r={rho:.4f}, aeff={aeff:.3f}")
        elif epoch % 10 == 0 or epoch == total_epochs:
            frozen_str = ' [FROZEN]' if frozen else ''
            print(f"    [{branch_name}]{frozen_str} Ep {epoch}: acc={ta:.4f}")

    return log


# ---------------------------------------------------------------------------
# Main experiment
# ---------------------------------------------------------------------------

def run_experiment(args):
    device = torch.device(f'cuda:{args.gpu}' if torch.cuda.is_available() else 'cpu')
    print(f"Using device: {device}")
    os.makedirs(args.output_dir, exist_ok=True)
    torch.manual_seed(args.seed); np.random.seed(args.seed)
    torch.cuda.manual_seed_all(args.seed)

    train_loader, test_loader = get_cifar10(args.batch_size)
    input_dim = 32 * 32 * 3; L = args.num_blocks; d = args.d_hidden

    # ----------------------------------------------------------------
    # Step 1: Train DFA and capture checkpoint at t0
    # ----------------------------------------------------------------
    print(f"\n{'='*60}\nTraining DFA baseline (checkpoint at t0={args.t0})\n{'='*60}")
    model_dfa = ResidualMLP(input_dim, d, 10, L).to(device)
    Bs, ckpt = train_dfa_get_checkpoint(
        model_dfa, train_loader, test_loader, device,
        args.epochs, args.t0, args.lr, args.wd)
    print(f"  Checkpoint acc at t0={args.t0}: {ckpt['acc']:.4f}")

    # ----------------------------------------------------------------
    # Step 2: Estimate per-block target RMS for norm-matched zero_target
    #         Run random_trainable for 10 epochs from checkpoint, collect stats.
    # ----------------------------------------------------------------
    print(f"\n{'='*60}\nCalibrating norm-matching (random_trainable 10 epochs)\n{'='*60}")
    normmatched_target_rms = estimate_normmatched_gammas(
        model_init_state=ckpt['model'],
        Bs=ckpt['Bs'],
        train_loader=train_loader,
        test_loader=test_loader,
        device=device,
        t0=args.t0,
        total_epochs=args.epochs,
        alpha=args.alpha,
        lr=args.lr,
        lr_fb=args.lr_fb,
        wd=args.wd,
        M=args.M,
        collect_epochs=10,
        input_dim=input_dim,
        d=d,
        L=L)
    print(f"  Per-block target RMS: {[round(r, 5) for r in normmatched_target_rms]}")

    # ----------------------------------------------------------------
    # Step 3: Define and run all 8 branches
    # ----------------------------------------------------------------
    VEC_SEED = args.seed + 7777

    def make_vec():
        torch.manual_seed(VEC_SEED)
        return VectorCreditNet(d_hidden=d, s_dim=10).to(device)

    def make_perlayer():
        torch.manual_seed(VEC_SEED)
        return PerLayerVector(d_hidden=d, num_blocks=L).to(device)

    # (name, branch_type, aux_factory, freeze_after, use_normmatched)
    branches = [
        ('continue_DFA',                  'dfa',                  lambda: None,       None, False),
        ('blend_random_trainable',         'blend_trainable',      make_vec,            None, False),
        ('blend_zero_target_trainable',    'blend_zero_target',    make_vec,            None, False),
        ('blend_zero_target_normmatched',  'blend_zero_normmatched', make_vec,          None, True),
        ('blend_perlayer_vector',          'blend_perlayer_vector', make_perlayer,       None, False),
        ('blend_random_freeze_after_1',    'blend_freeze_after_k', make_vec,            1,    False),
        ('blend_random_freeze_after_5',    'blend_freeze_after_k', make_vec,            5,    False),
        ('blend_random_freeze_after_10',   'blend_freeze_after_k', make_vec,            10,   False),
    ]

    all_results = {}
    for bname, btype, aux_factory, freeze_after, use_normmatched in branches:
        print(f"\n{'='*60}\n{bname}\n{'='*60}")
        model_b = ResidualMLP(input_dim, d, 10, L).to(device)
        model_b.load_state_dict(ckpt['model'])
        aux_net_b = aux_factory()
        nm_rms = normmatched_target_rms if use_normmatched else None

        log = run_branch(
            model_b, aux_net_b, ckpt['Bs'],
            train_loader, test_loader, device,
            args.t0, args.epochs, btype,
            args.alpha, args.lr, args.lr_fb, args.wd, args.M,
            branch_name=bname,
            freeze_after=freeze_after,
            normmatched_target_rms=nm_rms)
        all_results[bname] = log
        print(f"  {bname} final acc: {log['test_acc'][-1]:.4f}")

    # ----------------------------------------------------------------
    # Step 4: Summary table
    # ----------------------------------------------------------------
    dfa_final = all_results['continue_DFA']['test_acc'][-1]

    print(f"\n{'='*95}")
    print("SUMMARY — Phase 10A.7: Minimal Auxiliary Compression")
    print(f"{'='*95}")
    print(f"{'Branch':<36} {'@20':>6} {'final':>7} {'diff':>7} "
          f"{'mG_5:15':>9} {'mr_5:15':>9} {'aeff':>7}")
    print("-" * 83)

    for bname, log in all_results.items():
        accs = log['test_acc']
        idx20 = max(0, 20 - args.t0 - 1)
        acc20 = accs[idx20] if len(accs) > idx20 else accs[-1]
        final = accs[-1]
        diff = final - dfa_final
        gammas_e = [g for e, g in log['gamma'] if args.t0 < e <= args.t0 + 15]
        rhos_e = [r for e, r in log['rho'] if args.t0 < e <= args.t0 + 15]
        aeffs_e = [a for e, a in log['alpha_eff'] if args.t0 < e <= args.t0 + 15]
        mg = float(np.mean(gammas_e)) if gammas_e else float('nan')
        mr = float(np.mean(rhos_e)) if rhos_e else float('nan')
        mae = float(np.mean(aeffs_e)) if aeffs_e else float('nan')
        print(f"{bname:<36} {acc20:>6.4f} {final:>7.4f} {diff:>+7.4f} "
              f"{mg:>9.4f} {mr:>9.4f} {mae:>7.3f}")

    # ----------------------------------------------------------------
    # Step 5: Save results
    # ----------------------------------------------------------------
    save_data = {
        'args': vars(args),
        'dfa_ckpt_acc': float(ckpt['acc']),
        'normmatched_target_rms': normmatched_target_rms,
    }
    for bname, log in all_results.items():
        save_data[bname] = {
            'test_acc': log['test_acc'],
            'train_loss': log['train_loss'],
            'gamma': log['gamma'],
            'rho': log['rho'],
            'alpha_eff': log['alpha_eff'],
        }
    out_path = os.path.join(args.output_dir,
                            f'minimal_aux_compression_t{args.t0}_s{args.seed}.json')
    with open(out_path, 'w') as f:
        json.dump(save_data, f, indent=2, default=float)
    print(f"\nSaved to {out_path}")

    # ----------------------------------------------------------------
    # Step 6: Judgment
    # ----------------------------------------------------------------
    print(f"\n{'='*60}\nJUDGMENT\n{'='*60}")
    r = {bname: log['test_acc'][-1] for bname, log in all_results.items()}
    dfa = r['continue_DFA']
    rt  = r.get('blend_random_trainable', float('nan'))
    zt  = r.get('blend_zero_target_trainable', float('nan'))
    znm = r.get('blend_zero_target_normmatched', float('nan'))
    plv = r.get('blend_perlayer_vector', float('nan'))
    f1  = r.get('blend_random_freeze_after_1', float('nan'))
    f5  = r.get('blend_random_freeze_after_5', float('nan'))
    f10 = r.get('blend_random_freeze_after_10', float('nan'))

    print(f"  DFA={dfa:.4f}  rt={rt:.4f}  zt={zt:.4f}  znm={znm:.4f}  "
          f"plv={plv:.4f}  f1={f1:.4f}  f5={f5:.4f}  f10={f10:.4f}")

    thr = 0.003

    # Norm matching diagnosis
    if abs(znm - rt) < thr:
        print("  -> zero_normmatched ≈ random_trainable: "
              "gain is primarily norm/step-size effect, not directional signal")
    elif znm > zt + thr and znm < rt - thr:
        print("  -> zero_normmatched is between zt and rt: "
              "norm helps but signal direction still adds value")
    elif znm > rt - thr:
        print("  -> zero_normmatched ≈ random_trainable: "
              "norm matching fully recovers the blend gain (direction irrelevant)")

    # Per-layer vector vs full network
    if abs(plv - rt) < thr:
        print("  -> perlayer_vector ≈ random_trainable: "
              "per-block scalar direction sufficient; no input conditioning needed")
    elif plv > dfa + thr:
        print(f"  -> perlayer_vector improves over DFA (+{plv-dfa:.4f}): "
              "even input-agnostic per-block direction helps")
    else:
        print(f"  -> perlayer_vector does NOT improve over DFA: "
              "per-block vectors alone insufficient without input conditioning")

    # Freeze timing analysis
    print(f"\n  Freeze timing: f1={f1:.4f}  f5={f5:.4f}  f10={f10:.4f}  rt={rt:.4f}")
    if f10 > f1 + thr:
        print("  -> More training before freeze is better: "
              "Vec needs time to learn useful direction")
    if abs(f10 - rt) < thr:
        print("  -> freeze_after_10 ≈ random_trainable: "
              "most learning happens in first 10 epochs; subsequent training marginal")
    elif rt > f10 + thr:
        print("  -> continuous training > freeze_after_10: "
              "ongoing adaptation to evolving forward network matters")


def main():
    parser = argparse.ArgumentParser(
        description='Phase 10A.7: Minimal Auxiliary Compression')
    parser.add_argument('--num_blocks', type=int, default=4)
    parser.add_argument('--d_hidden', type=int, default=256)
    parser.add_argument('--batch_size', type=int, default=128)
    parser.add_argument('--epochs', type=int, default=100)
    parser.add_argument('--t0', type=int, default=5)
    parser.add_argument('--alpha', type=float, default=0.75)
    parser.add_argument('--lr', type=float, default=1e-3)
    parser.add_argument('--lr_fb', type=float, default=1e-3)
    parser.add_argument('--wd', type=float, default=0.01)
    parser.add_argument('--M', type=int, default=4)
    parser.add_argument('--seed', type=int, default=42)
    parser.add_argument('--gpu', type=int, default=0)
    parser.add_argument('--output_dir', type=str, default='results/minimal_aux_compression')
    args = parser.parse_args()
    run_experiment(args)


if __name__ == '__main__':
    main()