"""
gate_residual_recovery.py
-------------------------
Initial recovery for the gravity-consistency gate (to be stress-tested).

The gate (from GEOMETRIC_SUBSTANCE_CANONICAL.md): the model points to
"gravity does not exist" ONLY in the flat-and-torsionless case — braid residual
R_ij ≡ 0. Otherwise the geometry gravitates: teleparallel if the curvature part
is flat (gravity carried by torsion), Einstein-Cartan if it is not.

This computes a faithful PROTOTYPE of the residual (not the real BraidField
codec): R_ij = B_ij - (B_i + B_j) on a phase-circle model.

    B_i  = cos(theta_i)                      (single-strand polarity)
    B_ij = cos(theta_i + theta_j)            (symmetric joint  -> torsion)
         + eps * chi_ij * sin(theta_i-theta_j)  (antisymmetric -> curvature)
    R_ij = B_ij - B_i - B_j                  (the connected/interaction residual)

Symmetric part of R = torsion-like; antisymmetric part = curvature-like.
Achiral (eps=0) => curvature is EXACTLY zero => teleparallel.

Structured input = a Mian-Chowla Sidon set (all pairwise sums distinct).
The Sidon-vs-random separation is measured on the defining property itself
(pairwise-sum collisions), which is also what makes the residual non-degenerate.
"""

import numpy as np

rng = np.random.default_rng(0)


def mian_chowla(n):
    """Greedy B_2 (Sidon) sequence: all pairwise sums distinct, by construction."""
    seq, sums, c = [1], {2}, 2
    while len(seq) < n:
        new, ok = [], True
        for s in seq:
            v = s + c
            if v in sums or v in new:
                ok = False; break
            new.append(v)
        if ok and (c + c) not in sums and (c + c) not in new:
            for s in seq:
                sums.add(s + c)
            sums.add(c + c)
            seq.append(c)
        c += 1
    return np.array(seq)


def sum_collisions(s):
    """Defining Sidon failure count: #(pairs) - #(distinct pairwise sums)."""
    s = np.asarray(s)
    pairs = [s[i] + s[j] for i in range(len(s)) for j in range(i + 1, len(s))]
    return len(pairs) - len(set(pairs))


def residual(s, N, eps=0.0, chi=None):
    th = 2 * np.pi * np.asarray(s, float) / N
    Bi = np.cos(th)
    sym = np.cos(th[:, None] + th[None, :])                 # torsion-bearing
    if eps != 0.0 and chi is not None:
        asym = eps * chi * np.sin(th[:, None] - th[None, :])  # curvature-bearing
    else:
        asym = 0.0
    Bij = sym + asym
    return Bij - Bi[:, None] - Bi[None, :]


def split(R):
    Rs = 0.5 * (R + R.T)          # symmetric  -> torsion
    Ra = 0.5 * (R - R.T)          # antisymmetric -> curvature
    fro = np.linalg.norm(R)
    return np.linalg.norm(Rs), np.linalg.norm(Ra), fro


if __name__ == "__main__":
    n = 16
    S = mian_chowla(n)
    M = int(S.max())
    N = 2 * M + 1                 # phase modulus chosen so sums never wrap-collide
    print(f"Sidon set (Mian-Chowla, n={n}): {S.tolist()}")
    print(f"max={M}, phase modulus N={N}")
    print(f"sum-collisions (structured): {sum_collisions(S)}  (0 == perfectly Sidon)\n")

    # --- Gate condition 1: torsion present, curvature class -------------------
    # achiral
    R0 = residual(S, N, eps=0.0)
    tor0, cur0, fro0 = split(R0)
    # chiral (imbalanced handedness): fixed left-handed crossings
    chi = np.ones((n, n))
    Rc = residual(S, N, eps=0.5, chi=chi)
    torc, curc, froc = split(Rc)

    print("RESIDUAL / CURVATURE-TORSION SPLIT")
    print(f"  achiral:  ||R||={fro0:.4f}  torsion(sym)={tor0:.4f}  curvature(asym)={cur0:.2e}")
    print(f"            -> curvature {'== 0  => TELEPARALLEL' if cur0 < 1e-9 else '!= 0'}")
    print(f"  chiral :  ||R||={froc:.4f}  torsion(sym)={torc:.4f}  curvature(asym)={curc:.4f}")
    print(f"            -> curvature {'!= 0  => EINSTEIN-CARTAN' if curc > 1e-9 else '== 0'}")
    print(f"  (gravity-denying would require ||R|| == 0; here ||R|| = {fro0:.4f})\n")

    # --- Gate condition 2: structured vs random separation -------------------
    TRIALS = 5000
    rand_coll = np.empty(TRIALS)
    for t in range(TRIALS):
        r = rng.choice(np.arange(1, M + 1), size=n, replace=False)
        rand_coll[t] = sum_collisions(r)
    mu, sd = rand_coll.mean(), rand_coll.std()
    z = (sum_collisions(S) - mu) / sd
    frac_random_clean = np.mean(rand_coll == 0)

    print(f"SIDON SEPARATION vs {TRIALS} random {n}-subsets of [1,{M}]")
    print(f"  structured sum-collisions : 0")
    print(f"  random  mean +/- std      : {mu:.2f} +/- {sd:.2f}")
    print(f"  z-score                   : {z:.2f}")
    print(f"  random sets that are Sidon: {100*frac_random_clean:.2f}%\n")

    # --- Gate verdict --------------------------------------------------------
    flat = fro0 < 1e-9
    print("GATE VERDICT (initial recovery, to be stress-tested)")
    print(f"  flat-and-torsionless (gravity-denying)? {flat}")
    if not flat:
        cls = "TELEPARALLEL (curvature 0, gravity in torsion)" if cur0 < 1e-9 else "Einstein-Cartan"
        print(f"  geometry gravitates; achiral class: {cls}")
        print(f"  structure beats random by z = {z:.1f}  -> 'something rather than nothing'")
        print("  => model AFFIRMS gravity; does NOT point to gravity not existing. GATE OPENS.")
    else:
        print("  => residual flat; gate stays CLOSED.")