1import math, random, json
  2from pathlib import Path
  3import numpy as np
  4from scipy.linalg import expm
  5
  6SEED=3154
  7np.random.seed(SEED); random.seed(SEED)
  8
  9def hamiltonian(a,b,q,r):
 10    return np.array([[a, -b*b/r], [-q, -a]], dtype=float)
 11
 12def feedback(a,b,q,r,pf,T):
 13    M=expm(hamiltonian(a,b,q,r)*T)
 14    D=float(M[0,0]+M[0,1]*pf)
 15    P=float((M[1,0]+M[1,1]*pf)/D) if abs(D)>1e-15 else np.nan
 16    return D,P,M
 17
 18def math_check():
 19    # Direct determinant verification of the formula in the idea.
 20    samples=[(-1.,1.,1.,.2),(-1.,1.,.1,1.),(-1.,1.,1.,1.),(.3,2.,.7,.4)]
 21    rows=[]
 22    for a,b,q,r in samples:
 23        H=hamiltonian(a,b,q,r)
 24        stated=b*b*q/r-a*a
 25        actual=float(np.linalg.det(H))
 26        rows.append({'a':a,'b':b,'q':q,'r':r,'stated_det':stated,'actual_det':actual,
 27                     'identity_error':abs(actual-(-a*a-b*b*q/r))})
 28    # For valid q>=0,r>0 dynamics, roots can still occur in the hyperbolic case.
 29    a,b,q,r,pf=-1.,1.,1.,.2,1.
 30    H=hamiltonian(a,b,q,r)
 31    disc=a*a+b*b*q/r
 32    # exp(H T)=cosh(kT)I+sinh(kT)H/k, k^2=disc.
 33    k=math.sqrt(disc)
 34    coeff=(a-pf*b*b/r)/k
 35    # D=cosh(kT)+coeff*sinh(kT); solve tanh(kT)=-1/coeff if possible.
 36    predicted=None
 37    z=-1./coeff
 38    if coeff < -1 and 0 < z < 1:
 39        predicted=math.atanh(z)/k
 40    Ts=np.linspace(0,4,4001)
 41    Ds=np.array([feedback(a,b,q,r,pf,float(t))[0] for t in Ts])
 42    crossings=np.where(Ds[:-1]*Ds[1:]<=0)[0]
 43    numerical=float(Ts[crossings[0]]) if len(crossings) else None
 44    return {'determinant_rows':rows,
 45            'all_determinant_identities_correct':all(x['identity_error']<1e-12 for x in rows),
 46            'claimed_elliptic_det_for_first_case':rows[0]['stated_det'],
 47            'actual_det_for_first_case':rows[0]['actual_det'],
 48            'hyperbolic_predicted_conjugate':predicted,
 49            'hyperbolic_numerical_crossing_grid':numerical,
 50            'conjugate_point_observed':numerical is not None,
 51            'claim_formula_consistent':all(abs(x['stated_det']-x['actual_det'])<1e-12 for x in rows)}
 52
 53def quadratic_run(kind, steps=250, lr=.08, T=.18):
 54    A=np.diag([.5,1.,3.,8.]); x=np.array([2.,-1.5,.8,-.5],dtype=float)
 55    losses=[]; horizons=[]; gains=[]; denoms=[]
 56    for _ in range(steps):
 57        grad=A@x
 58        if kind=='sgd':
 59            x=x-lr*grad; horizons.append(0.); gains.append(1.); denoms.append(1.)
 60        else:
 61            gs=[]; ds=[]; ts=[]
 62            for lam in np.diag(A):
 63                r=.2; t=T
 64                if kind=='adaptive':
 65                    while True:
 66                        d,p,_=feedback(-lam,1.,lam,r,lam,t)
 67                        if d>.05 and np.isfinite(p): break
 68                        t*=.8
 69                        if t<1e-5: r*=2.; t=T
 70                    d,p,_=feedback(-lam,1.,lam,r,lam,t)
 71                else: d,p,_=feedback(-lam,1.,lam,r,lam,T)
 72                if not np.isfinite(p) or d<=.005: p=np.sign(p)*min(abs(p) if np.isfinite(p) else 1e6,100.)
 73                gs.append(p/r); ds.append(d); ts.append(t)
 74            x=x-lr*np.array(gs)*x
 75            horizons.append(float(np.mean(ts))); gains.append(float(np.max(np.abs(gs)))); denoms.append(float(np.min(ds)))
 76        losses.append(.5*float(x@A@x))
 77    return {'final_loss':losses[-1],'min_loss':min(losses),'max_gain':max(gains),
 78            'min_D':min(denoms),'final_x_norm':float(np.linalg.norm(x)),
 79            'mean_horizon':float(np.mean(horizons)),
 80            'diverged':bool(not np.isfinite(losses[-1]) or losses[-1]>1e8)}
 81
 82def mlp_run(kind, epochs=18):
 83    try:
 84        import torch
 85        from sklearn.datasets import load_digits
 86        from sklearn.model_selection import train_test_split
 87        torch.manual_seed(SEED)
 88        dev='cuda' if torch.cuda.is_available() else 'cpu'
 89        X,y=load_digits(return_X_y=True); X=X.astype('float32')/16.; y=y.astype('int64')
 90        Xtr,Xte,ytr,yte=train_test_split(X,y,test_size=.25,random_state=SEED,stratify=y)
 91        Xtr=torch.tensor(Xtr,device=dev); ytr=torch.tensor(ytr,device=dev); Xte=torch.tensor(Xte,device=dev); yte=torch.tensor(yte,device=dev)
 92        model=torch.nn.Sequential(torch.nn.Linear(64,32),torch.nn.Tanh(),torch.nn.Linear(32,10)).to(dev)
 93        opt=torch.optim.SGD(model.parameters(),lr=.08); losses=[]
 94        for _ in range(epochs):
 95            for ix in torch.randperm(len(Xtr),device=dev).split(128):
 96                opt.zero_grad(); loss=torch.nn.functional.cross_entropy(model(Xtr[ix]),ytr[ix]); loss.backward()
 97                if kind=='hh':
 98                    d,p,_=feedback(-1.,1.,1.,.5,1.,.2); mult=min(max(p/.5,.25),4.)
 99                    for par in model.parameters():
100                        if par.grad is not None: par.grad.mul_(mult)
101                opt.step()
102            with torch.no_grad(): losses.append(float(torch.nn.functional.cross_entropy(model(Xtr),ytr).cpu()))
103        with torch.no_grad(): acc=float((model(Xte).argmax(1)==yte).float().mean().cpu())
104        return {'final_loss':losses[-1],'test_accuracy':acc,'device':dev,
105                'feedback_multiplier':(feedback(-1,1,1,.5,1,.2)[1]/.5 if kind=='hh' else 1.)}
106    except Exception as e:
107        return {'error':type(e).__name__+': '+str(e)}
108
109if __name__=='__main__':
110    result={'math_check':math_check(),'quadratic':{},'digits_mlp':{}}
111    for k in ('sgd','fixed_horizon','adaptive'): result['quadratic'][k]=quadratic_run(k)
112    for k in ('sgd','hh'): result['digits_mlp'][k]=mlp_run(k)
113    Path('results.json').write_text(json.dumps(result,indent=2)); print(json.dumps(result,indent=2))