integrators.py

"""Shared nodal integrators and explicitly configured additive forcing.

The unforced nodal law is ``dEPI/dt = vf * dnfr``. The optional Gamma registry
implements the broader configured equation:

    ∂EPI/∂t = νf · ΔNFR(t) + Γi(R)

The extended equation includes:
  - Base term: νf · ΔNFR(t) - canonical structural evolution
  - Additive term: Γi - optional Kuramoto-based or harmonic forcing

Gamma is an additional rate source, not a derivation from the nodal product.
At positive capacity it can be written as an effective pressure ``dnfr+Gamma/vf``;
this representation is singular at zero capacity. Nonzero Gamma can therefore
move EPI at ``vf=0``. Unforced certificates and Silence claims require Gamma
to vanish, or separate evidence for the extended model. No units or phase
synchronization theorem follows merely from calling this integrator.

The base term is ``vf * dnfr`` with stored frequency and pressure held fixed
during each call. For the built-in time-only forcing at fixed phases, ``rk4``
is fourth-order quadrature. It does not reevaluate a state-dependent pressure
law at Runge-Kutta stages. Optional clipping can alter the unconstrained ODE
trajectory, so the order claim applies while clipping is inactive.

These solvers evolve the signed real scalar EPI chart, including its exact
uniform-real BEPI embedding. Richer BEPI values require their own dynamics;
their magnitude read-out must not silently replace the evolved state.
"""

from __future__ import annotations

import math
import sys
from abc import ABC, abstractmethod
from collections.abc import Callable, Iterable, Iterator, Mapping
from concurrent.futures import ProcessPoolExecutor
from contextlib import contextmanager
from fractions import Fraction
from multiprocessing import get_context
from typing import Any, Literal, TypeAlias, cast

import networkx as nx

from .._exact_time import finite_represented_real
from ..alias import get_attr, get_attr_str, set_attr, set_attr_str
from ..config.parsing import parse_bool
from ..constants import DEFAULTS
from ..constants.aliases import (
    ALIAS_D2EPI,
    ALIAS_DEPI,
    ALIAS_DNFR,
    ALIAS_EPI,
    ALIAS_EPI_KIND,
    ALIAS_THETA,
    ALIAS_VF,
)
from ..constants.canonical import (
    INTEGRATORS_DNFR_BOUNDS_CANONICAL,
    INTEGRATORS_HALF_STEP_CANONICAL,
    INTEGRATORS_RK4_SIXTH_CANONICAL,
    INTEGRATORS_SIGMOID_OFFSET_CANONICAL,
)
from ..errors.contextual import NetworkConfigError, TNFRUserError, TNFRValueError
from ..gamma import (
    _get_gamma_spec,
    _uses_builtin_gamma,
    eval_gamma,
    eval_gamma_vectorized,
)
from ..mathematics._exact_weighted import exact_weighted_sum_ratio
from ..mathematics._neighbor_differences import mean_neighbor_difference
from ..mathematics.unified_numerical import np
from ..types import NodeId, TNFRGraph, require_finite_real_scalar_epi
from ..utils import resolve_chunk_size
from ._euler_kernel import euler_update
from .structural_clip import resolve_clip_policy, structural_clip, structural_clip_array

__all__ = (
    "AbstractIntegrator",
    "DefaultIntegrator",
    "prepare_integration_params",
    "rk4_coupled_step",
    "update_epi_via_nodal_equation",
)

GammaMap: TypeAlias = dict[NodeId, float]
"""Γ evaluation cache keyed by node identifier."""

NodeIncrements: TypeAlias = dict[NodeId, tuple[float, ...]]
"""Mapping of nodes to staged integration increments."""

NodalUpdate: TypeAlias = dict[NodeId, tuple[float, float, float]]
"""Mapping of nodes to ``(EPI, dEPI/dt, ∂²EPI/∂t²)`` tuples."""

IntegratorMethod: TypeAlias = Literal["euler", "rk4"]
"""Supported explicit integration schemes for nodal updates."""

ClipPolicy: TypeAlias = tuple[float, float, Literal["hard", "soft"], float]

_PARALLEL_GRAPH: TNFRGraph | None = None


def _read_scalar_epi(nd: dict[str, Any]) -> float:
    """Read the authoritative scalar chart without a magnitude projection."""
    return cast(
        float,
        get_attr(
            nd,
            ALIAS_EPI,
            0.0,
            strict=True,
            conv=require_finite_real_scalar_epi,
        ),
    )


def _finite_output(value: Any, parameter: str = "derivative") -> float:
    """Admit represented-real solver scalars before coercion or commit."""
    try:
        return finite_represented_real(value, parameter)[0]
    except (TypeError, ValueError) as exc:
        raise NetworkConfigError(
            parameter=parameter, value=value, reason=str(exc)
        ) from exc


def _integration_clip_policy(graph: TNFRGraph) -> ClipPolicy:
    """Validate the active boundary policy before any solver-side evaluation."""
    try:
        return resolve_clip_policy(graph.graph)
    except ValueError as exc:
        raise NetworkConfigError(
            parameter="EPI clipping policy",
            value={
                key: graph.graph[key]
                for key in ("EPI_MIN", "EPI_MAX", "CLIP_MODE", "CLIP_SOFT_K")
                if key in graph.graph
            },
            reason=str(exc),
        ) from exc


def _validate_clock_grid(t0: float, dt_step: float, steps: int, method: str) -> float:
    """Validate the repeated-add clock actually used, without storing a grid."""
    current = t0
    for _ in range(steps):
        following = current + dt_step
        if not math.isfinite(following) or (dt_step > 0 and following <= current):
            raise NetworkConfigError(
                parameter="t",
                value=following,
                reason="Every positive solver substep must advance a finite represented clock",
            )
        if method == "rk4" and dt_step > 0:
            midpoint = current + dt_step / INTEGRATORS_HALF_STEP_CANONICAL
            if not math.isfinite(midpoint):
                raise NetworkConfigError(
                    parameter="t",
                    value=midpoint,
                    reason="RK4 stage times must be finite",
                )
        current = following
    return current


