1import json
  2import math
  3import numpy as np
  4from scipy.special import lambertw
  5from scipy.optimize import brentq
  6
  7SEED = 7
  8rng = np.random.default_rng(SEED)
  9
 10def roots(g, tau_r=1.0, tau_d=0.0, branches=range(-30,31)):
 11    """Roots of tau_r*lambda=-1+g exp(-lambda*tau_d), scalar real mode."""
 12    if tau_d == 0:
 13        return np.array([(-1.0 + g) / tau_r], dtype=complex)
 14    q = 1.0 / tau_r
 15    z = g * tau_d / tau_r * np.exp(tau_d / tau_r)
 16    return np.array([lambertw(z, k) / tau_d - q for k in branches])
 17
 18def dominant(g, tau_r=1.0, tau_d=0.0):
 19    r = roots(g, tau_r, tau_d)
 20    return r[np.argmax(r.real)]
 21
 22def discrete_growth(g, tau_r=1.0, tau_d=0.5, dt=0.002, duration=18.0):
 23    """Euler growth test initialized with the analytical dominant mode."""
 24    delay = int(round(tau_d / dt)); n = int(duration / dt)
 25    lam = dominant(g, tau_r, tau_d)
 26    times = (np.arange(n + delay + 1) - delay) * dt
 27    x = np.exp(lam * times).astype(complex)
 28    for i in range(delay, delay+n):
 29        x[i+1] = x[i] + dt/tau_r * (-x[i] + g*x[i-delay])
 30    t = np.arange(n) * dt
 31    sel = t > duration*0.25
 32    return float(np.polyfit(t[sel], np.log(np.abs(x[delay:delay+n][sel])+1e-30), 1)[0])
 33
 34def hopf_prediction(tau_r=1.0, tau_d=0.5):
 35    # Negative feedback first Hopf crossing: omega*tau_d + atan(omega*tau_r)=pi.
 36    f = lambda om: om*tau_d + math.atan(om*tau_r) - math.pi
 37    om = brentq(f, 1e-9, 30.0)
 38    return om, math.sqrt(1.0 + (om*tau_r)**2)
 39
 40def toy_verification():
 41    # Prediction 1: at zero delay the boundary is exactly g=1.
 42    gs = np.array([0.85, 0.98, 1.00, 1.02, 1.15])
 43    observed = [float(dominant(g).real) for g in gs]
 44    boundary = brentq(lambda g: dominant(g).real, .5, 1.5)
 45
 46    # Prediction 2: for finite delay, instability begins at the Hopf gain.
 47    om, gcrit = hopf_prediction()
 48    sweep = np.linspace(-gcrit-0.12, -gcrit+0.12, 13)
 49    reals = np.array([dominant(g, tau_d=.5).real for g in sweep])
 50    cross = brentq(lambda g: dominant(g, tau_d=.5).real, -gcrit-.2, -gcrit+.2)
 51
 52    # Prediction 3: dominant growth changes with gain according to the root;
 53    # validate against a direct delayed Euler rollout at three gains.
 54    gs3 = np.array([-gcrit-.04, -gcrit-.16, -gcrit-.30])
 55    predicted = np.array([dominant(g, tau_d=.5).real for g in gs3])
 56    measured = np.array([discrete_growth(g) for g in gs3])
 57    return {
 58        "zero_delay": {"predicted_boundary": 1.0, "observed_boundary": float(boundary),
 59                       "gains": gs.tolist(), "root_real_parts": observed},
 60        "delayed_hopf": {"predicted_omega": om, "predicted_gain": gcrit,
 61                          "observed_crossing": float(cross), "sweep_gains": sweep.tolist(),
 62                          "sweep_root_real": reals.tolist()},
 63        "growth_scaling": {"gains": gs3.tolist(), "predicted_real_root": predicted.tolist(),
 64                            "measured_euler_log_growth": measured.tolist(),
 65                            "mean_abs_error": float(np.mean(np.abs(predicted-measured)))}
 66    }
 67
 68def stabilized_ring_demo(N=64, tau_r=1.0, tau_d=.5, dt=.002, duration=12.0):
 69    # Build a ring with selected Fourier gains 1.35 (unstable for this delay),
 70    # and an idea version that clips those gains to 0.82. This changes no parameter count.
 71    selected = [1,2,3,N-3,N-2,N-1]
 72    wh_base = np.zeros(N, dtype=complex)
 73    for k in selected: wh_base[k] = 1.35
 74    # small harmless modes make the kernel nontrivial
 75    wh_base[0] = .15
 76    wh_base[4] = wh_base[N-4] = .12
 77    wh_idea = wh_base.copy()
 78    for k in selected: wh_idea[k] *= .82/1.35
 79    # Fourier-domain DDE: each mode is independent; initialize broadband perturbation.
 80    def rollout(wh):
 81        delay = int(round(tau_d/dt)); n=int(duration/dt)
 82        z=np.zeros((n+delay+1,N), dtype=complex)
 83        z[:delay+1] = rng.normal(size=(delay+1,N)) + 1j*rng.normal(size=(delay+1,N))
 84        z[:delay+1] *= 1e-3 / np.sqrt(np.mean(np.abs(z[:delay+1])**2))
 85        for i in range(delay, delay+n):
 86            z[i+1] = z[i] + dt/tau_r*(-z[i] + wh*z[i-delay])
 87        norms=np.sqrt(np.mean(np.abs(z[delay:delay+n])**2,axis=1))
 88        return norms
 89    b=rollout(wh_base); s=rollout(wh_idea)
 90    def fit(norms):
 91        t=np.arange(len(norms))*dt
 92        sel=t>duration*.35
 93        return float(np.polyfit(t[sel], np.log(norms[sel]+1e-30),1)[0])
 94    return {"selected_modes": selected, "baseline_gain":1.35, "idea_gain":.82,
 95            "baseline_predicted_max_root":float(dominant(1.35,tau_d=tau_d).real),
 96            "idea_predicted_max_root":float(dominant(.82,tau_d=tau_d).real),
 97            "baseline_rollout_growth":fit(b), "idea_rollout_growth":fit(s),
 98            "baseline_final_over_initial":float(b[-1]/b[0]),
 99            "idea_final_over_initial":float(s[-1]/s[0]),
100            "unstable_baseline":bool(b[-1] > 10*b[0]), "stable_idea":bool(s[-1] < 10*b[0])}
101
102def main():
103    result={"seed":SEED, "math_verification":toy_verification(), "ring_comparison":stabilized_ring_demo()}
104    with open("results.json","w") as f: json.dump(result,f,indent=2)
105    print(json.dumps(result, indent=2))
106
107if __name__ == "__main__": main()