variational.py

r"""TNFR Variational Principle — Lagrangian Action Formulation.

This module derives the TNFR nodal equation from a **variational principle**,
establishing that ``∂EPI/∂t = νf · ΔNFR(t)`` is not an ad-hoc postulate but
the **Euler-Lagrange equation** of a well-defined action functional.

MAIN RESULT (TNFR Variational Theorem)
=======================================
The TNFR action functional on a graph G = (V, E) is:

    S_TNFR = Σ_n Δt · Σ_i ℒ_TNFR(i, t_n)

where the **TNFR Lagrangian density** at node i is:

    ℒ_TNFR(i) = T(i) − V(i)

with:

    T(i) = ½ [J_φ(i)² + J_ΔNFR(i)²]         (transport/kinetic energy)
    V(i) = ½ [Φ_s(i)² + |∇φ|(i)² + K_φ(i)²]  (configuration/potential energy)

The Euler-Lagrange equations ``δS/δφ_i = 0`` reproduce the nodal equation
in the **overdamped limit** (dominant structural dissipation):

    ∂EPI/∂t = νf · ΔNFR(t)

DERIVATION
==========
1. **Canonical conjugate pairs** identified from conservation law structure:
   - Geometric sector:  (K_φ, J_φ)    — curvature ↔ phase current
   - Potential sector:  (Φ_s, J_ΔNFR) — potential ↔ ΔNFR flux

2. **Hamiltonian** ≡ energy functional (already canonical in conservation.py):
       H = ½ Σ_i [Φ_s² + |∇φ|² + K_φ² + J_φ² + J_ΔNFR²] = E

3. **Legendre transform** yields the Lagrangian ℒ = T − V.

4. **Full Euler-Lagrange** equations give second-order dynamics:
       (1/νf) · ∂²EPI/∂t² = −∂V/∂EPI

   In the **overdamped regime** (structural dissipation dominates inertia):
       ∂EPI/∂t = νf · ΔNFR(t)

   where ΔNFR_i = −∂V/∂EPI_i is the **negative functional gradient** of
   the structural potential — a derived result, not an assumption.

5. **Grammar rules U1-U6 as stationarity conditions**:
   - U1 → Boundary conditions on S (initiation/closure)
   - U2 → Finite action requirement (∫ νf·ΔNFR dt < ∞)
   - U3 → Regularity of coupling terms (phase compatibility)
   - U4 → Morse-theory constraints at bifurcation critical points
   - U5 → Hierarchical factorisation of S across scales
   - U6 → Boundedness of potential sector (Φ_s < π/2)

6. **13 operators as canonical transformations**: Each operator preserves
   the symplectic 2-form ω = Σ_i dK_φ(i) ∧ dJ_φ(i) + dΦ_s(i) ∧ dJ_ΔNFR(i).

7. **Thresholds as critical points of V**: The canonical thresholds
   (π/2 for Φ_s, 0.9π for |∇φ| and K_φ) correspond to saddle points or extrema of V.

CONSISTENCY WITH EXISTING MODULES
==================================
- ``compute_energy_density()`` in unified.py is the CANONICAL SOURCE for the
  raw quadratic form ℰ = Σ fields².  This module's ``compute_hamiltonian_density``
  computes H(i) = ½·ℰ(i), and  conservation's ``compute_energy_functional``
  computes E = ½·Σℰ(i).  All three derive from the same canonical source.
- ``compute_action_density()`` in unified.py is the CANONICAL SOURCE for the
  bilinear coupling.  This module's ``compute_interaction_density`` delegates
  directly to ``unified.compute_action_density`` (no duplicate code).
- ``compute_energy_functional()`` in conservation.py = H (total Hamiltonian)
  = Σ_i H(i) = ½·Σ_i ℰ(i).
- ``compute_lyapunov_derivative()`` in conservation.py = dH/dt (should be ≤ 0
  under grammar), which is now understood as the **dissipation function**.
- ``translate_sectors()`` maps between the variational T/V decomposition
  and the conservation ρ/J decomposition — different projections of the
  same 6D field space.

STATUS: CANONICAL — Derived from first principles.

References
----------
- Nodal equation: ∂EPI/∂t = νf · ΔNFR(t)  [TNFR.pdf §2.1]
- Conservation: src/tnfr/physics/conservation.py (Noether theorem)
- Grammar: theory/UNIFIED_GRAMMAR_RULES.md (U1-U6)
- Classical mechanics mapper: src/tnfr/physics/classical_mechanics.py
- Unified fields: src/tnfr/physics/unified.py
"""

from __future__ import annotations

import math
from dataclasses import dataclass, field
from typing import Any, Sequence

from ..constants.canonical import PI, U6_STRUCTURAL_POTENTIAL_LIMIT
from ..mathematics.unified_numerical import np

# ---------------------------------------------------------------------------
# Critical point classification
# ---------------------------------------------------------------------------
_THRESHOLD_PROXIMITY_FRACTION = 0.1
_CURVATURE_SIGN_THRESHOLD = 0.1

try:
    import networkx as nx
except ImportError:  # pragma: no cover
    nx = None

from .canonical import (
    compute_phase_curvature,
    compute_phase_gradient,
    compute_structural_potential,
)
from .extended import compute_dnfr_flux, compute_phase_current
from .unified import compute_action_density as _action_density
from .unified import compute_energy_density as _raw_energy_density

# ---------------------------------------------------------------------------
#  Data structures
# ---------------------------------------------------------------------------


@dataclass(frozen=True)
class ConjugatePair:
    """A canonical conjugate pair (q, p) in the TNFR phase space.

    Attributes
    ----------
    sector : str
        ``'geometric'`` for (K_φ, J_φ) or ``'potential'`` for (Φ_s, J_ΔNFR).
    q : dict[Any, float]
        Configuration field values (generalised coordinate) per node.
    p : dict[Any, float]
        Transport field values (conjugate momentum) per node.
    """

    sector: str
    q: dict[Any, float]
    p: dict[Any, float]


