"""Finite commutant comparisons on two independently chosen matrix models.

The benchmark retains the historical Camino-7 calculations but narrows their
scope. On a complete graph it checks that four explicitly listed matrices
``f(A, L)`` commute with the enumerated permutation automorphisms and preserve
the invariant subspace from one symmetric seed. On a separate tensor-product
model it checks that lifts ``M (x) I_2`` commute with sampled transformations
``I (x) U``, scalar phase matrices commute, and selected SU(2) matrices need
not commute.

These checks establish only that the tested finite matrix family is a subset of
the relevant commutants. They neither characterize the engine's reachable set
as an entire commutant nor construct a map from a Riemann or Yang--Mills target
into either finite space. The complementary vector used below is synthetic; it
is not ``S(T)``. The colour model and its U(1)/SU(2) groups are external choices,
not a gauge selected or generated by TNFR. Similar-looking orthogonal splittings
therefore do not identify the two Millennium-problem obstructions.

These finite identities are not a derivation from the nodal equation.

Run:
    python benchmarks/commutant_bridge.py

Status: RESEARCH benchmark; finite algebra and negative identification result.
"""

from __future__ import annotations

import os
import sys

import networkx as nx
import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
# Robust fallback so the harness also runs without PYTHONPATH=src preset.
sys.path.insert(
    0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "src")
)
from composition_arithmetic import automorphism_matrices  # noqa: E402

# Optional source of an externally specified ``log(p)`` diagonal.
try:  # pragma: no cover - exercised only when the package is importable
    from tnfr.dynamics.adelic import AdelicDynamics  # noqa: E402

    _HAVE_ADELIC = True
except Exception:  # pragma: no cover
    _HAVE_ADELIC = False

TOL = 1e-9
_TWO_PI = 2.0 * np.pi
PAULI = (
    np.array([[0, 1], [1, 0]], dtype=complex),
    np.array([[0, -1j], [1j, 0]], dtype=complex),
    np.array([[1, 0], [0, -1]], dtype=complex),
)


def _sieve(n):
    """Primes up to n (Sieve of Eratosthenes) -- fallback if adelic is absent."""
    flag = [True] * (n + 1)
    out = []
    for p in range(2, n + 1):
        if flag[p]:
            out.append(p)
            for k in range(p * p, n + 1, p):
                flag[k] = False
    return out


def selected_per_node_diagonal(n):
    """Return a selected node-distinct diagonal.

    The historical function name is preserved for imports. ``log(p)`` or
    ``diag(1..n)`` is supplied externally; this helper does not establish an RH
    role or a nodal derivation for the diagonal.
    """
    if _HAVE_ADELIC:
        eng = AdelicDynamics(max_prime=max(15, 4 * n))
        nu = np.asarray(eng.nu_f, dtype=float)[:n]
        if nu.size == n:
            return np.diag(nu), "nu_f = log p (adelic-model input)"
    primes = np.array(_sieve(max(15, 4 * n)), dtype=float)[:n]
    if primes.size == n:
        return np.diag(np.log(primes)), "nu_f = log p (sieve input)"
    return np.diag(np.arange(1, n + 1, dtype=float)), "diag(1..n) (selected input)"


# Historical API alias; the old name does not confer canonical status.
canonical_per_node_diagonal = selected_per_node_diagonal


# --------------------------------------------------------------------------- #
# Graph matrices used by the finite model. On the regular graphs below,
# ``D-A`` is a scalar multiple of the random-walk Laplacian and therefore has
# the same eigenspaces.
# --------------------------------------------------------------------------- #
def adjacency_laplacian(G, nodes):
    """Return adjacency and combinatorial-Laplacian matrices for ``G``."""
    A = nx.to_numpy_array(G, nodelist=nodes)
    L = np.diag(A.sum(axis=1)) - A
    return A, L


def _matrix_function(S, f):
    """Apply scalar f to a symmetric matrix S via its spectral decomposition."""
    w, V = np.linalg.eigh(S)
    return (V * f(w)) @ V.T


def catalog_operators(A, L):
    """Return the four matrices tested by this finite benchmark.

    This explicit set is not asserted to equal the engine's reachable set.
    ``exp(-L/2)`` is an ordinary finite matrix heat kernel here.
    """
    return {
        "A": A,
        "L = D - A": L,
        "L^2": L @ L,
        "exp(-L/2)": _matrix_function(L, lambda x: np.exp(-0.5 * x)),
    }


