Add Phase 4 diagnostic dissection: frozen credit recovery, online shallow scan, vector field pilot

Key findings:
- Frozen CIFAR: estimators CAN recover credit (SB best, CB 20x > DFA)
- Online shallow: cb_eT wr=0.2 tgw=1.0 achieves S1>0, S2 marginal
- Vector credit field: 0.91-0.96 Gamma/rho on synthetic (vs 0.34 scalar CB)
- Direct vector field avoids scalar V curvature problem entirely

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

3 files changed, 2118 insertions, 0 deletions

diff --git a/experiments/cifar_frozen_credit_recovery.py b/experiments/cifar_frozen_credit_recovery.py
new file mode 100644
index 0000000..5d39308
--- /dev/null
+++ b/experiments/cifar_frozen_credit_recovery.py
@@ -0,0 +1,693 @@
+"""
+Phase A: Frozen CIFAR Credit Recovery.
+
+Goal: Separate "estimator problem" from "forward exploitability problem".
+1. Train a BP reference network to convergence, freeze it.
+2. On frozen features, train credit estimators (state bridge, scalar CB with eT/deltaL).
+3. Evaluate Gamma, rho, nudging per layer.
+
+This answers: can the credit estimator recover useful local credit from fixed representations?
+"""
+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
+import time
+
+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 ValueNet, SinusoidalTimeEmbed, create_ema_model, update_ema
+from models.state_bridge import StateBridgeNet
+from metrics.credit_metrics import (
+    cosine_similarity_batch, perturbation_correlation, nudging_test,
+    offline_bp_cosine
+)
+
+
+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)
+    train_loader = DataLoader(trainset, batch_size=batch_size, shuffle=True, num_workers=4, pin_memory=True)
+    test_loader = DataLoader(testset, batch_size=batch_size, shuffle=False, num_workers=4, pin_memory=True)
+    return train_loader, test_loader
+
+
+def evaluate(model, test_loader, device):
+    model.eval()
+    correct, total = 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)
+            logits = model(x)
+            correct += (logits.argmax(1) == y).sum().item()
+            total += x.size(0)
+    return correct / total
+
+
+# =============================================================================
+# Step 1: Train BP reference network
+# =============================================================================
+def train_bp_reference(model, train_loader, test_loader, device, epochs=100, lr=1e-3, wd=0.01):
+    """Train BP reference to convergence."""
+    optimizer = optim.AdamW(model.parameters(), lr=lr, weight_decay=wd)
+    scheduler = optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=epochs)
+
+    for epoch in range(1, epochs + 1):
+        model.train()
+        total_loss, correct, total = 0, 0, 0
+        for x, y in train_loader:
+            x = x.view(x.size(0), -1).to(device)
+            y = y.to(device)
+            logits = model(x)
+            loss = F.cross_entropy(logits, y)
+            optimizer.zero_grad()
+            loss.backward()
+            optimizer.step()
+            total_loss += loss.item() * x.size(0)
+            correct += (logits.argmax(1) == y).sum().item()
+            total += x.size(0)
+        scheduler.step()
+        if epoch % 10 == 0 or epoch == 1:
+            test_acc = evaluate(model, test_loader, device)
+            print(f"  [BP ref] Epoch {epoch}: loss={total_loss/total:.4f}, "
+                  f"train_acc={correct/total:.4f}, test_acc={test_acc:.4f}")
+
+    test_acc = evaluate(model, test_loader, device)
+    print(f"  [BP ref] Final test accuracy: {test_acc:.4f}")
+    return test_acc
+
+
+# =============================================================================
+# Step 2: Train estimators on frozen features
+# =============================================================================
+
+def train_state_bridge_frozen(model, train_loader, device, args):
+    """Train state bridge on frozen BP features."""
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = 10
+
+    state_pred = StateBridgeNet(
+        d_hidden=d, s_dim=num_classes, time_embed_dim=32,
+        hidden_dim=256, num_layers=3
+    ).to(device)
+    state_opt = optim.Adam(state_pred.parameters(), lr=args.lr_fb)
+
+    model.eval()
+    for epoch in range(1, args.estimator_epochs + 1):
+        state_pred.train()
+        total_loss = 0
+        n = 0
+        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():
+                logits, hiddens = model(x, return_hidden=True)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+                s = e_T.detach()
+                hL_det = hiddens[-1].detach()
+
+            # Train state predictor
+            state_loss = 0.0
+            for l in range(L):
+                h_l_det = hiddens[l].detach()
+                t_l = torch.full((batch,), l / L, device=device)
+                pred_hL = state_pred(h_l_det, t_l, s)
+                target_norm = hL_det.norm(dim=-1, keepdim=True).clamp(min=1.0)
+                state_loss = state_loss + (((pred_hL - hL_det) / target_norm) ** 2).sum(dim=-1).mean()
+            state_loss = state_loss / L
+
+            state_opt.zero_grad()
+            state_loss.backward()
+            state_opt.step()
+            total_loss += state_loss.item() * batch
+            n += batch
+
+        if epoch % 20 == 0 or epoch == 1:
+            print(f"    [SB] Epoch {epoch}: state_loss={total_loss/n:.6f}")
+
+    return state_pred
+
+
+def train_scalar_cb_frozen(model, train_loader, device, args, s_type='eT'):
+    """
+    Train scalar credit bridge on frozen BP features.
+    s_type: 'eT' (softmax error, dim=10) or 'deltaL' (grad_{h_L} CE, dim=d_hidden)
+    """
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = 10
+
+    if s_type == 'eT':
+        s_dim = num_classes
+    elif s_type == 'deltaL':
+        s_dim = d
+    else:
+        raise ValueError(f"Unknown s_type: {s_type}")
+
+    value_net = ValueNet(
+        d_hidden=d, s_dim=s_dim, time_embed_dim=32,
+        hidden_dim=256, num_layers=3
+    ).to(device)
+    value_net_ema = create_ema_model(value_net)
+    value_opt = optim.Adam(value_net.parameters(), lr=args.lr_fb)
+
+    lam = args.lam
+    K_samples = args.K
+    sigma_bridge = args.sigma_bridge
+    ema_momentum = args.ema_momentum
+    term_grad_weight = args.term_grad_weight
+
+    model.eval()
+    for epoch in range(1, args.estimator_epochs + 1):
+        value_net.train()
+        total_vloss = 0
+        total_term = 0
+        total_tgrad = 0
+        total_bridge = 0
+        n = 0
+
+        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():
+                logits, hiddens = model(x, return_hidden=True)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+                true_loss = F.cross_entropy(logits, y, reduction='none').detach()
+
+            hL_det = hiddens[-1].detach()
+
+            # Compute s (conditioning code)
+            if s_type == 'eT':
+                s = e_T.detach()
+            elif s_type == 'deltaL':
+                # delta_L = grad_{h_L} CE (output-layer-local, allowed)
+                hL_req = hL_det.clone().requires_grad_(True)
+                logits_for_s = model.out_head(model.out_ln(hL_req))
+                ce_for_s = F.cross_entropy(logits_for_s, y, reduction='sum')
+                delta_L = torch.autograd.grad(ce_for_s, hL_req, create_graph=False)[0].detach()
+                s = delta_L
+
+            # Terminal boundary
+            t_L = torch.ones(batch, device=device)
+            V_terminal = value_net(hL_det, t_L, s)
+            loss_term = ((V_terminal - true_loss) ** 2).mean()
+
+            # Terminal gradient matching
+            loss_tgrad = torch.tensor(0.0, device=device)
+            if term_grad_weight > 0:
+                hL_req = hL_det.clone().requires_grad_(True)
+                V_at_L = value_net(hL_req, t_L, s)
+                grad_V_L = torch.autograd.grad(V_at_L.sum(), hL_req, create_graph=True)[0]
+                # Exact terminal gradient (output-layer-local)
+                hL_req2 = hL_det.clone().requires_grad_(True)
+                logits_tgt = model.out_head(model.out_ln(hL_req2))
+                ce_loss = F.cross_entropy(logits_tgt, y, reduction='sum')
+                a_L_exact = torch.autograd.grad(ce_loss, hL_req2, create_graph=False)[0].detach()
+                loss_tgrad = ((grad_V_L - a_L_exact) ** 2).sum(dim=-1).mean()
+
+            # Bridge consistency
+            loss_bridge = 0.0
+            for l in range(L):
+                h_l_det = hiddens[l].detach()
+                t_l = torch.full((batch,), l / L, device=device)
+                t_l_next = torch.full((batch,), (l + 1) / L, device=device)
+                V_l = value_net(h_l_det, t_l, s)
+
+                with torch.no_grad():
+                    h_next_det = hiddens[l + 1].detach()
+                    log_terms = []
+                    for k in range(K_samples):
+                        noise = sigma_bridge * torch.randn_like(h_next_det)
+                        V_next = value_net_ema(h_next_det + noise, t_l_next, s)
+                        log_terms.append(-V_next / lam)
+                    log_stack = torch.stack(log_terms, dim=-1)
+                    V_target = -lam * (torch.logsumexp(log_stack, dim=-1) - np.log(K_samples))
+
+                loss_bridge = loss_bridge + ((V_l - V_target.detach()) ** 2).mean()
+            loss_bridge = loss_bridge / L
+
+            value_loss = loss_term + loss_bridge + term_grad_weight * loss_tgrad
+
+            value_opt.zero_grad()
+            value_loss.backward()
+            torch.nn.utils.clip_grad_norm_(value_net.parameters(), 1.0)
+            value_opt.step()
+            update_ema(value_net, value_net_ema, ema_momentum)
+
+            total_vloss += value_loss.item() * batch
+            total_term += loss_term.item() * batch
+            total_tgrad += loss_tgrad.item() * batch
+            total_bridge += (loss_bridge.item() if isinstance(loss_bridge, torch.Tensor) else loss_bridge) * batch
+            n += batch
+
+        if epoch % 20 == 0 or epoch == 1:
+            print(f"    [CB_{s_type}] Epoch {epoch}: vloss={total_vloss/n:.6f}, "
+                  f"term={total_term/n:.6f}, tgrad={total_tgrad/n:.6f}, bridge={total_bridge/n:.6f}")
+
+    return value_net, value_net_ema
+
+
+# =============================================================================
+# Step 3: Evaluate credit quality on frozen features
+# =============================================================================
+
+def evaluate_credits(model, test_loader, device, estimators, args):
+    """
+    Evaluate credit quality for all estimators on frozen BP features.
+
+    Args:
+        estimators: dict of {name: {'type': 'sb'/'cb', 'net': ..., 's_type': ...}}
+    Returns:
+        dict of {name: {per-layer metrics}}
+    """
+    model.eval()
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = 10
+
+    # Accumulate over multiple test batches for robust statistics
+    all_results = {}
+    for name in estimators:
+        all_results[name] = {
+            'bp_cosine': [[] for _ in range(L)],
+            'perturbation_rho': [0.0] * L,
+            'nudging_0.001': [0.0] * L,
+            'nudging_0.003': [0.0] * L,
+            'nudging_0.01': [0.0] * L,
+        }
+
+    # Also add DFA baseline
+    dfa_Bs = [torch.randn(d, num_classes, device=device) / np.sqrt(num_classes) for _ in range(L)]
+    all_results['dfa'] = {
+        'bp_cosine': [[] for _ in range(L)],
+        'perturbation_rho': [0.0] * L,
+        'nudging_0.001': [0.0] * L,
+        'nudging_0.003': [0.0] * L,
+        'nudging_0.01': [0.0] * L,
+    }
+
+    n_batches_diag = min(10, len(test_loader))  # Use multiple batches
+    batch_idx = 0
+
+    for x, y in test_loader:
+        if batch_idx >= n_batches_diag:
+            break
+        batch_idx += 1
+
+        x = x.view(x.size(0), -1).to(device)
+        y = y.to(device)
+        batch = x.size(0)
+
+        # Get BP gradients (ground truth for Gamma)
+        # Temporarily enable grad on model params for BP gradient computation
+        for p in model.parameters():
+            p.requires_grad_(True)
+        model.zero_grad()
+        logits_bp, hiddens_bp = model(x, return_hidden=True)
+        for l in range(L + 1):
+            hiddens_bp[l].retain_grad()
+        loss_bp = F.cross_entropy(logits_bp, y)
+        loss_bp.backward()
+        bp_grads = {l: hiddens_bp[l].grad.detach().clone() for l in range(L + 1)}
+        # Re-freeze model
+        for p in model.parameters():
+            p.requires_grad_(False)
+
+        # Clean forward
+        with torch.no_grad():
+            logits, hiddens = model(x, return_hidden=True)
+            e_T = logits.softmax(dim=-1)
+            e_T[torch.arange(batch), y] -= 1
+            s_eT = e_T.detach()
+
+        hL_det = hiddens[-1].detach()
+
+        # Compute delta_L for deltaL conditioning
+        hL_req = hL_det.clone().requires_grad_(True)
+        logits_for_delta = model.out_head(model.out_ln(hL_req))
+        ce_for_delta = F.cross_entropy(logits_for_delta, y, reduction='sum')
+        delta_L = torch.autograd.grad(ce_for_delta, hL_req, create_graph=False)[0].detach()
+
+        for l in range(L):
+            h_l = hiddens[l].detach()
+            t_l = torch.full((batch,), l / L, device=device)
+
+            # Forward function for perturbation and nudging
+            def make_fwd_fn(start_l):
+                def fwd_fn(h):
+                    with torch.no_grad():
+                        curr = h
+                        for i in range(start_l, L):
+                            curr = curr + model.blocks[i](curr)
+                        out = model.out_head(model.out_ln(curr))
+                        return F.cross_entropy(out, y, reduction='none')
+                return fwd_fn
+
+            fwd_fn = make_fwd_fn(l)
+
+            # --- DFA credit ---
+            a_dfa = (s_eT @ dfa_Bs[l].T).detach()
+            bp_cos_dfa = cosine_similarity_batch(a_dfa, bp_grads[l])
+            all_results['dfa']['bp_cosine'][l].append(bp_cos_dfa)
+
+            if batch_idx == 1:  # Only compute rho/nudging on first batch (expensive)
+                rho_dfa = perturbation_correlation(h_l, a_dfa, fwd_fn, epsilon=1e-3, M=32)
+                all_results['dfa']['perturbation_rho'][l] = rho_dfa
+                for eta in [0.001, 0.003, 0.01]:
+                    nud = nudging_test(h_l, a_dfa, fwd_fn, eta=eta)
+                    all_results['dfa'][f'nudging_{eta}'][l] = nud
+
+            # --- Estimator credits ---
+            for name, est in estimators.items():
+                if est['type'] == 'sb':
+                    net = est['net']
+                    net.eval()
+                    h_l_req = h_l.clone().requires_grad_(True)
+                    pred_hL = net(h_l_req, t_l, s_eT)
+                    pred_logits = model.out_head(model.out_ln(pred_hL))
+                    pred_loss = F.cross_entropy(pred_logits, y, reduction='sum')
+                    a_l = torch.autograd.grad(pred_loss, h_l_req, create_graph=False)[0].detach()
+
+                elif est['type'] == 'cb':
+                    net = est['net']
+                    net.eval()
+                    s_type = est['s_type']
+                    if s_type == 'eT':
+                        s = s_eT
+                    elif s_type == 'deltaL':
+                        s = delta_L
+                    else:
+                        raise ValueError(f"Unknown s_type: {s_type}")
+
+                    h_l_req = h_l.clone().requires_grad_(True)
+                    V_l = net(h_l_req, t_l, s)
+                    a_l = torch.autograd.grad(V_l.sum(), h_l_req, create_graph=False)[0].detach()
+                else:
+                    raise ValueError(f"Unknown estimator type: {est['type']}")
+
+                bp_cos = cosine_similarity_batch(a_l, bp_grads[l])
+                all_results[name]['bp_cosine'][l].append(bp_cos)
+
+                if batch_idx == 1:
+                    rho = perturbation_correlation(h_l, a_l, fwd_fn, epsilon=1e-3, M=32)
+                    all_results[name]['perturbation_rho'][l] = rho
+                    for eta in [0.001, 0.003, 0.01]:
+                        nud = nudging_test(h_l, a_l, fwd_fn, eta=eta)
+                        all_results[name][f'nudging_{eta}'][l] = nud
+
+    # Average bp_cosine over batches
+    for name in all_results:
+        for l in range(L):
+            vals = all_results[name]['bp_cosine'][l]
+            all_results[name]['bp_cosine'][l] = float(np.mean(vals)) if vals else 0.0
+
+    return all_results
+
+
+def evaluate_state_bridge_pred_error(model, state_pred, test_loader, device):
+    """Evaluate state bridge's terminal state prediction error."""
+    model.eval()
+    state_pred.eval()
+    L = model.num_blocks
+
+    total_error = [0.0] * L
+    n = 0
+    for x, y in test_loader:
+        x = x.view(x.size(0), -1).to(device)
+        y = y.to(device)
+        batch = x.size(0)
+
+        with torch.no_grad():
+            logits, hiddens = model(x, return_hidden=True)
+            e_T = logits.softmax(dim=-1)
+            e_T[torch.arange(batch), y] -= 1
+            s = e_T.detach()
+            hL = hiddens[-1]
+
+            for l in range(L):
+                h_l = hiddens[l]
+                t_l = torch.full((batch,), l / L, device=x.device)
+                pred_hL = state_pred(h_l, t_l, s)
+                error = ((pred_hL - hL) ** 2).sum(dim=-1).mean().item()
+                total_error[l] += error * batch
+        n += batch
+
+    return [e / n for e in total_error]
+
+
+# =============================================================================
+# 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(batch_size=args.batch_size)
+    input_dim = 32 * 32 * 3
+    num_classes = 10
+
+    # ----- Step 1: Train BP reference -----
+    print(f"\n{'='*60}")
+    print(f"Step 1: Train BP reference (L={args.num_blocks}, d={args.d_hidden})")
+    print(f"{'='*60}")
+
+    bp_ckpt_path = os.path.join(args.output_dir, f'bp_ref_L{args.num_blocks}_d{args.d_hidden}_s{args.seed}.pt')
+
+    model = ResidualMLP(input_dim, args.d_hidden, num_classes, args.num_blocks).to(device)
+
+    if os.path.exists(bp_ckpt_path) and not args.retrain_bp:
+        print(f"  Loading BP reference from {bp_ckpt_path}")
+        model.load_state_dict(torch.load(bp_ckpt_path, map_location=device))
+        bp_acc = evaluate(model, test_loader, device)
+        print(f"  BP reference test accuracy: {bp_acc:.4f}")
+    else:
+        bp_acc = train_bp_reference(model, train_loader, test_loader, device,
+                                     epochs=args.bp_epochs, lr=args.lr, wd=args.wd)
+        torch.save(model.state_dict(), bp_ckpt_path)
+        print(f"  Saved BP reference to {bp_ckpt_path}")
+
+    # Freeze the model completely
+    model.eval()
+    for p in model.parameters():
+        p.requires_grad_(False)
+
+    # ----- Step 2: Train estimators -----
+    print(f"\n{'='*60}")
+    print(f"Step 2: Train estimators ({args.estimator_epochs} epochs each)")
+    print(f"{'='*60}")
+
+    estimators = {}
+
+    # 2a. State Bridge with s=eT
+    print("\n--- State Bridge (s=eT) ---")
+    torch.manual_seed(args.seed + 1000)
+    sb = train_state_bridge_frozen(model, train_loader, device, args)
+    estimators['sb_eT'] = {'type': 'sb', 'net': sb, 's_type': 'eT'}
+
+    # 2b. Scalar CB with s=eT
+    print("\n--- Scalar CB (s=eT) ---")
+    torch.manual_seed(args.seed + 2000)
+    cb_eT, cb_eT_ema = train_scalar_cb_frozen(model, train_loader, device, args, s_type='eT')
+    estimators['cb_eT'] = {'type': 'cb', 'net': cb_eT, 's_type': 'eT'}
+
+    # 2c. Scalar CB with s=deltaL
+    print("\n--- Scalar CB (s=deltaL) ---")
+    torch.manual_seed(args.seed + 3000)
+    cb_dL, cb_dL_ema = train_scalar_cb_frozen(model, train_loader, device, args, s_type='deltaL')
+    estimators['cb_deltaL'] = {'type': 'cb', 'net': cb_dL, 's_type': 'deltaL'}
+
+    # ----- Step 3: Evaluate -----
+    print(f"\n{'='*60}")
+    print(f"Step 3: Evaluate credit quality")
+    print(f"{'='*60}")
+
+    results = evaluate_credits(model, test_loader, device, estimators, args)
+
+    # State bridge prediction error
+    sb_pred_error = evaluate_state_bridge_pred_error(model, sb, test_loader, device)
+
+    # ----- Print results -----
+    L = args.num_blocks
+    print(f"\n{'='*60}")
+    print(f"RESULTS: Frozen CIFAR Credit Recovery (L={L}, d={args.d_hidden}, seed={args.seed})")
+    print(f"BP reference test accuracy: {bp_acc:.4f}")
+    print(f"{'='*60}")
+
+    # Summary table
+    methods = ['dfa', 'sb_eT', 'cb_eT', 'cb_deltaL']
+    method_labels = {
+        'dfa': 'DFA (random)',
+        'sb_eT': 'State Bridge (eT)',
+        'cb_eT': 'Scalar CB (eT)',
+        'cb_deltaL': 'Scalar CB (deltaL)',
+    }
+
+    print(f"\n{'Method':<25} {'mean Gamma':>12} {'mean rho':>12} {'mean nudge':>12}")
+    print("-" * 65)
+
+    summary = {}
+    for m in methods:
+        r = results[m]
+        mean_gamma = np.mean(r['bp_cosine'])
+        mean_rho = np.mean(r['perturbation_rho'])
+        mean_nudge = np.mean(r['nudging_0.003'])
+        summary[m] = {
+            'mean_gamma': float(mean_gamma),
+            'mean_rho': float(mean_rho),
+            'mean_nudge': float(mean_nudge),
+        }
+        print(f"{method_labels[m]:<25} {mean_gamma:>12.4f} {mean_rho:>12.4f} {mean_nudge:>12.6f}")
+
+    # Per-layer detail
+    print(f"\n--- Per-layer Gamma ---")
+    header = f"{'Layer':<8}"
+    for m in methods:
+        header += f" {method_labels[m]:>16}"
+    print(header)
+    for l in range(L):
+        row = f"  {l:<6}"
+        for m in methods:
+            row += f" {results[m]['bp_cosine'][l]:>16.4f}"
+        print(row)
+
+    print(f"\n--- Per-layer rho ---")
+    print(header)
+    for l in range(L):
+        row = f"  {l:<6}"
+        for m in methods:
+            row += f" {results[m]['perturbation_rho'][l]:>16.4f}"
+        print(row)
+
+    print(f"\n--- Per-layer nudge (eta=0.003) ---")
+    print(header)
+    for l in range(L):
+        row = f"  {l:<6}"
+        for m in methods:
+            row += f" {results[m]['nudging_0.003'][l]:>16.6f}"
+        print(row)
+
+    print(f"\n--- State Bridge prediction error per layer ---")
+    for l in range(L):
+        print(f"  Layer {l}: {sb_pred_error[l]:.6f}")
+
+    # ----- Save all results -----
+    save_data = {
+        'config': {
+            'num_blocks': args.num_blocks,
+            'd_hidden': args.d_hidden,
+            'seed': args.seed,
+            'bp_epochs': args.bp_epochs,
+            'estimator_epochs': args.estimator_epochs,
+            'lr_fb': args.lr_fb,
+            'lam': args.lam,
+            'K': args.K,
+            'sigma_bridge': args.sigma_bridge,
+            'ema_momentum': args.ema_momentum,
+            'term_grad_weight': args.term_grad_weight,
+        },
+        'bp_acc': float(bp_acc),
+        'summary': summary,
+        'per_layer': {},
+        'sb_pred_error': sb_pred_error,
+    }
+
+    for m in methods:
+        save_data['per_layer'][m] = {
+            'bp_cosine': results[m]['bp_cosine'],
+            'perturbation_rho': results[m]['perturbation_rho'],
+            'nudging_0.001': results[m]['nudging_0.001'],
+            'nudging_0.003': results[m]['nudging_0.003'],
+            'nudging_0.01': results[m]['nudging_0.01'],
+        }
+
+    out_path = os.path.join(args.output_dir,
+                            f'frozen_L{args.num_blocks}_d{args.d_hidden}_s{args.seed}.json')
+    with open(out_path, 'w') as f:
+        json.dump(save_data, f, indent=2)
+    print(f"\nResults saved to {out_path}")
+
+    # ----- Judgment -----
+    print(f"\n{'='*60}")
+    print("JUDGMENT")
+    print(f"{'='*60}")
+
+    best_cb = max(summary['cb_eT']['mean_rho'], summary['cb_deltaL']['mean_rho'])
+    dfa_rho = summary['dfa']['mean_rho']
+    best_cb_gamma = max(summary['cb_eT']['mean_gamma'], summary['cb_deltaL']['mean_gamma'])
+    dfa_gamma = summary['dfa']['mean_gamma']
+
+    if best_cb > dfa_rho + 0.02 and best_cb_gamma > dfa_gamma:
+        print("POSITIVE: Scalar CB recovers credit that is clearly better than DFA.")
+        print("  -> Bottleneck is in forward exploitability / local update, not estimator.")
+        print("  -> Next: Phase B (online shallow CIFAR).")
+    elif best_cb > 0.02:
+        print("MARGINAL: Scalar CB shows some signal but not clearly better than DFA.")
+        print("  -> Need more investigation before concluding estimator is the bottleneck.")
+    else:
+        print("NEGATIVE: Scalar CB cannot recover useful credit even on frozen features.")
+        print("  -> Estimator parameterization is the bottleneck.")
+        print("  -> Next: Phase C (direct vector field pilot).")
+
+    return save_data
+
+
+def main():
+    parser = argparse.ArgumentParser(description='Frozen CIFAR Credit Recovery')
+    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('--bp_epochs', type=int, default=100,
+                        help='Epochs to train BP reference')
+    parser.add_argument('--estimator_epochs', type=int, default=100,
+                        help='Epochs to train each estimator on frozen features')
+    parser.add_argument('--lr', type=float, default=1e-3, help='LR for BP reference')
+    parser.add_argument('--lr_fb', type=float, default=1e-3, help='LR for estimators')
+    parser.add_argument('--wd', type=float, default=0.01)
+    parser.add_argument('--lam', type=float, default=0.1)
+    parser.add_argument('--K', type=int, default=4)
+    parser.add_argument('--sigma_bridge', type=float, default=0.05)
+    parser.add_argument('--ema_momentum', type=float, default=0.995)
+    parser.add_argument('--term_grad_weight', type=float, default=1.0)
+    parser.add_argument('--seed', type=int, default=42)
+    parser.add_argument('--gpu', type=int, default=2)
+    parser.add_argument('--output_dir', type=str, default='results/frozen_cifar')
+    parser.add_argument('--retrain_bp', action='store_true')
+    args = parser.parse_args()
+    run_experiment(args)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/experiments/cifar_online_shallow_conditioning.py b/experiments/cifar_online_shallow_conditioning.py
new file mode 100644
index 0000000..cf5444b
--- /dev/null
+++ b/experiments/cifar_online_shallow_conditioning.py
@@ -0,0 +1,717 @@
+"""
+Phase B: Online shallow CIFAR with better conditioning.
+
+Goal: Find a minimal positive-signal regime on real CIFAR-10 with shallow depth.
+Sweep L={4,6}, d={256,512}, methods={DFA, CB_eT, CB_deltaL, SB_eT},
+warmup_ratio={0.0, 0.05, 0.2}, term_grad_weight={1.0, 4.0}.
+
+Single-seed smoke test first. Only expand to 3 seeds for configs with S1>0 and S2>0.
+"""
+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
+import time
+
+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 ValueNet, create_ema_model, update_ema
+from models.state_bridge import StateBridgeNet
+from metrics.credit_metrics import (
+    cosine_similarity_batch, perturbation_correlation, nudging_test,
+    offline_bp_cosine
+)
+
+
+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)
+    train_loader = DataLoader(trainset, batch_size=batch_size, shuffle=True, num_workers=4, pin_memory=True)
+    test_loader = DataLoader(testset, batch_size=batch_size, shuffle=False, num_workers=4, pin_memory=True)
+    return train_loader, test_loader
+
+
+def evaluate(model, test_loader, device):
+    model.eval()
+    correct, total = 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)
+            logits = model(x)
+            correct += (logits.argmax(1) == y).sum().item()
+            total += x.size(0)
+    return correct / total
+
+
+# =============================================================================
+# Training methods
+# =============================================================================
+
+def train_dfa(model, train_loader, test_loader, device, epochs, lr, wd):
+    """DFA training."""
+    d = model.d_hidden
+    num_classes = 10
+    L = model.num_blocks
+
+    Bs = [torch.randn(d, num_classes, device=device) / np.sqrt(num_classes) for _ in range(L)]
+
+    block_opts = [optim.AdamW(block.parameters(), lr=lr, weight_decay=wd) for block 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
+    )
+    all_schedulers = ([optim.lr_scheduler.CosineAnnealingLR(o, T_max=epochs) for o in block_opts]
+                      + [optim.lr_scheduler.CosineAnnealingLR(embed_opt, T_max=epochs),
+                         optim.lr_scheduler.CosineAnnealingLR(head_opt, T_max=epochs)])
+
+    log = {'train_loss': [], 'test_acc': []}
+
+    for epoch in range(1, epochs + 1):
+        model.train()
+        total_loss, correct, total = 0, 0, 0
+        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():
+                logits, hiddens = model(x, return_hidden=True)
+                loss_val = F.cross_entropy(logits, y)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+
+            hL_det = hiddens[-1].detach()
+            logits_out = model.out_head(model.out_ln(hL_det))
+            loss_out = F.cross_entropy(logits_out, y)
+            head_opt.zero_grad()
+            loss_out.backward()
+            head_opt.step()
+
+            for l in range(L):
+                h_l = hiddens[l].detach()
+                a_dfa = (e_T @ Bs[l].T).detach()
+                rms = (a_dfa ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                a_norm = a_dfa / rms
+                f_l = model.blocks[l](h_l)
+                local_loss = (f_l * a_norm).sum(dim=-1).mean()
+                block_opts[l].zero_grad()
+                local_loss.backward()
+                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
+                block_opts[l].step()
+
+            a_0 = (e_T @ Bs[0].T).detach()
+            rms_0 = (a_0 ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+            h0 = model.embed(x)
+            embed_loss = (h0 * (a_0 / rms_0)).sum(dim=-1).mean()
+            embed_opt.zero_grad()
+            embed_loss.backward()
+            embed_opt.step()
+
+            total_loss += loss_val.item() * batch
+            correct += (logits.argmax(1) == y).sum().item()
+            total += batch
+
+        for s in all_schedulers:
+            s.step()
+
+        train_loss = total_loss / total
+        test_acc = evaluate(model, test_loader, device)
+        log['train_loss'].append(train_loss)
+        log['test_acc'].append(test_acc)
+        if epoch % 20 == 0 or epoch == 1:
+            print(f"    [DFA] Ep {epoch}: loss={train_loss:.4f}, test={test_acc:.4f}")
+
+    return log, Bs
+
+
+def train_state_bridge_online(model, train_loader, test_loader, device, epochs, lr, lr_fb, wd):
+    """State bridge online training."""
+    d = model.d_hidden
+    num_classes = 10
+    L = model.num_blocks
+
+    state_pred = StateBridgeNet(
+        d_hidden=d, s_dim=num_classes, time_embed_dim=32, hidden_dim=256, num_layers=3
+    ).to(device)
+
+    block_opts = [optim.AdamW(block.parameters(), lr=lr, weight_decay=wd) for block 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
+    )
+    state_opt = optim.Adam(state_pred.parameters(), lr=lr_fb)
+    all_schedulers = ([optim.lr_scheduler.CosineAnnealingLR(o, T_max=epochs) for o in block_opts]
+                      + [optim.lr_scheduler.CosineAnnealingLR(embed_opt, T_max=epochs),
+                         optim.lr_scheduler.CosineAnnealingLR(head_opt, T_max=epochs)])
+
+    log = {'train_loss': [], 'test_acc': [], 'state_pred_error': []}
+
+    for epoch in range(1, epochs + 1):
+        model.train()
+        state_pred.train()
+        total_loss, correct, total = 0, 0, 0
+        total_se = 0
+
+        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():
+                logits, hiddens = model(x, return_hidden=True)
+                loss_val = F.cross_entropy(logits, y)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+                s = e_T.detach()
+
+            hL_det = hiddens[-1].detach()
+
+            # Train state predictor
+            state_loss = 0.0
+            for l in range(L):
+                h_l_det = hiddens[l].detach()
+                t_l = torch.full((batch,), l / L, device=device)
+                pred_hL = state_pred(h_l_det, t_l, s)
+                target_norm = hL_det.norm(dim=-1, keepdim=True).clamp(min=1.0)
+                state_loss = state_loss + (((pred_hL - hL_det) / target_norm) ** 2).sum(dim=-1).mean()
+            state_loss = state_loss / L
+            state_opt.zero_grad()
+            state_loss.backward()
+            state_opt.step()
+            total_se += state_loss.item() * batch
+
+            # Compute credits
+            credits = []
+            for l in range(L):
+                h_l_det = hiddens[l].detach().requires_grad_(True)
+                t_l = torch.full((batch,), l / L, device=device)
+                pred_hL = state_pred(h_l_det, t_l, s)
+                pred_logits = model.out_head(model.out_ln(pred_hL))
+                pred_loss = F.cross_entropy(pred_logits, y, reduction='sum')
+                a_l = torch.autograd.grad(pred_loss, h_l_det, create_graph=False)[0]
+                credits.append(a_l.detach())
+
+            # Update output head
+            logits_out = model.out_head(model.out_ln(hL_det))
+            loss_out = F.cross_entropy(logits_out, y)
+            head_opt.zero_grad()
+            loss_out.backward()
+            head_opt.step()
+
+            # Update blocks
+            for l in range(L):
+                h_l = hiddens[l].detach()
+                a = credits[l]
+                rms = (a ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                a_norm = a / rms
+                f_l = model.blocks[l](h_l)
+                local_loss = (f_l * a_norm).sum(dim=-1).mean()
+                block_opts[l].zero_grad()
+                local_loss.backward()
+                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
+                block_opts[l].step()
+
+            # Update embedding
+            a_0 = credits[0]
+            rms_0 = (a_0 ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+            h0 = model.embed(x)
+            embed_loss = (h0 * (a_0 / rms_0)).sum(dim=-1).mean()
+            embed_opt.zero_grad()
+            embed_loss.backward()
+            embed_opt.step()
+
+            total_loss += loss_val.item() * batch
+            correct += (logits.argmax(1) == y).sum().item()
+            total += batch
+
+        for sch in all_schedulers:
+            sch.step()
+
+        train_loss = total_loss / total
+        test_acc = evaluate(model, test_loader, device)
+        se = total_se / total
+        log['train_loss'].append(train_loss)
+        log['test_acc'].append(test_acc)
+        log['state_pred_error'].append(se)
+        if epoch % 20 == 0 or epoch == 1:
+            print(f"    [SB] Ep {epoch}: loss={train_loss:.4f}, test={test_acc:.4f}, se={se:.4f}")
+
+    return log, state_pred
+
+
+def train_credit_bridge_online(model, train_loader, test_loader, device,
+                                epochs, lr, lr_fb, wd, s_type='eT',
+                                warmup_ratio=0.2, term_grad_weight=1.0,
+                                lam=0.1, K=4, sigma_bridge=0.05, ema_momentum=0.995):
+    """Credit bridge online training with configurable s_type, warmup, tgw."""
+    d = model.d_hidden
+    num_classes = 10
+    L = model.num_blocks
+    warmup_epochs = max(1, int(epochs * warmup_ratio))
+
+    s_dim = num_classes if s_type == 'eT' else d
+
+    value_net = ValueNet(
+        d_hidden=d, s_dim=s_dim, time_embed_dim=32, hidden_dim=256, num_layers=3
+    ).to(device)
+    value_net_ema = create_ema_model(value_net)
+
+    # DFA fallback for warmup
+    Bs_fallback = [torch.randn(d, num_classes, device=device) / np.sqrt(num_classes) for _ in range(L)]
+
+    block_opts = [optim.AdamW(block.parameters(), lr=lr, weight_decay=wd) for block 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
+    )
+    value_opt = optim.Adam(value_net.parameters(), lr=lr_fb)
+
+    all_schedulers = ([optim.lr_scheduler.CosineAnnealingLR(o, T_max=epochs) for o in block_opts]
+                      + [optim.lr_scheduler.CosineAnnealingLR(embed_opt, T_max=epochs),
+                         optim.lr_scheduler.CosineAnnealingLR(head_opt, T_max=epochs)])
+
+    log = {'train_loss': [], 'test_acc': [], 'value_loss': []}
+
+    for epoch in range(1, epochs + 1):
+        model.train()
+        value_net.train()
+        total_loss, correct, total = 0, 0, 0
+        total_vloss = 0
+
+        if epoch <= warmup_epochs:
+            credit_blend = 0.0
+        else:
+            credit_blend = min(1.0, (epoch - warmup_epochs) / max(1, warmup_epochs))
+
+        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():
+                logits, hiddens = model(x, return_hidden=True)
+                loss_val = F.cross_entropy(logits, y)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+                true_loss = F.cross_entropy(logits, y, reduction='none').detach()
+
+            hL_det = hiddens[-1].detach()
+
+            # Compute s
+            if s_type == 'eT':
+                s = e_T.detach()
+            elif s_type == 'deltaL':
+                hL_req = hL_det.clone().requires_grad_(True)
+                logits_for_s = model.out_head(model.out_ln(hL_req))
+                ce_for_s = F.cross_entropy(logits_for_s, y, reduction='sum')
+                delta_L = torch.autograd.grad(ce_for_s, hL_req, create_graph=False)[0].detach()
+                s = delta_L
+
+            # Train value net
+            t_L = torch.ones(batch, device=device)
+            V_terminal = value_net(hL_det, t_L, s)
+            loss_term = ((V_terminal - true_loss) ** 2).mean()
+
+            loss_tgrad = torch.tensor(0.0, device=device)
+            if term_grad_weight > 0:
+                hL_req = hL_det.clone().requires_grad_(True)
+                V_at_L = value_net(hL_req, t_L, s)
+                grad_V_L = torch.autograd.grad(V_at_L.sum(), hL_req, create_graph=True)[0]
+                hL_req2 = hL_det.clone().requires_grad_(True)
+                logits_tgt = model.out_head(model.out_ln(hL_req2))
+                ce_loss = F.cross_entropy(logits_tgt, y, reduction='sum')
+                a_L_exact = torch.autograd.grad(ce_loss, hL_req2, create_graph=False)[0].detach()
+                loss_tgrad = ((grad_V_L - a_L_exact) ** 2).sum(dim=-1).mean()
+
+            loss_bridge = 0.0
+            for l in range(L):
+                h_l_det = hiddens[l].detach()
+                t_l = torch.full((batch,), l / L, device=device)
+                t_l_next = torch.full((batch,), (l + 1) / L, device=device)
+                V_l = value_net(h_l_det, t_l, s)
+                with torch.no_grad():
+                    h_next_det = hiddens[l + 1].detach()
+                    log_terms = []
+                    for k in range(K):
+                        noise = sigma_bridge * torch.randn_like(h_next_det)
+                        V_next = value_net_ema(h_next_det + noise, t_l_next, s)
+                        log_terms.append(-V_next / lam)
+                    log_stack = torch.stack(log_terms, dim=-1)
+                    V_target = -lam * (torch.logsumexp(log_stack, dim=-1) - np.log(K))
+                loss_bridge = loss_bridge + ((V_l - V_target.detach()) ** 2).mean()
+            loss_bridge = loss_bridge / L
+
+            value_loss = loss_term + loss_bridge + term_grad_weight * loss_tgrad
+            value_opt.zero_grad()
+            value_loss.backward()
+            torch.nn.utils.clip_grad_norm_(value_net.parameters(), 1.0)
+            value_opt.step()
+            update_ema(value_net, value_net_ema, ema_momentum)
+            total_vloss += value_loss.item() * batch
+
+            # Compute CB credits
+            cb_credits = []
+            for l in range(L):
+                h_l_det = hiddens[l].detach().requires_grad_(True)
+                t_l = torch.full((batch,), l / L, device=device)
+                V_l = value_net(h_l_det, t_l, s)
+                a_l = torch.autograd.grad(V_l.sum(), h_l_det, create_graph=False)[0]
+                cb_credits.append(a_l.detach())
+
+            dfa_credits = [(e_T @ Bs_fallback[l].T).detach() for l in range(L)]
+
+            # Blend
+            credits = []
+            for l in range(L):
+                if credit_blend >= 1.0:
+                    a = cb_credits[l]
+                elif credit_blend <= 0.0:
+                    a = dfa_credits[l]
+                else:
+                    cb_rms = (cb_credits[l] ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                    dfa_rms = (dfa_credits[l] ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                    a = credit_blend * (cb_credits[l] / cb_rms) + (1 - credit_blend) * (dfa_credits[l] / dfa_rms)
+                credits.append(a)
+
+            # Update head
+            logits_out = model.out_head(model.out_ln(hL_det))
+            loss_out = F.cross_entropy(logits_out, y)
+            head_opt.zero_grad()
+            loss_out.backward()
+            head_opt.step()
+
+            # Update blocks
+            for l in range(L):
+                h_l = hiddens[l].detach()
+                a = credits[l]
+                rms = (a ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                a_norm = a / rms
+                f_l = model.blocks[l](h_l)
+                local_loss = (f_l * a_norm).sum(dim=-1).mean()
+                block_opts[l].zero_grad()
+                local_loss.backward()
+                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
+                block_opts[l].step()
+
+            # Update embedding
+            a_0 = credits[0]
+            rms_0 = (a_0 ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+            h0 = model.embed(x)
+            embed_loss = (h0 * (a_0 / rms_0)).sum(dim=-1).mean()
+            embed_opt.zero_grad()
+            embed_loss.backward()
+            embed_opt.step()
+
+            total_loss += loss_val.item() * batch
+            correct += (logits.argmax(1) == y).sum().item()
+            total += batch
+
+        for sch in all_schedulers:
+            sch.step()
+
+        train_loss = total_loss / total
+        test_acc = evaluate(model, test_loader, device)
+        vloss = total_vloss / total
+        log['train_loss'].append(train_loss)
+        log['test_acc'].append(test_acc)
+        log['value_loss'].append(vloss)
+        if epoch % 20 == 0 or epoch == 1:
+            phase = "warmup" if epoch <= warmup_epochs else f"blend={credit_blend:.2f}"
+            print(f"    [CB_{s_type}] Ep {epoch} ({phase}): loss={train_loss:.4f}, test={test_acc:.4f}")
+
+    return log, value_net, value_net_ema
+
+
+# =============================================================================
+# Diagnostics
+# =============================================================================
+def compute_diagnostics(model, method_name, test_loader, device,
+                        value_net=None, state_pred=None, dfa_Bs=None, s_type='eT'):
+    """Compute Gamma, rho, nudging per layer."""
+    model.eval()
+    if value_net is not None:
+        value_net.eval()
+    if state_pred is not None:
+        state_pred.eval()
+
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = 10
+
+    # Get one batch
+    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 gradients (evaluation only)
+    logits_bp, hiddens_bp = model(x, return_hidden=True)
+    for l in range(L + 1):
+        hiddens_bp[l].retain_grad()
+    loss_bp = F.cross_entropy(logits_bp, y)
+    loss_bp.backward()
+    bp_grads = {l: hiddens_bp[l].grad.detach().clone() for l in range(L + 1)}
+
+    # Clean forward
+    with torch.no_grad():
+        logits, hiddens = model(x, return_hidden=True)
+        e_T = logits.softmax(dim=-1)
+        e_T[torch.arange(batch), y] -= 1
+        s_eT = e_T.detach()
+
+    hL_det = hiddens[-1].detach()
+
+    # delta_L for deltaL conditioning
+    hL_req = hL_det.clone().requires_grad_(True)
+    logits_for_delta = model.out_head(model.out_ln(hL_req))
+    ce_for_delta = F.cross_entropy(logits_for_delta, y, reduction='sum')
+    delta_L = torch.autograd.grad(ce_for_delta, hL_req, create_graph=False)[0].detach()
+
+    results = {
+        'bp_cosine': [],
+        'perturbation_rho': [],
+        'nudging': {'0.001': [], '0.003': [], '0.01': []},
+    }
+
+    for l in range(L):
+        h_l = hiddens[l].detach()
+        t_l = torch.full((batch,), l / L, device=device)
+
+        if method_name == 'dfa':
+            a_l = (s_eT @ dfa_Bs[l].T).detach()
+        elif method_name == 'state_bridge':
+            h_l_req = h_l.clone().requires_grad_(True)
+            pred_hL = state_pred(h_l_req, t_l, s_eT)
+            pred_logits = model.out_head(model.out_ln(pred_hL))
+            pred_loss = F.cross_entropy(pred_logits, y, reduction='sum')
+            a_l = torch.autograd.grad(pred_loss, h_l_req, create_graph=False)[0].detach()
+        elif method_name.startswith('cb_'):
+            s = s_eT if s_type == 'eT' else delta_L
+            h_l_req = h_l.clone().requires_grad_(True)
+            V_l = value_net(h_l_req, t_l, s)
+            a_l = torch.autograd.grad(V_l.sum(), h_l_req, create_graph=False)[0].detach()
+        else:
+            raise ValueError(f"Unknown method: {method_name}")
+
+        bp_cos = cosine_similarity_batch(a_l, bp_grads[l])
+        results['bp_cosine'].append(float(bp_cos))
+
+        def make_fwd_fn(start_l):
+            def fwd_fn(h):
+                with torch.no_grad():
+                    curr = h
+                    for i in range(start_l, L):
+                        curr = curr + model.blocks[i](curr)
+                    out = model.out_head(model.out_ln(curr))
+                    return F.cross_entropy(out, y, reduction='none')
+            return fwd_fn
+
+        fwd_fn = make_fwd_fn(l)
+        rho = perturbation_correlation(h_l, a_l, fwd_fn, epsilon=1e-3, M=16)
+        results['perturbation_rho'].append(float(rho))
+
+        for eta in [0.001, 0.003, 0.01]:
+            nud = nudging_test(h_l, a_l, fwd_fn, eta=eta)
+            results['nudging'][str(eta)].append(float(nud))
+
+    return results
+
+
+# =============================================================================
+# Single config runner
+# =============================================================================
+def run_config(L, d, method, seed, train_loader, test_loader, device,
+               epochs=100, lr=1e-3, lr_fb=1e-3, wd=0.01,
+               warmup_ratio=0.2, term_grad_weight=1.0,
+               lam=0.1, K=4, sigma_bridge=0.05, ema_momentum=0.995):
+    """Run a single (L, d, method, seed) config and return results."""
+    input_dim = 32 * 32 * 3
+    num_classes = 10
+
+    torch.manual_seed(seed)
+    np.random.seed(seed)
+    torch.cuda.manual_seed_all(seed)
+
+    model = ResidualMLP(input_dim, d, num_classes, L).to(device)
+
+    config_str = f"L={L}, d={d}, method={method}, seed={seed}"
+    if 'cb_' in method:
+        config_str += f", wr={warmup_ratio}, tgw={term_grad_weight}"
+    print(f"\n  --- {config_str} ---")
+
+    if method == 'dfa':
+        log, Bs = train_dfa(model, train_loader, test_loader, device, epochs, lr, wd)
+        diag = compute_diagnostics(model, 'dfa', test_loader, device, dfa_Bs=Bs)
+    elif method == 'sb_eT':
+        log, sp = train_state_bridge_online(model, train_loader, test_loader, device,
+                                             epochs, lr, lr_fb, wd)
+        diag = compute_diagnostics(model, 'state_bridge', test_loader, device, state_pred=sp)
+    elif method == 'cb_eT':
+        log, vnet, _ = train_credit_bridge_online(
+            model, train_loader, test_loader, device, epochs, lr, lr_fb, wd,
+            s_type='eT', warmup_ratio=warmup_ratio, term_grad_weight=term_grad_weight,
+            lam=lam, K=K, sigma_bridge=sigma_bridge, ema_momentum=ema_momentum
+        )
+        diag = compute_diagnostics(model, 'cb_eT', test_loader, device, value_net=vnet, s_type='eT')
+    elif method == 'cb_deltaL':
+        log, vnet, _ = train_credit_bridge_online(
+            model, train_loader, test_loader, device, epochs, lr, lr_fb, wd,
+            s_type='deltaL', warmup_ratio=warmup_ratio, term_grad_weight=term_grad_weight,
+            lam=lam, K=K, sigma_bridge=sigma_bridge, ema_momentum=ema_momentum
+        )
+        diag = compute_diagnostics(model, 'cb_deltaL', test_loader, device, value_net=vnet, s_type='deltaL')
+    else:
+        raise ValueError(f"Unknown method: {method}")
+
+    result = {
+        'method': method,
+        'L': L,
+        'd_hidden': d,
+        'seed': seed,
+        'warmup_ratio': warmup_ratio,
+        'term_grad_weight': term_grad_weight,
+        'test_acc': log['test_acc'][-1],
+        'mean_gamma': float(np.mean(diag['bp_cosine'])),
+        'mean_rho': float(np.mean(diag['perturbation_rho'])),
+        'mean_nudge': float(np.mean(diag['nudging']['0.003'])),
+        'per_layer_gamma': diag['bp_cosine'],
+        'per_layer_rho': diag['perturbation_rho'],
+        'per_layer_nudge': diag['nudging']['0.003'],
+    }
+
+    print(f"    Result: acc={result['test_acc']:.4f}, Gamma={result['mean_gamma']:.4f}, "
+          f"rho={result['mean_rho']:.4f}, nudge={result['mean_nudge']:.6f}")
+
+    return result
+
+
+# =============================================================================
+# Main
+# =============================================================================
+def main():
+    parser = argparse.ArgumentParser(description='Phase B: Online shallow CIFAR conditioning')
+    parser.add_argument('--depths', type=int, nargs='+', default=[4, 6])
+    parser.add_argument('--widths', type=int, nargs='+', default=[256, 512])
+    parser.add_argument('--methods', type=str, nargs='+',
+                        default=['dfa', 'sb_eT', 'cb_eT', 'cb_deltaL'])
+    parser.add_argument('--warmup_ratios', type=float, nargs='+', default=[0.0, 0.05, 0.2])
+    parser.add_argument('--tgws', type=float, nargs='+', default=[1.0, 4.0])
+    parser.add_argument('--epochs', type=int, default=100)
+    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('--lam', type=float, default=0.1)
+    parser.add_argument('--K', type=int, default=4)
+    parser.add_argument('--sigma_bridge', type=float, default=0.05)
+    parser.add_argument('--ema_momentum', type=float, default=0.995)
+    parser.add_argument('--seed', type=int, default=42)
+    parser.add_argument('--batch_size', type=int, default=128)
+    parser.add_argument('--gpu', type=int, default=2)
+    parser.add_argument('--output_dir', type=str, default='results/online_shallow')
+    args = parser.parse_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)
+
+    train_loader, test_loader = get_cifar10(batch_size=args.batch_size)
+
+    all_results = []
+
+    for L in args.depths:
+        for d in args.widths:
+            for method in args.methods:
+                if method in ['dfa', 'sb_eT']:
+                    # No warmup/tgw sweep for DFA and SB
+                    result = run_config(
+                        L, d, method, args.seed, train_loader, test_loader, device,
+                        epochs=args.epochs, lr=args.lr, lr_fb=args.lr_fb, wd=args.wd
+                    )
+                    all_results.append(result)
+                else:
+                    # Sweep warmup and tgw for CB methods
+                    for wr in args.warmup_ratios:
+                        for tgw in args.tgws:
+                            result = run_config(
+                                L, d, method, args.seed, train_loader, test_loader, device,
+                                epochs=args.epochs, lr=args.lr, lr_fb=args.lr_fb, wd=args.wd,
+                                warmup_ratio=wr, term_grad_weight=tgw,
+                                lam=args.lam, K=args.K, sigma_bridge=args.sigma_bridge,
+                                ema_momentum=args.ema_momentum
+                            )
+                            all_results.append(result)
+
+    # Summary table
+    print(f"\n{'='*80}")
+    print("SUMMARY")
+    print(f"{'='*80}")
+
+    # Find DFA baselines for S1, S2 computation
+    dfa_baselines = {}
+    for r in all_results:
+        if r['method'] == 'dfa':
+            dfa_baselines[(r['L'], r['d_hidden'])] = r
+
+    print(f"\n{'Method':<20} {'L':>3} {'d':>4} {'wr':>5} {'tgw':>5} {'Acc':>6} "
+          f"{'Gamma':>7} {'rho':>7} {'nudge':>10} {'S1':>7} {'S2':>7}")
+    print("-" * 95)
+
+    positive_configs = []
+    for r in all_results:
+        key = (r['L'], r['d_hidden'])
+        dfa_ref = dfa_baselines.get(key)
+        S1 = r['mean_gamma'] - (dfa_ref['mean_gamma'] if dfa_ref else 0)
+        S2 = r['mean_rho'] - (dfa_ref['mean_rho'] if dfa_ref else 0)
+        wr_str = f"{r.get('warmup_ratio', '-'):>5.2f}" if r['method'].startswith('cb_') else "    -"
+        tgw_str = f"{r.get('term_grad_weight', '-'):>5.1f}" if r['method'].startswith('cb_') else "    -"
+        print(f"{r['method']:<20} {r['L']:>3} {r['d_hidden']:>4} {wr_str} {tgw_str} "
+              f"{r['test_acc']:>6.4f} {r['mean_gamma']:>7.4f} {r['mean_rho']:>7.4f} "
+              f"{r['mean_nudge']:>10.6f} {S1:>7.4f} {S2:>7.4f}")
+
+        if r['method'].startswith('cb_') and S1 > 0 and S2 > 0:
+            nudge_better = r['mean_nudge'] < (dfa_ref['mean_nudge'] if dfa_ref else 0)
+            positive_configs.append({**r, 'S1': S1, 'S2': S2, 'nudge_better': nudge_better})
+
+    if positive_configs:
+        print(f"\nPOSITIVE CONFIGS (S1>0 AND S2>0):")
+        for pc in positive_configs:
+            print(f"  {pc['method']} L={pc['L']} d={pc['d_hidden']} wr={pc.get('warmup_ratio','-')} "
+                  f"tgw={pc.get('term_grad_weight','-')}: S1={pc['S1']:.4f} S2={pc['S2']:.4f} "
+                  f"nudge_better={pc['nudge_better']}")
+    else:
+        print(f"\nNO POSITIVE CONFIGS FOUND. All CB variants have S1<=0 or S2<=0.")
+
+    # Save
+    out_path = os.path.join(args.output_dir, f'scan_s{args.seed}.json')
+    with open(out_path, 'w') as f:
+        json.dump(all_results, f, indent=2)
+    print(f"\nResults saved to {out_path}")
+
+
+if __name__ == '__main__':
+    main()
diff --git a/experiments/synth_vector_credit.py b/experiments/synth_vector_credit.py
new file mode 100644
index 0000000..14e28e2
--- /dev/null
+++ b/experiments/synth_vector_credit.py
@@ -0,0 +1,708 @@
+"""
+Phase C: Direct Vector Credit Field Pilot.
+
+Compare scalar credit bridge vs direct vector credit field on synthetic best regime.
+Vector field: a_phi(h_l, t_l, s) -> R^d, trained with symmetric finite-difference
+directional targets (no hidden BP anchor).
+
+Loss:
+  L_proj = (1/M) sum_j ( <a_phi(h_l, t_l, s), v_j> - g_j )^2
+  where g_j = [ loss(h_l + eps*v_j) - loss(h_l - eps*v_j) ] / (2*eps)
+
+  L_term = || a_phi(h_L, 1, s) - delta_L ||^2
+
+  L_total = L_term + beta * L_proj
+"""
+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
+import copy
+import time
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+
+from models.value_net import ValueNet, SinusoidalTimeEmbed, create_ema_model, update_ema
+from metrics.credit_metrics import (
+    cosine_similarity_batch, perturbation_correlation, nudging_test,
+    offline_bp_cosine
+)
+
+
+# =============================================================================
+# Synthetic teacher-student (from synth_nonlinearity_ladder.py)
+# =============================================================================
+class TeacherNet(nn.Module):
+    """Fixed teacher with controllable nonlinearity."""
+    def __init__(self, d_hidden, num_classes, num_blocks, alpha=1.0, seed=0):
+        super().__init__()
+        self.d_hidden = d_hidden
+        self.num_blocks = num_blocks
+        self.alpha = alpha
+        rng = torch.Generator().manual_seed(seed)
+        self.Ws = nn.ParameterList()
+        for _ in range(num_blocks):
+            W = torch.randn(d_hidden, d_hidden, generator=rng) * 0.3 / (d_hidden ** 0.5)
+            U, S, Vh = torch.linalg.svd(W, full_matrices=False)
+            S_clamped = S.clamp(max=0.3)
+            W = U @ torch.diag(S_clamped) @ Vh
+            self.Ws.append(nn.Parameter(W, requires_grad=False))
+        self.U = nn.Parameter(torch.randn(num_classes, d_hidden, generator=rng) / (d_hidden ** 0.5),
+                              requires_grad=False)
+
+    def phi(self, z):
+        return (1 - self.alpha) * z + self.alpha * torch.tanh(z)
+
+    def forward(self, x):
+        h = x
+        for W in self.Ws:
+            h = h + self.phi(h @ W.T)
+        return h @ self.U.T
+
+
+class StudentBlock(nn.Module):
+    """Student block with pre-LayerNorm."""
+    def __init__(self, d_hidden, alpha=1.0):
+        super().__init__()
+        self.ln = nn.LayerNorm(d_hidden)
+        self.w = nn.Linear(d_hidden, d_hidden, bias=False)
+        nn.init.normal_(self.w.weight, std=0.01)
+        self.alpha = alpha
+
+    def phi(self, z):
+        return (1 - self.alpha) * z + self.alpha * torch.tanh(z)
+
+    def forward(self, h):
+        return self.w(self.phi(self.ln(h)))
+
+
+class StudentNet(nn.Module):
+    """Student network."""
+    def __init__(self, d_hidden, num_classes, num_blocks, alpha=1.0):
+        super().__init__()
+        self.blocks = nn.ModuleList([StudentBlock(d_hidden, alpha) for _ in range(num_blocks)])
+        self.out_head = nn.Linear(d_hidden, num_classes)
+        self.d_hidden = d_hidden
+        self.num_blocks = num_blocks
+
+    def forward(self, x, return_hidden=False):
+        h = x
+        hiddens = [h] if return_hidden else None
+        for block in self.blocks:
+            f = block(h)
+            h = h + f
+            if return_hidden:
+                hiddens.append(h)
+        logits = self.out_head(h)
+        if return_hidden:
+            return logits, hiddens
+        return logits
+
+    def forward_from_layer(self, h, start_layer):
+        for i in range(start_layer, self.num_blocks):
+            f = self.blocks[i](h)
+            h = h + f
+        return self.out_head(h)
+
+
+# =============================================================================
+# Vector Credit Field Network
+# =============================================================================
+class VectorCreditNet(nn.Module):
+    """
+    Direct vector credit field: a_phi(h_l, t_l, s) -> R^d.
+    Output is d-dimensional credit vector directly.
+    """
+    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):
+        """Returns credit vector (batch, d_hidden)."""
+        h_normed = self.ln(h)
+        t_emb = self.time_embed(t)
+        inp = torch.cat([h_normed, t_emb, s], dim=-1)
+        return self.net(inp)
+
+
+# =============================================================================
+# Training functions
+# =============================================================================
+def generate_batch(teacher, d_hidden, num_classes, batch_size, device):
+    """Generate synthetic data from teacher."""
+    x = torch.randn(batch_size, d_hidden, device=device)
+    with torch.no_grad():
+        teacher_logits = teacher(x)
+        y = teacher_logits.argmax(dim=-1)
+    return x, y
+
+
+def train_scalar_cb(model, teacher, device, args):
+    """Train scalar credit bridge (current method, as baseline)."""
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = args.num_classes
+
+    value_net = ValueNet(d_hidden=d, s_dim=num_classes, time_embed_dim=32,
+                         hidden_dim=256, num_layers=3).to(device)
+    value_net_ema = create_ema_model(value_net)
+
+    Bs = [torch.randn(d, num_classes, device=device) / np.sqrt(num_classes) for _ in range(L)]
+
+    block_opts = [optim.AdamW(b.parameters(), lr=args.lr, weight_decay=0.01) for b in model.blocks]
+    head_opt = optim.AdamW(model.out_head.parameters(), lr=args.lr, weight_decay=0.01)
+    value_opt = optim.Adam(value_net.parameters(), lr=args.lr_fb)
+
+    warmup_epochs = max(1, int(args.epochs * args.warmup_ratio))
+    log = {'train_loss': [], 'test_acc': []}
+
+    for epoch in range(1, args.epochs + 1):
+        model.train()
+        value_net.train()
+
+        if epoch <= warmup_epochs:
+            credit_blend = 0.0
+        else:
+            credit_blend = min(1.0, (epoch - warmup_epochs) / max(1, warmup_epochs))
+
+        total_loss, correct, total = 0, 0, 0
+        for _ in range(args.steps_per_epoch):
+            x, y = generate_batch(teacher, d, num_classes, args.batch_size, device)
+            batch = x.size(0)
+
+            with torch.no_grad():
+                logits, hiddens = model(x, return_hidden=True)
+                loss_val = F.cross_entropy(logits, y)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+                s = e_T.detach()
+                true_loss = F.cross_entropy(logits, y, reduction='none').detach()
+
+            hL_det = hiddens[-1].detach()
+
+            # Train value net
+            t_L = torch.ones(batch, device=device)
+            V_term = value_net(hL_det, t_L, s)
+            loss_term = ((V_term - true_loss) ** 2).mean()
+
+            # Terminal gradient matching
+            hL_req = hL_det.clone().requires_grad_(True)
+            V_at_L = value_net(hL_req, t_L, s)
+            grad_V_L = torch.autograd.grad(V_at_L.sum(), hL_req, create_graph=True)[0]
+            hL_req2 = hL_det.clone().requires_grad_(True)
+            logits_tgt = model.out_head(hL_req2)
+            ce = F.cross_entropy(logits_tgt, y, reduction='sum')
+            a_L_exact = torch.autograd.grad(ce, hL_req2, create_graph=False)[0].detach()
+            loss_tgrad = ((grad_V_L - a_L_exact) ** 2).sum(dim=-1).mean()
+
+            # Bridge consistency
+            loss_bridge = 0.0
+            for l in range(L):
+                h_l_det = hiddens[l].detach()
+                t_l = torch.full((batch,), l / L, device=device)
+                t_next = torch.full((batch,), (l + 1) / L, device=device)
+                V_l = value_net(h_l_det, t_l, s)
+                with torch.no_grad():
+                    h_next = hiddens[l + 1].detach()
+                    log_terms = []
+                    for k in range(args.K):
+                        noise = args.sigma_bridge * torch.randn_like(h_next)
+                        V_next = value_net_ema(h_next + noise, t_next, s)
+                        log_terms.append(-V_next / args.lam)
+                    log_stack = torch.stack(log_terms, dim=-1)
+                    V_target = -args.lam * (torch.logsumexp(log_stack, dim=-1) - np.log(args.K))
+                loss_bridge += ((V_l - V_target.detach()) ** 2).mean()
+            loss_bridge /= L
+
+            vloss = loss_term + loss_bridge + args.term_grad_weight * loss_tgrad
+            value_opt.zero_grad()
+            vloss.backward()
+            torch.nn.utils.clip_grad_norm_(value_net.parameters(), 1.0)
+            value_opt.step()
+            update_ema(value_net, value_net_ema, args.ema_momentum)
+
+            # Compute credits
+            cb_credits = []
+            for l in range(L):
+                h_l_det = hiddens[l].detach().requires_grad_(True)
+                t_l = torch.full((batch,), l / L, device=device)
+                V_l = value_net(h_l_det, t_l, s)
+                a_l = torch.autograd.grad(V_l.sum(), h_l_det, create_graph=False)[0]
+                cb_credits.append(a_l.detach())
+
+            dfa_credits = [(e_T @ Bs[l].T).detach() for l in range(L)]
+
+            credits = []
+            for l in range(L):
+                if credit_blend >= 1.0:
+                    credits.append(cb_credits[l])
+                elif credit_blend <= 0.0:
+                    credits.append(dfa_credits[l])
+                else:
+                    cb_rms = (cb_credits[l] ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                    dfa_rms = (dfa_credits[l] ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                    credits.append(credit_blend * cb_credits[l] / cb_rms + (1 - credit_blend) * dfa_credits[l] / dfa_rms)
+
+            # Update head
+            logits_out = model.out_head(hL_det)
+            loss_out = F.cross_entropy(logits_out, y)
+            head_opt.zero_grad()
+            loss_out.backward()
+            head_opt.step()
+
+            # Update blocks
+            for l in range(L):
+                h_l = hiddens[l].detach()
+                a = credits[l]
+                rms = (a ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                f_l = model.blocks[l](h_l)
+                local_loss = (f_l * (a / rms)).sum(dim=-1).mean()
+                block_opts[l].zero_grad()
+                local_loss.backward()
+                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
+                block_opts[l].step()
+
+            total_loss += loss_val.item() * batch
+            correct += (logits.argmax(1) == y).sum().item()
+            total += batch
+
+        if epoch % 10 == 0 or epoch == 1:
+            acc = correct / total
+            print(f"    [scalar_cb] Ep {epoch}: loss={total_loss/total:.4f}, acc={acc:.4f}")
+
+    return value_net
+
+
+def train_vector_field(model, teacher, device, args, M=4):
+    """
+    Train direct vector credit field with perturbation-based targets.
+    No hidden BP anchor.
+    """
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = args.num_classes
+
+    vector_net = VectorCreditNet(d_hidden=d, s_dim=num_classes, time_embed_dim=32,
+                                  hidden_dim=256, num_layers=3).to(device)
+
+    Bs = [torch.randn(d, num_classes, device=device) / np.sqrt(num_classes) for _ in range(L)]
+
+    block_opts = [optim.AdamW(b.parameters(), lr=args.lr, weight_decay=0.01) for b in model.blocks]
+    head_opt = optim.AdamW(model.out_head.parameters(), lr=args.lr, weight_decay=0.01)
+    vec_opt = optim.Adam(vector_net.parameters(), lr=args.lr_fb)
+
+    warmup_epochs = max(1, int(args.epochs * args.warmup_ratio))
+    eps = args.pert_eps
+    beta = args.pert_beta
+
+    for epoch in range(1, args.epochs + 1):
+        model.train()
+        vector_net.train()
+
+        if epoch <= warmup_epochs:
+            credit_blend = 0.0
+        else:
+            credit_blend = min(1.0, (epoch - warmup_epochs) / max(1, warmup_epochs))
+
+        total_loss, correct, total = 0, 0, 0
+        total_vloss = 0
+
+        for _ in range(args.steps_per_epoch):
+            x, y = generate_batch(teacher, d, num_classes, args.batch_size, device)
+            batch = x.size(0)
+
+            with torch.no_grad():
+                logits, hiddens = model(x, return_hidden=True)
+                loss_val = F.cross_entropy(logits, y)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+                s = e_T.detach()
+
+            hL_det = hiddens[-1].detach()
+
+            # Terminal matching: a_phi(h_L, 1, s) = delta_L
+            t_L = torch.ones(batch, device=device)
+            a_terminal = vector_net(hL_det, t_L, s)
+            # delta_L = grad_{h_L} CE (output-layer-local)
+            hL_req = hL_det.clone().requires_grad_(True)
+            logits_tgt = model.out_head(hL_req)
+            ce = F.cross_entropy(logits_tgt, y, reduction='sum')
+            delta_L = torch.autograd.grad(ce, hL_req, create_graph=False)[0].detach()
+            loss_term = ((a_terminal - delta_L) ** 2).sum(dim=-1).mean()
+
+            # Perturbation-based directional targets for all layers
+            loss_proj = torch.tensor(0.0, device=device)
+            for l in range(L):
+                h_l_det = hiddens[l].detach()
+                t_l = torch.full((batch,), l / L, device=device)
+                a_l = vector_net(h_l_det, t_l, s)
+
+                # Compute directional targets using symmetric finite difference
+                layer_proj_loss = 0.0
+                for _ in range(M):
+                    v = torch.randn_like(h_l_det)
+                    v = v / (v.norm(dim=-1, keepdim=True) + 1e-8)
+
+                    # Forward from perturbed h_l
+                    with torch.no_grad():
+                        logits_plus = model.forward_from_layer(h_l_det + eps * v, l)
+                        loss_plus = F.cross_entropy(logits_plus, y, reduction='none')
+                        logits_minus = model.forward_from_layer(h_l_det - eps * v, l)
+                        loss_minus = F.cross_entropy(logits_minus, y, reduction='none')
+                        g_j = (loss_plus - loss_minus) / (2 * eps)  # (batch,)
+
+                    # Predicted directional derivative
+                    pred_j = (a_l * v).sum(dim=-1)  # (batch,)
+                    layer_proj_loss = layer_proj_loss + ((pred_j - g_j.detach()) ** 2).mean()
+
+                loss_proj = loss_proj + layer_proj_loss / M
+            loss_proj = loss_proj / L
+
+            vec_loss = loss_term + beta * loss_proj
+            vec_opt.zero_grad()
+            vec_loss.backward()
+            torch.nn.utils.clip_grad_norm_(vector_net.parameters(), 1.0)
+            vec_opt.step()
+            total_vloss += vec_loss.item() * batch
+
+            # Compute credits for block updates
+            with torch.no_grad():
+                vec_credits = []
+                for l in range(L):
+                    h_l_det = hiddens[l].detach()
+                    t_l = torch.full((batch,), l / L, device=device)
+                    a_l = vector_net(h_l_det, t_l, s)
+                    vec_credits.append(a_l.detach())
+
+            dfa_credits = [(e_T @ Bs[l].T).detach() for l in range(L)]
+
+            credits = []
+            for l in range(L):
+                if credit_blend >= 1.0:
+                    credits.append(vec_credits[l])
+                elif credit_blend <= 0.0:
+                    credits.append(dfa_credits[l])
+                else:
+                    vc_rms = (vec_credits[l] ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                    dfa_rms = (dfa_credits[l] ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                    credits.append(credit_blend * vec_credits[l] / vc_rms + (1 - credit_blend) * dfa_credits[l] / dfa_rms)
+
+            # Update head
+            logits_out = model.out_head(hL_det)
+            loss_out = F.cross_entropy(logits_out, y)
+            head_opt.zero_grad()
+            loss_out.backward()
+            head_opt.step()
+
+            # Update blocks
+            for l in range(L):
+                h_l = hiddens[l].detach()
+                a = credits[l]
+                rms = (a ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                f_l = model.blocks[l](h_l)
+                local_loss = (f_l * (a / rms)).sum(dim=-1).mean()
+                block_opts[l].zero_grad()
+                local_loss.backward()
+                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
+                block_opts[l].step()
+
+            total_loss += loss_val.item() * batch
+            correct += (logits.argmax(1) == y).sum().item()
+            total += batch
+
+        if epoch % 10 == 0 or epoch == 1:
+            acc = correct / total
+            print(f"    [vec_M={M}] Ep {epoch}: loss={total_loss/total:.4f}, acc={acc:.4f}, "
+                  f"vloss={total_vloss/total:.6f}")
+
+    return vector_net
+
+
+def train_dfa(model, teacher, device, args):
+    """DFA baseline for comparison."""
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = args.num_classes
+
+    Bs = [torch.randn(d, num_classes, device=device) / np.sqrt(num_classes) for _ in range(L)]
+
+    block_opts = [optim.AdamW(b.parameters(), lr=args.lr, weight_decay=0.01) for b in model.blocks]
+    head_opt = optim.AdamW(model.out_head.parameters(), lr=args.lr, weight_decay=0.01)
+
+    for epoch in range(1, args.epochs + 1):
+        model.train()
+        total_loss, correct, total = 0, 0, 0
+        for _ in range(args.steps_per_epoch):
+            x, y = generate_batch(teacher, d, num_classes, args.batch_size, device)
+            batch = x.size(0)
+
+            with torch.no_grad():
+                logits, hiddens = model(x, return_hidden=True)
+                loss_val = F.cross_entropy(logits, y)
+                e_T = logits.softmax(dim=-1)
+                e_T[torch.arange(batch), y] -= 1
+
+            hL_det = hiddens[-1].detach()
+            logits_out = model.out_head(hL_det)
+            loss_out = F.cross_entropy(logits_out, y)
+            head_opt.zero_grad()
+            loss_out.backward()
+            head_opt.step()
+
+            for l in range(L):
+                h_l = hiddens[l].detach()
+                a = (e_T @ Bs[l].T).detach()
+                rms = (a ** 2).mean(dim=-1, keepdim=True).sqrt() + 1e-6
+                f_l = model.blocks[l](h_l)
+                local_loss = (f_l * (a / rms)).sum(dim=-1).mean()
+                block_opts[l].zero_grad()
+                local_loss.backward()
+                torch.nn.utils.clip_grad_norm_(model.blocks[l].parameters(), 1.0)
+                block_opts[l].step()
+
+            total_loss += loss_val.item() * batch
+            correct += (logits.argmax(1) == y).sum().item()
+            total += batch
+
+        if epoch % 10 == 0 or epoch == 1:
+            print(f"    [DFA] Ep {epoch}: loss={total_loss/total:.4f}, acc={correct/total:.4f}")
+
+    return Bs
+
+
+# =============================================================================
+# Diagnostics
+# =============================================================================
+def compute_diagnostics(model, teacher, device, method_name, args,
+                        value_net=None, vector_net=None, dfa_Bs=None):
+    """Compute Gamma, rho, nudging per layer."""
+    model.eval()
+    if value_net is not None:
+        value_net.eval()
+    if vector_net is not None:
+        vector_net.eval()
+
+    d = model.d_hidden
+    L = model.num_blocks
+    num_classes = args.num_classes
+
+    x, y = generate_batch(teacher, d, num_classes, 512, device)
+    batch = x.size(0)
+
+    # BP gradients (evaluation only)
+    h = x.detach().requires_grad_(True)
+    hiddens_bp = [h]
+    for block in model.blocks:
+        f = block(hiddens_bp[-1])
+        h_next = hiddens_bp[-1] + f
+        hiddens_bp.append(h_next)
+    logits_bp = model.out_head(hiddens_bp[-1])
+    loss_bp = F.cross_entropy(logits_bp, y)
+    grads = torch.autograd.grad(loss_bp, hiddens_bp, retain_graph=False)
+    bp_grads = {l: grads[l].detach().clone() for l in range(L + 1)}
+
+    with torch.no_grad():
+        logits, hiddens = model(x, return_hidden=True)
+        e_T = logits.softmax(dim=-1)
+        e_T[torch.arange(batch), y] -= 1
+        s = e_T.detach()
+
+    results = {'bp_cosine': [], 'perturbation_rho': [], 'nudging': []}
+
+    for l in range(L):
+        h_l = hiddens[l].detach()
+        t_l = torch.full((batch,), l / L, device=device)
+
+        if method_name == 'dfa':
+            a_l = (s @ dfa_Bs[l].T).detach()
+        elif method_name == 'scalar_cb':
+            h_l_req = h_l.clone().requires_grad_(True)
+            V_l = value_net(h_l_req, t_l, s)
+            a_l = torch.autograd.grad(V_l.sum(), h_l_req, create_graph=False)[0].detach()
+        elif method_name.startswith('vector'):
+            a_l = vector_net(h_l, t_l, s).detach()
+        else:
+            raise ValueError(f"Unknown: {method_name}")
+
+        bp_cos = cosine_similarity_batch(a_l, bp_grads[l])
+        results['bp_cosine'].append(float(bp_cos))
+
+        def make_fwd_fn(start_l):
+            def fwd_fn(h):
+                with torch.no_grad():
+                    logits = model.forward_from_layer(h, start_l)
+                    return F.cross_entropy(logits, y, reduction='none')
+            return fwd_fn
+
+        fwd_fn = make_fwd_fn(l)
+        rho = perturbation_correlation(h_l, a_l, fwd_fn, epsilon=1e-3, M=32)
+        results['perturbation_rho'].append(float(rho))
+
+        nud = nudging_test(h_l, a_l, fwd_fn, eta=0.003)
+        results['nudging'].append(float(nud))
+
+    return results
+
+
+# =============================================================================
+# Main
+# =============================================================================
+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)
+
+    all_results = []
+
+    for L in args.depths:
+        for seed in args.seeds:
+            print(f"\n{'='*60}")
+            print(f"L={L}, seed={seed}")
+            print(f"{'='*60}")
+
+            torch.manual_seed(seed)
+            np.random.seed(seed)
+            torch.cuda.manual_seed_all(seed)
+
+            teacher = TeacherNet(args.d_hidden, args.num_classes, L,
+                                 alpha=args.alpha, seed=seed * 1000).to(device)
+
+            # --- DFA ---
+            print("\n  --- DFA ---")
+            torch.manual_seed(seed)
+            model_dfa = StudentNet(args.d_hidden, args.num_classes, L, alpha=args.alpha).to(device)
+            Bs = train_dfa(model_dfa, teacher, device, args)
+            diag_dfa = compute_diagnostics(model_dfa, teacher, device, 'dfa', args, dfa_Bs=Bs)
+            r_dfa = {
+                'method': 'dfa', 'L': L, 'seed': seed,
+                'mean_gamma': float(np.mean(diag_dfa['bp_cosine'])),
+                'mean_rho': float(np.mean(diag_dfa['perturbation_rho'])),
+                'mean_nudge': float(np.mean(diag_dfa['nudging'])),
+                'per_layer_gamma': diag_dfa['bp_cosine'],
+                'per_layer_rho': diag_dfa['perturbation_rho'],
+            }
+            print(f"    Result: Gamma={r_dfa['mean_gamma']:.4f}, rho={r_dfa['mean_rho']:.4f}")
+            all_results.append(r_dfa)
+
+            # --- Scalar CB ---
+            print("\n  --- Scalar CB ---")
+            torch.manual_seed(seed)
+            model_cb = StudentNet(args.d_hidden, args.num_classes, L, alpha=args.alpha).to(device)
+            vnet = train_scalar_cb(model_cb, teacher, device, args)
+            diag_cb = compute_diagnostics(model_cb, teacher, device, 'scalar_cb', args, value_net=vnet)
+            r_cb = {
+                'method': 'scalar_cb', 'L': L, 'seed': seed,
+                'mean_gamma': float(np.mean(diag_cb['bp_cosine'])),
+                'mean_rho': float(np.mean(diag_cb['perturbation_rho'])),
+                'mean_nudge': float(np.mean(diag_cb['nudging'])),
+                'per_layer_gamma': diag_cb['bp_cosine'],
+                'per_layer_rho': diag_cb['perturbation_rho'],
+            }
+            print(f"    Result: Gamma={r_cb['mean_gamma']:.4f}, rho={r_cb['mean_rho']:.4f}")
+            all_results.append(r_cb)
+
+            # --- Vector Field M=4 ---
+            for M in args.M_values:
+                print(f"\n  --- Vector Field M={M} ---")
+                torch.manual_seed(seed)
+                model_vec = StudentNet(args.d_hidden, args.num_classes, L, alpha=args.alpha).to(device)
+                vec_net = train_vector_field(model_vec, teacher, device, args, M=M)
+                diag_vec = compute_diagnostics(model_vec, teacher, device, f'vector_M{M}', args,
+                                                vector_net=vec_net)
+                r_vec = {
+                    'method': f'vector_M{M}', 'L': L, 'seed': seed, 'M': M,
+                    'mean_gamma': float(np.mean(diag_vec['bp_cosine'])),
+                    'mean_rho': float(np.mean(diag_vec['perturbation_rho'])),
+                    'mean_nudge': float(np.mean(diag_vec['nudging'])),
+                    'per_layer_gamma': diag_vec['bp_cosine'],
+                    'per_layer_rho': diag_vec['perturbation_rho'],
+                }
+                print(f"    Result: Gamma={r_vec['mean_gamma']:.4f}, rho={r_vec['mean_rho']:.4f}")
+                all_results.append(r_vec)
+
+    # Summary
+    print(f"\n{'='*80}")
+    print("SUMMARY")
+    print(f"{'='*80}")
+    print(f"{'Method':<20} {'L':>3} {'seed':>5} {'Gamma':>8} {'rho':>8} {'nudge':>10}")
+    print("-" * 60)
+    for r in all_results:
+        print(f"{r['method']:<20} {r['L']:>3} {r['seed']:>5} {r['mean_gamma']:>8.4f} "
+              f"{r['mean_rho']:>8.4f} {r['mean_nudge']:>10.6f}")
+
+    # Save
+    out_path = os.path.join(args.output_dir, 'results.json')
+    with open(out_path, 'w') as f:
+        json.dump(all_results, f, indent=2)
+    print(f"\nResults saved to {out_path}")
+
+    # Compare vector field vs scalar CB
+    print(f"\n{'='*60}")
+    print("COMPARISON: Vector Field vs Scalar CB")
+    print(f"{'='*60}")
+    for L in args.depths:
+        for seed in args.seeds:
+            cb_r = [r for r in all_results if r['method'] == 'scalar_cb' and r['L'] == L and r['seed'] == seed]
+            if not cb_r:
+                continue
+            cb_r = cb_r[0]
+            for M in args.M_values:
+                vec_r = [r for r in all_results if r['method'] == f'vector_M{M}' and r['L'] == L and r['seed'] == seed]
+                if not vec_r:
+                    continue
+                vec_r = vec_r[0]
+                delta_gamma = vec_r['mean_gamma'] - cb_r['mean_gamma']
+                delta_rho = vec_r['mean_rho'] - cb_r['mean_rho']
+                print(f"  L={L} seed={seed} M={M}: delta_Gamma={delta_gamma:+.4f}, delta_rho={delta_rho:+.4f}")
+                if delta_rho >= 0.05 or delta_gamma >= 0.05:
+                    print(f"    -> SIGNIFICANT IMPROVEMENT")
+                elif delta_rho > 0 and delta_gamma > 0:
+                    print(f"    -> Modest improvement")
+                else:
+                    print(f"    -> No clear improvement")
+
+
+def main():
+    parser = argparse.ArgumentParser(description='Phase C: Vector Credit Field Pilot')
+    parser.add_argument('--d_hidden', type=int, default=128)
+    parser.add_argument('--num_classes', type=int, default=10)
+    parser.add_argument('--alpha', type=float, default=1.0)
+    parser.add_argument('--depths', type=int, nargs='+', default=[4, 8])
+    parser.add_argument('--M_values', type=int, nargs='+', default=[4, 8])
+    parser.add_argument('--epochs', type=int, default=80)
+    parser.add_argument('--steps_per_epoch', type=int, default=50)
+    parser.add_argument('--batch_size', type=int, default=256)
+    parser.add_argument('--lr', type=float, default=1e-3)
+    parser.add_argument('--lr_fb', type=float, default=1e-3)
+    parser.add_argument('--warmup_ratio', type=float, default=0.05)
+    parser.add_argument('--term_grad_weight', type=float, default=1.0)
+    parser.add_argument('--lam', type=float, default=0.1)
+    parser.add_argument('--K', type=int, default=4)
+    parser.add_argument('--sigma_bridge', type=float, default=0.05)
+    parser.add_argument('--ema_momentum', type=float, default=0.995)
+    parser.add_argument('--pert_eps', type=float, default=1e-3)
+    parser.add_argument('--pert_beta', type=float, default=1.0)
+    parser.add_argument('--seeds', type=int, nargs='+', default=[42, 123, 456])
+    parser.add_argument('--gpu', type=int, default=3)
+    parser.add_argument('--output_dir', type=str, default='results/vector_credit_pilot')
+    args = parser.parse_args()
+    run_experiment(args)
+
+
+if __name__ == '__main__':
+    main()