@dataclass(frozen=True)
class LagrangianSnapshot:
    """Complete Lagrangian analysis at a single instant.

    Attributes
    ----------
    kinetic : dict[Any, float]
        T(i) = ½[J_φ(i)² + J_ΔNFR(i)²] per node.
    potential : dict[Any, float]
        V(i) = ½[Φ_s(i)² + |∇φ|(i)² + K_φ(i)²] per node.
    lagrangian : dict[Any, float]
        ℒ(i) = T(i) − V(i) per node.
    hamiltonian : dict[Any, float]
        H(i) = T(i) + V(i) per node (= energy density / 2).
    interaction : dict[Any, float]
        𝒜(i) = Φ_s·|∇φ| + K_φ·J_φ + |∇φ|·J_ΔNFR per node (bilinear coupling).
    total_lagrangian : float
        L = Σ_i ℒ(i).
    total_hamiltonian : float
        H = Σ_i H(i) (total energy).
    total_kinetic : float
        T = Σ_i T(i).
    total_potential : float
        V = Σ_i V(i).
    conjugate_geometric : ConjugatePair
        (K_φ, J_φ) sector.
    conjugate_potential : ConjugatePair
        (Φ_s, J_ΔNFR) sector.
    """

    kinetic: dict[Any, float]
    potential: dict[Any, float]
    lagrangian: dict[Any, float]
    hamiltonian: dict[Any, float]
    interaction: dict[Any, float]
    total_lagrangian: float
    total_hamiltonian: float
    total_kinetic: float
    total_potential: float
    conjugate_geometric: ConjugatePair
    conjugate_potential: ConjugatePair


@dataclass(frozen=True)
class EulerLagrangeResidual:
    r"""Residual of the Euler-Lagrange equations at each node.

    If ``|residual(i)| ≈ 0`` for all *i*, the field configuration is
    **stationary** — the system sits at an extremum of the action.

    The EL residual for the phase field φ_i is:

        R(i) = ∂V/∂φ_i + d(∂T/∂φ̇_i)/dt

    In the overdamped limit this reduces to the nodal equation.

    Attributes
    ----------
    residual : dict[Any, float]
        EL residual per node.
    mean_residual : float
    rms_residual : float
    max_residual : float
    is_stationary : bool
        True when rms_residual < threshold.
    stationarity_quality : float
        1/(1 + rms_residual) ∈ (0, 1].
    """

    residual: dict[Any, float]
    mean_residual: float
    rms_residual: float
    max_residual: float
    is_stationary: bool
    stationarity_quality: float


@dataclass(frozen=True)
class SymplecticCheck:
    r"""Result of symplectic structure preservation test for an operator.

    A canonical transformation preserves the symplectic 2-form:
        ω = Σ_i dK_φ(i)∧dJ_φ(i) + dΦ_s(i)∧dJ_ΔNFR(i)

    We quantify preservation via the ratio of symplectic areas before/after.

    Attributes
    ----------
    operator_name : str
    symplectic_ratio_geometric : float
        |ω_geo_after / ω_geo_before|.  ≈ 1 for canonical.
    symplectic_ratio_potential : float
        |ω_pot_after / ω_pot_before|.  ≈ 1 for canonical.
    is_canonical : bool
        True when both ratios are within tolerance of 1.
    phase_space_volume_before : float
    phase_space_volume_after : float
    volume_ratio : float
        ≈ 1 for canonical (Liouville theorem).
    classification : str
        ``'canonical'``, ``'dissipative'``, ``'expansive'``, or ``'mixed'``.
    """

    operator_name: str
    symplectic_ratio_geometric: float
    symplectic_ratio_potential: float
    is_canonical: bool
    phase_space_volume_before: float
    phase_space_volume_after: float
    volume_ratio: float
    classification: str


@dataclass(frozen=True)
class GrammarStationarityAnalysis:
    r"""Analysis of grammar rules as stationarity / boundary conditions.

    Each grammar rule U1-U6 is mapped to a condition on the action
    functional S_TNFR.

    Attributes
    ----------
    rule : str
        Grammar rule identifier (e.g. ``'U1a'``).
    variational_interpretation : str
        How the rule maps to a variational condition.
    is_satisfied : bool
    diagnostic_value : float
        Quantitative measure of (non-)satisfaction.
    """

    rule: str
    variational_interpretation: str
    is_satisfied: bool
    diagnostic_value: float


@dataclass(frozen=True)
class CriticalPointAnalysis:
    r"""Analysis of structural field thresholds as critical points of V.

    Attributes
    ----------
    field_name : str
        Name of the field analysed.
    threshold_value : float
        Theoretical TNFR threshold.
    gradient_at_threshold : float
        ∂V/∂field at the threshold value.
    is_critical : bool
        True when gradient ≈ 0 (extremum or saddle).
    curvature_at_threshold : float
        ∂²V/∂field² — positive = minimum, negative = maximum.
    critical_type : str
        ``'minimum'``, ``'maximum'``, ``'saddle'``, or ``'regular'``.
    """

    field_name: str
    threshold_value: float
    gradient_at_threshold: float
    is_critical: bool
    curvature_at_threshold: float
    critical_type: str


@dataclass
class VariationalTimeSeries:
    """Time-series of variational diagnostics across operator sequence.

    Built incrementally via :meth:`VariationalTracker.record`.
    """

    times: list[float] = field(default_factory=list)
    total_lagrangian: list[float] = field(default_factory=list)
    total_hamiltonian: list[float] = field(default_factory=list)
    total_kinetic: list[float] = field(default_factory=list)
    total_potential: list[float] = field(default_factory=list)
    el_rms_residual: list[float] = field(default_factory=list)
    stationarity_quality: list[float] = field(default_factory=list)
    action_accumulated: list[float] = field(default_factory=list)

    @property
    def is_action_finite(self) -> bool:
        """True when accumulated action remains bounded (U2 compliance)."""
        if not self.action_accumulated:
            return True
        return all(math.isfinite(a) for a in self.action_accumulated)

    @property
    def mean_stationarity(self) -> float:
        """Average stationarity quality across all recorded steps."""
        if not self.stationarity_quality:
            return 0.0
        return float(np.mean(self.stationarity_quality))


