1import json
 2import numpy as np
 3
 4
 5def spectral_gate(A, c, tau=1e-3):
 6    c = np.asarray(c, float); A = np.asarray(A, float)
 7    lam = float(c @ A @ c / (c @ c))
 8    residual = float(np.linalg.norm(c @ A - lam * c) / np.linalg.norm(c))
 9    return residual, lam, (residual > tau or lam < 0)
10
11
12def powers_and_cov(A, P0, Q, H):
13    """Return propagated covariances P_j, including P_0."""
14    P = np.asarray(P0, float).copy(); out = [P.copy()]
15    for _ in range(H):
16        P = A @ P @ A.T + Q
17        out.append(P.copy())
18    return out
19
20
21def verify_math(rng):
22    d = 4
23    # An exact left eigenvector, and a non-eigenvector rotating readout.
24    A = np.diag([0.82, 0.61, 0.45, 0.30])
25    c = np.array([1., 0., 0., 0.])
26    lam = .82
27    max_power_err = 0.
28    for j in range(7):
29        max_power_err = max(max_power_err, np.linalg.norm(c @ np.linalg.matrix_power(A, j) - lam**j*c))
30    r, l, skip = spectral_gate(A, c)
31    # Covariance recursion versus the closed finite sum.
32    P0 = np.diag([.2, .1, .15, .05]); Q = .01*np.eye(d); H = 6
33    rec = powers_and_cov(A, P0, Q, H)
34    closed = []
35    for j in range(H+1):
36        Aj = np.linalg.matrix_power(A, j)
37        Pj = Aj @ P0 @ Aj.T
38        for i in range(j):
39            Ai = np.linalg.matrix_power(A, i)
40            Pj += Ai @ Q @ Ai.T
41        closed.append(Pj)
42    cov_err = max(np.linalg.norm(rec[j]-closed[j]) for j in range(H+1))
43    # A rotation must have a clearly nonzero residual and activate the gate.
44    th = .55; R = .98*np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
45    rr, ll, active = spectral_gate(R, np.array([1., 0.]), tau=1e-3)
46    return dict(power_identity_error=max_power_err, covariance_identity_error=cov_err,
47                aligned_residual=r, aligned_lambda=l, aligned_skips=not skip,
48                rotating_residual=rr, rotating_lambda=ll, rotating_activates=active)
49
50
51def rollout_decision(A, x, c, delta, H):
52    """Alarm if current or any of H predicted readouts crosses delta."""
53    y = np.asarray(x, float).copy(); best = float(c @ y)
54    for _ in range(H):
55        y = A @ y
56        best = max(best, float(c @ y))
57    return best >= delta
58
59
60def benchmark(rng, n=3000, H=8, delta=.8, tau=0.02):
61    th = .55
62    rotation = .98*np.array([[np.cos(th), -np.sin(th)], [np.sin(th), np.cos(th)]])
63    aligned = np.diag([.82, .60])
64    c = np.array([1., 0.])
65    systems = [("aligned_positive", aligned), ("rotating", rotation)]
66    rows = []
67    for name, A in systems:
68        # Random current states. Deterministic dynamics make the structural result visible.
69        X = rng.normal(size=(n, 2))
70        truth = np.array([rollout_decision(A, x, c, delta, H) for x in X])
71        current = (X @ c >= delta)
72        always = np.array([rollout_decision(A, x, c, delta, H) for x in X])
73        residual, lam, activates = spectral_gate(A, c, tau)
74        gated = current if not activates else always
75        # Cost is number of matrix-vector rollout steps, with current readout free in this comparison.
76        gate_steps = (H if activates else 0) * n
77        always_steps = H*n
78        def scores(pred):
79            tp = np.sum(pred & truth); fp = np.sum(pred & ~truth); fn = np.sum(~pred & truth)
80            return dict(recall=float(tp/max(1,tp+fn)), false_positive_rate=float(fp/max(1,np.sum(~truth))))
81        rows.append(dict(system=name, residual=residual, lambda_hat=lam,
82                         gate_activates=activates, truth_rate=float(truth.mean()),
83                         current=scores(current), always_rollout=scores(always), gated=scores(gated),
84                         always_rollout_steps=always_steps, gated_rollout_steps=gate_steps,
85                         rollout_step_reduction=float(1-gate_steps/always_steps)))
86    # Mixed workload: equal systems, showing selective compute and overall recall.
87    return rows
88
89
90def main():
91    rng = np.random.default_rng(12345)
92    math = verify_math(rng)
93    rows = benchmark(rng)
94    result = {"math_verification": math, "benchmark": rows}
95    with open("results.json", "w") as f: json.dump(result, f, indent=2)
96    print(json.dumps(result, indent=2))
97
98if __name__ == '__main__':
99    main()