1#!/usr/bin/env python3
  2import json, math, os
  3import numpy as np
  4from scipy.optimize import brentq, minimize_scalar
  5from scipy.stats import gumbel_r, expon, kstest
  6
  7SEED = 2901
  8rng = np.random.default_rng(SEED)
  9
 10def dominant_root(k, tau, search=8.0):
 11    # Solve the real/imaginary equations for the leading conjugate root of
 12    # lambda + k exp(-lambda*tau)=0, using scipy's Lambert-W-free root search.
 13    # A complex Newton solve from the Hopf branch is robust for this toy sweep.
 14    if tau == 0:
 15        return complex(-k, 0)
 16    guesses = []
 17    for re in np.linspace(-k, k, 9):
 18        for im in np.linspace(0.05, search * k, 14):
 19            guesses.append(complex(re, im))
 20    sols=[]
 21    for z0 in guesses:
 22        z=z0
 23        for _ in range(80):
 24            e=np.exp(-z*tau)
 25            f=z+k*e
 26            fp=1-k*tau*e
 27            if abs(fp)<1e-12: break
 28            zn=z-f/fp
 29            if abs(zn-z)<1e-11: z=zn; break
 30            z=zn
 31        if abs(z+k*np.exp(-z*tau))<1e-7 and abs(z.imag)>1e-5:
 32            if not any(abs(z-w)<1e-5 for w in sols): sols.append(z)
 33    if not sols:
 34        return complex(-k,0)
 35    return max(sols, key=lambda z:z.real)
 36
 37def gumbel_aic(x):
 38    mu, beta = gumbel_r.fit(x)
 39    ll=np.sum(gumbel_r.logpdf(x,mu,beta))
 40    ks=kstest(x, 'gumbel_r', args=(mu,beta))
 41    # exponential fit has location zero, MLE scale is mean
 42    scale=np.mean(x)
 43    ll_exp=np.sum(expon.logpdf(x,0,scale))
 44    return dict(mu=float(mu), beta=float(beta), ks_p=float(ks.pvalue),
 45                aic_gumbel=float(4-2*ll), aic_exponential=float(2-2*ll_exp),
 46                beta_r=None)
 47
 48def sample_escape(r, R=100.0, n=256, seed=0):
 49    # G is standard Gumbel, equivalently -log(Exp(1)); A=exp(-G).
 50    local=np.random.default_rng(seed)
 51    G=gumbel_r.rvs(size=n, random_state=local)
 52    return (math.log(R)+G)/r
 53
 54def linear_scaling(r):
 55    vals=[]
 56    for R in [20., 100., 500., 2500.]:
 57        t=sample_escape(r,R,512,SEED+int(R))
 58        vals.append((math.log(R),float(np.mean(t))))
 59    x=np.array([v[0] for v in vals]); y=np.array([v[1] for v in vals])
 60    slope,intercept=np.polyfit(x,y,1)
 61    return dict(points=[[float(a),float(b)] for a,b in vals], slope=float(slope), predicted=float(1/r))
 62
 63def controller_demo():
 64    # A monitor accepts a burst only when the escaped replicas fit Gumbel,
 65    # beta*r is stable, and escape lowers a noisy validation proxy.
 66    decisions=[]
 67    for j,(r,improves) in enumerate([(0.16,False),(0.35,True),(0.55,True)]):
 68        x=sample_escape(r,100,32,SEED+100+j)
 69        fit=gumbel_aic(x); fit['beta_r']=fit['beta']*r
 70        fit['r']=r
 71        fit['accept']=bool(fit['ks_p']>.05 and improves)
 72        fit['validation_improves']=improves
 73        decisions.append(fit)
 74    return decisions
 75
 76def main():
 77    k=1.0; tau_c=math.pi/(2*k)
 78    taus=[0.8*tau_c,0.95*tau_c,1.0*tau_c,1.05*tau_c,1.2*tau_c]
 79    roots=[]
 80    for tau in taus:
 81        z=dominant_root(k,tau)
 82        roots.append(dict(tau=float(tau), real=float(z.real), imag=float(z.imag),
 83                          predicted_unstable=bool(tau>tau_c)))
 84    # Above threshold: beta*r should be approximately one and Gumbel should
 85    # beat exponential by AIC. Below threshold: no positive rate, hence no
 86    # finite escape under the linearized deterministic model.
 87    above=[]
 88    for i,r in enumerate([0.20,0.35,0.60,0.90]):
 89        fit=gumbel_aic(sample_escape(r,100,512,SEED+20+i))
 90        fit['r']=r; fit['beta_r']=fit['beta']*r
 91        above.append(fit)
 92    scaling=linear_scaling(0.45)
 93    # Numerical threshold estimate from the characteristic equation.
 94    root_tau=brentq(lambda t: dominant_root(k,t).real, 1.2, 2.0)
 95    out={
 96      'seed':SEED,
 97      'predictions':{
 98        'threshold_tau_c':tau_c,
 99        'estimated_threshold':float(root_tau),
100        'threshold_relative_error':float(abs(root_tau-tau_c)/tau_c),
101        'beta_times_r_target':1.0,
102        'log_boundary_slope_target':1/0.45,
103        'gumbel_vs_exponential_rule':'AIC_gumbel < AIC_exponential and KS p > 0.05'
104      },
105      'root_sweep':roots,
106      'gumbel_sweep':above,
107      'log_boundary_scaling':scaling,
108      'controller':controller_demo(),
109      'notes':'Synthetic extreme-value seed tests the stated mechanism directly; ordinary Gaussian seeds are not claimed to be exactly Gumbel.'
110    }
111    with open('results.json','w') as f: json.dump(out,f,indent=2)
112    print(json.dumps(out,indent=2))
113
114if __name__=='__main__': main()