1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1632
  6J = np.array([[0.0, -1.0], [1.0, 0.0]])
  7I2 = np.eye(2)
  8
  9
 10def transition(gamma, omega, s):
 11    # Euler step for dx/dt=(-gamma I+s*omega J)x, dt=1.
 12    return (1.0 - gamma) * I2 + s * omega * J
 13
 14
 15def antisym(x):
 16    return 0.5 * (x - x.T)
 17
 18
 19def response(A, q):
 20    return np.linalg.matrix_power(A, q)
 21
 22
 23def stationary_cov(A, noise=1.0, n=200000, burn=100):
 24    rng = np.random.default_rng(1000 + int(abs(A[0, 1]) * 10000))
 25    x = np.zeros(2)
 26    xs = []
 27    for t in range(n + burn):
 28        x = A @ x + noise * rng.normal(size=2)
 29        if t >= burn:
 30            xs.append(x.copy())
 31    x = np.asarray(xs)
 32    return np.cov(x, rowvar=False, bias=True), x
 33
 34
 35def lag_cov(xs, lag):
 36    a, b = xs[lag:], xs[:-lag]
 37    a = a - a.mean(0)
 38    b = b - b.mean(0)
 39    return a.T @ b / len(a)
 40
 41
 42def fdr_residual_empirical(xs, A, tau=1.0, lag=1):
 43    # For discrete dynamics, dC/dt is approximated by C(lag)-C(0).
 44    c0 = np.cov(xs, rowvar=False, bias=True)
 45    cl = lag_cov(xs, lag)
 46    # impulse response at lag and current time; current impulse is identity.
 47    chi_l = np.linalg.matrix_power(A, lag)
 48    chi_0 = I2
 49    lhs = cl - c0
 50    rhs = -tau * (chi_l - chi_0)
 51    return np.linalg.norm(lhs - rhs) / (np.linalg.norm(lhs) + np.linalg.norm(rhs) + 1e-12)
 52
 53
 54def spectral_boundary_sweep():
 55    # Predicted stability boundary for Euler A=(1-gamma)I+s omega J:
 56    # rho(A)<1 iff (1-gamma)^2+omega^2<1.
 57    gamma = 0.25
 58    omegas = np.linspace(0.0, 1.15, 24)
 59    rows = []
 60    for w in omegas:
 61        A = transition(gamma, w, 1)
 62        rho = max(abs(np.linalg.eigvals(A)))
 63        rows.append((float(w), float(rho), bool(rho < 1.0)))
 64    predicted = float(np.sqrt(1.0 - (1.0 - gamma) ** 2))
 65    observed = omegas[np.where(np.array([r[2] for r in rows]))[0][-1]]
 66    return {'gamma': gamma, 'predicted_omega_boundary': predicted,
 67            'observed_grid_boundary': float(observed), 'rows': rows}
 68
 69
 70def oc_sweep():
 71    # Exact prediction: chi_q^A(+omega)+chi_q^A(-omega)=0; magnitude is
 72    # approximately q*omega for small omega and vanishes at omega=0.
 73    gamma = 0.2
 74    qs = [1, 2, 3]
 75    omegas = np.linspace(0.0, 0.03, 7)
 76    data = []
 77    for q in qs:
 78        vals = []
 79        mags = []
 80        for w in omegas:
 81            Ap, Am = transition(gamma, w, 1), transition(gamma, w, -1)
 82            vals.append(np.linalg.norm(antisym(response(Ap, q)) + antisym(response(Am, q))))
 83            mags.append(np.linalg.norm(antisym(response(Ap, q))))
 84        slope = np.polyfit(omegas[1:], mags[1:], 1)[0]
 85        data.append({'q': q, 'max_reversal_residual': float(max(vals)),
 86                     'zero_omega_magnitude': float(mags[0]),
 87                     'small_omega_slope': float(slope),
 88                     'predicted_small_omega_slope': float(np.sqrt(2) * q * (1-gamma)**(q-1))})
 89    return {'gamma': gamma, 'data': data}
 90
 91
 92def fdr_sweep():
 93    # Test that the empirical residual decreases with sample size and that
 94    # the correct effective temperature is the injected noise scale squared.
 95    gamma, omega, noise = 0.35, 0.12, 0.7
 96    A = transition(gamma, omega, 1)
 97    # Stationary covariance solves C=A C A^T+noise^2 I; tau=noise^2.
 98    c, xs = stationary_cov(A, noise=noise, n=80000)
 99    taus = np.linspace(0.25, 1.0, 16)