@contextmanager
def _restore_solver_outputs_on_error(
    graph: TNFRGraph, keys: tuple[str, ...]
) -> Iterator[None]:
    """Restore owned node outputs after a later scalar/extended substep fails.

    Callback effects, graph caches and external state are not transactional.
    """
    saved = {
        node: {key: nd[key] for key in keys if key in nd}
        for node, nd in graph.nodes(data=True)
    }
    try:
        yield
    except Exception:
        for node, original in saved.items():
            if node not in graph:
                continue
            nd = graph.nodes[node]
            for key in keys:
                nd.pop(key, None)
            nd.update(original)
        raise


def _gamma_worker_init(graph: TNFRGraph) -> None:
    """Initialise process-local graph reference for Γ evaluation."""

    global _PARALLEL_GRAPH
    _PARALLEL_GRAPH = graph


def _gamma_worker(task: tuple[list[NodeId], float]) -> list[tuple[NodeId, float]]:
    """Evaluate Γ for ``task`` chunk using process-local graph."""

    chunk, t = task
    if _PARALLEL_GRAPH is None:
        raise RuntimeError("Parallel Γ worker initialised without graph reference")
    return [
        (node, float(eval_gamma(_PARALLEL_GRAPH, node, t, strict=True)))
        for node in chunk
    ]


def _normalise_jobs(n_jobs: int | None, total: int) -> int | None:
    """Return an effective worker count respecting serial fallbacks."""

    if n_jobs is None:
        return None
    try:
        workers = int(n_jobs)
    except (TypeError, ValueError):
        return None
    if workers <= 1 or total <= 1:
        return None
    return max(1, min(workers, total))


def _chunk_nodes(nodes: list[NodeId], chunk_size: int) -> Iterable[list[NodeId]]:
    """Yield deterministic chunks from ``nodes`` respecting insertion order."""

    for idx in range(0, len(nodes), chunk_size):
        yield nodes[idx : idx + chunk_size]


def _exact_rk4_update(epi: float, dt_step: float, stages: tuple[float, ...]) -> float:
    """Round a range-limited quadrature after exact represented-input arithmetic."""
    before = _finite_output(epi, "EPI")
    step = _finite_output(dt_step, "dt")
    rates = tuple(_finite_output(stage, "RK4 stage") for stage in stages)
    numerator, denominator = exact_weighted_sum_ratio(
        (1.0, 2.0, 2.0, 1.0), rates, normalize=True
    )
    exact = Fraction(before) + Fraction(step) * Fraction(numerator, denominator)
    try:
        result = float(exact)
    except OverflowError as exc:
        raise NetworkConfigError(
            parameter="EPI", value="out of range", reason="Solver output must be finite"
        ) from exc
    return _finite_output(result, "EPI")


def _rk4_update(epi: Any, dt_step: float, stages: tuple[Any, ...]) -> Any:
    """Apply one shared 1:2:2:1 quadrature to scalar or array supplied rates.

    Ordinary values retain the existing separate arithmetic order. If h/6 is
    subnormal or a weighted sum/proposal overflows, reuse the exact weighted
    reducer and round the final represented-input expression for that row.
    This avoids losing a finite response to an intermediate range limit; it
    does not recompute state-dependent stages or add an evolution law.
    """
    k1, k2, k3, k4 = stages
    factor = dt_step / INTEGRATORS_RK4_SIXTH_CANONICAL
    factor_underflow = dt_step != 0.0 and abs(factor) < sys.float_info.min
    if np is not None and isinstance(epi, np.ndarray):
        with np.errstate(over="ignore", invalid="ignore", under="ignore"):
            weighted = k1 + 2 * k2 + 2 * k3 + k4
            result = epi + factor * weighted
        exceptional = ~np.isfinite(weighted) | ~np.isfinite(result)
        if factor_underflow:
            exceptional = np.ones_like(result, dtype=bool)
        for index in np.flatnonzero(exceptional):
            result[index] = _exact_rk4_update(
                epi[index], dt_step, tuple(stage[index] for stage in stages)
            )
        return result

    weighted = k1 + 2 * k2 + 2 * k3 + k4
    result = epi + factor * weighted
    if factor_underflow or not math.isfinite(weighted) or not math.isfinite(result):
        return _exact_rk4_update(epi, dt_step, stages)
    return result


def rk4_coupled_step(
    state: tuple[float, ...],
    rhs: Callable[[float, tuple[float, ...]], tuple[float, ...]],
    *,
    t: float,
    dt: float,
) -> tuple[float, ...]:
    """Advance one explicitly supplied complete vector field by classical RK4.

    Every stage receives an immutable, independently admitted full state and
    reevaluates ``rhs(time, state)``. This opt-in kernel does not change the
    graph integrator's held-pressure semantics. The caller owns the complete
    law, pressure refresh, phase lifts, Gamma admission, boundary conditions,
    numerical error certificate and association with graph state or history.

    Only nonempty tuples are accepted for state and rates. Scalars use the
    shared represented-real admission; dt is positive and all three distinct
    stage clocks must remain finite and strictly ordered after rounding.
    No clipping, wrapping, hidden substeps or graph writes occur. Callback
    effects outside its immutable arguments are not rolled back on failure.
    """
    if not callable(rhs):
        raise TypeError("rhs must be callable")

    def vector(values, label, *, size=None):
        if type(values) is not tuple or not values:
            raise TypeError(f"{label} must be a nonempty tuple")
        if size is not None and len(values) != size:
            raise ValueError(f"{label} dimension differs from the complete state")
        return tuple(
            finite_represented_real(value, f"{label}[{index}]")[0]
            for index, value in enumerate(values)
        )

    initial = vector(state, "state")
    start, _ = finite_represented_real(t, "t")
    step, exact_step = finite_represented_real(dt, "dt")
    if step <= 0:
        raise ValueError("dt must be positive")
    half, _ = finite_represented_real(exact_step / 2, "RK4 half-step")
    midpoint, _ = finite_represented_real(start + half, "RK4 midpoint time")
    endpoint, _ = finite_represented_real(start + step, "RK4 endpoint time")
    if not start < midpoint < endpoint:
        raise ValueError("RK4 stage clocks must remain strictly ordered")

    def evaluate(time, values):
        return vector(rhs(time, values), "rhs", size=len(initial))

    def stage(rates, duration):
        result = []
        for index, (value, rate) in enumerate(zip(initial, rates, strict=True)):
            updated = finite_represented_real(
                euler_update(value, duration, rate), f"RK4 stage[{index}]"
            )[0]
            if (
                updated == 0.0
                and (Fraction(value) + Fraction(duration) * Fraction(rate)) != 0
            ):
                raise ValueError("RK4 stage loses a nonzero represented-input result")
            result.append(updated)
        return tuple(result)

    first = evaluate(start, initial)
    second = evaluate(midpoint, stage(first, half))
    third = evaluate(midpoint, stage(second, half))
    fourth = evaluate(endpoint, stage(third, step))
    result = []
    for index, (value, *rates) in enumerate(
        zip(initial, first, second, third, fourth, strict=True)
    ):
        updated = finite_represented_real(
            _rk4_update(value, step, tuple(rates)), f"RK4 endpoint[{index}]"
        )[0]
        if (
            updated == 0.0
            and (
                Fraction(value)
                + exact_step
                * (
                    Fraction(rates[0])
                    + 2 * Fraction(rates[1])
                    + 2 * Fraction(rates[2])
                    + Fraction(rates[3])
                )
                / 6
            )
            != 0
        ):
            raise ValueError("RK4 endpoint loses a nonzero represented-input result")
        result.append(updated)
    return tuple(result)


