1"""Delay-Budget Controller MVP and numerical check of its stated DDE math."""
 2import json, math
 3import numpy as np
 4from scipy.optimize import brentq, root
 5
 6
 7def omega_from_G(G, rates):
 8    rates = np.asarray(rates, dtype=float)
 9    dc = float(np.prod(rates))
10    if G <= dc:
11        return None
12    f = lambda w: float(np.prod(rates * rates + w*w)) - G*G
13    hi = max(1., math.sqrt(G) + float(np.max(rates)))
14    while f(hi) < 0: hi *= 2
15    return brentq(f, 0., hi)
16
17
18def tau_crit(G, rates, q=0):
19    w = omega_from_G(G, rates)
20    if w is None: return math.inf
21    phase = sum(math.atan(w/r) for r in rates)
22    return ((2*q+1)*math.pi-phase)/w
23
24
25def characteristic(z, G, rates, tau):
26    return np.prod(z + np.asarray(rates)) + G*np.exp(-z*tau)
27
28
29def max_root_real_part(G, rates, tau):
30    """Find the dominant roots by seeded complex Newton solves."""
31    rates=np.asarray(rates,float); scale=max(1.,float(np.max(rates)),G)
32    found=[]
33    for re in np.linspace(-2*scale, 1.0*scale, 13):
34      for im in np.linspace(-8*scale, 8*scale, 25):
35        def fun(a):
36          z=a[0]+1j*a[1]; h=characteristic(z,G,rates,tau)
37          dz=np.prod(z+rates)*sum(1/(z+r) for r in rates)-G*tau*np.exp(-z*tau)
38          # complex Newton represented as real residual
39          return [h.real,h.imag]
40        sol=root(fun,[re,im])
41        if sol.success and np.linalg.norm(fun(sol.x))<1e-6:
42          z=sol.x[0]+1j*sol.x[1]
43          if not any(abs(z-w)<1e-4 for w in found): found.append(z)
44    return float(max(z.real for z in found)) if found else float('nan')
45
46
47def controller(G, rates, tau, threshold=.8):
48    rates=np.asarray(rates,float); margin=float(np.prod(rates)); tc=tau_crit(G,rates)
49    # This is the literal actionable logic, with an explicit undefined case.
50    if G < margin:
51        return {'branch':'no-Hopf-from-formula','tau':tau,'rates':rates.tolist(),'tau_crit':tc}
52    if tau > threshold*tc:
53        return {'branch':'delay','tau':threshold*tc,'rates':rates.tolist(),'tau_crit':tc}
54    new=rates.copy(); new[1]*=2.5; new[3]*=2.5
55    return {'branch':'filter','tau':tau,'rates':new.tolist(),'tau_crit':tc}
56
57
58def simulate(G,rates,tau,T=80.,dt=.002,seed=3):
59    rng=np.random.default_rng(seed); rates=np.asarray(rates,float)
60    muX,kX,muY,kY=rates; n=int(T/dt)+1; z=np.zeros((n,4)); z[0]=rng.normal(0,.02,4)
61    lag=max(0,int(round(tau/dt))); coupling=G/(kX*muY*kY)
62    for i in range(n-1):
63      j=max(0,i-lag); x,p,y,v=z[i]
64      d=np.array([-muX*x-coupling*z[j,3], kX*(x-p), muY*(p-y), kY*(y-v)])
65      z[i+1]=z[i]+dt*d
66      if not np.all(np.isfinite(z[i+1])) or np.max(abs(z[i+1]))>1e8:
67        return {'tail_rms':float('inf'),'max':float('inf'),'diverged':True}
68    return {'tail_rms':float(np.sqrt(np.mean(z[int(.7*n):]**2))), 'max':float(np.max(abs(z))), 'diverged':False}
69
70
71def main():
72    rates=np.ones(4); G=1.2; tc=tau_crit(G,rates)
73    mathcheck={'omega':omega_from_G(G,rates),'tau_crit':tc,
74      'magnitude_relative_error':abs(math.sqrt(np.prod(rates*rates+omega_from_G(G,rates)**2))-G)/G,
75      'phase_error':abs(omega_from_G(G,rates)*tc+sum(math.atan(omega_from_G(G,rates)/r) for r in rates)-math.pi),
76      'margin':float(np.prod(rates))}
77    cases=[]
78    for tau in [0., .8*tc, 1.2*tc, 2.*tc]:
79      c=controller(G,rates,tau)
80      cases.append({'G':G,'tau':tau,'controller':c,'dominant_real_root':max_root_real_part(G,rates,tau),
81                    'simulation':simulate(G,rates,tau),
82                    'controlled_simulation':simulate(G,np.array(c['rates']),c['tau'])})
83    # Below-margin case: no positive-frequency crossing predicted, even at large delay.
84    Glo=.8; taulo=20.; cases.append({'G':Glo,'tau':taulo,'controller':controller(Glo,rates,taulo),
85      'dominant_real_root':max_root_real_part(Glo,rates,taulo),'simulation':simulate(Glo,rates,taulo)})
86    out={'math_check':mathcheck,'cases':cases}
87    with open('results.json','w') as f: json.dump(out,f,indent=2)
88    print(json.dumps(out,indent=2))
89
90if __name__=='__main__': main()