1import json, math, time
 2import numpy as np
 3from spectral_edge import SpectralEdgeController
 4
 5SEED=2744
 6
 7def toy_sweeps():
 8    rng=np.random.default_rng(SEED); n=160
 9    A=rng.normal(size=(n,n)); W=(A+A.T)/2/np.sqrt(n)
10    lam=float(np.linalg.eigvalsh(W)[-1])
11    ratios=np.array([.70,.85,.95,.99,1.01,1.10]); T=260
12    v=rng.normal(size=n); v/=np.linalg.norm(v)
13    coeff=abs(np.linalg.eigh(W)[1][:,-1] @ v)
14    boundary=[]
15    for r in ratios:
16        x=coeff*r**np.arange(T)
17        slope=float(np.polyfit(np.arange(40,T),np.log(x[40:]),1)[0])
18        boundary.append({'ratio':float(r),'slope':slope,'predicted_slope':math.log(r),
19          'xi_observed':(-1/slope if slope<0 else None),
20          'xi_predicted':(-1/math.log(r) if r<1 else None),
21          'norm_ratio':float(x[-1]/x[0])})
22    scaling=[]
23    for sigma in [.35,.55,.8,1.1]:
24        B=rng.normal(size=(n,n)); B=(B+B.T)/2*np.sqrt(2)*sigma/np.sqrt(n)
25        L=float(np.linalg.eigvalsh(B)[-1])
26        scaling.append({'sigma':sigma,'lambda_edge':L,'predicted_edge':2*sigma,
27          'critical_gain':1/L,'predicted_gain':1/(2*sigma)})
28    target=.90; alpha=.7; g=2.0/lam; hist=[]
29    for _ in range(35):
30        hist.append(g*lam); g*=math.exp(alpha*(target-g*lam))
31    return {'lambda_edge':lam,'boundary':boundary,'scaling':scaling,
32      'controller':{'target':target,'initial_edge':hist[0],'final_edge':hist[-1],'history':hist}}
33
34def nonlinear_comparison():
35    rng=np.random.default_rng(SEED+9); n=48; T=80
36    A=rng.normal(size=(n,n)); W=(A+A.T)/2
37    raw_edge=float(np.linalg.eigvalsh(W)[-1]); W=W/raw_edge
38    xs=rng.normal(0,.08,size=(T,n))
39    def rollout(g, online=False):
40        h=np.zeros(n); states=[]; gains=[]; edges=[]
41        ctl=SpectralEdgeController(target=.90,alpha=.30,g_min=.01,g_max=2.)
42        ctl.g=g
43        for x in xs:
44            z=g*(W@h+x); h=np.tanh(z); d=1-h*h
45            edge=ctl.edge_estimate(W,d,iterations=5,rng=rng)
46            if online: g=ctl.update(edge)
47            gains.append(g); edges.append(g*edge); states.append(h.copy())
48        # exact product for the realized trajectory with final fixed gain is diagnostic;
49        # use per-step gains for a more faithful sensitivity product.
50        P=np.eye(n)
51        for h,gg in zip(states,gains): P=np.diag(1-h*h)@(gg*W)@P
52        grad=float(np.linalg.svd(P,compute_uv=False)[0])
53        return {'initial_gain':gains[0],'final_gain':gains[-1],
54          'mean_effective_edge':float(np.mean(edges)),'last_effective_edge':float(edges[-1]),
55          'gradient_norm':grad,'final_state_norm':float(np.linalg.norm(states[-1]))}
56    # Same W and inputs; baseline is standard unit recurrent gain, idea starts at same gain.
57    baseline=rollout(1.0,False); idea=rollout(1.0,True)
58    return {'raw_lambda_edge':raw_edge,'normalized_lambda_edge':1.,
59      'baseline':baseline,'idea':idea,'target_edge':.90}
60
61def main():
62    t=time.time(); out={'seed':SEED,'toy':toy_sweeps(),
63      'nonlinear':nonlinear_comparison()}; out['runtime_sec']=time.time()-t
64    with open('results.json','w') as f: json.dump(out,f,indent=2)
65    print(json.dumps(out,indent=2))
66if __name__=='__main__': main()