"""Tests for auxiliary U(1) field-coordinate diagnostics.

The specified rotations preserve selected coordinate norms and transform other
fields covariantly. The supplied connection A = d(arg Psi) is an exact vertex-
argument difference, with analytic cycle closure; the Yang-Mills-named outputs
measure numerical closure residuals of this construction. Tests also cover the
configured interaction-regime readouts on fresh seeded graph fixtures.
These algebraic checks do not establish that the coordinate rotations are
realizable as nodal operators, that their flow is autonomous or that the
regimes identify physical forces.
"""

from __future__ import annotations

import math
import os
import sys

import networkx as nx
import numpy as np
import pytest

sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "..", "src"))

from tests.diagnostic_graph_fixtures import (
    make_diagnostic_field_graph as _make_tnfr_graph,
)
from tnfr.constants import inject_defaults
from tnfr.physics.canonical import compute_phase_curvature, compute_phase_gradient
from tnfr.physics.extended import compute_dnfr_flux, compute_phase_current
from tnfr.physics.gauge import (
    GAUGE_CLOSURE_TOLERANCE,
    N_REGIMES,
    REGIME_ACTIVITY_SHARE,
    BianchiIdentityResult,
    YangMillsFieldEquations,
    apply_gauge_transformation,
    capture_gauge_snapshot,
    classify_interaction_regime,
    classify_interaction_regime_formal,
    classify_network_regimes,
    compute_chirality_norm,
    compute_covariant_derivative,
    compute_covariant_derivative_magnitude,
    compute_dual_chirality,
    compute_dual_topological_charge,
    compute_gauge_connection,
    compute_gauge_coupling_constant,
    compute_gauge_curvature,
    compute_gauge_energy_decomposition,
    compute_gauss_law_residual,
    compute_matter_current,
    compute_network_interaction_profile,
    compute_topological_norm,
    compute_yang_mills_action,
    compute_yang_mills_equations,
    verify_bianchi_identity,
    verify_gauge_invariance,
)
from tnfr.physics.unified import (
    compute_chirality_field,
    compute_complex_geometric_field,
    compute_topological_charge,
)

# ---------------------------------------------------------------------------
# Fixtures
# ---------------------------------------------------------------------------


@pytest.fixture
def ws_graph():
    return _make_tnfr_graph(30, "watts_strogatz", seed=42)


@pytest.fixture
def random_alpha(ws_graph):
    """Random gauge parameters for each node."""
    rng = np.random.default_rng(123)
    return {n: float(rng.uniform(0, 2 * math.pi)) for n in ws_graph.nodes()}


# ===================================================================
# 1. Gauge Transformation Mechanics
# ===================================================================


@pytest.fixture(scope="module")
def invariance_report():
    graph = _make_tnfr_graph(30, seed=42)
    rng = np.random.default_rng(123)
    angles = {node: float(rng.uniform(0, 2 * math.pi)) for node in graph}
    return verify_gauge_invariance(graph, angles)


@pytest.fixture(scope="module")
def field_equations():
    return compute_yang_mills_equations(_make_tnfr_graph(30, seed=42))


@pytest.fixture(scope="module")
def interaction_profile():
    return compute_network_interaction_profile(_make_tnfr_graph(30, seed=42))


@pytest.fixture(scope="module")
def formal_regimes():
    graph = _make_tnfr_graph(30, seed=42)
    return {node: classify_interaction_regime_formal(graph, node) for node in graph}


class TestGaugeTransformation:
    """Validate the rotation Ψ → e^{iα}Ψ produces correct K_φ'/J_φ'."""

    def test_zero_alpha_identity(self, ws_graph):
        """α = 0 everywhere → fields unchanged."""
        alpha_zero = {n: 0.0 for n in ws_graph.nodes()}
        result = apply_gauge_transformation(ws_graph, alpha_zero)

        k_phi = compute_phase_curvature(ws_graph)
        j_phi = compute_phase_current(ws_graph)

        for n in ws_graph.nodes():
            assert abs(result["k_phi"][n] - k_phi[n]) < 1e-12
            assert abs(result["j_phi"][n] - j_phi[n]) < 1e-12

    def test_pi_half_rotation(self, ws_graph):
        """α = π/2 → K_φ' = −J_φ, J_φ' = K_φ."""
        alpha = {n: math.pi / 2 for n in ws_graph.nodes()}
        result = apply_gauge_transformation(ws_graph, alpha)

        k_phi = compute_phase_curvature(ws_graph)
        j_phi = compute_phase_current(ws_graph)

        for n in ws_graph.nodes():
            assert abs(result["k_phi"][n] - (-j_phi[n])) < 1e-10
            assert abs(result["j_phi"][n] - k_phi[n]) < 1e-10

    def test_pi_rotation(self, ws_graph):
        """α = π → K_φ' = −K_φ, J_φ' = −J_φ (sign flip)."""
        alpha = {n: math.pi for n in ws_graph.nodes()}
        result = apply_gauge_transformation(ws_graph, alpha)

        k_phi = compute_phase_curvature(ws_graph)
        j_phi = compute_phase_current(ws_graph)

        for n in ws_graph.nodes():
            assert abs(result["k_phi"][n] + k_phi[n]) < 1e-10
            assert abs(result["j_phi"][n] + j_phi[n]) < 1e-10

    def test_psi_magnitude_preserved(self, ws_graph, random_alpha):
        """Local U(1) preserves |Ψ| at every node."""
        psi_before = compute_complex_geometric_field(ws_graph)
        result = apply_gauge_transformation(ws_graph, random_alpha)

        for n in ws_graph.nodes():
            assert abs(abs(result["psi"][n]) - abs(psi_before[n])) < 1e-10

    @pytest.mark.parametrize(
        "alpha",
        [None, {0: True}, {0: math.nan}, {0: math.inf}, {0: "0.1"}],
    )
    def test_transformation_rejects_invalid_angle_inputs(self, alpha):
        """Rotations reject invalid mappings and non-finite angles."""
        graph = nx.path_graph(2)
        inject_defaults(graph)
        value = [] if alpha is None else alpha
        with pytest.raises((TypeError, ValueError)):
            apply_gauge_transformation(graph, value)

    def test_double_transform_composes(self, ws_graph):
        """Two successive transforms α₁, α₂ equal one transform α₁ + α₂."""
        rng = np.random.default_rng(77)
        alpha1 = {n: float(rng.uniform(0, math.pi)) for n in ws_graph.nodes()}
        alpha2 = {n: float(rng.uniform(0, math.pi)) for n in ws_graph.nodes()}
        alpha_sum = {n: alpha1[n] + alpha2[n] for n in ws_graph.nodes()}

        # Apply sum in one step
        r_sum = apply_gauge_transformation(ws_graph, alpha_sum)

        # Apply sequentially: first α₁, then α₂ on transformed k_phi/j_phi
        r1 = apply_gauge_transformation(ws_graph, alpha1)

        # Manually apply second rotation on r1 results
        for n in ws_graph.nodes():
            a2 = alpha2[n]
            kp = r1["k_phi"][n]
            jp = r1["j_phi"][n]
            kp2 = kp * math.cos(a2) - jp * math.sin(a2)
            jp2 = kp * math.sin(a2) + jp * math.cos(a2)
            assert abs(kp2 - r_sum["k_phi"][n]) < 1e-10
            assert abs(jp2 - r_sum["j_phi"][n]) < 1e-10