def commutator_norm(M, N):
    """Frobenius norm ||M N - N M||."""
    return float(np.linalg.norm(M @ N - N @ M))


def symmetric_projector(mats):
    """Reynolds projector Pi = (1/|G|) sum_g P_g onto Fix(G) = V^G."""
    n = mats[0].shape[0]
    Pi = np.zeros((n, n))
    for M in mats:
        Pi += M
    return Pi / len(mats)


# --------------------------------------------------------------------------- #
# External tensor-product model: lifts, group actions and matrix products
# --------------------------------------------------------------------------- #
def su2_generators():
    """T_a = sigma_a / 2: traceless, Hermitian, [T_a, T_b] = i eps_abc T_c."""
    return [p / 2.0 for p in PAULI]


def su2_element(theta, axis):
    """exp(i theta n.sigma) = cos(theta) I + i sin(theta) n.sigma (genuine SU(2))."""
    axis = np.asarray(axis, dtype=float)
    axis = axis / np.linalg.norm(axis)
    nsig = axis[0] * PAULI[0] + axis[1] * PAULI[1] + axis[2] * PAULI[2]
    return np.cos(theta) * np.eye(2, dtype=complex) + 1j * np.sin(theta) * nsig


def u1_element(phase):
    """Scalar phase matrix chosen for the external U(1) comparison."""
    return np.exp(1j * phase) * np.eye(2, dtype=complex)


def color_scalar_part(X):
    """Projection of a colour matrix onto the scalar-matrix subspace."""
    d = X.shape[0]
    return (np.trace(X) / d) * np.eye(d, dtype=complex)


def lift(M, d):
    """Embed ``M`` as ``M (x) I_d`` in the selected tensor-product model."""
    return np.kron(M, np.eye(d, dtype=complex))


def gauge_transform(n_sites, U):
    """Embed ``U`` as the fibre action ``I_V (x) U``."""
    return np.kron(np.eye(n_sites), U)


def holonomy(edge_unitaries):
    """Ordered product of edge unitaries around a loop."""
    H = np.eye(edge_unitaries[0].shape[0], dtype=complex)
    for U in edge_unitaries:
        H = H @ U
    return H


# --------------------------------------------------------------------------- #
# TEST 1 -- selected graph matrices and an invariant seed
# --------------------------------------------------------------------------- #
def test_graph_commutant_subset():
    print("=" * 78)
    print("TEST 1 -- COMPLETE-GRAPH MODEL: tested matrices commute with S_5")
    print("=" * 78)
    G = nx.complete_graph(5)  # every vertex permutation is an automorphism
    nodes = list(G.nodes())
    n = len(nodes)
    mats = automorphism_matrices(G, nodes)
    A, L = adjacency_laplacian(G, nodes)
    ops = catalog_operators(A, L)
    eye = np.eye(n)

    # (a) The tested finite family is a subset of the commutant.
    e_worst = max(commutator_norm(M, P) for M in ops.values() for P in mats)
    # (b) Check preservation from one invariant seed and one synthetic target.
    Pi = symmetric_projector(mats)
    fix_dim = int(round(np.trace(Pi)))
    perp_dim = n - fix_dim
    v = Pi @ eye[:, 0]
    leak = max(float(np.linalg.norm((eye - Pi) @ (M @ v))) for M in ops.values())
    w = (eye - Pi) @ eye[:, 0]
    overlap = max(abs(float(w @ (M @ v))) for M in ops.values())

    ok = e_worst < TOL and leak < TOL and overlap < TOL and perp_dim >= 1
    print(f"  (a) tested subset commutator: max ||[M, P_s]|| = {e_worst:.2e}")
    print(
        f"  (b) V split              : dim Fix(S_n) = {fix_dim}, "
        f"dim Fix(S_n)^perp = {perp_dim}"
    )
    print(f"      symmetric seed leak  : max ||(I-Pi) M v|| = {leak:.2e}")
    print(f"      synthetic overlap    : max |<w, M v>| = {overlap:.2e}")
    print(
        f"  VERDICT: {'PASS' if ok else 'FAIL'} -- this four-matrix family is a "
        "commuting subset and preserves the selected invariant seed"
    )
    print("  No map from S(T) to the synthetic vector w is defined or inferred.")
    print()
    return ok


