1import json
  2import numpy as np
  3from scipy.linalg import expm
  4from pathlib import Path
  5
  6
  7def block(gf, gs, k):
  8    return np.array([[-gf, k], [k, -gs]], dtype=float)
  9
 10
 11def euler(A, h, n, x0=None, forcing=None):
 12    x = np.zeros((n + 1, A.shape[0]))
 13    if x0 is not None:
 14        x[0] = x0
 15    for t in range(n):
 16        u = 0.0 if forcing is None else forcing[t]
 17        b = np.ones(A.shape[0]) * u
 18        x[t + 1] = x[t] + h * (A @ x[t] + b)
 19    return x
 20
 21
 22def frequency_transfer(A, w, b=None, c=None):
 23    b = np.ones(A.shape[0]) if b is None else np.asarray(b)
 24    c = np.ones(A.shape[0]) if c is None else np.asarray(c)
 25    z = 1j * w
 26    return c @ np.linalg.solve(z * np.eye(A.shape[0]) - A, b)
 27
 28
 29def decay_rates_from_impulse(A, dt=0.002, duration=12.0):
 30    # Independently simulate exact sampled propagation, then fit each modal
 31    # amplitude. This tests decay rather than simply reporting eig(A).
 32    vals, vecs = np.linalg.eigh(A)
 33    x = np.ones(2)
 34    n = int(duration / dt)
 35    modal = np.empty((n + 1, 2))
 36    modal[0] = np.linalg.solve(vecs, x)
 37    P = expm(A * dt)
 38    for i in range(n):
 39        x = P @ x
 40        modal[i + 1] = np.linalg.solve(vecs, x)
 41    t = np.arange(n + 1) * dt
 42    fitted = []
 43    for j in range(2):
 44        # Avoid the final numerical underflow tail and fit log magnitude.
 45        keep = (np.abs(modal[:, j]) > 1e-10) & (t > 0.1) & (t < duration * .8)
 46        slope = np.polyfit(t[keep], np.log(np.abs(modal[keep, j])), 1)[0]
 47        fitted.append(-slope)
 48    return np.sort(-vals), np.sort(np.asarray(fitted))
 49
 50
 51def least_squares_reconstruction(gf=10.0, gs=1.0, k=0.15, seed=0):
 52    rng = np.random.default_rng(seed)
 53    n_train, n_test = 3000, 1500
 54    t = np.arange(n_train + n_test)
 55    y = (np.sin(2 * np.pi * 0.025 * t) + 0.7 * np.sin(2 * np.pi * 0.22 * t + .4))
 56    x = y + 0.03 * rng.standard_normal(len(t))
 57    # Discrete stable dynamics; each state is a causal filtered view of x.
 58    h = 0.03
 59    Ad = np.eye(2) + h * block(gf, gs, k)
 60    states = np.zeros((len(x) + 1, 2))
 61    for i, u in enumerate(x):
 62        states[i + 1] = Ad @ states[i] + h * np.array([u, u])
 63    # Baseline: ordinary one-timescale slow state, matched parameter count as
 64    # closely as possible by fitting its readout on the same training segment.
 65    sb = np.zeros(len(x) + 1)
 66    ad = 1 - h * gs
 67    for i, u in enumerate(x):
 68        sb[i + 1] = ad * sb[i] + h * u
 69    def fit_predict(features):
 70        F = np.column_stack([features, np.ones(len(features))])
 71        coef = np.linalg.lstsq(F[:n_train], y[:n_train], rcond=None)[0]
 72        return F[n_train:] @ coef
 73    pred_dual = fit_predict(states[1:])
 74    pred_single = fit_predict(sb[1:, None])
 75    return float(np.mean((pred_single - y[n_train:]) ** 2)), float(np.mean((pred_dual - y[n_train:]) ** 2))
 76
 77
 78def main():
 79    out = {}
 80    # Prediction 1: symmetric-part stability changes sign at k=sqrt(gf*gs).
 81    gf, gs = 10.0, 1.0
 82    kcrit = np.sqrt(gf * gs)
 83    ks = np.linspace(0, 1.5 * kcrit, 301)
 84    lmax = np.array([np.linalg.eigvalsh((block(gf, gs, k) + block(gf, gs, k).T) / 2).max() for k in ks])
 85    crossing = ks[np.argmin(np.abs(lmax))]
 86    out['stability_prediction'] = {
 87        'predicted_kcrit': float(kcrit), 'observed_kcrit_grid': float(crossing),
 88        'relative_error': float(abs(crossing-kcrit)/kcrit),
 89        'stable_at_0.8kcrit': bool(np.linalg.eigvalsh(block(gf, gs, .8*kcrit)).max() <= 1e-12),
 90        'unstable_at_1.2kcrit': bool(np.linalg.eigvalsh(block(gf, gs, 1.2*kcrit)).max() > 1e-12),
 91    }
 92    # Prediction 2: a gap yields two decay rates, with timescale ratio tracking
 93    # the imposed gap (small exchange perturbation).
 94    gap_rows = []
 95    for ratio in [2, 5, 10, 20, 40]:
 96        gs0, gf0 = 1.0, float(ratio)
 97        kk = .05 * np.sqrt(gf0 * gs0)
 98        rates, fitted = decay_rates_from_impulse(block(gf0, gs0, kk))
 99        gap_rows.append({'ratio': ratio, 'predicted_rates': rates.tolist(), 'measured_rates': fitted.tolist(),
100                         'measured_rate_ratio': float(fitted[1]/fitted[0])})
101    out['gap_prediction'] = gap_rows
102    # Prediction 3: slow-only reduction loses the fast response around w~gf.
103    A = block(10., 1., .15)
104    freqs = np.logspace(-2, 2, 300)
105    full = np.array([abs(frequency_transfer(A, w)) for w in freqs])
106    slow = np.array([abs(1 / (1j*w + 1)) for w in freqs])
107    relerr = np.abs(full-slow) / np.maximum(full, 1e-12)
108    low = float(np.mean(relerr[freqs < 0.3]))
109    high = float(np.mean(relerr[freqs > 10.]))
110    near = float(np.mean(relerr[(freqs > 7) & (freqs < 14)]))
111    out['frequency_prediction'] = {'fast_rate': 10.0, 'observed_peak_error_frequency': float(freqs[np.argmax(relerr)]),
112                                   'mean_relative_error_low_w': low, 'mean_relative_error_near_fast_rate': near,
113                                   'mean_relative_error_high_w': high, 'error_ratio_near_vs_low': near/max(low,1e-12)}
114    single_mse, dual_mse = least_squares_reconstruction()
115    out['toy_reconstruction'] = {'single_timescale_mse': single_mse, 'dual_band_mse': dual_mse,
116                                 'relative_mse_change': (single_mse-dual_mse)/single_mse}
117    Path('results.json').write_text(json.dumps(out, indent=2))
118    print(json.dumps(out, indent=2))
119
120if __name__ == '__main__':
121    main()