1import json
 2import numpy as np
 3
 4# Controlled stationary sweep for rho(x|a)=N(0,1/a), with Phi=a*x^2/2+log(a)/2.
 5# For adot=v, u(x,a)=-(v/(2a))*x solves div(u)-u*d_x Phi=v*d_a Phi.
 6
 7def covariance_check(seed=7, n=100000):
 8    rng=np.random.default_rng(seed); out=[]
 9    for a in (0.5, 1.0, 2.0):
10        x=rng.normal(size=n)/np.sqrt(a)
11        dphi=.5*(x*x-1/a)
12        cov=np.mean((x*x-np.mean(x*x))*(dphi-np.mean(dphi)))
13        out.append((a, float(-cov), float(-1/a**2), float(abs(-cov+1/a**2))))
14    return out
15
16def sweep(v, controlled, seed, dt=.01, n=10000, a0=.5, a1=2.0):
17    rng=np.random.default_rng(seed); x=rng.normal(size=n)/np.sqrt(a0)
18    a=a0; records=[]; steps=int(np.ceil((a1-a0)/(v*dt)))
19    for k in range(steps+1):
20        if k % max(1,steps//100)==0 or k==steps: records.append((a,np.mean(x*x),1/a))
21        if k==steps: break
22        rate=a+(v/(2*a) if controlled else 0.)
23        decay=np.exp(-rate*dt); variance=(1-decay*decay)/rate
24        x=decay*x+np.sqrt(variance)*rng.normal(size=n); a=min(a1,a+v*dt)
25    z=np.asarray(records); e=z[1:,1]-z[1:,2]
26    return z,float(np.mean(e)),float(np.sqrt(np.mean(e*e)))
27
28def deterministic_error(v, controlled, dt=.01, a0=.5, a1=2.0):
29    # Exact expected q=<x^2> recurrence for the frozen-coefficient OU update.
30    a=a0; q=1/a; steps=int(np.ceil((a1-a0)/(v*dt))); errors=[]
31    for _ in range(steps):
32        rate=a+(v/(2*a) if controlled else 0.)
33        d=np.exp(-rate*dt); q=d*d*q+(1-d*d)/rate
34        a=min(a1,a+v*dt); errors.append(q-1/a)
35    return float(np.mean(errors)),float(np.sqrt(np.mean(np.square(errors))))
36
37def main():
38    rates=np.array([.01,.02,.04,.08,.16])
39    sim=[]; det=[]
40    for i,v in enumerate(rates):
41        _,su,ru=sweep(v,False,100+i); _,sc,rc=sweep(v,True,200+i)
42        mu,ru_det=deterministic_error(v,False); mc,rc_det=deterministic_error(v,True)
43        sim.append((v,su,sc,ru,rc)); det.append((v,mu,mc,ru_det,rc_det))
44    sim=np.asarray(sim); det=np.asarray(det)
45    # RMS deterministic errors test the predicted O(v) versus O(v^2) scaling.
46    slope_u=float(np.polyfit(np.log(rates),np.log(det[:,3]),1)[0])
47    slope_c=float(np.polyfit(np.log(rates),np.log(det[:,4]),1)[0])
48    rng=np.random.default_rng(123); a=1.3; v=.17; x=rng.normal(size=100000)/np.sqrt(a)
49    lhs=-v/(2*a)-(-(v/(2*a))*x)*(a*x); rhs=v*.5*(x*x-1/a)
50    result={
51      'model':'rho=N(0,1/a), A=x^2, a: 0.5 -> 2.0',
52      'predictions':[
53        'covariance response is -1/a^2',
54        'the analytic control has zero continuity residual',
55        'uncontrolled tracking error scales O(v), controlled error O(v^2) near the adiabatic regime'
56      ],
57      'covariance_check':covariance_check(),
58      'simulation_rates_v_signed_uncontrolled_controlled_rms_uncontrolled_controlled':sim.tolist(),
59      'deterministic_rates_v_mean_uncontrolled_controlled_rms_uncontrolled_controlled':det.tolist(),
60      'deterministic_loglog_rms_slope_uncontrolled':slope_u,
61      'deterministic_loglog_rms_slope_controlled':slope_c,
62      'continuity_residual_rms':float(np.sqrt(np.mean((lhs-rhs)**2))),
63      'fastest_rate_controlled_to_uncontrolled_rms_ratio':float(det[-1,4]/det[-1,3])
64    }
65    with open('results.json','w') as f: json.dump(result,f,indent=2)
66    print(json.dumps(result,indent=2))
67
68if __name__=='__main__': main()