# --------------------------------------------------------------------------- #
# TEST 2 -- external tensor-product and U(1)/SU(2) matrix comparison
# --------------------------------------------------------------------------- #
def test_external_colour_model():
    print("=" * 78)
    print("TEST 2 -- EXTERNAL COLOUR MODEL: scalar lifts and selected U(1)/SU(2)")
    print("=" * 78)
    d = 2  # SU(2) colour
    G = nx.cycle_graph(6)
    nodes = list(G.nodes())
    n = len(nodes)
    A, L = adjacency_laplacian(G, nodes)
    ops = catalog_operators(A, L)
    rng = np.random.default_rng(7)

    # (a) Sample the tensor-factor commutation identity.
    g_worst = 0.0
    for _ in range(8):
        U = su2_element(rng.uniform(0.3, 1.2), rng.normal(size=3))
        Ug = gauge_transform(n, U)
        g_worst = max(
            g_worst, max(commutator_norm(lift(M, d), Ug) for M in ops.values())
        )
    ab_ok = g_worst < TOL

    # (b) Scalar U(1) matrices commute.
    u1_P = holonomy([u1_element(rng.uniform(0, _TWO_PI)) for _ in range(4)])
    u1_Q = holonomy([u1_element(rng.uniform(0, _TWO_PI)) for _ in range(4)])
    u1_comm = commutator_norm(u1_P, u1_Q)
    u1_ok = u1_comm < TOL

    # (c) Selected SU(2) products need not commute; [T_a,T_b] is traceless.
    tx, ty, _tz = su2_generators()
    f_na = tx @ ty - ty @ tx
    f_scalar = color_scalar_part(f_na)
    f_norm = float(np.linalg.norm(f_na))
    f_reach = float(np.linalg.norm(f_scalar))
    su2_P = holonomy(
        [su2_element(rng.uniform(0.3, 1.2), rng.normal(size=3)) for _ in range(4)]
    )
    su2_Q = holonomy(
        [su2_element(rng.uniform(0.3, 1.2), rng.normal(size=3)) for _ in range(4)]
    )
    su2_comm = commutator_norm(su2_P, su2_Q)
    na_ok = su2_comm > 1e-3 and f_norm > 1e-3 and f_reach < TOL

    ok = ab_ok and u1_ok and na_ok
    print(
        f"  (a) sampled tensor-factor commutators  : "
        f"max ||[f(A,L)(x)I, I(x)U]|| = {g_worst:.2e}  (8 sampled U)"
    )
    print(
        f"  (b) selected U(1) products commute     : " f"||[H_P, H_Q]|| = {u1_comm:.2e}"
    )
    print(
        f"  (c) SU(2) holonomies do NOT commute   : " f"||[H_P, H_Q]|| = {su2_comm:.3f}"
    )
    print(
        f"      generator commutator [T_a,T_b]    : ||F|| = {f_norm:.3f}, "
        f"scalar part ||P(F)|| = {f_reach:.2e}"
    )
    print(
        f"  VERDICT: {'PASS' if ok else 'FAIL'} -- the selected scalar and "
        "traceless matrix relations hold in this finite colour model"
    )
    print("  The model does not define a TNFR-derived gauge or a mass-gap target.")
    print()
    return ok


# --------------------------------------------------------------------------- #
# TEST 3 -- compare two independently defined orthogonal decompositions
# --------------------------------------------------------------------------- #
def test_independent_orthogonal_splits():
    print("=" * 78)
    print("TEST 3 -- TWO FINITE ORTHOGONAL SPLITS ON DIFFERENT SPACES")
    print("=" * 78)
    # Complete-graph permutation representation.
    G = nx.complete_graph(5)
    nodes = list(G.nodes())
    n = len(nodes)
    mats = automorphism_matrices(G, nodes)
    Pi = symmetric_projector(mats)
    graph_fix = int(round(np.trace(Pi)))
    graph_rest = n - graph_fix
    graph_orth = float(np.linalg.norm(Pi @ (np.eye(n) - Pi)))

    # Independently chosen scalar/traceless Hermitian matrix split.
    d = 2
    tx, ty, tz = su2_generators()
    rest_basis = [tx, ty, tz]  # su(2): traceless, dim d^2 - 1 = 3
    i_d = np.eye(d, dtype=complex)
    scalar_dim = 1
    traceless_dim = d * d - scalar_dim
    colour_orth = max(abs(complex(np.trace(i_d.conj().T @ T))) for T in rest_basis)

    ok = (
        graph_rest >= 1
        and traceless_dim >= 1
        and graph_orth < TOL
        and colour_orth < TOL
        and traceless_dim == len(rest_basis)
    )
    print("  graph model (S_n)    : V        = Fix(S_n) (+) orthogonal rest")
    print(
        f"                         dims     = {graph_fix} (+) {graph_rest}   "
        f"block orthogonality ||Pi(I-Pi)|| = {graph_orth:.2e}"
    )
    print("  colour model (SU(2)) : Hermitian = scalar I (+) traceless basis")
    print(
        f"                         dims     = {scalar_dim} (+) {traceless_dim}   "
        f"block orthogonality max|<I_d, T_a>| = {colour_orth:.2e}"
    )
    print(
        f"  VERDICT: {'PASS' if ok else 'FAIL'} -- both selected finite splits "
        "are orthogonal, but they act on different spaces"
    )
    print("  This shared algebraic pattern does not identify RH and YM obstructions.")
    print()
    return ok