# ===================================================================
# 2. Gauge Invariance of Physical Quantities
# ===================================================================


class TestGaugeInvariance:
    """Validate the invariant norms and record variant legacy quantities."""

    def test_energy_density_invariant(self, invariance_report):
        """Energy density ℰ(i) exactly invariant under local U(1)."""
        result = invariance_report
        assert result.energy_max_deviation < 1e-10

    def test_magnitude_invariant(self, invariance_report):
        """|Ψ(i)| exactly invariant under local U(1)."""
        result = invariance_report
        assert result.magnitude_max_deviation < 1e-10

    def test_topological_norm_invariant(self, invariance_report):
        """|𝒯|² = 𝒬² + 𝒬̃² invariant under local U(1)."""
        result = invariance_report
        assert result.topological_norm_max_deviation < 1e-10

    def test_chirality_norm_invariant(self, invariance_report):
        """|𝒳|² = χ² + χ̃² invariant under local U(1)."""
        result = invariance_report
        assert result.chirality_norm_max_deviation < 1e-10

    def test_symmetry_breaking_NOT_invariant(self, invariance_report):
        """𝒮 = (|∇φ|² − K_φ²) + (J_φ² − J_ΔNFR²) is NOT gauge-invariant.

        K_φ² and J_φ² individually change under rotation even though
        K_φ² + J_φ² = |Ψ|² is preserved.  For non-trivial α the
        deviation must be non-zero.
        """
        result = invariance_report
        assert (
            result.symmetry_breaking_max_deviation > 1e-6
        ), "Expected 𝒮 to change under gauge transform"

    def test_coherence_invariant(self, invariance_report):
        """C(t) invariant (external to Ψ internal rotation)."""
        result = invariance_report
        assert result.coherence_deviation < 1e-12

    def test_noether_charge_NOT_invariant(self, invariance_report):
        """Q = Σ(Φ_s + K_φ) is NOT gauge-invariant (K_φ rotates)."""
        result = invariance_report
        assert result.details["has_nontrivial_alpha"]
        # Q should change, but all *true* gauge-invariant quantities hold
        assert (
            result.noether_charge_deviation > 1e-6
        ), "Expected Q to change under gauge transform"
        assert result.is_invariant

    def test_verify_supports_signed_seed(self, ws_graph):
        """Signed seeds follow the repository deterministic convention."""
        first = verify_gauge_invariance(ws_graph, seed=-42)
        second = verify_gauge_invariance(ws_graph, seed=-42)
        assert first == second

    @pytest.mark.parametrize("tolerance", [0.0, -1.0, math.nan, math.inf, True])
    def test_verify_rejects_invalid_tolerance(self, ws_graph, tolerance):
        """An invariance threshold must be finite and strictly positive."""
        with pytest.raises((TypeError, ValueError)):
            verify_gauge_invariance(ws_graph, tolerance=tolerance)

    def test_invariance_factorisation(self, ws_graph):
        """Verify |𝒯|² = |Ψ|²·|Ω|² factorisation identity.

        For any node: 𝒬² + 𝒬̃² = |Ψ|² · (|∇φ|² + J_ΔNFR²).
        """
        psi = compute_complex_geometric_field(ws_graph)
        grad_phi = compute_phase_gradient(ws_graph)
        j_dnfr = compute_dnfr_flux(ws_graph)
        topo_norm = compute_topological_norm(ws_graph)

        for n in ws_graph.nodes():
            psi_sq = abs(psi[n]) ** 2
            omega_sq = grad_phi.get(n, 0.0) ** 2 + j_dnfr.get(n, 0.0) ** 2
            expected = psi_sq * omega_sq
            assert (
                abs(topo_norm[n] - expected) < 1e-10
            ), f"Node {n}: |𝒯|²={topo_norm[n]:.6e} ≠ |Ψ|²·|Ω|²={expected:.6e}"


