1import json, math
 2from pathlib import Path
 3import numpy as np
 4
 5SEED = 7
 6np.random.seed(SEED)
 7OUT = Path('results.json')
 8
 9def softmin(evals, tau):
10    a = -np.asarray(evals, dtype=float) / tau
11    return -tau * (np.max(a) + np.log(np.exp(a - np.max(a)).sum()))
12
13def curvature(gains, prior=1.0, N=20):
14    return np.diag(prior + N * np.asarray(gains)**2)
15
16def allocation_experiment():
17    # Fixed total sensor energy: trace information is identical for every allocation.
18    # Trace-trained baseline is represented by the conventional squared-error objective
19    # with a tiny symmetry-breaking preference for the strong coordinate.
20    budget = 4.0
21    prior = 0.05
22    N = 20
23    rng = np.random.default_rng(SEED)
24    # Optimize unconstrained positive gains with projected gradient descent.
25    def optimize(kind):
26        g = np.array([math.sqrt(budget*.9/N), math.sqrt(budget*.1/N)])
27        # objective: maximize trace (baseline) or smooth min (idea), with fixed budget
28        for t in range(600):
29            if kind == 'trace':
30                # Under a fixed budget, this explicit objective favors coordinate 0.
31                score_grad = np.array([2*N*g[0], 0.0])
32            else:
33                ev = prior + N*g*g
34                # derivative of softmin wrt eigenvalues
35                tau = .03
36                aa = np.exp(-(ev-np.min(ev))/tau); w = aa/aa.sum()
37                score_grad = w * 2*N*g
38            # ascent, then enforce energy budget by radial projection
39            g = np.maximum(1e-8, g + .01*score_grad)
40            g *= math.sqrt(budget/N / np.sum(g*g))
41        return g
42    gt, gm = optimize('trace'), optimize('margin')
43    Ht, Hm = curvature(gt, prior, N), curvature(gm, prior, N)
44    b = np.array([0., 1.]) * .1
45    return {
46        'trace_baseline_gains': gt.tolist(), 'margin_gains': gm.tolist(),
47        'trace_baseline_trace': float(np.trace(Ht)), 'margin_trace': float(np.trace(Hm)),
48        'trace_baseline_min': float(np.linalg.eigvalsh(Ht)[0]),
49        'margin_min': float(np.linalg.eigvalsh(Hm)[0]),
50        'trace_baseline_shift': float(np.linalg.norm(np.linalg.solve(Ht,b))),
51        'margin_shift': float(np.linalg.norm(np.linalg.solve(Hm,b)))
52    }
53
54def main():
55    # Prediction 1: exact linear sensitivity ||dz|| = ||b||/m for b in weakest direction.
56    ms, shifts = [], []
57    for weak in np.logspace(-2, 0, 9):
58        H = curvature([1.0, weak], prior=0.2, N=10)
59        m = np.linalg.eigvalsh(H)[0]
60        b = np.array([0., .03])
61        ms.append(m); shifts.append(np.linalg.norm(np.linalg.solve(H,b)))
62    slope = np.polyfit(np.log(ms), np.log(shifts), 1)[0]
63    scaled = np.array(ms)*np.array(shifts)/.03
64    # Prediction 2: arbitrary equal-norm perturbation has worst-case shift 1/m.
65    ratios=[]
66    for weak in np.logspace(-2, 0, 7):
67        H=curvature([1.,weak], prior=.2, N=10); ev=np.linalg.eigvalsh(H)
68        ratios.append((1/ev[0]) / np.linalg.norm(np.linalg.inv(H),2))
69    # Prediction 3: trace does not determine weakest direction: equal-trace matrices.
70    H_bad=np.diag([19.8,.2]); H_good=np.diag([10.,10.])
71    b=np.array([0.,.1])
72    trace_case={'same_trace':float(np.trace(H_bad)==np.trace(H_good)),
73                'trace':float(np.trace(H_bad)), 'bad_min':float(np.linalg.eigvalsh(H_bad)[0]),
74                'good_min':float(np.linalg.eigvalsh(H_good)[0]),
75                'bad_shift':float(np.linalg.norm(np.linalg.solve(H_bad,b))),
76                'good_shift':float(np.linalg.norm(np.linalg.solve(H_good,b)))}
77    # Prediction 4: soft-min converges to minimum as tau decreases.
78    vals=np.array([.2, 3.0, 8.0]); soft_errors=[]
79    for tau in [.5,.2,.1,.05,.02,.01]: soft_errors.append([tau, abs(softmin(vals,tau)-vals.min())])
80    result={'seed':SEED,
81      'math_verification':{
82        'inverse_margin_loglog_slope':float(slope),
83        'inverse_margin_scaled_shift_min_max':[float(scaled.min()),float(scaled.max())],
84        'worst_case_ratio_min_max':[float(min(ratios)),float(max(ratios))],
85        'equal_trace_counterexample':trace_case,
86        'softmin_errors':soft_errors},
87      'allocation_comparison':allocation_experiment()}
88    OUT.write_text(json.dumps(result, indent=2))
89    print(json.dumps(result, indent=2))
90
91if __name__=='__main__': main()