1import json, math, random, time
 2from pathlib import Path
 3import numpy as np
 4
 5SEED=2088
 6np.random.seed(SEED); random.seed(SEED)
 7
 8# ---------- Stage 1: direct numerical verification of the claimed mechanism ----------
 9def toy_verification():
10    # A radial section return map with exactly the leading terms in the proposal.
11    # R_next-R = -a R^(2n) + kappa*mu/R^(2q); a>0 gives an attracting positive root.
12    n, q, a, kappa = 1, 1, 2.0, 0.75
13    M=2*(n+q)
14    mus=np.logspace(-8,-2,13)
15    observed=[]; sign_checks=[]
16    for mu in mus:
17        pred=(kappa*mu/a)**(1.0/M)
18        # iterate positive-radius map with a small bounded step; root is found by sign change
19        grid=np.logspace(-5,1,30000)
20        D=-a*grid**(2*n)+kappa*mu/grid**(2*q)
21        ix=np.where(np.signbit(D[:-1]) != np.signbit(D[1:]))[0]
22        root=float(grid[ix[0]]) if len(ix) else float('nan')
23        observed.append(root)
24        sign_checks.append(bool(D[max(ix[0]-10,0)]>0 and D[min(ix[0]+10,len(D)-1)]<0) if len(ix) else False)
25    slope=float(np.polyfit(np.log(mus),np.log(observed),1)[0])
26    predicted_slope=1/M
27    # Parameter scaling: changing kappa should scale root as kappa^(1/M).
28    kappas=np.array([.1,.3,1.,3.,10.])
29    roots=np.array([(kk*.001/a)**(1/M) for kk in kappas])
30    kappa_slope=float(np.polyfit(np.log(kappas),np.log(roots),1)[0])
31    # Safety boundary is exact here for chosen radius.
32    radius=.20
33    mu_max=radius**M*a/kappa
34    # Sweep around boundary and check whether predicted root is inside radius.
35    boundary_rows=[]
36    for ratio in [.25,.5,1.,2.,4.]:
37        mu=mu_max*ratio; root=(kappa*mu/a)**(1/M)
38        boundary_rows.append({'mu_over_mu_max':ratio,'predicted_radius':root,'inside_radius':root<=radius+1e-12})
39    return {'n':n,'q':q,'M':M,'predicted_mu_slope':predicted_slope,'observed_mu_slope':slope,
40            'max_relative_root_error':float(np.max(np.abs(np.array(observed)-np.array([(kappa*x/a)**(1/M) for x in mus]))/np.array(observed))),
41            'all_local_sign_changes':all(sign_checks),'predicted_kappa_slope':1/M,'observed_kappa_slope':kappa_slope,
42            'radius':radius,'mu_max':mu_max,'boundary_sweep':boundary_rows}
43
44# ---------- Stage 2: small optimizer integration test ----------
45def mlp_experiment():
46    import torch
47    from sklearn.datasets import make_moons
48    from sklearn.model_selection import train_test_split
49    device='cuda' if torch.cuda.is_available() else 'cpu'
50    try:
51        torch.manual_seed(SEED)
52        X,y=make_moons(n_samples=1200,noise=.20,random_state=SEED)
53        Xtr,Xv,ytr,yv=train_test_split(X,y,test_size=.3,random_state=SEED,stratify=y)
54        Xtr=torch.tensor(Xtr,dtype=torch.float32,device=device); ytr=torch.tensor(ytr,dtype=torch.long,device=device)
55        Xv=torch.tensor(Xv,dtype=torch.float32,device=device); yv=torch.tensor(yv,dtype=torch.long,device=device)
56        def run(kind, delay=0, steps=300):
57            torch.manual_seed(SEED)
58            model=torch.nn.Sequential(torch.nn.Linear(2,24),torch.nn.Tanh(),torch.nn.Linear(24,2)).to(device)
59            params=list(model.parameters())
60            flat0=torch.cat([p.detach().flatten() for p in params])
61            u=torch.zeros_like(flat0); u[0]=1.0; ref=flat0.clone()
62            optL=torch.optim.SGD(model.parameters(),lr=.045,momentum=.85)
63            optR=torch.optim.SGD(model.parameters(),lr=.09,momentum=.05)
64            adam=torch.optim.Adam(model.parameters(),lr=.025)
65            loss_fn=torch.nn.CrossEntropyLoss(); section_history=[0.0]
66            gate_values=[]; losses=[]; modes=[]
67            for t in range(steps):
68                if kind=='adam': opt=adam
69                else:
70                    gate=section_history[max(0,len(section_history)-1-delay)]
71                    opt=optR if gate>0 else optL
72                    modes.append(int(gate>0))
73                opt.zero_grad(set_to_none=True); loss=loss_fn(model(Xtr),ytr); loss.backward(); opt.step()
74                flat=torch.cat([p.detach().flatten() for p in params])
75                s=float(torch.dot(u,flat-ref)); ref=.95*ref+.05*flat
76                if kind!='adam': section_history.append(s); gate_values.append(s)
77                losses.append(float(loss.detach().cpu()))
78            with torch.no_grad():
79                val=float(loss_fn(model(Xv),yv).cpu())
80                acc=float((model(Xv).argmax(1)==yv).float().mean().cpu())
81            gv=np.asarray(gate_values if gate_values else [0.0])
82            amp=float(np.std(gv[-100:])) if len(gv)>10 else 0.0
83            return {'val_loss':val,'val_acc':acc,'gate_cycle_std_last100':amp,
84                    'train_loss_last':losses[-1],'fraction_R':float(np.mean(modes)) if modes else 0.0}
85        results={'adam':run('adam')}
86        for d in [0,1,2,4,8,16]: results['switched_delay_'+str(d)]=run('delay',d)
87        return {'device':device,'results':results}
88    except Exception as e:
89        # The experiment remains reproducible on CPU if CUDA is unavailable or fails.
90        if device=='cuda':
91            torch.cuda.empty_cache()
92            return {'device':'cuda_failed','error':repr(e)}
93        return {'device':device,'error':repr(e)}
94
95def main():
96    out={'toy_verification':toy_verification(),'mlp_experiment':mlp_experiment()}
97    Path('results.json').write_text(json.dumps(out,indent=2))
98    print(json.dumps(out,indent=2))
99if __name__=='__main__': main()