# ===================================================================
# 3. Multi-Topology Validation
# ===================================================================


class TestMultiTopology:
    """Gauge invariance holds across different network topologies."""

    @pytest.mark.parametrize("topology", ["watts_strogatz", "barabasi_albert", "grid"])
    def test_gauge_invariance_cross_topology(self, topology):
        """ℰ invariance verified on WS, BA, and Grid graphs."""
        n = 25 if topology == "grid" else 30
        G = _make_tnfr_graph(n, topology, seed=42)
        result = verify_gauge_invariance(G, seed=99)
        assert result.is_invariant
        assert result.energy_max_deviation < 1e-10


# ===================================================================
# 4. Gauge Connection
# ===================================================================


class TestGaugeConnection:
    """Validate the exact connection A_ij = d(arg Ψ)_ij."""

    def test_antisymmetry(self, ws_graph):
        """A_ji = −A_ij for undirected graphs."""
        conn = compute_gauge_connection(ws_graph)
        for u, v in ws_graph.edges():
            a_uv = conn.get((u, v), 0.0)
            a_vu = conn.get((v, u), 0.0)
            assert abs(a_uv + a_vu) < 1e-12, f"A({u},{v}) + A({v},{u}) ≠ 0"

    def test_connection_range(self, ws_graph):
        """Oriented representatives lie in [−π, π], including stored reverses."""
        conn = compute_gauge_connection(ws_graph)
        for edge, val in conn.items():
            assert -math.pi <= val <= math.pi, f"A{edge} = {val} outside [−π, π]"

    def test_connection_zero_for_uniform_psi(self):
        """If Ψ has the same phase everywhere, A_ij = 0."""
        G = nx.complete_graph(5)
        inject_defaults(G)
        # Set all phases and delta_nfr identically
        for n in G.nodes():
            G.nodes[n]["phase"] = 1.0
            G.nodes[n]["delta_nfr"] = 0.3
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        conn = compute_gauge_connection(G)
        for val in conn.values():
            assert abs(val) < 1e-6, f"Non-zero connection {val} for uniform Ψ"

    def test_half_turn_preserves_oriented_antisymmetry(self, monkeypatch):
        monkeypatch.setattr(
            "tnfr.physics.gauge.compute_complex_geometric_field",
            lambda graph: {0: 1.0 + 0j, 1: -1.0 + 0j},
        )
        connection = compute_gauge_connection(nx.path_graph(2))
        assert connection[0, 1] == -math.pi
        assert connection[1, 0] == math.pi


# ===================================================================
# 5. Covariant Derivative
# ===================================================================


class TestCovariantDerivative:
    """Validate D_ij Ψ = Ψ(j) − e^{iA_ij}Ψ(i)."""

    def test_magnitude_gauge_invariant(self, ws_graph, random_alpha, monkeypatch):
        """Rebuild the connection after independent local field rotations."""
        psi = compute_complex_geometric_field(ws_graph)
        before = compute_covariant_derivative(ws_graph)
        rotations = {
            node: complex(math.cos(alpha), math.sin(alpha))
            for node, alpha in random_alpha.items()
        }
        rotated = {node: rotations[node] * value for node, value in psi.items()}
        monkeypatch.setattr(
            "tnfr.physics.gauge.compute_complex_geometric_field", lambda graph: rotated
        )
        after = compute_covariant_derivative(ws_graph)
        for (u, v), value in before.items():
            assert after[u, v] == pytest.approx(rotations[v] * value, abs=1e-12)
            assert abs(after[u, v]) == pytest.approx(abs(abs(psi[v]) - abs(psi[u])))

    @pytest.mark.parametrize("exponent", [-60, -500, -1022])
    @pytest.mark.parametrize("target_amplitude", [1.0, 2.0])
    def test_small_nonzero_fields_retain_phase(
        self, monkeypatch, exponent, target_amplitude
    ):
        """Compare scaled differences so absolute tolerances cannot hide lost phase."""
        graph = nx.path_graph(2)
        scale = math.ldexp(1.0, exponent)
        field = {0: complex(scale, 0.0), 1: complex(0.0, target_amplitude * scale)}
        monkeypatch.setattr(
            "tnfr.physics.gauge.compute_complex_geometric_field", lambda graph: field
        )
        assert compute_gauge_connection(graph)[0, 1] == pytest.approx(math.pi / 2)
        before = compute_covariant_derivative(graph)
        expected = target_amplitude - 1.0
        assert before[0, 1] / scale == pytest.approx(complex(0.0, expected), abs=1e-15)

        # Local rotations by +pi/2 and -pi/2; reconstruct A rather than rotating D.
        field = {0: complex(0.0, scale), 1: complex(target_amplitude * scale, 0.0)}
        after = compute_covariant_derivative(graph)
        assert after[0, 1] / scale == pytest.approx(complex(expected, 0.0), abs=1e-15)
        assert after[0, 1] / scale == pytest.approx(
            -1j * before[0, 1] / scale, abs=1e-15
        )

    @pytest.mark.parametrize("zero_node", [0, 1])
    def test_exact_zero_keeps_deterministic_phase(self, monkeypatch, zero_node):
        graph = nx.path_graph(2)
        scale = math.ldexp(1.0, -60)
        field = {zero_node: complex(-0.0, -0.0), 1 - zero_node: complex(0.0, scale)}
        monkeypatch.setattr(
            "tnfr.physics.gauge.compute_complex_geometric_field", lambda graph: field
        )
        expected_angle = math.pi / 2 if zero_node == 0 else -math.pi / 2
        assert compute_gauge_connection(graph)[0, 1] == pytest.approx(expected_angle)
        for magnitude in compute_covariant_derivative_magnitude(graph).values():
            assert magnitude / scale == pytest.approx(1.0)

    def test_zero_for_parallel_transport(self):
        """If Ψ is parallel-transported (constant magnitude, phase follows
        connection), then D_ij Ψ ≈ 0.

        For a pair of connected nodes with identical Ψ values,
        D = Ψ(j) − e^{i·0}·Ψ(i) = 0.
        """
        G = nx.path_graph(3)
        inject_defaults(G)
        # Set identical phases and delta_nfr to get identical Ψ
        for n in G.nodes():
            G.nodes[n]["phase"] = 2.0
            G.nodes[n]["delta_nfr"] = 0.4
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        cov_mag = compute_covariant_derivative_magnitude(G)
        for val in cov_mag.values():
            assert val < 1e-6, f"Non-zero |D| = {val} for parallel Ψ"


