r"""Structural diffusion and scoped diagnostics for the nodal EPI channel. For a fixed, nonnegative adjacency W, let D contain its row strengths and L_rw = I - D^-1 W, with a zero row at an isolated node. The canonical EPI pressure is DeltaNFR_epi = -L_rw EPI. With the capacity vector nu_f, the isolated channel therefore evolves as dEPI/dt = -A EPI, A = diag(nu_f) L_rw. The full pressure also includes phase, frequency and topology channels. Wrapped circular phase averaging is nonlinear and is not this scalar Laplacian identity. These diagnostics hold the graph and capacity fixed; they do not integrate a changing canonical operator sequence. Transport and spectral scope ---------------------------- For symmetric adjacency and a common positive capacity nu, each mode decays at nu*lambda_k. A connected graph relaxes to a uniform field and conserves sum(d_i*EPI_i). With disconnected components, stationarity need not be globally uniform. With positive heterogeneous capacity the conserved weights are d_i/nu_i and the decay rates are eigenvalues of A; replacing capacity by its mean changes the dynamics. Zero capacity freezes a node even when its pressure is nonzero. Directed transport uses the actual nonsymmetric A; real eigenvalue parts describe asymptotic damping, without certifying normality or absence of transient growth. The symmetric normalized Laplacian provides an orthonormal geometry basis only for symmetric adjacency. Its zero mode is proportional to sqrt(d) on a connected component, corresponding to a uniform EPI field after the degree-coordinate transformation. Finite dimension gives a finite spectrum; a nodal-domain upper bound does not imply monotonic domain counts for every graph or every basis of a degenerate eigenspace. Adding a scalar reaction rate r gives real growth rates r - Re(eigenvalues(A)). On a connected homogeneous network the first nonuniform threshold is nu*lambda_2. This is a spatial-mode threshold, not a bound on all evolution: the uniform mode already grows for any r>0. The correspondence is a linear diagnostic, not a derivation of grammar U2 for arbitrary sequences. Drift, waves and certificates ---------------------------- Under held pressure F the nodal equation is the first-order mobility law q_dot = nu_f*F. The damped graph-wave model q_ddot + gamma*q_dot + L*q = 0 has a slow diffusion limit with mobility 1/gamma. The wave-to-diffusion checks below concern that specified graph-wave model. They do not prove that the isotropic Hamiltonian implemented by symplectic_substrate, or all 13 engine operators, induces that graph wave. verify_structural_diffusion checks the canonical pressure on a replica and samples the actual frozen-capacity flow. It reports global uniformity and degree-weighted conservation separately from stationarity and the complete left-nullspace invariants of A. Finite-time residuals are diagnostics, not proofs of eventual convergence. Random walks and currents ------------------------ P = I - L_rw is row-stochastic, with an absorbing self-transition at an isolated node. For symmetric adjacency the degree distribution is stationary; convergence of a discrete walk additionally requires aperiodicity within an irreducible component. Resistance geometry requires symmetric nonnegative conductance. Disconnected pairs have infinite resistance and commute time; finite commute times use the volume of the pair's component. For symmetric W, edge currents J_ij = W_ij*(EPI_i-EPI_j) have divergence (D-W)*EPI. At nodes with positive degree and capacity, the nodal continuity balance is (d_i/nu_i)*dEPI_i/dt + div(J)_i = 0. This EPI-channel identity is distinct from the tetrad charge and currents in physics.conservation. The same conductance model gives the Dirichlet energy E_D = (1/4)*sum_ij W_ij*(EPI_i-EPI_j)^2. Its gradient is (D-W)*EPI, and the mobility is diag(nu_i/d_i), zero at isolates. Hence the frozen EPI-only flow satisfies dE_D/dt = -sum_i mobility_i*gradient_i^2 <= 0. compute_diffusion_energy reports this balance without evolving the graph. This energy is distinct from the tetrad potential in physics.variational. References within the implementation: dynamics.dnfr (canonical pressure), directed_diffusion (nonsymmetric transport), symplectic_substrate (specified Hamiltonian), and physics.conservation (tetrad diagnostics). """ from __future__ import annotations import math import struct import sys from collections.abc import Mapping from dataclasses import dataclass, field from fractions import Fraction from functools import lru_cache from typing import Any from ..alias import get_attr from ..constants.aliases import ALIAS_DNFR, ALIAS_EPI, ALIAS_VF from ..mathematics._weight_normalization import normalize_weights from ..mathematics.unified_numerical import np from ..types import real_scalar_epi from ..utils._structural_signature import ( proof_stamps_are_identical, structural_proof_signature, ) from ._conductance import ConductanceSnapshot, read_conductance from ._exact_linear_algebra import exact_matrix_inverse as _exact_matrix_inverse from ._exact_linear_algebra import ( exact_symmetric_semidefinite as _exact_symmetric_semidefinite, ) from ._helpers import finite_real_scalar __all__ = [ "DiffusionEnergyBalance", "HeterogeneousDiffusionStabilityCertificate", "SwitchingDiffusionStabilityCertificate", "EulerRelaxationWindowDiagnostic", "TimeVaryingDiffusionStabilityBound", "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", "compute_diffusion_energy", "derive_time_varying_diffusion_stability_bound", "diagnose_euler_relaxation_window", "verify_heterogeneous_diffusion_stability", "verify_switching_diffusion_stability", "structural_eigenvalues", "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", ] # Two binary64 significands make the rational bisection finer than any # representable rate near a well-scaled endpoint. Stopping early can only make # the returned lower bound more conservative. _EXACT_QUOTIENT_BISECTION_STEPS = 2 * sys.float_info.mant_dig # Refinement is optional because the inverse-norm endpoint is already a proof. # Limiting exact LDL bisection to four quotient coordinates keeps larger graph # certificates practical while retaining tight bounds in small theorem tests. _EXACT_QUOTIENT_REFINEMENT_MAX_DIMENSION = 4 def _reject_boolean_numeric(value: Any, name: str) -> None: """Reject booleans before NumPy can coerce them to zero or one.""" if isinstance(value, (bool, np.bool_)): raise ValueError(f"{name} must be numeric, not boolean") def _readonly_float_array(value: Any) -> Any: """Return a detached binary64 array whose ordinary writes are disabled.""" result = np.array(value, dtype=float, copy=True) result.setflags(write=False) return result def _finite_float_signature(value: Any) -> tuple[str, ...]: """Encode one finite array exactly enough for an in-memory proof stamp.""" array = np.asarray(value, dtype=float) if not np.all(np.isfinite(array)): raise ValueError("proof data must remain finite") return tuple(float(item).hex() for item in array.flat) def _fixed_flow_proof_stamp( nodes: Any, metric_weights: Any, exact_gap: Fraction, certified_rate: float, is_certified: bool, preserves_consensus_subspace: bool, preserves_uniform_fixed_points: bool, preserves_weighted_mean: bool, ) -> tuple[Any, ...]: """Snapshot the fields on which hybrid fixed-flow composition relies.""" return ( "fixed_heterogeneous_diffusion_v1", structural_proof_signature(tuple(nodes)), np.asarray(metric_weights).shape, _finite_float_signature(metric_weights), Fraction(exact_gap), float(certified_rate).hex(), bool(is_certified), bool(preserves_consensus_subspace), bool(preserves_uniform_fixed_points), bool(preserves_weighted_mean), ) def _switching_flow_proof_stamp( nodes: Any, regime_count: int, reference_metric_weights: Any, normalized_metric_weights: Any, exact_gaps: Any, certified_rates: Any, certified_rate: float, shares_exact_common_metric: bool, supports_exact_theorem: bool, preserves_consensus_by_regime: Any, preserves_consensus_subspace: bool, preserves_uniform_fixed_points_by_regime: Any, preserves_uniform_fixed_points: bool, preserves_weighted_mean_by_regime: Any, preserves_weighted_mean: bool, ) -> tuple[Any, ...]: """Snapshot the fields on which hybrid switching composition relies.""" return ( "exact_common_metric_switching_diffusion_v1", structural_proof_signature(tuple(nodes)), int(regime_count), np.asarray(reference_metric_weights).shape, _finite_float_signature(reference_metric_weights), np.asarray(normalized_metric_weights).shape, _finite_float_signature(normalized_metric_weights), tuple(Fraction(value) for value in exact_gaps), np.asarray(certified_rates).shape, _finite_float_signature(certified_rates), float(certified_rate).hex(), bool(shares_exact_common_metric), bool(supports_exact_theorem), tuple(bool(value) for value in preserves_consensus_by_regime), bool(preserves_consensus_subspace), tuple(bool(value) for value in preserves_uniform_fixed_points_by_regime), bool(preserves_uniform_fixed_points), tuple(bool(value) for value in preserves_weighted_mean_by_regime), bool(preserves_weighted_mean), ) def _fraction_lower_float(value: Fraction) -> float: """Largest nearby binary64 value that is known not to exceed ``value``.""" if value < 0: raise ValueError("a lower-rounded fraction must be nonnegative") try: rounded = float(value) except OverflowError: return math.nextafter(float("inf"), 0.0) if math.isinf(rounded): return math.nextafter(rounded, 0.0) if Fraction.from_float(rounded) > value: rounded = math.nextafter(rounded, 0.0) return rounded def _fraction_upper_float(value: Fraction) -> float: """Smallest nearby binary64 value that is known not to be below ``value``.""" if value < 0: raise ValueError("an upper-rounded fraction must be nonnegative") try: rounded = float(value) except OverflowError: return float("inf") if math.isinf(rounded): return rounded if Fraction.from_float(rounded) < value: rounded = math.nextafter(rounded, float("inf")) return rounded def _fraction_sqrt_upper_float(value: Fraction) -> float: """Return a proved binary64 upper bound on ``sqrt(value)``. The search compares squared binary64 candidates with the rational input, so the result does not rely on a directed-rounding guarantee from libm. Positive binary64 values have the same order as their unsigned bit patterns, which makes a fixed 63-step binary search sufficient. """ if value < 0: raise ValueError("a square-root bound must be nonnegative") if value == 0: return 0.0 low_bits = 0 high_bits = 0x7FF0000000000000 # positive infinity while low_bits + 1 < high_bits: middle_bits = (low_bits + high_bits) // 2 candidate = struct.unpack(">d", middle_bits.to_bytes(8, byteorder="big"))[0] if math.isinf(candidate) or Fraction.from_float(candidate) ** 2 >= value: high_bits = middle_bits else: low_bits = middle_bits return struct.unpack(">d", high_bits.to_bytes(8, byteorder="big"))[0] def _ordered_nodes(G: Any) -> list: """Stable node ordering for the matrix representation.""" return list(G.nodes()) def _nodal_frequencies(G: Any, nodes: list | None = None) -> Any: """Read the actual capacity vector; zero capacity freezes a node.""" if nodes is None: nodes = _ordered_nodes(G) try: raw = [ get_attr( G.nodes[node], ALIAS_VF, 0.0, conv=lambda value: value, strict=True ) for node in nodes ] except (TypeError, ValueError, OverflowError) as exc: raise ValueError( "Structural frequency values must be finite real numbers" ) from exc if any(isinstance(value, (bool, np.bool_)) for value in raw): raise ValueError( "Structural frequency must be a finite real scalar, not boolean" ) try: frequency = np.array( [finite_real_scalar(value, "Structural frequency") for value in raw], dtype=float, ) except (TypeError, ValueError, OverflowError) as exc: raise ValueError( "Structural frequency values must be finite real numbers" ) from exc if np.any(frequency < 0.0): raise ValueError("Structural frequency must be finite and nonnegative") return frequency 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). Finite effective edge weights are normalized in scaled coordinates, so a raw row sum need not fit in a float. Raw degree-weighted quantities retain their separate representability requirements. Parameters ---------- G : TNFRGraph Returns ------- (nodes, L_rw) : tuple[list, np.ndarray] The node ordering and the N×N random-walk Laplacian. """ conductance = read_conductance(G) probability, scale, _ = conductance.normalization() laplacian = conductance.dense(-probability) diagonal = np.diag_indices(len(conductance.nodes)) laplacian[diagonal] += scale > 0.0 return conductance.nodes, laplacian 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}. For symmetric adjacency, L_sym shares the spectrum of the 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. Asymmetric adjacency is rejected: it cannot be passed to a symmetric eigensolver. Directed damping rates are available through :func:`relaxation_spectrum`. Scaled conductance avoids overflowing raw row sums; square-root ratios retain representable symmetric coefficients at very unequal strengths. 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. """ conductance = read_conductance(G, nodes, symmetric=True) normalized, positive = conductance.symmetric_normalized_weights() lap = conductance.dense(-normalized) lap[np.diag_indices(len(conductance.nodes))] += positive return conductance.nodes, lap def _finite_scalar_epi(value: Any) -> float: """Read raw real EPI or an exact uniform-real BEPI embedding.""" is_bepi_representation = isinstance(value, Mapping) or all( hasattr(value, attribute) for attribute in ("f_continuous", "a_discrete", "x_grid") ) if not is_bepi_representation: return finite_real_scalar(value, "EPI") try: scalar = real_scalar_epi(value) except (KeyError, TypeError, ValueError, OverflowError) as exc: raise ValueError( "EPI must be a finite real scalar or uniform real BEPI embedding" ) from exc if scalar is None: raise ValueError( "EPI must be a finite real scalar or uniform real BEPI embedding" ) return scalar def structural_field(G: Any, nodes: list | None = None) -> Any: r"""Return the strict scalar EPI field aligned with ``nodes``. Raw finite real values and exact uniform-real BEPI embeddings share the signed scalar channel. Nonuniform or complex BEPI elements are rejected because their magnitude projection is not an equivalent diffusion state. """ if nodes is None: nodes = _ordered_nodes(G) raw = [ get_attr(G.nodes[node], ALIAS_EPI, 0.0, conv=lambda value: value, strict=True) for node in nodes ] return np.array([_finite_scalar_epi(value) for value in raw], dtype=float) def structural_diffusivity(G: Any) -> float: r"""Mean structural frequency, a descriptive capacity summary. In ∂EPI/∂t = −νf·L_rw·EPI, νf plays the role of the diffusivity: the larger the structural frequency, the faster the form spreads. This mean is a diffusion coefficient only when all nodal frequencies are equal; heterogeneous transport uses diag(νf)·L_rw, never mean(νf)·L_rw. """ vf = _nodal_frequencies(G) return float(np.mean(vf)) if len(vf) else 0.0 def degree_weighted_total(G: Any) -> float: r"""The row-strength-weighted total Σ_i d_i·EPI(i). For symmetric adjacency and a common frequency this total is conserved. For fixed positive heterogeneous frequencies the conserved weights are d_i/νf_i instead; a directed graph requires the generator's left nullspace. This helper keeps its literal weighted-total meaning in every case. Nonfinite EPI or an unrepresentable final total raises ValueError. Edge-based summation allows finite cancellation even when an intermediate degree or product exceeds float range. Under the stated symmetric homogeneous restriction this is an exact **EPI-channel** conserved quantity. It is distinct from the Noether-like diagnostic ``Q = Σ(Φ_s + K_φ)`` (:func:`tnfr.physics.conservation.compute_noether_charge`). Grammar U1–U6 does not generally conserve that diagnostic; its drift must be evaluated on the supplied trajectory. """ conductance = read_conductance(G) field = structural_field(G, conductance.nodes) if not np.all(np.isfinite(field)): raise ValueError("Structural transport requires finite scalar EPI") total = conductance.weighted_total(field) if not np.isfinite(total): raise ValueError("Degree-weighted total exceeds finite floating-point range") return total def _read_edge_flux(G: Any) -> tuple[ConductanceSnapshot, Any, Any]: """Share edge differences between current, divergence and energy gradient. Capacity is deliberately absent: constitutive current can remain nonzero at a frozen node. Effective zero edges are removed before subtraction. """ conductance = read_conductance(G, symmetric=True) field = structural_field(G, conductance.nodes) if not np.all(np.isfinite(field)): raise ValueError("Structural transport requires finite scalar EPI") try: with np.errstate(over="raise", under="raise", invalid="raise"): difference = field[conductance.source] - field[conductance.target] flux = conductance.weight * difference except FloatingPointError as exc: raise ValueError( "Structural current is outside finite floating-point range" ) from exc return conductance, difference, flux @dataclass(frozen=True) class DiffusionEnergyBalance: """Instantaneous balance for the fixed symmetric EPI-only channel. Arrays follow ``nodes`` and are detached from the graph. ``gradient`` is the Euclidean derivative of ``energy``; ``epi_rate`` includes the actual nodal capacities. Zero mobility is permitted and does not imply zero pressure. This is not a certificate for the full tetrad energy, changing topology, or arbitrary finite integration steps. """ nodes: list energy: float gradient: Any mobility: Any epi_rate: Any energy_rate: float @dataclass(frozen=True) class HeterogeneousDiffusionStabilityCertificate: """Lyapunov certificate for connected symmetric EPI diffusion. ``metric_weights`` is the binary64 evaluation of ``d_i / nu_f_i`` and ``lyapunov_value`` is the corresponding distance to consensus. The eigensolver fields are estimates. ``exact_quotient_gap_lower_bound`` is obtained from the actual represented generator and metric by rational arithmetic; ``certified_exponential_rate_lower_bound`` is its downward-rounded energy rate. Hybrid composition uses only the latter. ``exact_consensus_subspace_preservation`` says the represented generator maps constant vectors back into that subspace. ``exact_uniform_fixed_point_preservation`` is the stronger canonical diffusion requirement that it annihilate them; rounded Laplacian rows need not satisfy this automatically. ``exact_weighted_mean_preservation`` separately records whether this represented metric's mean is conserved exactly. ``conserved_total`` and ``equilibrium_value`` retain their compatibility names but are snapshot values of ``h.T x`` and its projection center unless that flag is true. The certificate concerns the frozen EPI-only channel; it does not certify changing topology or an operator word. """ nodes: tuple[Any, ...] equilibrium_value: float metric_weights: Any conserved_total: float lyapunov_value: float lyapunov_derivative: float generalized_gap: float exponential_rate: float derivative_upper_bound: float balance_residual: float is_consensus: bool is_equilibrium: bool is_certified: bool exact_quotient_gap_lower_bound: Fraction = Fraction(0) certified_exponential_rate_lower_bound: float = 0.0 exact_consensus_subspace_preservation: bool = False exact_uniform_fixed_point_preservation: bool = False exact_weighted_mean_preservation: bool = False _proof_stamp: tuple[Any, ...] = field(default=(), repr=False, compare=False) def _proof_fields_are_intact(self) -> bool: """Detect ordinary construction, replacement, or payload mutation. This private consistency stamp is not an authentication boundary; deliberate reconstruction of private fields is outside its scope. """ try: expected = _fixed_flow_proof_stamp( self.nodes, self.metric_weights, self.exact_quotient_gap_lower_bound, self.certified_exponential_rate_lower_bound, self.is_certified, self.exact_consensus_subspace_preservation, self.exact_uniform_fixed_point_preservation, self.exact_weighted_mean_preservation, ) observed = object.__getattribute__(self, "_proof_stamp") except BaseException: return False return proof_stamps_are_identical(observed, expected) @dataclass(frozen=True) class SwitchingDiffusionStabilityCertificate: """Common-metric certificate for a finite family of graph topologies. Each regime may change conductance and structural frequency. The exact common-Lyapunov theorem requires identical normalized weights ``d_i/nu_f_i`` in every regime and requires every represented generator to annihilate constant fields exactly. Numerical agreement within a caller chosen tolerance is reported separately and is not promoted to that theorem. ``equilibrium_value``, ``lyapunov_value``, and ``derivative_upper_bound`` are binary64 diagnostics for the first snapshot in the displayed normalized metric. The theorem is anchored instead in ``reference_metric_weights`` and the downward-rounded certified rate; callers evaluating other snapshots must recenter them independently. """ nodes: tuple[Any, ...] regime_count: int reference_metric_weights: Any normalized_metric_weights: Any metric_mismatches: Any generalized_gaps: Any equilibrium_value: float lyapunov_value: float uniform_exponential_rate: float derivative_upper_bound: float common_metric_residual: float shares_exact_common_metric: bool common_metric_within_tolerance: bool numerical_hypotheses_pass: bool supports_exact_switching_theorem: bool tolerance: float scope: str exact_quotient_gap_lower_bounds: tuple[Fraction, ...] = () certified_exponential_rate_lower_bounds: Any = field( default_factory=lambda: _readonly_float_array(()) ) certified_uniform_exponential_rate_lower_bound: float = 0.0 exact_consensus_subspace_preservation_by_regime: tuple[bool, ...] = () exact_common_consensus_subspace_preservation: bool = False exact_uniform_fixed_point_preservation_by_regime: tuple[bool, ...] = () exact_common_uniform_fixed_point_preservation: bool = False exact_weighted_mean_preservation_by_regime: tuple[bool, ...] = () exact_common_weighted_mean_preservation: bool = False _proof_stamp: tuple[Any, ...] = field(default=(), repr=False, compare=False) @property def shares_common_metric(self) -> bool: """Compatibility alias for exact common-metric equality.""" return self.shares_exact_common_metric @property def is_certified(self) -> bool: """Compatibility alias for the exact switching-theorem flag.""" return self.supports_exact_switching_theorem def _proof_fields_are_intact(self) -> bool: """Detect ordinary construction, replacement, or payload mutation. This private consistency stamp is not an authentication boundary; deliberate reconstruction of private fields is outside its scope. """ try: expected = _switching_flow_proof_stamp( self.nodes, self.regime_count, self.reference_metric_weights, self.normalized_metric_weights, self.exact_quotient_gap_lower_bounds, self.certified_exponential_rate_lower_bounds, self.certified_uniform_exponential_rate_lower_bound, self.shares_exact_common_metric, self.supports_exact_switching_theorem, self.exact_consensus_subspace_preservation_by_regime, self.exact_common_consensus_subspace_preservation, self.exact_uniform_fixed_point_preservation_by_regime, self.exact_common_uniform_fixed_point_preservation, self.exact_weighted_mean_preservation_by_regime, self.exact_common_weighted_mean_preservation, ) observed = object.__getattribute__(self, "_proof_stamp") except BaseException: return False return proof_stamps_are_identical(observed, expected) @dataclass(frozen=True) class TimeVaryingDiffusionStabilityBound: """Conditional exact-real bound induced by represented conductances. The supplied frequency arrays are assumptions on every instant of the schedule, not observations made by this read-only function. Certification interprets each effective binary64 conductance materialized by the shared reader, and each frequency bound, as an exact real number; it forms degrees and the Laplacian in rational arithmetic and proves a positive quotient gap on the Euclidean disagreement subspace. ``materialized_binary64_uniform_fixed_point_preservation`` separately says whether the ordinary binary64 ``diag(strength)-adjacency`` construction annihilates the uniform field; it is diagnostic and does not define the exact-real theorem. ``combinatorial_gap``, ``minimum_mobility``, ``maximum_mobility``, ``dirichlet_energy``, and the fields prefixed by ``spectral_`` are compatibility estimates. ``energy_decay_rate``, ``energy_derivative_upper_bound``, and ``consensus_distance_bound`` use only rational proof quantities and directed binary64 enclosures. An infinite consensus-distance bound records safe abstention. The spectral estimates may be nonfinite when the ordinary diagnostic path exceeds binary64 range. ``exact_real_continuous_time_model_certified`` records the rational theorem; ``operational_binary64_rate_available`` and its compatibility alias ``is_certified`` additionally require a positive public binary64 rate. ``runtime_integration_certified`` remains false because no numerical integrator or future schedule is observed here. """ nodes: tuple[Any, ...] frequency_lower_bounds: Any frequency_upper_bounds: Any combinatorial_gap: float minimum_mobility: float maximum_mobility: float dirichlet_energy: float energy_decay_rate: float energy_derivative_upper_bound: float consensus_distance_bound: float is_certified: bool spectral_energy_decay_rate_estimate: float spectral_consensus_distance_estimate: float exact_combinatorial_gap_lower_bound: Fraction certified_combinatorial_gap_lower_bound: float exact_minimum_mobility_lower_bound: Fraction certified_minimum_mobility_lower_bound: float exact_maximum_mobility_upper_bound: Fraction certified_maximum_mobility_upper_bound: float exact_energy_decay_rate_lower_bound: Fraction exact_dirichlet_energy: Fraction dirichlet_energy_lower_bound: float dirichlet_energy_upper_bound: float mobility_lower_bounds: Any mobility_upper_bounds: Any certified_mobility_lower_bounds: Any certified_mobility_upper_bounds: Any exact_real_uniform_fixed_point_preservation: bool materialized_binary64_uniform_fixed_point_preservation: bool exact_real_continuous_time_model_certified: bool operational_binary64_rate_available: bool runtime_integration_certified: bool scope: str @dataclass(frozen=True) class EulerRelaxationWindowDiagnostic: """Modal explicit-Euler relaxation for one frozen symmetric network. ``modal_steps`` counts integration steps, whereas ``policy_window`` counts operator positions. Both are reported for comparison and are not treated as interchangeable semantics. """ dt: float target_fraction: float spectral_relative_tolerance: float spectral_zero_threshold: float decay_rates: Any modal_multipliers: Any slowest_decay_rate: float fastest_decay_rate: float euler_stability_limit: float maximum_modal_factor: float modal_steps: int | None policy_window: int is_euler_stable: bool scope: str def _connected_symmetric_transport( G: Any, nodes: list | None = None, ) -> tuple[ConductanceSnapshot, Any, Any]: """Return one validated positive, connected symmetric transport snapshot.""" conductance = read_conductance(G, nodes, symmetric=True) if len(conductance.nodes) < 2: raise ValueError("Stability certification requires at least two nodes") strength = conductance.strength if np.any(strength <= 0.0): raise ValueError("Stability certification requires positive row strength") adjacency = conductance.dense() reached = {0} frontier = [0] while frontier: source = frontier.pop() for target in np.flatnonzero(adjacency[source] > 0.0): target = int(target) if target != source and target not in reached: reached.add(target) frontier.append(target) if len(reached) != len(conductance.nodes): raise ValueError( "Stability certification requires connected positive conductance" ) return conductance, adjacency, strength def _exact_projected_quadratic( matrix: tuple[tuple[Fraction, ...], ...], metric_weights: tuple[Fraction, ...], ) -> tuple[tuple[Fraction, ...], ...]: """Restrict a quadratic form to ``h.T y=0`` in a rational basis.""" dimension = len(metric_weights) pivot = max(range(dimension), key=metric_weights.__getitem__) free = tuple(index for index in range(dimension) if index != pivot) # Column a is e_free[a] - (h_free[a]/h_pivot)e_pivot. Choosing the # largest metric entry keeps every ratio at most one. Expanding this # two-sparse basis avoids an O(n^4) generic matrix multiplication. ratios = tuple(metric_weights[index] / metric_weights[pivot] for index in free) return tuple( tuple( ( matrix[free[left]][free[right]] - ratios[left] * matrix[pivot][free[right]] - ratios[right] * matrix[free[left]][pivot] + ratios[left] * ratios[right] * matrix[pivot][pivot] ) for right in range(dimension - 1) ) for left in range(dimension - 1) ) def _exact_inverse_norm_gap_lower_bound( dissipation: tuple[tuple[Fraction, ...], ...], metric: tuple[tuple[Fraction, ...], ...], ) -> Fraction: r"""Bound ``min z'Kz/z'Gz`` with exact rational arithmetic. For symmetric positive-definite ``K`` and ``G``, ``lambda_min(K,G) >= 1/(||K^-1||_inf ||G||_inf)``. That norm bound initializes a safe lower endpoint. Coordinate Rayleigh quotients give an upper endpoint, and exact semidefiniteness tests then refine the lower endpoint by rational bisection. Every returned value is a proved lower bound; floating-point eigensolver output is never used. """ if not _exact_symmetric_semidefinite(dissipation, strict=True): return Fraction(0) inverse = _exact_matrix_inverse(dissipation) inverse_norm = max( sum((abs(value) for value in row), Fraction(0)) for row in inverse ) metric_norm = max(sum((abs(value) for value in row), Fraction(0)) for row in metric) if inverse_norm <= 0 or metric_norm <= 0: return Fraction(0) lower = Fraction(1) / (inverse_norm * metric_norm) upper = min( dissipation[index][index] / metric[index][index] for index in range(len(dissipation)) ) if upper <= lower: return lower def candidate_is_valid(candidate: Fraction) -> bool: shifted = tuple( tuple( dissipation[i][j] - candidate * metric[i][j] for j in range(len(dissipation)) ) for i in range(len(dissipation)) ) return _exact_symmetric_semidefinite(shifted, strict=False) if candidate_is_valid(upper): return upper if len(dissipation) > _EXACT_QUOTIENT_REFINEMENT_MAX_DIMENSION: return lower for _ in range(_EXACT_QUOTIENT_BISECTION_STEPS): midpoint = (lower + upper) / 2 if candidate_is_valid(midpoint): lower = midpoint else: upper = midpoint return lower @lru_cache(maxsize=128) def _exact_real_laplacian_gap_lower_bound( laplacian: tuple[tuple[Fraction, ...], ...], ) -> tuple[Fraction, bool]: """Certify the Euclidean disagreement gap of one rational Laplacian. The caller forms degrees exactly from the rational interpretations of the stored binary64 conductances. Symmetry and annihilation of the uniform field are still checked before the quotient proof, so a malformed internal matrix can only cause safe abstention. """ dimension = len(laplacian) is_symmetric = all( laplacian[i][j] == laplacian[j][i] for i in range(dimension) for j in range(i) ) preserves_uniform_fixed_point = is_symmetric and all( sum(row, Fraction(0)) == 0 for row in laplacian ) if not preserves_uniform_fixed_point: return Fraction(0), False unit_weights = (Fraction(1),) * dimension identity = tuple( tuple(Fraction(i == j) for j in range(dimension)) for i in range(dimension) ) restricted_laplacian = _exact_projected_quadratic(laplacian, unit_weights) restricted_identity = _exact_projected_quadratic(identity, unit_weights) gap = _exact_inverse_norm_gap_lower_bound(restricted_laplacian, restricted_identity) return gap, True def _exact_real_dirichlet_energy( laplacian: tuple[tuple[Fraction, ...], ...], field_values: Any, ) -> Fraction: """Evaluate ``x.T @ B @ x / 2`` for a rational model exactly.""" exact_field = tuple( Fraction.from_float(float(value)) for value in np.asarray(field_values, dtype=float) ) return ( sum( ( exact_field[i] * laplacian[i][j] * exact_field[j] for i in range(len(exact_field)) for j in range(len(exact_field)) ), Fraction(0), ) / 2 ) @lru_cache(maxsize=128) def _exact_flow_gap_from_rationals( laplacian: tuple[tuple[Fraction, ...], ...], mobility: tuple[Fraction, ...], metric: tuple[Fraction, ...], ) -> tuple[Fraction, bool, bool, bool]: """Cached exact quotient proof for one represented frozen generator.""" generator = tuple( tuple(mobility[i] * laplacian[i][j] for j in range(len(metric))) for i in range(len(metric)) ) mapped_consensus = tuple(sum(row, Fraction(0)) for row in generator) preserves_consensus_subspace = all( value == mapped_consensus[0] for value in mapped_consensus[1:] ) preserves_uniform_fixed_points = all(value == 0 for value in mapped_consensus) weighted_generator_row = tuple( sum( (metric[i] * generator[i][j] for i in range(len(metric))), Fraction(0), ) for j in range(len(metric)) ) preserves_weighted_mean = all(value == 0 for value in weighted_generator_row) if not preserves_consensus_subspace: return ( Fraction(0), preserves_weighted_mean, False, preserves_uniform_fixed_points, ) weighted_mobility = tuple( metric[index] * mobility[index] for index in range(len(metric)) ) symmetric_dissipation = tuple( tuple( ( weighted_mobility[i] * laplacian[i][j] + laplacian[i][j] * weighted_mobility[j] ) / 2 for j in range(len(metric)) ) for i in range(len(metric)) ) metric_matrix = tuple( tuple(metric[i] if i == j else Fraction(0) for j in range(len(metric))) for i in range(len(metric)) ) restricted_dissipation = _exact_projected_quadratic(symmetric_dissipation, metric) restricted_metric = _exact_projected_quadratic(metric_matrix, metric) gap = _exact_inverse_norm_gap_lower_bound(restricted_dissipation, restricted_metric) return ( gap, preserves_weighted_mean, True, preserves_uniform_fixed_points, ) def _exact_flow_gap_lower_bound( laplacian: Any, mobility: Any, displayed_metric_weights: Any, ) -> tuple[Fraction, bool, bool, bool]: r"""Certify contraction in the displayed binary64 metric. ``fl(d/nu)`` need not be the exact reciprocal of ``fl(nu/d)``. Exact rational arithmetic therefore starts from the binary64 Laplacian and mobility already materialized by the verifier, including all degree-sum, subtraction, and division rounding. It forms the quotient dissipation of that represented generator directly. Exact LDL verifies positivity, and an inverse-norm bound supplies a conservative generalized gap. Small quotients receive an optional exact bisection refinement. The accompanying Booleans separately record preservation of the consensus subspace, uniform fixed points, and the displayed metric's weighted mean. """ metric = tuple( Fraction.from_float(float(value)) for value in np.asarray(displayed_metric_weights, dtype=float) ) exact_laplacian = tuple( tuple(Fraction.from_float(float(value)) for value in row) for row in np.asarray(laplacian, dtype=float) ) exact_mobility = tuple( Fraction.from_float(float(value)) for value in np.asarray(mobility, dtype=float) ) return _exact_flow_gap_from_rationals(exact_laplacian, exact_mobility, metric) def _exact_represented_flow_derivative_is_zero( laplacian: Any, mobility: Any, field: Any, ) -> bool: """Test ``diag(mobility) @ laplacian @ field == 0`` exactly.""" exact_field = tuple( Fraction.from_float(float(value)) for value in np.asarray(field) ) exact_mobility = tuple( Fraction.from_float(float(value)) for value in np.asarray(mobility) ) exact_laplacian = tuple( tuple(Fraction.from_float(float(value)) for value in row) for row in np.asarray(laplacian) ) return all( exact_mobility[i] * sum( (exact_laplacian[i][j] * exact_field[j] for j in range(len(exact_field))), Fraction(0), ) == 0 for i in range(len(exact_field)) ) def _float_flow_quotient_gap( laplacian: Any, frequency: Any, strength: Any, metric_weights: Any, ) -> float: """Estimate the actual flow quotient gap in the displayed metric.""" dimension = len(metric_weights) pivot = int(np.argmax(metric_weights)) free = np.array([index for index in range(dimension) if index != pivot], dtype=int) basis = np.zeros((dimension, dimension - 1), dtype=float) basis[pivot, :] = -metric_weights[free] / metric_weights[pivot] basis[free, np.arange(dimension - 1)] = 1.0 # Match the represented theorem and the derivative path: materialize # mobility first, then multiply by the displayed metric. Reassociating # these binary64 operations can change the diagnostic by one ulp. weighted_mobility = metric_weights * (frequency / strength) symmetric_dissipation = 0.5 * ( weighted_mobility[:, None] * laplacian + laplacian * weighted_mobility[None, :] ) restricted_dissipation = basis.T @ symmetric_dissipation @ basis restricted_metric = basis.T @ (metric_weights[:, None] * basis) factor = np.linalg.cholesky(restricted_metric) left_solved = np.linalg.solve(factor, restricted_dissipation) normalized = np.linalg.solve(factor, left_solved.T).T normalized = 0.5 * (normalized + normalized.T) rates = np.linalg.eigvalsh(normalized) if not np.all(np.isfinite(rates)): raise ValueError("flow quotient rate estimate must be finite") return float(max(rates[0], 0.0)) def _aligned_frequency_bound(values: Any, nodes: list, name: str) -> Any: """Align one scalar, mapping or vector capacity bound with ``nodes``.""" _reject_boolean_numeric(values, name) if np.isscalar(values): result = np.full(len(nodes), float(values), dtype=float) elif hasattr(values, "keys"): try: aligned = [values[node] for node in nodes] except (KeyError, TypeError) as exc: raise ValueError(f"{name} must define every graph node") from exc if any(isinstance(value, (bool, np.bool_)) for value in aligned): raise ValueError(f"{name} must contain numeric values, not booleans") try: result = np.array(aligned, dtype=float) except (TypeError, ValueError) as exc: raise ValueError(f"{name} must define every graph node") from exc else: try: unconverted = np.asarray(values, dtype=object) except (TypeError, ValueError) as exc: raise ValueError(f"{name} must be scalar or node-aligned") from exc if any(isinstance(value, (bool, np.bool_)) for value in unconverted.flat): raise ValueError(f"{name} must contain numeric values, not booleans") try: result = np.asarray(values, dtype=float) except (TypeError, ValueError) as exc: raise ValueError(f"{name} must be scalar or node-aligned") from exc if result.shape != (len(nodes),): raise ValueError(f"{name} must contain one value per graph node") if not np.all(np.isfinite(result)) or np.any(result <= 0.0): raise ValueError(f"{name} must be finite and positive") return result def compute_diffusion_energy(G: Any) -> DiffusionEnergyBalance: r"""Read the Dirichlet gradient-flow balance on symmetric conductance. With B = D-W and x = EPI, E_D = x^T B x / 2 and grad(E_D) = Bx. The canonical pure EPI channel is x' = -M grad(E_D), where M_ii = nu_f_i/d_i for positive row strength and zero otherwise. Thus E_D' = -grad(E_D)^T M grad(E_D) <= 0. For positive capacity and strength this is the gradient flow in metric diag(d_i/nu_f_i); zero capacity instead gives a degenerate positive-semidefinite mobility. Nonnegative finite weights, finite scalar EPI and nonnegative finite capacities are required. Parallel edges and loops follow the shared diffusion adjacency convention; loops add strength but no energy. Asymmetric adjacency has no such symmetric Dirichlet identity and is rejected. No graph attributes, pressure callbacks or caches are changed. Overflow and inexact underflow in fluxes or balance terms raise ValueError; a nonzero contribution must not silently become a zero-energy diagnostic. A row strength itself may exceed float range when the returned mobility, gradient, energy and rates remain representable in scaled coordinates. """ try: with np.errstate(over="raise", under="raise", invalid="raise", divide="raise"): conductance, difference, flux = _read_edge_flux(G) nodes = conductance.nodes frequency = _nodal_frequencies(G, nodes) mobility = conductance.divide_by_strength(frequency) gradient = conductance.divergence(flux) energy = float(0.25 * np.sum(flux * difference)) epi_rate = -mobility * gradient # Use elementwise products so NumPy checks underflow as well as # overflow; a BLAS dot product may bypass these floating exceptions. energy_rate = float(np.sum(gradient * epi_rate)) except FloatingPointError as exc: raise ValueError( "Diffusion energy balance is outside finite floating-point range" ) from exc return DiffusionEnergyBalance( nodes, energy, gradient, mobility, epi_rate, energy_rate ) def verify_heterogeneous_diffusion_stability( G: Any, *, tolerance: float = 1e-10, ) -> HeterogeneousDiffusionStabilityCertificate: r"""Certify exponential convergence of the frozen heterogeneous EPI channel. Let ``W`` be a fixed symmetric nonnegative conductance on a connected graph. The verifier materializes binary64 ``B = fl(D-W)`` and ``M = diag(fl(nu_f_i/d_i))``; the represented channel is ``x'=-M B x``. All strengths and capacities must be positive. The published metric is ``H = diag(fl(d_i/nu_f_i))``; its rounding means that ``H M`` need not be exactly the identity, and rounded rows of ``B`` need not sum to exact zero. A canonical diffusion certificate therefore requires ``M B 1=0`` exactly. With ``c = (h^T x)/(h^T 1)`` and ``y = x-c*1``, the weighted energy ``V = 1/2 y^T H y`` has ``V' = -y^T S y``, where ``S = (H M B + B M H)/2``. The implementation restricts ``S`` and ``H`` to ``h^T y=0``, proves positivity using exact rational LDL elimination, and bounds the generalized quotient with exact matrix norms. Consequently the disagreement energy decays at the returned certified rate. Exact preservation of ``h^T x`` is reported separately; only then is the initial value of ``c`` itself the certified consensus limit. Without that extra condition the theorem establishes disagreement decay, while the scalar consensus coordinate requires its own conclusion. This is an analytic certificate for continuous-time pure EPI diffusion. It makes no claim about the other pressure channels, finite-step integrator stability, changing capacities/topology, or grammar U2 sufficiency. """ _reject_boolean_numeric(tolerance, "tolerance") if not np.isfinite(tolerance) or tolerance <= 0.0: raise ValueError("tolerance must be finite and positive") conductance, adjacency, strength = _connected_symmetric_transport(G) nodes = conductance.nodes field = structural_field(G, nodes) frequency = _nodal_frequencies(G, nodes) if not np.all(np.isfinite(field)): raise ValueError("Structural transport requires finite scalar EPI") if np.any(frequency <= 0.0): raise ValueError( "Stability certification requires positive structural frequency" ) try: with np.errstate(over="raise", invalid="raise", divide="raise", under="ignore"): metric_weights = strength / frequency if not np.all(np.isfinite(metric_weights)) or np.any(metric_weights <= 0.0): raise ValueError("Stability metric weights must be finite and positive") metric_total = float(np.sum(metric_weights)) conserved_total = float(metric_weights @ field) equilibrium = conserved_total / metric_total centered = field - equilibrium lyapunov_value = float(0.5 * np.sum(metric_weights * centered**2)) laplacian = np.diag(strength) - adjacency generalized_gap = _float_flow_quotient_gap( laplacian, frequency, strength, metric_weights ) exponential_rate = 2.0 * generalized_gap mobility = frequency / strength diffusion_rate = -mobility * (laplacian @ field) lyapunov_derivative = float( np.sum(metric_weights * centered * diffusion_rate) ) weighted_mobility = metric_weights * mobility symmetric_dissipation = 0.5 * ( weighted_mobility[:, None] * laplacian + laplacian * weighted_mobility[None, :] ) quadratic_derivative = -float(centered @ symmetric_dissipation @ centered) balance_residual = abs(lyapunov_derivative - quadratic_derivative) except ( FloatingPointError, OverflowError, np.linalg.LinAlgError, ValueError, ) as exc: raise ValueError( "Stability certificate exceeds finite floating-point range" ) from exc ( exact_quotient_gap, exact_mean_preservation, exact_consensus_preservation, exact_uniform_fixed_points, ) = _exact_flow_gap_lower_bound(laplacian, mobility, metric_weights) certified_rate = _fraction_lower_float(2 * exact_quotient_gap) try: derivative_upper_bound = -certified_rate * lyapunov_value except OverflowError as exc: raise ValueError( "Stability certificate exceeds finite floating-point range" ) from exc if not np.isfinite(derivative_upper_bound): raise ValueError("Stability certificate exceeds finite floating-point range") positive_gap = certified_rate > 0.0 is_consensus = bool(np.max(np.abs(centered)) <= tolerance) is_equilibrium = _exact_represented_flow_derivative_is_zero( laplacian, mobility, field ) # Canonical pure-EPI diffusion must leave every uniform field fixed. certified = bool(positive_gap and exact_uniform_fixed_points) node_tuple = tuple(nodes) frozen_metric = _readonly_float_array(metric_weights) proof_stamp = _fixed_flow_proof_stamp( node_tuple, frozen_metric, exact_quotient_gap, certified_rate, certified, exact_consensus_preservation, exact_uniform_fixed_points, exact_mean_preservation, ) return HeterogeneousDiffusionStabilityCertificate( nodes=node_tuple, equilibrium_value=equilibrium, metric_weights=frozen_metric, conserved_total=conserved_total, lyapunov_value=lyapunov_value, lyapunov_derivative=lyapunov_derivative, generalized_gap=generalized_gap, exponential_rate=exponential_rate, derivative_upper_bound=derivative_upper_bound, balance_residual=balance_residual, is_consensus=is_consensus, is_equilibrium=is_equilibrium, is_certified=certified, exact_quotient_gap_lower_bound=exact_quotient_gap, certified_exponential_rate_lower_bound=certified_rate, exact_consensus_subspace_preservation=exact_consensus_preservation, exact_uniform_fixed_point_preservation=exact_uniform_fixed_points, exact_weighted_mean_preservation=exact_mean_preservation, _proof_stamp=proof_stamp, ) def verify_switching_diffusion_stability( graphs: Any, *, tolerance: float = 1e-10, ) -> SwitchingDiffusionStabilityCertificate: r"""Certify pure-EPI convergence under arbitrary switching of topology. Consider a finite family of fixed, connected, symmetric conductance graphs on one node set. Regime ``r`` has materialized binary64 ``B_r``, mobility ``M_r=diag(fl(nu_r/d_r))``, and metric weights ``h_r=fl(d_r/nu_r)``. If every represented ``h_r`` is a positive scalar multiple of one vector ``p``, then ``V=1/2 sum_i p_i (x_i-c)^2`` is common to all regimes, provided every represented generator also obeys ``M_r B_r 1=0`` exactly. Exact rational quotient tests certify its decay under each represented generator, so arbitrary switching satisfies ``V(t)<=exp(-r_* t)V(0)`` at the returned certified rate ``r_*``. Preservation of ``p^T x`` is tested separately; when it holds in every regime, ``c=p^T x/p^T 1`` is also conserved under arbitrary switching. The theorem requires exact proportionality. The certificate additionally reports a within-tolerance diagnostic, but a large caller tolerance cannot turn a nonzero metric mismatch into an exact common Lyapunov function. Metric normalization is performed after max-scaling. A positive component that would underflow to zero is rejected because losing it can make two distinct metrics appear exactly equal. The certificate allows no node insertion/removal, directed edge, isolate, zero capacity, nonlinear pressure channel or reset at a switch. Failure of the common-metric condition is a negative result for this theorem, not evidence of instability. """ _reject_boolean_numeric(tolerance, "tolerance") if not np.isfinite(tolerance) or tolerance <= 0.0: raise ValueError("tolerance must be finite and positive") regimes = list(graphs) if len(regimes) < 2: raise ValueError("switching stability requires at least two regimes") nodes = list(regimes[0]) node_set = set(nodes) raw_metrics = [] normalized_metrics = [] gaps = [] exact_gap_bounds: list[Fraction] = [] exact_consensus_flags: list[bool] = [] exact_uniform_fixed_flags: list[bool] = [] exact_mean_flags: list[bool] = [] for graph in regimes: if len(graph) != len(nodes) or set(graph) != node_set: raise ValueError("switching regimes must share exactly one node set") conductance, adjacency, strength = _connected_symmetric_transport(graph, nodes) frequency = _nodal_frequencies(graph, nodes) if np.any(frequency <= 0.0): raise ValueError( "switching stability requires positive structural frequency" ) try: with np.errstate(over="raise", invalid="raise", divide="raise"): metric = strength / frequency if not np.all(np.isfinite(metric)) or np.any(metric <= 0.0): raise ValueError( "switching metric weights must be finite and positive" ) raw_metrics.append(metric.copy()) normalized_metric, _, _ = normalize_weights(metric) if not np.all(np.isfinite(normalized_metric)) or np.any( normalized_metric <= 0.0 ): raise ValueError( "switching metric normalization exceeds floating-point " "dynamic range" ) normalized_metrics.append(normalized_metric) laplacian = np.diag(strength) - adjacency mobility = frequency / strength gaps.append( _float_flow_quotient_gap(laplacian, frequency, strength, metric) ) ( exact_gap_bound, exact_mean, exact_consensus, exact_uniform_fixed, ) = _exact_flow_gap_lower_bound(laplacian, mobility, metric) exact_gap_bounds.append(exact_gap_bound) exact_mean_flags.append(exact_mean) exact_consensus_flags.append(exact_consensus) exact_uniform_fixed_flags.append(exact_uniform_fixed) except (FloatingPointError, OverflowError, np.linalg.LinAlgError) as exc: raise ValueError( "switching certificate exceeds finite floating-point range" ) from exc reference = normalized_metrics[0] mismatches = np.array( [np.max(np.abs(metric - reference)) for metric in normalized_metrics], dtype=float, ) common_metric_residual = float(np.max(mismatches)) within_tolerance = bool(common_metric_residual <= tolerance) # Normalizing before comparison can collapse two distinct binary64 ratios # to the same rounded vector. Cross-products of exact float ratios test # proportionality of the represented positive metric vectors directly. def exactly_proportional(left: Any, right: Any) -> bool: left_ratio = tuple(Fraction.from_float(float(value)) for value in left) right_ratio = tuple(Fraction.from_float(float(value)) for value in right) return all( lhs * right_ratio[0] == rhs * left_ratio[0] for lhs, rhs in zip(left_ratio, right_ratio) ) raw_reference = raw_metrics[0] exact_common = all( exactly_proportional(metric, raw_reference) for metric in raw_metrics[1:] ) field = structural_field(regimes[0], nodes) if not np.all(np.isfinite(field)): raise ValueError("Structural transport requires finite scalar EPI") try: with np.errstate(over="raise", invalid="raise", divide="raise"): field_scale = float(np.max(np.abs(field), initial=0.0)) equilibrium = ( 0.0 if field_scale == 0.0 else float(field_scale * (reference @ (field / field_scale))) ) centered = field - equilibrium value = float(0.5 * np.sum(reference * centered**2)) uniform_rate = 2.0 * min(gaps) exact_uniform_gap_bound = min(exact_gap_bounds) certified_uniform_rate = _fraction_lower_float(2 * exact_uniform_gap_bound) derivative_upper_bound = -certified_uniform_rate * value except (FloatingPointError, OverflowError) as exc: raise ValueError( "switching Lyapunov diagnostics exceed finite floating-point range" ) from exc if not all( np.isfinite(number) for number in (equilibrium, value, uniform_rate, derivative_upper_bound) ): raise ValueError( "switching Lyapunov diagnostics exceed finite floating-point range" ) certified_rates = _readonly_float_array( [_fraction_lower_float(2 * bound) for bound in exact_gap_bounds] ) numerical_hypotheses_pass = bool( within_tolerance and certified_uniform_rate > 0.0 and all(exact_uniform_fixed_flags) ) exact_theorem = bool( exact_common and certified_uniform_rate > 0.0 and all(exact_uniform_fixed_flags) ) exact_common_consensus = bool(exact_common and all(exact_consensus_flags)) exact_common_uniform_fixed = bool(exact_common and all(exact_uniform_fixed_flags)) exact_common_mean = bool(exact_common and all(exact_mean_flags)) node_tuple = tuple(nodes) frozen_reference = _readonly_float_array(raw_reference) frozen_normalized = _readonly_float_array(reference) frozen_mismatches = _readonly_float_array(mismatches) frozen_gaps = _readonly_float_array(gaps) proof_stamp = _switching_flow_proof_stamp( node_tuple, len(regimes), frozen_reference, frozen_normalized, tuple(exact_gap_bounds), certified_rates, certified_uniform_rate, exact_common, exact_theorem, tuple(exact_consensus_flags), exact_common_consensus, tuple(exact_uniform_fixed_flags), exact_common_uniform_fixed, tuple(exact_mean_flags), exact_common_mean, ) return SwitchingDiffusionStabilityCertificate( nodes=node_tuple, regime_count=len(regimes), reference_metric_weights=frozen_reference, normalized_metric_weights=frozen_normalized, metric_mismatches=frozen_mismatches, generalized_gaps=frozen_gaps, equilibrium_value=equilibrium, lyapunov_value=value, uniform_exponential_rate=uniform_rate, derivative_upper_bound=derivative_upper_bound, common_metric_residual=common_metric_residual, shares_exact_common_metric=exact_common, common_metric_within_tolerance=within_tolerance, numerical_hypotheses_pass=numerical_hypotheses_pass, supports_exact_switching_theorem=exact_theorem, tolerance=float(tolerance), scope="finite fixed-node connected symmetric pure-EPI switching family", exact_quotient_gap_lower_bounds=tuple(exact_gap_bounds), certified_exponential_rate_lower_bounds=certified_rates, certified_uniform_exponential_rate_lower_bound=certified_uniform_rate, exact_consensus_subspace_preservation_by_regime=tuple(exact_consensus_flags), exact_common_consensus_subspace_preservation=exact_common_consensus, exact_uniform_fixed_point_preservation_by_regime=tuple( exact_uniform_fixed_flags ), exact_common_uniform_fixed_point_preservation=(exact_common_uniform_fixed), exact_weighted_mean_preservation_by_regime=tuple(exact_mean_flags), exact_common_weighted_mean_preservation=exact_common_mean, _proof_stamp=proof_stamp, ) def derive_time_varying_diffusion_stability_bound( G: Any, frequency_lower_bounds: Any, frequency_upper_bounds: Any, ) -> TimeVaryingDiffusionStabilityBound: r"""Derive a common-Lyapunov bound for time-varying nodal capacities. Each effective binary64 conductance materialized by the shared reader and each supplied binary64 capacity bound is interpreted as an exact real coefficient. Exact rational row sums define ``D`` and ``B=D-W``. Assume a continuous-time capacity schedule satisfying ``lower_i <= nu_i(t) <= upper_i`` at every instant. For ``E_D=x^T Bx/2`` the induced pure-EPI model obeys ``E_D' = -(Bx)^T diag(nu_i(t)/d_i) (Bx)``. If ``mu = min_i lower_i/d_i`` and a rational quotient proof establishes ``lambda_2(B)>0`` on ``1_perp``, then ``E_D' <= -2*mu*lambda_2(B)*E_D``. The public operational decay rate is a downward-rounded lower bound on that exact rational rate. The derivative and current consensus-distance fields use corresponding directed enclosures. The ordinary binary64 eigengap, energy, and products remain explicitly labelled estimates for compatibility. Finite upper capacity bounds make the state velocity integrable in this exact-real continuous model, so the field converges to a uniform value. That value is generally schedule-dependent; no fixed weighted mean is claimed unless all capacities share a common scalar modulation. The routine cannot observe a future schedule and does not certify a numerical integration path. Its result is therefore a conditional analytic theorem induced by the represented input coefficients. """ conductance, adjacency, strength = _connected_symmetric_transport(G) nodes = conductance.nodes lower = _aligned_frequency_bound( frequency_lower_bounds, nodes, "frequency_lower_bounds" ) upper = _aligned_frequency_bound( frequency_upper_bounds, nodes, "frequency_upper_bounds" ) if np.any(lower > upper): raise ValueError("frequency lower bounds cannot exceed upper bounds") field = structural_field(G, nodes) if not np.all(np.isfinite(field)): raise ValueError("Structural transport requires finite scalar EPI") # These are compatibility diagnostics for the ordinary binary64 matrix # path. They never establish the theorem Boolean, so loss of range in an # estimate must not suppress a valid exact-real certificate. laplacian = np.diag(strength) - adjacency try: with np.errstate( over="ignore", invalid="ignore", divide="ignore", under="ignore" ): eigenvalues = np.linalg.eigvalsh(laplacian) combinatorial_gap = float(max(eigenvalues[1], 0.0)) except np.linalg.LinAlgError: combinatorial_gap = float("nan") with np.errstate(over="ignore", invalid="ignore", divide="ignore", under="ignore"): mobility_lower_estimates = lower / strength mobility_upper_estimates = upper / strength minimum_mobility = float(np.min(mobility_lower_estimates)) maximum_mobility = float(np.max(mobility_upper_estimates)) energy = float(0.5 * field @ laplacian @ field) spectral_decay_rate = float( np.multiply(np.multiply(2.0, minimum_mobility), combinatorial_gap) ) spectral_distance_estimate = ( float( np.sqrt( np.maximum( np.divide(np.multiply(2.0, energy), combinatorial_gap), 0.0, ) ) ) if combinatorial_gap > 0.0 else float("inf") ) exact_adjacency = tuple( tuple(Fraction.from_float(float(value)) for value in row) for row in np.asarray(adjacency, dtype=float) ) exact_strength = tuple(sum(row, Fraction(0)) for row in exact_adjacency) exact_laplacian = tuple( tuple( (exact_strength[i] if i == j else Fraction(0)) - exact_adjacency[i][j] for j in range(len(nodes)) ) for i in range(len(nodes)) ) exact_gap, exact_real_uniform_fixed = _exact_real_laplacian_gap_lower_bound( exact_laplacian ) exact_materialized_laplacian = tuple( tuple(Fraction.from_float(float(value)) for value in row) for row in np.asarray(laplacian, dtype=float) ) materialized_binary64_uniform_fixed = all( sum(row, Fraction(0)) == 0 for row in exact_materialized_laplacian ) exact_lower_frequency = tuple(Fraction.from_float(float(value)) for value in lower) exact_upper_frequency = tuple(Fraction.from_float(float(value)) for value in upper) exact_lower_mobility = tuple( exact_lower_frequency[i] / exact_strength[i] for i in range(len(nodes)) ) exact_upper_mobility = tuple( exact_upper_frequency[i] / exact_strength[i] for i in range(len(nodes)) ) exact_minimum_mobility = min(exact_lower_mobility) exact_maximum_mobility = max(exact_upper_mobility) certified_mobility_lower = _readonly_float_array( [_fraction_lower_float(value) for value in exact_lower_mobility] ) certified_mobility_upper = _readonly_float_array( [_fraction_upper_float(value) for value in exact_upper_mobility] ) certified_gap = _fraction_lower_float(exact_gap) certified_minimum_mobility = _fraction_lower_float(exact_minimum_mobility) certified_maximum_mobility = _fraction_upper_float(exact_maximum_mobility) exact_decay_rate = 2 * exact_minimum_mobility * exact_gap decay_rate = _fraction_lower_float(exact_decay_rate) exact_energy = _exact_real_dirichlet_energy(exact_laplacian, field) exact_real_model_certified = bool( exact_real_uniform_fixed and exact_gap > 0 and exact_minimum_mobility > 0 and exact_maximum_mobility > 0 and exact_energy >= 0 ) operational_rate_available = bool(exact_real_model_certified and decay_rate > 0.0) if not operational_rate_available: decay_rate = 0.0 if exact_energy >= 0: energy_lower_bound = _fraction_lower_float(exact_energy) energy_upper_bound = _fraction_upper_float(exact_energy) else: energy_lower_bound = float("-inf") energy_upper_bound = float("inf") if operational_rate_available: exact_derivative_magnitude = Fraction.from_float(decay_rate) * exact_energy derivative_bound = -_fraction_lower_float(exact_derivative_magnitude) else: derivative_bound = 0.0 if exact_gap > 0 and exact_energy >= 0: distance_bound = _fraction_sqrt_upper_float(2 * exact_energy / exact_gap) else: distance_bound = float("inf") return TimeVaryingDiffusionStabilityBound( nodes=tuple(nodes), frequency_lower_bounds=_readonly_float_array(lower), frequency_upper_bounds=_readonly_float_array(upper), combinatorial_gap=combinatorial_gap, minimum_mobility=minimum_mobility, maximum_mobility=maximum_mobility, dirichlet_energy=energy, energy_decay_rate=decay_rate, energy_derivative_upper_bound=derivative_bound, consensus_distance_bound=distance_bound, is_certified=operational_rate_available, spectral_energy_decay_rate_estimate=spectral_decay_rate, spectral_consensus_distance_estimate=spectral_distance_estimate, exact_combinatorial_gap_lower_bound=exact_gap, certified_combinatorial_gap_lower_bound=certified_gap, exact_minimum_mobility_lower_bound=exact_minimum_mobility, certified_minimum_mobility_lower_bound=(certified_minimum_mobility), exact_maximum_mobility_upper_bound=exact_maximum_mobility, certified_maximum_mobility_upper_bound=(certified_maximum_mobility), exact_energy_decay_rate_lower_bound=exact_decay_rate, exact_dirichlet_energy=exact_energy, dirichlet_energy_lower_bound=energy_lower_bound, dirichlet_energy_upper_bound=energy_upper_bound, mobility_lower_bounds=_readonly_float_array(mobility_lower_estimates), mobility_upper_bounds=_readonly_float_array(mobility_upper_estimates), certified_mobility_lower_bounds=certified_mobility_lower, certified_mobility_upper_bounds=certified_mobility_upper, exact_real_uniform_fixed_point_preservation=(exact_real_uniform_fixed), materialized_binary64_uniform_fixed_point_preservation=( materialized_binary64_uniform_fixed ), exact_real_continuous_time_model_certified=(exact_real_model_certified), operational_binary64_rate_available=operational_rate_available, runtime_integration_certified=False, scope=( "conditional exact-real continuous pure-EPI schedule induced by " "effective binary64 conductances materialized by the shared reader " "and binary64 capacity bounds" ), ) def diagnose_euler_relaxation_window( G: Any, *, dt: float, target_fraction: float | None = None, tolerance: float = 1e-12, ) -> EulerRelaxationWindowDiagnostic: r"""Compute the graph-specific modal relaxation steps for explicit Euler. For frozen connected symmetric pure EPI diffusion, continuous decay rates are the eigenvalues of ``diag(nu_f)L_rw``. Explicit Euler maps each nonstationary mode by ``q_k=1-dt*lambda_k``. The worst modal amplitude falls below ``target_fraction`` after the returned number of steps exactly when ``max_k |q_k| < 1``. Otherwise ``modal_steps`` is ``None``. The canonical U4 value is included as ``policy_window`` only to expose the calibration gap. It counts operator positions, not solver timesteps, and this diagnostic does not modify or validate grammar U4. ``tolerance`` is a dimensionless relative cutoff against the fastest decay rate; it is not an absolute rate in ``Hz_str``. Symmetry concerns effective conductance, not the graph container type. Reciprocal directed representations are admitted under the same positive, connected transport check; asymmetric effective weights are rejected. """ import math _reject_boolean_numeric(dt, "dt") if target_fraction is not None: _reject_boolean_numeric(target_fraction, "target_fraction") _reject_boolean_numeric(tolerance, "tolerance") if not np.isfinite(dt) or dt <= 0.0: raise ValueError("dt must be finite and positive") if target_fraction is None: target_fraction = 1.0 / (math.pi + 1.0) if not np.isfinite(target_fraction) or not 0.0 < target_fraction < 1.0: raise ValueError("target_fraction must lie strictly between zero and one") if not np.isfinite(tolerance) or not 0.0 < tolerance < 1.0: raise ValueError( "tolerance must be finite and lie strictly between zero and one" ) _connected_symmetric_transport(G) frequency = _nodal_frequencies(G) if np.any(frequency <= 0.0): raise ValueError( "Euler modal diagnostic requires positive structural frequency" ) rates = relaxation_spectrum(G) fastest_rate = float(np.max(rates, initial=0.0)) if not np.isfinite(fastest_rate) or fastest_rate <= 0.0: raise ValueError("Euler modal diagnostic requires a positive decay scale") zero_threshold = tolerance * fastest_rate positive = rates[rates > zero_threshold] zero_count = len(rates) - len(positive) if zero_count != 1 or not len(positive): raise ValueError("Euler modal diagnostic requires exactly one stationary mode") multipliers = 1.0 - dt * positive maximum_factor = float(np.max(np.abs(multipliers))) fastest = float(positive[-1]) stability_limit = 2.0 / fastest stable = maximum_factor < 1.0 steps = None if stable: if maximum_factor == 0.0: steps = 1 else: steps = max( 1, int(math.floor(math.log(target_fraction) / math.log(maximum_factor))) + 1, ) while maximum_factor**steps >= target_fraction: steps += 1 from ..config.physics_derivation import derive_bifurcation_window_from_physics return EulerRelaxationWindowDiagnostic( dt=float(dt), target_fraction=float(target_fraction), spectral_relative_tolerance=float(tolerance), spectral_zero_threshold=float(zero_threshold), decay_rates=positive.copy(), modal_multipliers=multipliers.copy(), slowest_decay_rate=float(positive[0]), fastest_decay_rate=fastest, euler_stability_limit=stability_limit, maximum_modal_factor=maximum_factor, modal_steps=steps, policy_window=derive_bifurcation_window_from_physics(), is_euler_stable=stable, scope="frozen connected symmetric pure EPI explicit-Euler modes", ) def relaxation_spectrum(G: Any) -> Any: r"""Real decay rates of diag(νf)·L_rw, sorted ascending. Equal frequencies reduce to νf·λ_k and reuse the cached symmetric geometry on an undirected graph. With heterogeneous frequency, the generator is diag(νf)·L_rw. Directed rates are real parts of its possibly complex eigenvalues; they do not specify oscillation or transient amplification. Returns ------- np.ndarray The decay rates sorted ascending (real parts). """ frequency = _nodal_frequencies(G) if not len(frequency): return np.empty(0, dtype=float) if not G.is_directed(): if np.all(frequency == frequency[0]): eig = _cached_eigenvalues(G) return frequency[0] * eig _, lap_sym = symmetric_normalized_laplacian(G) root_frequency = np.sqrt(frequency) # AB and BA share eigenvalues even if a frequency is zero. rates = np.linalg.eigvalsh( root_frequency[:, None] * lap_sym * root_frequency[None, :] ) else: _, lap = structural_diffusion_operator(G) rates = np.linalg.eigvals(frequency[:, None] * lap).real return np.maximum(np.sort(rates), 0.0) 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. """ _reject_boolean_numeric(decimals, "decimals") _, 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, a capacity summary (a diffusion coefficient only if uniform). spectral_gap : float Second-smallest real part of the geometry spectrum of L_rw. slowest_relaxation_rate : float Smallest positive real decay rate of diag(νf)·L_rw. 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). invariant_weighted_conserved : bool or None All left-nullspace coordinates of diag(νf)·L_rw remain constant. max_invariant_conservation_drift : float Maximum absolute drift in the orthonormal invariant coordinates. reaches_stationary_state : bool or None Final nodal derivative is below the sampled-flow threshold. final_stationarity_residual : float Maximum absolute final nodal derivative; distinct from ΔNFR pressure. """ 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 invariant_weighted_conserved: bool | None = None max_invariant_conservation_drift: float = 0.0 reaches_stationary_state: bool | None = None final_stationarity_residual: float = 0.0 @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 if self.invariant_weighted_conserved is None else self.invariant_weighted_conserved ) and ( self.relaxes_to_uniform if self.reaches_stationary_state is None else self.reaches_stationary_state ) ) 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"mean capacity νf={self.diffusivity:.4f}, " f"spectral gap λ₂={self.spectral_gap:.4f}, " f"slowest positive nodal rate={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}), " f"nodal invariants conserved={self.invariant_weighted_conserved}, " f"stationary={self.reaches_stationary_state}" ) 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. Checks the EPI pressure identity and samples the actual frozen-frequency nodal flow ``e' = -diag(νf) L_rw e``. Conservation is checked against all left-nullspace invariants of that generator. Degree-weighted conservation and global uniformity are reported separately: they need not hold with heterogeneous frequencies, disconnected components or frozen nodes. ``reaches_stationary_state`` tests the final nodal derivative, not pressure equilibrium: zero capacity can freeze a nonuniform field. This is a finite sampled-flow diagnostic, not a proof of eventual convergence for every graph. The explicit step must preserve the nonnegative diffusion transition. 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 """ _reject_boolean_numeric(dt, "dt") _reject_boolean_numeric(steps, "steps") _reject_boolean_numeric(tolerance, "tolerance") if not np.isfinite(dt) or dt <= 0.0: raise ValueError("dt must be finite and positive") if not isinstance(steps, (int, np.integer)) or steps < 0: raise ValueError("steps must be a nonnegative integer") if not np.isfinite(tolerance) or tolerance <= 0.0: raise ValueError("tolerance must be finite and positive") nodes, lap = structural_diffusion_operator(G) n = len(nodes) epi = structural_field(G, nodes) if not np.all(np.isfinite(epi)): raise ValueError("The diffusion certificate requires a finite EPI field") frequency = _nodal_frequencies(G, nodes) generator = frequency[:, None] * lap if dt * float(np.max(np.diag(generator), initial=0.0)) > 1.0: raise ValueError("dt exceeds the nonnegative explicit diffusion step bound") # (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), initial=0.0)) is_laplacian = residual < max(tolerance, 1e-12) except Exception: residual = float("nan") is_laplacian = False # diffusivity and spectrum diffusivity = float(np.mean(frequency)) if n else 0.0 eig = np.linalg.eigvals(lap).real eig.sort() spectral_gap = max(0.0, float(eig[1])) if n > 1 else 0.0 rates = relaxation_spectrum(G) slowest_rate = next((float(rate) for rate in rates if rate > tolerance), 0.0) # degree vector for the conserved weighted total deg = read_conductance(G, nodes).strength # Columns of U beyond rank(A) span ker(A.T), hence their scalar products # with e are conserved. On an undirected positive-frequency graph this # includes the analytic weights d_i/νf_i, component by component. left_vectors, singular_values, _ = np.linalg.svd(generator, full_matrices=True) rank_tolerance = ( n * np.finfo(float).eps * float(np.max(singular_values, initial=0.0)) ) rank = int(np.count_nonzero(singular_values > rank_tolerance)) invariants = left_vectors[:, rank:].T # (2)+(3) integrate the nodal flow at the specified physical time step. e = epi.copy() w0 = float(deg @ e) max_drift = 0.0 initial_invariants = invariants @ e invariant_drift = 0.0 for _ in range(steps): e = e - dt * (generator @ e) max_drift = max(max_drift, abs(float(deg @ e) - w0)) invariant_drift = max( invariant_drift, float( np.max( np.abs(invariants @ e - initial_invariants), initial=0.0, ) ), ) conserved = max_drift < max(tolerance, 1e-9 * (abs(w0) + 1e-12)) invariant_conserved = invariant_drift <= tolerance * max( 1.0, float(np.max(np.abs(initial_invariants), initial=0.0)), ) final_std = float(np.std(e)) if n else 0.0 initial_std = float(np.std(epi)) if n else 0.0 relaxes = final_std < max(1e-3, 1e-2 * (initial_std + 1e-12)) stationarity_residual = float(np.max(np.abs(generator @ e), initial=0.0)) stationary = stationarity_residual < max( 1e-3, 1e-2 * float( np.max(np.abs(generator @ epi), initial=0.0), ), ) 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, invariant_weighted_conserved=bool(invariant_conserved), max_invariant_conservation_drift=invariant_drift, reaches_stationary_state=bool(stationary), final_stationarity_residual=stationarity_residual, ) # --------------------------------------------------------------------------- # 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. This held-force diagnostic does not establish a projection from that Hamiltonian. 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 for name, value in ( ("nu_f", nu_f), ("pressure", pressure), ("dt", dt), ("steps", steps), ("tolerance", tolerance), ): _reject_boolean_numeric(value, name) # 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 limit of a specified damped graph wave # to the dissipative structural diffusion # --------------------------------------------------------------------------- def _damped_wave_roots(lambdas: Any, gamma: float) -> tuple[Any, Any]: """Return both complex roots without subtracting near-equal real terms. Form the larger-magnitude root directly and recover the other from their product, lambda. Zero modes include the gamma=0 double root. The spectral inputs and finite arithmetic range retain their callers' scope; this helper does not assert a nonlinear or runtime wave law. """ values = np.asarray(lambdas, dtype=float) root = np.sqrt(gamma * gamma - 4.0 * values + 0j) if gamma >= 0.0: fast = -gamma / 2.0 - root / 2.0 slow = np.divide(values, fast, out=np.zeros_like(fast), where=fast != 0.0) else: # Preserve the historical plus/minus ordering for negative damping. slow = -gamma / 2.0 + root / 2.0 fast = np.divide(values, slow, out=np.zeros_like(slow), where=slow != 0.0) return slow, fast 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 specified graph-wave model has the 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. """ _reject_boolean_numeric(gamma, "gamma") _, lap = structural_diffusion_operator(G) lambdas = np.sort(np.linalg.eigvals(lap).real) lambdas = np.clip(lambdas, 0.0, None) s_slow, s_fast = _damped_wave_roots(lambdas, gamma) return lambdas, s_slow.real, s_fast.real @dataclass(frozen=True) class OverdampedProjectionCertificate: r"""Verification of the specified graph-wave-to-diffusion limit. The model starts from the graph wave q̈ = −L q. 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). This does not establish a projection of the different isotropic Hamiltonian implemented in :mod:`tnfr.physics.symplectic_substrate`. 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 Euclidean error in nodal EPI coordinates between the damped-wave trajectory and exp(−L_rw t/γ)·q₀ over the sampled 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 specified damped graph wave approaches 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 graph 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 the overdamped diffusion limit of the specified graph 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. Initial EPI is transformed by sqrt(row strength) before modal projection, and both trajectories are decoded back to nodal coordinates before the relative-error comparison. A common scale avoids raw degree overflow; isolated coordinates use the identity transform and remain stationary from rest. This is a binary64 diagnostic of the specified held graph wave, not a runtime certificate. Unrepresentable transformed initial values or decoded endpoints raise ``ValueError`` instead of producing a trajectory certificate. 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 """ _reject_boolean_numeric(gamma, "gamma") _reject_boolean_numeric(n_time_samples, "n_time_samples") _reject_boolean_numeric(tolerance, "tolerance") 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 s_slow = _damped_wave_roots(lambdas, gamma)[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 = _damped_wave_roots(lam2, gamma)[0].real slow_gap = float(-s_gap) else: slow_gap = 0.0 diff_gap = lam2 / gamma # (ii) L_rw trajectories through the degree-coordinate similarity to L_sym. 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) conductance = read_conductance(G, sym_nodes, symmetric=True) _, row_scale, row_total = conductance.normalization() positive = row_scale > 0.0 degree_root = np.ones(n, dtype=float) if np.any(positive): roots = np.sqrt(row_scale[positive]) * np.sqrt(row_total[positive]) degree_root[positive] = roots / np.max(roots) scaled_q0 = degree_root * q0 if np.any((q0 != 0.0) & (scaled_q0 == 0.0)): raise ValueError("Wave degree coordinates are below floating-point range") c0 = V.T @ scaled_q0 ss, sf = _damped_wave_roots(w, gamma) 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: try: with np.errstate(over="raise", invalid="raise", divide="raise"): q_wave = ( V @ (a * np.exp(ss * t) + b * np.exp(sf * t)).real ) / degree_root q_diff = (V @ (c0 * np.exp(-w * t / gamma))) / degree_root except FloatingPointError as exc: raise ValueError( "Decoded wave coordinates exceed finite floating-point range" ) from exc 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 graph wave is the undamped standing-wave spectrum √λ_k (the discrete modes). The overdamped projection (γ→∞, :func:`verify_overdamped_projection`) collapses the damped graph 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 graph 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 """ _reject_boolean_numeric(gamma, "gamma") _reject_boolean_numeric(tolerance, "tolerance") 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) s_plus = _damped_wave_roots(lambdas, gamma)[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]: """Compatibility wrapper for the canonical normalized Laplacian. Reuse its zero-row convention for isolated nodes so the cached spectrum and the nodal diffusion operator have the same stationary modes. """ return symmetric_normalized_laplacian(G) 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. """ conductance = read_conductance(G) edges = tuple( sorted( zip( conductance.source.tolist(), conductance.target.tolist(), conductance.weight.tolist(), ) ) ) return (G.is_directed(), G.is_multigraph(), tuple(conductance.nodes), edges) def _cached_eigenvalues(G: Any) -> Any: """Cache only the spectrum when a read-out does not require mode shapes. A full decomposition, when already cached, supplies the same values. Otherwise eigvalsh avoids computing and retaining an unused N-by-N basis. The full-spectrum solver and its Laplacian are still dense. """ sig = _topology_signature(G) cache = G.graph.get("_tnfr_spectrum_cache") if isinstance(cache, dict) and cache.get("sig") == sig: return cache["vals"] _, lap = _symmetric_normalized_laplacian(G) values = np.clip(np.linalg.eigvalsh(lap), 0.0, None) G.graph["_tnfr_spectrum_cache"] = {"sig": sig, "vals": values} return values 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 cache -- L_sym and L_rw have the same spectrum. A value-only cache is upgraded on the first request for actual eigenvectors. """ sig = _topology_signature(G) cache = G.graph.get("_tnfr_spectrum_cache") if isinstance(cache, dict) and cache.get("sig") == sig and "vecs" in cache: 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_eigenvalues(G: Any) -> Any: """Return only the canonical structural eigenvalues. Unlike :func:`structural_eigenmodes`, this value-only readout neither computes nor retains eigenvectors unless a previous mode request already populated them. Asymmetric adjacency is rejected by the shared canonical Laplacian constructor. """ return _cached_eigenvalues(G).copy() 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. For symmetric adjacency the zero mode is proportional to sqrt(row strength), corresponding to uniform EPI in degree coordinates. Asymmetric adjacency is rejected. Returned arrays are independent copies. 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. """ eigenvalues, eigenvectors = _cached_eigh(G) return eigenvalues.copy(), eigenvectors.copy() def nodal_domain_count(mode: Any) -> int: r"""Legacy ordered sign-change statistic for a supplied mode vector. Count sign changes after discarding entries with absolute value at most ``1e-12``. The result depends on the supplied node ordering. No adjacency is read, so it does not count connected sign domains on a general graph or certify Courant ordering. The historical function name is retained. 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"""Read the auxiliary graph-wave spectrum without time integration. The separately specified graph-wave model assigns standing-wave frequencies :math:`\omega_k = \sqrt{\lambda_k}` to structural modes (:func:`verify_undamped_limit` / :func:`structural_eigenmodes`). Their positive differences give possible beat frequencies in that model's time convention. The first-order nodal equation alone does not derive this second-order wave, its physical time scale, or an observed vibration. Neither zero pressure nor spectral multiplicity proves a beat event or fractal structure. The legacy ``vibration_energy`` key stores half the spectral trace. An actual wave state's energy also depends on its modal amplitudes and velocities: in normalized orthonormal modal coordinates it is :math:`\frac12\sum_k(\dot q_k^2+\lambda_k q_k^2)`. The returned statistic equals that energy for the selected initialization ``q_k=1, qdot_k=0``; this function does not extract or assert that initialization from ``G``. Parameters ---------- G : TNFRGraph n_modes : int Nonnegative number of leading resonant frequencies to report. Zero leaves the spectrum list empty; the other statistics still use all available modes, including the returned total ``n_modes`` count. Returns ------- dict ``resonant_spectrum`` (leading :math:`\omega_k = \sqrt{\lambda_k}`), ``fundamental`` (the smallest retained frequency with eigenvalue strictly above the selected cutoff ``1e-9``), ``dominant_beat`` (the smallest adjacent retained frequency gap strictly above ``1e-9``), ``spectral_multiplicity`` (largest multiplicity after the implemented eigenvalue rounding), ``vibration_energy`` (legacy key for :math:`\tfrac12\sum\lambda_k`, not measured state energy), ``n_modes`` (count above the eigenvalue cutoff). The eigenvalue and frequency-gap cuts have different units in the supplied wave model; they are numerical read-out policies. A zero reported fundamental or beat can mean that positive values were excluded, not that all physical motion is absent. """ from numbers import Integral _reject_boolean_numeric(n_modes, "n_modes") if not isinstance(n_modes, Integral) or n_modes < 0: raise ValueError("n_modes must be a nonnegative integer") eigvals = _cached_eigenvalues(G) eigvals = np.asarray(eigvals, dtype=float) omega = np.sqrt(np.clip(eigvals, 0.0, None)) nz = np.sort(omega[eigvals > 1e-9]) # selected positive-mode resolution 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()) if counts.size else 0, "vibration_energy": float(0.5 * np.sum(eigvals)), "n_modes": int(nz.size), } def compute_nodal_pulse(G: Any) -> dict[str, Any]: r"""Read per-NFR capacity and phase synchrony without time integration. The nodal equation :math:`\partial\mathrm{EPI}_i/\partial t=\nu_{f,i}\,\Delta\mathrm{NFR}_i` does not by itself derive an oscillator or identify phase speed with structural capacity. This read-out summarizes the stored capacities and phases; it does not observe a period or prove a phase-evolution law. Invalid capacities/phases and unavailable synchrony calculations raise; this includes disabled affinity on a nonempty graph. They are not interpreted as structural silence or zero alignment. Local phase synchrony measures neighbor alignment, global Kuramoto order ``R`` measures collective alignment, and the U3 phase bound reports admissibility. These are distinct from the auxiliary graph-wave spectrum returned by :func:`compute_emergent_pulse`. A high ``R`` does not identify an engine trajectory with that wave model. See ``theory/FORCED_SUPPORT_BALANCE.md`` section 23. Parameters ---------- G : TNFRGraph Returns ------- dict ``mean_frequency`` / ``frequency_spread`` (structural capacities nu_f: mean and std in Hz_str), ``phase_coherence`` (the collective Kuramoto ``R`` in ``[0, 1]`` -- alignment of the current phases), ``mean_local_resonance`` (mean per-NFR local phase synchrony in ``[0, 1]``), ``resonance_gate`` (the canonical default U3 bound Delta phi_max, not a configured live admission), ``n_pulsing`` (legacy name for NFRs with nu_f > 0, not an observed oscillation count), ``n_nodes``. """ from ..constants.canonical import DELTA_PHI_MAX from ..gamma import kuramoto_R_psi from ..metrics.coherence import coherence_matrix, local_phase_sync_weighted from ..metrics.common import finite_mean_absolute, finite_population_std from ._helpers import get_phase 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 = _nodal_frequencies(G, nodes) for node in nodes: get_phase(G, node) mean_frequency = finite_mean_absolute(vf, name="structural frequency") spread = finite_population_std(vf, name="structural frequency") # Collective alignment of the stored phases (Kuramoto order parameter). phase_coherence = finite_real_scalar(kuramoto_R_psi(G)[0], "phase coherence") if phase_coherence < 0.0: raise ValueError("phase coherence must be nonnegative") # per-NFR resonance: mean local phase synchrony (one shared matrix build) order, W = coherence_matrix(G, _record_history=False) if order is None or W is None: raise ValueError("local phase synchrony is unavailable: coherence is disabled") local = tuple( finite_real_scalar( local_phase_sync_weighted(G, k, nodes_order=order, W_row=W), f"local phase synchrony[{k!r}]", ) for k in order ) if any(value < 0.0 for value in local): raise ValueError("local phase synchrony must be nonnegative") mean_local = finite_mean_absolute(local, name="local phase synchrony") return { "mean_frequency": mean_frequency, "frequency_spread": spread, "phase_coherence": phase_coherence, "mean_local_resonance": mean_local, "resonance_gate": gate, "n_pulsing": int(np.sum(vf > 0.0)), "n_nodes": n, } @dataclass(frozen=True) class DiscreteModeCertificate: r"""Finite structural-spectrum and orthonormal-mode diagnostics. The mode basis belongs to the symmetric normalized Laplacian. Its spectrum is compared with the random-walk Laplacian; the legacy ordered sign-change heuristic is reported separately and is not a validity gate. 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 Legacy name for the numerical λ_1 = 0 check. On an irregular connected graph the L_sym zero mode is proportional to sqrt(row strength), corresponding to uniform EPI after the degree-coordinate transform. spectral_gap : float The second eigenvalue λ_2; it can be zero on disconnected support. matches_diffusion_spectrum : bool The L_sym spectrum equals the diffusion operator (L_rw) spectrum. nodal_domains_grow : bool Legacy name for the ordered-sign heuristic: the first mode has no ordered sign changes and the last has more (True for at most one mode). It does not test intermediate modes or graph nodal domains. 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: """Combine the finite spectrum, orthonormality and zero-mode checks.""" return ( self.spectrum_is_discrete and self.modes_orthonormal and self.has_uniform_zero_mode and self.matches_diffusion_spectrum ) 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"zero eigenvalue λ₁=0={self.has_uniform_zero_mode}, " f"spectral gap λ₂={self.spectral_gap:.4f}, " f"matches diffusion spectrum={self.matches_diffusion_spectrum}, " f"ordered-sign heuristic={self.nodal_domains_grow}; " f"ω_k=√λ_k=[{freqs}]" ) def verify_discrete_modes( G: Any, *, tolerance: float = 1e-9 ) -> DiscreteModeCertificate: r"""Check a finite structural spectrum and its orthonormal mode basis. Compare the symmetric normalized spectrum with the random-walk diffusion spectrum and check a zero eigenvalue. The legacy ``nodal_domains_grow`` field compares ordered sign changes of the first and last modes. This heuristic depends on enumeration and the chosen degenerate eigenbasis; it is neither a graph nodal-domain theorem nor part of mode validity. Parameters ---------- G : TNFRGraph tolerance : float Numerical tolerance for the orthonormality / spectrum checks. Returns ------- DiscreteModeCertificate """ _reject_boolean_numeric(tolerance, "tolerance") 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)) # Retain the legacy first/last ordered-sign heuristic as a separate readout. 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"""Real growth rates r − Re(eigenvalues(diag(νf)·L_rw)). 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: nonnegative damping, possibly multiple frozen modes). Parameters ---------- G : TNFRGraph reaction_rate : float, optional A supplied finite local growth/decay rate r added to every mode. The operator catalog and grammar do not derive this auxiliary reaction law or identify an arbitrary operator with a change in r. Returns ------- np.ndarray Growth rates sorted by ascending nodal decay rate. """ rate = finite_real_scalar(reaction_rate, "reaction_rate") return rate - relaxation_spectrum(G) def instability_threshold(G: Any) -> float: r"""Second nodal decay rate, the first nonuniform reaction threshold. On a connected symmetric graph with common positive frequency this is νf·λ_2 and its geometry mode is Fiedler. Heterogeneous capacity uses the actual nodal generator. Multiple zero modes give a zero threshold. For any positive r the constant mode already grows, so this is not a global boundedness certificate or an operator-sequence U2 test. Parameters ---------- G : TNFRGraph Returns ------- float Second-smallest real decay rate (0 with fewer than two modes). """ rates = relaxation_spectrum(G) return float(rates[1]) if len(rates) > 1 else 0.0 def fiedler_partition(G: Any) -> tuple[list, list]: r"""Partition nodes by the sign of the second normalized spatial mode. On a connected symmetric nonnegative graph this is a Fiedler sign partition, useful as a candidate regional decomposition. It is not a general minimum-cut theorem or evidence that a persistent pattern forms. Disconnected graphs can have a second zero mode; degenerate eigenvalues make the returned partition dependent on the selected eigenbasis. Parameters ---------- G : TNFRGraph Returns ------- (part_a, part_b) : tuple[list, list] Positive and nonpositive sign groups of the returned second mode. """ 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 certificate requires a connected symmetric graph and a common positive capacity. It does not certify grammar U2. 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. This homogeneous Fiedler certificate requires symmetric adjacency, a connected positive-weight support, at least two nodes, and one common positive frequency. For heterogeneous/directed damping use :func:`relaxation_spectrum` and :func:`dispersion_relation`. Parameters ---------- G : TNFRGraph tolerance : float Numerical tolerance. Returns ------- StructuralStabilityCertificate """ _reject_boolean_numeric(tolerance, "tolerance") if not np.isfinite(tolerance) or tolerance <= 0.0: raise ValueError("tolerance must be finite and positive") frequency = _nodal_frequencies(G) if ( len(frequency) < 2 or not np.all(frequency == frequency[0]) or frequency[0] <= 0.0 ): raise ValueError( "Fiedler stability requires at least two nodes and a common positive frequency" ) eigvals = _cached_eigenvalues(G) n = len(eigvals) nu = float(frequency[0]) gap = float(eigvals[1]) if gap <= tolerance: raise ValueError("Fiedler stability requires connected positive-weight support") 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 # Independently check the modal readout against the actual nodal matrix. _, lap = structural_diffusion_operator(G) rates = np.linalg.eigvals(frequency[:, None] * lap).real 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).""" conductance = read_conductance(G) return conductance.nodes, conductance.dense(), conductance.strength def _symmetric_adjacency_degree(G: Any) -> tuple[list, Any, Any]: """Require the reversible conductance model used by resistance geometry.""" conductance = read_conductance(G, symmetric=True) return conductance.nodes, conductance.dense(), conductance.strength def random_walk_matrix(G: Any) -> tuple[list, Any]: r"""Return the random-walk transition matrix P = I - L_rw. 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); nodes with zero outgoing strength are absorbing. This includes isolated nodes, directed sinks and nodes with only zero-weight edges. Parallel edge conductances are summed by the shared adjacency. Parameters ---------- G : TNFRGraph Returns ------- (nodes, P) : tuple[list, np.ndarray] """ conductance = read_conductance(G) return conductance.nodes, conductance.transition() def stationary_distribution(G: Any) -> tuple[list, Any]: r"""A stationary distribution of a walk with symmetric conductances. Uses row strengths ``π_i = d_i / Σ d`` when any conductance is present, and the uniform distribution when every node is absorbing. The empty graph returns an empty vector. This is a stationary measure, not a claim of uniqueness or convergence on disconnected or periodic walks. Asymmetric adjacency is rejected: its stationary measure is generally not the degree distribution. The separate directed-diffusion stationary solver supports the strictly positive measure on strongly connected directed networks. Self-loop conductance is counted once, as in P. Scaled strengths avoid overflowing row sums and total volume. Stationary probabilities below the floating-point range round to zero. Parameters ---------- G : TNFRGraph Returns ------- (nodes, pi) : tuple[list, np.ndarray] """ conductance = read_conductance(G, symmetric=True) nodes = conductance.nodes scaled = conductance.relative_strength() if np.any(scaled > 0.0): with np.errstate(under="ignore"): pi = scaled / scaled.sum() else: pi = np.full(len(nodes), 1.0 / len(nodes)) if nodes else scaled return nodes, pi def _resistance_geometry(G: Any) -> tuple[list, Any, Any, list]: """Return dimensionless resistance, component scales and exact-capable volumes. Work component by component so unreachable pairs stay infinite and the pseudoinverse's numerical rank is scaled locally. Connectivity is determined by positive conductance, including summed parallel edges. Non-loop weights set the scale: loops change holding time, not resistance. A volume too large for float is retained as a Fraction internally; its ratio with the conductance scale can still give a finite commute time. """ import math from fractions import Fraction import networkx as nx conductance = read_conductance(G, symmetric=True) nodes, adjacency = conductance.nodes, conductance.dense() resistance = np.full(adjacency.shape, np.inf) scales = np.ones(len(nodes)) volumes = [0.0] * len(nodes) for component in nx.connected_components(nx.from_numpy_array(adjacency)): indices = np.asarray(sorted(component), dtype=int) block = np.ix_(indices, indices) weights = adjacency[block] positive_weights = weights[weights > 0.0] try: volume = math.fsum(positive_weights) except OverflowError: volume = sum( (Fraction.from_float(float(weight)) for weight in positive_weights), Fraction(), ) for index in indices: volumes[index] = volume weights = weights.copy() np.fill_diagonal(weights, 0.0) scale = float(np.max(weights, initial=0.0)) or 1.0 positive = weights > 0.0 with np.errstate(under="ignore"): weights /= scale if np.any(positive & (weights == 0.0)): raise ValueError( "Resistance conductance ratios are below floating-point range" ) scales[indices] = scale laplacian = np.diag(weights.sum(axis=1)) - weights inverse = np.linalg.pinv(laplacian, hermitian=True) diagonal = np.diag(inverse) resistance[block] = np.maximum( diagonal[:, None] + diagonal[None, :] - 2.0 * inverse, 0.0, ) np.fill_diagonal(resistance, 0.0) return nodes, resistance, scales, volumes def _resistance_scale_ratios(numerators: Any, denominators: Any) -> list: """Retain an exact exceptional ratio until the requested matrix is formed.""" from fractions import Fraction result = [] for numerator, denominator in zip(numerators, denominators): if isinstance(numerator, Fraction): result.append(numerator / Fraction.from_float(float(denominator))) else: with np.errstate(over="ignore", under="ignore"): ratio = float(numerator) / float(denominator) if np.isfinite(ratio) and (ratio != 0.0 or numerator == 0.0): result.append(ratio) else: result.append( Fraction.from_float(float(numerator)) / Fraction.from_float(float(denominator)) ) return result def _commute_from_resistance(resistance: Any, volumes: Any) -> Any: """Apply finite or exact scale factors; only unreachable pairs stay infinite.""" finite = np.isfinite(resistance) try: with np.errstate(over="raise", invalid="raise", under="ignore"): factors = np.asarray(volumes, dtype=float) result = np.multiply( factors[:, None], resistance, out=np.full_like(resistance, np.inf), where=finite, ) if not np.all(np.isfinite(result[finite])): raise OverflowError except (FloatingPointError, OverflowError): from fractions import Fraction result = np.full_like(resistance, np.inf) for row, factor in enumerate(volumes): exact_factor = ( factor if isinstance(factor, Fraction) else Fraction.from_float(float(factor)) ) for column in np.flatnonzero(finite[row]): exact = ( Fraction.from_float(float(resistance[row, column])) * exact_factor ) try: result[row, column] = float(exact) except OverflowError as exc: raise ValueError( "Resistance read-out exceeds finite floating-point range" ) from exc if np.any(finite & (resistance > 0.0) & (result == 0.0)): raise ValueError("Resistance read-out is below floating-point range") return result 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 within the pair's connected component. R_eff is an extended metric: different positive-conductance components have infinite resistance, including distinct isolated nodes. Self-resistance is zero. The conductance matrix must be symmetric; directed asymmetric resistance is outside this formula's scope. Component scaling avoids raw degree/volume overflow and removes loops before forming the Laplacian. Unrepresentable requested finite entries raise ValueError; infinity is reserved for unreachable pairs. This remains a floating-point pseudoinverse, with its numerical rank-resolution limit. Parameters ---------- G : TNFRGraph Returns ------- (nodes, R) : tuple[list, np.ndarray] ``R[i, j]`` is the effective resistance between node i and node j. """ nodes, resistance, scales, _ = _resistance_geometry(G) return nodes, _commute_from_resistance( resistance, _resistance_scale_ratios(np.ones(len(nodes)), scales), ) def commute_time(G: Any) -> tuple[list, Any]: r"""Commute-time matrix C(i,j) = volume(component)·R_eff(i,j). The expected round-trip time of the structural random walk between two reachable nodes equals their component's total row strength times the effective resistance. For a connected, unweighted, loopless graph the volume is 2m. Self-loop conductances contribute once to row strengths and account for holding time. Unreachable pairs have infinite commute time; self-commute time is zero. Symmetric conductance is required. A common conductance scale cancels before raw resistance or volume is materialized, so their separate overflow need not prevent a finite result. A requested reachable commute time outside float range raises ValueError. Parameters ---------- G : TNFRGraph Returns ------- (nodes, C) : tuple[list, np.ndarray] """ nodes, resistance, scales, volumes = _resistance_geometry(G) return nodes, _commute_from_resistance( resistance, _resistance_scale_ratios(volumes, scales), ) @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). For symmetric conductance its stationary measure follows row strength, and resistance / commute time describe transport within components. Disconnected pairs have infinite distance. This certificate retains its historical field names; ``2m`` means the pair's component volume, and an edgeless graph uses uniform stationarity. Attributes ---------- n_nodes : int operator_is_walk_generator : bool L_rw = I − P exactly (the diffusion operator generates the walk). transition_row_stochastic : bool Every row of P sums to one, including absorbing zero-strength rows. stationary_is_degree : bool The normalized degree measure (uniform if edgeless) obeys π·P = π. resistance_is_metric : bool R_eff is symmetric, non-negative, and obeys the triangle inequality, allowing infinity between different components. max_resistance : float The largest pairwise effective resistance (transport diameter). commute_equals_2m_resistance : bool C(i,j) = component volume·R_eff(i,j) for reachable pairs; both quantities are infinite for unreachable pairs. 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 measure={self.stationary_is_degree}, " f"R_eff extended metric={self.resistance_is_metric} " f"(max R {self.max_resistance:.4f}), " f"commute=component volume·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 proportional to row strength (uniform if edgeless), that effective resistance is an extended transport metric, and that commute time uses each component's volume. Requires symmetric conductance. Stationarity does not assert convergence or uniqueness; triangle checks are sampled. Parameters ---------- G : TNFRGraph tolerance : float Numerical tolerance. Returns ------- RandomWalkCertificate """ _reject_boolean_numeric(tolerance, "tolerance") if not np.isfinite(tolerance) or tolerance <= 0.0: raise ValueError("tolerance must be finite and positive") nodes, scaled_resistance, scales, volumes = _resistance_geometry(G) n = len(nodes) r = _commute_from_resistance( scaled_resistance, _resistance_scale_ratios(np.ones(n), scales), ) _, lrw = structural_diffusion_operator(G) _, 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) # Every row is stochastic, including absorbing zero-strength rows. row_sums = p.sum(axis=1) row_stochastic = bool(np.all(np.abs(row_sums - 1.0) < tolerance)) # stationary distribution π = degree, π·P = π _, pi = stationary_distribution(G) stationary_ok = ( bool( np.allclose(pi @ p, pi, atol=tolerance, rtol=0.0) and np.all(pi >= 0.0) and abs(float(pi.sum()) - 1.0) < tolerance ) if n else True ) # effective resistance is a metric 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 # Reuse this call's admitted component geometry. The comparison still uses # separately materialized resistance, while commute cancels scales first. # The volume belongs to the pair's component, never to unreachable nodes. c = _commute_from_resistance( scaled_resistance, _resistance_scale_ratios(volumes, scales), ) commute_ok = bool( np.allclose( c, _commute_from_resistance(r, volumes), atol=tolerance, rtol=0.0, ) ) 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 edge-current matrix J_ij = W_ij (EPI_i − EPI_j). The transport carries a current: along each edge i∼j the diffusion flux is J_ij = W_ij (EPI_i − EPI_j) (flux from high to low concentration). The matrix is **antisymmetric** (J_ij = −J_ji) and supported on positive-conductance edges only. Symmetric conductance is required; parallel edges add their conductances and self-loops carry zero current. Asymmetric directed transport is outside this contract. 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. """ conductance, _, flux = _read_edge_flux(G) return conductance.nodes, conductance.dense(flux) 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. At positive row strength d_i and frequency νf_i, the nodal equation gives ``(d_i / νf_i) ∂EPI_i/∂t + div(J)_i = 0``. This constitutive flux does not itself include mobility: a frozen node can carry nonzero pressure/current while its nodal derivative vanishes. Parameters ---------- G : TNFRGraph Returns ------- (nodes, div) : tuple[list, np.ndarray] ``div[i]`` is the net current leaving node i. """ conductance, _, flux = _read_edge_flux(G) return conductance.nodes, conductance.divergence(flux) @dataclass(frozen=True) class StructuralFlowCertificate: r"""Verification of the structural flow (current, Kirchhoff, Ohm). Symmetric diffusion transport carries the current J_ij = W_ij (EPI_i − EPI_j). 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). The latter checks sampled reachable pairs and the injected current equation, not injections between disconnected components. 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), with L·V equal to the injected source/sink vector, for sampled reachable pairs. 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 = W_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). Conductance must be symmetric. Up to 20 distinct reachable pairs are sampled reproducibly; their voltage fields must satisfy L·V = b as well as the resistance drop. When no distinct pair is reachable, this sampled Ohm check is vacuous. Parameters ---------- G : TNFRGraph tolerance : float Numerical tolerance. Returns ------- StructuralFlowCertificate """ _reject_boolean_numeric(tolerance, "tolerance") if not np.isfinite(tolerance) or tolerance <= 0.0: raise ValueError("tolerance must be finite and positive") 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 = _symmetric_adjacency_degree(G) lap = np.diag(deg) - adj epi = structural_field(G, nodes) 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 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: _, resistance = effective_resistance(G) pairs = np.argwhere(np.triu(np.isfinite(resistance), k=1)) rng = np.random.default_rng(0) if len(pairs) > 20: pairs = pairs[rng.choice(len(pairs), size=20, replace=False)] inverses = {} for s, t in pairs: component = tuple(np.flatnonzero(np.isfinite(resistance[s]))) indices = np.asarray(component, dtype=int) if component not in inverses: inverses[component] = np.linalg.pinv( lap[np.ix_(indices, indices)], hermitian=True, ) b = np.zeros(n) b[s], b[t] = 1.0, -1.0 v = np.zeros(n) v[indices] = inverses[component] @ b[indices] drop = v[s] - v[t] if not ( np.allclose(lap @ v, b, atol=tolerance, rtol=0.0) and np.isclose(drop, resistance[s, t], atol=tolerance, rtol=tolerance) ): 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, )