1import json, math
 2import numpy as np
 3
 4# Two independent symmetric binary coordinates, each flipping at rate q.
 5# Exact CTMC endpoint probability after T is known analytically; tau-leap uses
 6# Bernoulli(1-exp(-qh)) at each interval. This isolates discretization bias.
 7SEED=2485
 8q=0.8
 9T=1.0
10N=400000
11
12def exact_p(T):
13    return 0.5*(1.0-math.exp(-2*q*T))
14
15def tau_p(h):
16    n=round(T/h)
17    # Starting at 0, parity of n Bernoulli flips.
18    a=1.0-2.0*(-math.expm1(-q*h))
19    return 0.5*(1.0-a**n)
20
21def collision_scaling():
22    # Two competing edges from one coordinate, rates q1,q2. Probability both
23    # indicators activate in one leap is exactly (1-e^-hq1)(1-e^-hq2) ~ q1q2 h^2.
24    q1,q2=0.7,1.1
25    rows=[]
26    for h in [0.002,0.005,0.01,0.02,0.05,0.1]:
27        p=(1-math.exp(-h*q1))*(1-math.exp(-h*q2))
28        rows.append({'h':h,'collision_prob':p,'p_over_h2':p/h**2,'leading_prediction':q1*q2})
29    return rows
30
31def main():
32    rows=[]
33    exact=exact_p(T)
34    for h in [0.002,0.005,0.01,0.02,0.05,0.1,0.2]:
35        pred=tau_p(h)
36        # Monte Carlo confirms the Bernoulli chain prediction.
37        rng=np.random.default_rng(SEED+round(10000*h))
38        x=np.zeros(N,dtype=np.int8)
39        for _ in range(round(T/h)):
40            x ^= (rng.random(N)<(-np.expm1(-q*h)))
41        obs=float(x.mean())
42        rows.append({'h':h,'exact_endpoint_p':exact,'tau_predicted_p':pred,
43                     'tau_mc_p':obs,'abs_bias':abs(pred-exact),
44                     'bias_over_h':abs(pred-exact)/h})
45    out={'q':q,'T':T,'rows':rows,'collision_scaling':collision_scaling()}
46    with open('controlled_results.json','w') as f: json.dump(out,f,indent=2)
47    print(json.dumps(out,indent=2))
48
49if __name__=='__main__': main()