1import json
 2import numpy as np
 3SEED=2058; T=80; A=np.diag([.9,.5])
 4def roll(u,b):
 5    x=np.zeros((len(u)+1,2)); B=np.array([1.,b])
 6    for t,v in enumerate(u): x[t+1]=A@x[t]+B*v
 7    return x
 8def features(x):
 9    z=x[:-1]; P=np.zeros((2*len(z),4)); P[0::2,:2]=z; P[1::2,2:]=z; return P
10def gram(u,b):
11    P=features(roll(u,b)); return P.T@P/len(u)
12def mineig(u,b): return float(np.linalg.eigvalsh(gram(u,b)+1e-12*np.eye(4))[0])
13def stacked(us,b):
14    P=np.concatenate([features(roll(u,b)) for u in us]); return float(np.linalg.eigvalsh(P.T@P/len(P)+1e-12*np.eye(4))[0])
15def fit(us,b,sig,seed):
16    r=np.random.default_rng(seed); X=[];Y=[]
17    for u in us:
18        x=roll(u,b); X.append(x[:-1]); Y.append(x[1:]+sig*r.normal(size=(len(u),2)))
19    X=np.concatenate(X);Y=np.concatenate(Y); ah=np.linalg.solve(X.T@X+1e-8*np.eye(2),X.T@Y).T
20    return float(np.linalg.norm(ah-A))
21def pulse(r, active=True):
22    if not active:
23        z=r.normal(size=T); z=np.convolve(z,np.ones(9)/9, mode='same')
24    else:
25        z=np.zeros(T); ix=r.choice(T,12,replace=False); z[ix]=r.choice([-1.,1.],12)
26    return z / max(np.linalg.norm(z),1e-12) * np.sqrt(T)
27def main():
28    r=np.random.default_rng(SEED); out={}
29    u=pulse(r); amps=np.array([.25,.5,1.,2.]); vals=np.array([mineig(a*u,.35) for a in amps]); out['scaling']={'amps':amps.tolist(),'observed_norm':(vals/vals[2]).tolist(),'predicted':((amps/amps[2])**2).tolist(),'max_abs_error':float(np.max(np.abs(vals/vals[2]-(amps/amps[2])**2)))}
30    bs=np.array([0,.02,.05,.1,.2,.35,.6,1.]); vv=np.array([mineig(u,b) for b in bs]); out['weak_channel']={'b':bs.tolist(),'observed':vv.tolist(),'predicted_zero_at_b0':float(vv[0]),'monotonic':bool(np.all(np.diff(vv)>=-1e-12))}
31    out['weak_channel']['normalized_b2']=(vv[1:]/vv[-1]).tolist(); out['weak_channel']['predicted_b2']=((bs[1:]/bs[-1])**2).tolist()
32    candidates=[pulse(r) for _ in range(500)]; chosen=[]; remaining=list(candidates)
33    for _ in range(6):
34        j=max(range(len(remaining)),key=lambda i:stacked(chosen+[remaining[i]],.35)); q=remaining.pop(j); chosen.append(q)
35    passive=[pulse(r,False) for _ in range(6)]; active=chosen
36    pg=stacked(passive,.35); ag=stacked(active,.35)
37    out['design']={'passive_gramian':pg,'active_gramian':ag,'gain':ag/max(pg,1e-12),'candidate_pool':500,'greedy_steps':6}
38    for sig in [.02,.1]: out['identification_'+str(sig)]={'passive':float(np.mean([fit(passive,.35,sig,SEED+i) for i in range(20)])),'active':float(np.mean([fit(active,.35,sig,SEED+i) for i in range(20)]))}
39    print(json.dumps(out,indent=2))
40if __name__=='__main__': main()