1import json, math
  2import numpy as np
  3
  4SEED = 2817
  5
  6def telegraph_run(gamma, lambdas, vs, dt=0.01, n=60000, burn=6000, seed=0):
  7    rng = np.random.default_rng(seed)
  8    lambdas, vs = np.asarray(lambdas, float), np.asarray(vs, float)
  9    s = rng.choice([-1., 1.], size=len(vs))
 10    u = 0.0
 11    a = math.exp(-gamma*dt)
 12    b = (1-a)/gamma
 13    p = -np.expm1(-lambdas*dt)
 14    out = np.empty(n-burn)
 15    maxabs = 0.
 16    j = 0
 17    for i in range(n):
 18        s[rng.random(len(s)) < p] *= -1
 19        u = a*u + b*float(np.dot(vs, s))
 20        maxabs = max(maxabs, abs(u))
 21        if i >= burn:
 22            out[j] = u; j += 1
 23    return out, maxabs
 24
 25def variance_pred(gamma, lambdas, vs):
 26    return float(np.sum(np.asarray(vs)**2/(gamma*(gamma+2*np.asarray(lambdas)))))
 27
 28def cov_correct(gamma, lam, v, tau):
 29    # Correct covariance of du=-gamma*u+v*s, E[s(t)s(t+tau)]=exp(-2 lam tau).
 30    tau = np.asarray(tau, float)
 31    if abs(gamma-2*lam) < 1e-8:
 32        # continuous limiting form
 33        return v*v*np.exp(-gamma*tau)*(1/(2*gamma) + tau/2)
 34    return v*v*(gamma*np.exp(-2*lam*tau)-2*lam*np.exp(-gamma*tau))/(gamma*(gamma**2-4*lam**2))
 35
 36def cov_as_given(gamma, lam, v, tau):
 37    tau=np.asarray(tau,float)
 38    return v*v*(np.exp(-2*lam*tau)-np.exp(-gamma*tau))/(gamma**2-4*lam**2)
 39
 40def acf(x, maxlag):
 41    x=x-x.mean(); den=np.dot(x,x)
 42    return np.array([np.dot(x[:-k] if k else x, x[k:] if k else x)/den for k in range(maxlag+1)])
 43
 44def excess_kurtosis(x):
 45    z=x-x.mean(); m2=np.mean(z*z)
 46    return float(np.mean(z**4)/m2**2-3)
 47
 48def quadratic_opt(noise_kind, dim=20, steps=5000, lr=.08, alpha=1., gamma=1., lam=.7, K=4, seed=1):
 49    rng=np.random.default_rng(seed); theta=rng.normal(2., .5, dim)
 50    u=np.zeros(dim); s=rng.choice([-1.,1.], size=(dim, K))
 51    # Match stationary per-coordinate telegraph variance with OU noise.
 52    v=1./math.sqrt(K)
 53    var=K*v*v/(gamma*(gamma+2*lam))
 54    ou=np.zeros(dim); oo=math.exp(-gamma)
 55    losses=[]
 56    for t in range(steps):
 57        g=theta
 58        if noise_kind=='telegraph':
 59            s[rng.random((dim,K)) < -np.expm1(-lam)] *= -1
 60            u=oo*u+(1-oo)*(v*np.sum(s,axis=1)/gamma)
 61            force=u
 62        elif noise_kind=='ou':
 63            ou=oo*ou+math.sqrt(var*(1-oo*oo))*rng.normal(size=dim)
 64            force=ou
 65        else: force=0.
 66        theta -= lr*(g+alpha*force)
 67        losses.append(float(np.mean(theta*theta)))
 68    return float(np.mean(losses[-500:])), float(np.mean(theta*theta))
 69
 70def main():
 71    gamma=1.; dt=.01
 72    rows=[]
 73    # Prediction 1: hard bound independent of lambda and time.
 74    for lam in [.1, .7, 2., 8.]:
 75        x,m=telegraph_run(gamma,[lam]*4,[.25]*4,dt=dt,seed=int(lam*100))
 76        rows.append({'lambda':lam,'bound':1.0,'max_abs':m,'ratio':m})
 77    # Prediction 2: exact variance scaling in lambda.
 78    variance=[]
 79    for lam in [.1,.3,.7,1.5,3.,8.]:
 80        x,_=telegraph_run(gamma,[lam]*4,[.25]*4,dt=dt,seed=10+int(lam*10))
 81        pred=variance_pred(gamma,[lam]*4,[.25]*4)
 82        variance.append({'lambda':lam,'emp_var':float(np.var(x)),'pred_var':pred,'rel_err':float(np.var(x)/pred-1)})
 83    # Prediction 3: covariance is sum of exponentials; compare corrected derivation and supplied expression.
 84    lam=.7; x,_=telegraph_run(gamma,[lam],[1.],dt=dt,seed=77)
 85    lags=np.array([0,10,25,50,100,200,400]); empirical=acf(x,int(lags[-1]))[lags]
 86    empirical_cov=empirical*np.var(x)
 87    corrected=cov_correct(gamma,lam,1.,lags*dt)
 88    supplied=cov_as_given(gamma,lam,1.,lags*dt)
 89    cov={'lags':(lags*dt).tolist(),'empirical':empirical_cov.tolist(),'corrected':corrected.tolist(),'supplied':supplied.tolist(),
 90         'rmse_corrected':float(np.sqrt(np.mean((empirical_cov-corrected)**2))),
 91         'rmse_supplied':float(np.sqrt(np.mean((empirical_cov-supplied)**2)))}
 92    # Prediction 4: many independent sources reduce standardized excess kurtosis.
 93    kurt=[]
 94    for K in [1,2,4,8,16]:
 95        x,_=telegraph_run(gamma,[.7]*K,[1/math.sqrt(K)]*K,dt=dt,seed=100+K)
 96        kurt.append({'K':K,'excess_kurtosis':excess_kurtosis(x),'variance':float(np.var(x)),
 97                     'pred_var':variance_pred(gamma,[.7]*K,[1/math.sqrt(K)]*K)})
 98    # Small secondary optimization sanity check (one telegraph source per coordinate).
 99    opt={k:quadratic_opt(k,seed=200) for k in ['none','telegraph','ou']}
100    result={'config':{'gamma':gamma,'dt':dt,'steps':60000,'burn':30000},'bound_sweep':rows,
101            'variance_sweep':variance,'covariance_check':cov,'kurtosis_sweep':kurt,'quadratic_optimizer':opt,
102            'note':'The covariance formula in the prompt is tested literally; it is zero at tau=0 and is not a valid covariance. corrected is derived from the same SDE.'}
103    with open('results.json','w') as f: json.dump(result,f,indent=2)
104    print(json.dumps(result,indent=2))
105
106if __name__=='__main__': main()