1import json
 2from pathlib import Path
 3import numpy as np
 4from scipy.linalg import expm
 5
 6SEED = 2273
 7
 8
 9def oscillator_matrix(damping=0.08, omega=1.7):
10    return np.array([[damping, -omega], [omega, damping]], dtype=float)
11
12
13def delay_row(A, c, t):
14    return np.asarray(c @ expm(-A * float(t))).reshape(-1)
15
16
17def sigma_min(O):
18    return float(np.linalg.svd(O, compute_uv=False)[-1]) if O.shape[0] >= O.shape[1] else 0.0
19
20
21def greedy_delays(A, c, candidates, budget, delta=0.0, initial=(0.0,)):
22    selected = list(map(float, initial))
23    O = np.stack([delay_row(A, c, t) for t in selected])
24    margins = [sigma_min(O)]
25    gains = []
26    for _ in range(max(0, budget - len(selected))):
27        remaining = [float(t) for t in candidates if all(abs(t-u) > 1e-10 for u in selected)]
28        if not remaining:
29            break
30        scored = [(sigma_min(np.vstack((O, delay_row(A, c, t)))), t) for t in remaining]
31        best, t = max(scored, key=lambda x: x[0])
32        gain = best - margins[-1]
33        if gain < delta:
34            break
35        selected.append(t)
36        O = np.vstack((O, delay_row(A, c, t)))
37        margins.append(float(best)); gains.append(float(gain))
38    return np.array(selected), O, np.array(margins), np.array(gains)
39
40
41def uniform_delays(tmax, count):
42    return np.linspace(0.0, float(tmax), count)
43
44
45def reconstruction_mse(O, noise_std, trials=12000, seed=991):
46    rng = np.random.default_rng(seed)
47    x = rng.normal(size=(O.shape[1], trials))
48    y = O @ x + noise_std * rng.normal(size=(O.shape[0], trials))
49    xhat = np.linalg.pinv(O) @ y
50    return float(np.mean((xhat-x)**2))
51
52
53def main():
54    c = np.array([1.0, 0.0])
55    candidates = np.unique(np.r_[np.linspace(0.0, 4.0, 161), np.geomspace(0.01, 12.0, 100)])
56    rows = []
57    A = oscillator_matrix(); selected, O, margins, _ = greedy_delays(A, c, candidates, 8)
58    rows.append({'prediction':'row append monotonicity for every greedy step', 'observed_min_increment':float(np.min(np.diff(margins))), 'predicted':'>= 0', 'confirmed':bool(np.all(np.diff(margins) >= -1e-12))})
59
60    # Frequency sweep: informative quadrature delay is pi/(2 omega).
61    freq_sweep = []
62    for omega in (0.8, 1.2, 1.7, 2.4, 3.2):
63        A = oscillator_matrix(omega=omega)
64        ts, _, _, _ = greedy_delays(A, c, candidates, 2)
65        predicted = np.pi/(2*omega)
66        observed = float(ts[1])
67        freq_sweep.append({'omega':omega, 'predicted_quarter':predicted, 'observed_first_delay':observed, 'absolute_error':abs(observed-predicted)})
68    rows.append({'prediction':'first greedy delay follows quarter-period pi/(2 omega)', 'sweep':freq_sweep, 'max_absolute_error':max(x['absolute_error'] for x in freq_sweep), 'confirmed':max(x['absolute_error'] for x in freq_sweep) < 0.15})
69
70    # Noise sweep: fixed linear inverse predicts MSE proportional to noise variance.
71    A = oscillator_matrix(); _, O4, margins4, _ = greedy_delays(A, c, candidates, 4)
72    noise_sweep = []
73    for noise in (0.005, 0.01, 0.02, 0.04, 0.08):
74        mse = reconstruction_mse(O4, noise)
75        noise_sweep.append({'noise_std':noise, 'mse':mse, 'mse_over_noise2':mse/(noise*noise)})
76    ratios = np.array([x['mse_over_noise2'] for x in noise_sweep])
77    rows.append({'prediction':'reconstruction MSE scales as noise_std^2 / sigma_min^2', 'sigma_min':float(margins4[-1]), 'sweep':noise_sweep, 'ratio_cv':float(np.std(ratios)/np.mean(ratios)), 'confirmed':bool(np.std(ratios)/np.mean(ratios) < 0.04)})
78
79    # Margin sweep and inverse-margin relation across retained tap counts.
80    checks=[]
81    for budget in range(2,9):
82        _, OO, mm, _ = greedy_delays(A, c, candidates, budget)
83        mse = reconstruction_mse(OO, 0.03)
84        checks.append((float(mm[-1]), mse, float(1/(mm[-1]**2 + 1e-12))))
85    margins3=np.array([x[0] for x in checks]); mses=np.array([x[1] for x in checks]); inv=np.array([x[2] for x in checks])
86    corr=float(np.corrcoef(mses, inv)[0,1])
87    rows.append({'prediction':'error decreases with inverse stable margin squared', 'budgets':list(range(2,9)), 'sigma_min':margins3.tolist(), 'mse':mses.tolist(), 'inv_sigma2':inv.tolist(), 'correlation':corr, 'confirmed':corr > 0.85})
88
89    compare=[]
90    for budget in (4, 8, 16):
91        gts, gO, gm, _ = greedy_delays(A,c,candidates,budget)
92        uts=uniform_delays(4.0,budget); uO=np.stack([delay_row(A,c,t) for t in uts])
93        compare.append({'budget':budget, 'greedy_delays':gts.tolist(), 'uniform_delays':uts.tolist(), 'greedy_sigma_min':sigma_min(gO), 'uniform_sigma_min':sigma_min(uO), 'greedy_mse':reconstruction_mse(gO,.03), 'uniform_mse':reconstruction_mse(uO,.03)})
94    result={'seed':SEED,'system':{'damping':.08,'omega':1.7},'checks':rows,'comparison':compare}
95    Path('results.json').write_text(json.dumps(result, indent=2))
96    print(json.dumps(result, indent=2))
97
98if __name__ == '__main__': main()