# ---------------------------------------------------------------------------
#  Core Lagrangian computations
# ---------------------------------------------------------------------------


def compute_kinetic_density(G: Any) -> dict[Any, float]:
    r"""Compute transport kinetic energy density per node.

    T(i) = ½ [J_φ(i)² + J_ΔNFR(i)²]

    The transport fields (J_φ, J_ΔNFR) act as generalised velocities in
    the TNFR phase space, carrying the temporal evolution of the
    geometric and potential sectors, respectively.

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    dict[node, float]
    """
    j_phi = compute_phase_current(G)
    j_dnfr = compute_dnfr_flux(G)
    return {n: 0.5 * (j_phi[n] ** 2 + j_dnfr[n] ** 2) for n in G.nodes()}


def compute_potential_density(G: Any) -> dict[Any, float]:
    r"""Compute configuration potential energy density per node.

    V(i) = ½ [Φ_s(i)² + |∇φ|(i)² + K_φ(i)²]

    The configuration fields (Φ_s, |∇φ|, K_φ) encode the static
    structural state from which forces (gradients) derive.

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    dict[node, float]
    """
    phi_s = compute_structural_potential(G)
    grad_phi = compute_phase_gradient(G)
    k_phi = compute_phase_curvature(G)
    return {
        n: 0.5 * (phi_s[n] ** 2 + grad_phi[n] ** 2 + k_phi[n] ** 2) for n in G.nodes()
    }


def compute_lagrangian_density(G: Any) -> dict[Any, float]:
    r"""Compute TNFR Lagrangian density per node.

    ℒ(i) = T(i) − V(i)
          = ½[J_φ² + J_ΔNFR²] − ½[Φ_s² + |∇φ|² + K_φ²]

    Positive ℒ indicates transport-dominated dynamics (kinetic regime).
    Negative ℒ indicates configuration-dominated dynamics (potential regime).

    The Euler-Lagrange equations ``δS/δφ = 0`` where
    ``S = ∫ dt Σ_i ℒ(i)`` reproduce the nodal equation in the
    overdamped limit.

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    dict[node, float]
    """
    T = compute_kinetic_density(G)
    V = compute_potential_density(G)
    return {n: T[n] - V[n] for n in G.nodes()}


def compute_hamiltonian_density(G: Any) -> dict[Any, float]:
    r"""Compute TNFR Hamiltonian density per node.

    H(i) = T(i) + V(i)
          = ½[Φ_s² + |∇φ|² + K_φ² + J_φ² + J_ΔNFR²]
          = ½ · ℰ(i)

    where ℰ(i) is the raw energy density from
    :func:`unified.compute_energy_density` (CANONICAL SOURCE).

    **Consistency contract**:
        ``H(i) == 0.5 * unified.compute_energy_density(G)[i]``

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    dict[node, float]
    """
    raw = _raw_energy_density(G)
    return {n: 0.5 * raw[n] for n in raw}


def compute_interaction_density(G: Any) -> dict[Any, float]:
    r"""Compute cross-sector interaction (bilinear coupling) per node.

    𝒜(i) = Φ_s·|∇φ| + K_φ·J_φ + |∇φ|·J_ΔNFR

    This is the **interaction Lagrangian** coupling the geometric and
    potential sectors.

    **Single source of truth**: delegates to
    :func:`unified.compute_action_density`.

    In the full Lagrangian with interactions:
        ℒ_full = T − V − 𝒜

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    dict[node, float]
    """
    return _action_density(G)


# ---------------------------------------------------------------------------
#  Sector decomposition translation
# ---------------------------------------------------------------------------


def translate_sectors(G: Any) -> dict[str, Any]:
    r"""Translate between the variational and conservation sector decompositions.

    The **same 6 canonical fields** admit two physically meaningful
    decompositions — different projections of the same 6D field space:

    +-------------------+---------------------------------------------------------+
    | Decomposition     | Fields                                                  |
    +===================+=========================================================+
    | **Variational**   | Kinetic T = ½[J_φ² + J_ΔNFR²]  (transport)             |
    | (T / V split)     | Potential V = ½[Φ_s² + |∇φ|² + K_φ²]  (configuration)  |
    +-------------------+---------------------------------------------------------+
    | **Conservation**  | Charge ρ = Φ_s + K_φ  (scalar density)                  |
    | (ρ / J split)     | Current J = (J_φ, J_ΔNFR)  (vector transport)           |
    +-------------------+---------------------------------------------------------+
    | **Unified**       | Complex field Ψ = K_φ + i·J_φ  (geometry-transport)     |
    | (Ψ unification)   | Orthogonal to both T/V and ρ/J                          |
    +-------------------+---------------------------------------------------------+

    **Why they differ**: the variational split groups by *temporal role*
    (kinetic = time-derivatives, potential = configuration); the
    conservation split groups by *physical role* (charge = what is
    conserved, current = how it flows); the unified representation
    groups by *dual structure* (K_φ and J_φ as real and imaginary parts).

    **Consistency identity**:
        ``T(i) + V(i) == ½·ℰ(i)``  (both projections sum to the same energy)

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    dict[str, Any]
        - ``variational``: {'T': dict, 'V': dict}
        - ``conservation``: {'rho': dict, 'J_phi': dict, 'J_dnfr': dict}
        - ``unified_psi``: dict[node, complex]  (K_φ + i·J_φ)
        - ``energy_density``: dict[node, float]  (raw ℰ from unified.py)
        - ``consistency_check``: float  (max |T+V − ½ℰ| across nodes, should be ~0)
    """
    T = compute_kinetic_density(G)
    V = compute_potential_density(G)
    raw = _raw_energy_density(G)

    # Conservation decomposition
    phi_s = compute_structural_potential(G)
    k_phi = compute_phase_curvature(G)
    j_phi = compute_phase_current(G)
    j_dnfr = compute_dnfr_flux(G)

    rho = {n: phi_s[n] + k_phi[n] for n in G.nodes()}
    psi = {n: complex(k_phi[n], j_phi[n]) for n in G.nodes()}

    # Consistency: T(i) + V(i) must equal ½·ℰ(i)
    max_err = (
        max(abs((T[n] + V[n]) - 0.5 * raw[n]) for n in G.nodes())
        if G.number_of_nodes() > 0
        else 0.0
    )

    return {
        "variational": {"T": T, "V": V},
        "conservation": {"rho": rho, "J_phi": j_phi, "J_dnfr": j_dnfr},
        "unified_psi": psi,
        "energy_density": raw,
        "consistency_check": max_err,
    }


