1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 1419
  6RHO = np.array([1.0, 1.25, 1.5, 2.0, 3.0, 4.0])
  7
  8
  9def orthonormal_kernel(C):
 10    # SVD gives an orthonormal basis of ker(C), with C shape (p,n).
 11    u, s, vh = np.linalg.svd(C, full_matrices=True)
 12    rank = int(np.sum(s > 1e-10))
 13    return vh[rank:].T
 14
 15
 16def projected_sweep(J, C, G0, P):
 17    U = orthonormal_kernel(C)
 18    rows = []
 19    for rho in RHO:
 20        H = U.T @ (J - rho * G0 @ C) @ U
 21        S = P @ H + H.T @ P
 22        rows.append({
 23            "rho": float(rho),
 24            "H": H.tolist(),
 25            "lambda_max_sym": float(np.linalg.eigvalsh((S + S.T) / 2).max()),
 26        })
 27    return U, rows
 28
 29
 30def scheduler(rows, alpha=0.05, P=None):
 31    # The stated formula uses -2 alpha lambda_min(P).
 32    pmin = 1.0 if P is None else float(np.linalg.eigvalsh(P).min())
 33    threshold = -2 * alpha * pmin
 34    for r in rows:
 35        if r["lambda_max_sym"] <= threshold:
 36            return r["rho"]
 37    return float(rows[-1]["rho"])
 38
 39
 40def noisy_observer(J, C, G0, rho_policy, steps=3000, dt=0.01,
 41                   obs_every=5, noise=0.08, seed=0):
 42    rng = np.random.default_rng(seed)
 43    n = J.shape[0]
 44    x = np.array([1.0, 1.0])
 45    z = np.array([-1.0, 0.0])
 46    errors = []
 47    used = []
 48    for t in range(steps):
 49        x = x + dt * (J @ x)
 50        if t % obs_every == 0:
 51            y = C @ x + noise * rng.normal(size=(C.shape[0],))
 52            rho = float(rho_policy(z, J, C, G0))
 53            z = z + dt * (J @ z + rho * G0 @ (y - C @ z))
 54            used.append(rho)
 55        else:
 56            z = z + dt * (J @ z)
 57        errors.append(float(np.linalg.norm(x - z)))
 58    tail = np.asarray(errors[-1000:])
 59    return {"rmse_tail": float(np.sqrt(np.mean(tail ** 2))),
 60            "mean_rho": float(np.mean(used)), "max_error": float(np.max(errors))}
 61
 62
 63def main():
 64    rng = np.random.default_rng(SEED)
 65    # C observes x1; K is exactly the x2 axis. G0 injects the measurement into x1.
 66    C = np.array([[1.0, 0.0]])
 67    G0 = np.array([[1.0], [0.35]])
 68    P = np.eye(1)
 69    U = orthonormal_kernel(C)
 70
 71    # Prediction 1: C U_K=0, hence H(rho) is invariant in rho to machine precision.
 72    J_stable = np.array([[-0.5, 0.0], [0.0, -0.20]])
 73    _, stable_rows = projected_sweep(J_stable, C, G0, P)
 74    Hs = np.asarray([r["H"] for r in stable_rows]).reshape(len(RHO), -1)
 75    invariant_error = float(np.max(np.abs(Hs - Hs[0])))
 76
 77    # Prediction 2: the projected boundary is gamma=0.2 (for P=I, alpha=0),
 78    # independent of gain; the scheduler therefore cannot move it.
 79    gammas = np.linspace(-0.5, 0.5, 101)
 80    boundary_records = []
 81    for gamma in gammas:
 82        J = np.diag([-0.5, gamma])
 83        _, rows = projected_sweep(J, C, G0, P)
 84        boundary_records.append((gamma, rows[0]["lambda_max_sym"]))
 85    vals = np.array([v for _, v in boundary_records])
 86    crossing = gammas[np.where(np.sign(vals[:-1]) != np.sign(vals[1:]))[0][0]]
 87
 88    # Prediction 3: for projected-stable dynamics, scheduler always chooses rho=1,
 89    # because every candidate has the same certified margin.
 90    _, stable_rows2 = projected_sweep(J_stable, C, G0, P)
 91    selected = scheduler(stable_rows2, alpha=0.05, P=P)
 92    all_selected = []
 93    for gamma in [-0.5, -0.2, -0.05, 0.05, 0.5]:
 94        J = np.diag([-0.5, gamma])
 95        _, rows = projected_sweep(J, C, G0, P)
 96        all_selected.append({"gamma": gamma, "selected_rho": scheduler(rows, alpha=0.0, P=P)})
 97
 98    # Comparison: fixed low/high and stated scheduler on the same stable toy system.
 99    def fixed(rho):
100        return lambda z, J, C_, G: rho
101    def scheduled(z, J, C_, G):
102        _, rows = projected_sweep(J, C_, G, P)
103        return scheduler(rows, alpha=0.05, P=P)
104    comparison = {
105        "fixed_low_rho_1": noisy_observer(J_stable, C, G0, fixed(1.0), seed=SEED),
106        "fixed_high_rho_4": noisy_observer(J_stable, C, G0, fixed(4.0), seed=SEED),
107        "scheduled": noisy_observer(J_stable, C, G0, scheduled, seed=SEED),
108    }
109
110    # Also show that an unstable projected mode cannot be certified at any rho.
111    J_unstable = np.diag([-0.5, 0.10])
112    _, unstable_rows = projected_sweep(J_unstable, C, G0, P)
113    unstable_pick = scheduler(unstable_rows, alpha=0.05, P=P)
114
115    result = {
116        "seed": SEED,
117        "U_kernel": U.tolist(),
118        "prediction_1_invariance_max_abs": invariant_error,
119        "prediction_1_expected": 0.0,
120        "prediction_2_predicted_boundary_gamma": 0.0,
121        "prediction_2_observed_grid_crossing_left_endpoint": float(crossing),
122        "prediction_2_boundary_error_grid": float(abs(crossing - 0.0)),
123        "prediction_3_stable_selected_rho": selected,
124        "prediction_3_expected_rho": 1.0,
125        "prediction_3_sweep": all_selected,
126        "unstable_projected_mode_scheduler_fallback_rho": unstable_pick,
127        "comparison": comparison,
128        "interpretation": "C U_K=0 makes the proposed projected gain term identically zero; gain scheduling cannot alter latent projected stability.",
129    }
130    Path("results.json").write_text(json.dumps(result, indent=2))
131    print(json.dumps(result, indent=2))
132
133
134if __name__ == "__main__":
135    main()