1import json, math, random
 2from pathlib import Path
 3import numpy as np
 4
 5SEED = 2096
 6np.random.seed(SEED); random.seed(SEED)
 7OUT = Path('results.json')
 8
 9def toy_scaling():
10    # g(t)=a+b*t+c*t^2 gives exact linear-extrapolation residual 2|c|h^2.
11    hs = np.logspace(-3, -0.3, 12)
12    a, b, c = 0.7, -0.4, 1.3
13    vals = []
14    for h in hs:
15        old2 = a + b*(-h) + c*(-h)**2
16        old1 = a
17        current = a + b*h + c*h*h
18        predicted = old1 + (old1-old2)
19        vals.append(abs(current-predicted))
20    slope = float(np.polyfit(np.log(hs), np.log(vals), 1)[0])
21    C, eps = 2.0*abs(c), 0.05
22    claimed_target = math.sqrt(eps/C)
23    # This is exactly the formula in the idea: h_next=sqrt(eps/||r||),
24    # with projection and ratio limiting. Its fixed point is fourth-root,
25    # because ||r||=C*h^2; this is tested separately from the claimed target.
26    h, hmin, hmax, q, delta = 0.8, 1e-4, 2.0, 1.35, 1e-12
27    trace=[]
28    for _ in range(100):
29        residual=C*h*h
30        raw=math.sqrt(eps/(residual+delta))
31        nh=float(np.clip(raw,hmin,hmax)); nh=float(np.clip(nh,h/q,q*h))
32        trace.append((h,residual,nh)); h=nh
33    observed_h=float(np.mean([z[0] for z in trace[-10:]]))
34    observed_r=float(np.mean([z[1] for z in trace[-10:]]))
35    formula_fixed_point=(eps/C)**0.25
36    rng=np.random.default_rng(SEED); h=.2; noisy_ratios=[]
37    for _ in range(100):
38        measured=max(1e-9,C*h*h*(1+.15*rng.normal()))
39        nh=float(np.clip(math.sqrt(eps/measured),hmin,hmax))
40        nh=float(np.clip(nh,h/q,q*h)); noisy_ratios.append(nh/h); h=nh
41    allrat=np.array([z[2]/z[0] for z in trace]+noisy_ratios)
42    return {'residual_scaling': {'observed_slope':slope,'predicted_slope':2.0,'absolute_error':abs(slope-2)},
43      'scheduler_target': {'C':C,'epsilon':eps,'claimed_sqrt_target':claimed_target,
44        'observed_h':observed_h,'observed_residual':observed_r,
45        'formula_fixed_point_fourth_root':formula_fixed_point,
46        'observed_vs_fourth_root_relative_error':abs(observed_h-formula_fixed_point)/formula_fixed_point,
47        'claimed_target_relative_error':abs(observed_h-claimed_target)/claimed_target,
48        'note':'The implemented formula lacks the multiplicative h_k present in the paper-style ratio controller.'},
49      'ratio_sweep': {'observed_min':float(allrat.min()),'observed_max':float(allrat.max()),'predicted_bounds':[1/q,q]}}
50
51def benchmark():
52    import torch
53    torch.manual_seed(SEED); np.random.seed(SEED)
54    device='cuda' if torch.cuda.is_available() else 'cpu'
55    try:
56        n=1200; rng=np.random.default_rng(SEED)
57        x=np.concatenate([rng.normal([-1.,-1.],.75,(n//2,2)),rng.normal([1.,1.],.75,(n//2,2))]).astype('float32')
58        y=np.concatenate([np.zeros(n//2),np.ones(n//2)]).astype('int64'); p=rng.permutation(n); x=x[p]; y=y[p]
59        X=torch.tensor(x,device=device); Y=torch.tensor(y,device=device)
60        def run(kind):
61            torch.manual_seed(SEED+3)
62            model=torch.nn.Sequential(torch.nn.Linear(2,24),torch.nn.Tanh(),torch.nn.Linear(24,2)).to(device)
63            lossfn=torch.nn.CrossEntropyLoss(); params=list(model.parameters()); prev=[]; intervals=[.08,.08]
64            h=.08; eps=.003; hmin=.015; hmax=.20; q=1.25; beta=.8; R=None; losses=[]; hs=[]; accs=[]
65            for k in range(260):
66                ix=torch.randint(0,n,(64,),device=device); model.zero_grad(set_to_none=True)
67                loss=lossfn(model(X[ix]),Y[ix]); loss.backward(); g=torch.cat([p.grad.detach().reshape(-1) for p in params])
68                if kind=='adaptive' and len(prev)>=2:
69                    pred=prev[-1]+(h/intervals[-1])*(prev[-1]-prev[-2]); res=torch.linalg.vector_norm(g-pred).item()
70                    R=res if R is None else beta*R+(1-beta)*res
71                    nh=float(np.clip(math.sqrt(eps/(R+1e-8)),hmin,hmax)); nh=float(np.clip(nh,h/q,q*h))
72                else: nh=h
73                with torch.no_grad():
74                    for z in params: z -= nh*z.grad
75                prev.append(g); prev=prev[-3:]; intervals.append(nh); intervals=intervals[-3:]
76                h=nh if kind=='adaptive' else (.08*(.5*(1+math.cos(math.pi*(k+1)/260)) if kind=='cosine' else 1.))
77                losses.append(float(loss)); hs.append(nh)
78                with torch.no_grad(): accs.append(float((model(X).argmax(1)==Y).float().mean()))
79            return {'final_loss':float(np.mean(losses[-20:])),'final_accuracy':float(np.mean(accs[-20:])),'best_loss':float(min(losses)),'mean_interval':float(np.mean(hs)),'interval_min':float(min(hs)),'interval_max':float(max(hs))}
80        out={k:run(k) for k in ['fixed','cosine','adaptive']}; out['device']=device; return out
81    except Exception as e:
82        if device=='cuda':
83            torch.cuda.empty_cache(); return {'cuda_error':str(e),'fallback':'run with CUDA_VISIBLE_DEVICES=-1'}
84        return {'error':str(e)}
85
86def main():
87    data={'seed':SEED,'mechanism':toy_scaling()}
88    try: data['benchmark']=benchmark()
89    except Exception as e: data['benchmark_error']=repr(e)
90    OUT.write_text(json.dumps(data,indent=2)); print(json.dumps(data,indent=2))
91if __name__=='__main__': main()