1import json
 2import numpy as np
 3from hidden_diffusion_monitor import (analytic_spectrum, paper_sigma, simulate,
 4                                      forward_reverse_action_gap, autocorr, expm)
 5
 6
 7def run():
 8    rng = np.random.default_rng(20250308)
 9    A = np.array([[1., 1.], [-1., 0.]])  # A_yy=0, stable, coupled
10    Dx = Dm = 1.; dt = .03; n = 120000
11    dcs = [-.75, 0., .75]
12    omega = np.logspace(-2, 2, 250)
13    spectra = []; acs = []; rows = []
14    for dc in dcs:
15        D = np.array([[Dx, dc], [dc, Dm]])
16        spectra.append(analytic_spectrum(omega, A, D))
17        path, C, F, Q = simulate(A, D, dt, n, rng)
18        ac = autocorr(path, 80)[1:]
19        acs.append(ac)
20        sigma = paper_sigma(1., 1., Dx, Dm, dc)
21        gap = forward_reverse_action_gap(path, F, Q, dt)
22        exact = np.array([(expm(-A*t) @ C)[0, 0]/C[0, 0]
23                          for t in np.arange(1, 81)*dt])
24        rows.append({'D_xy': dc, 'spectrum_rel_to_baseline': float(
25            np.max(np.abs(spectra[-1]-spectra[0]) / np.maximum(spectra[0], 1e-12))),
26            'sigma_formula': float(sigma), 'monitor_gap': float(gap),
27            'autocorrelation_max_error': float(np.max(np.abs(ac-exact)))})
28    # Baseline has only observed x: its distinguishability is spectral/AC variation.
29    baseline_pairwise_ac = max(float(np.max(np.abs(acs[i]-acs[j])))
30                               for i in range(3) for j in range(i))
31    baseline_pairwise_spectrum = max(float(np.max(np.abs(spectra[i]-spectra[j]) /
32                                                  np.maximum(spectra[j], 1e-12)))
33                                     for i in range(3) for j in range(i))
34    monitor_gap_spread = max(r['monitor_gap'] for r in rows)-min(r['monitor_gap'] for r in rows)
35    formula_spread = max(r['sigma_formula'] for r in rows)-min(r['sigma_formula'] for r in rows)
36    result = {'rows': rows, 'baseline_observed_AC_pairwise_max': baseline_pairwise_ac,
37              'baseline_observed_spectrum_pairwise_max': baseline_pairwise_spectrum,
38              'monitor_gap_spread': monitor_gap_spread, 'formula_sigma_spread': formula_spread}
39    print(json.dumps(result, indent=2))
40    with open('results.json', 'w') as f: json.dump(result, f, indent=2)
41
42if __name__ == '__main__': run()