1import json, math
  2from pathlib import Path
  3import numpy as np
  4from scipy.linalg import expm
  5from scipy.stats import spearmanr
  6
  7SEED = 2272
  8rng = np.random.default_rng(SEED)
  9
 10# Stable damped oscillator: x(t)=exp(A t)x(0), y=c^T x.
 11A = np.array([[-0.08, -2.0], [2.0, -0.08]], dtype=float)
 12c = np.array([1.0, 0.0])
 13
 14def observability(delays):
 15    return np.stack([c @ expm(-A * float(t)) for t in delays])
 16
 17def metrics(O):
 18    s = np.linalg.svd(O, compute_uv=False)
 19    return float(s[-1]), float(s[0]), float(s[0]/s[-1])
 20
 21def reconstruct(O, x, sigma, reps=3000):
 22    # Gaussian measurement noise; average squared state error.
 23    pinv = np.linalg.pinv(O)
 24    noise = rng.normal(0.0, sigma, size=(reps, O.shape[0]))
 25    errs = (noise @ pinv.T)
 26    rmse = float(np.sqrt(np.mean(np.sum(errs*errs, axis=1))))
 27    return rmse
 28
 29def main():
 30    out = {"seed": SEED, "A": A.tolist(), "c": c.tolist(), "predictions": {}, "designs": [], "noise_sweep": [], "spearman": {}}
 31
 32    # Prediction 1: isotropic LS RMSE is sigma*sqrt(trace((O^T O)^-1)),
 33    # hence linear in noise level (and bounded by sqrt(n)*sigma/smin).
 34    delays = np.array([0.0, 0.30, 0.63, 0.95, 1.30])
 35    O = observability(delays)
 36    smin, smax, cond = metrics(O)
 37    trace_factor = math.sqrt(float(np.trace(np.linalg.inv(O.T @ O))))
 38    noise_levels = [0.002, 0.005, 0.01, 0.02, 0.05, 0.1]
 39    vals = []
 40    for sig in noise_levels:
 41        actual = reconstruct(O, None, sig)
 42        predicted = sig * trace_factor
 43        bound = math.sqrt(2.0) * sig / smin
 44        vals.append({"sigma": sig, "observed_rmse": actual, "predicted_rmse": predicted, "worst_case_bound": bound})
 45    slopes = np.polyfit(np.log(noise_levels), np.log([v['observed_rmse'] for v in vals]), 1)[0]
 46    out["predictions"]["noise_linear_scaling"] = {"predicted_slope": 1.0, "observed_loglog_slope": float(slopes), "design_smin": smin, "design_condition": cond}
 47    out["noise_sweep"] = vals
 48
 49    # Prediction 2: for fixed white noise, MSE^0.5 follows the exact inverse Gram factor,
 50    # and should be strongly monotone with ill-conditioning across delay sets.
 51    designs = {
 52        "uniform_short": np.linspace(0, 0.35, 5),
 53        "uniform_medium": np.linspace(0, 1.30, 5),
 54        "uniform_long": np.linspace(0, 2.60, 5),
 55        "irregular_good": np.array([0.0, 0.22, 0.61, 1.07, 1.56]),
 56        "irregular_bad_clustered": np.array([0.0, 0.02, 0.04, 0.07, 0.11]),
 57        "irregular_multiscale": np.array([0.0, 0.08, 0.40, 1.30, 2.35]),
 58    }
 59    sigma = 0.03
 60    factors, observed = [], []
 61    for name, d in designs.items():
 62        oo = observability(d)
 63        sm, sx, kk = metrics(oo)
 64        factor = math.sqrt(float(np.trace(np.linalg.inv(oo.T @ oo))))
 65        obs = reconstruct(oo, None, sigma)
 66        factors.append(factor); observed.append(obs)
 67        out["designs"].append({"name": name, "delays": d.tolist(), "smin": sm, "smax": sx, "condition": kk, "predicted_rmse": sigma*factor, "observed_rmse": obs})
 68    rho, p = spearmanr(factors, observed)
 69    out["spearman"] = {"rho_predicted_factor_vs_observed_rmse": float(rho), "p_value": float(p)}
 70    out["predictions"]["inverse_singular_value_noise_amplification"] = {"expected": "positive monotonic relation", "observed_spearman": float(rho), "threshold": 0.7}
 71
 72    # Prediction 3: adding a non-redundant delay can improve smin; sweep a fifth delay.
 73    # Compare a clustered baseline against candidate irregular placement.
 74    base = np.array([0.0, 0.25, 0.55, 0.85])
 75    candidates = np.linspace(0.0, 3.0, 301)
 76    smins = np.array([metrics(observability(np.r_[base, t]))[0] for t in candidates])
 77    best_i = int(np.argmax(smins))
 78    out["predictions"]["delay_design_sweep"] = {
 79        "predicted": "non-redundant delay increases smallest singular value",
 80        "base_delays": base.tolist(), "base_smin": metrics(observability(base))[0],
 81        "best_added_delay": float(candidates[best_i]), "best_smin": float(smins[best_i]),
 82        "improvement_ratio": float(smins[best_i]/metrics(observability(base))[0]),
 83        "candidate_grid_step": 0.01
 84    }
 85
 86    # Mini baseline comparison: actual delays are Poisson-irregular. The baseline
 87    # decodes with a fixed uniform-delay observation matrix, while the idea uses
 88    # the observed timestamps to build O for each history.
 89    nominal = np.array([0.0, 0.30, 0.60, 0.90, 1.20])
 90    baseline_O = observability(nominal)
 91    poisson_rows = []
 92    for mean_gap in [0.08, 0.15, 0.30, 0.50]:
 93        base_errors, conditioned_errors = [], []
 94        for _ in range(1200):
 95            # Past gaps, newest observation at delay zero.
 96            gaps = rng.exponential(mean_gap, size=4)
 97            actual_delays = np.r_[0.0, np.cumsum(gaps)]
 98            actual_O = observability(actual_delays)
 99            noise = rng.normal(0.0, 0.03, size=5)