def _apply_increment_chunk(
    chunk: list[tuple[NodeId, float, float, tuple[float, ...]]],
    dt_step: float,
    method: str,
) -> list[tuple[NodeId, tuple[float, float, float]]]:
    """Compute updated states for ``chunk`` using scalar arithmetic."""

    results: list[tuple[NodeId, tuple[float, float, float]]] = []
    dt_nonzero = dt_step != 0

    for node, epi_i, dEPI_prev, ks in chunk:
        if method == "rk4":
            epi = _rk4_update(epi_i, dt_step, ks)
            dEPI_dt = ks[3]
        else:
            (k1,) = ks
            epi = euler_update(epi_i, dt_step, k1)
            dEPI_dt = k1
        d2epi = (dEPI_dt - dEPI_prev) / dt_step if dt_nonzero else 0.0
        results.append(
            (
                node,
                (
                    _finite_output(epi, "EPI"),
                    _finite_output(dEPI_dt),
                    _finite_output(d2epi),
                ),
            )
        )

    return results


def _evaluate_gamma_map(
    G: TNFRGraph,
    nodes: list[NodeId],
    t: float,
    *,
    n_jobs: int | None = None,
) -> GammaMap:
    """Return Γ evaluations for ``nodes`` at time ``t`` respecting parallelism."""

    # Arbitrary callbacks may read graph state or keep invocation-side evidence;
    # execute them once on the live staged path, not on process-local copies.
    builtin = _uses_builtin_gamma(_get_gamma_spec(G, strict=True))
    workers = _normalise_jobs(n_jobs, len(nodes)) if builtin else None
    if workers is None:
        return {n: float(eval_gamma(G, n, t, strict=True)) for n in nodes}

    approx_chunk = math.ceil(len(nodes) / (workers * 4)) if workers > 0 else None
    chunk_size = resolve_chunk_size(
        approx_chunk,
        len(nodes),
        minimum=1,
    )
    mp_ctx = get_context("spawn")
    tasks = ((chunk, t) for chunk in _chunk_nodes(nodes, chunk_size))

    results: GammaMap = {}
    with ProcessPoolExecutor(
        max_workers=workers,
        mp_context=mp_ctx,
        initializer=_gamma_worker_init,
        initargs=(G,),
    ) as executor:
        futures = [executor.submit(_gamma_worker, task) for task in tasks]
        for fut in futures:
            for node, value in fut.result():
                results[node] = value
    return results


def prepare_integration_params(
    G: TNFRGraph,
    dt: float | None = None,
    t: float | None = None,
    method: Literal["euler", "rk4"] | None = None,
) -> tuple[float, int, float, Literal["euler", "rk4"]]:
    """Validate and normalise ``dt``, ``t`` and ``method`` for integration.

    Explicit and graph-default timesteps must be finite nonnegative reals.
    Boolean/text values and nonzero inputs lost in materialization reject for
    the step, initial clock and minimum step alike.
    Invalid parameters raise :class:`NetworkConfigError`. When ``dt``
    exceeds a positive ``DT_MIN`` stored on ``G`` the span is deterministically
    subdivided into integer steps so that the resulting ``dt_step`` never falls
    below that minimum threshold.

    Returns ``(dt_step, steps, t0, method)`` where ``dt_step`` is the effective
    step, ``steps`` the number of substeps and ``t0`` the prepared initial
    time.
    """
    if dt is None:
        dt = G.graph.get("DT", DEFAULTS.get("DT", 0.1))
    dt = _finite_output(dt, "dt")
    if not math.isfinite(dt) or dt < 0:
        raise NetworkConfigError(
            parameter="dt", value=dt, reason="Time step must be finite and non-negative"
        )

    t = _finite_output(G.graph.get("_t", 0.0) if t is None else t, "t")
    if not math.isfinite(t):
        raise NetworkConfigError(
            parameter="t", value=t, reason="Initial time must be finite"
        )

    method_value = (
        G.graph.get("INTEGRATOR_METHOD", DEFAULTS.get("INTEGRATOR_METHOD", "euler"))
        if method is None
        else method
    )
    method_value = (
        method_value.lower() if isinstance(method_value, str) else method_value
    )
    if method_value not in ("euler", "rk4"):
        raise NetworkConfigError(
            parameter="method",
            value=method_value,
            reason="Integration method must be 'euler' or 'rk4'",
        )

    dt_min = _finite_output(
        G.graph.get("DT_MIN", DEFAULTS.get("DT_MIN", 0.0)), "DT_MIN"
    )
    if not math.isfinite(dt_min) or dt_min < 0:
        raise NetworkConfigError(
            parameter="DT_MIN",
            value=dt_min,
            reason="Minimum time step must be finite and non-negative",
        )
    steps = 1
    if dt_min > 0 and dt > dt_min:
        ratio = dt / dt_min
        if not math.isfinite(ratio):
            raise NetworkConfigError(
                parameter="DT_MIN",
                value=dt_min,
                reason="Time-step subdivision ratio must be finite",
            )
        steps = max(1, int(math.floor(ratio)))
    dt_step = dt / steps if steps else 0.0
    _validate_clock_grid(t, dt_step, steps, method_value)

    return dt_step, steps, t, cast(Literal["euler", "rk4"], method_value)


