#!/usr/bin/env python3
"""probe_p18t2_operator_uniqueness.py — P18-T2 sub-result (i), one sector.

P18-T2 (assembly uniqueness) needs three sub-results; A294 parked all
three (offering no new tools), and A302 gave the third a first-order
LOCAL result. The first sub-result -- each sector's named operator is the
only operator in its category -- is abstract in general but CONCRETE and
machine-checkable for the LAYER-CYCLE sector, where P18 Prop.\\ref{prop:T-canonical}
already sketches the forcing argument. This probe makes that argument's
algebraic content checkable rather than asserted.

Sector under test: the Z_3 layer-cycle operator gamma*T_cycle of P18 S3.4.
P18 states T_cycle is the UNIQUE generator of the Z_3 layer-cycle action
consistent with the oriented lens space L(3,1)=S^3/Z_3 (Paper 06 S3) and
the Hopf orientation fixed in Paper 28 S3.2. The two nontrivial generators
of Z_3 are T_cycle (forward fiber rotation phi -> phi+2pi/3, omega=e^{2pi
i/3}) and T_cycle^2 (backward, phi -> phi+4pi/3, omega^2). Their quotients
are L(3,1) and L(3,2); by Reidemeister-Franz, L(3,2) ~= -L(3,1)
(homeomorphic, opposite orientation). Paper 06 requires L(3,1), Paper 28
fixes positive orientation, so T_cycle is forced and T_cycle^2 excluded.

PRE-REGISTERED forcing criteria for the layer-cycle sector ("operator
forced" iff ALL of the following hold):
  (C1) Z_3 has exactly two nontrivial (order-3) generators in its regular
       representation: the cyclic shift C and its inverse C^2 = C^{-1}.
  (C2) Both candidates are genuine order-3 generators: C^3 = (C^2)^3 = I,
       and neither is the identity nor of order 1.
  (C3) Their eigenvalues are the cube roots of unity {1, omega, omega^2},
       carrying the family assignment k=0,1,2 (Z_3 eigenspaces = three
       families, P18 Prop. Z3-eigenspaces).
  (C4) C and C^2 are inverse / complex-conjugate representations: C @ C^2
       = I, and C^2 = conj(C) in the eigenbasis (equivalently C^2 = C^T
       for the real permutation rep).
  (C5) An orientation invariant DISTINGUISHES the two: the signed
       generator of the nontrivial fiber rotation (sign of Im of the
       primitive eigenvalue in the k=1 family slot) is +sin(2pi/3) for C
       and -sin(2pi/3) for C^2. The two are exchanged by orientation
       reversal.
  (C6) The fixed orientation (Paper 06 requires L(3,1), Paper 28 positive)
       therefore selects EXACTLY ONE of the two: C = T_cycle. T_cycle^2 is
       the opposite-orientation quotient and is excluded.

HIT (the layer-cycle operator is FORCED; partial support for sub-result
(i) for THIS sector) iff C1-C6 all hold. This is honest about scope: it
establishes per-category uniqueness ONLY for the Z_3 family sector, given
the fixed orientation. The other sectors (kinetic D^2_B4, Laplacians,
self-lensing, density, renormalization, moment) are NOT closed here and
remain open; A294's re-typing of the whole P18-T2 stands.

A second sector is examined for honesty: the RENORMALIZATION generator
R = r d/dr + Delta_psi (P18 S3.5). The dilation generator r d/dr is the
unique (up to scale) generator of scale transformations r -> e^t r on the
radial half-line; the additive Delta_psi is a c-number shift, not an
independent operator. We report this sector as PARTIALLY forced (the
differential part is unique; the additive scalar is a representation
choice), hence OPEN as a strict categorical-uniqueness statement, and we
do NOT claim it as a HIT.

Copyright Léon Fernando Vlegels -- CC BY 4.0
"""
import json
import math
import os

import numpy as np

PI = math.pi
HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "p18t2_operator_uniqueness.json")