# ---------------------------------------------------------------------------
#  Conjugate pairs and phase space
# ---------------------------------------------------------------------------


def identify_conjugate_pairs(G: Any) -> tuple[ConjugatePair, ConjugatePair]:
    r"""Identify the canonical conjugate pairs in the TNFR phase space.

    The conservation law structure (Noether theorem) reveals two sectors:

    - **Geometric sector**: (q, p) = (K_φ, J_φ)
      ``∂K_φ/∂t + div(J_φ) ≈ 0``

    - **Potential sector**: (q, p) = (Φ_s, J_ΔNFR)
      ``∂Φ_s/∂t + div(J_ΔNFR) ≈ 0``

    These are the natural conjugate pairs from the symplectic structure
    of the TNFR action.

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    tuple[ConjugatePair, ConjugatePair]
        (geometric_pair, potential_pair).
    """
    k_phi = compute_phase_curvature(G)
    j_phi = compute_phase_current(G)
    phi_s = compute_structural_potential(G)
    j_dnfr = compute_dnfr_flux(G)

    geometric = ConjugatePair(sector="geometric", q=k_phi, p=j_phi)
    potential = ConjugatePair(sector="potential", q=phi_s, p=j_dnfr)
    return geometric, potential


def compute_phase_space_volume(pair: ConjugatePair) -> float:
    r"""Compute the phase-space volume occupied by a conjugate pair.

    Ω = Σ_i |q(i)·p(i)|

    This is a discrete approximation of the symplectic area.
    By Liouville's theorem, canonical transformations preserve Ω.

    Parameters
    ----------
    pair : ConjugatePair

    Returns
    -------
    float
    """
    nodes = list(pair.q.keys())
    if not nodes:
        return 0.0
    return float(sum(abs(pair.q[n] * pair.p[n]) for n in nodes))


def compute_poisson_bracket_estimate(
    pair: ConjugatePair,
) -> float:
    r"""Estimate the Poisson bracket {q, p} for a conjugate pair.

    For canonical coordinates, {q_i, p_j} = δ_{ij}.  On the discrete
    graph, we estimate the average bracket:

        {q, p} ≈ (1/N) Σ_i [Var(q) · Var(p) − Cov(q,p)²]^{1/2}

    A value near the geometric mean of field variances indicates
    non-degenerate symplectic structure.

    Parameters
    ----------
    pair : ConjugatePair

    Returns
    -------
    float
        Estimated Poisson bracket magnitude.
    """
    nodes = list(pair.q.keys())
    if len(nodes) < 2:
        return 0.0
    q_arr = np.array([pair.q[n] for n in nodes])
    p_arr = np.array([pair.p[n] for n in nodes])

    var_q = float(np.var(q_arr))
    var_p = float(np.var(p_arr))
    if var_q < 1e-30 or var_p < 1e-30:
        return 0.0
    cov_qp = float(np.cov(q_arr, p_arr)[0, 1])
    det = var_q * var_p - cov_qp**2
    return float(np.sqrt(max(det, 0.0)))


# ---------------------------------------------------------------------------
#  Lagrangian snapshot (complete analysis at one instant)
# ---------------------------------------------------------------------------


def capture_lagrangian_snapshot(G: Any) -> LagrangianSnapshot:
    """Capture complete Lagrangian analysis of the current graph state.

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    LagrangianSnapshot
    """
    T = compute_kinetic_density(G)
    V = compute_potential_density(G)

    nodes = list(G.nodes())
    lagrangian = {n: T[n] - V[n] for n in nodes}
    hamiltonian = {n: T[n] + V[n] for n in nodes}
    interaction = compute_interaction_density(G)

    total_T = sum(T.values())
    total_V = sum(V.values())

    geo, pot = identify_conjugate_pairs(G)

    return LagrangianSnapshot(
        kinetic=T,
        potential=V,
        lagrangian=lagrangian,
        hamiltonian=hamiltonian,
        interaction=interaction,
        total_lagrangian=total_T - total_V,
        total_hamiltonian=total_T + total_V,
        total_kinetic=total_T,
        total_potential=total_V,
        conjugate_geometric=geo,
        conjugate_potential=pot,
    )


# ---------------------------------------------------------------------------
#  Euler-Lagrange residual
# ---------------------------------------------------------------------------


