1import json, random
 2from pathlib import Path
 3import numpy as np
 4import torch
 5
 6SEED = 2291
 7np.random.seed(SEED); random.seed(SEED); torch.manual_seed(SEED)
 8torch.set_num_threads(4)
 9
10def laplacian(a, b):
11    # Directed path 1 -> 2 -> 3, weighted in-degree Laplacian.
12    return np.array([[0., 0., 0.], [-a, a, 0.], [0., -b, b]])
13
14def spectrum(a, b):
15    return np.linalg.eig(laplacian(a, b))
16
17def gap(a, b):
18    ev, _ = spectrum(a, b)
19    nz = [z for z in ev if abs(z) > 1e-7]
20    return float(min(abs(nz[i]-nz[j])/(1+abs(nz[i])+abs(nz[j]))
21                    for i in range(len(nz)) for j in range(i)))
22
23def cond(a, b):
24    return float(np.linalg.cond(spectrum(a, b)[1]))
25
26def discriminant_nonzero(a, b):
27    # q(s)=(s-a)(s-b), so Disc(q)=(a-b)^2.
28    return float((a-b)**2)
29
30def analytic_gap(a,b):
31    return abs(a-b)/(1+a+b)
32
33def math_sweep():
34    rows=[]; s=1.0
35    for eps in [1e-4,3e-4,1e-3,3e-3,1e-2,3e-2,1e-1]:
36        a,b=s+eps,s-eps
37        rows.append({'eps':eps,'gap':gap(a,b),'analytic_gap':analytic_gap(a,b),
38                     'discriminant':discriminant_nonzero(a,b),'condition':cond(a,b)})
39    x=np.log([r['eps'] for r in rows[:5]])
40    gs=float(np.polyfit(x,np.log([r['gap'] for r in rows[:5]]),1)[0])
41    ds=float(np.polyfit(x,np.log([r['discriminant'] for r in rows[:5]]),1)[0])
42    cp=[r['condition']*r['eps'] for r in rows[:5]]
43    L=laplacian(1.,1.)
44    ev=np.linalg.eigvals(L)
45    # At a=b, algebraic multiplicity of eigenvalue 1 is two but nullity of
46    # L-I is one: this is the exact defective boundary.
47    nullity=int(3-np.linalg.matrix_rank(L-np.eye(3),tol=1e-10))
48    return {'rows':rows,'observed_gap_log_slope':gs,
49            'observed_discriminant_log_slope':ds,'condition_times_epsilon':cp,
50            'exact_collision':{'eigenvalues':ev.tolist(),
51                'disc_q':discriminant_nonzero(1.,1.),
52                'mult_1_algebraic':2,'nullity_L_minus_I':nullity,
53                'defective':bool(nullity<2)},
54            'predictions':{'gap_slope':1.,'discriminant_slope':2.,
55                           'condition_scaling':'condition*epsilon approximately constant'}}
56
57def train(gated, steps=700, gmin=.03, beta=5000.):
58    torch.manual_seed(SEED)
59    theta=torch.tensor([.4,-.4],dtype=torch.float64,requires_grad=True)
60    opt=torch.optim.Adam([theta],lr=.025); omin=1e-3; hist=[]
61    for k in range(steps):
62        a,b=omin+torch.nn.functional.softplus(theta)
63        task=(a-1.)**2+(b-1.)**2
64        gg=torch.abs(a-b)/(1.+a+b)
65        barrier=beta*torch.relu(torch.tensor(gmin,dtype=torch.float64)-gg)**2 if gated else 0.
66        loss=task+barrier
67        opt.zero_grad(); loss.backward(); opt.step()
68        if k in (0,99,299,699):
69            hist.append({'step':k+1,'a':float(a.detach()),'b':float(b.detach()),
70                         'task':float(task.detach()),'gap':float(gg.detach()),
71                         'barrier':float(barrier.detach()) if gated else 0.})
72    af,bf=[float(x.detach()) for x in (a,b)]
73    return {'a':af,'b':bf,'gap':gap(af,bf),'condition':cond(af,bf),
74            'task':(af-1)**2+(bf-1)**2,'history':hist}
75
76def main():
77    math_result=math_sweep()
78    base=train(False)
79    sweep={str(g):train(True,gmin=g) for g in [1e-4,1e-3,1e-2,3e-2]}
80    out={'math':math_result,'training':{'baseline_unconstrained_positive':base,
81          'gap_gated_gmin_sweep':sweep,'beta':5000.},'seed':SEED}
82    Path('results.json').write_text(json.dumps(out,indent=2))
83    print(json.dumps(out,indent=2))
84if __name__=='__main__': main()