#!/usr/bin/env python3
"""semisymmetry_discriminator.py — Milestone §0: covariant semi-symmetry test.

Places the drift-perturbed Fisher–Rao geometry on the classical ladder
(docs/reviews/CONJECTURE_UPGRADE_ROADMAP.md §0):

    locally symmetric   ∇R = 0                    (Cartan)
    semisymmetric       R·R = 0, ∇R ≠ 0           (Szabó)
    pseudosymmetric     R·R = L·Q(g,R)            (Deszcz)

Geometry. The open simplex Δ_{m−1} = {p ∈ ℝ^m_{>0} : Σp = 1} in reduced
coordinates p₁..p_{m−1} (p_m = 1 − Σ) carries the Fisher–Rao metric

    g_ij = δ_ij / p_i + 1 / p_m            (positive-definite, (m−1, 0)).

0a  BASELINE: g has constant sectional curvature 1/4 (isometric to an
    orthant of the radius-2 sphere via p ↦ 2√p), hence ∇R ≡ 0 — fully
    symmetric, holonomy SO(m−1). Verified symbolically below.

0b  DRIFT PERTURBATION: BraidStateN.lean's rossbyDriftFromChirality
    assigns per-strand weights (left = +1, right = −1, scarred = +1/2,
    achiral = 0). Mean-centering w gives a tangent drift vector β
    (Σβ = 0). The *drift-flipped* metric

        g' = g − 2 (β♭ ⊗ β♭) / g(β, β),        β♭ = g β,

    is the reflection of g along β: it flips the sign of g exactly on
    span(β), so sig(g') = (1, m−2) with the TIME-LIKE direction = the
    drift direction. This is the concrete geometrization of "the
    Kelvin/Rossby directionality supplies the (1, m−2) signature"
    (§0-0b). β is extended as a constant field in the reduced chart —
    a modelling choice, stated explicitly.

0c  DISCRIMINATOR: compute R', ∇R', R'·R' and Q(g',R') for g' and test
    the ladder. Conventions (Kobayashi–Nomizu):

        R^ρ_{σμν} = ∂_μ Γ^ρ_{νσ} − ∂_ν Γ^ρ_{μσ}
                    + Γ^ρ_{μλ} Γ^λ_{νσ} − Γ^ρ_{νλ} Γ^λ_{μσ}
        R_{ρσμν} = g_{ρλ} R^λ_{σμν}

    Derivation action of an endomorphism A on a (0,4) tensor T:

        (A·T)_{ijkl} = −A^m_i T_{mjkl} − A^m_j T_{imkl}
                       −A^m_k T_{ijml} − A^m_l T_{ijkm}

    R·R uses A = R(∂a,∂b) (components A^m_i = R^m_{iab});
    Q(g,R) uses A = ∂a ∧_g ∂b (components A^m_i = δ^m_a g_{bi} − δ^m_b g_{ai}).
    Pseudosymmetry asks for a *function* L with R·R = L·Q(g,R).

All arithmetic exact (sympy Rational); point evaluations at rational
points. A nonzero tensor at one rational point is a PROOF of ≠ 0;
symbolic identities are proven over the whole chart.

Usage:
    python3 python/semisymmetry_discriminator.py            # m = 4 (Δ₃)
    python3 python/semisymmetry_discriminator.py --m 8      # Δ₇ (slow)
"""

from __future__ import annotations

import argparse
import itertools
import sys

import sympy as sp


# ── tensor machinery (exact, chart-based) ─────────────────────────────

def fisher_metric(m: int, ps):
    """Fisher–Rao on Δ_{m−1} in reduced coordinates."""
    pm = 1 - sum(ps)
    n = m - 1
    g = sp.zeros(n, n)
    for i in range(n):
        for j in range(n):
            g[i, j] = (1 / ps[i] if i == j else 0) + 1 / pm
    return sp.Matrix(g)


def christoffel(g: sp.Matrix, coords):
    n = len(coords)
    ginv = g.inv()
    Gamma = [[[sp.S.Zero] * n for _ in range(n)] for _ in range(n)]
    dg = [[[sp.diff(g[i, j], coords[k]) for k in range(n)]
           for j in range(n)] for i in range(n)]
    for l in range(n):
        for i in range(n):
            for j in range(n):
                s = sp.S.Zero
                for k in range(n):
                    s += ginv[l, k] * (dg[k][i][j] + dg[k][j][i] - dg[i][j][k])
                Gamma[l][i][j] = sp.together(s / 2)
    return Gamma