def _apply_increments(
    G: TNFRGraph,
    dt_step: float,
    increments: NodeIncrements,
    *,
    method: str,
    n_jobs: int | None = None,
) -> NodalUpdate:
    """Combine precomputed increments to update node states."""

    nodes: list[NodeId] = list(G.nodes)
    if not nodes:
        return {}

    epi_initial: list[float] = []
    dEPI_prev: list[float] = []
    ordered_increments: list[tuple[float, ...]] = []

    for node in nodes:
        nd = G.nodes[node]
        _, _, dEPI_dt_prev, epi_i = _node_state(nd)
        epi_initial.append(float(epi_i))
        dEPI_prev.append(float(dEPI_dt_prev))
        ordered_increments.append(increments[node])

    if np is not None:
        epi_arr = np.asarray(epi_initial, dtype=float)
        dEPI_prev_arr = np.asarray(dEPI_prev, dtype=float)
        k_arr = np.asarray(ordered_increments, dtype=float)

        if method == "rk4":
            if k_arr.ndim != 2 or k_arr.shape[1] != 4:
                raise TNFRUserError(
                    message="RK4 integration requires four staged increments",
                    suggestion="Check integrator implementation logic",
                    context={"shape": str(k_arr.shape)},
                )
            k1 = k_arr[:, 0]
            k2 = k_arr[:, 1]
            k3 = k_arr[:, 2]
            k4 = k_arr[:, 3]
            epi = _rk4_update(epi_arr, dt_step, (k1, k2, k3, k4))
            dEPI_dt = k4
        else:
            if k_arr.ndim == 1:
                k1 = k_arr
            else:
                k1 = k_arr[:, 0]
            epi = euler_update(epi_arr, dt_step, k1)
            dEPI_dt = k1

        if dt_step != 0:
            d2epi = (dEPI_dt - dEPI_prev_arr) / dt_step
        else:
            d2epi = np.zeros_like(dEPI_dt)

        results: NodalUpdate = {}
        for idx, node in enumerate(nodes):
            results[node] = (
                _finite_output(epi[idx], "EPI"),
                _finite_output(dEPI_dt[idx]),
                _finite_output(d2epi[idx]),
            )
        return results

    payload: list[tuple[NodeId, float, float, tuple[float, ...]]] = list(
        zip(nodes, epi_initial, dEPI_prev, ordered_increments)
    )

    workers = _normalise_jobs(n_jobs, len(nodes))
    if workers is None:
        return dict(_apply_increment_chunk(payload, dt_step, method))

    approx_chunk = math.ceil(len(nodes) / (workers * 4)) if workers > 0 else None
    chunk_size = resolve_chunk_size(
        approx_chunk,
        len(nodes),
        minimum=1,
    )
    mp_ctx = get_context("spawn")

    results: NodalUpdate = {}
    with ProcessPoolExecutor(max_workers=workers, mp_context=mp_ctx) as executor:
        futures = [
            executor.submit(
                _apply_increment_chunk,
                chunk,
                dt_step,
                method,
            )
            for chunk in _chunk_nodes(payload, chunk_size)
        ]
        for fut in futures:
            for node, value in fut.result():
                results[node] = value

    return {node: results[node] for node in nodes}


def _collect_nodal_increments(
    G: TNFRGraph,
    gamma_maps: tuple[GammaMap, ...],
    *,
    method: str,
) -> NodeIncrements:
    """Combine node base state with staged Γ contributions.

    Implements the canonical TNFR nodal equation in two parts:

    1. **Base term** (canonical equation):
       base = vf * dnfr  →  ∂EPI/∂t = νf · ΔNFR(t)

       This is the fundamental TNFR equation where:
         - vf (νf): structural frequency in Hz_str
         - dnfr (ΔNFR): nodal gradient (reorganization operator)
         - base: instantaneous rate of EPI evolution

    2. **Network coupling term**:
       Γi(R) from gamma_maps - optional Kuramoto order parameter

    The full extended equation is: ∂EPI/∂t = νf·ΔNFR(t) + Γi(R)

    Args:
        G: TNFR graph with node attributes vf and dnfr
        gamma_maps: Staged Γ evaluations (1 for Euler, 4 for RK4)
        method: Integration method ('euler' or 'rk4')

    Returns:
        Mapping of nodes to staged integration increments

    Pressure units are EPI units divided by capacity and the declared clock.
    This numerical combination does not establish an operator closure theorem.
    """

    nodes: list[NodeId] = list(G.nodes())
    if not nodes:
        return {}

    if method == "rk4":
        expected_maps = 4
    elif method == "euler":
        expected_maps = 1
    else:
        raise TNFRValueError(
            "method must be 'euler' or 'rk4'",
            context={"method": method, "available": ["euler", "rk4"]},
            suggestion="Use 'euler' or 'rk4' as the integration method.",
        )

    if len(gamma_maps) != expected_maps:
        raise TNFRValueError(
            f"{method} integration requires {expected_maps} gamma maps",
            context={
                "method": method,
                "required": expected_maps,
                "provided": len(gamma_maps),
            },
            suggestion=f"Provide exactly {expected_maps} gamma maps for {method} integration.",
        )

    if np is not None:
        states = [_node_state(G.nodes[node]) for node in nodes]
        vf = np.asarray([state[0] for state in states], dtype=float)
        dnfr = np.asarray([state[1] for state in states], dtype=float)
        # CANONICAL TNFR EQUATION: ∂EPI/∂t = νf · ΔNFR(t)
        # This implements the fundamental nodal equation explicitly
        base = vf * dnfr

        gamma_arrays = [
            np.fromiter((gm.get(n, 0.0) for n in nodes), float, count=len(nodes))
            for gm in gamma_maps
        ]
        if gamma_arrays:
            gamma_stack = np.stack(gamma_arrays, axis=1)
            combined = base[:, None] + gamma_stack
        else:
            combined = base[:, None]

        return {
            node: tuple(float(value) for value in combined[idx])
            for idx, node in enumerate(nodes)
        }

    increments: NodeIncrements = {}
    for node in nodes:
        nd = G.nodes[node]
        vf, dnfr, *_ = _node_state(nd)
        # CANONICAL TNFR EQUATION: ∂EPI/∂t = νf · ΔNFR(t)
        # Scalar implementation of the fundamental nodal equation
        base = vf * dnfr
        gammas = [gm.get(node, 0.0) for gm in gamma_maps]

        if method == "rk4":
            k1, k2, k3, k4 = gammas
            increments[node] = (
                base + k1,
                base + k2,
                base + k3,
                base + k4,
            )
        else:
            (k1,) = gammas
            increments[node] = (base + k1,)

    return increments