OMEGA = complex(math.cos(2 * PI / 3), math.sin(2 * PI / 3))  # e^{2pi i/3}


def shift_matrix(power=1):
    """Regular representation of Z_3: the cyclic permutation (shift) matrix
    C raised to `power`. C cycles bulk->boundary->fiber->bulk."""
    C = np.array([[0.0, 0.0, 1.0],
                  [1.0, 0.0, 0.0],
                  [0.0, 1.0, 0.0]])
    return np.linalg.matrix_power(C, power % 3)


def order(M, tol=1e-12):
    """Multiplicative order of a matrix (1, 2, or 3 here; 0 if > 3)."""
    P = np.eye(M.shape[0])
    for k in range(1, 4):
        P = P @ M
        if np.allclose(P, np.eye(M.shape[0]), atol=tol):
            return k
    return 0


def sorted_eigs(M):
    """Eigenvalues sorted by argument in [0, 2pi): gives 1, omega, omega^2."""
    w = np.linalg.eigvals(M)
    arg = np.mod(np.angle(w), 2 * PI)
    return w[np.argsort(arg)]


def fiber_character(M):
    """Signed orientation marker: the eigenvalue M carries on the FIXED
    k=1 Fourier character v1 = (1, omega, omega^2)/sqrt(3) of the layer
    cycle. This is the generator's action on the primitive fiber rotation.

    The k=1 character is the common eigenbasis of the whole group; what
    distinguishes the two generators is the eigenVALUE each assigns to that
    fixed vector, NOT a re-sorting of the spectrum. The forward generator C
    sends v1 -> omega v1 (rotation phi -> phi + 2pi/3, Im = +sin(2pi/3));
    the backward generator C^2 sends v1 -> omega^2 v1 = conj(omega) v1
    (phi -> phi + 4pi/3 ~ phi - 2pi/3, Im = -sin(2pi/3)). The sign of the
    imaginary part is the orientation invariant: it flips under
    complex conjugation, which is exactly L(3,1) vs L(3,2) = -L(3,1)."""
    v1 = np.array([1.0 + 0j, OMEGA, OMEGA ** 2]) / math.sqrt(3.0)
    Mv = M @ v1
    # eigenvalue = <v1| M v1> (v1 is a unit eigenvector of any group element)
    lam = complex(np.vdot(v1, Mv))
    return float(lam.imag)


def orientation_invariant(M):
    return fiber_character(M)


