1"""Gauge-covariant Wilson-loop toy verification.
  2
  3Run with: python wilson_experiment.py
  4Outputs results.json and prints quantitative predictions versus observations.
  5Quaternion q=(w,x,y,z) represents SU(2); Re Tr(U)/2 equals w.
  6"""
  7import json
  8import math
  9import numpy as np
 10
 11SEED = 2472
 12rng = np.random.default_rng(SEED)
 13
 14
 15def qmul(a, b):
 16    a, b = np.asarray(a), np.asarray(b)
 17    w = a[..., 0]*b[..., 0] - np.sum(a[..., 1:]*b[..., 1:], axis=-1)
 18    v = (a[..., 0, None]*b[..., 1:] + b[..., 0, None]*a[..., 1:]
 19         + np.cross(a[..., 1:], b[..., 1:]))
 20    return np.concatenate([w[..., None], v], axis=-1)
 21
 22
 23def qconj(a):
 24    a = np.asarray(a).copy()
 25    a[..., 1:] *= -1
 26    return a
 27
 28
 29def qnorm(a):
 30    return a / np.linalg.norm(a, axis=-1, keepdims=True)
 31
 32
 33def random_q(n=1, angle=None):
 34    if angle is None:
 35        x = rng.normal(size=(n, 4))
 36        return qnorm(x)
 37    axis = rng.normal(size=(n, 3))
 38    axis /= np.linalg.norm(axis, axis=1, keepdims=True)
 39    a = np.broadcast_to(np.asarray(angle), (n,))
 40    return np.concatenate([np.cos(a)[:, None], np.sin(a)[:, None]*axis], axis=1)
 41
 42
 43def loop(uij, ujk, uki):
 44    return qmul(qmul(uij, ujk), uki)
 45
 46
 47def compatibility(qloops):
 48    return float(np.mean(qloops[:, 0]))
 49
 50
 51def gauge_transform(uij, hi, hj):
 52    return qmul(qmul(hi, uij), qconj(hj))
 53
 54
 55def analytic_mean(theta):
 56    # For independently isotropic fixed-angle edge errors, scalar loop mean is cos(theta)^3.
 57    return math.cos(theta)**3
 58
 59
 60def math_checks():
 61    n = 2000
 62    max_err = 0.0
 63    for _ in range(n):
 64        u, v, w = random_q(3)
 65        hi, hj, hk = random_q(3)
 66        h = loop(u, v, w)
 67        hp = loop(gauge_transform(u, hi, hj),
 68                  gauge_transform(v, hj, hk),
 69                  gauge_transform(w, hk, hi))
 70        expected = qmul(qmul(hi, h), qconj(hi))
 71        max_err = max(max_err, float(np.max(np.abs(hp - expected))))
 72    # A pure gauge connection U_ij=h_i h_j^dagger has exactly identity loops.
 73    flat = []
 74    for _ in range(1000):
 75        hi, hj, hk = random_q(3)
 76        flat.append(loop(qmul(hi, qconj(hj)), qmul(hj, qconj(hk)),
 77                         qmul(hk, qconj(hi))))
 78    flat = np.asarray(flat)
 79    return {"conjugation_max_abs_error": max_err,
 80            "pure_gauge_mean_M": compatibility(flat),
 81            "pure_gauge_max_energy": float(np.max(1-flat[:, 0]))}
 82
 83
 84def disorder_sweep():
 85    # Triangles are independent here; this isolates the predicted edge-disorder law.
 86    sigmas = np.linspace(0.0, 0.9, 10)
 87    rows = []
 88    for sigma in sigmas:
 89        vals = []
 90        for _ in range(300):
 91            angles = rng.normal(0.0, sigma, size=(300, 3))
 92            # Each edge has an independent random axis, but its scalar part is cos(angle).
 93            edges = [random_q(300, angles[:, i]) for i in range(3)]
 94            vals.append(np.mean(loop(edges[0], edges[1], edges[2])[:, 0]))
 95        vals = np.asarray(vals)
 96        observed = float(np.mean(vals))
 97        # E cos(theta) = exp(-sigma^2/2), so E M = exp(-3 sigma^2/2).
 98        predicted = math.exp(-1.5*sigma*sigma)
 99        rows.append({"sigma": float(sigma), "observed_M": observed,
100                     "predicted_M": predicted,
101                     "abs_error": abs(observed-predicted),
102                     "susceptibility": float(300*np.var(vals))})
103    return rows
104
105
106def regularization_toy():
107    """Optimize scalar edge angles to fit noisy target angles.
108
109    Baseline fits each of three directed edge angles independently. Wilson adds
110    lambda*(sum(edge angles))^2, the small-angle form of 1-ReTr(H)/2.
111    This is a local coordinate chart of the SU(2) loop near identity.
112    """
113    targets = np.array([0.55, -0.40, 0.35])
114    noise = np.array([0.30, -0.15, 0.20])
115    y = targets + noise
116    lr, steps = 0.08, 300
117    out = []
118    for lam in [0.0, 0.1, 0.5, 1.0]:
119        x = np.zeros(3)
120        for _ in range(steps):
121            # task MSE + lambda*(1-cos(sum)); exact Wilson scalar in this chart
122            s = np.sum(x)
123            grad = 2*(x-y)/3.0 + lam*np.sin(s)
124            x -= lr*grad
125        task = float(np.mean((x-targets)**2))
126        wilson = float(1-math.cos(np.sum(x)))
127        out.append({"lambda": lam, "task_mse": task, "loop_energy": wilson,
128                    "learned_angles": x.tolist()})
129    return out
130
131
132def main():
133    checks = math_checks()
134    sweep = disorder_sweep()
135    opt = regularization_toy()
136    # Quantitative acceptance tests: exact covariance, flatness, and disorder scaling.
137    max_sweep_error = max(r["abs_error"] for r in sweep)
138    result = {"seed": SEED, "math_checks": checks,
139              "disorder_sweep": sweep, "regularization_toy": opt,
140              "summary": {"max_disorder_prediction_error": max_sweep_error,
141                          "disorder_prediction_tolerance": 0.015}}
142    with open("results.json", "w") as f:
143        json.dump(result, f, indent=2)
144    print(json.dumps(result, indent=2))
145
146
147if __name__ == "__main__":
148    main()