1import json
 2from pathlib import Path
 3import numpy as np
 4from carrier_probe import H, DT, Q_TRUE, carrier, gramians, rollout, fit_q
 5
 6
 7def main():
 8    rng = np.random.default_rng(123)
 9    a = 0.16
10    delta = 1e-8
11    alpha = 0.1
12    candidates = []
13    for w in [0.5, 1.0, 2.0, 4.0, 8.0]:
14        candidates.append((f"sin_w{w}", carrier("sin", w)))
15    t = np.arange(H) * DT
16    candidates.append(("binary_w2", carrier("binary", 2.0)))
17    # A short orthogonal-ish pulse, normalized to unit magnitude.
18    pulse = np.zeros(H); pulse[H//3: H//3 + H//10] = 1.0
19    candidates.append(("short_pulse", pulse))
20    rows = []
21    for name, m in candidates:
22        Wo, Wr, eo, er = gramians(a, m)
23        score = np.linalg.slogdet(Wo + delta*np.eye(2))[1] + alpha*np.linalg.slogdet(Wr + delta*np.eye(2))[1]
24        rows.append({"candidate": name, "score": float(score),
25                     "lambda_min_Wo": float(eo[0]), "lambda_min_Wr": float(er[0])})
26    rows.sort(key=lambda x: x["score"], reverse=True)
27
28    # Recovery curve under identical observation noise scale.
29    recovery = []
30    for aa in [0.0, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32]:
31        m = carrier("sin", 2.0)
32        obs = rollout(aa, m)[:, 0] + rng.normal(0, 0.01, H+1)
33        qhat = fit_q(obs, m, aa)
34        recovery.append({"a": aa, "q_hat": qhat, "abs_error": abs(qhat-Q_TRUE)})
35    out = {"amplitude": a, "alpha": alpha, "candidates_by_score": rows,
36           "recovery_curve": recovery,
37           "selected": rows[0]["candidate"]}
38    Path("candidate_results.json").write_text(json.dumps(out, indent=2))
39    print(json.dumps(out, indent=2))
40
41if __name__ == "__main__":
42    main()