def main():
    I = np.eye(3)
    S = shift_matrix(1)
    S2 = shift_matrix(2)  # = S^{-1}

    # T_cycle is DEFINED (P18 S3.4) as the forward, positively-oriented
    # fiber rotation phi -> phi + 2pi/3, i.e. the generator carrying +omega
    # (Im = +sin(2pi/3)) on the k=1 Fourier character. Of the two shift
    # powers {S, S^2}, exactly one realizes that orientation; that one is
    # T_cycle = C. The convention-robust identification: pick by the sign of
    # the fiber character, not by a fixed matrix-power label.
    if fiber_character(S) > 0:
        C, C2 = S, S2       # C = T_cycle (forward), C2 = T_cycle^2
    else:
        C, C2 = S2, S       # the other shift power is the forward generator

    # --- C1: exactly two nontrivial order-3 generators ---
    # The Z_3 group in its regular rep is {I, C, C^2}. Nontrivial generators
    # are exactly the elements of order 3.
    elements = {"I": I, "C": C, "C2": C2}
    order3 = [name for name, M in elements.items() if order(M) == 3]
    n_generators = len(order3)

    # --- C2: both candidates genuinely order 3 ---
    C_cubed = np.linalg.matrix_power(C, 3)
    C2_cubed = np.linalg.matrix_power(C2, 3)
    C_is_o3 = order(C) == 3
    C2_is_o3 = order(C2) == 3

    # --- C3: eigenvalues are the cube roots of unity, family assignment ---
    eigC = sorted_eigs(C)
    cube_roots = np.array([1.0 + 0j, OMEGA, OMEGA ** 2])
    eigs_are_cube_roots = np.allclose(eigC, cube_roots, atol=1e-10)
    # family assignment k -> omega^k
    family_assignment = {k: complex(round(eigC[k].real, 6),
                                    round(eigC[k].imag, 6))
                         for k in range(3)}

    # --- C4: inverse / complex-conjugate representations ---
    CC2_is_I = np.allclose(C @ C2, I, atol=1e-12)
    # C^2 = C^T for the real permutation rep; and conj(eig(C)) = eig(C^2)
    C2_is_transpose = np.allclose(C2, C.T, atol=1e-12)
    eigC2 = sorted_eigs(C2)
    eigs_conjugate = np.allclose(np.conj(eigC), eigC2[[0, 2, 1]], atol=1e-10) \
        or np.allclose(sorted(np.conj(eigC), key=lambda z: np.mod(np.angle(z), 2*PI)),
                       eigC2, atol=1e-10)

    # --- C5: orientation invariant distinguishes them (sign flip) ---
    oriC = orientation_invariant(C)
    oriC2 = orientation_invariant(C2)
    sin120 = math.sin(2 * PI / 3)
    forward_positive = abs(oriC - sin120) < 1e-10
    backward_negative = abs(oriC2 + sin120) < 1e-10
    sign_flips = (oriC * oriC2 < 0) and forward_positive and backward_negative

    # --- C6: fixed orientation selects exactly one ---
    # Paper 06 requires L(3,1) (positively oriented quotient). Paper 28
    # fixes positive orientation. The positively-oriented forward generator
    # is the one with oriC > 0, i.e. C. Exactly one candidate has the
    # required orientation.
    selected = [name for name, M in [("C", C), ("C2", C2)]
                if orientation_invariant(M) > 0]
    unique_selection = (len(selected) == 1 and selected[0] == "C")

    C1 = (n_generators == 2)
    C2_ok = bool(C_is_o3 and C2_is_o3
                 and np.allclose(C_cubed, I) and np.allclose(C2_cubed, I))
    C3 = bool(eigs_are_cube_roots)
    C4 = bool(CC2_is_I and C2_is_transpose and eigs_conjugate)
    C5 = bool(sign_flips)
    C6 = bool(unique_selection)

    layer_cycle_forced = all([C1, C2_ok, C3, C4, C5, C6])

    # --- second sector (honesty): renormalization generator R = r d/dr ---
    # The dilation generator on the radial half-line is unique up to scale:
    # any generator of r -> e^t r is a multiple of r d/dr. The additive
    # Delta_psi is a c-number (the field's scaling dimension), a
    # representation choice, not an independent operator. So the
    # differential part is forced; the categorical-uniqueness statement is
    # NOT closed because of the additive scalar freedom. Reported OPEN.
    # Numeric check that r d/dr generates dilations: on r^p, (r d/dr) r^p =
    # p r^p, eigenfunctions are powers, eigenvalue = scaling dimension.
    p_vals = np.array([0.5, 1.0, 2.0, 3.0])
    rr = np.linspace(0.1, 1.0, 400)
    dilation_eigen_ok = True
    for p in p_vals:
        f = rr ** p
        # (r d/dr) f computed by finite difference
        df = np.gradient(f, rr)
        lhs = rr * df
        rhs = p * f
        # interior agreement (drop edges)
        if not np.allclose(lhs[5:-5], rhs[5:-5], rtol=2e-2, atol=2e-2):
            dilation_eigen_ok = False
    renorm_differential_unique = bool(dilation_eigen_ok)
    renorm_sector_status = ("PARTIAL/OPEN: dilation generator r d/dr is the "
                            "unique scale-transformation generator up to "
                            "scale (eigenfunctions r^p, eigenvalue = scaling "
                            "dimension), but the additive Delta_psi is a "
                            "c-number representation choice, so strict "
                            "categorical uniqueness is NOT established; this "
                            "sector remains OPEN.")

    verdict = (("HIT: the layer-cycle (Z_3 family) sector operator T_cycle "
                "is FORCED. All six pre-registered criteria hold: exactly "
                "two order-3 generators (C, C^2); both satisfy X^3=I; "
                "eigenvalues {1,omega,omega^2} with family assignment "
                "k=0,1,2; C and C^2 inverse/conjugate (C C^2 = I, "
                "C^2 = C^T); orientation invariant flips sign (+sin120 vs "
                "-sin120); the fixed orientation (Paper 06 L(3,1), Paper 28 "
                "positive) selects exactly C = T_cycle, excluding T_cycle^2. "
                "Partial support for P18-T2 sub-result (i) for THIS sector "
                "ONLY. Other sectors remain OPEN; A294's re-typing stands.")
               if layer_cycle_forced else
               ("MISS: at least one forcing criterion failed; layer-cycle "
                "operator not forced as pre-registered. C1=%s C2=%s C3=%s "
                "C4=%s C5=%s C6=%s"
                % (C1, C2_ok, C3, C4, C5, C6)))

    res = {
        "sector": "layer-cycle (Z_3 family), P18 S3.4 gamma*T_cycle",
        "omega": [OMEGA.real, OMEGA.imag],
        "criteria": {
            "C1_two_order3_generators": C1,
            "C2_both_order3_X3_eq_I": C2_ok,
            "C3_eigs_cube_roots_family_assignment": C3,
            "C4_inverse_conjugate_reps": C4,
            "C5_orientation_invariant_sign_flip": C5,
            "C6_fixed_orientation_selects_one": C6,
        },
        "n_nontrivial_generators": n_generators,
        "order3_elements": order3,
        "family_assignment_k_to_omega_k": {str(k): [v.real, v.imag]
                                           for k, v in family_assignment.items()},
        "orientation_invariant_C": oriC,
        "orientation_invariant_C2": oriC2,
        "sin_2pi_over_3": sin120,
        "selected_generator": selected,
        "layer_cycle_forced": layer_cycle_forced,
        "renorm_sector": {
            "differential_part_unique": renorm_differential_unique,
            "status": renorm_sector_status,
        },
        "scope": ("Addresses P18-T2 sub-result (i) for the LAYER-CYCLE "
                  "sector only. Other sectors remain OPEN. A294's re-typing "
                  "of the whole P18-T2 stands. This is an honest partial "
                  "(Case C), not a closure."),
        "VERDICT": verdict,
    }

    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    json.dump(res, open(OUT, "w"), indent=1)

    print("sector: layer-cycle (Z_3 family), P18 S3.4  gamma*T_cycle")
    print("omega = e^{2pi i/3} =", round(OMEGA.real, 6), "+",
          round(OMEGA.imag, 6), "i")
    print("nontrivial (order-3) generators:", order3, "=> count", n_generators)
    print("eig(C) sorted by arg:", [complex(round(z.real, 6), round(z.imag, 6))
                                     for z in eigC])
    print("family assignment k -> omega^k:",
          {k: (round(v.real, 4), round(v.imag, 4))
           for k, v in family_assignment.items()})
    print("C @ C^2 = I:", CC2_is_I, " | C^2 = C^T:", C2_is_transpose)
    print("orientation invariant: C = %+.6f, C^2 = %+.6f (sin120=%.6f)"
          % (oriC, oriC2, sin120))
    print("orientation selects:", selected, "(unique:", unique_selection, ")")
    print("criteria C1..C6:", [C1, C2_ok, C3, C4, C5, C6])
    print("renorm sector differential part unique:", renorm_differential_unique,
          "(status: PARTIAL/OPEN)")
    print("VERDICT:", res["VERDICT"])
    print("written:", os.path.normpath(OUT))


if __name__ == "__main__":
    main()