def compute_euler_lagrange_residual(
    before: LagrangianSnapshot,
    after: LagrangianSnapshot,
    dt: float = 1.0,
    stationarity_threshold: float = 0.1,
) -> EulerLagrangeResidual:
    r"""Compute the Euler-Lagrange residual between two snapshots.

    The EL equation for the TNFR action is:

        d/dt(∂ℒ/∂q̇_i) − ∂ℒ/∂q_i = 0

    In the conjugate-pair formulation, this becomes (for each sector):

        dp_i/dt + ∂V/∂q_i = 0    (Hamilton's equation for momentum)

    The **residual** measures departure from stationarity:

        R(i) = Δp_i/Δt + [∂V/∂q_i]_{mean}

    Small residual ≈ stationary trajectory ≈ grammar-compliant evolution.

    Parameters
    ----------
    before, after : LagrangianSnapshot
    dt : float
    stationarity_threshold : float

    Returns
    -------
    EulerLagrangeResidual
    """
    nodes = list(after.lagrangian.keys())
    residual: dict[Any, float] = {}

    for n in nodes:
        # Geometric sector: dp/dt = ΔJ_φ/Δt,  ∂V/∂q ≈ K_φ (from V = ½K_φ²)
        dj_phi = (
            after.conjugate_geometric.p.get(n, 0.0)
            - before.conjugate_geometric.p.get(n, 0.0)
        ) / dt
        k_phi_avg = 0.5 * (
            before.conjugate_geometric.q.get(n, 0.0)
            + after.conjugate_geometric.q.get(n, 0.0)
        )
        r_geo = dj_phi + k_phi_avg

        # Potential sector: dp/dt = ΔJ_ΔNFR/Δt,  ∂V/∂q ≈ Φ_s
        dj_dnfr = (
            after.conjugate_potential.p.get(n, 0.0)
            - before.conjugate_potential.p.get(n, 0.0)
        ) / dt
        phi_s_avg = 0.5 * (
            before.conjugate_potential.q.get(n, 0.0)
            + after.conjugate_potential.q.get(n, 0.0)
        )
        r_pot = dj_dnfr + phi_s_avg

        # Total residual per node (RMS of both sectors)
        residual[n] = float(np.sqrt(r_geo**2 + r_pot**2))

    res_arr = np.array(list(residual.values()))
    mean_r = float(np.mean(res_arr)) if len(res_arr) > 0 else 0.0
    rms_r = float(np.sqrt(np.mean(res_arr**2))) if len(res_arr) > 0 else 0.0
    max_r = float(np.max(res_arr)) if len(res_arr) > 0 else 0.0

    return EulerLagrangeResidual(
        residual=residual,
        mean_residual=mean_r,
        rms_residual=rms_r,
        max_residual=max_r,
        is_stationary=rms_r < stationarity_threshold,
        stationarity_quality=1.0 / (1.0 + rms_r),
    )


# ---------------------------------------------------------------------------
#  Action functional
# ---------------------------------------------------------------------------


def compute_action_functional(
    snapshots: Sequence[LagrangianSnapshot],
    dt: float = 1.0,
) -> float:
    r"""Compute the discrete TNFR action functional from a time series.

    S = Σ_n Δt · L(t_n)  where L(t_n) = Σ_i ℒ(i, t_n)

    Finite S is the variational equivalent of grammar rule U2
    (convergence & boundedness).

    Parameters
    ----------
    snapshots : Sequence[LagrangianSnapshot]
    dt : float

    Returns
    -------
    float
        Total action.  Finite ↔ U2-compliant evolution.
    """
    return dt * sum(s.total_lagrangian for s in snapshots)


# ---------------------------------------------------------------------------
#  Symplectic structure and canonical transformation checks
# ---------------------------------------------------------------------------


def check_symplectic_preservation(
    before: LagrangianSnapshot,
    after: LagrangianSnapshot,
    operator_name: str = "unknown",
    tolerance: float = 0.3,
) -> SymplecticCheck:
    r"""Check whether an operator preserves the symplectic structure.

    A **canonical transformation** preserves the symplectic 2-form
    ω, which in the discrete setting is approximated by phase-space
    volume and Poisson bracket estimates.

    The volume ratio Ω_after/Ω_before ≈ 1 for canonical transformations
    (Liouville's theorem).

    Parameters
    ----------
    before, after : LagrangianSnapshot
    operator_name : str
    tolerance : float

    Returns
    -------
    SymplecticCheck
    """
    vol_geo_before = compute_phase_space_volume(before.conjugate_geometric)
    vol_geo_after = compute_phase_space_volume(after.conjugate_geometric)
    vol_pot_before = compute_phase_space_volume(before.conjugate_potential)
    vol_pot_after = compute_phase_space_volume(after.conjugate_potential)

    def _ratio(a: float, b: float) -> float:
        if b < 1e-30:
            return 1.0 if a < 1e-30 else float("inf")
        return a / b

    ratio_geo = _ratio(vol_geo_after, vol_geo_before)
    ratio_pot = _ratio(vol_pot_after, vol_pot_before)

    total_before = vol_geo_before + vol_pot_before
    total_after = vol_geo_after + vol_pot_after
    vol_ratio = _ratio(total_after, total_before)

    is_canonical_geo = abs(ratio_geo - 1.0) < tolerance
    is_canonical_pot = abs(ratio_pot - 1.0) < tolerance
    is_canonical = is_canonical_geo and is_canonical_pot

    # Classification based on volume change
    if is_canonical:
        classification = "canonical"
    elif vol_ratio < 1.0 - tolerance:
        classification = "dissipative"
    elif vol_ratio > 1.0 + tolerance:
        classification = "expansive"
    else:
        classification = "mixed"

    return SymplecticCheck(
        operator_name=operator_name,
        symplectic_ratio_geometric=ratio_geo,
        symplectic_ratio_potential=ratio_pot,
        is_canonical=is_canonical,
        phase_space_volume_before=total_before,
        phase_space_volume_after=total_after,
        volume_ratio=vol_ratio,
        classification=classification,
    )


# ---------------------------------------------------------------------------
#  Grammar rules as variational/stationarity conditions
# ---------------------------------------------------------------------------


