1import itertools, json
  2import numpy as np
  3
  4SEED = 123
  5ORDERS = list(itertools.permutations(range(4)))
  6
  7def rho(a):
  8    return float(np.max(np.abs(np.linalg.eigvals(a))))
  9
 10def block_maps(lam, gamma):
 11    # A local error model: each update changes one coordinate block. The directed
 12    # sensitivities are deliberately nonsymmetric, as neural block Jacobians need not
 13    # be symmetric. At lam=0 all blocks are decoupled and order must be irrelevant.
 14    c = np.array([[0, 1.0, -.6, .3], [-.8, 0, .9, -.4],
 15                  [.5, -.7, 0, 1.1], [-.3, .8, -.9, 0]], dtype=float)
 16    out = []
 17    for a in range(4):
 18        t = np.eye(4)
 19        t[a, :] -= gamma * (np.eye(4)[a, :] + lam * c[a, :])
 20        out.append(t)
 21    return out
 22
 23def composed(lam, gamma, order):
 24    j = np.eye(4)
 25    for b in order:                         # order[0] acts first
 26        j = block_maps(lam, gamma)[b] @ j
 27    return j
 28
 29def metrics(lam, gamma):
 30    vals = np.array([rho(composed(lam, gamma, o)) for o in ORDERS])
 31    return vals
 32
 33def observed_ratio(j, e0=None, n=100):
 34    e = np.ones(j.shape[0]) if e0 is None else np.asarray(e0, float).copy()
 35    rs = []
 36    for _ in range(n):
 37        en = j @ e
 38        rs.append(np.linalg.norm(en) / np.linalg.norm(e))
 39        e = en
 40    return float(np.mean(rs[-20:]))
 41
 42def crossing(lam, order, criterion, lo=.001, hi=2.0):
 43    # Find the gamma boundary predicted by rho(J)=1. For this family rho is
 44    # monotonic over the relevant interval; return None if no crossing exists.
 45    f = lambda g: criterion(composed(lam, g, order)) - 1.0
 46    if f(lo) >= 0: return lo
 47    if f(hi) < 0: return None
 48    for _ in range(60):
 49        mid = (lo + hi) / 2
 50        if f(mid) < 0: lo = mid
 51        else: hi = mid
 52    return (lo + hi) / 2
 53
 54def toy_verification():
 55    # Prediction 1: at zero coupling, all order radii coincide exactly.
 56    # Prediction 2: increasing directed coupling increases order spread.
 57    coupling = []
 58    for lam in [0, .1, .2, .3, .4, .5, .6, .7, .8]:
 59        vals = metrics(lam, .42)
 60        coupling.append({'lambda': lam, 'best_rho': float(vals.min()),
 61                         'worst_rho': float(vals.max()),
 62                         'order_gap': float(vals.max()-vals.min()),
 63                         'best_order': list(ORDERS[int(vals.argmin())])})
 64    # Prediction 3: rho=1 is the stability boundary; compare predicted crossing
 65    # against direct finite-trajectory growth classification.
 66    stability = []
 67    for lam in [.2, .4, .6, .8]:
 68        vals = metrics(lam, .42)
 69        order = ORDERS[int(vals.argmin())]
 70        pred = crossing(lam, order, rho)
 71        test_gammas = np.linspace(.05, 1.5, 30)
 72        observed = None
 73        for g in test_gammas:
 74            j = composed(lam, g, order)
 75            e = np.ones(4)
 76            for _ in range(80): e = j @ e
 77            if np.linalg.norm(e) > 1.05: # finite-run growth, empirical boundary
 78                observed = float(g); break
 79        stability.append({'lambda': lam, 'order': list(order),
 80                          'predicted_rho1_gamma': pred,
 81                          'observed_first_growth_gamma': observed})
 82    # Direct asymptotic-ratio check for the selected composition.
 83    lam, gamma = .5, .42
 84    vals = metrics(lam, gamma); order = ORDERS[int(vals.argmin())]
 85    j = composed(lam, gamma, order)
 86    ratio = {'lambda': lam, 'gamma': gamma, 'order': list(order),
 87             'predicted_rho': rho(j),
 88             'observed_tail_ratio': observed_ratio(j, [.8, -.4, .6, 1.1])}
 89    return {'coupling_sweep': coupling, 'stability_sweep': stability,
 90            'ratio_check': ratio}
 91
 92def neural_experiment():
 93    import torch
 94    torch.manual_seed(SEED); np.random.seed(SEED)
 95    device = 'cuda' if torch.cuda.is_available() else 'cpu'
 96    try:
 97        x = torch.linspace(-2, 2, 256, device=device).unsqueeze(1)
 98        y = torch.sin(3*x) + .15*x*x
 99        def run(mode):
100            torch.manual_seed(SEED)
101            net = torch.nn.Sequential(torch.nn.Linear(1,16), torch.nn.Tanh(),
102                                      torch.nn.Linear(16,1)).to(device)
103            blocks = [[net[0].weight, net[0].bias], [net[2].weight, net[2].bias]]
104            order = [0, 1]; losses=[]; orders=[]; lr=.035
105            for step in range(180):
106                if mode == 'adaptive' and step % 15 == 0:
107                    # A cheap empirical cross-block proxy: test each candidate by
108                    # one block update and score the resulting loss.
109                    scores=[]
110                    for candidate in ([0,1],[1,0]):
111                        saved=[p.detach().clone() for p in net.parameters()]
112                        bi=candidate[0]; net.zero_grad(); q=((net(x)-y)**2).mean(); q.backward()
113                        with torch.no_grad():
114                            for p in blocks[bi]: p -= lr*p.grad
115                        scores.append(float(((net(x)-y)**2).mean()))
116                        with torch.no_grad():
117                            for p,s in zip(net.parameters(),saved): p.copy_(s)
118                    order = [0,1] if scores[0] <= scores[1] else [1,0]
119                if mode == 'sim':
120                    net.zero_grad(); loss=((net(x)-y)**2).mean(); loss.backward()
121                    with torch.no_grad():
122                        for p in net.parameters(): p -= lr*p.grad
123                else:
124                    for bi in order:
125                        net.zero_grad(); loss=((net(x)-y)**2).mean(); loss.backward()
126                        with torch.no_grad():
127                            for p in blocks[bi]: p -= lr*p.grad
128                orders.append(tuple(order)); losses.append(float(((net(x)-y)**2).mean()))
129            return losses, orders
130        out={m:run(m) for m in ['sim','fixed','adaptive']}
131        return {'device':device, 'final_loss':{k:v[0][-1] for k,v in out.items()},
132                'loss_at_30':{k:v[0][29] for k,v in out.items()},
133                'adaptive_order_counts':{str(o):out['adaptive'][1].count(o) for o in [(0,1),(1,0)]}}
134    except Exception as e:
135        return {'device':'cpu-fallback','error':repr(e)}
136
137def main():
138    result={'seed':SEED, 'toy':toy_verification(), 'neural':neural_experiment()}
139    with open('results.json','w') as f: json.dump(result,f,indent=2)
140    print(json.dumps(result,indent=2))
141if __name__ == '__main__': main()