def riemann(g: sp.Matrix, Gamma, coords):
    """R^ρ_{σμν} (up) and R_{ρσμν} (down)."""
    n = len(coords)
    Rup = [[[[sp.S.Zero] * n for _ in range(n)] for _ in range(n)] for _ in range(n)]
    for rho in range(n):
        for sig in range(n):
            for mu in range(n):
                for nu in range(mu + 1, n):  # antisymmetry in (μ,ν)
                    term = sp.diff(Gamma[rho][nu][sig], coords[mu]) \
                         - sp.diff(Gamma[rho][mu][sig], coords[nu])
                    for lam in range(n):
                        term += Gamma[rho][mu][lam] * Gamma[lam][nu][sig] \
                              - Gamma[rho][nu][lam] * Gamma[lam][mu][sig]
                    term = sp.together(term)
                    Rup[rho][sig][mu][nu] = term
                    Rup[rho][sig][nu][mu] = -term
    Rdn = [[[[sp.S.Zero] * n for _ in range(n)] for _ in range(n)] for _ in range(n)]
    for rho in range(n):
        for sig in range(n):
            for mu in range(n):
                for nu in range(mu + 1, n):
                    s = sp.S.Zero
                    for lam in range(n):
                        s += g[rho, lam] * Rup[lam][sig][mu][nu]
                    s = sp.together(s)
                    Rdn[rho][sig][mu][nu] = s
                    Rdn[rho][sig][nu][mu] = -s
    return Rup, Rdn


def cov_deriv_riemann(Rdn, Gamma, coords):
    """∇_a R_{ijkl} — the local-symmetry obstruction."""
    n = len(coords)
    out = [[[[[sp.S.Zero] * n for _ in range(n)] for _ in range(n)]
            for _ in range(n)] for _ in range(n)]
    for a in range(n):
        for i in range(n):
            for j in range(n):
                for k in range(n):
                    for l in range(n):
                        s = sp.diff(Rdn[i][j][k][l], coords[a])
                        for mth in range(n):
                            s -= Gamma[mth][a][i] * Rdn[mth][j][k][l]
                            s -= Gamma[mth][a][j] * Rdn[i][mth][k][l]
                            s -= Gamma[mth][a][k] * Rdn[i][j][mth][l]
                            s -= Gamma[mth][a][l] * Rdn[i][j][k][mth]
                        out[a][i][j][k][l] = sp.together(s)
    return out


def derivation_action(A, Rdn, n):
    """(A·R)_{ijkl} for endomorphism components A[m][i] = A^m_i."""
    out = [[[[sp.S.Zero] * n for _ in range(n)] for _ in range(n)] for _ in range(n)]
    for i in range(n):
        for j in range(n):
            for k in range(n):
                for l in range(n):
                    s = sp.S.Zero
                    for mth in range(n):
                        s -= A[mth][i] * Rdn[mth][j][k][l]
                        s -= A[mth][j] * Rdn[i][mth][k][l]
                        s -= A[mth][k] * Rdn[i][j][mth][l]
                        s -= A[mth][l] * Rdn[i][j][k][mth]
                    out[i][j][k][l] = s
    return out


def rr_and_q(g, Rup, Rdn, n):
    """R·R and Q(g,R): 6-index tensors indexed [a][b] → (0,4) blocks."""
    RR = {}
    Q = {}
    for a in range(n):
        for b in range(a + 1, n):
            A_R = [[Rup[mth][i][a][b] for i in range(n)] for mth in range(n)]
            A_w = [[(g[b, i] if mth == a else 0) - (g[a, i] if mth == b else 0)
                    for i in range(n)] for mth in range(n)]
            RR[(a, b)] = derivation_action(A_R, Rdn, n)
            Q[(a, b)] = derivation_action(A_w, Rdn, n)
    return RR, Q


# ── evaluation helpers ────────────────────────────────────────────────

def eval_tensor_at(t, subsmap, depth_idx):
    """Exact-substitute every component; return list of nonzero (idx, val)."""
    nz = []
    for idx in depth_idx:
        expr = t
        for i in idx:
            expr = expr[i]
        v = sp.nsimplify(sp.together(expr).subs(subsmap))
        v = sp.cancel(v)
        if v != 0:
            nz.append((idx, v))
    return nz


