import json, math
import numpy as np
from scipy.linalg import solve_discrete_lyapunov, eigvalsh
from event_triggered_preconditioner import make_hessians, lmi_margin


def fixed_run(H, K, h, sigma, steps=80, max_interval=10000):
    x=np.array([2.0,-1.5]); last=x.copy(); events=0; since=0; norms=[]
    for t in range(steps):
        x=x-h*K@(H@last); since+=1
        d=last-x
        if d@d > sigma*sigma*(x@x+1e-12) or since>=max_interval:
            last=x.copy(); events+=1; since=0
        norms.append(np.linalg.norm(x))
        if not np.isfinite(norms[-1]) or norms[-1]>1e10: return False, events, norms
    return True,events,norms


def main():
    Hs=make_hessians(); H=Hs[0]
    Hbar=sum(Hs)/2
    Ks={'diag':np.diag(.72/np.diag(Hbar)), 'dense':.72*np.linalg.inv(Hbar)}
    print('PREDICTION_1 discrete fixed-H boundary: h_boundary=2/rho(KH)')
    for name,K in Ks.items():
        rho=max(abs(np.linalg.eigvals(K@H)))
        hb=2/rho
        rows=[]
        for mult in [.90,.99,1.01,1.10]:
            ok,_,norms=fixed_run(H,K,mult*hb,0,steps=100,max_interval=1)
            rows.append({'mult':mult,'h':mult*hb,'stable_100_steps':ok,'final_norm':float(norms[-1])})
        print(json.dumps({'method':name,'rho':float(rho),'predicted_boundary':float(hb),'observed':rows}))

    print('PREDICTION_2 LMI threshold: margin crosses zero as sigma increases')
    for name,K in Ks.items():
        lo,hi=0.,10.
        # Fixed P and lambda as in the implementation; locate the numerical certificate boundary.
        while lmi_margin(Hs,K,.02,hi,lam=1)[0]<0: hi*=2
        for _ in range(60):
            mid=(lo+hi)/2
            if lmi_margin(Hs,K,.02,mid,lam=1)[0]<0: lo=mid
            else: hi=mid
        vals=[]
        for s in [0,.25*lo,.5*lo,.75*lo,lo,1.05*lo]:
            mar,_=lmi_margin(Hs,K,.02,s,lam=1)
            vals.append({'sigma':float(s),'margin':float(mar),'certified':bool(mar<0)})
        print(json.dumps({'method':name,'predicted_sigma_boundary':float(lo),'observed':vals}))

    print('PREDICTION_3 finite-horizon trigger scaling: larger sigma lowers events, with bounded error')
    H=Hs[0]; K=Ks['dense']; h=.25
    # This h is safely below the fixed-H synchronous boundary and exposes stale updates.
    for sigma in [0,.02,.05,.10,.20,.40]:
        ok,ev,norms=fixed_run(H,K,h,sigma,steps=40,max_interval=1000)
        sync=fixed_run(H,K,h,0,steps=40,max_interval=1)[2][-1]
        print(json.dumps({'sigma':sigma,'events':ev,'event_fraction':ev/40,
            'norm_at_40':float(norms[-1]),'ratio_to_sync_norm':float(norms[-1]/(sync+1e-30)),
            'stable':ok}))

if __name__=='__main__': main()