1import itertools, json
 2import numpy as np
 3
 4SEED = 7
 5
 6
 7def subset_stats(theta, m):
 8    theta = np.asarray(theta, dtype=float)
 9    d = len(theta)
10    subs = np.array(list(itertools.combinations(range(d), m)), dtype=int)
11    scores = theta[subs].sum(axis=1)
12    scores -= scores.max()
13    p = np.exp(scores); p /= p.sum()
14    X = np.zeros((len(subs), d))
15    X[np.arange(len(subs))[:, None], subs] = 1.0
16    mu = p @ X
17    cov = (X * p[:, None]).T @ X - np.outer(mu, mu)
18    return mu, cov, subs, p
19
20
21def resistance(cov, i, j):
22    pinv = np.linalg.pinv(cov, rcond=1e-11)
23    z = np.zeros(cov.shape[0]); z[i] = 1; z[j] = -1
24    return float(z @ pinv @ z)
25
26
27def bound_sweep():
28    rows=[]
29    for d,m in [(6,2),(8,4),(12,2),(12,4),(12,6)]:
30        worst=0.; worst_pair=None; nullerr=0.; min_psd=1e9
31        for a in np.linspace(0, 12, 13):
32            theta=np.zeros(d); theta[0]=a; theta[1]=-a/2
33            mu,cov,_,_=subset_stats(theta,m)
34            nullerr=max(nullerr, float(np.max(np.abs(cov@np.ones(d)))))
35            v=np.diag(cov); V=v.sum(); B=.5*(np.diag(v)-np.outer(v,v)/V)
36            min_psd=min(min_psd, float(np.linalg.eigvalsh(cov-B).min()))
37            for i in range(d):
38                for j in range(i+1,d):
39                    r=resistance(cov,i,j)/(1/v[i]+1/v[j])
40                    if r>worst: worst=float(r); worst_pair=(a,i,j)
41        rows.append(dict(d=d,m=m,max_resistance_ratio=worst,worst_case=worst_pair,
42                         max_null_residual=nullerr,min_cov_minus_bound_eigenvalue=min_psd))
43    return rows
44
45
46def trust_scaling():
47    d,m=12,4; theta=np.linspace(2,-2,d); mu,cov,_,_=subset_stats(theta,m)
48    g=np.linspace(-1,1,d); g-=g.mean()
49    u=np.linalg.pinv(cov, rcond=1e-11)@g; u-=u.mean()
50    v=np.diag(cov)
51    raw=max(abs(u[i]-u[j])/np.sqrt(1/v[i]+1/v[j]) for i in range(d) for j in range(i+1,d))
52    out=[]
53    for rho in [0.01,0.03,0.1,0.3,1.0]:
54        scale=min(1.,rho/raw); delta=-u*scale
55        observed=max(abs(delta[i]-delta[j])/np.sqrt(1/v[i]+1/v[j]) for i in range(d) for j in range(i+1,d))
56        out.append(dict(rho=rho,raw_normalized_step=raw,observed=observed,predicted=min(rho,raw)))
57    return out
58
59
60def optimize_compare():
61    d,m=12,4
62    target,_c,_,_=subset_stats(np.array([3.,2.,1.,.5,0,0,0,0,0,0,-1.,-2.]),m)
63    results=[]
64    for method in ['vanilla','natural_trust']:
65        theta=np.zeros(d); losses=[]; cvs=[]
66        for t in range(80):
67            mu,cov,_,_=subset_stats(theta,m)
68            err=mu-target; loss=.5*float(err@err)
69            losses.append(loss); cvs.append(float(np.std(mu)/np.mean(mu)))
70            g=cov@err
71            if method=='vanilla':
72                delta=-1.4*g
73            else:
74                u=np.linalg.pinv(cov,rcond=1e-11)@g; u-=u.mean()
75                v=np.diag(cov); rho=.35
76                mx=max(abs(u[i]-u[j])/np.sqrt(1/v[i]+1/v[j]) for i in range(d) for j in range(i+1,d))
77                delta=-1.4*u*min(1.,rho/mx) if mx>0 else np.zeros(d)
78            theta += delta
79        results.append(dict(method=method,initial_loss=losses[0],final_loss=losses[-1],
80                            loss_step_10=losses[10],loss_step_40=losses[40],
81                            final_load_cv=cvs[-1],min_loss=float(min(losses))))
82    return results
83
84
85def main():
86    report={'seed':SEED,
87      'predictions':{
88        'resistance_bound':'max normalized resistance ratio <= 1 for every field and pair',
89        'covariance_lower_bound':'Sigma - 0.5(D-vv^T/V) is PSD',
90        'trust_scaling':'normalized max pairwise update equals min(rho, raw normalized step)'},
91      'bound_sweep':bound_sweep(),'trust_sweep':trust_scaling(),
92      'optimization':optimize_compare()}
93    print(json.dumps(report,indent=2))
94
95if __name__=='__main__': main()