100            x_true = rng.normal(size=2)
101            y = actual_O @ x_true + noise
102            base_errors.append(np.linalg.norm(np.linalg.pinv(baseline_O) @ y - x_true))
103            conditioned_errors.append(np.linalg.norm(np.linalg.pinv(actual_O) @ y - x_true))
104        poisson_rows.append({
105            "mean_gap": mean_gap,
106            "baseline_fixed_uniform_error": float(np.mean(base_errors)),
107            "conditioned_actual_delay_error": float(np.mean(conditioned_errors)),
108            "relative_improvement": float(1.0 - np.mean(conditioned_errors)/np.mean(base_errors))
109        })
110    out["poisson_baseline_comparison"] = poisson_rows
111
112    # A direct bound check over random noise vectors: ||pinv(O)e|| <= ||e||/smin.
113    bound_ratios = []
114    for _ in range(10000):
115        e = rng.normal(size=O.shape[0])
116        bound_ratios.append(np.linalg.norm(np.linalg.pinv(O) @ e) / (np.linalg.norm(e)/smin))
117    out["bound_check"] = {"max_ratio": float(max(bound_ratios)), "mean_ratio": float(np.mean(bound_ratios)), "claim": "ratio <= 1"}
118
119    # Pass criteria: mechanism checks, not merely a baseline win.
120    out["worked_checks"] = {
121        "noise_slope_within_0.10": bool(abs(slopes-1.0) <= 0.10),
122        "spearman_above_0.7": bool(rho >= 0.7),
123        "bound_holds": bool(max(bound_ratios) <= 1.0 + 1e-10),
124        "delay_improves_smin": bool(smins[best_i] > metrics(observability(base))[0] * 1.05),
125    }
126    out["worked"] = all(out["worked_checks"].values())
127    Path("results.json").write_text(json.dumps(out, indent=2))
128    print(json.dumps(out, indent=2))
129
130if __name__ == '__main__':
131    main()