structural_diffusion.py

r"""TNFR Structural Diffusion — the transport content of the nodal equation.

This module makes explicit, and verifies, that the TNFR nodal equation

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

is **structurally a diffusion equation on the network**.  This is not an
analogy imported from another paradigm: it is the literal content of the
canonical ΔNFR computation.

THE NODAL EQUATION IS GRAPH DIFFUSION
=====================================
The canonical ΔNFR (:func:`tnfr.dynamics.default_compute_delta_nfr`) is a
weighted sum of *neighbour-mean-minus-self* gradients, one per structural
channel (see :mod:`tnfr.dynamics.dnfr`):

    g_epi(i)   = mean_{j∈N(i)} EPI(j) − EPI(i)
    g_phase(i) = −angle_diff(θ(i), mean θ neighbours) / π
    g_vf(i)    = mean νf(neighbours) − νf(i)
    g_topo(i)  = mean deg(neighbours) − deg(i)

Each ``neighbour-mean − self`` term is exactly the action of the
**random-walk graph Laplacian** L_rw = I − D⁻¹W on that field:

    g_epi = −(L_rw · EPI)   (verified to machine precision).

So the EPI channel of the nodal equation is

    ∂EPI/∂t = νf · ΔNFR_epi = −νf · L_rw · EPI,

i.e. the **discrete diffusion (heat) equation** with diffusivity νf.  The
structural form EPI spreads across the network exactly as heat or a
concentration diffuses; ΔNFR is the diffusive gradient (the structural
pressure) driving the flux, and νf is the mobility / diffusivity.

WHAT EMERGES (empirically-grounded, in TNFR's own terms)
========================================================
- **Structural diffusion** (EPI channel): the form relaxes to a uniform
  field; each Laplacian eigenmode decays as exp(−νf·λ_k·t); the slowest
  rate is set by the spectral gap λ₂ (the Fiedler value).
- **Conserved structural total**: the random-walk Laplacian conserves the
  **degree-weighted total** Σ_i deg(i)·EPI(i) (its left null vector is the
  degree vector), the analogue of the conserved amount of diffusing
  substance.
- **Equilibrium ⟺ no gradients**: ΔNFR = 0 ⟺ the field is uniform across
  neighbourhoods — the diffusive steady state.
- **Synchronization** (phase channel): the phase term aligns θ to the
  neighbour mean, driving Kuramoto-type synchronization (R → 1).

These are the registers whose existence is established by the strictest
empirical method — diffusion (Fourier 1822, Fick 1855, Einstein 1905) and
synchronization (Kuramoto; observed in fireflies, pacemaker cells, neurons,
Josephson junctions).  They are reproduced here as the **same mathematics**
(the graph Laplacian is the discrete diffusion operator), not as a
metaphor.

THE MECHANICAL REGIME IS OVERDAMPED DRIFT (not inertial)
========================================================
Because the nodal equation is **first order in time**, the mechanical
regime it produces directly is the **overdamped drift law**, not Newtonian
inertia.  Reading EPI as a position-like coordinate q and ΔNFR as the
structural pressure F, the nodal equation is

    q̇ = νf · F,

i.e. **velocity proportional to applied force**, with νf the **mobility**.
Under a sustained structural pressure the field drifts at *constant*
velocity (linear in time), it does **not** accelerate.  This is the
empirically-demonstrated mobility / drift law — Stokes drag (1851),
Einstein's mobility relation (1905), terminal velocity, sedimentation,
electrophoresis — where νf is the mobility, NOT an inverse inertial mass.

The **inertial** Newtonian regime (second order, q̈ = F/m, oscillation)
is a *different* structure: it lives in the conservative **symplectic
substrate** Hamiltonian flow (:mod:`tnfr.physics.symplectic_substrate`,
where the flow is q̈ = −q per conjugate pair).  The bare nodal equation is
the **overdamped projection** of that substrate flow.  So:

    bare nodal equation (1st order)  →  overdamped drift  v = νf·F
    symplectic substrate (2nd order) →  inertial oscillation  q̈ = −∂V/∂q

both empirically grounded, but distinct regimes — a single first-order
nodal equation cannot, by itself, be Newton's second law.

DISCRETE MODES ARE THE BOUNDED-MANIFOLD STANDING WAVES
======================================================
On a **bounded** structural manifold (a finite graph) the diffusion
operator has a **discrete** spectrum of eigenmodes — the same structure as
the discrete harmonics of a bounded vibrating medium.  The symmetric
normalized Laplacian L_sym = I − D^{-1/2} W D^{-1/2} shares the diffusion
operator L_rw's spectrum {λ_k} but has **orthonormal** eigenvectors v_k:

- **Discrete spectrum**: a finite manifold supports a finite, discrete set
  of eigenvalues {λ_k} (not a continuum) — the structural origin of
  "discrete modes".  λ_1 = 0 is the uniform mode (the conserved diffusion
  mode); λ_2 (the spectral gap) is the first non-trivial mode.
- **Standing-wave shapes**: the eigenvectors v_k are orthonormal standing
  waves.  On a path graph they are exactly the cosine standing waves of a
  vibrating string (overlap 1.0 to machine precision).
- **Nodal-domain ordering** (Courant): the number of sign changes (nodal
  domains) grows with the mode index k — the structural "mode number" k
  emerges from the bounded geometry, not from a postulate.
- **Two time-regimes, same modes**: under diffusion (first order) mode k
  relaxes as exp(−νf·λ_k·t); under the wave/substrate flow (second order)
  it oscillates at the standing-wave frequency ω_k = √λ_k.

This is the discrete-harmonic structure of a bounded elastic medium —
vibrating strings (Pythagoras), Chladni plate modes (1787), molecular
vibrational spectra — all established by the strictest empirical method.
The discreteness is a consequence of the **bounded structural geometry**,
not an imported quantum postulate.

STRUCTURAL STABILITY: THE DISPERSION RELATION
=============================================
The growth or decay of each structural eigenmode under diffusion plus a
local reaction rate r is governed by the **dispersion relation**

    σ_k = r − νf · λ_k,

the universal linear-stability law (the same tool that governs every
instability and pattern-forming system — convective instability, the
onset of pattern formation).  From the canonical Laplacian spectrum {λ_k}:

- **Pure diffusion** (r = 0): σ_k = −νf·λ_k ≤ 0, so every non-uniform mode
  *decays* — structural equilibrium is stable, the integral ∫νf·ΔNFR dt
  converges.  This is the linear-stability content of "diffusion relaxes
  to uniform".
- **Structural instability threshold** r_c = νf·λ_2 (the spectral gap, the
  Fiedler value, times the diffusivity).  For 0 < r < r_c only the uniform
  mode grows (global amplification, no spatial structure); for r > r_c the
  **Fiedler mode** (k = 1) also grows — the first *structural* pattern.
- **The first structural pattern is the Fiedler partition**: the Fiedler
  eigenvector splits the network along its **weakest structural cut** (the
  two most weakly-connected communities) — the empirically-validated
  spectral-clustering result.
- **U2 grammar, spectrally**: a destabilizing reaction raises r, a
  stabilizer lowers it; bounded evolution (U2) ⟺ keeping r below r_c.
  Above r_c the Fiedler mode grows unboundedly → fragmentation (the
  U2-violation the grammar prevents).

The reaction rate r is a generic local rate; in TNFR the operators supply
it (stabilizers lower r, destabilizers raise it).  A *two-channel*
structural diffusion with **differential diffusivity** and an
activator–inhibitor coupling supports a finite-wavelength (Turing)
instability — the empirically-demonstrated pattern-formation mechanism
(Belousov–Zhabotinsky, morphogenesis); those kinetics are a model input,
not TNFR-derived, so only the dispersion-relation mechanism is certified
here.

THE STRUCTURAL RANDOM WALK AND RESISTANCE GEOMETRY
==================================================
The diffusion operator is **literally the generator of a random walk** on
the network: L_rw = I − D^{-1}W = I − P, where P = D^{-1}W is the
random-walk transition matrix (verified exactly).  So the structural
transport is **Brownian motion on the network** — the empirically-
demonstrated random walk (Einstein 1905, Perrin 1908, the proof of atoms):

- **Stationary distribution ∝ degree**: the random walk converges to
  π_i = deg(i) / Σ deg — exactly the **degree-weighted total** the
  diffusion conserves.  The conserved quantity *is* the equilibrium
  measure.
- **Effective resistance** (Ohm's law): treating the network as a
  resistor network (the combinatorial Laplacian L = D − W is the
  conductance matrix — Kirchhoff 1847), the effective resistance
  R_eff(i,j) = L⁺_ii + L⁺_jj − 2L⁺_ij (L⁺ the pseudoinverse) is a
  **transport metric** (symmetric, non-negative, triangle inequality) —
  the structural "difficulty of transport" between two nodes.
- **Commute time = 2m·R_eff**: the expected round-trip time of the random
  walk between two nodes equals 2m times the effective resistance (m the
  number of edges) — the exact link between the diffusion random walk and
  the resistance geometry (Chandra et al. 1996), confirmed against
  Monte-Carlo walks.

These are the same mathematics as Brownian motion (random walk) and
electrical networks (Ohm/Kirchhoff resistance) — both established by the
strictest empirical method.

THE STRUCTURAL FLOW: CURRENT, KIRCHHOFF, AND CONTINUITY
======================================================
The transport carries a **structural current**: the diffusion edge current
J_ij = EPI_i − EPI_j (Fick's law — flux from high to low, antisymmetric).
Its node-level balance is **Kirchhoff's current law**, which *is* the
discrete continuity equation:

    div(J)(i) = Σ_{j∼i} J_ij = (L·EPI)(i),

so the net outflow at a node equals the combinatorial Laplacian acting on
EPI.  Hence the diffusion continuity equation ∂EPI/∂t + div(J) = 0 holds,
and for a closed network (no sources) the total flux balances, Σ_i div(J)
= 0 (L has zero column sums — the structural-conservation analogue here).
Under an injected unit current from s to t the induced potential drop is
the **effective resistance** R_eff(s,t) (Ohm's law) — tying the current to
the resistance geometry above.  These are Fick diffusion, Kirchhoff's
circuit laws, and Ohm's law — all empirically demonstrated.

This is the EPI-channel current; it complements the tetrad-field
continuity of :mod:`tnfr.physics.conservation` (which tracks the charge
ρ = Φ_s + K_φ and the current J = (J_φ, J_ΔNFR)).

HONEST SCOPE
============
- The identity ΔNFR_epi = −L_rw·EPI is EXACT (machine precision), a
  mathematical fact about the canonical ΔNFR.
- The full ΔNFR is multi-channel: EPI **diffusion** + phase
  **synchronization** + νf/topology **homogenization**.  This module
  isolates and certifies the diffusion (EPI) channel and reports the
  synchronization channel qualitatively.
- **λ₂ is topological, NOT tied to the canonical constants.** The spectral
  gap λ₂ (which governs relaxation, stability, and the instability
  threshold) is a purely spectral/topological quantity — determined by N,
  degree, and connectivity (e.g. ring λ₂ = 1 − cos(2π/N), complete-graph
  λ₂ = n/(n−1)).  The overlay threshold scales of the tetrad fields do
  **not** enter the Laplacian spectrum; any numerical proximity is
  coincidental (the 2π/N in a ring is a geometric polygon angle, not a
  structural scale).  Measured
  negative result — do not assert a λ₂ ↔ constant relation.
- This characterises the transport content of the nodal dynamics; it does
  not, by itself, resolve any open program (Riemann G4, Navier–Stokes).

References
----------
- :mod:`tnfr.dynamics.dnfr` — the canonical ΔNFR neighbour-mean gradients
- :func:`tnfr.observers.kuramoto_order` — the synchronization order R
- :mod:`tnfr.physics.conservation` — the structural continuity equation
- AGENTS.md §"Foundational Physics" — the nodal equation
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Any

from ..alias import get_attr
from ..constants.aliases import ALIAS_DNFR, ALIAS_EPI, ALIAS_VF
from ..mathematics.unified_numerical import np

__all__ = [
    "StructuralDiffusionCertificate",
    "OverdampedRegimeCertificate",
    "OverdampedProjectionCertificate",
    "UndampedLimitCertificate",
    "DiscreteModeCertificate",
    "StructuralStabilityCertificate",
    "RandomWalkCertificate",
    "StructuralFlowCertificate",
    "structural_diffusion_operator",
    "symmetric_normalized_laplacian",
    "structural_field",
    "structural_diffusivity",
    "relaxation_spectrum",
    "structural_frequency_rank",
    "degree_weighted_total",
    "structural_eigenmodes",
    "nodal_domain_count",
    "compute_emergent_pulse",
    "compute_nodal_pulse",
    "dispersion_relation",
    "instability_threshold",
    "fiedler_partition",
    "random_walk_matrix",
    "stationary_distribution",
    "effective_resistance",
    "commute_time",
    "structural_current",
    "current_divergence",
    "verify_structural_diffusion",
    "verify_overdamped_regime",
    "damped_wave_rates",
    "verify_overdamped_projection",
    "verify_undamped_limit",
    "verify_discrete_modes",
    "verify_structural_stability",
    "verify_structural_random_walk",
    "verify_structural_flow",
]


def _ordered_nodes(G: Any) -> list:
    """Stable node ordering for the matrix representation."""
    return list(G.nodes())


def structural_diffusion_operator(G: Any) -> tuple[list, Any]:
    r"""Return the random-walk graph Laplacian L_rw = I − D⁻¹W.

    This is the operator whose action on a field is exactly the canonical
    ΔNFR ``neighbour-mean − self`` gradient: g = −L_rw·field.  Built from
    the (optionally weighted) adjacency; isolated nodes (degree 0) get a
    zero row (no diffusion).

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (nodes, L_rw) : tuple[list, np.ndarray]
        The node ordering and the N×N random-walk Laplacian.
    """
    nodes = _ordered_nodes(G)
    index = {n: i for i, n in enumerate(nodes)}
    n = len(nodes)
    lap = np.zeros((n, n), dtype=float)
    for node in nodes:
        i = index[node]
        neigh = list(G.neighbors(node))
        if not neigh:
            continue
        # weighted degree (weight defaults to 1.0 when absent)
        weights = [float(G[node][m].get("weight", 1.0)) for m in neigh]
        deg = sum(weights)
        if deg <= 0.0:
            continue
        lap[i, i] = 1.0
        for m, w in zip(neigh, weights):
            lap[i, index[m]] -= w / deg
    return nodes, lap


def symmetric_normalized_laplacian(
    G: Any, nodes: list | None = None
) -> tuple[list, Any]:
    r"""Return the symmetric normalized Laplacian L_sym = I − D^{-1/2} W D^{-1/2}.

    L_sym shares the spectrum of the canonical diffusion operator
    L_rw = I − D⁻¹W (:func:`structural_diffusion_operator`) but is symmetric, so
    it has an orthonormal eigenbasis and real eigenvalues — the canonical choice
    for the relaxation spectrum (its λ₂ is the structural ``diffusion_gap``).
    Isolated nodes (degree 0) get a zero row.

    Parameters
    ----------
    G : TNFRGraph
    nodes : list, optional
        Node ordering; defaults to the stable ``list(G.nodes())`` order.

    Returns
    -------
    (nodes, L_sym) : tuple[list, np.ndarray]
        The node ordering and the N×N symmetric normalized Laplacian.
    """
    if nodes is None:
        nodes = _ordered_nodes(G)
    index = {nd: i for i, nd in enumerate(nodes)}
    n = len(nodes)
    deg = np.zeros(n, dtype=float)
    for node in nodes:
        deg[index[node]] = sum(
            float(G[node][m].get("weight", 1.0)) for m in G.neighbors(node)
        )
    # Compute D^{-1/2} only on connected nodes; isolated nodes (deg = 0) keep
    # d_inv_sqrt = 0. Masked assignment avoids evaluating 1/sqrt(0) (which the
    # np.where form does for every entry before selecting, emitting a warning).
    d_inv_sqrt = np.zeros(n, dtype=float)
    positive = deg > 0.0
    d_inv_sqrt[positive] = 1.0 / np.sqrt(deg[positive])
    lap = np.zeros((n, n), dtype=float)
    for node in nodes:
        i = index[node]
        if deg[i] <= 0.0:
            continue
        lap[i, i] = 1.0
        for m in G.neighbors(node):
            j = index[m]
            w = float(G[node][m].get("weight", 1.0))
            lap[i, j] -= w * d_inv_sqrt[i] * d_inv_sqrt[j]
    return nodes, lap


def structural_field(G: Any, nodes: list | None = None) -> Any:
    r"""Return the EPI field as a vector aligned with ``nodes``."""
    if nodes is None:
        nodes = _ordered_nodes(G)
    return np.array(
        [float(get_attr(G.nodes[n], ALIAS_EPI, 0.0)) for n in nodes],
        dtype=float,
    )


def structural_diffusivity(G: Any) -> float:
    r"""Mean structural frequency νf — the diffusion coefficient (mobility).

    In ∂EPI/∂t = −νf·L_rw·EPI, νf plays the role of the diffusivity: the
    larger the structural frequency, the faster the form spreads.
    """
    nodes = _ordered_nodes(G)
    vf = [float(get_attr(G.nodes[n], ALIAS_VF, 0.0)) for n in nodes]
    return float(np.mean(vf)) if vf else 0.0


def degree_weighted_total(G: Any) -> float:
    r"""The conserved structural total Σ_i deg(i)·EPI(i).

    The random-walk Laplacian conserves the degree-weighted total (its left
    null vector is the degree vector), the analogue of the conserved amount
    of a diffusing substance.

    This is the **EPI-channel** conserved quantity (of the diffusion
    ∂EPI/∂t = −νf·L_rw·EPI).  It is **distinct** from the tetrad Noether charge
    Q = Σ(Φ_s + K_φ)
    (:func:`tnfr.physics.conservation.compute_noether_charge`), conserved under
    grammar U1–U6: TNFR carries two distinct conservation laws, on the EPI
    field and on the tetrad fields respectively (see
    STRUCTURAL_CONSERVATION_THEOREM §8.7).
    """
    nodes = _ordered_nodes(G)
    total = 0.0
    for node in nodes:
        neigh = list(G.neighbors(node))
        deg = sum(float(G[node][m].get("weight", 1.0)) for m in neigh)
        total += deg * float(get_attr(G.nodes[node], ALIAS_EPI, 0.0))
    return float(total)


def relaxation_spectrum(G: Any) -> Any:
    r"""Diffusion relaxation rates νf·λ_k (sorted ascending).

    The eigenvalues λ_k of the random-walk Laplacian L_rw scaled by the
    diffusivity νf give the decay rates of the diffusion eigenmodes:
    mode k relaxes as exp(−νf·λ_k·t).  λ₁ = 0 (the conserved uniform mode);
    λ₂ (the spectral gap / Fiedler value) sets the slowest relaxation.

    Returns
    -------
    np.ndarray
        The rates νf·λ_k sorted ascending (real parts).
    """
    # L_sym and L_rw share the spectrum; reuse the cached symmetric
    # eigendecomposition (ascending, clipped >= 0). The rates are nu_f*lambda_k.
    eig, _ = _cached_eigh(G)
    return structural_diffusivity(G) * eig


def structural_frequency_rank(G: Any, decimals: int = 8) -> int:
    r"""Number of distinct structural frequencies (the structural rank).

    The distinct eigenvalues of the canonical random-walk Laplacian L_rw are
    the network's structural frequencies (the relaxation rates of
    ∂EPI/∂t = −νf·L_rw·EPI, up to the νf scale). This returns their count s(G)
    — the size of the distinct-frequency spectrum — complementing
    ``relaxation_spectrum`` (which returns the rates themselves).

    For a connected graph, s distinct eigenvalues bound the diameter by s−1,
    and s = 2 iff the graph is complete (regular case). On arithmetic Cayley
    networks the rank is a primality / cyclotomy diagnostic (see
    :mod:`tnfr.mathematics.number_theory`): the quadratic-residue network on an
    odd prime has rank 3, and the k-th power residue network on a prime p has
    rank ``gcd(k, p-1) + 1``.

    Note
    ----
    For large or dense graphs the distinct-eigenvalue count is sensitive to the
    ``decimals`` rounding (floating-point noise in ``eigvals`` can split truly
    equal eigenvalues). For arithmetic residue networks the exact multiplicative
    :func:`tnfr.mathematics.number_theory.quadratic_residue_annotated_rank` is
    the robust closed-form object; this scalar count agrees with it for small
    moduli.

    Parameters
    ----------
    G : TNFRGraph
    decimals : int
        Rounding applied to the real and imaginary parts before counting
        distinct values (the spectrum may be complex for directed graphs).

    Returns
    -------
    int
        The number of distinct eigenvalues of L_rw.
    """
    _, lap = structural_diffusion_operator(G)
    eig = np.linalg.eigvals(lap)
    rounded = np.round(eig.real, decimals) + 1j * np.round(eig.imag, decimals)
    return int(np.unique(rounded).size)


@dataclass(frozen=True)
class StructuralDiffusionCertificate:
    r"""Verification that the nodal equation's EPI channel is graph diffusion.

    Attributes
    ----------
    n_nodes : int
    dnfr_is_graph_laplacian : bool
        The canonical ΔNFR (EPI channel) equals −L_rw·EPI.
    max_laplacian_residual : float
        Max |ΔNFR_epi − (−L_rw·EPI)| over the nodes (≈ 0).
    diffusivity : float
        Mean νf (the diffusion coefficient / mobility).
    spectral_gap : float
        λ₂ of L_rw (the Fiedler value); sets the slowest relaxation.
    slowest_relaxation_rate : float
        νf·λ₂ — the slowest diffusion decay rate.
    degree_weighted_conserved : bool
        Σ deg·EPI is conserved under the diffusion flow.
    max_conservation_drift : float
        Max drift of the degree-weighted total over the sampled flow.
    relaxes_to_uniform : bool
        The field relaxes to a spatially uniform diffusive equilibrium.
    final_field_std : float
        Std of the field after the sampled diffusion flow (≈ 0).
    """

    n_nodes: int
    dnfr_is_graph_laplacian: bool
    max_laplacian_residual: float
    diffusivity: float
    spectral_gap: float
    slowest_relaxation_rate: float
    degree_weighted_conserved: bool
    max_conservation_drift: float
    relaxes_to_uniform: bool
    final_field_std: float

    @property
    def is_valid_diffusion(self) -> bool:
        """True when the nodal EPI channel verifies as graph diffusion."""
        return (
            self.dnfr_is_graph_laplacian
            and self.degree_weighted_conserved
            and self.relaxes_to_uniform
        )

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_valid_diffusion else "INVALID"
        return (
            f"Structural diffusion [{ok}]: "
            f"ΔNFR_epi = −L_rw·EPI={self.dnfr_is_graph_laplacian} "
            f"(res {self.max_laplacian_residual:.1e}), "
            f"diffusivity νf={self.diffusivity:.4f}, "
            f"spectral gap λ₂={self.spectral_gap:.4f}, "
            f"slowest rate νf·λ₂={self.slowest_relaxation_rate:.4f}, "
            f"deg-weighted conserved={self.degree_weighted_conserved} "
            f"(drift {self.max_conservation_drift:.1e}), "
            f"relaxes to uniform={self.relaxes_to_uniform} "
            f"(final std {self.final_field_std:.1e})"
        )


def _dnfr_epi_channel(G: Any, nodes: list) -> Any:
    r"""Canonical ΔNFR restricted to the EPI channel, on a clean replica.

    Isolates the EPI diffusion channel by computing the canonical ΔNFR with
    weights (phase=0, epi=1, vf=0, topo=0) on a minimal structural replica
    (nodes + edges + EPI/θ/νf only), so the caller's graph is never mutated
    and the non-copyable runtime caches are not duplicated.
    """
    from ..dynamics import default_compute_delta_nfr

    g2 = G.__class__()
    for node in nodes:
        data = G.nodes[node]
        g2.add_node(
            node,
            EPI=float(get_attr(data, ALIAS_EPI, 0.0)),
            theta=float(data.get("theta", 0.0)),
            nu_f=float(get_attr(data, ALIAS_VF, 0.0)),
        )
    for u, v, data in G.edges(data=True):
        g2.add_edge(u, v, weight=float(data.get("weight", 1.0)))
    g2.graph["DNFR_WEIGHTS"] = {
        "phase": 0.0,
        "epi": 1.0,
        "vf": 0.0,
        "topo": 0.0,
    }
    default_compute_delta_nfr(g2)
    return np.array(
        [float(get_attr(g2.nodes[n], ALIAS_DNFR, 0.0)) for n in nodes],
        dtype=float,
    )


def verify_structural_diffusion(
    G: Any,
    *,
    dt: float = 0.1,
    steps: int = 400,
    tolerance: float = 1e-9,
) -> StructuralDiffusionCertificate:
    r"""Verify the nodal equation's EPI channel is graph diffusion.

    Confirms (1) the canonical ΔNFR EPI channel equals −L_rw·EPI to machine
    precision, (2) the degree-weighted total is conserved under the
    diffusion flow, and (3) the field relaxes to a uniform diffusive
    equilibrium; and reports the diffusivity νf and the relaxation spectrum.

    The caller's graph is never mutated (the ΔNFR check runs on a copy).

    Parameters
    ----------
    G : TNFRGraph
    dt : float
        Forward-Euler step for the diffusion-flow checks.
    steps : int
        Number of diffusion steps for the relaxation / conservation checks.
    tolerance : float
        Maximum allowed Laplacian residual and conservation drift.

    Returns
    -------
    StructuralDiffusionCertificate
    """
    nodes, lap = structural_diffusion_operator(G)
    n = len(nodes)
    epi = structural_field(G, nodes)

    # (1) ΔNFR (epi channel) == −L_rw·EPI ?
    try:
        dnfr_epi = _dnfr_epi_channel(G, nodes)
        residual = float(np.max(np.abs(dnfr_epi - (-(lap @ epi)))))
        is_laplacian = residual < max(tolerance, 1e-12)
    except Exception:
        residual = float("nan")
        is_laplacian = False

    # diffusivity and spectrum
    diffusivity = structural_diffusivity(G)
    eig = np.linalg.eigvals(lap).real
    eig.sort()
    spectral_gap = float(eig[1]) if n > 1 else 0.0
    slowest_rate = diffusivity * spectral_gap

    # degree vector for the conserved weighted total
    deg = np.array(
        [
            sum(float(G[node][m].get("weight", 1.0)) for m in G.neighbors(node))
            for node in nodes
        ],
        dtype=float,
    )

    # (2)+(3) integrate the pure diffusion flow e ← e − dt·L_rw·e
    e = epi.copy()
    w0 = float(deg @ e)
    max_drift = 0.0
    for _ in range(steps):
        e = e - dt * (lap @ e)
        max_drift = max(max_drift, abs(float(deg @ e) - w0))
    conserved = max_drift < max(tolerance, 1e-9 * (abs(w0) + 1e-12))
    final_std = float(np.std(e))
    relaxes = final_std < max(1e-3, 1e-2 * float(np.std(epi) + 1e-12))

    return StructuralDiffusionCertificate(
        n_nodes=n,
        dnfr_is_graph_laplacian=is_laplacian,
        max_laplacian_residual=residual,
        diffusivity=diffusivity,
        spectral_gap=spectral_gap,
        slowest_relaxation_rate=slowest_rate,
        degree_weighted_conserved=conserved,
        max_conservation_drift=max_drift,
        relaxes_to_uniform=relaxes,
        final_field_std=final_std,
    )


# ---------------------------------------------------------------------------
# The overdamped drift regime: the bare nodal equation is first-order
# ---------------------------------------------------------------------------


@dataclass(frozen=True)
class OverdampedRegimeCertificate:
    r"""Verification that the bare nodal equation is the overdamped drift law.

    The nodal equation ∂EPI/∂t = νf·ΔNFR is **first order in time**, so —
    reading EPI as a position q and ΔNFR as the structural pressure F — it
    is the **mobility / drift law** q̇ = νf·F: velocity proportional to
    force, with νf the mobility.  Under a sustained pressure the field
    drifts at *constant* velocity (linear in time), it does not accelerate.

    This is the empirically-demonstrated overdamped regime (Stokes 1851,
    Einstein 1905, terminal velocity, sedimentation, electrophoresis).  The
    inertial Newtonian regime (q̈ = F/m, second order) is the separate
    :mod:`tnfr.physics.symplectic_substrate` Hamiltonian flow; the nodal
    equation is its overdamped projection.

    Attributes
    ----------
    drift_velocity : float
        v = νf·F evaluated at the reference (νf, F).
    velocity_is_constant : bool
        Under sustained pressure, dEPI/dt is constant (first-order/drift).
    max_velocity_variation : float
        Max |dEPI/dt − v| over the held-pressure integration (≈ 0).
    position_linear_in_time : bool
        EPI(t) grows linearly (slope = drift), not quadratically.
    position_slope : float
        Measured slope of EPI(t) (= the drift velocity).
    mobility_linear_in_nu_f : bool
        v ∝ νf (the mobility law): v/νf is constant across νf.
    drift_linear_in_pressure : bool
        v ∝ F: v/F is constant across F.
    is_second_order : bool
        Whether the bare equation is second order (always False — it is the
        overdamped, first-order regime).
    """

    drift_velocity: float
    velocity_is_constant: bool
    max_velocity_variation: float
    position_linear_in_time: bool
    position_slope: float
    mobility_linear_in_nu_f: bool
    drift_linear_in_pressure: bool
    is_second_order: bool

    @property
    def is_overdamped_drift(self) -> bool:
        """True when the bare nodal equation verifies as overdamped drift."""
        return (
            self.velocity_is_constant
            and self.position_linear_in_time
            and self.mobility_linear_in_nu_f
            and self.drift_linear_in_pressure
            and not self.is_second_order
        )

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_overdamped_drift else "INVALID"
        return (
            f"Overdamped drift regime [{ok}]: "
            f"q̇ = νf·F = {self.drift_velocity:.4f} "
            f"(mobility law); "
            f"velocity constant={self.velocity_is_constant} "
            f"(var {self.max_velocity_variation:.1e}), "
            f"position linear={self.position_linear_in_time} "
            f"(slope {self.position_slope:.4f}), "
            f"v∝νf={self.mobility_linear_in_nu_f}, "
            f"v∝F={self.drift_linear_in_pressure}, "
            f"second-order={self.is_second_order}"
        )


def verify_overdamped_regime(
    *,
    nu_f: float = 0.7,
    pressure: float = 1.3,
    dt: float = 0.01,
    steps: int = 300,
    tolerance: float = 1e-9,
) -> OverdampedRegimeCertificate:
    r"""Verify the bare nodal equation is the overdamped drift law q̇ = νf·F.

    Integrates the canonical nodal equation
    (:func:`tnfr.dynamics.canonical.compute_canonical_nodal_derivative`)
    under a *sustained* structural pressure and measures that the EPI
    coordinate drifts at constant velocity v = νf·F (first-order, mobility
    law), linear in νf (mobility) and in the pressure F.  Uses the canonical
    nodal-equation function — no formula is re-implemented here.

    Parameters
    ----------
    nu_f : float
        Reference structural frequency (mobility).
    pressure : float
        Sustained structural pressure ΔNFR (= F).
    dt : float
        Integration step.
    steps : int
        Number of integration steps.
    tolerance : float
        Maximum allowed velocity variation / linearity residual.

    Returns
    -------
    OverdampedRegimeCertificate
    """
    from ..dynamics.canonical import compute_canonical_nodal_derivative

    # integrate the bare nodal equation under a held pressure
    epi = 0.0
    velocities = []
    positions = []
    for _ in range(steps):
        v = compute_canonical_nodal_derivative(nu_f, pressure).derivative
        epi = epi + dt * v
        velocities.append(v)
        positions.append(epi)
    vel = np.array(velocities, dtype=float)
    pos = np.array(positions, dtype=float)

    drift = float(vel[0])
    vel_var = float(np.max(np.abs(vel - drift)))
    vel_constant = vel_var < tolerance

    # position grows linearly with slope = drift (first-order, not quadratic)
    t = np.arange(steps, dtype=float) * dt
    slope, _ = np.polyfit(t, pos, 1)
    quad = np.polyfit(t, pos, 2)[0]  # leading quadratic coefficient ≈ 0
    pos_linear = abs(float(slope) - drift) < max(tolerance, 1e-6 * abs(drift)) and abs(
        float(quad)
    ) < max(tolerance, 1e-6 * abs(drift) + 1e-9)

    # mobility law: v ∝ νf (v/νf constant across νf)
    ratios_nu = [
        compute_canonical_nodal_derivative(nf, pressure).derivative / nf
        for nf in (0.2, 0.5, 1.0, 1.5)
    ]
    mobility_linear = float(np.std(ratios_nu)) < tolerance

    # drift ∝ F (v/F constant across F)
    ratios_f = [
        compute_canonical_nodal_derivative(nu_f, f).derivative / f
        for f in (0.3, 0.8, 1.3, 2.0)
    ]
    pressure_linear = float(np.std(ratios_f)) < tolerance

    return OverdampedRegimeCertificate(
        drift_velocity=drift,
        velocity_is_constant=vel_constant,
        max_velocity_variation=vel_var,
        position_linear_in_time=pos_linear,
        position_slope=float(slope),
        mobility_linear_in_nu_f=mobility_linear,
        drift_linear_in_pressure=pressure_linear,
        is_second_order=False,
    )


# ---------------------------------------------------------------------------
# Overdamped projection: the bridge from the conservative substrate wave
# to the dissipative structural diffusion
# ---------------------------------------------------------------------------


def damped_wave_rates(G: Any, gamma: float) -> tuple[Any, Any, Any]:
    r"""Per-mode slow/fast rates of the damped graph wave q̈ + γq̇ + Lq = 0.

    The conservative symplectic substrate carries the graph **wave**
    equation q̈ = −L q (second order, mode k oscillating at √λ_k — the
    standing-wave ``discrete modes`` of :func:`verify_discrete_modes`).
    Adding a damping γ gives the damped oscillator q̈ + γq̇ + L q = 0, whose
    per-mode characteristic equation is

        s² + γ s + λ_k = 0   ⟹   s± = ½(−γ ± √(γ² − 4λ_k)).

    For γ² > 4λ_k (overdamped per mode) both roots are real: a **slow** root
    s₋ → −λ_k/γ (the diffusion rate) and a **fast** root s₊ → −γ (a
    transient that dies immediately).  This is the spectral content of the
    overdamped projection: after the fast transient, mode k relaxes at
    λ_k/γ = ν_f·λ_k with ν_f = 1/γ.

    Parameters
    ----------
    G : TNFRGraph
    gamma : float
        Damping coefficient.  Its inverse is the effective diffusivity
        (mobility) ν_f = 1/γ.

    Returns
    -------
    (lambdas, s_slow, s_fast) : tuple[np.ndarray, np.ndarray, np.ndarray]
        Sorted Laplacian eigenvalues and the (real-part) slow/fast roots.
    """
    _, lap = structural_diffusion_operator(G)
    lambdas = np.sort(np.linalg.eigvals(lap).real)
    lambdas = np.clip(lambdas, 0.0, None)
    disc = gamma * gamma - 4.0 * lambdas + 0j
    root = np.sqrt(disc)
    s_slow = ((-gamma + root) / 2.0).real
    s_fast = ((-gamma - root) / 2.0).real
    return lambdas, s_slow, s_fast


@dataclass(frozen=True)
class OverdampedProjectionCertificate:
    r"""Verification that structural diffusion is the overdamped projection
    of the conservative symplectic-substrate wave flow.

    The conservative substrate (:mod:`tnfr.physics.symplectic_substrate`)
    carries the graph wave q̈ = −L q (second order).  Damping it and taking
    the strong-damping (Smoluchowski) limit collapses it onto the
    first-order structural diffusion q̇ = −(1/γ) L q, with the
    identification **ν_f = 1/γ** (structural frequency = inverse damping =
    mobility).  Both endpoints are canonical TNFR objects; this certificate
    measures the bridge between them.

    Attributes
    ----------
    n_nodes : int
    gamma : float
        Damping coefficient used for the projection.
    nu_f_effective : float
        The effective diffusivity 1/γ recovered by the projection.
    spectral_gap : float
        λ₂ of L_rw (the Fiedler value).
    lambda_max : float
        Largest Laplacian eigenvalue (sets the bridge error scale).
    max_rate_rel_error : float
        Max over modes of |s_slow + λ_k/γ| / (λ_k/γ): how far the damped
        slow rate is from the diffusion rate.
    rate_error_times_gamma_sq : float
        ``max_rate_rel_error · γ²`` — converges to ≈ λ_max, confirming the
        bridge error scales as O(λ_max/γ²).
    slowest_slow_rate : float
        Overdamped slow rate of the Fiedler mode, −s_slow(λ₂).
    slowest_diffusion_rate : float
        The diffusion spectral gap ν_f·λ₂ = λ₂/γ.
    trajectory_max_rel_error : float
        Max relative L² error between the damped-wave trajectory and the
        diffusion trajectory exp(−L t/γ)·q₀ over an overdamped time window.
    projects_to_diffusion : bool
        Whether both the rate and trajectory errors fall within tolerance.
    """

    n_nodes: int
    gamma: float
    nu_f_effective: float
    spectral_gap: float
    lambda_max: float
    max_rate_rel_error: float
    rate_error_times_gamma_sq: float
    slowest_slow_rate: float
    slowest_diffusion_rate: float
    trajectory_max_rel_error: float
    projects_to_diffusion: bool

    @property
    def is_valid_projection(self) -> bool:
        """True when the damped substrate wave projects onto diffusion."""
        return self.projects_to_diffusion

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_valid_projection else "INVALID"
        return (
            f"Overdamped projection [{ok}]: damped substrate wave "
            f"projects onto diffusion with nu_f=1/gamma="
            f"{self.nu_f_effective:.4f}; rate error "
            f"{self.max_rate_rel_error:.2e} (x gamma^2="
            f"{self.rate_error_times_gamma_sq:.3f} ~ lambda_max="
            f"{self.lambda_max:.3f}), slow gap {self.slowest_slow_rate:.5f} "
            f"vs diffusion gap {self.slowest_diffusion_rate:.5f}, "
            f"trajectory error {self.trajectory_max_rel_error:.2e}"
        )


def verify_overdamped_projection(
    G: Any,
    *,
    gamma: float = 50.0,
    n_time_samples: int = 40,
    tolerance: float = 1e-2,
) -> OverdampedProjectionCertificate:
    r"""Verify diffusion is the overdamped projection of the substrate wave.

    Builds the canonical random-walk Laplacian L_rw, forms the damped graph
    wave q̈ + γq̇ + L q = 0, and measures two things in the strong-damping
    limit: (i) the per-mode **slow rate** converges to the diffusion rate
    λ_k/γ = ν_f·λ_k (ν_f = 1/γ), with error scaling as O(λ_max/γ²); and
    (ii) the damped-wave **trajectory** (from q₀ at rest) collapses onto the
    structural-diffusion trajectory exp(−L t/γ)·q₀.  No field formula is
    re-implemented — L_rw comes from
    :func:`structural_diffusion_operator` and the orthonormal eigenbasis
    from the symmetric normalized Laplacian.

    Parameters
    ----------
    G : TNFRGraph
    gamma : float
        Damping coefficient (effective diffusivity ν_f = 1/γ).  Must satisfy
        γ² > 4·λ_max for every mode to be overdamped.
    n_time_samples : int
        Time samples in the overdamped window for the trajectory check.
    tolerance : float
        Maximum relative error (rate and trajectory) for a valid projection.

    Returns
    -------
    OverdampedProjectionCertificate
    """
    nodes, lap = structural_diffusion_operator(G)
    n = len(nodes)
    lambdas = np.sort(np.linalg.eigvals(lap).real)
    lambdas = np.clip(lambdas, 0.0, None)
    lam_max = float(lambdas[-1]) if n else 0.0
    nonzero = lambdas[lambdas > 1e-9]
    lam2 = float(nonzero[0]) if nonzero.size else 0.0
    nu_f = 1.0 / gamma

    # (i) per-mode slow rate vs diffusion rate
    disc = gamma * gamma - 4.0 * lambdas + 0j
    s_slow = ((-gamma + np.sqrt(disc)) / 2.0).real
    diff_rate = lambdas / gamma  # = nu_f * lambda_k
    mask = lambdas > 1e-9
    if np.any(mask):
        rel = np.abs(s_slow[mask] + diff_rate[mask]) / diff_rate[mask]
        max_rate_rel = float(np.max(rel))
    else:
        max_rate_rel = 0.0
    # Fiedler (slowest) mode.  Complex-safe: an under-damped gamma (gamma^2 <
    # 4*lam2, reachable when a caller fits gamma from oscillatory data) yields
    # a complex root whose real part -gamma/2 is the envelope decay rate.
    if lam2 > 0.0:
        s_gap = (
            (-gamma + np.sqrt(gamma * gamma - 4.0 * lam2 + 0j)) / 2.0
        ).real
        slow_gap = float(-s_gap)
    else:
        slow_gap = 0.0
    diff_gap = lam2 / gamma

    # (ii) trajectory: damped wave vs diffusion in the orthonormal eigenbasis
    sym_nodes, lsym = _symmetric_normalized_laplacian(G)
    w, V = np.linalg.eigh(lsym)
    w = np.clip(w, 0.0, None)
    q0 = structural_field(G, sym_nodes)
    c0 = V.T @ q0
    disc_w = gamma * gamma - 4.0 * w + 0j
    root_w = np.sqrt(disc_w)
    ss = (-gamma + root_w) / 2.0
    sf = (-gamma - root_w) / 2.0
    denom = sf - ss
    safe = np.abs(denom) > 1e-12
    a = np.where(safe, c0 * sf / np.where(safe, denom, 1.0), c0)
    b = np.where(safe, -c0 * ss / np.where(safe, denom, 1.0), 0.0)
    horizon = 3.0 * gamma / lam2 if lam2 > 0.0 else gamma
    ts = np.linspace(0.05 * horizon, horizon, max(2, n_time_samples))
    traj_err = 0.0
    for t in ts:
        q_wave = V @ (a * np.exp(ss * t) + b * np.exp(sf * t)).real
        q_diff = V @ (c0 * np.exp(-w * t / gamma))
        denom_t = float(np.linalg.norm(q_diff)) + 1e-12
        err = float(np.linalg.norm(q_wave - q_diff)) / denom_t
        traj_err = max(traj_err, err)

    projects = (max_rate_rel < tolerance) and (traj_err < tolerance)

    return OverdampedProjectionCertificate(
        n_nodes=n,
        gamma=float(gamma),
        nu_f_effective=float(nu_f),
        spectral_gap=lam2,
        lambda_max=lam_max,
        max_rate_rel_error=max_rate_rel,
        rate_error_times_gamma_sq=float(max_rate_rel * gamma * gamma),
        slowest_slow_rate=slow_gap,
        slowest_diffusion_rate=float(diff_gap),
        trajectory_max_rel_error=traj_err,
        projects_to_diffusion=bool(projects),
    )


@dataclass(frozen=True)
class UndampedLimitCertificate:
    r"""Verification that the γ→0 limit of the damped substrate wave is the
    undamped standing-wave spectrum √λ_k (the discrete modes).

    The overdamped projection (γ→∞, :func:`verify_overdamped_projection`)
    collapses the damped substrate wave onto structural diffusion.  Its
    opposite end, γ→0, is the **conservative** limit: the roots of
    s² + γs + λ_k = 0 become the pure-imaginary pair s = ±i√λ_k, so every
    mode oscillates undamped at the standing-wave frequency ω_k = √λ_k —
    exactly the discrete modes of :func:`verify_discrete_modes`.  γ is the
    single dial between the two regimes: γ→∞ diffusion, γ→0 standing waves.

    Attributes
    ----------
    n_nodes : int
    gamma : float
        Small damping used to probe the conservative limit.
    max_decay_rate : float
        Max |Re(s)| over the non-uniform modes (→ 0 as γ → 0; equals γ/2
        for underdamped modes — the envelope decay).
    max_freq_rel_error : float
        Max over modes of |Im(s) − √λ_k| / √λ_k: how close the damped
        oscillation frequency is to the undamped standing-wave frequency.
    freq_error_times_inv_gamma_sq : float
        ``max_freq_rel_error / γ²`` — converges to a constant, confirming
        the frequency error scales as O(γ²).
    matches_discrete_modes : bool
        Whether the γ→0 frequencies match the standing-wave spectrum √λ_k.
    standing_wave_frequencies : tuple[float, ...]
        The lowest few undamped frequencies ω_k = √λ_k (in ascending
        eigenvalue order, starting from the uniform mode ω₀ ≈ 0 — the same
        convention as :func:`verify_discrete_modes`).
    """

    n_nodes: int
    gamma: float
    max_decay_rate: float
    max_freq_rel_error: float
    freq_error_times_inv_gamma_sq: float
    matches_discrete_modes: bool
    standing_wave_frequencies: tuple[float, ...]

    @property
    def is_valid_undamped_limit(self) -> bool:
        """True when the γ→0 wave recovers the standing-wave spectrum."""
        return self.matches_discrete_modes

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_valid_undamped_limit else "INVALID"
        return (
            f"Undamped limit [{ok}]: gamma->0 damped wave -> standing waves "
            f"s=+-i*sqrt(lambda_k); at gamma={self.gamma:.3g} decay "
            f"{self.max_decay_rate:.3e} (=gamma/2), freq error "
            f"{self.max_freq_rel_error:.2e} (/gamma^2="
            f"{self.freq_error_times_inv_gamma_sq:.3f}), "
            f"matches discrete modes={self.matches_discrete_modes}"
        )


def verify_undamped_limit(
    G: Any,
    *,
    gamma: float = 1e-3,
    tolerance: float = 1e-2,
) -> UndampedLimitCertificate:
    r"""Verify the γ→0 limit of the damped substrate wave is standing waves.

    Forms the damped graph wave q̈ + γq̇ + L q = 0 at a *small* damping γ and
    measures that each underdamped mode's complex root tends to the
    pure-imaginary standing-wave value s = ±i√λ_k: the envelope decay
    Re(s) = −γ/2 → 0, and the oscillation frequency Im(s) → √λ_k (the
    discrete-mode frequency of :func:`verify_discrete_modes`).  This is the
    conservative end of the same γ-dial whose γ→∞ end is the overdamped
    projection onto structural diffusion.

    Parameters
    ----------
    G : TNFRGraph
    gamma : float
        Small damping (γ² < 4·λ₂ keeps the slow modes underdamped /
        oscillatory).
    tolerance : float
        Maximum relative frequency error for a valid undamped limit.

    Returns
    -------
    UndampedLimitCertificate
    """
    nodes, lap = structural_diffusion_operator(G)
    n = len(nodes)
    lambdas, s_slow, _ = damped_wave_rates(G, gamma)
    omega = np.sqrt(np.clip(lambdas, 0.0, None))  # standing-wave frequencies

    # complex roots of s^2 + gamma s + lambda = 0 (full, not just real part)
    disc = gamma * gamma - 4.0 * lambdas + 0j
    root = np.sqrt(disc)
    s_plus = (-gamma + root) / 2.0
    decay = np.abs(s_plus.real)  # envelope decay = gamma/2 (underdamped)
    freq = np.abs(s_plus.imag)  # oscillation frequency

    mask = lambdas > 1e-9
    if np.any(mask):
        rel = np.abs(freq[mask] - omega[mask]) / omega[mask]
        max_freq_rel = float(np.max(rel))
        max_decay = float(np.max(decay[mask]))
    else:
        max_freq_rel = 0.0
        max_decay = 0.0

    freqs = tuple(float(omega[k]) for k in range(min(6, n)))
    matches = max_freq_rel < tolerance

    return UndampedLimitCertificate(
        n_nodes=n,
        gamma=float(gamma),
        max_decay_rate=max_decay,
        max_freq_rel_error=max_freq_rel,
        freq_error_times_inv_gamma_sq=float(max_freq_rel / (gamma * gamma)),
        matches_discrete_modes=bool(matches),
        standing_wave_frequencies=freqs,
    )


# ---------------------------------------------------------------------------
# Discrete modes: the standing waves of the bounded structural manifold
# ---------------------------------------------------------------------------


def _symmetric_normalized_laplacian(G: Any) -> tuple[list, Any]:
    r"""Return the symmetric normalized Laplacian L_sym.

    L_sym = I − D^{-1/2} W D^{-1/2} is symmetric (orthonormal eigenvectors)
    and shares the spectrum of the random-walk diffusion operator L_rw used
    by ΔNFR.
    """
    nodes = _ordered_nodes(G)
    index = {n: i for i, n in enumerate(nodes)}
    n = len(nodes)
    deg = np.zeros(n, dtype=float)
    adj = np.zeros((n, n), dtype=float)
    for node in nodes:
        i = index[node]
        for m in G.neighbors(node):
            w = float(G[node][m].get("weight", 1.0))
            adj[i, index[m]] = w
            deg[i] += w
    # Compute D^{-1/2} only on connected nodes; isolated nodes (deg = 0) keep
    # dinv = 0 (no division by zero -- the np.where form still evaluates
    # 1/sqrt(0) for every entry before selecting, emitting a spurious warning).
    dinv = np.zeros(n, dtype=float)
    positive = deg > 0.0
    dinv[positive] = 1.0 / np.sqrt(deg[positive])
    lap = np.eye(n) - (dinv[:, None] * adj * dinv[None, :])
    lap = 0.5 * (lap + lap.T)  # symmetrise residual numerical asymmetry
    return nodes, lap


def _topology_signature(G: Any) -> tuple:
    """A hashable signature of the graph topology (nodes + weighted edges).

    The diffusion Laplacian depends only on the node set and the weighted edge
    set -- never on node state (EPI, nu_f, theta). When this signature is
    unchanged the spectrum is unchanged, so the eigendecomposition can be
    reused across evolution steps and across read-outs.
    """
    nodes = _ordered_nodes(G)
    idx = {nd: i for i, nd in enumerate(nodes)}
    edges = tuple(
        sorted(
            (idx[u], idx[v], float(d.get("weight", 1.0)))
            for u, v, d in G.edges(data=True)
        )
    )
    return (len(nodes), tuple(nodes), edges)


def _cached_eigh(G: Any) -> tuple[Any, Any]:
    """Eigendecomposition of L_sym, memoized on the topology signature.

    Returns ``(eigenvalues, eigenvectors)`` ascending and clipped ``>= 0``.
    The decomposition is invariant under evolution on a fixed graph, so the
    O(N^3) ``eigh`` runs once per topology; the cache lives in ``G.graph`` and
    self-invalidates when the signature changes. ``structural_eigenmodes`` (the
    rhythm path) and ``relaxation_spectrum`` (the spectrum / NFR path) share
    this one decomposition -- L_sym and L_rw have the same spectrum.
    """
    sig = _topology_signature(G)
    cache = G.graph.get("_tnfr_spectrum_cache")
    if isinstance(cache, dict) and cache.get("sig") == sig:
        return cache["vals"], cache["vecs"]
    _, lap = _symmetric_normalized_laplacian(G)
    eigvals, eigvecs = np.linalg.eigh(lap)
    eigvals = np.clip(eigvals, 0.0, None)
    G.graph["_tnfr_spectrum_cache"] = {
        "sig": sig,
        "vals": eigvals,
        "vecs": eigvecs,
    }
    return eigvals, eigvecs


def structural_eigenmodes(G: Any) -> tuple[Any, Any]:
    r"""Return the discrete eigenmodes of the bounded structural manifold.

    Computes the eigenvalues {λ_k} and orthonormal eigenvectors {v_k} of the
    symmetric normalized Laplacian L_sym (same spectrum as the diffusion
    operator L_rw).  The eigenvalues are the discrete mode "energies"; the
    eigenvectors are the standing-wave mode shapes (orthonormal), sorted by
    ascending λ_k.  λ_1 = 0 is the uniform mode.

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (eigenvalues, eigenvectors) : tuple[np.ndarray, np.ndarray]
        ``eigenvalues`` shape ``(N,)`` ascending; ``eigenvectors`` shape
        ``(N, N)`` with column ``k`` the k-th standing-wave mode shape.
    """
    return _cached_eigh(G)


def nodal_domain_count(mode: Any) -> int:
    r"""Number of sign changes (nodal domains − 1) of a standing-wave mode.

    The structural "mode number": the k-th standing wave has k sign changes
    on a 1D manifold (Courant's nodal-domain ordering).  Near-zero entries
    are ignored to avoid spurious sign flips.

    Parameters
    ----------
    mode : np.ndarray
        A mode shape (eigenvector).

    Returns
    -------
    int
        The number of sign changes along the node ordering.
    """
    v = np.asarray(mode, dtype=float)
    sig = np.sign(v[np.abs(v) > 1e-12])
    if sig.size < 2:
        return 0
    return int(np.sum(np.abs(np.diff(sig)) > 0))


def compute_emergent_pulse(G: Any, n_modes: int = 8) -> dict[str, Any]:
    r"""The emergent pulse: the rhythm the substrate plays [CANONICAL].

    The conservative face of the nodal dynamics is a *sustained vibration*:
    every structural mode oscillates at the standing-wave frequency
    :math:`\omega_k = \sqrt{\lambda_k}` (the discrete modes of
    :func:`verify_undamped_limit` / :func:`structural_eigenmodes`).  The rhythm
    is the interference of those resonances -- beats at the differences
    :math:`\omega_j - \omega_k` -- and the equilibria (the ``dNFR = 0``
    coherence states) are the beats the vibration passes through.  This unifies
    the scattered conservative machinery (the spectrum, the standing-wave
    frequencies, the self-similar decimation) into one *pulse* read-out,
    computed closed-form from the structural spectrum (no time integration) --
    the conservative twin of the dissipative coherence read-out.

    Parameters
    ----------
    G : TNFRGraph
    n_modes : int
        Number of leading resonant frequencies to report.

    Returns
    -------
    dict
        ``resonant_spectrum`` (leading :math:`\omega_k = \sqrt{\lambda_k}`),
        ``fundamental`` (the slowest non-uniform resonance), ``dominant_beat``
        (the slowest beat = smallest positive :math:`\omega_j - \omega_k`),
        ``spectral_multiplicity`` (the largest eigenvalue multiplicity = the
        self-similar / fractal signature), ``vibration_energy``
        (:math:`\tfrac12\sum\lambda_k`, the conserved structural-pressure
        energy of the vibration), ``n_modes``.
    """
    eigvals, _ = structural_eigenmodes(G)
    eigvals = np.asarray(eigvals, dtype=float)
    omega = np.sqrt(np.clip(eigvals, 0.0, None))
    nz = np.sort(omega[eigvals > 1e-9])  # exclude the uniform (lambda~0) mode
    beats = np.diff(nz) if nz.size > 1 else np.asarray([])
    pos = beats[beats > 1e-9] if beats.size else beats
    _, counts = np.unique(np.round(eigvals, 9), return_counts=True)
    return {
        "resonant_spectrum": [float(x) for x in nz[:n_modes]],
        "fundamental": float(nz[0]) if nz.size else 0.0,
        "dominant_beat": float(pos.min()) if pos.size else 0.0,
        "spectral_multiplicity": int(counts.max()),
        "vibration_energy": float(0.5 * np.sum(eigvals)),
        "n_modes": int(nz.size),
    }


def compute_nodal_pulse(G: Any) -> dict[str, Any]:
    r"""The per-NFR pulse and its resonance [CANONICAL].

    The collective rhythm (:func:`compute_emergent_pulse`) is what the NFR
    *bricks* produce: every NFR is itself a phase oscillator -- the
    single-node reduction of the nodal equation
    :math:`\partial\mathrm{EPI}_i/\partial t=\nu_{f,i}\,\Delta\mathrm{NFR}_i`
    -- pulsing at its own structural frequency :math:`\nu_{f,i}` with phase
    :math:`\varphi_i`.  *Resonance* couples those pulses: the local
    phase-synchrony (:func:`~tnfr.metrics.coherence.local_phase_sync`)
    measures how phase-locked each NFR is with its neighbours, the global
    Kuramoto order ``R`` (:func:`~tnfr.gamma.kuramoto_R_psi`) the collective
    locking, and the U3 gate :data:`~tnfr.constants.canonical.DELTA_PHI_MAX`
    sets admissibility.  The collective pulse emerges as these per-NFR pulses
    resonate (``R -> 1``).  This is the *local* face of the rhythm -- the
    pulsing NFRs that generate the network rhythm -- read from canonical
    per-node quantities (no time integration).

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    dict
        ``mean_frequency`` / ``frequency_spread`` (the per-NFR pulse rates
        nu_f: mean and std in Hz_str), ``phase_coherence`` (the collective
        Kuramoto ``R`` in ``[0, 1]`` -- how locked the pulses are),
        ``mean_local_resonance`` (mean per-NFR local phase synchrony in
        ``[0, 1]``), ``resonance_gate`` (the U3 admissibility bound
        Delta phi_max), ``n_pulsing`` (NFRs with nu_f > 0), ``n_nodes``.
    """
    from ..constants.canonical import DELTA_PHI_MAX

    nodes = list(G.nodes())
    n = len(nodes)
    gate = float(DELTA_PHI_MAX)
    if not n:
        return {
            "mean_frequency": 0.0,
            "frequency_spread": 0.0,
            "phase_coherence": 0.0,
            "mean_local_resonance": 0.0,
            "resonance_gate": gate,
            "n_pulsing": 0,
            "n_nodes": 0,
        }
    vf = np.asarray(
        [float(get_attr(G.nodes[k], ALIAS_VF, 0.0) or 0.0) for k in nodes],
        dtype=float,
    )
    # collective resonance of the per-NFR pulses (Kuramoto order parameter)
    try:
        from ..gamma import kuramoto_R_psi

        phase_coherence = float(kuramoto_R_psi(G)[0])
    except Exception:
        phase_coherence = 0.0
    # per-NFR resonance: mean local phase synchrony (one shared matrix build)
    try:
        from ..metrics.coherence import (
            coherence_matrix,
            local_phase_sync_weighted,
        )

        order, W = coherence_matrix(G)
        if order is None:
            mean_local = 0.0
        else:
            mean_local = float(
                np.mean(
                    [
                        local_phase_sync_weighted(
                            G, k, nodes_order=order, W_row=W
                        )
                        for k in order
                    ]
                )
            )
    except Exception:
        mean_local = 0.0
    return {
        "mean_frequency": float(np.mean(vf)),
        "frequency_spread": float(np.std(vf)),
        "phase_coherence": phase_coherence,
        "mean_local_resonance": mean_local,
        "resonance_gate": gate,
        "n_pulsing": int(np.sum(vf > 1e-9)),
        "n_nodes": n,
    }


@dataclass(frozen=True)
class DiscreteModeCertificate:
    r"""Verification of the discrete standing-wave modes of a bounded manifold.

    A bounded structural manifold (finite graph) supports a discrete
    spectrum of orthonormal standing-wave eigenmodes — the same structure
    as the discrete harmonics of a vibrating string (Pythagoras), a Chladni
    plate, or a molecular vibrational spectrum.

    Attributes
    ----------
    n_modes : int
        Number of discrete modes (= number of nodes; finite/discrete).
    spectrum_is_discrete : bool
        The manifold has a finite, discrete eigenvalue spectrum.
    modes_orthonormal : bool
        The standing-wave mode shapes are orthonormal.
    max_orthonormality_residual : float
        Max |⟨v_i, v_j⟩ − δ_ij| (≈ 0).
    has_uniform_zero_mode : bool
        λ_1 = 0 (the uniform mode / conserved diffusion mode).
    spectral_gap : float
        λ_2 — the first non-trivial mode.
    matches_diffusion_spectrum : bool
        The L_sym spectrum equals the diffusion operator (L_rw) spectrum.
    nodal_domains_grow : bool
        The nodal-domain count grows from the lowest to the highest mode
        (Courant ordering; structural mode number).
    standing_wave_frequencies : tuple
        The first few standing-wave frequencies ω_k = √λ_k.
    """

    n_modes: int
    spectrum_is_discrete: bool
    modes_orthonormal: bool
    max_orthonormality_residual: float
    has_uniform_zero_mode: bool
    spectral_gap: float
    matches_diffusion_spectrum: bool
    nodal_domains_grow: bool
    standing_wave_frequencies: tuple

    @property
    def is_valid_discrete_modes(self) -> bool:
        """True when the manifold verifies as discrete standing waves."""
        return (
            self.spectrum_is_discrete
            and self.modes_orthonormal
            and self.has_uniform_zero_mode
            and self.matches_diffusion_spectrum
            and self.nodal_domains_grow
        )

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_valid_discrete_modes else "INVALID"
        freqs = ", ".join(f"{f:.3f}" for f in self.standing_wave_frequencies)
        return (
            f"Discrete standing-wave modes [{ok}]: "
            f"{self.n_modes} discrete modes, "
            f"orthonormal={self.modes_orthonormal} "
            f"(res {self.max_orthonormality_residual:.1e}), "
            f"uniform λ₁=0={self.has_uniform_zero_mode}, "
            f"spectral gap λ₂={self.spectral_gap:.4f}, "
            f"matches diffusion spectrum={self.matches_diffusion_spectrum}, "
            f"nodal domains grow={self.nodal_domains_grow}; "
            f"ω_k=√λ_k=[{freqs}]"
        )


def verify_discrete_modes(
    G: Any, *, tolerance: float = 1e-9
) -> DiscreteModeCertificate:
    r"""Verify the discrete standing-wave modes of the bounded manifold.

    Confirms that the finite manifold has a discrete spectrum of orthonormal
    standing-wave eigenmodes, with a uniform λ_1 = 0 mode, a spectrum
    matching the diffusion operator (L_rw), and nodal-domain counts growing
    with the mode index (Courant) — the structural origin of "discrete
    modes", the same as the discrete harmonics of a bounded elastic medium.

    Parameters
    ----------
    G : TNFRGraph
    tolerance : float
        Numerical tolerance for the orthonormality / spectrum checks.

    Returns
    -------
    DiscreteModeCertificate
    """
    eigvals, eigvecs = structural_eigenmodes(G)
    n = len(eigvals)

    discrete = n > 0 and np.all(np.isfinite(eigvals))

    gram = eigvecs.T @ eigvecs
    ortho_res = float(np.max(np.abs(gram - np.eye(n)))) if n else 0.0
    orthonormal = ortho_res < max(tolerance, 1e-9)

    uniform_zero = bool(abs(float(eigvals[0])) < 1e-6) if n else False
    gap = float(eigvals[1]) if n > 1 else 0.0

    # spectrum matches the random-walk diffusion operator L_rw
    _, lrw = structural_diffusion_operator(G)
    rw_spec = np.sort(np.linalg.eigvals(lrw).real)
    matches = bool(np.allclose(np.sort(eigvals), rw_spec, atol=1e-7))

    # nodal-domain counts grow from lowest to highest mode (Courant)
    counts = [nodal_domain_count(eigvecs[:, k]) for k in range(n)]
    grow = (counts[0] == 0 and counts[-1] > counts[0]) if n > 1 else True

    freqs = tuple(float(np.sqrt(eigvals[k])) for k in range(min(6, n)))

    return DiscreteModeCertificate(
        n_modes=n,
        spectrum_is_discrete=discrete,
        modes_orthonormal=orthonormal,
        max_orthonormality_residual=ortho_res,
        has_uniform_zero_mode=uniform_zero,
        spectral_gap=gap,
        matches_diffusion_spectrum=matches,
        nodal_domains_grow=grow,
        standing_wave_frequencies=freqs,
    )


# ---------------------------------------------------------------------------
# Structural stability: the dispersion relation and the instability threshold
# ---------------------------------------------------------------------------


def dispersion_relation(G: Any, reaction_rate: float = 0.0) -> Any:
    r"""Per-mode growth rates σ_k = reaction_rate − νf·λ_k (dispersion).

    The growth (σ_k > 0) or decay (σ_k < 0) rate of each structural
    eigenmode under diffusion plus a local reaction rate.  At
    ``reaction_rate = 0`` this is the negative of the relaxation spectrum
    (pure diffusion: every non-uniform mode decays).

    Parameters
    ----------
    G : TNFRGraph
    reaction_rate : float, optional
        A local growth/decay rate r added to every mode (the operators
        supply it: stabilizers lower r, destabilizers raise it).

    Returns
    -------
    np.ndarray
        The growth rates σ_k sorted by ascending λ_k.
    """
    eigvals, _ = structural_eigenmodes(G)
    nu = structural_diffusivity(G)
    return float(reaction_rate) - nu * eigvals


def instability_threshold(G: Any) -> float:
    r"""Structural instability threshold r_c = νf·λ_2 (spectral gap).

    The reaction rate above which the first *non-uniform* (Fiedler) mode
    becomes unstable.  For 0 < r < r_c only the uniform mode grows (global
    amplification); for r > r_c a structural pattern (the Fiedler partition)
    emerges.

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    float
        νf·λ_2 (0 if the graph has fewer than two modes).
    """
    eigvals, _ = structural_eigenmodes(G)
    nu = structural_diffusivity(G)
    return float(nu * eigvals[1]) if len(eigvals) > 1 else 0.0


def fiedler_partition(G: Any) -> tuple[list, list]:
    r"""The Fiedler-mode 2-partition — the first structural pattern.

    Splits the nodes by the sign of the Fiedler eigenvector (the mode with
    the smallest non-zero λ).  This is the network's **weakest structural
    cut** (the two most weakly-connected communities) — the empirically-
    validated spectral-clustering partition, and the first pattern to grow
    once the reaction rate crosses the instability threshold.

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (part_a, part_b) : tuple[list, list]
        Node lists for the two structural communities.
    """
    nodes = _ordered_nodes(G)
    eigvals, eigvecs = structural_eigenmodes(G)
    if eigvecs.shape[1] < 2:
        return list(nodes), []
    fiedler = eigvecs[:, 1]
    part_a = [nodes[i] for i in range(len(nodes)) if fiedler[i] > 0]
    part_b = [nodes[i] for i in range(len(nodes)) if fiedler[i] <= 0]
    return part_a, part_b


@dataclass(frozen=True)
class StructuralStabilityCertificate:
    r"""Verification of the structural linear-stability / dispersion relation.

    The dispersion relation σ_k = r − νf·λ_k governs the growth/decay of
    every structural eigenmode.  Pure diffusion (r = 0) decays every
    non-uniform mode (stable equilibrium); the threshold r_c = νf·λ_2
    separates uniform amplification from structural pattern formation (the
    Fiedler partition).  This is the spectral form of U2 grammar.

    Attributes
    ----------
    n_nodes : int
    diffusivity : float
        νf (the diffusion coefficient).
    spectral_gap : float
        λ_2 (the Fiedler value).
    instability_threshold : float
        r_c = νf·λ_2.
    pure_diffusion_stable : bool
        At r = 0 every non-uniform mode decays (σ_k < 0 for k ≥ 1).
    max_nonuniform_growth_at_zero : float
        Max σ_k over k ≥ 1 at r = 0 (should be < 0).
    dispersion_matches_relaxation : bool
        The r = 0 dispersion equals the negative relaxation spectrum.
    first_unstable_mode : int
        The first non-uniform mode to go unstable above threshold (= 1,
        the Fiedler mode).
    fiedler_partition_sizes : tuple
        Sizes (|A|, |B|) of the Fiedler 2-partition (the first pattern).
    """

    n_nodes: int
    diffusivity: float
    spectral_gap: float
    instability_threshold: float
    pure_diffusion_stable: bool
    max_nonuniform_growth_at_zero: float
    dispersion_matches_relaxation: bool
    first_unstable_mode: int
    fiedler_partition_sizes: tuple

    @property
    def is_valid_stability(self) -> bool:
        """True when the structural-stability picture verifies."""
        return (
            self.pure_diffusion_stable
            and self.dispersion_matches_relaxation
            and self.first_unstable_mode == 1
        )

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_valid_stability else "INVALID"
        a, b = self.fiedler_partition_sizes
        return (
            f"Structural stability [{ok}]: "
            f"νf={self.diffusivity:.4f}, spectral gap λ₂="
            f"{self.spectral_gap:.4f}, "
            f"instability threshold r_c=νf·λ₂="
            f"{self.instability_threshold:.4f}, "
            f"pure diffusion stable={self.pure_diffusion_stable} "
            f"(max non-uniform growth "
            f"{self.max_nonuniform_growth_at_zero:.1e}), "
            f"first unstable mode k={self.first_unstable_mode} (Fiedler), "
            f"Fiedler partition {a}/{b}"
        )


def verify_structural_stability(
    G: Any, *, tolerance: float = 1e-9
) -> StructuralStabilityCertificate:
    r"""Verify the structural linear-stability / dispersion relation.

    Confirms that pure diffusion decays every non-uniform mode (stable
    equilibrium), that the r = 0 dispersion equals the negative relaxation
    spectrum, that the instability threshold is r_c = νf·λ_2, and that the
    first structural mode to go unstable above threshold is the Fiedler
    mode (whose eigenvector gives the weakest-cut partition).

    Parameters
    ----------
    G : TNFRGraph
    tolerance : float
        Numerical tolerance.

    Returns
    -------
    StructuralStabilityCertificate
    """
    eigvals, eigvecs = structural_eigenmodes(G)
    n = len(eigvals)
    nu = structural_diffusivity(G)
    gap = float(eigvals[1]) if n > 1 else 0.0
    r_c = nu * gap

    # pure diffusion (r=0): non-uniform modes (k>=1) decay
    sigma0 = dispersion_relation(G, 0.0)
    max_nonuniform = float(np.max(sigma0[1:])) if n > 1 else 0.0
    stable = max_nonuniform < tolerance

    # dispersion at r=0 == negative relaxation spectrum
    rates = relaxation_spectrum(G)
    matches = bool(np.allclose(np.sort(-sigma0), np.sort(rates), atol=1e-7))

    # first non-uniform mode to go unstable just above threshold
    sigma = dispersion_relation(G, r_c + max(1e-3, 1e-3 * r_c))
    if n > 1 and np.any(sigma[1:] > tolerance):
        first_unstable = int(np.argmax(sigma[1:] > tolerance)) + 1
    else:
        first_unstable = 0

    # Fiedler partition sizes
    part_a, part_b = fiedler_partition(G)
    sizes = (len(part_a), len(part_b))

    return StructuralStabilityCertificate(
        n_nodes=n,
        diffusivity=nu,
        spectral_gap=gap,
        instability_threshold=r_c,
        pure_diffusion_stable=stable,
        max_nonuniform_growth_at_zero=max_nonuniform,
        dispersion_matches_relaxation=matches,
        first_unstable_mode=first_unstable,
        fiedler_partition_sizes=sizes,
    )


# ---------------------------------------------------------------------------
# The structural random walk: Brownian motion and resistance geometry
# ---------------------------------------------------------------------------


def _adjacency_degree(G: Any) -> tuple[list, Any, Any]:
    """Return (nodes, weighted adjacency W, weighted degree vector)."""
    nodes = _ordered_nodes(G)
    index = {n: i for i, n in enumerate(nodes)}
    n = len(nodes)
    adj = np.zeros((n, n), dtype=float)
    for node in nodes:
        i = index[node]
        for m in G.neighbors(node):
            adj[i, index[m]] = float(G[node][m].get("weight", 1.0))
    deg = adj.sum(axis=1)
    return nodes, adj, deg


def random_walk_matrix(G: Any) -> tuple[list, Any]:
    r"""Return the random-walk transition matrix P = D^{-1}W.

    The diffusion operator is L_rw = I − P, so P is exactly the
    random-walk the structural diffusion generates.  Row-stochastic (each
    row sums to 1); isolated nodes get a zero row.

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (nodes, P) : tuple[list, np.ndarray]
    """
    nodes, adj, deg = _adjacency_degree(G)
    with np.errstate(divide="ignore", invalid="ignore"):
        dinv = np.where(deg > 0.0, 1.0 / deg, 0.0)
    return nodes, dinv[:, None] * adj


def stationary_distribution(G: Any) -> tuple[list, Any]:
    r"""Stationary distribution π_i = deg(i) / Σ deg of the random walk.

    The random walk converges to the degree distribution — exactly the
    degree-weighted total the diffusion conserves.  The conserved quantity
    is the equilibrium measure.

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (nodes, pi) : tuple[list, np.ndarray]
    """
    nodes, _, deg = _adjacency_degree(G)
    total = float(deg.sum())
    pi = deg / total if total > 0 else deg
    return nodes, pi


def _laplacian_pinv(G: Any) -> tuple[list, Any, float]:
    """Return (nodes, pinv(L) of the combinatorial Laplacian, n_edges)."""
    nodes, adj, deg = _adjacency_degree(G)
    lap = np.diag(deg) - adj
    return nodes, np.linalg.pinv(lap), float(deg.sum()) / 2.0


def effective_resistance(G: Any) -> tuple[list, Any]:
    r"""Effective-resistance matrix R_eff(i,j) (Ohm's law).

    Treating the network as a resistor network (the combinatorial
    Laplacian L = D − W is the conductance matrix — Kirchhoff 1847), the
    effective resistance between nodes is

        R_eff(i,j) = L⁺_ii + L⁺_jj − 2·L⁺_ij,

    with L⁺ the Moore–Penrose pseudoinverse.  R_eff is a metric
    (symmetric, non-negative, triangle inequality) — the structural
    "difficulty of transport" between two nodes.

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (nodes, R) : tuple[list, np.ndarray]
        ``R[i, j]`` is the effective resistance between node i and node j.
    """
    nodes, lp, _ = _laplacian_pinv(G)
    diag = np.diag(lp)
    r = diag[:, None] + diag[None, :] - 2.0 * lp
    np.fill_diagonal(r, 0.0)
    return nodes, np.maximum(r, 0.0)


def commute_time(G: Any) -> tuple[list, Any]:
    r"""Commute-time matrix C(i,j) = 2m·R_eff(i,j).

    The expected round-trip time of the structural random walk between two
    nodes equals 2m times the effective resistance (m the number of
    edges) — the exact link between the diffusion random walk and the
    resistance geometry (Chandra et al. 1996).

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (nodes, C) : tuple[list, np.ndarray]
    """
    nodes, r = effective_resistance(G)
    _, _, deg = _adjacency_degree(G)
    m_edges = float(deg.sum()) / 2.0
    return nodes, 2.0 * m_edges * r


@dataclass(frozen=True)
class RandomWalkCertificate:
    r"""Verification of the structural random walk and resistance geometry.

    The diffusion operator is the generator of a random walk (Brownian
    motion on the network); its stationary distribution is the degree, and
    the effective resistance / commute time give the transport geometry.

    Attributes
    ----------
    n_nodes : int
    operator_is_walk_generator : bool
        L_rw = I − P exactly (the diffusion operator generates the walk).
    transition_row_stochastic : bool
        P = D⁻¹W is row-stochastic.
    stationary_is_degree : bool
        π = deg / Σ deg is the stationary distribution (π·P = π).
    resistance_is_metric : bool
        R_eff is symmetric, non-negative, and obeys the triangle
        inequality.
    max_resistance : float
        The largest pairwise effective resistance (transport diameter).
    commute_equals_2m_resistance : bool
        C(i,j) = 2m·R_eff(i,j) (random walk ↔ resistance identity).
    max_walk_generator_residual : float
        Max |L_rw − (I − P)| (≈ 0).
    """

    n_nodes: int
    operator_is_walk_generator: bool
    transition_row_stochastic: bool
    stationary_is_degree: bool
    resistance_is_metric: bool
    max_resistance: float
    commute_equals_2m_resistance: bool
    max_walk_generator_residual: float

    @property
    def is_valid_random_walk(self) -> bool:
        """True when the structural random-walk picture verifies."""
        return (
            self.operator_is_walk_generator
            and self.transition_row_stochastic
            and self.stationary_is_degree
            and self.resistance_is_metric
            and self.commute_equals_2m_resistance
        )

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_valid_random_walk else "INVALID"
        return (
            f"Structural random walk [{ok}]: "
            f"L_rw = I−P={self.operator_is_walk_generator} "
            f"(res {self.max_walk_generator_residual:.1e}), "
            f"P row-stochastic={self.transition_row_stochastic}, "
            f"stationary π=degree={self.stationary_is_degree}, "
            f"R_eff metric={self.resistance_is_metric} "
            f"(max R {self.max_resistance:.4f}), "
            f"commute=2m·R_eff={self.commute_equals_2m_resistance}"
        )


def verify_structural_random_walk(
    G: Any, *, tolerance: float = 1e-9
) -> RandomWalkCertificate:
    r"""Verify the structural random walk and resistance geometry.

    Confirms that the diffusion operator is the random-walk generator
    (L_rw = I − P), that P is row-stochastic with stationary distribution
    π = degree, that the effective resistance is a transport metric, and
    that the commute time equals 2m·R_eff (random walk ↔ resistance).

    Parameters
    ----------
    G : TNFRGraph
    tolerance : float
        Numerical tolerance.

    Returns
    -------
    RandomWalkCertificate
    """
    nodes, lrw = structural_diffusion_operator(G)
    n = len(nodes)
    _, p = random_walk_matrix(G)

    # L_rw = I − P
    gen_res = float(np.max(np.abs(lrw - (np.eye(n) - p)))) if n else 0.0
    is_generator = gen_res < max(tolerance, 1e-12)

    # P row-stochastic (rows of connected nodes sum to 1)
    row_sums = p.sum(axis=1)
    deg_nonzero = np.array([sum(1 for _ in G.neighbors(nd)) > 0 for nd in nodes])
    row_stochastic = (
        bool(np.all(np.abs(row_sums[deg_nonzero] - 1.0) < tolerance)) if n else True
    )

    # stationary distribution π = degree, π·P = π
    _, pi = stationary_distribution(G)
    stationary_ok = bool(np.allclose(pi @ p, pi, atol=1e-7)) if n else True

    # effective resistance is a metric
    _, r = effective_resistance(G)
    symmetric = bool(np.allclose(r, r.T))
    nonneg = bool(np.all(r >= -1e-9))
    # triangle inequality on a sample of triples
    triangle = True
    if n >= 3:
        rng = np.random.default_rng(0)
        for _ in range(200):
            i, j, k = rng.integers(0, n, size=3)
            if r[i, k] > r[i, j] + r[j, k] + 1e-7:
                triangle = False
                break
    is_metric = symmetric and nonneg and triangle
    max_r = float(np.max(r)) if n else 0.0

    # commute time = 2m·R_eff
    _, c = commute_time(G)
    _, _, deg = _adjacency_degree(G)
    m_edges = float(deg.sum()) / 2.0
    commute_ok = bool(np.allclose(c, 2.0 * m_edges * r, atol=1e-7))

    return RandomWalkCertificate(
        n_nodes=n,
        operator_is_walk_generator=is_generator,
        transition_row_stochastic=row_stochastic,
        stationary_is_degree=stationary_ok,
        resistance_is_metric=is_metric,
        max_resistance=max_r,
        commute_equals_2m_resistance=commute_ok,
        max_walk_generator_residual=gen_res,
    )


# ---------------------------------------------------------------------------
# The structural flow: current, Kirchhoff's law, and continuity
# ---------------------------------------------------------------------------
def structural_current(G: Any) -> tuple[list, Any]:
    r"""Structural (diffusion) edge-current matrix J_ij = EPI_i − EPI_j.

    The transport carries a current: along each edge i∼j the diffusion
    flux is J_ij = EPI_i − EPI_j (Fick's law — flux from high to low
    concentration).  The matrix is **antisymmetric** (J_ij = −J_ji) and
    supported on edges only (J_ij = 0 when i and j are not adjacent).

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (nodes, J) : tuple[list, np.ndarray]
        ``J[i, j]`` is the current from node i to node j across edge i∼j.
    """
    from ..alias import get_attr
    from ..constants.aliases import ALIAS_EPI

    nodes, adj, _ = _adjacency_degree(G)
    epi = np.array(
        [float(get_attr(G.nodes[n], ALIAS_EPI, 0.0)) for n in nodes],
        dtype=float,
    )
    mask = adj != 0.0
    # J_ij = EPI_i − EPI_j on edges, zero elsewhere
    j = (epi[:, None] - epi[None, :]) * mask
    return nodes, j


def current_divergence(G: Any) -> tuple[list, Any]:
    r"""Net current outflow per node div(J)(i) = Σ_{j∼i} J_ij = (L·EPI)(i).

    Kirchhoff's current law: the net outflow at a node equals the
    combinatorial Laplacian L = D − W acting on EPI.  This is the discrete
    continuity equation div(J) = L·EPI, so ∂EPI/∂t + div(J) = 0 for the
    diffusion dynamics (∂EPI/∂t = −L_rw·EPI carries the same content
    degree-normalized).

    Parameters
    ----------
    G : TNFRGraph

    Returns
    -------
    (nodes, div) : tuple[list, np.ndarray]
        ``div[i]`` is the net current leaving node i.
    """
    nodes, j = structural_current(G)
    return nodes, j.sum(axis=1)


@dataclass(frozen=True)
class StructuralFlowCertificate:
    r"""Verification of the structural flow (current, Kirchhoff, Ohm).

    The diffusion transport carries a structural current J_ij = EPI_i −
    EPI_j (Fick's law).  Its node balance is Kirchhoff's current law —
    the discrete continuity equation div(J) = L·EPI — and under an
    injected current the potential drop is the effective resistance
    (Ohm's law).

    Attributes
    ----------
    n_nodes : int
    current_antisymmetric : bool
        J_ij = −J_ji (the current is a directed edge flow).
    kirchhoff_holds : bool
        Net outflow Σ_j J_ij = (L·EPI)_i (current law = continuity).
    total_flux_balances : bool
        Σ_i div(J)_i = 0 (closed network, no sources/sinks).
    equilibrium_zero_current : bool
        A uniform EPI field produces zero current everywhere.
    ohm_law_holds : bool
        An injected unit current s→t induces a potential drop equal to
        R_eff(s,t).
    max_kirchhoff_residual : float
        Max |Σ_j J_ij − (L·EPI)_i| (≈ 0).
    """

    n_nodes: int
    current_antisymmetric: bool
    kirchhoff_holds: bool
    total_flux_balances: bool
    equilibrium_zero_current: bool
    ohm_law_holds: bool
    max_kirchhoff_residual: float

    @property
    def is_valid_flow(self) -> bool:
        """True when the structural-flow picture verifies."""
        return (
            self.current_antisymmetric
            and self.kirchhoff_holds
            and self.total_flux_balances
            and self.equilibrium_zero_current
            and self.ohm_law_holds
        )

    def summary(self) -> str:
        """Human-readable one-line verdict."""
        ok = "VALID" if self.is_valid_flow else "INVALID"
        return (
            f"Structural flow [{ok}]: "
            f"current antisymmetric={self.current_antisymmetric}, "
            f"Kirchhoff div(J)=L·EPI={self.kirchhoff_holds} "
            f"(res {self.max_kirchhoff_residual:.1e}), "
            f"total flux balances={self.total_flux_balances}, "
            f"equilibrium zero current={self.equilibrium_zero_current}, "
            f"Ohm drop=R_eff={self.ohm_law_holds}"
        )


def verify_structural_flow(
    G: Any, *, tolerance: float = 1e-9
) -> StructuralFlowCertificate:
    r"""Verify the structural flow: current, Kirchhoff's law, Ohm's law.

    Confirms that the diffusion edge current J_ij = EPI_i − EPI_j is
    antisymmetric, that Kirchhoff's current law div(J) = L·EPI holds (the
    discrete continuity equation), that the total flux balances on a closed
    network, that a uniform EPI field carries zero current, and that an
    injected unit current induces a potential drop equal to the effective
    resistance (Ohm's law).

    Parameters
    ----------
    G : TNFRGraph
    tolerance : float
        Numerical tolerance.

    Returns
    -------
    StructuralFlowCertificate
    """
    from ..alias import get_attr
    from ..constants.aliases import ALIAS_EPI

    nodes, j = structural_current(G)
    n = len(nodes)

    # current antisymmetry J_ij = −J_ji
    antisym = bool(np.allclose(j, -j.T, atol=tolerance)) if n else True

    # Kirchhoff: net outflow = (L·EPI) with L the combinatorial Laplacian
    _, adj, deg = _adjacency_degree(G)
    lap = np.diag(deg) - adj
    epi = np.array(
        [float(get_attr(G.nodes[nd], ALIAS_EPI, 0.0)) for nd in nodes],
        dtype=float,
    )
    net_out = j.sum(axis=1)
    kirchhoff_res = float(np.max(np.abs(net_out - lap @ epi))) if n else 0.0
    kirchhoff_ok = kirchhoff_res < max(tolerance, 1e-9)

    # total flux balances (L has zero column sums) — closed network
    total_balances = (
        bool(abs(float(net_out.sum())) < max(tolerance, 1e-9)) if n else True
    )

    # equilibrium: a uniform EPI field carries zero current
    uniform = np.ones(n)
    j_uniform = (uniform[:, None] - uniform[None, :]) * (adj != 0.0)
    equilibrium_ok = bool(np.allclose(j_uniform, 0.0, atol=tolerance)) if n else True

    # Ohm's law: injected unit current s→t induces drop V_s − V_t = R_eff
    ohm_ok = True
    if n >= 2:
        lp = np.linalg.pinv(lap)
        diag = np.diag(lp)
        rng = np.random.default_rng(0)
        for _ in range(min(20, n)):
            s, t = rng.integers(0, n, size=2)
            if s == t:
                continue
            b = np.zeros(n)
            b[s], b[t] = 1.0, -1.0
            v = lp @ b
            drop = v[s] - v[t]
            r_eff = diag[s] + diag[t] - 2.0 * lp[s, t]
            if not np.isclose(drop, r_eff, atol=1e-7):
                ohm_ok = False
                break

    return StructuralFlowCertificate(
        n_nodes=n,
        current_antisymmetric=antisym,
        kirchhoff_holds=kirchhoff_ok,
        total_flux_balances=total_balances,
        equilibrium_zero_current=equilibrium_ok,
        ohm_law_holds=ohm_ok,
        max_kirchhoff_residual=kirchhoff_res,
    )