1import json
  2from pathlib import Path
  3import numpy as np
  4from scipy.optimize import least_squares
  5
  6DT = 0.02
  7H = 250
  8Q_TRUE = 1.0
  9
 10
 11def carrier(kind="sin", omega=2.0, phase=0.0):
 12    t = np.arange(H) * DT
 13    if kind == "sin":
 14        return np.sin(omega * t + phase)
 15    if kind == "binary":
 16        return np.sign(np.sin(omega * t + phase))
 17    raise ValueError(kind)
 18
 19
 20def rollout(a, m, q=Q_TRUE, u=None, z0=None):
 21    if u is None:
 22        u = np.zeros(H)
 23    if z0 is None:
 24        z0 = np.zeros(2)
 25    z = np.zeros((H + 1, 2))
 26    z[0] = z0
 27    for t in range(H):
 28        inp = u[t] + a * m[t]
 29        z[t + 1, 0] = z[t, 0] + DT * (-z[t, 0] + q * z[t, 1] ** 2 + inp)
 30        z[t + 1, 1] = z[t, 1] + DT * (-0.1 * z[t, 1] + inp)
 31    return z
 32
 33
 34def gramians(a, m, q=Q_TRUE, delta=1e-12):
 35    z = rollout(a, m, q)
 36    # C=[1,0], and A is the exact Jacobian of one Euler step.
 37    Phi = np.eye(2)
 38    Wo = np.zeros((2, 2))
 39    # Reachability is computed by propagating each input impulse forward.
 40    Wr = np.zeros((2, 2))
 41    for t in range(H):
 42        A = np.array([[1 - DT, 2 * DT * q * z[t, 1]],
 43                      [0, 1 - 0.1 * DT]])
 44        B = np.array([DT, DT])
 45        C = np.array([1.0, 0.0])
 46        Wo += (Phi.T @ np.outer(C, C) @ Phi) * DT
 47        Psi = np.eye(2)
 48        for j in range(t + 1, H):
 49            Aj = np.array([[1 - DT, 2 * DT * q * z[j, 1]],
 50                           [0, 1 - 0.1 * DT]])
 51            Psi = Aj @ Psi
 52        v = Psi @ B
 53        Wr += np.outer(v, v) * DT
 54        Phi = A @ Phi
 55    ev_o = np.linalg.eigvalsh(Wo)
 56    ev_r = np.linalg.eigvalsh(Wr)
 57    return Wo, Wr, ev_o, ev_r
 58
 59
 60def fit_q(observations, m, a):
 61    # Fit q from observed z1 only, with known input and known initial state.
 62    def residual(x):
 63        pred = rollout(a, m, float(x[0]))[:, 0]
 64        return (pred - observations) / 0.01
 65    return float(least_squares(residual, [0.2], bounds=(-3, 3)).x[0])
 66
 67
 68def main():
 69    rng = np.random.default_rng(123)
 70    m = carrier("sin", omega=2.0)
 71    amplitudes = np.array([0.0, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.6])
 72    rows = []
 73    for a in amplitudes:
 74        Wo, Wr, eo, er = gramians(a, m)
 75        rows.append({"a": float(a), "lambda_min_Wo": float(eo[0]),
 76                     "lambda_min_Wr": float(er[0]),
 77                     "logdet_Wo": float(np.linalg.slogdet(Wo + 1e-12*np.eye(2))[1])})
 78
 79    # Prediction 1: passive hidden observability is zero (up to numerical precision).
 80    passive = rows[0]["lambda_min_Wo"]
 81    # Prediction 2: small-amplitude hidden observability is quadratic in a.
 82    small = np.array([r for r in rows if 0.01 <= r["a"] <= 0.08])
 83    coef = float(np.polyfit(np.log(small[:, 0]) if False else np.log([r["a"] for r in small]),
 84                            np.log([r["lambda_min_Wo"] for r in small]), 1)[0])
 85    ratios = [r["lambda_min_Wo"] / (r["a"] ** 2) for r in small]
 86    # Prediction 3: a practical threshold is where the hidden eigenvalue exceeds 100x passive.
 87    threshold = next((r["a"] for r in rows if r["lambda_min_Wo"] > max(100*passive, 1e-8)), None)
 88
 89    # Secondary task check: estimate q from noisy observed trajectories.
 90    noise = rng.normal(0, 0.01, H + 1)
 91    passive_obs = rollout(0.0, m)[:, 0] + noise
 92    probed_a = 0.16
 93    probed_obs = rollout(probed_a, m)[:, 0] + rng.normal(0, 0.01, H + 1)
 94    q_passive = fit_q(passive_obs, m, 0.0)
 95    q_probed = fit_q(probed_obs, m, probed_a)
 96    result = {"dt": DT, "horizon": H*DT, "rows": rows,
 97              "predictions": {
 98                  "passive_lambda_hidden_expected": 0.0,
 99                  "passive_lambda_observed": passive,
100                  "quadratic_exponent_expected": 2.0,
101                  "quadratic_exponent_observed": coef,
102                  "quadratic_ratio_mean": float(np.mean(ratios)),
103                  "quadratic_ratio_cv": float(np.std(ratios)/np.mean(ratios)),
104                  "threshold_rule": "lambda_min(Wo) > max(100*passive, 1e-8)",
105                  "threshold_observed_a": threshold
106              },
107              "fit_q": {"true": Q_TRUE, "passive": q_passive, "probed": q_probed}}
108    Path("results.json").write_text(json.dumps(result, indent=2))
109    print(json.dumps(result, indent=2))
110
111
112if __name__ == "__main__":
113    main()