1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4from scipy.special import lambertw
  5from delay_controller import critical_delay, max_admissible_gain
  6
  7# Delay-aware plug-and-play residual capacity: analytic checks and a tiny residual toy.
  8# The DDE is integrated by explicit Euler with linear interpolation-free integer delay.
  9def tau_critical(a, G):
 10    if G <= a: return float('inf')
 11    return math.acos(-a/G) / math.sqrt(G*G-a*a)
 12
 13def gmax_for_tau(a, tau):
 14    if tau <= 0: return float('inf')
 15    # tau_c(G) decreases from infinity at G=a to zero at infinity.
 16    lo, hi = a*(1+1e-10), max(2*a, 1.0/tau + a)
 17    while tau_critical(a, hi) > tau: hi *= 2
 18    for _ in range(100):
 19        mid=(lo+hi)/2
 20        if tau_critical(a, mid) > tau: lo=mid
 21        else: hi=mid
 22    return (lo+hi)/2
 23
 24def simulate_dde(a,G,tau,dt=0.0005,T=20.0):
 25    # history e(t)=1 for t<=0; method of steps, integer delay approximation.
 26    n=int(T/dt)+1; d=max(1,int(round(tau/dt)))
 27    x=np.ones(n); growth=[]
 28    for i in range(1,n):
 29        delayed = 1.0 if i-d < 0 else x[i-d]
 30        x[i]=x[i-1] + dt*(-a*x[i-1]-G*delayed)
 31        if i > int(.5*n): growth.append(abs(x[i]))
 32    # ratio of late RMS to early RMS; unstable cases are strongly larger.
 33    early=np.sqrt(np.mean(x[int(.15*n):int(.30*n)]**2))+1e-12
 34    late=np.sqrt(np.mean(x[int(.70*n):int(.95*n)]**2))
 35    return float(late/early), float(np.max(np.abs(x)))
 36
 37def _char(z, a, G, tau):
 38    lam=complex(float(z[0]), float(z[1]))
 39    q=lam+a+G*np.exp(-lam*tau)
 40    return [q.real,q.imag]
 41
 42def spectral_abscissa(a,G,tau):
 43    if tau <= 1e-12: return -(a+G)
 44    from scipy.optimize import root
 45    vals=[]
 46    for re in np.linspace(-3,2,7):
 47        for im in np.linspace(-10,10,25):
 48            sol=root(_char,[re,im],args=(a,G,tau))
 49            if sol.success and np.linalg.norm(_char(sol.x,a,G,tau)) < 1e-7:
 50                z=complex(*sol.x)
 51                if all(abs(z-w)>1e-5 for w in vals): vals.append(z)
 52    return max(z.real for z in vals)
 53
 54def classify(a,G,tau):
 55    ab=spectral_abscissa(a,G,tau)
 56    return ab > 1e-7, ab
 57
 58def main():
 59    np.random.seed(7); random.seed(7)
 60    a=1.0
 61    # Prediction 1: below a, no finite delay destabilizes the linear model.
 62    below=[]
 63    for G in [.25,.75,1.0]:
 64        for tau in [.2,1.,3.,8.]:
 65            unstable,r=classify(a,G,tau); below.append((G,tau,unstable,r))
 66    # Prediction 2: for G>a, transition is at the closed-form tau_c.
 67    rows=[]
 68    for G in [1.25,1.5,2.,3.]:
 69        pred=tau_critical(a,G)
 70        # Numerically locate the transition using the characteristic roots.
 71        lo,hi=pred*.5,pred*1.5
 72        for _ in range(24):
 73            mid=(lo+hi)/2
 74            if classify(a,G,mid)[0]: hi=mid
 75            else: lo=mid
 76        observed=(lo+hi)/2
 77        rows.append({'G':G,'predicted_tau_c':pred,'observed_first_unstable_grid':observed,
 78                     'relative_error':abs(observed-pred)/pred if math.isfinite(observed) else None,
 79                     'root_at_predicted':spectral_abscissa(a,G,pred)})
 80    # Prediction 3: identical modules have Nmax=floor(Gmax/g), checked at fixed latency.
 81    tau=.8; g=.35; gm=gmax_for_tau(a,tau)
 82    n_pred=math.floor(gm/g) # strict inequality means this can be conservative at exact integer.
 83    count_rows=[]
 84    for N in range(1, n_pred+3):
 85        G=N*g; unstable,r=classify(a,G,tau)
 86        count_rows.append({'N':N,'G':G,'predicted_stable':G<gm,'sim_unstable':unstable,'ratio':r})
 87    # Tiny residual-stack proxy: delayed feedback recurrence, controller bypasses modules at boundary.
 88    # Each module contributes g to aggregate gain; benefit is a saturating function of active count.
 89    def stack(active, controller):
 90        chosen=0
 91        for _ in range(12):
 92            if controller and (chosen+1)*g >= gm: break
 93            chosen += 1
 94        G=chosen*g
 95        # proxy task error decreases with useful modules, but instability causes explosive error.
 96        _,maxabs=simulate_dde(a,G,tau,T=8.)
 97        error=1/(1+chosen) + 0.02*maxabs
 98        return chosen,G,error
 99    baseline=stack(12,False); controlled=stack(12,True)
100    out={'parameters':{'a':a,'dt':.0005,'T':20.0},
101         'prediction_1_G_le_a_all_delays':below,
102         'prediction_2_boundary_sweep':rows,
103         'prediction_3_count_boundary':{'tau':tau,'g':g,'Gmax':gm,'predicted_Nmax':n_pred,'rows':count_rows},
104         'toy_residual_proxy':{'unconstrained':baseline,'delay_controller':controlled}}
105    Path('results.json').write_text(json.dumps(out,indent=2))
106    print(json.dumps(out,indent=2))
107
108if __name__=='__main__': main()