1import json, math, random
 2import numpy as np
 3from scipy.optimize import least_squares, minimize
 4from scipy.signal import residue
 5
 6SEED=7
 7np.random.seed(SEED); random.seed(SEED)
 8
 9def tf_step(t, T, z):
10    # Exact inverse Laplace response of (1+z s)/(s prod(1+T_i s)).
11    den=np.poly([-1.0/x for x in T])
12    num=np.array([z,1.0])
13    # append integrator for the step; residue gives y=sum(r exp(p t)).
14    p,r,k=residue(num, np.polymul(den,[1.,0.]))
15    y=np.zeros_like(t,dtype=float)
16    for pi,ri in zip(p,r): y += np.real(ri*np.exp(np.real(pi)*t))
17    return y+float(np.real(k[0])) if len(k) else y
18
19def moment_rho(T,z):
20    B=float(sum(T)-z)
21    C=float(T[0]*T[1]+T[0]*T[2]+T[1]*T[2]-z*sum(T)+z*z)
22    return B*B/C
23
24def twostate_response(t,p):
25    Tf,Ts,z=np.exp(p)
26    den=np.poly([-1/Tf,-1/Ts])
27    q,r,k=residue(np.array([z,1.]),np.polymul(den,[1.,0.]))
28    y=np.zeros_like(t,dtype=float)
29    for qi,ri in zip(q,r): y += np.real(ri*np.exp(np.real(qi)*t))
30    return y+(float(np.real(k[0])) if len(k) else 0.)
31
32def fit_window(t,y,init=None):
33    if init is None: init=np.log([.3,1.5,.1])
34    r=least_squares(lambda p:twostate_response(t,p)-y,init,max_nfev=70)
35    Tf,Ts,z=np.exp(r.x); a1=1/Tf+1/Ts; a0=1/(Tf*Ts)
36    return float(a1*a1/a0),r.x,float(np.mean(r.fun*r.fun))
37
38def find_case():
39    t=np.linspace(.01,10,121); best=None
40    for _ in range(180):
41        T=np.exp(np.random.uniform(np.log(.08),np.log(4),3)); z=np.random.uniform(-1,1)
42        rm=moment_rho(T,z)
43        if not 4.35<rm<4.65: continue
44        y=tf_step(t,T,z); rw,_,err=fit_window(t,y)
45        score=abs(rm-4.5)+abs(rw-3.31)
46        if best is None or score<best[0]: best=(score,T,z,rm,rw,err)
47    if best is None: raise RuntimeError('no case; increase search')
48    return t,best[1],best[2],best[3],best[4]
49
50def reduction_penalty(T,z,lam=.35,margin=.15):
51    ts=np.linspace(.01,10,81); rm=moment_rho(T,z); rw,_,_=fit_window(ts,tf_step(ts,T,z))
52    return abs(rm-rw)+lam*max(0,margin-abs(rm-4))+lam*max(0,margin-abs(rw-4)),rm,rw
53
54def optimize_model(t,y,T0,z0,regularized):
55    x0=np.r_[np.log(T0),z0]
56    def obj(x):
57        T=np.exp(x[:3]); z=x[3]; task=np.mean((tf_step(t,T,z)-y)**2)
58        if not regularized:return task
59        return task+.002*reduction_penalty(T,z)[0]
60    r=minimize(obj,x0,method='Nelder-Mead',options={'maxiter':70,'xatol':2e-5,'fatol':1e-9})
61    T=np.exp(r.x[:3]);z=float(r.x[3]); pen,rm,rw=reduction_penalty(T,z)
62    return {'task_mse':float(np.mean((tf_step(t,T,z)-y)**2)),'rho_moment':float(rm),'rho_window':float(rw),'disagreement':float(abs(rm-rw)),'objective':float(r.fun),'T':T.tolist(),'z':z}
63
64def main():
65    t,T,z,rm,rw=find_case(); y=tf_step(t,T,z)
66    check=[]
67    for rho in [2.,3.31,4.,4.5,7.]:
68        D=rho-4;check.append({'rho':rho,'Delta':D,'class':'complex' if D<0 else 'real'})
69    out={'seed':SEED,'source_system':{'T':T.tolist(),'z':float(z)},'initial_reductions':{'rho_moment':rm,'rho_window':rw,'disagreement':abs(rm-rw),'window_fit_mse':float(np.mean((tf_step(t,T,z)-tf_step(t,T,z))**2))},'discriminant_check':check,'baseline':optimize_model(t,y,T,z,False),'idea':optimize_model(t,y,T,z,True)}
70    with open('results.json','w') as f:json.dump(out,f,indent=2)
71    print(json.dumps(out,indent=2))
72if __name__=='__main__':main()