def signature_by_minors(gm: sp.Matrix):
    """Signature via Jacobi/Sylvester: signs of leading principal minors.
    Valid when all leading minors are nonzero (checked)."""
    n = gm.shape[0]
    minors = [gm[:k, :k].det() for k in range(1, n + 1)]
    assert all(mv != 0 for mv in minors), "degenerate leading minor"
    signs = [sp.sign(mv) for mv in minors]
    neg = 0
    prev = 1
    for s in signs:
        if s * prev < 0:
            neg += 1
        prev = s
    return (n - neg, neg), minors


# ── main ──────────────────────────────────────────────────────────────

def main(argv=None) -> int:
    ap = argparse.ArgumentParser()
    ap.add_argument("--m", type=int, default=4,
                    help="ambient outcomes (simplex dim = m−1); default 4")
    ap.add_argument("--fast", action="store_true",
                    help="defer simplification; point-evaluation only "
                         "(≠0 results remain exact proofs; ≡0 results are "
                         "sample-point evidence, confirm symbolically later)")
    args = ap.parse_args(argv)
    if args.fast:
        sp.together = lambda x: x  # simplification deferred to point-eval
    m = args.m
    n = m - 1

    ps = sp.symbols(f"p1:{m}", positive=True)
    coords = list(ps)
    g = fisher_metric(m, ps)

    # rossbyDriftFromChirality weights (BraidStateN.lean):
    #   left=+1, right=−1, scarred=+1/2, achiral=0 — cycled over strands.
    base = [sp.Integer(1), sp.Integer(-1), sp.Rational(1, 2), sp.Integer(0)]
    w = [base[i % 4] for i in range(m)]
    wbar = sp.Rational(sum(w), m)
    beta_amb = [wi - wbar for wi in w]          # Σ = 0: tangent
    beta = sp.Matrix(beta_amb[:n])              # reduced-chart components
    print(f"Δ_{n}: m={m} outcomes, chirality weights w={w}, β_reduced={list(beta)}")

    # ---- 0a: baseline ------------------------------------------------
    print("\n[0a] static Fisher–Rao baseline")
    Gamma = christoffel(g, coords)
    Rup, Rdn = riemann(g, Gamma, coords)
    K = sp.Rational(1, 4)
    csc_ok = True
    for i, j, k, l in itertools.product(range(n), repeat=4):
        expect = K * (g[i, k] * g[j, l] - g[i, l] * g[j, k])
        if sp.simplify(sp.together(Rdn[i][j][k][l] - expect)) != 0:
            csc_ok = False
            break
    print(f"  R = 1/4 (g∧g) (constant curvature 1/4, symbolic): {csc_ok}")
    if csc_ok:
        print("  ⇒ ∇R ≡ 0 (constant-curvature ⇒ locally symmetric): True")
        print("  ⇒ signature (m−1, 0) positive-definite, holonomy SO(m−1)")
    else:
        print("  UNEXPECTED — baseline is not constant curvature; abort")
        return 1

    # ---- 0b: drift-flipped metric -------------------------------------
    print("\n[0b] drift-flipped metric g' = g − 2 β♭⊗β♭ / g(β,β)")
    beta_flat = g * beta
    beta_norm2 = (beta.T * g * beta)[0, 0]
    gp = sp.Matrix(n, n, lambda i, j: sp.together(
        g[i, j] - 2 * beta_flat[i] * beta_flat[j] / beta_norm2))

    # signature at rational sample points (exact minors)
    pts = []
    centroid = {ps[i]: sp.Rational(1, m) for i in range(n)}
    pts.append(("centroid", centroid))
    off1 = {ps[i]: sp.Rational(i + 2, 2 * m * (m + 2)) for i in range(n)}
    pts.append(("offcenter1", off1))
    off2 = {ps[i]: sp.Rational(2 * i + 1, m * m + 3) for i in range(n)}
    pts.append(("offcenter2", off2))

    for name, pt in pts:
        gm = sp.Matrix(n, n, lambda i, j: sp.cancel(gp[i, j].subs(pt)))
        sig, _ = signature_by_minors(gm)
        # time-like = drift: g'(β,β) = −g(β,β) < 0
        bb = sp.cancel((beta.T * gm * beta)[0, 0])
        print(f"  {name}: signature {sig}, g'(β,β) = {bb} (<0 ⇒ β time-like)")

    # ---- 0c: discriminator on g' --------------------------------------
    print("\n[0c] discriminator on g'")
    Gp = christoffel(gp, coords)
    Rpu, Rpd = riemann(gp, Gp, coords)

    idx4 = list(itertools.product(range(n), repeat=4))
    idx5 = list(itertools.product(range(n), repeat=5))

    nablaR = cov_deriv_riemann(Rpd, Gp, coords)
    name0, pt0 = pts[0]
    nz = eval_tensor_at(nablaR, pt0, idx5)
    print(f"  ∇R' at {name0}: {len(nz)} nonzero components "
          f"{'⇒ NOT locally symmetric' if nz else '⇒ vanishes here'}")
    if nz:
        idx, v = nz[0]
        print(f"    e.g. (∇R')_{idx} = {v}")

    RR, Q = rr_and_q(gp, Rpu, Rpd, n)
    max_rr_nz = 0
    ratios = {}
    ok_pseudo = True
    L_val = None
    for (a, b) in RR:
        for name, pt in pts[:2]:
            rr_nz = {i: sp.cancel(sp.together(RRc).subs(pt))
                     for i, RRc in _iter4(RR[(a, b)], idx4)}
            q_nz = {i: sp.cancel(sp.together(Qc).subs(pt))
                    for i, Qc in _iter4(Q[(a, b)], idx4)}
            rr_nz = {i: v for i, v in rr_nz.items() if v != 0}
            q_nz = {i: v for i, v in q_nz.items() if v != 0}
            max_rr_nz = max(max_rr_nz, len(rr_nz))
            if not rr_nz:
                continue
            # pseudosymmetry: R·R = L · Q(g,R) componentwise
            for i, v in rr_nz.items():
                if i not in q_nz:
                    ok_pseudo = False
                    break
                r = sp.cancel(v / q_nz[i])
                ratios.setdefault((name, a, b), set()).add(r)
            if set(rr_nz) != set(q_nz):
                # Q has support where R·R doesn't (or vice versa) → L=0 forced there
                extra = set(q_nz) - set(rr_nz)
                if extra:
                    ratios.setdefault((name, a, b), set()).add(sp.S.Zero)

    if max_rr_nz == 0:
        print("  R'·R' = 0 at sample points (candidate SEMISYMMETRIC — "
              "confirm symbolically)")
        if args.fast:
            verdict = ("SEMISYMMETRIC (Szabó) at sample points — symbolic "
                       "≡ 0 confirmation pending (rerun without --fast)")
        else:
            allzero = all(sp.simplify(sp.together(c)) == 0
                          for (a, b) in RR for _, c in _iter4(RR[(a, b)], idx4))
            print(f"  R'·R' ≡ 0 symbolically: {allzero}")
            verdict = "SEMISYMMETRIC (Szabó)" if allzero and nz else "check"
    else:
        per_point = {}
        for (name, a, b), rs in ratios.items():
            per_point.setdefault(name, set()).update(rs)
        for name, rs in per_point.items():
            print(f"  L candidates at {name}: {rs}")
        consistent = all(len(rs) == 1 for rs in per_point.values()) and ok_pseudo
        if consistent:
            L_val = {name: next(iter(rs)) for name, rs in per_point.items()}
            verdict = f"PSEUDOSYMMETRIC (Deszcz), L = {L_val}"
        else:
            verdict = ("PROPER (neither semi- nor pseudo-symmetric at sample "
                       "points) — R·R ≠ L·Q(g,R)")
    print(f"\n[verdict] baseline: LOCALLY SYMMETRIC (∇R=0). "
          f"drift-perturbed: {'∇R ≠ 0; ' if nz else ''}{verdict}")

    # ---- 0d: drift-direction diagnostics ------------------------------
    print("\n[0d] obstruction vs drift direction (centroid)")
    ib = [sp.S.Zero] * 4
    # contract ∇R' with β in the derivative slot vs a g-orthogonal vector
    contr_beta = sp.S.Zero
    for a in range(n):
        for idx in idx4:
            comp = nablaR[a][idx[0]][idx[1]][idx[2]][idx[3]]
            contr_beta += (beta[a] * comp.subs(pt0)) ** 2
    contr_beta = sp.cancel(contr_beta)
    print(f"  Σ (β^a ∇_a R')² at centroid = {sp.nsimplify(contr_beta)} "
          f"({'drift direction carries obstruction' if contr_beta != 0 else 'obstruction ⟂ drift'})")
    return 0


def _iter4(t, idx4):
    for idx in idx4:
        yield idx, t[idx[0]][idx[1]][idx[2]][idx[3]]


if __name__ == "__main__":
    sys.exit(main())