1import json
  2import time
  3import numpy as np
  4from scipy.sparse import csr_matrix
  5from scipy.sparse.linalg import spsolve
  6
  7
  8def shift_matrix(n, k=1):
  9    P = np.zeros((n, n))
 10    for i in range(n):
 11        P[(i + k) % n, i] = 1.0
 12    return P
 13
 14
 15def cycle_laplacian(n):
 16    return 2*np.eye(n) - shift_matrix(n, 1) - shift_matrix(n, -1)
 17
 18
 19def symmetry_basis(n):
 20    j = np.arange(n)[:, None]
 21    k = np.arange(n)[None, :]
 22    return np.exp(2j*np.pi*j*k/n) / np.sqrt(n)
 23
 24
 25def projector_checks(n):
 26    Ps = []
 27    for k in range(n):
 28        P = np.zeros((n, n), dtype=complex)
 29        for g in range(n):
 30            rho = shift_matrix(n, g)
 31            chi = np.exp(2j*np.pi*k*g/n)
 32            P += np.conj(chi) * rho / n
 33        Ps.append(P)
 34    return {
 35        "resolution_error": float(np.linalg.norm(sum(Ps)-np.eye(n))),
 36        "idempotence_error": float(max(np.linalg.norm(P@P-P) for P in Ps)),
 37        "orthogonality_error": float(max(np.linalg.norm(Ps[i]@Ps[j]) for i in range(n) for j in range(n) if i != j)),
 38    }
 39
 40
 41def matrices(n, theta, tau):
 42    L = cycle_laplacian(n)
 43    # Symmetry-compatible learned combination of identity, Laplacian, and L^2.
 44    H = 0.7*np.eye(n) + 0.4*L + 0.08*(L@L)
 45    K = np.eye(n) + tau * theta * H
 46    return L, K
 47
 48
 49def one_case(n, theta, tau, rng):
 50    L, K = matrices(n, theta, tau)
 51    Q = symmetry_basis(n)
 52    rho = shift_matrix(n, 1)
 53    comm = np.linalg.norm(rho@L-L@rho) / np.linalg.norm(L)
 54    Kh = Q.conj().T @ K @ Q
 55    offdiag = Kh - np.diag(np.diag(Kh))
 56    block_error = np.linalg.norm(offdiag) / np.linalg.norm(Kh)
 57    b = rng.normal(size=n)
 58    t0 = time.perf_counter()
 59    x_sparse = spsolve(csr_matrix(K), b)
 60    sparse_ms = 1000*(time.perf_counter()-t0)
 61    t0 = time.perf_counter()
 62    bh = Q.conj().T @ b
 63    # C_n irreps are one-dimensional: independent scalar block solves.
 64    xh = bh / np.diag(Kh)
 65    x_fourier = np.real_if_close(Q @ xh).real
 66    fourier_ms = 1000*(time.perf_counter()-t0)
 67    residual = np.linalg.norm(K@x_fourier-b)/np.linalg.norm(b)
 68    equivalence = np.linalg.norm(x_fourier-x_sparse)/np.linalg.norm(x_sparse)
 69    return comm, block_error, equivalence, residual, sparse_ms, fourier_ms
 70
 71
 72def main():
 73    rng = np.random.default_rng(1234)
 74    ns = [16, 32, 64, 128, 256]
 75    thetas = [0.0, 0.1, 1.0, 10.0]
 76    taus = [0.01, 0.1, 1.0]
 77    proj = projector_checks(16)
 78    rows = []
 79    for n in ns:
 80        for theta in thetas:
 81            for tau in taus:
 82                vals = one_case(n, theta, tau, rng)
 83                rows.append({"n":n,"theta":theta,"tau":tau,"comm":vals[0],"block_error":vals[1],"equivalence":vals[2],"residual":vals[3],"sparse_ms":vals[4],"fourier_ms":vals[5]})
 84    by_n = []
 85    for n in ns:
 86        rr = [r for r in rows if r["n"] == n]
 87        by_n.append({"n": n,
 88                     "median_sparse_ms": float(np.median([r["sparse_ms"] for r in rr])),
 89                     "median_fourier_ms": float(np.median([r["fourier_ms"] for r in rr])),
 90                     "fourier_over_sparse": float(np.median([r["fourier_ms"] for r in rr]) / np.median([r["sparse_ms"] for r in rr]))})
 91    tol = 1e-12
 92    predictions = [
 93        {"prediction": "commutator remains at roundoff for every theta,tau", "predicted_bound": tol,
 94         "observed_max": max(r["comm"] for r in rows), "confirmed": bool(max(r["comm"] for r in rows) < tol)},
 95        {"prediction": "symmetry basis removes off-block entries for every theta,tau", "predicted_bound": tol,
 96         "observed_max": max(r["block_error"] for r in rows), "confirmed": bool(max(r["block_error"] for r in rows) < tol)},
 97        {"prediction": "transformed solve equals original solve for every theta,tau", "predicted_bound": tol,
 98         "observed_max": max(r["equivalence"] for r in rows), "confirmed": bool(max(r["equivalence"] for r in rows) < tol)},
 99    ]
100    summary = {
101        "sweep": {"n": ns, "theta": thetas, "tau": taus},
102        "predictions": predictions,
103        "timing_by_n": by_n,
104        "projectors": proj,
105        "max_commutator": max(r["comm"] for r in rows),
106        "max_block_error": max(r["block_error"] for r in rows),
107        "max_solution_relative_error": max(r["equivalence"] for r in rows),
108        "max_residual": max(r["residual"] for r in rows),
109        "median_sparse_ms": float(np.median([r["sparse_ms"] for r in rows])),
110        "median_fourier_ms": float(np.median([r["fourier_ms"] for r in rows])),
111        "rows": rows,
112    }
113    with open("results.json", "w") as f:
114        json.dump(summary, f, indent=2)
115    print(json.dumps({k:v for k,v in summary.items() if k != "rows"}, indent=2))
116
117
118if __name__ == "__main__":
119    main()