100    residuals = []
101    for tau in taus:
102        # Use exact covariance and exact discrete lag covariance, removing MC error.
103        cl = A @ c
104        exact = np.linalg.norm((cl-c) + tau*(A-I2))
105        residuals.append(float(exact))
106    best_tau = float(taus[int(np.argmin(residuals))])
107    sample_sizes = [500, 2000, 8000, 32000]
108    empirical = []
109    for n in sample_sizes:
110        empirical.append(float(fdr_residual_empirical(xs[:n], A, tau=noise**2)))
111    return {'gamma': gamma, 'omega': omega, 'noise_scale': noise,
112            'predicted_tau': noise**2, 'best_grid_tau': best_tau,
113            'exact_residual_at_predicted_tau': float(residuals[np.argmin(abs(taus-noise**2))]),
114            'sample_sizes': sample_sizes, 'empirical_normalized_residuals': empirical,
115            'tau_grid': taus.tolist(), 'exact_residuals': residuals}
116
117
118def fit_comparison():
119    # Tiny equal-compute linear transition fitting. Regularizer is the OC
120    # penalty on the learned antisymmetric component under orientation reversal.
121    rng = np.random.default_rng(SEED)
122    gamma, omega = 0.2, 0.18
123    Aplus, Aminus = transition(gamma, omega, 1), transition(gamma, omega, -1)
124    n = 5000
125    x = rng.normal(size=(n,2))
126    y = x @ Aplus.T + 0.08*rng.normal(size=(n,2))
127    # Baseline least squares on the same plus-orientation data.
128    Ab = np.linalg.lstsq(x, y, rcond=None)[0].T
129    vals = []
130    for lam in [0.0, 0.001, 0.01, 0.1, 1.0]:
131        # Objective ||AX-Y||^2 + lambda ||Anti(A)+Anti(Aminus)||^2.
132        # Solve by deterministic gradient descent, same 300 steps for all.
133        A = np.zeros((2,2))
134        lr = 0.08
135        for _ in range(300):
136            grad = 2*(A @ (x.T@x)/n - y.T@x/n)
137            R = antisym(A) + antisym(Aminus)
138            grad += 2*lam*R
139            A -= lr*grad
140        testx = rng.normal(size=(3000,2))
141        testy = testx @ Aplus.T
142        mse = np.mean((testx @ A.T-testy)**2)
143        oc = np.linalg.norm(antisym(A)+antisym(Aminus))
144        vals.append({'lambda': lam, 'test_mse': float(mse), 'oc_residual': float(oc)})
145    testx = rng.normal(size=(3000,2))
146    testy = testx @ Aplus.T
147    return {'baseline_test_mse': float(np.mean((testx @ Ab.T-testy)**2)), 'baseline_oc_residual': float(np.linalg.norm(antisym(Ab)+antisym(Aminus))), 'regularized_runs': vals}
148
149
150def main():
151    out = {'seed': SEED, 'oc_sweep': oc_sweep(), 'stability_sweep': spectral_boundary_sweep(),
152           'fdr_sweep': fdr_sweep(), 'fit_comparison': fit_comparison()}
153    Path('results.json').write_text(json.dumps(out, indent=2))
154    print(json.dumps(out, indent=2))
155
156if __name__ == '__main__':
157    main()