1import json, math, time
  2from pathlib import Path
  3import numpy as np
  4
  5
  6def wrap(x):
  7    return (x + np.pi) % (2*np.pi) - np.pi
  8
  9
 10def spectral_init(n, edges, phi, weights, ref=0, iters=30):
 11    C = np.zeros((n, n), dtype=np.complex128)
 12    for (i, j), p, w in zip(edges, phi, weights):
 13        C[i, j] += w * np.exp(1j*p)
 14        C[j, i] += w * np.exp(-1j*p)
 15    # power iteration, with Rayleigh-sign stabilization
 16    v = np.ones(n, dtype=np.complex128) / np.sqrt(n)
 17    for _ in range(iters):
 18        q = C @ v
 19        nq = np.linalg.norm(q)
 20        if nq == 0: break
 21        v = q / nq
 22    z = v / np.maximum(np.abs(v), 1e-10)
 23    theta = np.angle(z)
 24    return wrap(theta - theta[ref])
 25
 26
 27def gauge_rmse(est, truth):
 28    d = wrap(est - truth)
 29    # both are nominally gauged, but optimize residual global gauge for robustness
 30    shift = np.angle(np.mean(np.exp(1j*(est-truth))))
 31    return float(np.sqrt(np.mean(wrap(d-shift)**2)))
 32
 33
 34def edge_rmse(theta, truth, edges):
 35    return float(np.sqrt(np.mean([wrap((theta[i]-theta[j])-(truth[i]-truth[j]))**2 for i,j in edges])))
 36
 37
 38def centered_jacobian_margin(n, edges, weights):
 39    """Smallest nonzero singular value after removing the global phase gauge."""
 40    J = np.zeros((len(edges), n))
 41    for k, (i, j) in enumerate(edges):
 42        J[k, i] = weights[k]
 43        J[k, j] = -weights[k]
 44    P = np.eye(n) - np.ones((n, n)) / n
 45    sv = np.linalg.svd(J @ P, compute_uv=False)
 46    nz = sv[sv > 1e-8]
 47    return float(nz[-1]) if nz.size else 0.0
 48
 49
 50def gated_refine(init, edges, phi, weights, rho, eta, threshold=0.35):
 51    """Use spectral estimate directly in the low-eta regime; otherwise refine."""
 52    if eta <= threshold:
 53        return init.copy(), False
 54    return gd_refine(init, edges, phi, weights, lr=0.04, steps=100), True
 55
 56
 57def gd_refine(init, edges, phi, weights, lr=0.08, steps=100):
 58    # WLS on circular relative-phase residuals; node 0 fixed as the gauge.
 59    x = init.copy(); x[0] = 0.0
 60    for _ in range(steps):
 61        g = np.zeros_like(x)
 62        for (i,j), p, w in zip(edges, phi, weights):
 63            r = wrap(x[i]-x[j]-p)
 64            # derivative of 1-cos(r), robust and smooth
 65            s = w*np.sin(r)
 66            g[i] += s; g[j] -= s
 67        x[1:] -= lr*g[1:]
 68        x[0] = 0.0
 69    return wrap(x)
 70
 71
 72def trial(seed, n=20, degree=3, noise=0.35):
 73    rng=np.random.default_rng(seed)
 74    truth=wrap(rng.normal(0, 1.0, n)); truth -= truth[0]
 75    # connected ring plus random undirected edges
 76    edges={(i,(i+1)%n) for i in range(n)}
 77    target=n*degree//2
 78    while len(edges)<target:
 79        i,j=sorted(rng.choice(n,2,replace=False)); edges.add((i,j))
 80    edges=sorted(edges)
 81    weights=np.ones(len(edges))
 82    phi=np.array([wrap(truth[i]-truth[j]+rng.normal(0,noise)) for i,j in edges])
 83    spec=spectral_init(n,edges,phi,weights)
 84    rand=wrap(rng.uniform(-np.pi,np.pi,n)); rand -= rand[0]
 85    spec_loss=edge_rmse(spec,truth,edges); rand_loss=edge_rmse(rand,truth,edges)
 86    spec_final=gd_refine(spec,edges,phi,weights)
 87    rand_final=gd_refine(rand,edges,phi,weights)
 88    rho = centered_jacobian_margin(n, edges, weights)
 89    residual=np.array([wrap(spec[i]-spec[j]-p) for (i,j),p in zip(edges,phi)])
 90    mad=float(np.median(np.abs(residual-np.median(residual))))
 91    eta=mad/rho
 92    gated, did_refine = gated_refine(spec, edges, phi, weights, rho, eta)
 93
 94    # robust residual and normalized perturbation proxy
 95    return dict(rho=rho, eta=eta, gate_refined=did_refine, gated_rmse=gauge_rmse(gated,truth), spec_init_rmse=gauge_rmse(spec,truth), rand_init_rmse=gauge_rmse(rand,truth),
 96                spec_edge_rmse=spec_loss, rand_edge_rmse=rand_loss,
 97                spec_final_rmse=gauge_rmse(spec_final,truth), rand_final_rmse=gauge_rmse(rand_final,truth),
 98                spec_final_edge_rmse=edge_rmse(spec_final,truth,edges), rand_final_edge_rmse=edge_rmse(rand_final,truth,edges),
 99                mad=mad)
100
101
102def jacobian_check(n=8):
103    # Each edge observes theta_i-theta_j. J has one row per edge.
104    edges=[(i,(i+1)%n) for i in range(n)] + [(i,(i+2)%n) for i in range(n)]
105    J=np.zeros((len(edges),n))
106    for k,(i,j) in enumerate(edges): J[k,i]=1; J[k,j]=-1
107    one=np.ones(n)
108    gauge_norm=float(np.linalg.norm(J@one))
109    P=np.eye(n)-np.ones((n,n))/n
110    sv=np.linalg.svd(J@P,compute_uv=False)
111    rho=float(sv[-2] if sv.size>1 and sv[-1]<1e-8 else sv[-1])
112    # constrained singular values should equal the nonzero singular values of J
113    nonzero=np.linalg.svd(J,compute_uv=False); nonzero=nonzero[nonzero>1e-8]
114    return dict(gauge_null_norm=gauge_norm, rho=rho, nonzero_min=float(nonzero[-1]),
115                rho_matches=float(abs(rho-nonzero[-1])))
116
117
118def main():
119    t=time.time(); jac=jacobian_check()
120    allrows=[]
121    for noise in (0.15,0.35,0.70,1.10):
122        rows=[trial(1000+q,noise=noise) for q in range(40)]
123        summary={'noise':noise}
124        for key in rows[0]: summary[key+'_mean']=float(np.mean([r[key] for r in rows]))
125        summary['spec_better_init_rate']=float(np.mean([r['spec_init_rmse']<r['rand_init_rmse'] for r in rows]))
126        summary['spec_better_final_rate']=float(np.mean([r['spec_final_rmse']<r['rand_final_rmse'] for r in rows]))
127        allrows.append(summary)
128    out={'jacobian_check':jac,'results':allrows,'runtime_sec':time.time()-t}
129    Path('results.json').write_text(json.dumps(out,indent=2))
130    print(json.dumps(out,indent=2))
131
132if __name__=='__main__': main()