def _build_gamma_increments(
    G: TNFRGraph,
    dt_step: float,
    t_local: float,
    *,
    method: str,
    n_jobs: int | None = None,
) -> NodeIncrements:
    """Evaluate Γ contributions and merge them with ``νf·ΔNFR`` base terms."""

    if method == "rk4":
        gamma_count = 4
    elif method == "euler":
        gamma_count = 1
    else:
        raise TNFRValueError(
            "method must be 'euler' or 'rk4'",
            context={"method": method, "available": ["euler", "rk4"]},
            suggestion="Use 'euler' or 'rk4' as the integration method.",
        )

    gamma_spec = _get_gamma_spec(G, strict=True)

    gamma_type = ""
    if isinstance(gamma_spec, Mapping):
        gamma_type = str(gamma_spec.get("type", ""))

    if gamma_type == "none" and _uses_builtin_gamma(gamma_spec):
        gamma_maps: tuple[GammaMap, ...] = tuple(
            cast(GammaMap, {}) for _ in range(gamma_count)
        )
        return _collect_nodal_increments(G, gamma_maps, method=method)

    nodes: list[NodeId] = list(G.nodes)
    if not nodes:
        gamma_maps = tuple(cast(GammaMap, {}) for _ in range(gamma_count))
        return _collect_nodal_increments(G, gamma_maps, method=method)

    if method == "rk4":
        t_mid = t_local + dt_step / INTEGRATORS_HALF_STEP_CANONICAL
        t_end = t_local + dt_step
        g1_map = _evaluate_gamma_map(G, nodes, t_local, n_jobs=n_jobs)
        g_mid_map = _evaluate_gamma_map(G, nodes, t_mid, n_jobs=n_jobs)
        g4_map = _evaluate_gamma_map(G, nodes, t_end, n_jobs=n_jobs)
        gamma_maps = (g1_map, g_mid_map, g_mid_map, g4_map)
    else:  # method == "euler"
        gamma_maps = (_evaluate_gamma_map(G, nodes, t_local, n_jobs=n_jobs),)

    return _collect_nodal_increments(G, gamma_maps, method=method)


def _integrate_euler(
    G: TNFRGraph,
    dt_step: float,
    t_local: float,
    *,
    n_jobs: int | None = None,
) -> NodalUpdate:
    """One explicit Euler integration step."""
    increments = _build_gamma_increments(
        G,
        dt_step,
        t_local,
        method="euler",
        n_jobs=n_jobs,
    )
    return _apply_increments(
        G,
        dt_step,
        increments,
        method="euler",
        n_jobs=n_jobs,
    )


def _integrate_rk4(
    G: TNFRGraph,
    dt_step: float,
    t_local: float,
    *,
    n_jobs: int | None = None,
) -> NodalUpdate:
    """Fourth-order forcing quadrature with the stored nodal base held fixed."""
    increments = _build_gamma_increments(
        G,
        dt_step,
        t_local,
        method="rk4",
        n_jobs=n_jobs,
    )
    return _apply_increments(
        G,
        dt_step,
        increments,
        method="rk4",
        n_jobs=n_jobs,
    )