# ===================================================================
# 6. Gauge Curvature (Field Strength)
# ===================================================================


class TestGaugeCurvature:
    """Validate F_C = Σ_C A_ij closes for the exact connection."""

    def test_curvature_on_complete_graph(self):
        """Complete graph K_4 has many triangles; curvature computed."""
        G = nx.complete_graph(5)
        inject_defaults(G)
        rng = np.random.default_rng(42)
        for n in G.nodes():
            G.nodes[n]["phase"] = float(rng.uniform(0, 2 * math.pi))
            G.nodes[n]["delta_nfr"] = float(rng.uniform(-0.5, 0.5))
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        curv = compute_gauge_curvature(G)
        assert len(curv) > 0, "K_5 should have triangles"
        # A=d(arg Ψ), so every cycle closes modulo floating-point error.
        for cycle, f in curv.items():
            assert abs(f) <= GAUGE_CLOSURE_TOLERANCE

    @pytest.mark.parametrize("topology", ["watts_strogatz", "barabasi_albert", "grid"])
    def test_canonical_connection_is_pure_gauge_across_topologies(self, topology):
        """Every detected cycle closes for arbitrary canonical snapshots."""
        G = _make_tnfr_graph(25, topology, seed=171)
        curv = compute_gauge_curvature(G)
        assert all(abs(value) <= GAUGE_CLOSURE_TOLERANCE for value in curv.values())

    def test_four_cycle_with_heterogeneous_labels_keeps_traversal(self):
        """Cycle detection neither compares labels nor sorts away edge order."""
        labels = [0, "one", (2,), frozenset({3})]
        G = nx.cycle_graph(labels)
        inject_defaults(G)
        for index, node in enumerate(labels):
            G.nodes[node]["phase"] = 0.7 * index
            G.nodes[node]["delta_nfr"] = 0.2 * (index - 1)
            G.nodes[node]["frequency"] = 1.0
            G.nodes[node]["EPI"] = "epi"

        curv = compute_gauge_curvature(G)
        assert len(curv) == 1
        cycle = next(iter(curv))
        assert all(
            G.has_edge(cycle[index], cycle[(index + 1) % len(cycle)])
            for index in range(len(cycle))
        )
        assert abs(curv[cycle]) <= GAUGE_CLOSURE_TOLERANCE

    def test_curvature_zero_flat_connection(self):
        """Gauge-flat configuration: identical Ψ → F_C = 0 everywhere."""
        G = nx.complete_graph(4)
        inject_defaults(G)
        for n in G.nodes():
            G.nodes[n]["phase"] = 1.5
            G.nodes[n]["delta_nfr"] = 0.2
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        curv = compute_gauge_curvature(G)
        for cycle, f in curv.items():
            assert abs(f) < 1e-6, f"Non-zero F on {cycle} for flat connection"

    def test_yang_mills_action_is_only_closure_residual(self, ws_graph):
        """The legacy action is at roundoff for A=d(arg Ψ)."""
        s_ym = compute_yang_mills_action(ws_graph)
        n_cycles = len(compute_gauge_curvature(ws_graph))
        assert 0.0 <= s_ym <= 0.5 * n_cycles * GAUGE_CLOSURE_TOLERANCE**2

    @pytest.mark.parametrize("value", [2, 0, -1])
    def test_curvature_rejects_too_short_cycle_limit(self, value):
        """Cycle closure is defined only for cycles of length at least three."""
        with pytest.raises(ValueError):
            compute_gauge_curvature(nx.cycle_graph(4), max_cycle_length=value)

    @pytest.mark.parametrize("value", [3.0, True, "4"])
    def test_curvature_rejects_noninteger_cycle_limit(self, value):
        """Cycle limits are not silently coerced."""
        with pytest.raises(TypeError):
            compute_gauge_curvature(nx.cycle_graph(4), max_cycle_length=value)

    def test_yang_mills_zero_for_flat(self):
        """S_YM = 0 for a gauge-flat configuration."""
        G = nx.complete_graph(4)
        inject_defaults(G)
        for n in G.nodes():
            G.nodes[n]["phase"] = 1.0
            G.nodes[n]["delta_nfr"] = 0.3
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        assert compute_yang_mills_action(G) < 1e-10


