"""TNFR structural conservation diagnostics. The module records the balance Δρ/Δt + div J = S, ρ = Φ_s + K_φ, J = (J_φ, J_ΔNFR), between two observed graph states. ``S`` is the measured source/residual of that finite-interval balance. This identity follows from the definition of ``S``; grammar-valid operator labels alone do not prove that it vanishes or is small. Likewise, the non-negative quadratic structural energy exposed here is a Lyapunov *candidate*: monotonicity requires trajectory evidence or a model-specific proof. These sums use the engine's normalized numeric field convention. With raw units, pressure has EPI units, the inverse-square potential also carries inverse distance squared, and phase discrepancies are angles. They cannot be added without declared reference scales. The neighbor-difference divergence also supplies no time-rate factor: ``dt`` and the current normalization must share a declared time convention before this diagnostic can represent a physical continuity equation. The implementation does not derive those scales. The exact conservation result owned elsewhere is the degree-weighted EPI total for fixed symmetric pure diffusion under its stated capacity assumptions. It must not be conflated with the tetrad charge reported by this module. STATUS ====== CANONICAL DIAGNOSTIC INTERFACE. Residuals, charge drift, and energy change are observations of the supplied trajectory, not consequences of U1-U6 by label. References ---------- - Nodal identity and complete-law premises: theory/FUNDAMENTAL_THEORY.md - Grammar U1-U6: theory/UNIFIED_GRAMMAR_RULES.md - Structural fields: src/tnfr/physics/canonical.py - Extended fields: src/tnfr/physics/extended.py """ from __future__ import annotations import math from collections.abc import Mapping from copy import deepcopy from dataclasses import dataclass, field from dataclasses import fields as dataclass_fields from dataclasses import replace from fractions import Fraction from typing import Any, Sequence import networkx as nx # noqa: F401 - retain the existing module attribute from .._exact_time import finite_represented_real from ..constants.canonical import ( K_PHI_CANONICAL_THRESHOLD, PI, U6_STRUCTURAL_POTENTIAL_LIMIT, ) from ..mathematics._exact_weighted import exact_weighted_sum_ratio from ..mathematics._neighbor_differences import mean_neighbor_difference from ..mathematics.unified_numerical import np from ..metrics.common import ( finite_mean_absolute, finite_pearson_correlation, finite_population_std, ) from ._helpers import finite_real_scalar from .canonical import compute_phase_curvature, compute_structural_potential from .extended import compute_dnfr_flux, compute_phase_current from .unified import ( _capture_structural_fields, _total_charge_from_density, _total_energy_from_fields, ) # --------------------------------------------------------------------------- # Conservation diagnostic alert levels # --------------------------------------------------------------------------- _BALANCE_RMS_ALERT = 1.0 _SECTOR_IMBALANCE_RATIO = 1.5 # --------------------------------------------------------------------------- # Data structures # --------------------------------------------------------------------------- class ConservationAlertLevels(dict[str, float]): """Backward-compatible numeric alert mapping with explicit scope metadata. The historical API returned a plain ``dict[str, float]`` from :func:`compute_grammar_conservation_bounds`. This subclass preserves that mapping behaviour and its legacy keys while exposing why the values must not be read as mathematical bounds or as grammar-rule classifiers. """ scope = "finite_structural_balance_alerts" thresholds_are_proven_bounds = False grammar_validation_applicable = False applicable_grammar_rules: tuple[str, ...] = () u6_drift_requires_reference = True @property def metadata(self) -> dict[str, Any]: """Return detached scope metadata without changing numeric iteration.""" return { "scope": self.scope, "thresholds_are_proven_bounds": self.thresholds_are_proven_bounds, "grammar_validation_applicable": self.grammar_validation_applicable, "applicable_grammar_rules": self.applicable_grammar_rules, "u6_drift_requires_reference": self.u6_drift_requires_reference, "legacy_pi_scaled_alerts": True, } @dataclass(frozen=True) class ConservationSnapshot: """Single-time snapshot of structural balance fields at every node. Attributes ---------- charge_density : dict[Any, float] ρ(i) = Φ_s(i) + K_φ(i) at each node. phi_s : dict[Any, float] Structural potential Φ_s per node. k_phi : dict[Any, float] Phase curvature K_φ per node. j_phi : dict[Any, float] Phase current J_φ per node. j_dnfr : dict[Any, float] ΔNFR flux J_ΔNFR per node. grad_phi : dict[Any, float] Phase gradient |∇φ| per node. divergence : dict[Any, float] Discrete divergence div(J) at each node. divergence_phi, divergence_dnfr : dict or None Captured applications of that same graph operator to each current. Legacy manually constructed snapshots may omit them; their sector decomposition is then unavailable, rather than inferred from magnitudes. """ charge_density: dict[Any, float] phi_s: dict[Any, float] k_phi: dict[Any, float] j_phi: dict[Any, float] j_dnfr: dict[Any, float] grad_phi: dict[Any, float] divergence: dict[Any, float] divergence_phi: dict[Any, float] | None = None divergence_dnfr: dict[Any, float] | None = None def _observed_scalar(value: Any, name: str) -> float: """Admit represented evidence, retaining this module's ValueError API.""" try: return finite_represented_real(value, name)[0] except TypeError as exc: raise ValueError(str(exc)) from exc def _positive_interval(dt: float) -> float: value = _observed_scalar(dt, "dt") if value <= 0.0: raise ValueError("dt must be strictly positive") return value def _observed_secant(before: float, after: float, dt: float, name: str) -> float: """Compute a finite represented secant without false zero or range loss. Ordinary differences retain their historical rounding. Exceptional subtraction/division uses the shared exact difference kernel and rounds only the final quotient; a nonzero unrepresentable rate is unavailable. """ before = _observed_scalar(before, name) after = _observed_scalar(after, name) result = (after - before) / dt if math.isfinite(result) and (result != 0.0 or before == after): return result numerator, denominator = exact_weighted_sum_ratio( (1.0,), (after,), center=before, normalize=False ) time_num, time_den = dt.as_integer_ratio() return _observed_scalar( Fraction(numerator * time_den, denominator * time_num), name ) def _snapshot_nodes(snapshot: ConservationSnapshot, fields: Sequence[str]) -> tuple: """Admit complete finite maps without fabricating missing observations.""" if not isinstance(snapshot, ConservationSnapshot): raise TypeError("expected a ConservationSnapshot") nodes = tuple(snapshot.charge_density) if not nodes: raise ValueError("conservation observations require nonempty node support") support = set(nodes) for name in fields: values = getattr(snapshot, name) if not isinstance(values, Mapping) or set(values) != support: raise ValueError(f"snapshot {name} must match the complete node support") for node in nodes: _observed_scalar(values[node], f"snapshot {name}[{node!r}]") return nodes def _snapshot_pair_nodes(before, after, fields: Sequence[str]) -> tuple: nodes = _snapshot_nodes(before, fields) if set(nodes) != set(_snapshot_nodes(after, fields)): raise ValueError("snapshot comparison requires identical node support") return nodes def _rms(values: Sequence[float]) -> float: """Reduce finite residuals without squaring their absolute scale.""" scale = max(map(abs, values), default=0.0) if scale == 0.0: return 0.0 return scale * math.sqrt( math.fsum((value / scale) ** 2 for value in values) / len(values) ) @dataclass class ConservationBalance: """Result of the continuity equation verification across two snapshots. Uses **Crank-Nicolson (trapezoidal)** discretization for O(Δt²) accuracy: * ``residual[i] = Δρ(i)/Δt + ½[div J_before(i) + div J_after(i)]`` — a finite-interval balance diagnostic, not a grammar validator. * ``mean_residual``, ``max_residual`` — aggregate diagnostics. * ``conservation_quality`` — represented ``1/(1+RMS)`` in [0, 1]. Rounding can produce one for a small nonzero residual; inspect the RMS itself when testing represented zero. * ``grammar_violation_index`` — legacy name for the mean absolute residual; zero does not certify grammar compliance. """ residual: dict[Any, float] delta_rho: dict[Any, float] divergence_after: dict[Any, float] # trapezoidal average (before+after)/2 mean_residual: float std_residual: float max_residual: float rms_residual: float conservation_quality: float grammar_violation_index: float total_charge_before: float total_charge_after: float charge_drift: float diagnostic_scope: str = "two_snapshot_structural_balance" grammar_validation_applicable: bool = False assessed_grammar_rules: tuple[str, ...] = () @property def balance_alert_index(self) -> float: """Return the mean absolute residual under an accurate public name.""" return self.grammar_violation_index @dataclass class ConservationTimeSeries: """Full time-series of conservation diagnostics. Built incrementally via :meth:`ConservationTracker.record`. """ times: list[float] = field(default_factory=list) total_charge: list[float] = field(default_factory=list) mean_residuals: list[float] = field(default_factory=list) rms_residuals: list[float] = field(default_factory=list) conservation_quality: list[float] = field(default_factory=list) grammar_violation_index: list[float] = field(default_factory=list) charge_drift: list[float] = field(default_factory=list) @property def is_conserved(self) -> bool: """Return the legacy finite-series quality alert classification. This alias neither proves a conservation law nor validates grammar. Prefer :attr:`aggregate_balance_within_alert` in new code. """ return self.aggregate_balance_within_alert @property def aggregate_balance_within_alert(self) -> bool: """Whether measured interval quality reaches the legacy alert cut. The synthetic baseline stored for the first snapshot has no preceding interval and is excluded. """ if len(self.conservation_quality) <= 1: return False return self.sampled_mean_quality >= 0.9 @property def grammar_validation_applicable(self) -> bool: """Grammar cannot be inferred from this residual time series.""" return False @property def sampled_mean_quality(self) -> float: """Average quality over actual two-snapshot balance intervals.""" if len(self.conservation_quality) <= 1: return 0.0 return float(np.mean(self.conservation_quality[1:])) @property def mean_quality(self) -> float: """Compatibility name for :attr:`sampled_mean_quality`.""" return self.sampled_mean_quality @property def mean_quality_including_baseline(self) -> float: """Return the historical mean including the synthetic first value.""" if not self.conservation_quality: return 0.0 return float(np.mean(self.conservation_quality)) # --------------------------------------------------------------------------- # Core computation: structural charge density and divergence # --------------------------------------------------------------------------- def compute_charge_density(G: Any) -> dict[Any, float]: r"""Compute structural charge density ρ(i) = Φ_s(i) + K_φ(i). This is the Noether-like charge-density diagnostic used by TNFR: - Φ_s captures global structural potential (long-range coupling) - K_φ captures local geometric curvature (short-range confinement) Their sum combines global and local structural fields. Conservation of its total must be checked on the actual trajectory or derived for a specified auxiliary model. Parameters ---------- G : TNFRGraph Returns ------- dict[node, float] ρ(i) per node. """ phi_s = compute_structural_potential(G) k_phi = compute_phase_curvature(G) return _charge_density_from_fields(phi_s, k_phi) def _charge_density_from_fields( phi_s: dict[Any, float], k_phi: dict[Any, float] ) -> dict[Any, float]: """Shared charge definition for live and captured field maps.""" return {n: phi_s[n] + k_phi.get(n, 0.0) for n in phi_s} def compute_current_divergence(G: Any) -> dict[Any, float]: r"""Compute discrete divergence of structural current div J(i). The structural current is J = (J_φ, J_ΔNFR). On a graph, the stored quantity at node i uses the legacy neighbor-minus-center sign: div J(i) = (1/|N(i)|) Σ_{j∈N(i)} [ (J_φ(j) - J_φ(i)) + (J_ΔNFR(j) - J_ΔNFR(i)) ] Thus it is ``-L_rw`` applied to each current component for an unweighted graph, conventionally an inward rather than outward flux. The public name is retained for compatibility; all balance routines use this same sign. This is a normalized graph diagnostic, with no metric-length divisor or transport-rate factor. Summing its phase and pressure components assumes declared field reference scales, not equality of their raw physical units. Parameters ---------- G : TNFRGraph Returns ------- dict[node, float] div J(i) per node. """ j_phi = compute_phase_current(G) j_dnfr = compute_dnfr_flux(G) return _current_divergence_from_fields(G, j_phi, j_dnfr) def _current_divergence_from_fields( G: Any, j_phi: dict[Any, float], j_dnfr: dict[Any, float] ) -> dict[Any, float]: """Apply the existing neighbor-mean divergence to recorded currents.""" phase, pressure = _current_divergence_components_from_fields(G, j_phi, j_dnfr) return {node: phase[node] + pressure[node] for node in phase} def _current_divergence_components_from_fields(G, j_phi, j_dnfr): """Read each linear current channel on the actual graph neighborhood.""" nodes = tuple(G.nodes()) phase, pressure = {}, {} for i in nodes: neighbors = list(G.neighbors(i)) # Legacy neighbor-minus-center (inward-flux / -L_rw) convention. phase[i] = mean_neighbor_difference(j_phi[i], [j_phi[j] for j in neighbors]) pressure[i] = mean_neighbor_difference( j_dnfr[i], [j_dnfr[j] for j in neighbors] ) return phase, pressure # --------------------------------------------------------------------------- # Snapshot capture # --------------------------------------------------------------------------- def capture_conservation_snapshot(G: Any) -> ConservationSnapshot: """Capture structural charge, currents, and fields at one instant. This is a *read-only* operation that never mutates EPI. The caller must hold graph state fixed during the complete call, including the topology-dependent divergence. Returned maps are detached; capture is not atomic with concurrent evolution. Parameters ---------- G : TNFRGraph Network in its current state. Returns ------- ConservationSnapshot """ fields = _capture_structural_fields(G) phi_s, k_phi, grad_phi = fields.phi_s, fields.k_phi, fields.grad_phi j_phi, j_dnfr = fields.j_phi, fields.j_dnfr charge = _charge_density_from_fields(phi_s, k_phi) div_phi, div_dnfr = _current_divergence_components_from_fields(G, j_phi, j_dnfr) div_j = {node: div_phi[node] + div_dnfr[node] for node in div_phi} return ConservationSnapshot( charge_density=charge, phi_s=phi_s, k_phi=k_phi, j_phi=j_phi, j_dnfr=j_dnfr, grad_phi=grad_phi, divergence=div_j, divergence_phi=div_phi, divergence_dnfr=div_dnfr, ) # --------------------------------------------------------------------------- # Conservation balance (two-snapshot comparison) # --------------------------------------------------------------------------- def verify_conservation_balance( before: ConservationSnapshot, after: ConservationSnapshot, dt: float = 1.0, ) -> ConservationBalance: r"""Verify the structural continuity equation between two snapshots. Computes the residual of the continuity equation using a **Crank-Nicolson (trapezoidal)** discretization: Δρ(i)/Δt + ½[div J_before(i) + div J_after(i)] ≈ 0 For sufficiently smooth evolution on fixed support, the trapezoidal approximation has O(Δt²) accuracy, compared to the O(Δt) of a right-endpoint scheme. A small residual records agreement with this balance over the sampled interval; it does not prove sequence-wide conservation or grammar compliance. Large residuals can reflect numerical error, changing topology or node support, source terms, or a failure of the assumed balance. Their cause requires separate investigation. Requires nonempty identical node support, complete finite charge/divergence maps and a finite positive interval. A support change needs explicit source and correspondence accounting outside this observer; it is never filled with missing-node zeros. All inputs use a common declared normalization and time convention. The routine does not convert the raw units of potential, curvature, currents or ``dt``; a small numeric residual alone does not validate those units. A nonzero rate outside binary64 range raises instead of becoming a zero residual. Recoverable intermediate subtraction overflow uses an exact represented difference before the final quotient is rounded. Parameters ---------- before : ConservationSnapshot State before operator application. after : ConservationSnapshot State after operator application. dt : float Step in the declared diagnostic time coordinate (default 1.0). Returns ------- ConservationBalance Finite-interval structural-balance diagnostics. """ dt = _positive_interval(dt) nodes = _snapshot_pair_nodes(before, after, ("charge_density", "divergence")) delta_rho: dict[Any, float] = {} residual: dict[Any, float] = {} average_divergence: dict[Any, float] = {} for n in nodes: rho_before = before.charge_density[n] rho_after = after.charge_density[n] d_rho = _observed_secant(rho_before, rho_after, dt, "charge rate") delta_rho[n] = d_rho # Crank-Nicolson: trapezoidal average of divergence at both endpoints div_j = mean_neighbor_difference( 0.0, (before.divergence[n], after.divergence[n]) ) average_divergence[n] = div_j # Continuity: ∂ρ/∂t + div J = S → residual = ∂ρ/∂t + div J residual[n] = finite_real_scalar(d_rho + div_j, "balance residual") residual_vals = tuple(residual.values()) mean_res = mean_neighbor_difference(0.0, residual_vals) std_res = finite_population_std(residual_vals, name="balance residual") max_res = max(map(abs, residual_vals)) rms_res = _rms(residual_vals) # Conservation quality: 1/(1 + RMS) maps [0, ∞) → (0, 1] quality = 1.0 / (1.0 + rms_res) # Legacy field name. This is only the mean absolute balance residual and # has no rule-classification semantics. gvi = finite_mean_absolute(residual_vals, name="balance residual") # Total charge tracking q_before = _total_charge_from_density(before.charge_density) q_after = _total_charge_from_density(after.charge_density) charge_drift = finite_real_scalar(abs(q_after - q_before), "charge drift") return ConservationBalance( residual=residual, delta_rho=delta_rho, divergence_after=average_divergence, mean_residual=mean_res, std_residual=std_res, max_residual=max_res, rms_residual=rms_res, conservation_quality=quality, grammar_violation_index=gvi, total_charge_before=q_before, total_charge_after=q_after, charge_drift=charge_drift, ) # --------------------------------------------------------------------------- # Conservation tracker (multi-step) # --------------------------------------------------------------------------- class ConservationTracker: """Track finite structural-balance diagnostics across an operator sequence. Usage ----- >>> tracker = ConservationTracker(G) >>> tracker.record(t=0.0) # initial snapshot >>> Emission()(G, node) >>> tracker.record(t=1.0) # after operator >>> Coherence()(G, node) >>> tracker.record(t=2.0) >>> report = tracker.report() >>> print(f"Within alert: {report.aggregate_balance_within_alert}") """ def __init__(self, G: Any) -> None: self._G = G self._snapshots: list[tuple[float, ConservationSnapshot]] = [] self._series = ConservationTimeSeries() def record( self, t: float = 0.0, *, _require_energy: bool = False ) -> ConservationSnapshot: """Capture current state and compute balance against previous snapshot. Parameters ---------- t : float Finite timestamp, strictly greater than the previous timestamp. Equal timestamps do not imply an invented unit-duration step. Returns ------- ConservationSnapshot A detached copy of the captured snapshot. Editing it cannot alter the tracker's retained evidence. Rejected samples leave the retained snapshots and series unchanged. The private energy flag lets the SDK admit its combined balance/energy report before commit; pure balance tracking does not require energy. """ t = _observed_scalar(t, "snapshot timestamp") if self._snapshots: t_prev, snap_prev = self._snapshots[-1] dt = _positive_interval(t - t_prev) snap = capture_conservation_snapshot(self._G) _snapshot_nodes(snap, ("charge_density", "divergence")) if self._snapshots: balance = verify_conservation_balance(snap_prev, snap, dt=dt) if _require_energy: compute_lyapunov_derivative(snap_prev, snap, dt=dt) # Commit only after timestamp, capture and interval admission pass. self._snapshots.append((t, snap)) self._series.times.append(t) self._series.total_charge.append(balance.total_charge_after) self._series.mean_residuals.append(balance.mean_residual) self._series.rms_residuals.append(balance.rms_residual) self._series.conservation_quality.append(balance.conservation_quality) self._series.grammar_violation_index.append(balance.grammar_violation_index) self._series.charge_drift.append(balance.charge_drift) else: # First snapshot — record initial charge only if _require_energy: _energy_from_snapshot(snap) q_total = _total_charge_from_density(snap.charge_density) self._snapshots.append((t, snap)) self._series.times.append(t) self._series.total_charge.append(q_total) self._series.mean_residuals.append(0.0) self._series.rms_residuals.append(0.0) self._series.conservation_quality.append(1.0) self._series.grammar_violation_index.append(0.0) self._series.charge_drift.append(0.0) # Copy maps without copying node identifiers: labels may use object # identity, and they must still address the caller's original nodes. return replace( snap, **{ descriptor.name: dict(values) for descriptor in dataclass_fields(snap) if (values := getattr(snap, descriptor.name)) is not None }, ) def report(self) -> ConservationTimeSeries: """Return detached accumulated diagnostics, not editable retained evidence.""" return deepcopy(self._series) @property def latest_balance(self) -> ConservationBalance | None: """Return the most recent balance check, or None.""" if len(self._snapshots) < 2: return None t_prev, snap_prev = self._snapshots[-2] t_curr, snap_curr = self._snapshots[-1] dt = t_curr - t_prev return verify_conservation_balance(snap_prev, snap_curr, dt=dt) # --------------------------------------------------------------------------- # Decomposed analysis: which field component contributes most to residual # --------------------------------------------------------------------------- def decompose_conservation_residual( before: ConservationSnapshot, after: ConservationSnapshot, dt: float = 1.0, ) -> dict[str, dict[Any, float]]: r"""Decompose the continuity residual into Φ_s and K_φ contributions. The charge density ρ = Φ_s + K_φ, so: Δρ/Δt = ΔΦ_s/Δt + ΔK_φ/Δt This function separates potential-field and curvature-field contributions to the measured residual. Their magnitudes do not classify U2 or U3, and potential drift evaluates the U6 policy only when compared with the canonical two-snapshot threshold through a dedicated U6 checker. Divergence is evaluated using the **Crank-Nicolson (trapezoidal)** average of the before and after snapshots for O(Δt²) accuracy. Its sector components are captured independently on the graph. Local current magnitudes cannot recover their divergences or even their signs. Snapshots missing this component evidence must be recaptured. Returns ------- dict with keys: 'phi_s_drift' : per-node ΔΦ_s/Δt 'k_phi_drift' : per-node ΔK_φ/Δt 'j_phi_div' : per-node div(J_φ) contribution 'j_dnfr_div' : per-node div(J_ΔNFR) contribution 'potential_residual' : ΔΦ_s/Δt + div(J_ΔNFR) [potential sector] 'geometric_residual' : ΔK_φ/Δt + div(J_φ) [geometric sector] """ dt = _positive_interval(dt) if any( getattr(snapshot, field) is None for snapshot in (before, after) for field in ("divergence_phi", "divergence_dnfr") ): raise ValueError("sector divergences are unavailable; recapture both snapshots") nodes = _snapshot_pair_nodes( before, after, ("phi_s", "k_phi", "divergence_phi", "divergence_dnfr"), ) phi_s_drift: dict[Any, float] = {} k_phi_drift: dict[Any, float] = {} j_phi_div: dict[Any, float] = {} j_dnfr_div: dict[Any, float] = {} potential_residual: dict[Any, float] = {} geometric_residual: dict[Any, float] = {} for n in nodes: # Rate of change of each charge component d_phi_s = _observed_secant( before.phi_s[n], after.phi_s[n], dt, "potential rate" ) d_k_phi = _observed_secant( before.k_phi[n], after.k_phi[n], dt, "curvature rate" ) phi_s_drift[n] = d_phi_s k_phi_drift[n] = d_k_phi j_phi_div[n] = mean_neighbor_difference( 0.0, (before.divergence_phi[n], after.divergence_phi[n]) ) j_dnfr_div[n] = mean_neighbor_difference( 0.0, (before.divergence_dnfr[n], after.divergence_dnfr[n]) ) potential_residual[n] = finite_real_scalar( d_phi_s + j_dnfr_div[n], "potential residual" ) geometric_residual[n] = finite_real_scalar( d_k_phi + j_phi_div[n], "geometric residual" ) return { "phi_s_drift": phi_s_drift, "k_phi_drift": k_phi_drift, "j_phi_div": j_phi_div, "j_dnfr_div": j_dnfr_div, "potential_residual": potential_residual, "geometric_residual": geometric_residual, } # --------------------------------------------------------------------------- # Legacy pi-scaled conservation alert levels # --------------------------------------------------------------------------- def compute_grammar_conservation_bounds(G: Any) -> ConservationAlertLevels: r"""Compute legacy grammar-scaled diagnostic alert levels. The legacy alert construction combines: - U2's stabilizer/debt role as qualitative context; - the U6 drift policy ``ΔPhi_s < pi/2`` as a numeric scale; - the exact trigonometric bound ``|J_phi| <= 1``. These values combine historical pi-scaled policies into monitoring alert levels. They are not derived U2/U3/U6 bounds, do not evaluate any grammar rule, and must not be used as a conservation or grammar certificate. The returned object remains a ``dict`` subclass containing only numeric legacy keys; its :attr:`ConservationAlertLevels.metadata` property exposes this scope without breaking callers that iterate the numeric mapping. Returns ------- ConservationAlertLevels 'max_charge_density' : legacy policy-scaled alert level for |ρ| 'max_current_magnitude' : legacy alert level for |J| 'max_allowed_residual' : legacy alert level for |Δρ/Δt + div J| 'phi_s_confinement' : legacy key for the π/2 U6 *drift alert scale* 'k_phi_hotspot' : 0.9×π ≈ 2.8274 (curvature hotspot threshold) """ n_nodes = G.number_of_nodes() # Legacy key/value retained. U6 constrains two-snapshot mean |ΔΦ_s|; # it does not bound the magnitude |Φ_s| in one graph. phi_s_alert = U6_STRUCTURAL_POTENTIAL_LIMIT # K_φ is bounded by π (wrapped angle difference) k_phi_bound = PI # Historical composite alert, not a maximum charge-density theorem. max_charge = phi_s_alert + k_phi_bound # J_φ = mean(sin(Δθ)), bounded by 1 j_phi_bound = 1.0 # Legacy J_ΔNFR alert proxy. U2/U6 do not bound pressure spread. j_dnfr_alert = 2.0 * phi_s_alert # Maximum current magnitude max_current = math.sqrt(j_phi_bound**2 + j_dnfr_alert**2) # Historical divergence alert derived from the composite current scale. avg_degree = 2.0 * G.number_of_edges() / max(n_nodes, 1) max_div = 2.0 * max_current # Historical residual alert, not a maximum allowed by grammar. max_residual = max_charge + max_div return ConservationAlertLevels( { "max_charge_density": max_charge, "max_current_magnitude": max_current, "max_allowed_residual": max_residual, "phi_s_confinement": phi_s_alert, "k_phi_hotspot": K_PHI_CANONICAL_THRESHOLD, "average_degree": avg_degree, } ) # --------------------------------------------------------------------------- # Legacy grammar-named entry point for balance alerts # --------------------------------------------------------------------------- def detect_grammar_violations_from_conservation( balance: ConservationBalance, bounds: dict[str, float] | None = None, ) -> dict[str, Any]: r"""Report legacy-scaled balance alerts without classifying grammar. Residuals and charge drift contain no operator-history evidence and no direct phase-admissibility evidence. They therefore cannot classify U2, U3, or U6, and cannot validate U1-U6. Numerical discretization, topology changes, source terms, or an external pressure law can all produce the same values for grammatically different histories. The function name and the legacy ``violations_*`` keys remain for API compatibility. Those keys now report that no grammar verdict was made. New ``alerts_*`` keys contain the actual finite-balance alerts. Parameters ---------- balance : ConservationBalance Result from verify_conservation_balance. bounds : dict[str, float], optional Legacy alert levels from :func:`compute_grammar_conservation_bounds`. Returns ------- dict with: ``violations_detected`` is always ``False`` because grammar is not assessed; ``alerts_detected``, ``alert_types``, ``severity`` and ``nodes_alerted`` expose the legacy-scaled observations. Scope and applicability metadata are included explicitly. """ threshold = U6_STRUCTURAL_POTENTIAL_LIMIT if bounds is not None: threshold = bounds.get("max_allowed_residual", U6_STRUCTURAL_POTENTIAL_LIMIT) threshold = _observed_scalar(threshold, "max_allowed_residual") if threshold <= 0.0: raise ValueError("max_allowed_residual must be strictly positive") alert_types: list[str] = [] nodes_alerted: list[Any] = [] for node, res in balance.residual.items(): if abs(res) > threshold: nodes_alerted.append(node) if nodes_alerted: alert_types.append("balance_residual_above_legacy_alert") if balance.charge_drift > U6_STRUCTURAL_POTENTIAL_LIMIT: alert_types.append("charge_drift_above_legacy_pi_alert") if balance.rms_residual > _BALANCE_RMS_ALERT: alert_types.append("balance_rms_above_legacy_alert") if balance.max_residual > 2 * threshold: alert_types.append("balance_peak_above_legacy_alert") severity = min(1.0, balance.rms_residual / max(threshold, 1e-10)) return { # Backward-compatible grammar-shaped fields: no grammar assessment was # made, so these must never carry inferred U-rule failures. "violations_detected": False, "violation_count": 0, "violation_types": [], "nodes_violating": [], # Accurate finite-balance alert surface. "alerts_detected": bool(alert_types), "alert_count": len(alert_types), "alert_types": alert_types, "severity": severity, "nodes_alerted": nodes_alerted, "alert_threshold": threshold, # Applicability metadata. "diagnostic_scope": "two_snapshot_structural_balance_alerts", "grammar_validation_applicable": False, "grammar_validated": False, "grammar_rules_assessed": (), "thresholds_are_proven_bounds": False, "u6_drift_requires_phi_s_reference": True, } # --------------------------------------------------------------------------- # Structural charge (historical Noether-like name) # --------------------------------------------------------------------------- def compute_noether_charge(G: Any) -> float: r"""Compute the historically named tetrad charge candidate. ``Q = Σ_i ρ(i) = Σ_i [Φ_s(i) + K_φ(i)]``. The charge Q integrates global (potential) and local (geometric) structural information into a single scalar. Its drift must be measured along the actual trajectory; grammar compliance is neither sufficient nor inferred from a small drift. This tetrad charge candidate is distinct from the EPI-channel degree-weighted total ``Σ_i deg(i)·EPI(i)`` (:func:`tnfr.physics.structural_diffusion.degree_weighted_total`). The latter has an exact restricted conservation theorem for fixed symmetric random-walk diffusion; this function provides no corresponding theorem for the full tetrad dynamics. Parameters ---------- G : TNFRGraph Returns ------- float Total structural-charge candidate. """ charge = compute_charge_density(G) return _total_charge_from_density(charge) def compute_energy_functional(G: Any) -> float: r"""Compute the TNFR structural energy functional. E = (1/2) Σ_i ℰ(i) where ℰ(i) = Φ_s² + |∇φ|² + K_φ² + J_φ² + J_ΔNFR² is the raw energy density defined by :func:`unified.compute_energy_density`. **Single source of truth**: delegates to the shared normalized-total helper in ``unified`` for the same quadratic form and ½ normalization. It avoids materializing unrepresentable squares when the final total remains representable. **Mathematical normalization**: ``E = sum(variational.compute_hamiltonian_density(G).values())`` ``E = 0.5 * sum(unified.compute_energy_density(G).values())`` These equalities identify the real quadratic form, not bitwise equality after independent per-node rounding. Raw densities can underflow/overflow even when the normalized network total is representable. The result is non-negative. It is a Lyapunov candidate only: U2 validity alone does not prove ``dE/dt <= 0`` for an arbitrary pressure law or operator sequence. Parameters ---------- G : TNFRGraph Returns ------- float Total structural energy. See Also -------- unified.compute_energy_density : Raw ℰ(i) per node. variational.compute_hamiltonian_density : H(i) = ½·ℰ(i) per node. """ fields = _capture_structural_fields(G) return _total_energy_from_fields( fields.phi_s, fields.grad_phi, fields.k_phi, fields.j_phi, fields.j_dnfr ) def _energy_from_snapshot(snapshot: ConservationSnapshot) -> float: r"""E = ½ Σ_i (Φ_s² + |∇φ|² + K_φ² + J_φ² + J_ΔNFR²) from a captured snapshot. Snapshot-field counterpart of :func:`compute_energy_functional` (which reads the live graph): both evaluate the same canonical quadratic form, so the energy-density definition lives in one place rather than being inlined. """ _snapshot_nodes(snapshot, ("phi_s", "grad_phi", "k_phi", "j_phi", "j_dnfr")) return _total_energy_from_fields( snapshot.phi_s, snapshot.grad_phi, snapshot.k_phi, snapshot.j_phi, snapshot.j_dnfr, ) # --------------------------------------------------------------------------- # Sector coupling analysis (the deep physics insight) # --------------------------------------------------------------------------- def analyze_sector_coupling( before: ConservationSnapshot, after: ConservationSnapshot, dt: float = 1.0, ) -> dict[str, Any]: r"""Analyze coupling between potential and geometric conservation sectors. The measured balance decomposes into two residual sectors: **Potential sector** (global, ΔNFR-driven): ∂Φ_s/∂t + div(J_ΔNFR) ≈ 0 - Records pressure-field transport and source mismatch - Does not classify the U2 stabilizer/debt rule - U6 needs a two-snapshot mean |ΔΦ_s| comparison **Geometric sector** (local, phase-driven): ∂K_φ/∂t + div(J_φ) ≈ 0 - Records phase-curvature transport and source mismatch - Does not assess U3 edge-phase admissibility - The 0.9π hotspot cut is a selected alert inside the exact π wrap The cross-correlation summarizes co-variation between the two residual channels. It does not establish a causal coupling or derive the complex field ``Psi = K_phi + i*J_phi``. The legacy three-sample and 1e-15 dispersion cuts are retained as observation policies, not physical bounds. An unavailable correlation keeps its compatibility zero but is explicitly tagged; a nonrepresentable asymmetry ratio is returned as None. Parameters ---------- before, after : ConservationSnapshot Snapshots before and after an operator sequence. dt : float Time step. Returns ------- dict with: 'potential_sector_residual' : RMS residual of potential sector 'geometric_sector_residual' : RMS residual of geometric sector 'cross_coupling_strength' : Correlation between sector residuals 'dominant_sector' : 'potential' | 'geometric' | 'balanced' 'sector_asymmetry' : Ratio of dominant to subdominant residual """ decomp = decompose_conservation_residual(before, after, dt=dt) nodes = list(decomp["potential_residual"].keys()) pot_res = [float(decomp["potential_residual"][n]) for n in nodes] geo_res = [float(decomp["geometric_residual"][n]) for n in nodes] rms_pot = _rms(pot_res) rms_geo = _rms(geo_res) # Cross-coupling: correlation between sector residuals if len(nodes) <= 2: correlation_status = "insufficient_samples" cross_corr = None elif ( finite_population_std(pot_res) <= 1e-15 or finite_population_std(geo_res) <= 1e-15 ): correlation_status = "below_legacy_dispersion_cut" cross_corr = None else: cross_corr = finite_pearson_correlation(pot_res, geo_res) correlation_status = ( "available" if cross_corr is not None else "constant_sample" ) # Determine dominant sector if rms_pot > _SECTOR_IMBALANCE_RATIO * rms_geo: dominant = "potential" elif rms_geo > _SECTOR_IMBALANCE_RATIO * rms_pot: dominant = "geometric" else: dominant = "balanced" asymmetry = max(rms_pot, rms_geo) / (min(rms_pot, rms_geo) + 1e-15) return { "potential_sector_residual": rms_pot, "geometric_sector_residual": rms_geo, "cross_coupling_strength": cross_corr if cross_corr is not None else 0.0, "cross_coupling_available": cross_corr is not None, "cross_coupling_status": correlation_status, "dominant_sector": dominant, "sector_asymmetry": asymmetry if math.isfinite(asymmetry) else None, "sector_asymmetry_available": math.isfinite(asymmetry), } # --------------------------------------------------------------------------- # Ward-like diagnostics: per-operator finite-step signatures # --------------------------------------------------------------------------- @dataclass(frozen=True) class WardIdentity: """Conservation signature of a single operator application. This legacy-named object records finite before/after observables for one reported operator application. It does not prove a symmetry or validate the operator sequence. For operator O_k at step k: ΔQ_k/(N Δt_k) + mean_i[(div J_before + div J_after)/2]_i = mean_i[S_k(i)] Thus ``delta_charge`` is a total structural-charge-candidate change while ``mean_source`` is a per-node rate residual. They are distinct quantities. Attributes ---------- operator_name : str Name of the applied operator (e.g. "AL", "IL", "OZ"). delta_charge : float Total structural-charge-candidate change ΔQ = Q_after - Q_before. delta_energy : float Energy functional change ΔE = E_after - E_before. mean_source : float Network-averaged source term ⟨S⟩ (continuity residual). conservation_quality : float Balance quality for this single step. charge_character : str Finite-difference label: 'source' (ΔQ > ε), 'sink' (ΔQ < -ε), 'transport' (|ΔQ| < ε), or legacy 'exact' (both changes within ε). The last label means threshold-neutral, not exact conservation. energy_character : str 'dissipative' (ΔE < -ε), 'injective' (ΔE > ε), or 'neutral'. """ operator_name: str delta_charge: float delta_energy: float mean_source: float conservation_quality: float charge_character: str energy_character: str def compute_ward_identity( before: ConservationSnapshot, after: ConservationSnapshot, operator_name: str, G_before: Any = None, G_after: Any = None, dt: float = 1.0, threshold: float = 0.01, ) -> WardIdentity: r"""Compute a legacy-named Ward diagnostic for one observed step. Measures structural-charge change from the supplied snapshots. Energy also comes from those snapshots unless both optional graphs are supplied. The labels are finite threshold classifications and contain no grammar verdict. Parameters ---------- before, after : ConservationSnapshot Network state before and after operator application. operator_name : str Name/glyph of the operator that was applied (e.g. "AL", "IL"). G_before, G_after : TNFRGraph, optional If both are supplied, they provide the energy endpoints. Correspondence of those graph states to the charge snapshots is the caller's unverified responsibility; neither support nor field equality is checked. Supplying only one graph retains the legacy snapshot-energy fallback. Snapshot energy evaluates the same declared quadratic form, not an estimate of an independently supplied graph. dt : float Time step between snapshots (default 1.0). threshold : float Minimum |ΔQ| or |ΔE| to classify as non-neutral (default 0.01). Returns ------- WardIdentity """ if not isinstance(operator_name, str) or not operator_name.strip(): raise ValueError("operator_name must be a non-empty string") if isinstance(threshold, (bool, np.bool_)): raise TypeError("threshold must be a finite positive real number") threshold = float(threshold) if not np.isfinite(threshold) or threshold <= 0.0: raise ValueError("threshold must be finite and strictly positive") balance = verify_conservation_balance(before, after, dt=dt) q_before = balance.total_charge_before q_after = balance.total_charge_after delta_q = q_after - q_before # Energy from graphs if available, else estimate from snapshots if G_before is not None and G_after is not None: e_before = compute_energy_functional(G_before) e_after = compute_energy_functional(G_after) else: e_before = _energy_from_snapshot(before) e_after = _energy_from_snapshot(after) delta_e = e_after - e_before # Classify charge character if abs(delta_q) < threshold and abs(delta_e) < threshold: charge_char = "exact" elif delta_q > threshold: charge_char = "source" elif delta_q < -threshold: charge_char = "sink" else: charge_char = "transport" # Classify energy character if delta_e < -threshold: energy_char = "dissipative" elif delta_e > threshold: energy_char = "injective" else: energy_char = "neutral" return WardIdentity( operator_name=operator_name, delta_charge=delta_q, delta_energy=delta_e, mean_source=balance.mean_residual, conservation_quality=balance.conservation_quality, charge_character=charge_char, energy_character=energy_char, ) def verify_sequence_ward_identity( identities: Sequence[WardIdentity], *, alert_level: float | None = None, ) -> dict[str, Any]: r"""Measure the sequence Ward residual ``Σ_k ``. The legacy ``sequence_conserved`` flag applies a finite-sequence alert level. It is retained as an alias for ``aggregate_balance_within_alert``; it does not infer grammar validity or prove a conservation law. When no alert level is supplied, the historical pi-scaled value is retained for compatibility and is explicitly reported as a legacy alert, not a bound. Parameters ---------- identities : Sequence[WardIdentity] Ordered Ward identities for each operator in the sequence. alert_level : float, optional Finite aggregate-source alert. Defaults to the historical ``(pi/2) / n_steps`` scale. Returns ------- dict with: 'total_source' : float — Σ⟨S_k⟩ (should be ≈ 0) 'total_charge_change' : float — net ΔQ 'total_energy_change' : float — net ΔE 'sequence_conserved' : legacy alias for a threshold comparison 'operator_summary' : dict[str, int] — count by charge_character """ total_source = finite_real_scalar( math.fsum(finite_real_scalar(w.mean_source, "mean source") for w in identities), "total source", ) total_dq = finite_real_scalar( math.fsum( finite_real_scalar(w.delta_charge, "charge change") for w in identities ), "total charge change", ) total_de = finite_real_scalar( math.fsum( finite_real_scalar(w.delta_energy, "energy change") for w in identities ), "total energy change", ) summary: dict[str, int] = {} for w in identities: summary[w.charge_character] = summary.get(w.charge_character, 0) + 1 n_steps = max(len(identities), 1) if alert_level is None: threshold = U6_STRUCTURAL_POTENTIAL_LIMIT / n_steps else: if isinstance(alert_level, (bool, np.bool_)): raise TypeError("alert_level must be a finite positive real number") threshold = float(alert_level) if not np.isfinite(threshold) or threshold <= 0.0: raise ValueError("alert_level must be finite and strictly positive") within_alert = bool(identities) and abs(total_source) < threshold return { "total_source": total_source, "total_charge_change": total_dq, "total_energy_change": total_de, "sequence_conserved": within_alert, "aggregate_balance_within_alert": within_alert, "alert_threshold": threshold, "thresholds_are_proven_bounds": False, "grammar_validation_applicable": False, "grammar_validated": False, "diagnostic_scope": "finite_sequence_balance_aggregate", "sample_available": bool(identities), "operator_summary": summary, } # --------------------------------------------------------------------------- # Lyapunov stability analysis # --------------------------------------------------------------------------- @dataclass(frozen=True) class LyapunovResult: """Result of Lyapunov stability analysis for an operator step. The energy functional E = ½Σ(Φ_s² + |∇φ|² + K_φ² + J_φ² + J_ΔNFR²) is a Lyapunov candidate. This result reports whether one observed finite step decreases it; grammar validity alone does not determine that sign. Attributes ---------- energy_before : float E[G] before the operator application. energy_after : float E[G] after the operator application. energy_derivative : float (E_after - E_before) / dt — approximation of dE/dt. dissipation : float D[G] = max(0, -dE/dt) — structural dissipation rate. is_stable : bool Compatibility alert: dE/dt <= the configured nonnegative tolerance. Use ``energy_nonincreasing`` for the captured endpoint-energy comparison. is_strongly_stable : bool True when dE/dt < -ε (energy strictly decreasing). """ energy_before: float energy_after: float energy_derivative: float dissipation: float is_stable: bool is_strongly_stable: bool @property def energy_nonincreasing(self) -> bool: """Compare the captured endpoint energies, independently of rate rounding.""" return self.energy_after <= self.energy_before def compute_lyapunov_derivative( before: ConservationSnapshot, after: ConservationSnapshot, dt: float = 1.0, stability_threshold: float = 1e-6, ) -> LyapunovResult: r"""Compute the Lyapunov derivative dE/dt between two snapshots. The structural energy functional: E = ½ Σ_i [Φ_s(i)² + |∇φ|(i)² + K_φ(i)² + J_φ(i)² + J_ΔNFR(i)²] The structural dissipation readout is ``D[G] = max(0, -dE/dt)``. A non-increasing observation supports stability for that step; this function does not prove a general Lyapunov theorem from grammar labels. Parameters ---------- before, after : ConservationSnapshot Network state before and after operator application. dt : float Time step (default 1.0). stability_threshold : float Nonnegative numerical allowance for ``is_stable`` and the minimum decrease for ``is_strongly_stable`` (default 1e-6). It is not a theorem. Returns ------- LyapunovResult """ dt = _positive_interval(dt) stability_threshold = _observed_scalar(stability_threshold, "stability_threshold") if stability_threshold < 0.0: raise ValueError("stability_threshold must be nonnegative") _snapshot_pair_nodes(before, after, ("charge_density",)) e_before = _energy_from_snapshot(before) e_after = _energy_from_snapshot(after) de_dt = _observed_secant(e_before, e_after, dt, "energy derivative") dissipation = max(0.0, -de_dt) return LyapunovResult( energy_before=e_before, energy_after=e_after, energy_derivative=de_dt, dissipation=dissipation, is_stable=de_dt <= stability_threshold, is_strongly_stable=de_dt < -stability_threshold, ) # --------------------------------------------------------------------------- # Spectral conservation analysis (graph Laplacian decomposition) # --------------------------------------------------------------------------- @dataclass(frozen=True) class SpectralConservation: r"""Spectral decomposition of the conservation fields. Expands charge density ρ and current divergence in the eigenbasis of the symmetric normalized Laplacian L_sym = I − D^{-1/2} W D^{-1/2}, whose spectrum is that of the canonical TNFR diffusion operator L_rw = I − D⁻¹W (the EPI channel of the nodal equation): ρ(i) = Σ_k ρ̂_k ψ_k(i). The continuity equation mode-by-mode reads: dρ̂_k/dt + λ_k Ĵ_k = Ŝ_k Eigenvalue ordering supplies graph scales only. A static snapshot cannot infer conservation at low frequency, relaxation at high frequency, U5 compliance, or the omitted time derivative. Attributes ---------- eigenvalues : np.ndarray L_sym eigenvalues λ_k (sorted ascending) — the canonical L_rw relaxation spectrum. rho_spectrum : np.ndarray Charge density coefficients ρ̂_k in the eigenbasis. div_spectrum : np.ndarray Current divergence coefficients in the eigenbasis. modal_divergence_magnitude : np.ndarray Static activity ``|div_hat_k|`` of the already-computed divergence. No second Laplacian factor is applied. low_divergence_activity_modes : int Number of modes with divergence magnitude at or below its median. spectral_gap : float λ_1 — gap between zero mode and first non-trivial mode. node_order : tuple Node order used to build both field vectors and the Laplacian. """ eigenvalues: Any # np.ndarray rho_spectrum: Any # np.ndarray div_spectrum: Any # np.ndarray conservation_by_mode: Any # np.ndarray dominant_conservation_modes: int spectral_gap: float node_order: tuple[Any, ...] = () @property def modal_divergence_magnitude(self) -> Any: """Accurate name for the legacy stored per-mode activity field.""" return self.conservation_by_mode @property def low_divergence_activity_modes(self) -> int: """Accurate name for the legacy stored median-split count.""" return self.dominant_conservation_modes def compute_spectral_conservation( G: Any, snapshot: ConservationSnapshot | None = None, ) -> SpectralConservation: r"""Decompose conservation fields in the normalized Laplacian eigenbasis. Represents one charge-density and current-divergence snapshot in the ``L_sym`` eigenbasis. ``modal_divergence_magnitude=|div_hat_k|`` is a static activity readout. The divergence has already been computed in node space, so multiplying it by ``lambda_k`` again would apply an unintended second graph derivative. The compatibility field ``conservation_by_mode`` exposes the same array. With no ``d rho_hat / dt``, neither name is a per-mode conservation residual or a U5 assessment. Parameters ---------- G : TNFRGraph The TNFR network. snapshot : ConservationSnapshot, optional Pre-computed snapshot. If ``None``, captured from *G*. Returns ------- SpectralConservation """ if snapshot is None: snapshot = capture_conservation_snapshot(G) # Build the symmetric normalized Laplacian L_sym = I − D^{-1/2} W D^{-1/2}, # whose spectrum is that of the canonical TNFR diffusion operator # L_rw = I − D⁻¹W (the EPI channel; see ``structural_diffusion``). This is # The field vectors must use exactly the node order returned with L. from .structural_diffusion import symmetric_normalized_laplacian nodes, L = symmetric_normalized_laplacian(G) n = len(nodes) # Eigendecomposition (orthonormal eigenbasis of L_sym) eigvals, eigvecs = np.linalg.eigh(L) # Project charge density and divergence into eigenbasis rho_vec = np.array([snapshot.charge_density[nd] for nd in nodes]) div_vec = np.array([snapshot.divergence[nd] for nd in nodes]) rho_hat = eigvecs.T @ rho_vec # coefficients in eigenbasis div_hat = eigvecs.T @ div_vec # Static per-mode divergence magnitude. ``div_vec`` is already a graph # divergence, so an additional eigenvalue factor would apply L twice. divergence_activity = np.abs(div_hat) # Median split is descriptive only; it is not a conservation verdict. median_res = float(np.median(divergence_activity)) if n > 0 else 0.0 n_low_activity = int(np.sum(divergence_activity <= median_res + 1e-15)) # Spectral gap sorted_eigs = np.sort(eigvals) spectral_gap = float(sorted_eigs[1]) if n > 1 else 0.0 return SpectralConservation( eigenvalues=eigvals, rho_spectrum=rho_hat, div_spectrum=div_hat, conservation_by_mode=divergence_activity, dominant_conservation_modes=n_low_activity, spectral_gap=spectral_gap, node_order=tuple(nodes), ) # --------------------------------------------------------------------------- # Public API # --------------------------------------------------------------------------- __all__ = [ # Data structures "ConservationSnapshot", "ConservationBalance", "ConservationTimeSeries", "ConservationAlertLevels", "WardIdentity", "LyapunovResult", "SpectralConservation", # Core computations "compute_charge_density", "compute_current_divergence", "capture_conservation_snapshot", "verify_conservation_balance", # Tracking "ConservationTracker", # Analysis "decompose_conservation_residual", "analyze_sector_coupling", "compute_grammar_conservation_bounds", "detect_grammar_violations_from_conservation", # Conserved quantities "compute_noether_charge", "compute_energy_functional", # Ward identities "compute_ward_identity", "verify_sequence_ward_identity", # Lyapunov stability "compute_lyapunov_derivative", # Spectral conservation "compute_spectral_conservation", ]