def _integrate_vectorized_step(
    G: TNFRGraph,
    dt_step: float,
    steps: int,
    t0: float,
    method: str,
    np: Any,
    clip_policy: ClipPolicy,
    states: list[tuple[float, float, float, float]],
) -> float:
    """Perform full integration steps using vectorized operations.

    Returns the final time t_local.
    """
    from ..alias import get_theta_attr

    nodes = list(G.nodes)
    n_nodes = len(nodes)
    if n_nodes == 0:
        return _validate_clock_grid(t0, dt_step, steps, method)

    # 1. Extract state into arrays
    vf, dnfr, dEPI, epi = (
        np.asarray([state[column] for state in states], dtype=float)
        for column in range(4)
    )

    # The unforced row consumes no phase. Active built-ins validate raw aliases
    # before array coercion, matching the scalar trigonometric-cache owner.
    if _get_gamma_spec(G, strict=True).get("type", "none") == "none":
        theta = np.zeros(n_nodes)
    else:
        theta = np.asarray(
            [
                get_theta_attr(
                    G.nodes[node],
                    0.0,
                    strict=True,
                    conv=lambda raw: _finite_output(raw, "Gamma phase"),
                )
                for node in nodes
            ]
        )

    # Base term: dEPI/dt = vf * dnfr
    # Assumed constant during the step (dnfr doesn't change)
    base = vf * dnfr

    epi_min, epi_max, clip_mode, clip_k = clip_policy

    t_local = t0

    # Pre-allocate d2EPI
    d2EPI = np.zeros_like(dEPI)

    for _ in range(steps):
        epi_previous = epi.copy()
        dEPI_prev = dEPI.copy()

        if method == "rk4":
            # k1
            gamma1 = eval_gamma_vectorized(G, theta, t_local, np, strict=True)
            k1 = base + gamma1

            # k2
            gamma2 = eval_gamma_vectorized(
                G,
                theta,
                t_local + dt_step / INTEGRATORS_HALF_STEP_CANONICAL,
                np,
                strict=True,
            )
            k2 = base + gamma2

            # k3
            gamma3 = gamma2  # Same time point
            k3 = base + gamma3

            # k4
            gamma4 = eval_gamma_vectorized(G, theta, t_local + dt_step, np, strict=True)
            k4 = base + gamma4

            # Update
            epi = _rk4_update(epi, dt_step, (k1, k2, k3, k4))
            dEPI = k4

        else:  # Euler
            gamma = eval_gamma_vectorized(G, theta, t_local, np, strict=True)
            k1 = base + gamma
            epi = euler_update(epi, dt_step, k1)
            dEPI = k1

        # d2EPI
        if dt_step != 0:
            d2EPI = (dEPI - dEPI_prev) / dt_step
        else:
            d2EPI[:] = 0.0

        for values, parameter in (
            (epi, "EPI"),
            (dEPI, "dEPI_dt"),
            (d2EPI, "d2EPI_dt2"),
        ):
            if not np.isfinite(values).all():
                raise NetworkConfigError(
                    parameter=parameter,
                    value=values,
                    reason="Solver output must be finite",
                )

        # Boundary projection must not move a node with zero integrated change.
        changed = epi != epi_previous
        if np.any(changed):
            epi[changed] = structural_clip_array(
                epi[changed], lo=epi_min, hi=epi_max, mode=clip_mode, k=clip_k
            )

        t_local += dt_step

    # Preserve each mapping's authoritative spelling, as the scalar path does.
    # A bulk primary-key write would shadow existing legacy state aliases.
    for node, form, rate, acceleration in zip(nodes, epi, dEPI, d2EPI):
        nd = G.nodes[node]
        set_attr(nd, ALIAS_EPI, form)
        set_attr(nd, ALIAS_DEPI, rate)
        set_attr(nd, ALIAS_D2EPI, acceleration)

    return t_local


class AbstractIntegrator(ABC):
    """Abstract base class encapsulating nodal equation integration."""

    @abstractmethod
    def integrate(
        self,
        graph: TNFRGraph,
        *,
        dt: float | None,
        t: float | None,
        method: str | None,
        n_jobs: int | None,
    ) -> None:
        """Advance ``graph`` coherence states according to the nodal equation."""


class DefaultIntegrator(AbstractIntegrator):
    """Explicit integrator combining Euler and RK4 step implementations."""

    def integrate(
        self,
        graph: TNFRGraph,
        *,
        dt: float | None,
        t: float | None,
        method: str | None,
        n_jobs: int | None,
    ) -> None:
        """Integrate the nodal equation updating EPI, ΔEPI and Δ²EPI."""

        if not isinstance(
            graph, (nx.Graph, nx.DiGraph, nx.MultiGraph, nx.MultiDiGraph)
        ):
            raise TypeError("G must be a networkx graph instance")

        dt_step, steps, t0, resolved_method = prepare_integration_params(
            graph, dt, t, cast(IntegratorMethod | None, method)
        )

        if dt_step == 0.0:
            return

        clip_policy = _integration_clip_policy(graph)
        # Validate all nodal inputs before resolving forcing can populate caches.
        states = [_node_state(nd) for nd in graph.nodes.values()]
        builtin_gamma = _uses_builtin_gamma(_get_gamma_spec(graph, strict=True))
        if np is not None and builtin_gamma:
            t_final = _integrate_vectorized_step(
                graph, dt_step, steps, t0, resolved_method, np, clip_policy, states
            )
            graph.graph["_t"] = t_final
            return

        epi_min, epi_max, clip_mode, clip_k = clip_policy

        owned_keys = ALIAS_EPI + ALIAS_EPI_KIND + ALIAS_DEPI + ALIAS_D2EPI
        t_local = t0
        with _restore_solver_outputs_on_error(graph, owned_keys):
            for _ in range(steps):
                if resolved_method == "rk4":
                    updates = _integrate_rk4(graph, dt_step, t_local, n_jobs=n_jobs)
                else:
                    updates = _integrate_euler(graph, dt_step, t_local, n_jobs=n_jobs)

                # Stage every projected result before committing this substep.
                staged = {}
                for n, (epi, dEPI_dt, d2epi) in updates.items():
                    nd = graph.nodes[n]
                    epi_previous = _read_scalar_epi(nd)
                    epi_clipped = (
                        epi_previous
                        if epi == epi_previous
                        else structural_clip(
                            epi,
                            lo=epi_min,
                            hi=epi_max,
                            mode=clip_mode,
                            k=clip_k,
                            record_stats=False,
                        )
                    )
                    staged[n] = (
                        _finite_output(epi_clipped, "EPI"),
                        _finite_output(dEPI_dt),
                        _finite_output(d2epi),
                        get_attr_str(nd, ALIAS_EPI_KIND, ""),
                    )

                for n, (epi, dEPI_dt, d2epi, epi_kind) in staged.items():
                    nd = graph.nodes[n]
                    set_attr(nd, ALIAS_EPI, epi)
                    if epi_kind:
                        set_attr_str(nd, ALIAS_EPI_KIND, epi_kind)
                    set_attr(nd, ALIAS_DEPI, dEPI_dt)
                    set_attr(nd, ALIAS_D2EPI, d2epi)
                t_local += dt_step

        graph.graph["_t"] = t_local


