1import json
 2from pathlib import Path
 3import numpy as np
 4
 5
 6def phi(delta, eps):
 7    if delta == 0 and eps == 0:
 8        return 1.0
 9    return 2.0 * eps / (delta + np.sqrt(delta * delta + 4.0 * eps * eps))
10
11
12def eig_dist(A, At, norm='spectral'):
13    la = np.linalg.eigvalsh(A)[::-1]
14    lt = np.linalg.eigvalsh(At)[::-1]
15    d = np.abs(la - lt)
16    if norm in ('spectral', 2, np.inf):
17        return float(np.max(d))
18    if norm in ('fro', 'F',  'euclidean'):
19        return float(np.linalg.norm(d))
20    raise ValueError(norm)
21
22
23def two_by_two(gap, eps):
24    A = np.array([[gap, eps], [eps, 0.0]])
25    At = np.diag([gap, 0.0])
26    shift = eig_dist(A, At, 'spectral')
27    exact = (np.sqrt(gap * gap + 4 * eps * eps) - gap) / 2
28    cert = phi(gap, eps) * eps
29    return shift, exact, cert
30
31
32def rank_one_check(rng, n=5, m=4, trials=100):
33    worst = 0.0
34    violations = 0
35    for _ in range(trials):
36        h1 = np.diag(np.sort(rng.uniform(2, 8, n)))
37        h2 = np.diag(np.sort(rng.uniform(-4, 0, m)))
38        u = rng.normal(size=n); u /= np.linalg.norm(u)
39        v = rng.normal(size=m); v /= np.linalg.norm(v)
40        eps = rng.uniform(0.01, 2.0)
41        E = eps * np.outer(v, u)
42        A = np.block([[h1, E.T], [E, h2]])
43        At = np.block([[h1, np.zeros((n, m))], [np.zeros((m, n)), h2]])
44        eta = np.min(np.abs(h1.diagonal()[:, None] - h2.diagonal()[None, :]))
45        V = np.block([[np.zeros((n, n)), E.T], [E, np.zeros((m, m))]])
46        bound = phi(eta, eps) * np.linalg.norm(V, 'fro')
47        actual = eig_dist(A, At, 'fro')
48        ratio = actual / bound if bound > 0 else 0
49        worst = max(worst, ratio)
50        violations += int(actual > bound * (1 + 1e-10))
51    return {'trials': trials, 'max_actual_over_bound': worst, 'violations': violations}
52
53
54def adaptive_merge_demo():
55    gap = 1.0
56    tau = 0.08
57    rows = []
58    for eps in [0.01, 0.03, 0.1, 0.3, 1.0]:
59        _, _, c = two_by_two(gap, eps)
60        cert_f = phi(gap, eps) * np.sqrt(2.0) * eps
61        rows.append({'eps': eps, 'certificate_fro': cert_f,
62                     'decision': 'retain separate blocks' if cert_f <= tau else 'merge blocks'})
63    return {'gap': gap, 'threshold': tau, 'rows': rows}
64
65
66def main():
67    rng = np.random.default_rng(2802)
68    weak = []
69    for gap in [0.5, 1.0, 2.0, 4.0]:
70        eps = gap * 1e-5
71        actual, _, cert = two_by_two(gap, eps)
72        weak.append({'gap': gap, 'eps': eps, 'observed_shift_over_eps2': actual/(eps*eps),
73                     'predicted_limit_1_over_gap': 1/gap, 'certificate_ratio': actual/cert})
74    strong = []
75    for eps in [1, 3, 10, 30, 100]:
76        actual, _, cert = two_by_two(1.0, eps)
77        strong.append({'eps': eps, 'observed_shift_over_eps': actual/eps,
78                       'predicted_limit': 1.0, 'certificate_ratio': actual/cert})
79    zero_gap = []
80    for eps in [0.01, 0.1, 1.0, 10.0]:
81        actual, _, cert = two_by_two(0.0, eps)
82        zero_gap.append({'eps': eps, 'observed_shift': actual, 'predicted_eps': eps,
83                         'certificate_ratio': actual/cert})
84    result = {'weak_coupling_prediction': weak, 'strong_coupling_prediction': strong,
85              'zero_gap_prediction': zero_gap,
86              'rank_one_frobenius_check': rank_one_check(rng),
87              'adaptive_merge_demo': adaptive_merge_demo()}
88    Path('results.json').write_text(json.dumps(result, indent=2))
89    print(json.dumps(result, indent=2))
90
91
92if __name__ == '__main__':
93    main()