1import json
  2import numpy as np
  3
  4EPS = 1e-6
  5
  6
  7def admissibility(U, eps=EPS):
  8    U = np.asarray(U)
  9    D = U[..., 0]
 10    m = U[..., 1:-1]
 11    E = U[..., -1]
 12    return (D >= eps) & (E - np.sqrt(D * D + np.sum(m * m, axis=-1)) >= eps)
 13
 14
 15def g(U, s):
 16    U = np.asarray(U)
 17    s = np.asarray(s)
 18    # Supports a shared s=(d,) or batched s=(N,d).
 19    if s.ndim == 1:
 20        return U[..., -1] - s[0] * U[..., 0] - np.sum(s[1:] * U[..., 1:-1], axis=-1)
 21    return U[..., -1] - s[..., 0] * U[..., 0] - np.sum(s[..., 1:] * U[..., 1:-1], axis=-1)
 22
 23
 24def limit_residual(U0, U1, eps=EPS, iterations=60):
 25    """Exact scalar residual limiter via bisection on the relativistic margin.
 26
 27    The admissible set is convex, so with valid U0 and invalid U1 there is one
 28    boundary crossing on [0,1]. This is a robust reference implementation of
 29    the scalar limiter; theta is detached because this NumPy MVP has no grads.
 30    """
 31    U0 = np.asarray(U0, dtype=float)
 32    U1 = np.asarray(U1, dtype=float)
 33    delta = U1 - U0
 34    def margin(t):
 35        V = U0 + t[..., None] * delta
 36        return V[..., -1] - np.sqrt(np.sum(V[..., :-1] * V[..., :-1], axis=-1))
 37    theta = np.ones(U0.shape[:-1])
 38    bad = (U1[..., 0] < eps) | (margin(np.ones(U0.shape[:-1])) < eps)
 39    lo = np.zeros_like(theta)
 40    hi = np.ones_like(theta)
 41    # Only bad rows need a root; valid rows remain at theta=1.
 42    for _ in range(iterations):
 43        mid = 0.5 * (lo + hi)
 44        ok = margin(mid) >= eps
 45        lo = np.where(ok, mid, lo)
 46        hi = np.where(ok, hi, mid)
 47    theta = np.where(bad, lo, theta)
 48    # Density is linear and can be the active constraint independently.
 49    dbad = U1[..., 0] < eps
 50    td = np.clip((U0[..., 0] - eps) / np.maximum(U0[..., 0] - U1[..., 0], 1e-30), 0., 1.)
 51    theta = np.where(dbad, np.minimum(theta, td), theta)
 52    # Move infinitesimally inside both closed constraints to avoid roundoff.
 53    theta = np.where(theta < 1.0, np.maximum(0.0, theta * (1.0 - 1e-12)), theta)
 54    V = U0 + theta[..., None] * delta
 55    # A direct safety correction handles ulp-scale endpoint errors.
 56    V[..., 0] = np.maximum(V[..., 0], eps)
 57    return V, theta
 58
 59def component_clip(U, eps=EPS):
 60    """Control: independently clip density and energy, leaving momentum unchanged."""
 61    V = np.asarray(U, dtype=float).copy()
 62    V[..., 0] = np.maximum(V[..., 0], eps)
 63    V[..., -1] = np.maximum(V[..., -1], eps)
 64    return V
 65
 66
 67def main(seed=1234):
 68    rng = np.random.default_rng(seed)
 69
 70    # Verify min_{||s||<=1} g_s(U) = E - sqrt(D^2+||m||^2).
 71    x = rng.normal(size=(10000, 3))
 72    x /= np.linalg.norm(x, axis=1, keepdims=True)
 73    U = rng.uniform(0.1, 3.0, size=(10000, 4))
 74    U[:, 1:3] = rng.normal(0.0, 1.0, size=(10000, 2))
 75    exact = U[:, -1] - np.linalg.norm(U[:, :-1], axis=1)
 76    sampled = np.min(np.stack([g(U, x[i]) for i in range(len(x))], axis=1), axis=1)
 77    # The sampled minimum is only an approximation; test the analytic worst normal.
 78    w = U[:, :-1]
 79    sstar = w / np.linalg.norm(w, axis=1, keepdims=True)
 80    analytic = np.array([g(U[i], sstar[i]) for i in range(len(U))])
 81    identity_abs_err = float(np.max(np.abs(exact - analytic)))
 82    sampled_gap_mean = float(np.mean(sampled - exact))
 83
 84    # Valid baselines, then deliberately corrupt learned candidates.
 85    n = 20000
 86    U0 = rng.uniform(0.2, 2.0, size=(n, 4))
 87    U0[:, 1:3] = rng.normal(0.0, 0.7, size=(n, 2))
 88    U0[:, -1] = np.linalg.norm(U0[:, :-1], axis=1) + rng.uniform(0.05, 1.5, size=n)
 89    # Errors are large enough to produce both density and coupled energy/momentum violations.
 90    U1 = U0 + rng.normal(0.0, 0.9, size=U0.shape)
 91    U1[:, 0] -= rng.uniform(0.0, 1.3, size=n)
 92    U1[:, -1] -= rng.uniform(0.0, 0.9, size=n)
 93
 94    safe, theta = limit_residual(U0, U1)
 95    clipped = component_clip(U1)
 96    valid0 = admissibility(U0).mean()
 97    valid1 = admissibility(U1).mean()
 98    valid_safe = admissibility(safe).mean()
 99    valid_clip = admissibility(clipped).mean()
100    # Retention is the fraction of the learned residual preserved.
101    dnorm = np.linalg.norm(U1 - U0, axis=1)
102    retained = np.linalg.norm(safe - U0, axis=1) / np.maximum(dnorm, 1e-12)
103    clip_retained = np.linalg.norm(clipped - U0, axis=1) / np.maximum(dnorm, 1e-12)
104    result = {
105        "seed": seed,
106        "identity_max_abs_error": identity_abs_err,
107        "random_sphere_mean_gap_vs_exact": sampled_gap_mean,
108        "baseline_valid_rate": float(valid0),
109        "candidate_valid_rate": float(valid1),
110        "gql_safe_valid_rate": float(valid_safe),
111        "component_clip_valid_rate": float(valid_clip),
112        "gql_mean_theta": float(theta.mean()),
113        "gql_fraction_limited": float(np.mean(theta < 1.0 - 1e-12)),
114        "gql_mean_residual_retention": float(retained.mean()),
115        "component_clip_mean_residual_retention": float(clip_retained.mean()),
116        "gql_min_margin": float(np.min(safe[:, -1] - np.linalg.norm(safe[:, :-1], axis=1))),
117        "component_clip_min_margin": float(np.min(clipped[:, -1] - np.linalg.norm(clipped[:, :-1], axis=1))),
118    }
119    print(json.dumps(result, indent=2, sort_keys=True))
120    with open("gql_results.json", "w") as f:
121        json.dump(result, f, indent=2, sort_keys=True)
122
123
124if __name__ == "__main__":
125    main()
126