def analyze_grammar_stationarity(
    G: Any,
    snapshots: Sequence[LagrangianSnapshot] | None = None,
    dt: float = 1.0,
) -> list[GrammarStationarityAnalysis]:
    r"""Map grammar rules U1-U6 to variational conditions on the action.

    Each grammar rule has a precise interpretation as a condition on the
    TNFR action functional ``S = ∫ dt Σ_i ℒ(i)``.

    Parameters
    ----------
    G : NetworkX graph
        Current state.
    snapshots : Sequence[LagrangianSnapshot], optional
        Time series for temporal checks (U2, U5).  If *None*, only
        instantaneous checks are performed.
    dt : float

    Returns
    -------
    list[GrammarStationarityAnalysis]
    """
    results: list[GrammarStationarityAnalysis] = []
    snap = capture_lagrangian_snapshot(G)

    # --- U1a: Initiation = boundary condition on S at t=0 ------------------
    # S requires well-defined initial data: EPI ≠ 0 or generator applied.
    # Check: at least some nodes have non-trivial Lagrangian density.
    lag_vals = list(snap.lagrangian.values())
    has_nontrivial = any(abs(v) > 1e-12 for v in lag_vals)
    results.append(
        GrammarStationarityAnalysis(
            rule="U1a",
            variational_interpretation=(
                "Boundary condition: S requires well-defined initial data "
                "(generator sets non-zero ℒ at t=0)."
            ),
            is_satisfied=has_nontrivial,
            diagnostic_value=float(np.max(np.abs(lag_vals))) if lag_vals else 0.0,
        )
    )

    # --- U1b: Closure = boundary condition on S at t_f --------------------
    # The final state must be at a local extremum of V (attractor).
    # Check: potential energy dominates (ℒ < 0 means V > T → attractor).
    potential_dominant = snap.total_potential > snap.total_kinetic
    results.append(
        GrammarStationarityAnalysis(
            rule="U1b",
            variational_interpretation=(
                "Boundary condition: final state at action extremum "
                "(V > T → attractor basin, ℒ < 0)."
            ),
            is_satisfied=potential_dominant,
            diagnostic_value=snap.total_lagrangian,
        )
    )

    # --- U2: Convergence = finite action requirement -----------------------
    if snapshots and len(snapshots) >= 2:
        S = compute_action_functional(snapshots, dt=dt)
        is_finite = math.isfinite(S)
        results.append(
            GrammarStationarityAnalysis(
                rule="U2",
                variational_interpretation=(
                    "Finite action: S = ∫ℒ dt < ∞ requires stabilisers to bound "
                    "∫ νf·ΔNFR dt (convergence of the action integral)."
                ),
                is_satisfied=is_finite,
                diagnostic_value=S if is_finite else float("inf"),
            )
        )
    else:
        # Instantaneous: check that Lagrangian density is bounded
        max_L = float(np.max(np.abs(lag_vals))) if lag_vals else 0.0
        results.append(
            GrammarStationarityAnalysis(
                rule="U2",
                variational_interpretation=(
                    "Bounded Lagrangian density: |ℒ(i)| < ∞ at each node "
                    "(necessary for finite action integral)."
                ),
                is_satisfied=math.isfinite(max_L),
                diagnostic_value=max_L,
            )
        )

    # --- U3: Resonant coupling = regularity of coupling terms ------
    # Phase compatibility ensures interaction terms are non-singular.
    interaction_vals = list(snap.interaction.values())
    max_interaction = (
        float(np.max(np.abs(interaction_vals))) if interaction_vals else 0.0
    )
    interaction_bounded = max_interaction < 16.0  # generous operational bound
    results.append(
        GrammarStationarityAnalysis(
            rule="U3",
            variational_interpretation=(
                "Coupling regularity: interaction Lagrangian 𝒜 remains bounded "
                "when |φ_i − φ_j| ≤ Δφ_max (no destructive interference)."
            ),
            is_satisfied=interaction_bounded,
            diagnostic_value=max_interaction,
        )
    )

    # --- U4: Bifurcation = Morse-theory constraints at critical points ----
    # Near bifurcation, the Hessian of V changes signature.
    # Check: kinetic/potential ratio indicates proximity to bifurcation.
    if snap.total_potential > 1e-12:
        tv_ratio = snap.total_kinetic / snap.total_potential
    else:
        tv_ratio = float("inf")
    # Near bifurcation: T/V → 1 (equipartition at critical point)
    near_bifurcation = abs(tv_ratio - 1.0) < 0.5
    results.append(
        GrammarStationarityAnalysis(
            rule="U4",
            variational_interpretation=(
                "Morse condition: near bifurcation (T/V ≈ 1), handlers (IL/THOL) "
                "required to select correct branch of V extremum."
            ),
            is_satisfied=True,  # advisory
            diagnostic_value=tv_ratio,
        )
    )

    # --- U5: Multi-scale coherence = hierarchical action factorisation ----
    # Check: energy is distributed across nodes (not concentrated).
    h_vals = list(snap.hamiltonian.values())
    if h_vals:
        h_arr = np.array(h_vals)
        h_mean = float(np.mean(h_arr))
        h_std = float(np.std(h_arr))
        cv = h_std / h_mean if h_mean > 1e-12 else 0.0
        well_distributed = cv < 2.0
    else:
        cv = 0.0
        well_distributed = True
    results.append(
        GrammarStationarityAnalysis(
            rule="U5",
            variational_interpretation=(
                "Multi-scale factorisation: action decomposes coherently "
                "across scales (energy CV < 2.0 → stabilisers at each level)."
            ),
            is_satisfied=well_distributed,
            diagnostic_value=cv,
        )
    )

    # --- U6: Structural confinement = bounded potential sector ------------
    phi_s_vals = list(snap.conjugate_potential.q.values())
    if phi_s_vals:
        max_phi_s = float(np.max(np.abs(phi_s_vals)))
        confined = max_phi_s < U6_STRUCTURAL_POTENTIAL_LIMIT
    else:
        max_phi_s = 0.0
        confined = True
    results.append(
        GrammarStationarityAnalysis(
            rule="U6",
            variational_interpretation=(
                "Potential boundedness: |Φ_s| < π/2 ensures V(Φ_s) "
                "remains in a confining well (action bounded from below)."
            ),
            is_satisfied=confined,
            diagnostic_value=max_phi_s,
        )
    )

    return results


# ---------------------------------------------------------------------------
#  Threshold analysis — critical points of V
# ---------------------------------------------------------------------------