# ===================================================================
# 7. Dual Charges
# ===================================================================


class TestDualCharges:
    """Validate dual topological charge 𝒬̃ and dual chirality χ̃."""

    def test_dual_topological_charge_definition(self, ws_graph):
        """𝒬̃ = K_φ·|∇φ| + J_φ·J_ΔNFR computed correctly."""
        grad_phi = compute_phase_gradient(ws_graph)
        k_phi = compute_phase_curvature(ws_graph)
        j_phi = compute_phase_current(ws_graph)
        j_dnfr = compute_dnfr_flux(ws_graph)

        q_dual = compute_dual_topological_charge(ws_graph)

        for n in ws_graph.nodes():
            expected = k_phi.get(n, 0.0) * grad_phi.get(n, 0.0) + j_phi.get(
                n, 0.0
            ) * j_dnfr.get(n, 0.0)
            assert abs(q_dual[n] - expected) < 1e-12

    def test_dual_chirality_definition(self, ws_graph):
        """χ̃ = |∇φ|·J_φ + K_φ·J_ΔNFR computed correctly."""
        grad_phi = compute_phase_gradient(ws_graph)
        k_phi = compute_phase_curvature(ws_graph)
        j_phi = compute_phase_current(ws_graph)
        j_dnfr = compute_dnfr_flux(ws_graph)

        chi_dual = compute_dual_chirality(ws_graph)

        for n in ws_graph.nodes():
            expected = grad_phi.get(n, 0.0) * j_phi.get(n, 0.0) + k_phi.get(
                n, 0.0
            ) * j_dnfr.get(n, 0.0)
            assert abs(chi_dual[n] - expected) < 1e-12

    def test_topological_norm_from_q_and_q_dual(self, ws_graph):
        """Verify |𝒯|² = 𝒬² + 𝒬̃² by explicit computation."""
        topo_charge = compute_topological_charge(ws_graph)
        q_dual = compute_dual_topological_charge(ws_graph)
        topo_norm = compute_topological_norm(ws_graph)

        for n in ws_graph.nodes():
            expected = topo_charge.get(n, 0.0) ** 2 + q_dual[n] ** 2
            assert abs(topo_norm[n] - expected) < 1e-10

    def test_chirality_norm_from_chi_and_chi_dual(self, ws_graph):
        """Verify |𝒳|² = χ² + χ̃² by explicit computation."""
        chirality = compute_chirality_field(ws_graph)
        chi_dual = compute_dual_chirality(ws_graph)
        chiral_norm = compute_chirality_norm(ws_graph)

        for n in ws_graph.nodes():
            expected = chirality.get(n, 0.0) ** 2 + chi_dual[n] ** 2
            assert abs(chiral_norm[n] - expected) < 1e-10


# ===================================================================
# 8. Gauge Snapshot
# ===================================================================


class TestGaugeSnapshot:
    """Validate the GaugeSnapshot capture."""

    def test_snapshot_fields_present(self, ws_graph):
        """All fields are populated in the snapshot."""
        snap = capture_gauge_snapshot(ws_graph)
        assert len(snap.psi) == len(ws_graph)
        assert len(snap.psi_magnitude) == len(ws_graph)
        assert len(snap.psi_phase) == len(ws_graph)
        assert len(snap.connection) > 0
        assert len(snap.energy_density) == len(ws_graph)
        assert len(snap.topological_norm) == len(ws_graph)
        assert len(snap.chirality_norm) == len(ws_graph)

    def test_snapshot_magnitude_consistency(self, ws_graph):
        """snapshot.psi_magnitude matches |snapshot.psi|."""
        snap = capture_gauge_snapshot(ws_graph)
        for n in ws_graph.nodes():
            assert abs(snap.psi_magnitude[n] - abs(snap.psi[n])) < 1e-12

    def test_snapshot_declares_pure_gauge_scope(self, ws_graph):
        """Dataclass semantics expose the exact-connection limitation."""
        snap = capture_gauge_snapshot(ws_graph)
        assert snap.canonical_connection_is_pure_gauge
        assert snap.curvature_is_numerical_residual


# ===================================================================
# 9. Energy Decomposition
# ===================================================================


class TestEnergyDecomposition:
    """Validate gauge-theoretic energy sector decomposition."""

    def test_sectors_sum_to_total(self, ws_graph):
        """potential + gradient + gauge + flux = 2 × total_energy."""
        decomp = compute_gauge_energy_decomposition(ws_graph)
        sector_sum = (
            decomp["potential_sector"]
            + decomp["gradient_sector"]
            + decomp["gauge_sector"]
            + decomp["flux_sector"]
        )
        assert abs(sector_sum - 2.0 * decomp["total_energy"]) < 1e-10

    def test_fractions_sum_to_one(self, ws_graph):
        """Sector fractions sum to 1 (within numerical precision)."""
        decomp = compute_gauge_energy_decomposition(ws_graph)
        frac_sum = (
            decomp["potential_fraction"]
            + decomp["gradient_fraction"]
            + decomp["gauge_fraction"]
            + decomp["flux_fraction"]
        )
        assert abs(frac_sum - 1.0) < 1e-6


# ===================================================================
# 10. Interaction Regime Classification
# ===================================================================


