import json, math
import numpy as np

# Two independent symmetric binary coordinates, each flipping at rate q.
# Exact CTMC endpoint probability after T is known analytically; tau-leap uses
# Bernoulli(1-exp(-qh)) at each interval. This isolates discretization bias.
SEED=2485
q=0.8
T=1.0
N=400000

def exact_p(T):
    return 0.5*(1.0-math.exp(-2*q*T))

def tau_p(h):
    n=round(T/h)
    # Starting at 0, parity of n Bernoulli flips.
    a=1.0-2.0*(-math.expm1(-q*h))
    return 0.5*(1.0-a**n)

def collision_scaling():
    # Two competing edges from one coordinate, rates q1,q2. Probability both
    # indicators activate in one leap is exactly (1-e^-hq1)(1-e^-hq2) ~ q1q2 h^2.
    q1,q2=0.7,1.1
    rows=[]
    for h in [0.002,0.005,0.01,0.02,0.05,0.1]:
        p=(1-math.exp(-h*q1))*(1-math.exp(-h*q2))
        rows.append({'h':h,'collision_prob':p,'p_over_h2':p/h**2,'leading_prediction':q1*q2})
    return rows

def main():
    rows=[]
    exact=exact_p(T)
    for h in [0.002,0.005,0.01,0.02,0.05,0.1,0.2]:
        pred=tau_p(h)
        # Monte Carlo confirms the Bernoulli chain prediction.
        rng=np.random.default_rng(SEED+round(10000*h))
        x=np.zeros(N,dtype=np.int8)
        for _ in range(round(T/h)):
            x ^= (rng.random(N)<(-np.expm1(-q*h)))
        obs=float(x.mean())
        rows.append({'h':h,'exact_endpoint_p':exact,'tau_predicted_p':pred,
                     'tau_mc_p':obs,'abs_bias':abs(pred-exact),
                     'bias_over_h':abs(pred-exact)/h})
    out={'q':q,'T':T,'rows':rows,'collision_scaling':collision_scaling()}
    with open('controlled_results.json','w') as f: json.dump(out,f,indent=2)
    print(json.dumps(out,indent=2))

if __name__=='__main__': main()