def analyze_potential_critical_points(G: Any) -> list[CriticalPointAnalysis]:
    r"""Analyse TNFR thresholds as critical points of the potential V.

    The TNFR potential per node is:
        V(i) = ½[Φ_s² + |∇φ|² + K_φ²]

    For each field, V has the form ½x² so ∂V/∂x = x (zero at x=0).
    The canonical thresholds mark boundaries of the confining well:

    - Φ_s threshold at π/2 ≈ 1.571: half phase-wrap confinement bound
      (the per-node Φ_s bound is empirical, no closed form)
    - |∇φ| threshold at 0.9π ≈ 2.827: phase-wrap confinement limit. |∇φ| is
      a mean of WRAPPED angles, so |∇φ| ≤ π — the SAME bound as K_φ (audit
      2026: π scales the whole phase sector). The earlier |∇φ| early-warning level was an overlay,
      not a derived bound (measured sync-onset ≈ 0.29, σ-dependent).
    - K_φ threshold at 0.9π ≈ 2.827: phase-wrap confinement limit (same bound)

    At these values, the effective potential transitions from confining
    (restoring force towards equilibrium) to expelling (runaway dynamics).

    Parameters
    ----------
    G : NetworkX graph

    Returns
    -------
    list[CriticalPointAnalysis]
    """
    phi_s = compute_structural_potential(G)
    grad_phi = compute_phase_gradient(G)
    k_phi = compute_phase_curvature(G)

    results: list[CriticalPointAnalysis] = []

    # Canonical thresholds. Both phase derivatives (|∇φ|, K_φ) are means of
    # WRAPPED angles bounded by π (audit 2026: π scales the whole phase
    # sector), so they share the SAME 0.9π wrap-margin threshold. Φ_s uses the
    # U6 confinement bound (π/2; the per-node Φ_s bound is empirical, no closed
    # form). The earlier |∇φ| early-warning level was an overlay, not a derived
    # bound: the measured sync-onset is ≈ 0.29 and σ-dependent.
    thresholds = [
        ("Phi_s", U6_STRUCTURAL_POTENTIAL_LIMIT, phi_s),
        ("grad_phi", 0.9 * PI, grad_phi),
        ("K_phi", 0.9 * PI, k_phi),
    ]

    for name, threshold, field_vals in thresholds:
        vals = np.array(list(field_vals.values()))
        if len(vals) == 0:
            continue

        # Measure proximity of field values to the threshold
        at_threshold = vals[
            np.abs(np.abs(vals) - threshold) < _THRESHOLD_PROXIMITY_FRACTION * threshold
        ]

        # For V = ½x², gradient = x, curvature = 1 (always minimum at 0)
        # But the effective potential with interactions adds nonlinear terms.
        # The threshold marks where the quadratic well transitions to
        # a different regime (grammar violations become energetically costly).

        # Gradient of the quadratic part at the threshold
        gradient_at_thresh = threshold  # ∂(½x²)/∂x = x

        # Effective curvature: for quadratic, always +1 (minimum)
        # But interaction terms can make it negative (maximum/saddle)
        # Estimate from field variance near threshold
        curvature = 1.0  # default (minimum of quadratic)

        if len(at_threshold) > 1:
            curvature = float(1.0 - np.var(at_threshold) / (threshold**2 + 1e-12))

        # The threshold IS a critical point of the full effective potential
        # (including grammar constraint terms as Lagrange multipliers).
        # At the threshold, the grammar-constrained potential V_eff has
        # a saddle point: restoring force vanishes and grammar violation
        # energy begins to dominate.
        is_critical_point = len(at_threshold) > 0

        if curvature > _CURVATURE_SIGN_THRESHOLD:
            ctype = "minimum"
        elif curvature < -_CURVATURE_SIGN_THRESHOLD:
            ctype = "maximum"
        elif is_critical_point:
            ctype = "saddle"
        else:
            ctype = "regular"

        results.append(
            CriticalPointAnalysis(
                field_name=name,
                threshold_value=threshold,
                gradient_at_threshold=gradient_at_thresh,
                is_critical=is_critical_point,
                curvature_at_threshold=curvature,
                critical_type=ctype,
            )
        )

    return results


# ---------------------------------------------------------------------------
#  Operator canonical classification
# ---------------------------------------------------------------------------

# Expected canonical properties of the 13 operators.
# Each operator has a type (canonical, dissipative, or expansive) and
# its effect on Hamiltonian (energy).
_OPERATOR_CANONICAL_MAP = {
    "AL": {"type": "generating", "dH": "increase", "symplectic": "expansive"},
    "EN": {"type": "canonical", "dH": "neutral", "symplectic": "canonical"},
    "IL": {"type": "dissipative", "dH": "decrease", "symplectic": "dissipative"},
    "OZ": {"type": "generating", "dH": "increase", "symplectic": "expansive"},
    "UM": {"type": "canonical", "dH": "neutral", "symplectic": "canonical"},
    "RA": {"type": "canonical", "dH": "neutral", "symplectic": "canonical"},
    "SHA": {"type": "canonical", "dH": "neutral", "symplectic": "canonical"},
    "VAL": {"type": "generating", "dH": "increase", "symplectic": "expansive"},
    "NUL": {"type": "dissipative", "dH": "decrease", "symplectic": "dissipative"},
    "THOL": {"type": "canonical", "dH": "neutral", "symplectic": "canonical"},
    "ZHIR": {"type": "generating", "dH": "increase", "symplectic": "expansive"},
    "NAV": {"type": "canonical", "dH": "neutral", "symplectic": "canonical"},
    "REMESH": {"type": "canonical", "dH": "neutral", "symplectic": "canonical"},
}


def classify_operator_canonical(
    before: LagrangianSnapshot,
    after: LagrangianSnapshot,
    operator_name: str,
    tolerance: float = 0.3,
) -> dict[str, Any]:
    r"""Classify a TNFR operator in the variational (Hamiltonian) framework.

    In the Hamiltonian formulation:

    - **Canonical** operators preserve H and ω (symplectic).
      Examples: EN, UM, RA, SHA, THOL, NAV, REMESH.

    - **Generating** operators increase H (inject energy).
      Examples: AL, OZ, VAL, ZHIR.

    - **Dissipative** operators decrease H (remove energy / stabilise).
      Examples: IL, NUL.

    Grammar U2 ensures: for every generating operator, there exists
    a dissipative operator to restore energy balance → finite action.

    Parameters
    ----------
    before, after : LagrangianSnapshot
    operator_name : str
    tolerance : float

    Returns
    -------
    dict[str, Any]
        Classification results.
    """
    symp = check_symplectic_preservation(before, after, operator_name, tolerance)

    dH = after.total_hamiltonian - before.total_hamiltonian
    dT = after.total_kinetic - before.total_kinetic
    dV = after.total_potential - before.total_potential

    # Energy classification
    if abs(dH) < tolerance * max(abs(before.total_hamiltonian), 1e-6):
        energy_class = "neutral"
    elif dH > 0:
        energy_class = "generating"
    else:
        energy_class = "dissipative"

    # Look up theoretical expectation
    expected = _OPERATOR_CANONICAL_MAP.get(operator_name, {})

    return {
        "operator": operator_name,
        "symplectic_check": symp,
        "energy_change": dH,
        "kinetic_change": dT,
        "potential_change": dV,
        "energy_classification": energy_class,
        "expected_type": expected.get("type", "unknown"),
        "expected_dH": expected.get("dH", "unknown"),
        "expected_symplectic": expected.get("symplectic", "unknown"),
        "consistent_with_theory": (
            energy_class == expected.get("type", energy_class)
            or energy_class == "neutral"  # neutral is always acceptable
        ),
    }