def update_epi_via_nodal_equation(
    G: TNFRGraph,
    *,
    dt: float | None = None,
    t: float | None = None,
    method: Literal["euler", "rk4"] | None = None,
    n_jobs: int | None = None,
) -> None:
    """TNFR nodal equation with optional extended dynamics.

    Implements either:

    **Classical**: ∂EPI/∂t = νf · ΔNFR(t) + Γi(R)
      - EPI is the node's Primary Information Structure
      - νf is the node's structural frequency (Hz_str)
      - ΔNFR(t) is the nodal gradient (reorganisation need)
      - Γi(R) is optional network coupling via Kuramoto order

    **Extended**: Coupled system with flux fields (when use_extended_dynamics=True)
      - ∂EPI/∂t = νf · ΔNFR(t) + Γi(R) [Same live source registry]
      - ∂θ/∂t = f(νf, ΔNFR, J_φ) [Phase evolution with transport]
      - ∂ΔNFR/∂t = g(∇·J_ΔNFR) [Configured independent pressure response]

    The extended system consumes the shared diagnostic fields J_φ and
    J_ΔNFR through additional configured laws. It does not certify conservation,
    all invariants, or agreement with freshly recomputed canonical pressure.
    On regular undirected support its pressure block has a positive squared-
    Laplacian generator; it is not pressure diffusion. See the constitutive
    audit in theory/DIAGNOSTIC_AND_GRAMMAR_SCOPE.md, section 14.

    Args:
        G: TNFR graph with nodes containing structural attributes
        dt: Integration time step (uses graph default if None)
        t: Current time (uses graph default if None)
        method: Integration method ('euler' or 'rk4')
        n_jobs: Number of parallel jobs for integration

    Notes:
        - This wrapper selects the extension when use_extended_dynamics is True
        - Ordinary runtime dispatches its integrator directly; this flag alone
          does not replace that integrator or its later coordination substep
        - Extended dynamics require J_φ and J_ΔNFR fields (from physics module)
        - Zero flux leaves the forced EPI row and can leave a phase response
        - Extended system preserves backward compatibility (default: False)

    Examples:
        >>> # Classical dynamics (default)
        >>> update_epi_via_nodal_equation(G, dt=0.01)

        >>> # Extended dynamics with flux fields
        >>> G.graph['use_extended_dynamics'] = True
        >>> update_epi_via_nodal_equation(G, dt=0.01)
    """
    # Use the shared configuration parser: the text "false" must not select
    # a different nodal completion merely because it is a nonempty string.
    raw_extended = G.graph.get("use_extended_dynamics", False)
    try:
        use_extended = parse_bool(raw_extended)
    except (TypeError, ValueError) as exc:
        raise NetworkConfigError(
            parameter="use_extended_dynamics",
            value=raw_extended,
            reason=str(exc),
        ) from exc

    if use_extended:
        # Use extended nodal system with flux fields
        _update_extended_nodal_system(G, dt=dt, t=t, method=method, n_jobs=n_jobs)
    else:
        # Use classical TNFR dynamics
        DefaultIntegrator().integrate(
            G,
            dt=dt,
            t=t,
            method=method,
            n_jobs=n_jobs,
        )


def _node_state(nd: dict[str, Any]) -> tuple[float, float, float, float]:
    """Return common node state attributes for canonical equation evaluation.

    Extracts the fundamental TNFR variables from node data:
      - νf (vf): Structural frequency in Hz_str
      - ΔNFR (dnfr): Nodal gradient (reorganization operator)
      - dEPI/dt (previous): Last computed EPI derivative
      - EPI (current): Current Primary Information Structure

    These variables are used in the canonical nodal equation:
        ∂EPI/∂t = νf · ΔNFR(t)

    Args:
        nd: Node data dictionary containing TNFR attributes

    Returns:
        tuple of (vf, dnfr, dEPI_dt_prev, epi_i) with 0.0 defaults

    Notes:
        - vf alias maps to VF, frequency, or structural_frequency
        - dnfr alias maps to DNFR, delta_nfr, or reorganization_gradient
        - EPI must belong to the finite signed scalar chart; rich BEPI is rejected
    """

    from .canonical import validate_nodal_gradient, validate_structural_frequency

    epi_i = _read_scalar_epi(nd)
    vf = get_attr(nd, ALIAS_VF, 0.0, strict=True, conv=validate_structural_frequency)
    dnfr = get_attr(nd, ALIAS_DNFR, 0.0, strict=True, conv=validate_nodal_gradient)
    dEPI_dt_prev = get_attr(nd, ALIAS_DEPI, 0.0, strict=True, conv=_finite_output)
    return vf, dnfr, dEPI_dt_prev, epi_i