class TestInteractionRegimes:
    """Validate gauge-based interaction regime classification."""

    def test_regime_is_valid_string(self, ws_graph):
        """Regime classification returns a known regime type."""
        valid_regimes = {"em_like", "weak_like", "strong_like", "gravity_like"}
        for n in ws_graph.nodes():
            info = classify_interaction_regime(ws_graph, n)
            assert info["regime"] in valid_regimes

    def test_regime_scores_sum(self, ws_graph):
        """Regime scores should be non-negative."""
        for n in list(ws_graph.nodes())[:5]:
            info = classify_interaction_regime(ws_graph, n)
            for score in info["regime_scores"].values():
                assert score >= 0.0

    def test_network_regimes(self, ws_graph):
        """Full network classification returns valid structure."""
        result = classify_network_regimes(ws_graph)
        assert "per_node" in result
        assert "regime_distribution" in result
        assert "dominant_regime" in result
        assert "mean_gauge_curvature" in result
        assert "gauge_flatness" in result

        total = sum(result["regime_distribution"].values())
        assert total == len(ws_graph)

    def test_gauge_flatness_range(self, ws_graph):
        """Gauge flatness is a fraction in [0, 1]."""
        result = classify_network_regimes(ws_graph)
        assert 0.0 <= result["gauge_flatness"] <= 1.0


# ===================================================================
# 11. Reproducibility
# ===================================================================


class TestReproducibility:
    """Verify deterministic results under fixed seeds (Invariant #6)."""

    def test_gauge_invariance_reproducible(self):
        """Same graph + same seed → identical gauge invariance result."""
        G1 = _make_tnfr_graph(20, "watts_strogatz", seed=42)
        G2 = _make_tnfr_graph(20, "watts_strogatz", seed=42)

        r1 = verify_gauge_invariance(G1, seed=99)
        r2 = verify_gauge_invariance(G2, seed=99)

        assert r1.energy_max_deviation == r2.energy_max_deviation
        assert r1.noether_charge_deviation == r2.noether_charge_deviation
        assert r1.topological_norm_max_deviation == r2.topological_norm_max_deviation

    def test_snapshot_reproducible(self):
        """Same graph → same snapshot values."""
        G = _make_tnfr_graph(15, "watts_strogatz", seed=42)
        s1 = capture_gauge_snapshot(G)
        s2 = capture_gauge_snapshot(G)

        for n in G.nodes():
            assert s1.psi[n] == s2.psi[n]
            assert s1.energy_density[n] == s2.energy_density[n]


# ===================================================================
# 12. Edge Cases
# ===================================================================


class TestEdgeCases:
    """Test behaviour on minimal and degenerate graphs."""

    def test_pair_graph(self):
        """Two-node graph: gauge structure still valid."""
        G = nx.path_graph(2)
        inject_defaults(G)
        for n in G.nodes():
            G.nodes[n]["phase"] = float(n)
            G.nodes[n]["delta_nfr"] = 0.5
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        result = verify_gauge_invariance(G, seed=42)
        assert result.is_invariant

    def test_triangle_graph(self):
        """Triangle graph: single curvature plaquette."""
        G = nx.complete_graph(3)
        inject_defaults(G)
        rng = np.random.default_rng(42)
        for n in G.nodes():
            G.nodes[n]["phase"] = float(rng.uniform(0, 2 * math.pi))
            G.nodes[n]["delta_nfr"] = float(rng.uniform(-0.5, 0.5))
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        curv = compute_gauge_curvature(G)
        assert len(curv) == 1  # Exactly one triangle
        result = verify_gauge_invariance(G, seed=42)
        assert result.is_invariant

    def test_star_graph_no_triangles(self):
        """Star graph has no triangles → empty curvature."""
        G = nx.star_graph(5)
        inject_defaults(G)
        rng = np.random.default_rng(42)
        for n in G.nodes():
            G.nodes[n]["phase"] = float(rng.uniform(0, 2 * math.pi))
            G.nodes[n]["delta_nfr"] = float(rng.uniform(-0.5, 0.5))
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"

        curv = compute_gauge_curvature(G)
        assert len(curv) == 0  # No triangles in star graph
        assert compute_yang_mills_action(G) == 0.0


# ===================================================================
# Yang-Mills Field Equations & Formal Regime Classification Tests
# ===================================================================


class TestMatterCurrent:
    """Legacy invariant link-current residual J_matter(i,j)."""

    def test_antisymmetry(self, ws_graph):
        """J(j,i) = −J(i,j) for every oriented edge."""
        j_mat = compute_matter_current(ws_graph)
        for (u, v), val in j_mat.items():
            if (v, u) in j_mat:
                assert abs(j_mat[(v, u)] + val) < 1e-12

    def test_pure_gauge_link_current_is_numerical_zero(self, ws_graph):
        """The correctly signed invariant bilinear is real for A=d(arg Ψ)."""
        current = compute_matter_current(ws_graph)
        assert max(map(abs, current.values()), default=0.0) < 1e-12

    def test_uniform_psi_zero_current(self):
        """Spatially uniform Ψ → zero matter current."""
        G = nx.complete_graph(5)
        inject_defaults(G)
        for n in G.nodes():
            G.nodes[n]["phase"] = 0.0
            G.nodes[n]["delta_nfr"] = 0.5
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"
        j_mat = compute_matter_current(G)
        for val in j_mat.values():
            assert abs(val) < 1e-12

    def test_covers_all_edges(self, ws_graph):
        """Current dict contains both orientations for every edge."""
        j_mat = compute_matter_current(ws_graph)
        for u, v in ws_graph.edges():
            assert (u, v) in j_mat or (v, u) in j_mat