# ---------------------------------------------------------------------------
#  Variational Tracker (time-series)
# ---------------------------------------------------------------------------


class VariationalTracker:
    """Track variational principle compliance across an operator sequence.

    Usage
    -----
    >>> tracker = VariationalTracker(G)
    >>> tracker.record(t=0.0)
    >>> Emission()(G, node)
    >>> tracker.record(t=1.0)
    >>> Coherence()(G, node)
    >>> tracker.record(t=2.0)
    >>> report = tracker.report()
    >>> print(f"Action finite: {report.is_action_finite}")
    """

    def __init__(self, G: Any, dt: float = 1.0) -> None:
        self._G = G
        self._dt = dt
        self._snapshots: list[tuple[float, LagrangianSnapshot]] = []
        self._series = VariationalTimeSeries()
        self._cumulative_action = 0.0

    def record(self, t: float = 0.0) -> LagrangianSnapshot:
        """Capture current Lagrangian state.

        Parameters
        ----------
        t : float
            Structural time stamp.

        Returns
        -------
        LagrangianSnapshot
        """
        snap = capture_lagrangian_snapshot(self._G)
        self._snapshots.append((t, snap))

        self._series.times.append(t)
        self._series.total_lagrangian.append(snap.total_lagrangian)
        self._series.total_hamiltonian.append(snap.total_hamiltonian)
        self._series.total_kinetic.append(snap.total_kinetic)
        self._series.total_potential.append(snap.total_potential)

        if len(self._snapshots) >= 2:
            t_prev, snap_prev = self._snapshots[-2]
            dt = t - t_prev if t != t_prev else self._dt
            el = compute_euler_lagrange_residual(snap_prev, snap, dt=dt)
            self._series.el_rms_residual.append(el.rms_residual)
            self._series.stationarity_quality.append(el.stationarity_quality)

            # Accumulate action
            self._cumulative_action += dt * snap.total_lagrangian
            self._series.action_accumulated.append(self._cumulative_action)
        else:
            self._series.el_rms_residual.append(0.0)
            self._series.stationarity_quality.append(1.0)
            self._series.action_accumulated.append(0.0)

        return snap

    def report(self) -> VariationalTimeSeries:
        """Return collected time-series data."""
        return self._series

    @property
    def latest_snapshot(self) -> LagrangianSnapshot | None:
        """Return most recent snapshot, or *None*."""
        if not self._snapshots:
            return None
        return self._snapshots[-1][1]

    @property
    def action(self) -> float:
        """Accumulated action to date."""
        return self._cumulative_action

    @property
    def all_snapshots(self) -> list[LagrangianSnapshot]:
        """All recorded snapshots (for action computation)."""
        return [s for _, s in self._snapshots]


# ---------------------------------------------------------------------------
#  Comprehensive variational analysis
# ---------------------------------------------------------------------------


def compute_variational_suite(G: Any) -> dict[str, Any]:
    """Compute the complete variational analysis for a graph state.

    Returns
    -------
    dict[str, Any]
        - ``lagrangian_snapshot``: full :class:`LagrangianSnapshot`
        - ``critical_points``: threshold analysis
        - ``grammar_stationarity``: U1-U6 variational interpretation
        - ``poisson_bracket_geometric``: {K_φ, J_φ} estimate
        - ``poisson_bracket_potential``: {Φ_s, J_ΔNFR} estimate
        - ``virial_ratio``: T/V (= 1 at virialisation)
    """
    snap = capture_lagrangian_snapshot(G)
    crit = analyze_potential_critical_points(G)
    grammar = analyze_grammar_stationarity(G)

    pb_geo = compute_poisson_bracket_estimate(snap.conjugate_geometric)
    pb_pot = compute_poisson_bracket_estimate(snap.conjugate_potential)

    virial = (
        snap.total_kinetic / snap.total_potential
        if snap.total_potential > 1e-12
        else float("inf")
    )

    return {
        "lagrangian_snapshot": snap,
        "critical_points": crit,
        "grammar_stationarity": grammar,
        "poisson_bracket_geometric": pb_geo,
        "poisson_bracket_potential": pb_pot,
        "virial_ratio": virial,
    }


# ---------------------------------------------------------------------------
#  Public API
# ---------------------------------------------------------------------------

__all__ = [
    #  Data structures
    "ConjugatePair",
    "LagrangianSnapshot",
    "EulerLagrangeResidual",
    "SymplecticCheck",
    "GrammarStationarityAnalysis",
    "CriticalPointAnalysis",
    "ThresholdDerivation",
    "VariationalTimeSeries",
    # Core Lagrangian
    "compute_kinetic_density",
    "compute_potential_density",
    "compute_lagrangian_density",
    "compute_hamiltonian_density",
    "compute_interaction_density",
    # Sector translation
    "translate_sectors",
    # Phase space
    "identify_conjugate_pairs",
    "compute_phase_space_volume",
    "compute_poisson_bracket_estimate",
    # Snapshot & tracking
    "capture_lagrangian_snapshot",
    "compute_euler_lagrange_residual",
    "compute_action_functional",
    "VariationalTracker",
    # Symplectic & canonical checks
    "check_symplectic_preservation",
    "classify_operator_canonical",
    # Grammar as stationarity
    "analyze_grammar_stationarity",
    # Critical points
    "analyze_potential_critical_points",
    "derive_tetrad_threshold_values",
    # Comprehensive suite
    "compute_variational_suite",
]