1import json, math
 2import numpy as np
 3from scipy.linalg import eigvalsh, solve_discrete_lyapunov
 4
 5SEED = 2075
 6
 7def hessians():
 8    a=.42; R=np.array([[math.cos(a),-math.sin(a)],[math.sin(a),math.cos(a)]])
 9    return [np.diag([1.,8.]), R@np.diag([1.5,5.5])@R.T]
10
11def lmi_matrices(Hs,K,alpha,sigma,lam=1.):
12    n=K.shape[0]; Abar=-K@sum(Hs)/len(Hs)
13    P=solve_discrete_lyapunov(np.eye(n)+.02*Abar,np.eye(n)); P=P/np.trace(P)*n
14    out=[]
15    for H in Hs:
16        A=-K@H; B=A
17        M=np.block([[A.T@P+P@A+2*alpha*P+lam*sigma*sigma*np.eye(n),P@B],
18                    [B.T@P,-lam*np.eye(n)]])
19        out.append((M+M.T)/2)
20    return P,out
21
22def margin(Hs,K,alpha,sigma,lam=1.):
23    return max(float(eigvalsh(M).max()) for M in lmi_matrices(Hs,K,alpha,sigma,lam)[1])
24
25def max_certified_sigma(Hs,K,alpha=.02):
26    lo,hi=0.,1.
27    while margin(Hs,K,alpha,hi)<0 and hi<100: hi*=2
28    for _ in range(55):
29        mid=(lo+hi)/2
30        if margin(Hs,K,alpha,mid)<0: lo=mid
31        else: hi=mid
32    return lo
33
34def simulate(Hs,K,h,sigma,steps=1200,max_interval=10000):
35    x=np.array([2.,-1.5]); last=x.copy(); events=0; since=0; intervals=[]; norms=[]; losses=[]
36    for t in range(steps):
37        H=Hs[t%len(Hs)]; x=x-h*K@(H@last); since+=1
38        d=last-x
39        if sigma==0 or d@d>sigma*sigma*(x@x+1e-30) or since>=max_interval:
40            last=x.copy(); events+=1; intervals.append(since); since=0
41        norms.append(float(np.linalg.norm(x))); losses.append(float(.5*x@H@x))
42        if not np.all(np.isfinite(x)) or norms[-1]>1e12: break
43    return {'stable':len(norms)==steps and max(norms)<1e12,'events':events,'fraction':events/steps,'mean_interval':float(np.mean(intervals)) if intervals else steps,'final_norm':norms[-1],'tail_loss':float(np.mean(losses[-max(20,steps//10):])),'max_norm':max(norms)}
44
45def main():
46    Hs=hessians(); Hb=sum(Hs)/2
47    Kd=np.diag(.72/np.diag(Hb)); Kf=.72*np.linalg.inv(Hb)
48    print('HESSIANS',json.dumps([H.tolist() for H in Hs]))
49    print('LMI_AND_TRIGGER')
50    for name,K in [('diagonal',Kd),('dense',Kf)]:
51        sig=max_certified_sigma(Hs,K)
52        vals=[]
53        for s in [0,.5*sig,.9*sig,1.1*sig]: vals.append({'sigma':s,'margin':margin(Hs,K,.02,s),'run':simulate(Hs,K,.16,s)})
54        print(json.dumps({'method':name,'certified_sigma':sig,'sweep':vals}))
55    print('DISCRETE_BOUNDARY')
56    for name,K in [('diagonal',Kd),('dense',Kf)]:
57        rho=max(max(abs(np.linalg.eigvals(K@H))) for H in Hs); boundary=2/rho
58        rows=[]
59        for mult in [.9,1.02]:
60            r=simulate(Hs,K,mult*boundary,0,steps=200,max_interval=1)
61            rows.append({'mult':mult,'h':mult*boundary,'stable':r['stable'],'max_norm':r['max_norm']})
62        print(json.dumps({'method':name,'rho':float(rho),'boundary':boundary,'rows':rows}))
63    print('LMI_RANDOM_IMPLICATION')
64    rng=np.random.default_rng(SEED); K=Kf; sig=.8*max_certified_sigma(Hs,K); P,Ms=lmi_matrices(Hs,K,.02,sig)
65    worst=-1e9
66    for M in Ms:
67        for _ in range(20000):
68            e=rng.normal(size=2); d=rng.normal(size=2); d*=sig*np.linalg.norm(e)/(np.linalg.norm(d)+1e-30)
69            z=np.r_[e,d]; worst=max(worst,float(z@M@z))
70    print(json.dumps({'sigma':sig,'max_sampled_quadratic_form':worst,'all_negative':worst<0}))
71
72if __name__=='__main__': main()