1"""Small noisy snapshot comparison: unfiltered versus forward-compatible EDMD."""
 2import json
 3import numpy as np
 4from forward_intersection import make_K, principal_filter
 5
 6
 7def one_trial(seed, T=80, noise=0.04, tau=0.20, horizon=10):
 8    rng = np.random.default_rng(seed)
 9    Q0 = np.eye(5)[:, :4]
10    K = make_K(.60)
11    # Snapshot inputs lie in the dictionary; outputs are noisy ambient observations.
12    X = rng.normal(size=(4, T))
13    Y = K @ (Q0 @ X) + noise * rng.normal(size=(5, T))
14    Khat = Y @ np.linalg.pinv(Q0 @ X)
15    Q1, _ = principal_filter(Q0, Khat @ Q0, tau)
16    A0 = Q0.T @ Khat @ Q0
17    Af = Q1.T @ Khat @ Q1
18
19    # Held-out initial conditions contain only represented stable modes.
20    Z = rng.normal(size=(3, 30))
21    truth_diag = np.array([.9, .7, .8])
22    err0, err1 = [], []
23    for h in range(1, horizon + 1):
24        target = truth_diag ** h @ Z
25        pred0 = (A0[:3, :3] ** 0) @ Z  # overwritten below; explicit matrix power follows
26        pred0 = np.linalg.matrix_power(A0[:3, :3], h) @ Z
27        # The selected Q1 is normally the first three coordinate directions; use
28        # its ambient prediction and compare only to the true supported trajectory.
29        pred1 = Q1 @ np.linalg.matrix_power(Af, h) @ (Q1.T @ (Q0[:,:3] @ Z))
30        err0.append(np.mean((pred0 - target) ** 2))
31        err1.append(np.mean((pred1[:3] - target) ** 2))
32    return {
33        "dim": int(Q1.shape[1]),
34        "spectral_radius_baseline": float(max(abs(np.linalg.eigvals(A0)))),
35        "spectral_radius_idea": float(max(abs(np.linalg.eigvals(Af)))),
36        "rollout_mse_10_baseline": float(err0[-1]),
37        "rollout_mse_10_idea": float(err1[-1]),
38    }
39
40
41def run():
42    results = [one_trial(s) for s in range(10)]
43    keys = ["spectral_radius_baseline", "spectral_radius_idea",
44            "rollout_mse_10_baseline", "rollout_mse_10_idea"]
45    means = {k: float(np.mean([r[k] for r in results])) for k in keys}
46    return {"settings": {"trials": 10, "T": 80, "noise_std": .04,
47                          "tau": .20, "horizon": 10},
48            "mean": means,
49            "dimensions": [r["dim"] for r in results], "trials": results}
50
51
52if __name__ == "__main__":
53    print(json.dumps(run(), indent=2))