class TestYangMillsFieldEquations:
    """Legacy-named consistency residual on the pure-gauge surface."""

    def test_result_type(self, ws_graph, field_equations):
        """Returns YangMillsFieldEquations dataclass."""
        result = field_equations
        assert isinstance(result, YangMillsFieldEquations)
        assert result.canonical_connection_is_pure_gauge
        assert not result.is_dynamical_derivation

    def test_actions_non_negative(self, ws_graph, field_equations):
        """S_YM ≥ 0 and S_matter ≥ 0."""
        eq = field_equations
        assert eq.yang_mills_action >= 0.0
        assert eq.matter_action >= 0.0
        assert eq.total_action >= eq.yang_mills_action
        assert eq.total_action >= eq.matter_action

    def test_self_scale_is_zero_after_closure_tolerance(
        self, ws_graph, field_equations
    ):
        """No coupling is inferred from roundoff in an exact connection."""
        eq = field_equations
        assert eq.coupling_constant == 0.0
        assert eq.yang_mills_action == 0.0
        assert eq.max_residual < 1e-12

    def test_explicit_coupling(self, ws_graph):
        """User-specified diagnostic scale is preserved."""
        eq = compute_yang_mills_equations(ws_graph, coupling=1.0)
        assert eq.coupling_constant == pytest.approx(1.0, abs=1e-12)

    @pytest.mark.parametrize("coupling", [True, "1", -1.0, math.nan, math.inf])
    def test_invalid_explicit_coupling_rejected(self, ws_graph, coupling):
        """Invalid diagnostic scales cannot be silently floored or propagated."""
        with pytest.raises((TypeError, ValueError)):
            compute_yang_mills_equations(ws_graph, coupling=coupling)

    def test_residual_structure(self, ws_graph, field_equations):
        """Residuals are non-negative with correct mean/max."""
        eq = field_equations
        for r in eq.equation_residual.values():
            assert r >= 0.0
        assert eq.mean_residual >= 0.0
        assert eq.max_residual >= eq.mean_residual or len(eq.equation_residual) <= 1

    def test_matter_current_consistency(self, ws_graph, field_equations):
        """Matter current matches standalone computation."""
        eq = field_equations
        standalone = compute_matter_current(ws_graph)
        for edge, val in standalone.items():
            assert abs(eq.matter_current[edge] - val) < 1e-12

    def test_reproducibility(self):
        """Same seed → identical Yang-Mills equations."""
        G1 = _make_tnfr_graph(25, seed=99)
        G2 = _make_tnfr_graph(25, seed=99)
        eq1 = compute_yang_mills_equations(G1)
        eq2 = compute_yang_mills_equations(G2)
        assert eq1.yang_mills_action == pytest.approx(eq2.yang_mills_action, abs=1e-14)
        assert eq1.matter_action == pytest.approx(eq2.matter_action, abs=1e-14)
        assert eq1.coupling_constant == pytest.approx(eq2.coupling_constant, abs=1e-14)


class TestBianchiIdentity:
    """Legacy-named verification of exact-connection cycle closure."""

    def test_satisfied_ws(self, ws_graph):
        """All checked cycles close on a Watts-Strogatz graph."""
        result = verify_bianchi_identity(ws_graph)
        assert isinstance(result, BianchiIdentityResult)
        assert result.num_coboundaries_tested > 0
        assert result.num_cycles_tested == len(compute_gauge_curvature(ws_graph))
        assert result.max_residual <= GAUGE_CLOSURE_TOLERANCE
        assert result.is_satisfied

    def test_tree_graph_trivial(self):
        """Tree graph has no plaquettes → trivially satisfied."""
        G = nx.star_graph(6)
        inject_defaults(G)
        for n in G.nodes():
            G.nodes[n]["phase"] = 0.0
            G.nodes[n]["delta_nfr"] = 0.1
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"
        result = verify_bianchi_identity(G)
        assert result.is_satisfied
        assert result.num_coboundaries_tested == 0

    @pytest.mark.parametrize("tolerance", [True, "1e-8", -1.0, math.nan, math.inf])
    def test_invalid_tolerance_rejected(self, ws_graph, tolerance):
        """The closure comparison requires a finite non-negative tolerance."""
        with pytest.raises((TypeError, ValueError)):
            verify_bianchi_identity(ws_graph, tolerance=tolerance)


class TestGaussLawResidual:
    """Legacy-named divergence of the numerical link-current residual."""

    def test_residual_per_node(self, ws_graph):
        """Returns one residual per node, all non-negative."""
        residuals = compute_gauss_law_residual(ws_graph)
        assert len(residuals) == ws_graph.number_of_nodes()
        for val in residuals.values():
            assert val >= 0.0

    def test_uniform_psi_small_residual(self):
        """Spatially uniform Ψ → small Gauss residual."""
        G = nx.complete_graph(5)
        inject_defaults(G)
        for n in G.nodes():
            G.nodes[n]["phase"] = 0.0
            G.nodes[n]["delta_nfr"] = 0.5
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"
        residuals = compute_gauss_law_residual(G)
        for val in residuals.values():
            assert val < 1e-10


