1import json, math
  2from pathlib import Path
  3import numpy as np
  4from scipy.optimize import minimize
  5
  6# Two stable limit cycles are represented by stable radii r=1 and r=2,
  7# with a common angular velocity. The separatrix is r=1.5.
  8R1, RS, R2 = 1.0, 1.5, 2.0
  9OMEGA = 1.0
 10DT = 0.01
 11
 12def fr(r):
 13    return -(r - R1) * (r - RS) * (r - R2)
 14
 15def dfr(r):
 16    return -(3*r*r - 2*(R1+RS+R2)*r + (R1*RS + R1*R2 + RS*R2))
 17
 18def action_integral(a, b, n=100000):
 19    x = np.linspace(a, b, n)
 20    return float(np.trapz(-fr(x), x))
 21
 22def discrete_action(x, T):
 23    # B=1, so this is exactly the stated 1/4 integral discretization.
 24    h = T / (len(x)-1)
 25    v = np.diff(x) / h
 26    mid = 0.5*(x[:-1] + x[1:])
 27    return float(h * np.sum((v - fr(mid))**2) / 4.0)
 28
 29def optimize_action(a, b, T=20.0, n=81):
 30    # Fixed endpoints and a smooth monotone initial path. The optimizer can
 31    # discover the minimum-action uphill path for the finite endpoint problem.
 32    x0 = np.linspace(a, b, n)
 33    h = T/(n-1)
 34    def fun(y):
 35        x = np.r_[a, y, b]
 36        return discrete_action(x, T)
 37    ans = minimize(fun, x0[1:-1], method='L-BFGS-B', options={'maxiter': 1200, 'ftol': 1e-12})
 38    return float(ans.fun), bool(ans.success)
 39
 40def simulate_rates(D, paths=160, steps=125000, seed=0):
 41    # Vectorized first-passage counting. Hysteretic labels avoid recrossing noise.
 42    rng = np.random.default_rng(seed)
 43    r = np.full(paths, R1)
 44    state = np.ones(paths, dtype=np.int8)
 45    counts = np.zeros(2, dtype=np.int64)
 46    for _ in range(steps):
 47        r += fr(r)*DT + math.sqrt(2*D*DT)*rng.standard_normal(paths)
 48        # radial coordinate is reflected at zero (irrelevant for rare events)
 49        r = np.abs(r)
 50        hit2 = (state == 1) & (r >= RS)
 51        hit1 = (state == 2) & (r <= RS)
 52        counts[0] += np.count_nonzero(hit2)
 53        counts[1] += np.count_nonzero(hit1)
 54        state[hit2] = 2
 55        state[hit1] = 1
 56    total_time = paths * steps * DT
 57    return counts / total_time, counts
 58
 59def main():
 60    # Prediction 1: action barrier from cycle 1 to separatrix.
 61    barrier = action_integral(R1, RS)
 62    reverse_barrier = action_integral(R2, RS)
 63    # Prediction 2: transverse Floquet multiplier exp(f'(cycle)*period).
 64    period = 2*math.pi/OMEGA
 65    floquet1 = math.exp(dfr(R1)*period)
 66    floquet2 = math.exp(dfr(R2)*period)
 67    opt_action, opt_ok = optimize_action(1.02, 1.48)
 68    finite_barrier = action_integral(1.02, 1.48)
 69
 70    # Prediction 3: log k is affine in 1/D in the rare-event regime, while
 71    # sufficiently large D leaves that regime. Five-plus values are required.
 72    Ds = np.array([0.004, 0.006, 0.010, 0.015, 0.030, 0.080, 0.200])
 73    rates = []
 74    counts = []
 75    for i, D in enumerate(Ds):
 76        rr, cc = simulate_rates(float(D), paths=48, steps=50000, seed=100+i)
 77        rates.append(rr)
 78        counts.append(cc.tolist())
 79    rates = np.asarray(rates)
 80    counts = np.asarray(counts)
 81    # Fit only forward rates with enough events; report rare subset and all-data fits.
 82    valid = (counts[:,0] >= 8) & (rates[:,0] > 0)
 83    rare = valid & (Ds <= 0.015)
 84    high = valid & (Ds >= 0.080)
 85    def fit(mask):
 86        if np.count_nonzero(mask) < 2: return [float('nan')]*3
 87        p = np.polyfit(1/Ds[mask], np.log(rates[mask,0]), 1)
 88        pred = np.polyval(p, 1/Ds[mask])
 89        r2 = 1 - np.sum((np.log(rates[mask,0])-pred)**2)/np.sum((np.log(rates[mask,0])-np.mean(np.log(rates[mask,0])))**2)
 90        return [float(p[0]), float(p[1]), float(r2)]
 91    fit_rare = fit(rare)
 92    fit_all = fit(valid)
 93    fit_high = fit(high)
 94    crossover_D = (RS-R1)**2/(2*period)
 95    result = {
 96      'model': 'dr=-(r-1)(r-1.5)(r-2)dt + sqrt(2D)dW; dtheta=dt',
 97      'predictions': {
 98        'action_barrier_1_to_separatrix': {'predicted': barrier, 'estimated_discrete_path_1.02_to_1.48': opt_action, 'finite_endpoint_integral': finite_barrier, 'relative_error_vs_finite_integral': abs(opt_action-finite_barrier)/finite_barrier},
 99        'floquet_transverse': {'period': period, 'cycle_1_derivative': dfr(R1), 'cycle_2_derivative': dfr(R2), 'cycle_1_multiplier': floquet1, 'cycle_2_multiplier': floquet2, 'neutral_phase_multiplier': 1.0},
100        'arrhenius_slope': {'predicted_slope': -barrier, 'rare_fit_slope': fit_rare[0], 'rare_fit_intercept': fit_rare[1], 'rare_fit_R2': fit_rare[2], 'all_fit_slope': fit_all[0], 'all_fit_R2': fit_all[2], 'high_noise_fit_slope': fit_high[0], 'high_noise_fit_R2': fit_high[2], 'predicted_crossover_D_from_cycle_width': crossover_D}
101      },
102      'sweep': [{'D':float(D), 'forward_rate':float(rates[i,0]), 'reverse_rate':float(rates[i,1]), 'forward_count':int(counts[i,0]), 'reverse_count':int(counts[i,1]), 'x=1/D':float(1/D)} for i,D in enumerate(Ds)],
103      'settings': {'dt':DT, 'paths':48, 'steps':50000, 'total_time_each_D':48*50000*DT, 'seed_base':100},
104      'mechanism_checks': {'action_agreement_pass': abs(opt_action-finite_barrier)/finite_barrier < .20, 'floquet_stable_pass': floquet1 < 1 and floquet2 < 1, 'rare_arrhenius_linear_pass': fit_rare[2] > .90, 'large_noise_deviation_from_rare_slope': abs(fit_high[0]-fit_rare[0]) > .20*abs(fit_rare[0])}
105    }
106    Path('results.json').write_text(json.dumps(result, indent=2))
107    print(json.dumps(result, indent=2))
108
109if __name__ == '__main__':
110    main()