1import json
 2import numpy as np
 3from residual_screened_koopman import eig_residuals, screen_mask, spectral_forecast
 4
 5
 6def trajectory(lams, x0, T):
 7    z=np.zeros((T,len(lams))); z[0]=x0
 8    for t in range(T-1): z[t+1]=lams*z[t]
 9    return z
10
11
12def run():
13    # Prediction 1: exact eigenpairs have zero residual without corruption.
14    lams=np.array([0.8,0.98]); clean=trajectory(lams,[1.,1.],500)
15    r0=eig_residuals(np.diag(lams),clean)[2]
16
17    # Prediction 2: iid observation noise makes residual energy scale as sigma^2.
18    sigmas=np.array([0.,.01,.02,.04,.08]); obs=[]
19    for j,s in enumerate(sigmas):
20        a=[]
21        for q in range(40):
22            g=np.random.default_rng(1000+j*100+q)
23            y=clean+s*g.normal(size=clean.shape)
24            a.append(np.mean(eig_residuals(np.diag(lams),y)[2]**2))
25        obs.append(float(np.mean(a)))
26    slope=float(np.polyfit(sigmas[1:4]**2,np.array(obs[1:4]),1)[0])
27    # Theoretical small-noise slope, accounting for finite clean denominator.
28    preds=[]
29    n=len(clean)-1
30    for s in sigmas:
31        p=[]
32        for i,lam in enumerate(lams):
33            den=np.sum(clean[:-1,i]**2)
34            p.append(s*s*(1+lam*lam)*n/(den+s*s*n))
35        preds.append(float(np.mean(p)))
36
37    # Prediction 3: sparse corruption produces residual energy approximately
38    # proportional to corruption rate at fixed impulse amplitude.
39    rates=np.array([0.,.002,.005,.01,.02,.04]); impulse=.25; rate_obs=[]
40    for j,p in enumerate(rates):
41        a=[]
42        for q in range(40):
43            g=np.random.default_rng(5000+j*100+q)
44            hit=(g.random(clean.shape)<p)*g.choice([-1.,1.],size=clean.shape)
45            y=clean+impulse*hit
46            a.append(np.mean(eig_residuals(np.diag(lams),y)[2]**2))
47        rate_obs.append(float(np.mean(a)))
48    rate_slope=float(np.polyfit(rates[1:4],rate_obs[1:4],1)[0])
49    rate_preds=[]
50    n=len(clean)-1
51    for p in rates:
52        per=[]
53        for i,lam in enumerate(lams):
54            den=np.sum(clean[:-1,i]**2)
55            per.append(p*impulse**2*(1+lam*lam)*n/(den+p*impulse**2*n))
56        rate_preds.append(float(np.mean(per)))
57
58    # Held-out screening: mode 0 is reliable, mode 1 is corrupted in validation.
59    train=clean[:220].copy(); val=clean[220:].copy()
60    g=np.random.default_rng(77)
61    train[:,1]+=0.10*g.normal(size=len(train)); val[:,1]+=0.10*g.normal(size=len(val))
62    K=np.diag(lams)
63    rv=eig_residuals(K,val)[2]; mask,tau=screen_mask(rv,.5)
64    z0=val[0]; truth=clean[221:271]
65    full=spectral_forecast(K,z0,50)
66    screened=spectral_forecast(K,z0,50,mask)
67    mse_full=float(np.mean((full-truth)**2)); mse_screen=float(np.mean((screened-truth)**2))
68    result={'zero_residual':r0.tolist(),'variance':{'sigma':sigmas.tolist(),'observed':obs,'predicted':preds,'small_noise_slope':slope},'impulse_rate':{'rate':rates.tolist(),'observed':rate_obs,'predicted':rate_preds,'small_rate_slope':rate_slope},'screening':{'validation_residuals':rv.tolist(),'tau':float(tau),'mask':mask.tolist(),'retained':int(mask.sum()),'full_rollout_mse':mse_full,'screened_rollout_mse':mse_screen}}
69    print(json.dumps(result,indent=2))
70    with open('results.json','w') as f: json.dump(result,f,indent=2)
71
72if __name__=='__main__': run()