class TestGaugeCouplingConstant:
    """Legacy mean-squared closure statistic g² = ⟨F²⟩."""

    def test_tree_zero(self):
        """No plaquettes → g² = 0."""
        G = nx.star_graph(5)
        inject_defaults(G)
        for n in G.nodes():
            G.nodes[n]["phase"] = 0.0
            G.nodes[n]["delta_nfr"] = 0.1
            G.nodes[n]["frequency"] = 1.0
            G.nodes[n]["EPI"] = "epi"
        assert compute_gauge_coupling_constant(G) == 0.0


class TestRegimeActivityCriterion:
    """Legacy equal-share activity reporting convention."""

    def test_four_regimes(self):
        """Four retained labels imply a one-quarter reporting share."""
        assert N_REGIMES == 4
        assert REGIME_ACTIVITY_SHARE == pytest.approx(0.25, abs=1e-14)


class TestFormalInteractionRegimes:
    """Formalized historical four-label snapshot heuristic."""

    def test_order_params_bounded(self, formal_regimes):
        """All order parameters in [0, 1]."""
        for m in formal_regimes.values():
            assert 0.0 <= m.em_order_parameter <= 1.0 + 1e-12
            assert 0.0 <= m.weak_order_parameter <= 1.0 + 1e-12
            assert 0.0 <= m.strong_order_parameter  # can exceed 1 in principle
            assert 0.0 <= m.gravity_order_parameter <= 1.0 + 1e-12

    def test_strong_like_slot_does_not_promote_roundoff(self, formal_regimes):
        """Pure-gauge closure noise cannot masquerade as confinement."""
        for metrics in formal_regimes.values():
            assert metrics.strong_order_parameter == 0.0
            assert metrics.regime_scores["strong_like"] == 0.0

    def test_scores_sum_to_one(self, formal_regimes):
        """Normalised scores sum ≈ 1."""
        for m in formal_regimes.values():
            total = sum(m.regime_scores.values())
            assert total == pytest.approx(1.0, abs=1e-10)

    def test_dominant_has_max_score(self, formal_regimes):
        """Dominant regime has the highest score."""
        for m in tuple(formal_regimes.values())[:5]:
            max_regime = max(m.regime_scores, key=m.regime_scores.get)  # type: ignore
            assert m.dominant_regime == max_regime

    def test_threshold_consistency(self, formal_regimes):
        """above_threshold flags use the equipartition share on the scores."""
        for m in tuple(formal_regimes.values())[:5]:
            for regime in ("em_like", "weak_like", "strong_like", "gravity_like"):
                assert m.above_threshold[regime] == (
                    m.regime_scores[regime] > REGIME_ACTIVITY_SHARE
                )

    def test_mixing_angle_range(self, formal_regimes):
        """Mixing angle arg(Ψ) ∈ [−π, π]."""
        for m in formal_regimes.values():
            assert -math.pi - 1e-10 <= m.mixing_angle <= math.pi + 1e-10


class TestNetworkInteractionProfile:
    """Network-wide interaction regime aggregation."""

    def test_empty_graph_has_finite_zero_order_parameters(self):
        """Empty-network aggregation must not emit NaN telemetry."""
        profile = compute_network_interaction_profile(nx.Graph())
        assert profile.mean_order_parameters == {
            "em_like": 0.0,
            "weak_like": 0.0,
            "strong_like": 0.0,
            "gravity_like": 0.0,
        }
        assert all(math.isfinite(value) for value in profile.regime_fractions.values())
        assert math.isfinite(profile.mixing_entropy)

    def test_distribution_sums_to_n(self, ws_graph, interaction_profile):
        """Sum of regime counts equals number of nodes."""
        profile = interaction_profile
        assert sum(profile.regime_distribution.values()) == ws_graph.number_of_nodes()

    def test_fractions_sum_to_one(self, ws_graph, interaction_profile):
        """Regime fractions sum to 1."""
        profile = interaction_profile
        assert sum(profile.regime_fractions.values()) == pytest.approx(1.0, abs=1e-10)

    def test_entropy_bounds(self, ws_graph, interaction_profile):
        """Shannon entropy H ∈ [0, ln(4)]."""
        profile = interaction_profile
        assert profile.mixing_entropy >= 0.0
        assert profile.mixing_entropy <= math.log(4) + 1e-10

    def test_per_node_coverage(self, ws_graph, interaction_profile):
        """per_node contains every node."""
        profile = interaction_profile
        assert set(profile.per_node.keys()) == set(ws_graph.nodes())

    def test_gauge_metrics_non_negative(self, ws_graph, interaction_profile):
        """Gauge coupling, YM action, flatness are non-negative."""
        profile = interaction_profile
        assert profile.gauge_coupling_constant >= 0.0
        assert profile.yang_mills_action >= 0.0
        assert 0.0 <= profile.gauge_flatness <= 1.0

    def test_mean_order_parameters(self, ws_graph, interaction_profile):
        """Mean order parameters are within expected range."""
        profile = interaction_profile
        for key, val in profile.mean_order_parameters.items():
            assert val >= 0.0, f"Mean O_P should be non-negative: {key}={val}"

    def test_reproducibility(self):
        """Same seed → identical profile."""
        G1 = _make_tnfr_graph(20, seed=101)
        G2 = _make_tnfr_graph(20, seed=101)
        p1 = compute_network_interaction_profile(G1)
        p2 = compute_network_interaction_profile(G2)
        assert p1.mixing_entropy == pytest.approx(p2.mixing_entropy, abs=1e-14)
        assert p1.dominant_regime == p2.dominant_regime
        assert p1.regime_distribution == p2.regime_distribution