def _update_extended_nodal_system(
    G: TNFRGraph,
    *,
    dt: float | None = None,
    t: float | None = None,
    method: Literal["euler", "rk4"] | None = None,
    n_jobs: int | None = None,
) -> None:
    """Advance the coupled EPI/phase/pressure system by synchronous Euler.

    Each substep evaluates canonical fields and all nodal derivatives from the
    same graph state before writing any updates. This optional coupled system
    currently supports Euler only; unsupported methods are rejected explicitly.
    The EPI row reuses the ordinary integrator's live Gamma registry at every
    substep, including at zero capacity. Gamma is a supplied source, not an
    inferred contribution from the diagnostic fluxes.
    Clipping is a boundary policy, not a proof of numerical stability.
    Pressure is independently evolved, not refreshed from EPI/phase/capacity;
    consistency with that constitutive map requires a separate chain-rule test.
    """
    from ..physics.extended import compute_dnfr_flux, compute_phase_current
    from .canonical import compute_extended_nodal_system

    dt_step, steps, t_local, resolved_method = prepare_integration_params(
        G, dt, t, method
    )
    if resolved_method != "euler":
        raise NetworkConfigError(
            parameter="method",
            value=resolved_method,
            reason="Extended nodal dynamics supports only 'euler'",
        )
    if dt_step == 0.0:
        return

    epi_min, epi_max, clip_mode, clip_k = _integration_clip_policy(G)
    # Validate the full nodal input before field readers populate graph caches.
    for nd in G.nodes.values():
        _node_state(nd)
        get_attr(nd, ALIAS_THETA, 0.0, strict=True, conv=_finite_output)

    owned_keys = (
        ALIAS_EPI
        + ALIAS_THETA
        + ALIAS_DNFR
        + ALIAS_DEPI
        + ALIAS_D2EPI
        + ("dtheta_dt", "ddnfr_dt")
    )
    with _restore_solver_outputs_on_error(G, owned_keys):
        for _ in range(steps):
            # Zero flux is valid; it is never a missing-data signal.
            phase_current = compute_phase_current(G)
            pressure_flux = compute_dnfr_flux(G)
            divergences = compute_flux_divergence_vectorized(G, pressure_flux)
            increments = _build_gamma_increments(
                G, dt_step, t_local, method="euler", n_jobs=n_jobs
            )
            updates = {}
            for node in G:
                nd = G.nodes[node]
                vf, dnfr, previous_derivative, epi = _node_state(nd)
                theta = get_attr(nd, ALIAS_THETA, 0.0, strict=True, conv=_finite_output)
                result = compute_extended_nodal_system(
                    nu_f=vf,
                    delta_nfr=dnfr,
                    theta=theta,
                    j_phi=phase_current.get(node, 0.0),
                    j_dnfr_divergence=divergences.get(node, 0.0),
                    coupling_strength=_estimate_local_coupling_strength(G, node),
                    validate_units=False,
                )
                rate = _finite_output(increments[node][0], "dEPI_dt")
                phase_rate = _finite_output(result.phase_derivative, "dtheta_dt")
                pressure_rate = _finite_output(result.dnfr_derivative, "ddnfr_dt")
                d2epi = _finite_output(
                    (rate - previous_derivative) / dt_step, "d2EPI_dt2"
                )
                new_epi = _finite_output(euler_update(epi, dt_step, rate), "EPI")
                if new_epi != epi:
                    new_epi = structural_clip(
                        new_epi,
                        lo=epi_min,
                        hi=epi_max,
                        mode=clip_mode,
                        k=clip_k,
                    )
                new_theta = _finite_output(theta + phase_rate * dt_step, "theta") % (
                    2 * math.pi
                )
                new_dnfr = _finite_output(dnfr + pressure_rate * dt_step, "delta_nfr")
                if new_dnfr != dnfr:
                    new_dnfr = max(
                        -INTEGRATORS_DNFR_BOUNDS_CANONICAL,
                        min(INTEGRATORS_DNFR_BOUNDS_CANONICAL, new_dnfr),
                    )
                updates[node] = (
                    new_epi,
                    new_theta,
                    new_dnfr,
                    rate,
                    d2epi,
                    phase_rate,
                    pressure_rate,
                )

            for node, (
                epi,
                theta,
                dnfr,
                rate,
                d2epi,
                phase_rate,
                pressure_rate,
            ) in updates.items():
                nd = G.nodes[node]
                set_attr(nd, ALIAS_EPI, epi)
                set_attr(nd, ALIAS_THETA, theta)
                set_attr(nd, ALIAS_DNFR, dnfr)
                set_attr(nd, ALIAS_DEPI, rate)
                set_attr(nd, ALIAS_D2EPI, d2epi)
                nd["dtheta_dt"] = phase_rate
                nd["ddnfr_dt"] = pressure_rate
            t_local += dt_step
    G.graph["_t"] = t_local


# Centralized flux divergence computation


def _compute_flux_divergence_centralized(
    G: TNFRGraph, flux_dict: dict[NodeId, float], node: NodeId
) -> float:
    """
    Compute flux divergence using centralized finite difference method.

    This scalar neighbor contrast owns the arithmetic for the reader below.
    It is not an incidence divergence of oriented edge fluxes and carries no
    general conservation guarantee.
    """
    neighbors = list(G.neighbors(node))
    if not neighbors:
        return 0.0

    central_flux = _finite_output(flux_dict.get(node, 0.0), "flux")
    neighbor_fluxes = [
        _finite_output(flux_dict.get(neighbor, 0.0), "flux") for neighbor in neighbors
    ]
    # Retain differences before reduction. Averaging the absolute field first
    # can overflow on a constant field or erase a finite common-offset contrast.
    return mean_neighbor_difference(
        central_flux,
        neighbor_fluxes,
        coefficient=-math.sqrt(len(neighbors)),
    )


def compute_flux_divergence_vectorized(
    G: TNFRGraph, flux_dict: dict[NodeId, float]
) -> dict[NodeId, float]:
    """Evaluate the scalar unique-neighbor divergence for every graph size.

    The existing discretization is sqrt(k_i) * (J_i - mean_neighbor(J)).
    Neighbors are outgoing successors on directed graphs; parallel edges and
    self-loops follow G.neighbors semantics. Edge weights are not metric spacing
    in this diagnostic. Nodes without outgoing neighbors have zero divergence.
    Applied to the shared pressure contrast J=-L_U*p, it gives
    -diag(sqrt(k))*L_U**2*p. The optional pressure response negates this again,
    producing a positive squared-Laplacian generator, not pressure diffusion.
    Active readings use represented-real admission before shared neighbor-
    difference arithmetic. Unrepresentable final outputs reject; a large
    common offset alone neither creates a divergence nor makes it unavailable.
    """
    return {
        node: _compute_flux_divergence_centralized(G, flux_dict, node) for node in G
    }


def _estimate_local_coupling_strength(G: TNFRGraph, node: NodeId) -> float:
    """Evaluate the selected degree-sigmoid transport coefficient."""
    degree = G.degree(node)
    if degree == 0:
        return 0.0

    # Sigmoid coupling: stronger for well-connected nodes
    normalized_degree = min(degree / 10.0, 1.0)  # Saturation at degree 10
    # Configured sensitivity, not a coefficient derived from the nodal law.
    coupling_factor = 4.5
    coupling = 1.0 / (
        1.0
        + math.exp(
            -coupling_factor
            * (normalized_degree - INTEGRATORS_SIGMOID_OFFSET_CANONICAL)
        )
    )

    return coupling