1import numpy as np
  2from scipy.linalg import expm, solve_continuous_lyapunov
  3
  4
  5def stationary_cov(A, D):
  6    # A C + C A^T = 2D for dX=-A X dt + sqrt(2D)dW
  7    return -solve_continuous_lyapunov(A, -2.0 * D)
  8
  9
 10def exact_step(A, D, dt):
 11    C = stationary_cov(A, D)
 12    F = expm(-A * dt)
 13    Q = (C - F @ C @ F.T)
 14    Q = (Q + Q.T) / 2
 15    vals, vecs = np.linalg.eigh(Q)
 16    Q = vecs @ np.diag(np.maximum(vals, 1e-12)) @ vecs.T
 17    return C, F, Q
 18
 19
 20def simulate(A, D, dt, n, rng):
 21    C, F, Q = exact_step(A, D, dt)
 22    out = np.empty((n, 2))
 23    out[0] = rng.multivariate_normal(np.zeros(2), C)
 24    noise = rng.multivariate_normal(np.zeros(2), Q, size=n - 1)
 25    for k in range(n - 1):
 26        out[k + 1] = F @ out[k] + noise[k]
 27    return out, C, F, Q
 28
 29
 30def gaussian_logpdf(v, cov):
 31    sign, logdet = np.linalg.slogdet(cov)
 32    inv = np.linalg.inv(cov)
 33    return -.5 * (logdet + v @ inv @ v + 2*np.log(2*np.pi))
 34
 35
 36def forward_reverse_action_gap(path, F, Q, dt):
 37    """Discretized path-action gap using the forward transition both ways.
 38    This is the standard forward path versus its time-reversed ordering;
 39    unlike a fitted conditional-entropy difference it detects probability currents.
 40    """
 41    vals = []
 42    for x, y in zip(path[:-1], path[1:]):
 43        vals.append(gaussian_logpdf(y - F @ x, Q) -
 44                    gaussian_logpdf(x - F @ y, Q))
 45    return float(np.mean(vals) / dt)
 46
 47
 48def analytic_spectrum(omega, A, D):
 49    out = np.empty_like(omega, dtype=float)
 50    for i, w in enumerate(omega):
 51        H = np.linalg.inv(A + 1j * w * np.eye(2))
 52        out[i] = np.real(2 * (H @ D @ H.conj().T)[0, 0])
 53    return out
 54
 55
 56def paper_sigma(a, b, Dx, Dm, Dc):
 57    return (a * (Dc + Dm) + b * Dx) ** 2 / (a * (Dx * Dm - Dc * Dc))
 58
 59
 60def autocorr(path, maxlag):
 61    z = path[:, 0] - path[:, 0].mean()
 62    den = np.dot(z, z)
 63    return np.array([np.dot(z[:-k] if k else z, z[k:] if k else z) /
 64                     (den if k == 0 else np.sqrt(np.dot(z[:-k], z[:-k]) * np.dot(z[k:], z[k:])))
 65                     for k in range(maxlag + 1)])
 66
 67
 68def main():
 69    rng = np.random.default_rng(20250308)
 70    # Exact adaptation: A_yy=0 and nontrivial coupling. Stable since det(A)>0.
 71    a, b, g = 1.0, 1.0, -1.0
 72    A = np.array([[a, b], [g, 0.0]])
 73    Dx = Dm = 1.0
 74    dcs = [-0.75, 0.0, 0.75]
 75    omega = np.logspace(-2, 2, 300)
 76    reference = analytic_spectrum(omega, A, np.diag([Dx, Dm]))
 77    print('A=', A.tolist())
 78    print('Analytic blind-direction check (max relative PSD error):')
 79    rows = []
 80    for dc in dcs:
 81        D = np.array([[Dx, dc], [dc, Dm]])
 82        spec = analytic_spectrum(omega, A, D)
 83        rel = np.max(np.abs(spec - reference) / np.maximum(reference, 1e-12))
 84        predicted = paper_sigma(a, b, Dx, Dm, dc)
 85        path, C, F, Q = simulate(A, D, dt=.03, n=180000, rng=rng)
 86        gap = forward_reverse_action_gap(path, F, Q, .03)
 87        # Compare observed autocorrelation to its exact stationary prediction.
 88        lags = np.arange(1, 101) * .03
 89        exact_ac = np.array([(expm(-A*t) @ C)[0, 0] / C[0, 0] for t in lags])
 90        sample_ac = autocorr(path, 100)[1:]
 91        ac_err = float(np.max(np.abs(sample_ac - exact_ac)))
 92        rows.append((dc, rel, predicted, gap, ac_err))
 93    sigmas = np.array([r[2] for r in rows]); gaps = np.array([r[3] for r in rows])
 94    print(' dc     analytic_rel_err    sigma_formula    action_gap       max_AC_error')
 95    for r in rows:
 96        print('% .2f   % .3e          % .6f       % .6f        % .4f' % r)
 97    print('formula sigma range:', float(sigmas.min()), float(sigmas.max()))
 98    print('action-gap range:', float(gaps.min()), float(gaps.max()))
 99    print('spectrum remains blind:', max(r[1] for r in rows) < 1e-10)
100    print('autocorrelation agreement:', max(r[4] for r in rows) < .03)
101    print('gap spread / formula spread:', float(np.ptp(gaps)), float(np.ptp(sigmas)))
102
103if __name__ == '__main__':
104    main()