# --------------------------------------------------------------------------- #
# TEST 4 -- compare two externally supplied symmetry-breaking ingredients
# --------------------------------------------------------------------------- #
def test_supplied_symmetry_breakers():
    print("=" * 78)
    print("TEST 4 -- SELECTED DIAGONAL AND SU(2) GENERATORS: INDEPENDENT INPUTS")
    print("=" * 78)
    # A supplied node-distinct diagonal breaks the complete-graph symmetry.
    G = nx.complete_graph(5)
    nodes = list(G.nodes())
    n = len(nodes)
    mats = automorphism_matrices(G, nodes)
    Pi = symmetric_projector(mats)
    eye = np.eye(n)
    # The diagonal comes from an optional adelic model or a local fallback.
    p2, p2_label = selected_per_node_diagonal(n)
    graph_break = max(commutator_norm(p2, P) for P in mats)
    graph_leak = float(np.linalg.norm((eye - Pi) @ (p2 @ (Pi @ eye[:, 0]))))
    graph_breaks_symmetry = graph_break > 1e-3 and graph_leak > 1e-3

    # Independently supplied SU(2) generators do not commute.
    tx, ty, _tz = su2_generators()
    fibre_commutator = commutator_norm(tx, ty)
    fibre_is_noncommutative = fibre_commutator > 1e-3

    ok = graph_breaks_symmetry and fibre_is_noncommutative
    print(f"  graph input : D = diag({p2_label})")
    print(
        f"                breaks S_n (||[D,P_s]|| = {graph_break:.2f}, leak = "
        f"{graph_leak:.2f})"
    )
    print(
        "  colour input: non-commuting generators [T_x,T_y] != 0 "
        f"(||[T_x,T_y]|| = {fibre_commutator:.3f}) leave the scalar subspace"
    )
    print(
        f"  VERDICT: {'PASS' if ok else 'FAIL'} -- both supplied matrices break "
        "their respective finite symmetries"
    )
    print("  No common origin or shared Millennium-problem role is established.")
    print()
    return ok


# Historical callable names remain aliases for compatibility only.
test_rh_commutant_wall = test_graph_commutant_subset
test_ym_commutant_wall = test_external_colour_model
test_one_shape_two_groups = test_independent_orthogonal_splits
test_nonderivable_escape_contrast = test_supplied_symmetry_breakers


def main():
    print(__doc__)
    t1 = test_graph_commutant_subset()
    t2 = test_external_colour_model()
    t3 = test_independent_orthogonal_splits()
    t4 = test_supplied_symmetry_breakers()

    print("=" * 78)
    print("SUMMARY")
    print("=" * 78)
    print(f"  TEST 1 graph-family commutators          : {'PASS' if t1 else 'FAIL'}")
    print(f"  TEST 2 external colour-matrix model      : {'PASS' if t2 else 'FAIL'}")
    print(f"  TEST 3 two independent orthogonal splits : {'PASS' if t3 else 'FAIL'}")
    print(f"  TEST 4 two supplied symmetry breakers    : {'PASS' if t4 else 'FAIL'}")
    structural = t1 and t2 and t3 and t4
    print()
    print(f"  STRUCTURAL CHECKS: {'ALL PASS' if structural else 'SOME FAIL'}")
    print()
    print("  SCOPE VERDICT: the listed matrices satisfy the reported finite algebra.")
    print("  The tested family is only a subset of each commutant. No engine")
    print("  reachability theorem, map for S(T), TNFR-derived gauge, or identification")
    print("  between RH and Yang--Mills obstructions is supplied.")
    return 0 if structural else 1


if __name__ == "__main__":
    raise SystemExit(main())