"""ΔNFR structural reorganization pressure utilities and strategies.

This module provides helper functions to configure, cache and apply ΔNFR
components from phase, EPI, structural frequency and topology during
simulations.  The neighbour accumulation helpers reuse cached edge indices
and NumPy workspaces whenever available so cosine, sine, EPI, νf and topology
means remain faithful to the canonical ΔNFR reorganisation without redundant
allocations.
"""

from __future__ import annotations

import inspect
import math
import sys
from collections.abc import Callable, Iterator, Mapping, MutableMapping, Sequence
from concurrent.futures import ProcessPoolExecutor
from fractions import Fraction
from time import perf_counter
from types import ModuleType
from typing import Any, cast

from .._exact_time import finite_represented_real
from ..alias import get_attr, get_theta_attr, set_dnfr
from ..constants import DEFAULTS, get_param
from ..constants.aliases import ALIAS_EPI, ALIAS_VF
from ..mathematics._neighbor_differences import (
    _require_finite_pressure,
    edge_mean_differences,
    mean_neighbor_difference,
)
from ..mathematics._phase_midpoint import certified_two_neighbor_phase
from ..mathematics.unified_numerical import compute_phase_difference, np
from ..metrics.common import (
    _materialize_weight_mapping,
    merge_and_normalize_weights,
    merge_graph_weights,
)
from ..metrics.trig_cache import compute_theta_trig
from ..types import (
    DeltaNFRHook,
    DnfrCacheVectors,
    DnfrVectorMap,
    NeighborStats,
    NodeId,
    TNFRGraph,
    require_finite_real_scalar_epi,
)
from ..utils import (
    DNFR_PREP_STATE_KEY,
    CacheManager,
    DnfrCache,
    DnfrPrepState,
    _graph_cache_manager,
    angle_diff,
    cached_node_list,
    cached_nodes_and_A,
    new_dnfr_cache,
    normalize_weights,
    resolve_chunk_size,
)
from .canonical import validate_structural_frequency
from .fused_dnfr import compute_fused_gradients_symmetric

_SPARSE_DENSITY_THRESHOLD = 0.25
_DNFR_APPROX_BYTES_PER_EDGE = 48
_DNFR_CHANNELS = ("phase", "epi", "vf", "topo")


def _should_vectorize(G: TNFRGraph, np_module: ModuleType | None) -> bool:
    """Return ``True`` when NumPy is available unless the graph disables it."""

    if np_module is None:
        return False
    flag = G.graph.get("vectorized_dnfr")
    if flag is None:
        return True
    return bool(flag)


_NUMPY_CACHE_ATTRS = (
    "theta_np",
    "epi_np",
    "vf_np",
    "cos_theta_np",
    "sin_theta_np",
    "deg_array",
    "neighbor_x_np",
    "neighbor_y_np",
    "neighbor_epi_sum_np",
    "neighbor_vf_sum_np",
    "neighbor_count_np",
    "neighbor_deg_sum_np",
    "neighbor_inv_count_np",
    "neighbor_cos_avg_np",
    "neighbor_sin_avg_np",
    "neighbor_mean_tmp_np",
    "neighbor_mean_length_np",
    "neighbor_accum_np",
    "neighbor_edge_values_np",
    "dense_components_np",
    "dense_accum_np",
    "dense_degree_np",
)


def _profile_start_stop(
    profile: MutableMapping[str, Any] | None,
    *,
    keys: Sequence[str] = (),
) -> tuple[Callable[[], float], Callable[[str, float], None]]:
    """Return helpers to measure wall-clock durations for ``profile`` keys."""

    if profile is not None:
        for key in keys:
            profile.setdefault(key, 0.0)

        def _start() -> float:
            return perf_counter()

        def _stop(metric: str, start: float) -> None:
            profile[metric] = float(profile.get(metric, 0.0)) + (perf_counter() - start)

    else:

        def _start() -> float:
            return 0.0

        def _stop(
            metric: str, start: float
        ) -> None:  # noqa: ARG001 - uniform signature
            return None

    return _start, _stop


def _iter_chunk_offsets(total: int, jobs: int) -> Iterator[tuple[int, int]]:
    """Yield ``(start, end)`` offsets splitting ``total`` items across ``jobs``."""

    if total <= 0 or jobs <= 1:
        return

    jobs = max(1, min(int(jobs), total))
    base, extra = divmod(total, jobs)
    start = 0
    for i in range(jobs):
        size = base + (1 if i < extra else 0)
        if size <= 0:
            continue
        end = start + size
        yield start, end
        start = end


def _neighbor_sums_worker(
    start: int,
    end: int,
    neighbor_indices: Sequence[Sequence[int]],
    cos_th: Sequence[float],
    sin_th: Sequence[float],
    epi: Sequence[float],
    vf: Sequence[float],
    x_base: Sequence[float],
    y_base: Sequence[float],
    epi_base: Sequence[float],
    vf_base: Sequence[float],
    count_base: Sequence[float],
    deg_base: Sequence[float] | None,
    deg_list: Sequence[float] | None,
    degs_list: Sequence[float] | None,
) -> tuple[
    int,
    list[float],
    list[float],
    list[float],
    list[float],
    list[float],
    list[float] | None,
]:
    """Return partial neighbour sums for the ``[start, end)`` range."""

    chunk_x: list[float] = []
    chunk_y: list[float] = []
    chunk_epi: list[float] = []
    chunk_vf: list[float] = []
    chunk_count: list[float] = []
    chunk_deg: list[float] | None = [] if deg_base is not None else None

    for offset, idx in enumerate(range(start, end)):
        neighbors = neighbor_indices[idx]
        x_i = float(x_base[offset])
        y_i = float(y_base[offset])
        epi_i = float(epi_base[offset])
        vf_i = float(vf_base[offset])
        count_i = float(count_base[offset])
        if deg_base is not None and chunk_deg is not None:
            deg_i_acc = float(deg_base[offset])
        else:
            deg_i_acc = 0.0
        deg_i = float(degs_list[idx]) if degs_list is not None else 0.0

        for neighbor_idx in neighbors:
            x_i += float(cos_th[neighbor_idx])
            y_i += float(sin_th[neighbor_idx])
            epi_i += float(epi[neighbor_idx])
            vf_i += float(vf[neighbor_idx])
            count_i += 1.0
            if chunk_deg is not None:
                if deg_list is not None:
                    deg_i_acc += float(deg_list[neighbor_idx])
                else:
                    deg_i_acc += deg_i

        chunk_x.append(x_i)
        chunk_y.append(y_i)
        chunk_epi.append(epi_i)
        chunk_vf.append(vf_i)
        chunk_count.append(count_i)
        if chunk_deg is not None:
            chunk_deg.append(deg_i_acc)

    return (
        start,
        chunk_x,
        chunk_y,
        chunk_epi,
        chunk_vf,
        chunk_count,
        chunk_deg,
    )


def _dnfr_gradients_worker(
    start: int,
    end: int,
    nodes: Sequence[NodeId],
    theta: Sequence[float] | Any,
    epi: Sequence[float] | Any,
    vf: Sequence[float] | Any,
    th_bar: Sequence[float] | Any,
    epi_bar: Sequence[float] | Any,
    vf_bar: Sequence[float] | Any,
    deg_bar: Sequence[float] | Any | None,
    degs: Mapping[Any, float] | Sequence[float] | None,
    w_phase: float,
    w_epi: float,
    w_vf: float,
    w_topo: float,
    epi_gradient: Sequence[float] | None = None,
    vf_gradient: Sequence[float] | None = None,
    phase_overrides: Mapping[int, float] | None = None,
) -> tuple[int, list[float]]:
    """Return partial ΔNFR gradients for the ``[start, end)`` range."""

    chunk: list[float] = []
    for idx in range(start, end):
        n = nodes[idx]
        g_phase = (
            phase_overrides[idx]
            if phase_overrides and idx in phase_overrides
            else -angle_diff(theta[idx], th_bar[idx]) / math.pi
        )
        g_epi = (
            (epi_gradient[idx] if epi_gradient is not None else epi_bar[idx] - epi[idx])
            if w_epi
            else 0.0
        )
        g_vf = (
            (vf_gradient[idx] if vf_gradient is not None else vf_bar[idx] - vf[idx])
            if w_vf
            else 0.0
        )
        if w_topo != 0.0 and deg_bar is not None and degs is not None:
            if isinstance(degs, dict):
                deg_i = float(degs.get(n, 0))
            else:
                deg_i = float(degs[idx])
            g_topo = deg_bar[idx] - deg_i
        else:
            g_topo = 0.0
        chunk.append(
            w_phase * g_phase
            + (g_epi if epi_gradient is not None else w_epi * g_epi)
            + (g_vf if vf_gradient is not None else w_vf * g_vf)
            + w_topo * g_topo
        )
    return start, chunk


def _resolve_parallel_jobs(n_jobs: int | None, total: int) -> int | None:
    """Return an effective worker count for ``total`` items or ``None``."""

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


def _is_numpy_like(obj) -> bool:
    return (
        getattr(obj, "dtype", None) is not None
        and getattr(obj, "shape", None) is not None
    )


def _has_cached_numpy_buffers(
    data: MutableMapping[str, Any], cache: DnfrCache | None
) -> bool:
    for attr in _NUMPY_CACHE_ATTRS:
        arr = data.get(attr)
        if _is_numpy_like(arr):
            return True
    if cache is not None:
        for attr in _NUMPY_CACHE_ATTRS:
            arr = getattr(cache, attr, None)
            if _is_numpy_like(arr):
                return True
    A = data.get("A")
    if _is_numpy_like(A):
        return True
    return False


__all__ = (
    "default_compute_delta_nfr",
    "set_delta_nfr_hook",
    "dnfr_phase_only",
    "dnfr_epi_vf_mixed",
    "dnfr_laplacian",
    "compute_delta_nfr_hamiltonian",
)


def _write_dnfr_metadata(
    G, *, weights: Mapping[str, float], hook_name: str, note: str | None = None
) -> None:
    """Write a ``_DNFR_META`` block in ``G.graph`` with the mix and hook name.

    ``weights`` may include arbitrary components (phase/epi/vf/topo/etc.).
    """
    weights_norm = normalize_weights(weights, weights.keys())
    meta = {
        "hook": hook_name,
        "weights_raw": dict(weights),
        "weights_effective": dict(weights),
        "weights_norm": weights_norm,
        "components": [k for k, v in weights.items() if v != 0.0],
        "doc": "ΔNFR = Σ w_i·g_i",
    }
    if note:
        meta["note"] = str(note)
    G.graph["_DNFR_META"] = meta
    G.graph["_dnfr_hook_name"] = hook_name  # string friendly


def _resolve_dnfr_weights(G) -> Mapping[str, float]:
    """Read effective coefficients without mutating the graph.

    A legacy explicitly supplied private mix remains an override until the
    public configuration changes after its first pressure preparation. An
    engine-owned cached mix must never hide a subsequent public configuration
    edit. Detached observers use this same resolution rule.
    """
    cached = G.graph.get("_dnfr_weights")
    source = G.graph.get("_dnfr_weights_source")
    configured = _materialize_weight_mapping(
        merge_graph_weights(G, "DNFR_WEIGHTS", strict=True),
        _DNFR_CHANNELS,
        name="DNFR_WEIGHTS",
    )
    if cached is not None and (source is None or source == configured):
        return _materialize_weight_mapping(
            cached, _DNFR_CHANNELS, name="cached DNFR weights"
        )
    return merge_and_normalize_weights(G, "DNFR_WEIGHTS", _DNFR_CHANNELS, default=0.0)


def _configure_dnfr_weights(G) -> Mapping[str, float]:
    """Store the effective mix and a detached snapshot of its public source."""
    weights = _resolve_dnfr_weights(G)
    G.graph["_dnfr_weights"] = weights
    G.graph["_dnfr_weights_source"] = _materialize_weight_mapping(
        merge_graph_weights(G, "DNFR_WEIGHTS", strict=True),
        _DNFR_CHANNELS,
        name="DNFR_WEIGHTS",
    )
    return weights


def _read_scalar_epi(nd: Mapping[str, Any]) -> float:
    """Read the finite signed scalar chart required by this pressure model."""
    return get_attr(
        nd, ALIAS_EPI, 0.0, strict=True, conv=require_finite_real_scalar_epi
    )


def _read_capacity(nd: Mapping[str, Any], default: float = 0.0) -> float:
    """Apply nodal capacity admission before any alias or float coercion."""
    return get_attr(
        nd, ALIAS_VF, default, strict=True, conv=validate_structural_frequency
    )


def _read_phase(nd: Mapping[str, Any]) -> float:
    """Read the authoritative finite phase before building trigonometric data."""
    return get_theta_attr(
        nd,
        0.0,
        strict=True,
        conv=lambda raw: finite_represented_real(raw, "pressure phase")[0],
    )


def _pressure_phase_coordinate(value: float) -> float:
    """Read the principal phasor chart before a phase-pressure subtraction."""
    return (
        value if abs(value) <= math.pi else math.atan2(math.sin(value), math.cos(value))
    )


def _init_dnfr_cache(
    G: TNFRGraph,
    nodes: Sequence[NodeId],
    cache_or_manager: CacheManager | DnfrCache | None = None,
    checksum: Any | None = None,
    force_refresh: bool = False,
    *,
    manager: CacheManager | None = None,
) -> tuple[
    DnfrCache,
    dict[NodeId, int],
    list[float],
    list[float],
    list[float],
    list[float],
    list[float],
    bool,
]:
    """Initialise or reuse cached ΔNFR arrays.

    ``manager`` telemetry records cache hits, misses and timings. Legacy
    callers can still pass a ``cache`` instance as the
    third positional argument; this helper supports both signatures by seeding
    the manager-backed state with the provided cache when necessary.
    """

    if manager is None and isinstance(cache_or_manager, CacheManager):
        manager = cache_or_manager
        cache_or_manager = None

    if manager is None:
        manager = _graph_cache_manager(G.graph)

    graph = G.graph
    state = manager.get(DNFR_PREP_STATE_KEY)
    if not isinstance(state, DnfrPrepState):
        manager.clear(DNFR_PREP_STATE_KEY)
        state = manager.get(DNFR_PREP_STATE_KEY)

    if isinstance(cache_or_manager, DnfrCache):
        state.cache = cache_or_manager
        if checksum is None:
            checksum = cache_or_manager.checksum

    cache = state.cache
    reuse = (
        not force_refresh
        and isinstance(cache, DnfrCache)
        and cache.checksum == checksum
        and len(cache.theta) == len(nodes)
    )
    if reuse:
        manager.increment_hit(DNFR_PREP_STATE_KEY)
        graph["_dnfr_prep_cache"] = cache
        return (
            cache,
            cache.idx,
            cache.theta,
            cache.epi,
            cache.vf,
            cache.cos_theta,
            cache.sin_theta,
            False,
        )

    def _rebuild(current: DnfrPrepState | Any) -> DnfrPrepState:
        if not isinstance(current, DnfrPrepState):
            raise RuntimeError("ΔNFR prep state unavailable during rebuild")
        prev_cache = current.cache if isinstance(current.cache, DnfrCache) else None
        idx_local = {n: i for i, n in enumerate(nodes)}
        size = len(nodes)
        zeros = [0.0] * size
        cache_new = prev_cache if prev_cache is not None else new_dnfr_cache()
        cache_new.idx = idx_local
        cache_new.theta = zeros.copy()
        cache_new.epi = zeros.copy()
        cache_new.vf = zeros.copy()
        cache_new.cos_theta = [1.0] * size
        cache_new.sin_theta = [0.0] * size
        cache_new.neighbor_x = zeros.copy()
        cache_new.neighbor_y = zeros.copy()
        cache_new.neighbor_epi_sum = zeros.copy()
        cache_new.neighbor_vf_sum = zeros.copy()
        cache_new.neighbor_count = zeros.copy()
        cache_new.neighbor_deg_sum = zeros.copy() if size else []
        cache_new.degs = None
        cache_new.edge_src = None
        cache_new.edge_dst = None
        cache_new.checksum = checksum

        # Reset any numpy mirrors or aggregated buffers to avoid leaking
        # state across refresh cycles (e.g. switching between vectorised
        # and Python paths or reusing legacy caches).
        if prev_cache is not None:
            for attr in _NUMPY_CACHE_ATTRS:
                setattr(cache_new, attr, None)
            for attr in (
                "th_bar_np",
                "epi_bar_np",
                "vf_bar_np",
                "deg_bar_np",
                "grad_phase_np",
                "grad_epi_np",
                "grad_vf_np",
                "grad_topo_np",
                "grad_total_np",
            ):
                setattr(cache_new, attr, None)
            cache_new.edge_src = None
            cache_new.edge_dst = None
            cache_new.edge_signature = None
            cache_new.neighbor_accum_signature = None
        cache_new.degs = prev_cache.degs if prev_cache else None
        cache_new.checksum = checksum
        current.cache = cache_new
        graph["_dnfr_prep_cache"] = cache_new
        return current

    with manager.timer(DNFR_PREP_STATE_KEY):
        state = manager.update(DNFR_PREP_STATE_KEY, _rebuild)
    manager.increment_miss(DNFR_PREP_STATE_KEY)
    cache = state.cache
    if not isinstance(cache, DnfrCache):  # pragma: no cover - defensive guard
        raise RuntimeError("ΔNFR cache initialisation failed")
    return (
        cache,
        cache.idx,
        cache.theta,
        cache.epi,
        cache.vf,
        cache.cos_theta,
        cache.sin_theta,
        True,
    )


def _ensure_numpy_vectors(cache: DnfrCache) -> DnfrCacheVectors:
    """Ensure NumPy copies of cached vectors are initialised and up to date."""

    if cache is None:
        return (None, None, None, None, None)

    arrays: list[Any | None] = []
    size = len(cache.theta)
    for attr_np, source_attr in (
        ("theta_np", "theta"),
        ("epi_np", "epi"),
        ("vf_np", "vf"),
        ("cos_theta_np", "cos_theta"),
        ("sin_theta_np", "sin_theta"),
    ):
        arr = getattr(cache, attr_np)
        if arr is not None and getattr(arr, "shape", None) == (size,):
            arrays.append(arr)
            continue
        src = getattr(cache, source_attr)
        if src is None:
            setattr(cache, attr_np, None)
            arrays.append(None)
            continue
        arr = np.asarray(src, dtype=float)
        if getattr(arr, "shape", None) != (size,):
            arr = np.array(src, dtype=float)
        setattr(cache, attr_np, arr)
        arrays.append(arr)
    return tuple(arrays)


def _ensure_numpy_degrees(
    cache: DnfrCache,
    deg_list: Sequence[float] | None,
) -> np.ndarray | None:
    """Initialise/update NumPy array mirroring ``deg_list``.

    Deg_array reuse pattern:
    -------------------------
    The degree array (deg_array) is a cached NumPy buffer that stores node
    degrees for topology-based ΔNFR computations. The reuse pattern follows:

    1. **Allocation**: Created once when topology weight (w_topo) > 0 or when
       caching is enabled, sized to match the node count.

    2. **Reuse across steps**: When the graph topology is stable (no edge
       additions/removals), the same deg_array buffer is reused across
       multiple ΔNFR computation steps by updating in-place via np.copyto.

    3. **Count buffer optimization**: For undirected graphs where node degree
       equals neighbor count, deg_array can serve double duty as the count
       buffer (see _accumulate_neighbors_numpy lines 2185-2194), eliminating
       the need for an extra accumulator row.

    4. **Invalidation**: Cache is cleared when graph.edges changes or when
       _dnfr_prep_dirty flag is set, ensuring fresh allocation on next use.

    This pattern maintains ΔNFR computational accuracy (Invariant #1:
    Nodal Equation Integrity) while
    minimizing allocations for stable topologies.
    """

    if deg_list is None:
        if cache is not None:
            cache.deg_array = None
        return None
    if cache is None:
        return np.array(deg_list, dtype=float)
    arr = cache.deg_array
    if arr is None or len(arr) != len(deg_list):
        arr = np.array(deg_list, dtype=float)
    else:
        np.copyto(arr, deg_list, casting="unsafe")
    cache.deg_array = arr
    return arr


def _resolve_numpy_degree_array(
    data: MutableMapping[str, Any],
    count: np.ndarray | None,
    *,
    cache: DnfrCache | None,
) -> np.ndarray | None:
    """Return the vector of node degrees required for topology gradients."""

    if data["w_topo"] == 0.0:
        return None
    deg_array = data.get("deg_array")
    if deg_array is not None:
        return deg_array
    deg_list = data.get("deg_list")
    if deg_list is not None:
        deg_array = np.array(deg_list, dtype=float)
        data["deg_array"] = deg_array
        if cache is not None:
            cache.deg_array = deg_array
        return deg_array
    return count


def _ensure_cached_array(
    cache: DnfrCache | None,
    attr: str,
    shape: tuple[int, ...],
) -> np.ndarray:
    """Return a cached NumPy buffer with ``shape`` creating/reusing it."""

    if np is None:
        raise RuntimeError("NumPy is required to build cached arrays")
    arr = getattr(cache, attr) if cache is not None else None
    if arr is None or getattr(arr, "shape", None) != shape:
        arr = np.empty(shape, dtype=float)
        if cache is not None:
            setattr(cache, attr, arr)
    return arr


def _ensure_numpy_state_vectors(data: MutableMapping[str, Any]) -> DnfrVectorMap:
    """Synchronise list-based state vectors with their NumPy counterparts."""

    nodes = data.get("nodes") or ()
    size = len(nodes)
    cache: DnfrCache | None = data.get("cache")

    cache_arrays: DnfrCacheVectors = (None, None, None, None, None)
    if cache is not None:
        cache_arrays = _ensure_numpy_vectors(cache)

    result: dict[str, Any | None] = {}
    for plain_key, np_key, cached_arr, result_key in (
        ("theta", "theta_np", cache_arrays[0], "theta"),
        ("epi", "epi_np", cache_arrays[1], "epi"),
        ("vf", "vf_np", cache_arrays[2], "vf"),
        ("cos_theta", "cos_theta_np", cache_arrays[3], "cos"),
        ("sin_theta", "sin_theta_np", cache_arrays[4], "sin"),
    ):
        arr = data.get(np_key)
        if arr is None:
            arr = cached_arr
        if arr is None or getattr(arr, "shape", None) != (size,):
            src = data.get(plain_key)
            if src is None and cache is not None:
                src = getattr(cache, plain_key)
            if src is None:
                arr = None
            else:
                arr = np.asarray(src, dtype=float)
                if getattr(arr, "shape", None) != (size,):
                    arr = np.array(src, dtype=float)
        if arr is not None:
            data[np_key] = arr
            data[plain_key] = arr
            if cache is not None:
                setattr(cache, np_key, arr)
        else:
            data[np_key] = None
        result[result_key] = arr

    return result


def _build_edge_index_arrays(
    G: TNFRGraph,
    nodes: Sequence[NodeId],
    idx: Mapping[NodeId, int],
) -> tuple[Any, Any]:
    """Create (src, dst) index arrays for ``G`` respecting ``nodes`` order."""
    from ..mathematics.unified_numerical import np

    if np is None:
        return None, None
    if not nodes:
        empty = np.empty(0, dtype=np.intp)
        return empty, empty

    src = []
    dst = []
    append_src = src.append
    append_dst = dst.append
    for node in nodes:
        i = idx.get(node)
        if i is None:
            continue
        for neighbor in G.neighbors(node):
            j = idx.get(neighbor)
            if j is None:
                continue
            append_src(i)
            append_dst(j)
    if not src:
        empty = np.empty(0, dtype=np.intp)
        return empty, empty
    edge_src = np.asarray(src, dtype=np.intp)
    edge_dst = np.asarray(dst, dtype=np.intp)
    return edge_src, edge_dst


def _build_edge_weight_array(
    G: TNFRGraph,
    nodes: Sequence[NodeId],
    edge_src: Any,
    edge_dst: Any,
) -> Any:
    """Per-edge weights aligned to ``(edge_src, edge_dst)``, read live from G.

    The EPI channel realizes the random-walk Laplacian L_rw = I − D⁻¹W, so it
    must use the current edge weights.  Reading them fresh each call keeps the
    weighted diffusion correct even when edge weights mutate without a topology
    change.  Returns ``None`` when every weight is unity — the fused kernel then
    takes the exact unweighted path.
    """
    if np is None or edge_src is None or edge_dst is None:
        return None
    n_edges = int(getattr(edge_src, "shape", (0,))[0])
    if n_edges == 0:
        return None
    weights = None
    is_multigraph = G.is_multigraph()
    for k in range(n_edges):
        u = nodes[int(edge_src[k])]
        v = nodes[int(edge_dst[k])]
        edge_data = G[u][v]
        if is_multigraph:
            # NetworkX adjacency sums parallel edges; the EPI channel must
            # realize the same W while other channels count unique neighbors.
            w = float(sum(attrs.get("weight", 1.0) for attrs in edge_data.values()))
        else:
            w = float(edge_data.get("weight", 1.0))
        if w != 1.0:
            if weights is None:
                weights = np.ones(n_edges, dtype=float)
            weights[k] = w
    return weights


def _refresh_dnfr_vectors(
    G: TNFRGraph, nodes: Sequence[NodeId], cache: DnfrCache
) -> None:
    """Update cached angle and state vectors for ΔNFR."""
    np_module = np
    trig = compute_theta_trig((n, _read_phase(G.nodes[n])) for n in nodes)
    use_numpy = _should_vectorize(G, np_module)
    node_count = len(nodes)
    trig_theta = getattr(trig, "theta_values", None)
    trig_cos = getattr(trig, "cos_values", None)
    trig_sin = getattr(trig, "sin_values", None)
    np_ready = (
        use_numpy
        and np_module is not None
        and isinstance(trig_theta, getattr(np_module, "ndarray", tuple()))
        and isinstance(trig_cos, getattr(np_module, "ndarray", tuple()))
        and isinstance(trig_sin, getattr(np_module, "ndarray", tuple()))
        and getattr(trig_theta, "shape", None) == getattr(trig_cos, "shape", None)
        and getattr(trig_theta, "shape", None) == getattr(trig_sin, "shape", None)
        and (trig_theta.shape[0] if getattr(trig_theta, "ndim", 0) else 0) == node_count
    )

    if np_ready:
        if node_count:
            epi_arr = np_module.fromiter(
                (_read_scalar_epi(G.nodes[node]) for node in nodes),
                dtype=float,
                count=node_count,
            )
            vf_arr = np_module.fromiter(
                (_read_capacity(G.nodes[node]) for node in nodes),
                dtype=float,
                count=node_count,
            )
        else:
            epi_arr = np_module.empty(0, dtype=float)
            vf_arr = np_module.empty(0, dtype=float)

        theta_arr = np_module.asarray(trig_theta, dtype=float)
        cos_arr = np_module.asarray(trig_cos, dtype=float)
        sin_arr = np_module.asarray(trig_sin, dtype=float)

        def _sync_numpy(attr: str, source: Any) -> Any:
            dest = getattr(cache, attr)
            if dest is None or getattr(dest, "shape", None) != source.shape:
                dest = np_module.array(source, dtype=float)
            else:
                np_module.copyto(dest, source, casting="unsafe")
            setattr(cache, attr, dest)
            return dest

        _sync_numpy("theta_np", theta_arr)
        _sync_numpy("epi_np", epi_arr)
        _sync_numpy("vf_np", vf_arr)
        _sync_numpy("cos_theta_np", cos_arr)
        _sync_numpy("sin_theta_np", sin_arr)

        # Python mirrors remain untouched while the vectorised path is active.
        # They will be rebuilt the next time the runtime falls back to lists.
        if cache.theta is not None and len(cache.theta) != node_count:
            cache.theta = [0.0] * node_count
        if cache.epi is not None and len(cache.epi) != node_count:
            cache.epi = [0.0] * node_count
        if cache.vf is not None and len(cache.vf) != node_count:
            cache.vf = [0.0] * node_count
        if cache.cos_theta is not None and len(cache.cos_theta) != node_count:
            cache.cos_theta = [1.0] * node_count
        if cache.sin_theta is not None and len(cache.sin_theta) != node_count:
            cache.sin_theta = [0.0] * node_count
    else:
        for index, node in enumerate(nodes):
            i: int = int(index)
            node_id: NodeId = node
            nd = G.nodes[node_id]
            cache.theta[i] = trig.theta[node_id]
            cache.epi[i] = _read_scalar_epi(nd)
            cache.vf[i] = _read_capacity(nd)
            cache.cos_theta[i] = trig.cos[node_id]
            cache.sin_theta[i] = trig.sin[node_id]
        if use_numpy and np_module is not None:
            _ensure_numpy_vectors(cache)
        else:
            cache.theta_np = None
            cache.epi_np = None
            cache.vf_np = None
            cache.cos_theta_np = None
            cache.sin_theta_np = None


def _prepare_dnfr_data(
    G: TNFRGraph,
    *,
    cache_size: int | None = 128,
    profile: MutableMapping[str, Any] | None = None,
) -> dict[str, Any]:
    """Precompute common data for ΔNFR strategies.

    The helper decides between edge-wise and dense adjacency accumulation
    heuristically.  Graphs whose edge density exceeds
    ``_SPARSE_DENSITY_THRESHOLD`` receive a cached adjacency matrix so the
    dense path can be exercised; callers may also force the dense mode by
    setting ``G.graph['dnfr_force_dense']`` to a truthy value.

    Parameters
    ----------
    profile : MutableMapping[str, float] or None, optional
        Mutable mapping that accumulates wall-clock timings for ΔNFR
        preparation. When provided the helper increases the
        ``"dnfr_cache_rebuild"`` bucket with the time spent refreshing cached
        node vectors and associated NumPy workspaces.
    """
    start_timer, stop_timer = _profile_start_stop(
        profile,
        keys=("dnfr_cache_rebuild",),
    )

    graph = G.graph
    weights = _configure_dnfr_weights(G)

    result: dict[str, Any] = {
        "weights": weights,
        "cache_size": cache_size,
    }

    np_module = np
    use_numpy = _should_vectorize(G, np_module)

    nodes = cast(tuple[NodeId, ...], cached_node_list(G))
    edge_count = G.number_of_edges()

    # Centralized decision logic for sparse vs dense accumulation path.
    # This decision affects which accumulation strategy will be used:
    #   - "sparse": edge-based accumulation (_accumulate_neighbors_broadcasted)
    #   - "dense": matrix multiplication with adjacency matrix (_accumulate_neighbors_dense)
    # The decision is stored in dnfr_path_decision for telemetry and debugging.
    prefer_sparse = False
    dense_override = bool(G.graph.get("dnfr_force_dense"))
    dnfr_path_decision = "fallback"  # Default when numpy unavailable

    if use_numpy:
        # Heuristic: use sparse path when density <= _SPARSE_DENSITY_THRESHOLD (0.25)
        prefer_sparse = _prefer_sparse_accumulation(len(nodes), edge_count)

        if dense_override:
            # User explicitly requested dense mode
            prefer_sparse = False
            dnfr_path_decision = "dense_forced"
        elif not prefer_sparse:
            # Heuristic chose dense path (high density graph)
            dnfr_path_decision = "dense_auto"
        else:
            # Heuristic chose sparse path (low density graph)
            dnfr_path_decision = "sparse"

    nodes_cached, A_untyped = cached_nodes_and_A(
        G,
        cache_size=cache_size,
        require_numpy=False,
        prefer_sparse=prefer_sparse,
        nodes=nodes,
    )
    nodes = cast(tuple[NodeId, ...], nodes_cached)
    A: np.ndarray | None = A_untyped
    result["nodes"] = nodes
    result["A"] = A
    manager = _graph_cache_manager(G.graph)
    checksum = G.graph.get("_dnfr_nodes_checksum")
    dirty_flag = bool(G.graph.pop("_dnfr_prep_dirty", False))
    existing_cache = cast(DnfrCache | None, graph.get("_dnfr_prep_cache"))
    cache_timer = start_timer()
    cache, idx, theta, epi, vf, cos_theta, sin_theta, refreshed = _init_dnfr_cache(
        G,
        nodes,
        existing_cache,
        checksum,
        force_refresh=dirty_flag,
        manager=manager,
    )
    stop_timer("dnfr_cache_rebuild", cache_timer)
    dirty = dirty_flag or refreshed
    caching_enabled = cache is not None and (cache_size is None or cache_size > 0)
    result["cache"] = cache
    result["idx"] = idx
    result["theta"] = theta
    result["epi"] = epi
    result["vf"] = vf
    result["cos_theta"] = cos_theta
    result["sin_theta"] = sin_theta
    if cache is not None:
        _refresh_dnfr_vectors(G, nodes, cache)
        if np_module is None and not caching_enabled:
            for attr in (
                "neighbor_x_np",
                "neighbor_y_np",
                "neighbor_epi_sum_np",
                "neighbor_vf_sum_np",
                "neighbor_count_np",
                "neighbor_deg_sum_np",
                "neighbor_inv_count_np",
                "neighbor_cos_avg_np",
                "neighbor_sin_avg_np",
                "neighbor_mean_tmp_np",
                "neighbor_mean_length_np",
                "neighbor_accum_np",
                "neighbor_edge_values_np",
            ):
                setattr(cache, attr, None)
            cache.neighbor_accum_signature = None
            for attr in (
                "th_bar_np",
                "epi_bar_np",
                "vf_bar_np",
                "deg_bar_np",
                "grad_phase_np",
                "grad_epi_np",
                "grad_vf_np",
                "grad_topo_np",
                "grad_total_np",
            ):
                setattr(cache, attr, None)

    w_phase = float(weights.get("phase", 0.0))
    w_epi = float(weights.get("epi", 0.0))
    w_vf = float(weights.get("vf", 0.0))
    w_topo = float(weights.get("topo", 0.0))
    result["w_phase"] = w_phase
    result["w_epi"] = w_epi
    result["w_vf"] = w_vf
    result["w_topo"] = w_topo
    degree_map = cast(dict[NodeId, float] | None, cache.degs if cache else None)
    if cache is not None and dirty:
        cache.degs = None
        cache.deg_list = None
        cache.deg_array = None
        cache.edge_src = None
        cache.edge_dst = None
        cache.edge_signature = None
        cache.neighbor_accum_signature = None
        cache.neighbor_accum_np = None
        cache.neighbor_edge_values_np = None
        degree_map = None

    deg_list: list[float] | None = None
    degs: dict[NodeId, float] | None = None
    deg_array: np.ndarray | None = None

    if w_topo != 0.0 or caching_enabled:
        if degree_map is None or len(degree_map) != len(G):
            # Every channel uses unique outgoing neighbors. Total degree
            # includes incoming arcs and doubles loops, unlike the fused
            # kernel's counts and the scalar G.neighbors contract.
            degree_map = {
                cast(NodeId, node): float(sum(1 for _ in G.neighbors(node)))
                for node in nodes
            }
            if cache is not None:
                cache.degs = degree_map

        if (
            cache is not None
            and cache.deg_list is not None
            and not dirty
            and len(cache.deg_list) == len(nodes)
        ):
            deg_list = cache.deg_list
        else:
            deg_list = [float(degree_map.get(node, 0.0)) for node in nodes]
            if cache is not None:
                cache.deg_list = deg_list

        degs = degree_map

        if np_module is not None and deg_list is not None:
            if cache is not None:
                deg_array = _ensure_numpy_degrees(cache, deg_list)
            else:
                deg_array = np_module.array(deg_list, dtype=float)
        elif cache is not None:
            cache.deg_array = None
    elif cache is not None and dirty:
        cache.deg_list = None
        cache.deg_array = None

    G.graph["_dnfr_prep_dirty"] = False

    result["degs"] = degs
    result["deg_list"] = deg_list

    theta_np: np.ndarray | None
    epi_np: np.ndarray | None
    vf_np: np.ndarray | None
    cos_theta_np: np.ndarray | None
    sin_theta_np: np.ndarray | None
    edge_src: np.ndarray | None
    edge_dst: np.ndarray | None
    if use_numpy:
        theta_np, epi_np, vf_np, cos_theta_np, sin_theta_np = _ensure_numpy_vectors(
            cache
        )
        edge_src = None
        edge_dst = None
        if cache is not None:
            edge_src = cache.edge_src
            edge_dst = cache.edge_dst
            if edge_src is None or edge_dst is None or dirty:
                edge_src, edge_dst = _build_edge_index_arrays(G, nodes, idx)
                cache.edge_src = edge_src
                cache.edge_dst = edge_dst
        else:
            edge_src, edge_dst = _build_edge_index_arrays(G, nodes, idx)

        if cache is not None:
            for attr in ("neighbor_accum_np", "neighbor_edge_values_np"):
                arr = getattr(cache, attr, None)
                if arr is not None:
                    result[attr] = arr
        if edge_src is not None and edge_dst is not None:
            signature = (id(edge_src), id(edge_dst), len(nodes))
            result["edge_signature"] = signature
            if cache is not None:
                cache.edge_signature = signature
    else:
        theta_np = None
        epi_np = None
        vf_np = None
        cos_theta_np = None
        sin_theta_np = None
        edge_src = None
        edge_dst = None
        if cache is not None:
            cache.edge_src = None
            cache.edge_dst = None

    result.setdefault("neighbor_edge_values_np", None)
    if cache is not None and "edge_signature" not in result:
        result["edge_signature"] = cache.edge_signature

    result["theta_np"] = theta_np
    result["epi_np"] = epi_np
    result["vf_np"] = vf_np
    result["cos_theta_np"] = cos_theta_np
    result["sin_theta_np"] = sin_theta_np
    if theta_np is not None and getattr(theta_np, "shape", None) == (len(nodes),):
        result["theta"] = theta_np
    if epi_np is not None and getattr(epi_np, "shape", None) == (len(nodes),):
        result["epi"] = epi_np
    if vf_np is not None and getattr(vf_np, "shape", None) == (len(nodes),):
        result["vf"] = vf_np
    if cos_theta_np is not None and getattr(cos_theta_np, "shape", None) == (
        len(nodes),
    ):
        result["cos_theta"] = cos_theta_np
    if sin_theta_np is not None and getattr(sin_theta_np, "shape", None) == (
        len(nodes),
    ):
        result["sin_theta"] = sin_theta_np
    result["deg_array"] = deg_array
    result["edge_src"] = edge_src
    result["edge_dst"] = edge_dst
    result["edge_count"] = edge_count
    result["prefer_sparse"] = prefer_sparse
    result["dense_override"] = dense_override
    result["dnfr_path_decision"] = dnfr_path_decision
    result["stable_linear_gradients"] = True
    result.setdefault("neighbor_accum_np", None)
    result.setdefault("neighbor_accum_signature", None)

    return result


def _apply_dnfr_gradients(
    G: TNFRGraph,
    data: MutableMapping[str, Any],
    th_bar: Sequence[float] | np.ndarray,
    epi_bar: Sequence[float] | np.ndarray,
    vf_bar: Sequence[float] | np.ndarray,
    deg_bar: Sequence[float] | np.ndarray | None = None,
    degs: Mapping[Any, float] | Sequence[float] | np.ndarray | None = None,
    *,
    n_jobs: int | None = None,
    profile: MutableMapping[str, Any] | None = None,
) -> None:
    """Combine precomputed gradients and write ΔNFR to each node.

    Parameters
    ----------
    profile : MutableMapping[str, float] or None, optional
        Mutable mapping receiving aggregated timings for the gradient assembly
        (``"dnfr_gradient_assembly"``) and in-place writes
        (``"dnfr_inplace_write"``).
    """
    start_timer, stop_timer = _profile_start_stop(
        profile,
        keys=("dnfr_gradient_assembly", "dnfr_inplace_write"),
    )

    nodes = data["nodes"]
    theta = data["theta"]
    epi = data["epi"]
    vf = data["vf"]
    w_phase = data["w_phase"]
    w_epi = data["w_epi"]
    w_vf = data["w_vf"]
    w_topo = data["w_topo"]
    epi_gradient = data.get("epi_gradient")
    vf_gradient = data.get("vf_gradient")
    if degs is None:
        degs = data.get("degs")

    cache: DnfrCache | None = data.get("cache")

    theta_np = data.get("theta_np")
    epi_np = data.get("epi_np")
    vf_np = data.get("vf_np")
    deg_array = data.get("deg_array") if w_topo != 0.0 else None

    use_vector = (
        np is not None
        and theta_np is not None
        and epi_np is not None
        and vf_np is not None
        and isinstance(th_bar, np.ndarray)
        and isinstance(epi_bar, np.ndarray)
        and isinstance(vf_bar, np.ndarray)
    )
    if use_vector and w_topo != 0.0:
        use_vector = (
            deg_bar is not None
            and isinstance(deg_bar, np.ndarray)
            and isinstance(deg_array, np.ndarray)
        )

    grad_timer = start_timer()
    phase_overrides = {}
    if w_phase != 0.0:
        indices = data["idx"]
        for index, node in enumerate(nodes):
            neighbors = tuple(G.neighbors(node))
            if len(neighbors) == 2:
                certified = certified_two_neighbor_phase(
                    float(theta[index]),
                    float(theta[indices[neighbors[0]]]),
                    float(theta[indices[neighbors[1]]]),
                )
                if certified is not None:
                    phase_overrides[index] = certified.delta / math.pi

    if use_vector:
        grad_phase = _ensure_cached_array(cache, "grad_phase_np", theta_np.shape)
        grad_epi = _ensure_cached_array(cache, "grad_epi_np", epi_np.shape)
        grad_vf = _ensure_cached_array(cache, "grad_vf_np", vf_np.shape)
        grad_total = _ensure_cached_array(cache, "grad_total_np", theta_np.shape)
        grad_topo = None
        if w_topo != 0.0:
            grad_topo = _ensure_cached_array(cache, "grad_topo_np", deg_array.shape)

        center_phase = np.where(
            np.abs(theta_np) <= math.pi,
            theta_np,
            np.arctan2(data["sin_theta"], data["cos_theta"]),
        )
        mean_phase = np.where(
            np.abs(th_bar) <= math.pi,
            th_bar,
            np.arctan2(np.sin(th_bar), np.cos(th_bar)),
        )
        np.copyto(grad_phase, compute_phase_difference(mean_phase, center_phase))
        np.divide(grad_phase, math.pi, out=grad_phase)
        for index, value in phase_overrides.items():
            grad_phase[index] = value

        if w_epi == 0.0:
            grad_epi.fill(0.0)
        elif epi_gradient is not None:
            np.copyto(grad_epi, epi_gradient, casting="unsafe")
        else:
            np.subtract(epi_bar, epi_np, out=grad_epi)

        if w_vf == 0.0:
            grad_vf.fill(0.0)
        elif vf_gradient is not None:
            np.copyto(grad_vf, vf_gradient, casting="unsafe")
        else:
            np.subtract(vf_bar, vf_np, out=grad_vf)

        if grad_topo is not None and deg_bar is not None:
            np.copyto(grad_topo, deg_bar, casting="unsafe")
            grad_topo -= deg_array

        if w_phase != 0.0:
            np.multiply(grad_phase, w_phase, out=grad_total)
        else:
            grad_total.fill(0.0)
        if w_epi != 0.0:
            if w_epi != 1.0 and epi_gradient is None:
                np.multiply(grad_epi, w_epi, out=grad_epi)
            with np.errstate(over="ignore", invalid="ignore"):
                np.add(grad_total, grad_epi, out=grad_total)
        if w_vf != 0.0:
            if w_vf != 1.0 and vf_gradient is None:
                np.multiply(grad_vf, w_vf, out=grad_vf)
            with np.errstate(over="ignore", invalid="ignore"):
                np.add(grad_total, grad_vf, out=grad_total)
        if w_topo != 0.0 and grad_topo is not None:
            if w_topo != 1.0:
                np.multiply(grad_topo, w_topo, out=grad_topo)
            with np.errstate(over="ignore", invalid="ignore"):
                np.add(grad_total, grad_topo, out=grad_total)

        dnfr_values = grad_total
    else:
        effective_jobs = _resolve_parallel_jobs(n_jobs, len(nodes))
        if effective_jobs:
            chunk_results = []
            with ProcessPoolExecutor(max_workers=effective_jobs) as executor:
                futures = []
                for start, end in _iter_chunk_offsets(len(nodes), effective_jobs):
                    if start == end:
                        continue
                    futures.append(
                        executor.submit(
                            _dnfr_gradients_worker,
                            start,
                            end,
                            nodes,
                            [_pressure_phase_coordinate(value) for value in theta],
                            epi,
                            vf,
                            [_pressure_phase_coordinate(value) for value in th_bar],
                            epi_bar,
                            vf_bar,
                            deg_bar,
                            degs,
                            w_phase,
                            w_epi,
                            w_vf,
                            w_topo,
                            epi_gradient,
                            vf_gradient,
                            phase_overrides,
                        )
                    )
                for future in futures:
                    chunk_results.append(future.result())

            dnfr_values = [0.0] * len(nodes)
            for start, chunk in sorted(chunk_results, key=lambda item: item[0]):
                end = start + len(chunk)
                dnfr_values[start:end] = chunk
        else:
            dnfr_values = []
            for i, n in enumerate(nodes):
                g_phase = (
                    phase_overrides[i]
                    if i in phase_overrides
                    else angle_diff(
                        _pressure_phase_coordinate(th_bar[i]),
                        _pressure_phase_coordinate(theta[i]),
                    )
                    / math.pi
                )
                g_epi = (
                    (
                        epi_gradient[i]
                        if epi_gradient is not None
                        else epi_bar[i] - epi[i]
                    )
                    if w_epi
                    else 0.0
                )
                g_vf = (
                    (vf_gradient[i] if vf_gradient is not None else vf_bar[i] - vf[i])
                    if w_vf
                    else 0.0
                )
                if w_topo != 0.0 and deg_bar is not None and degs is not None:
                    if isinstance(degs, dict):
                        deg_i = float(degs.get(n, 0))
                    else:
                        deg_i = float(degs[i])
                    g_topo = deg_bar[i] - deg_i
                else:
                    g_topo = 0.0
                dnfr_values.append(
                    w_phase * g_phase
                    + (g_epi if epi_gradient is not None else w_epi * g_epi)
                    + (g_vf if vf_gradient is not None else w_vf * g_vf)
                    + w_topo * g_topo
                )

        if cache is not None:
            cache.grad_phase_np = None
            cache.grad_epi_np = None
            cache.grad_vf_np = None
            cache.grad_topo_np = None
            cache.grad_total_np = None

    stop_timer("dnfr_gradient_assembly", grad_timer)

    _require_finite_pressure(dnfr_values)
    write_timer = start_timer()
    for i, n in enumerate(nodes):
        set_dnfr(G, n, float(dnfr_values[i]))
    stop_timer("dnfr_inplace_write", write_timer)


def _init_bar_arrays(
    data: MutableMapping[str, Any],
    *,
    degs: Mapping[Any, float] | Sequence[float] | None = None,
) -> tuple[Sequence[float], Sequence[float], Sequence[float], Sequence[float] | None]:
    """Prepare containers for neighbour means.

    If NumPy is available, NumPy arrays are created; otherwise lists are used.
    ``degs`` is optional and only initialised when the topological term is
    active.
    """
    from ..mathematics.unified_numerical import np

    nodes = data["nodes"]
    theta = data["theta"]
    epi = data["epi"]
    vf = data["vf"]
    w_topo = data["w_topo"]
    cache: DnfrCache | None = data.get("cache")

    if np is not None:
        size = len(theta)
        if cache is not None:
            th_bar = cache.th_bar_np
            if th_bar is None or getattr(th_bar, "shape", None) != (size,):
                th_bar = np.array(theta, dtype=float)
            else:
                np.copyto(th_bar, theta, casting="unsafe")
            cache.th_bar_np = th_bar

            epi_bar = cache.epi_bar_np
            if epi_bar is None or getattr(epi_bar, "shape", None) != (size,):
                epi_bar = np.array(epi, dtype=float)
            else:
                np.copyto(epi_bar, epi, casting="unsafe")
            cache.epi_bar_np = epi_bar

            vf_bar = cache.vf_bar_np
            if vf_bar is None or getattr(vf_bar, "shape", None) != (size,):
                vf_bar = np.array(vf, dtype=float)
            else:
                np.copyto(vf_bar, vf, casting="unsafe")
            cache.vf_bar_np = vf_bar

            if w_topo != 0.0 and degs is not None:
                if isinstance(degs, dict):
                    deg_size = len(nodes)
                else:
                    deg_size = len(degs)
                deg_bar = cache.deg_bar_np
                if deg_bar is None or getattr(deg_bar, "shape", None) != (deg_size,):
                    if isinstance(degs, dict):
                        deg_bar = np.array(
                            [float(degs.get(node, 0.0)) for node in nodes],
                            dtype=float,
                        )
                    else:
                        deg_bar = np.array(degs, dtype=float)
                else:
                    if isinstance(degs, dict):
                        for i, node in enumerate(nodes):
                            deg_bar[i] = float(degs.get(node, 0.0))
                    else:
                        np.copyto(deg_bar, degs, casting="unsafe")
                cache.deg_bar_np = deg_bar
            else:
                deg_bar = None
                if cache is not None:
                    cache.deg_bar_np = None
        else:
            th_bar = np.array(theta, dtype=float)
            epi_bar = np.array(epi, dtype=float)
            vf_bar = np.array(vf, dtype=float)
            deg_bar = (
                np.array(degs, dtype=float)
                if w_topo != 0.0 and degs is not None
                else None
            )
    else:
        size = len(theta)
        if cache is not None:
            th_bar = cache.th_bar
            if th_bar is None or len(th_bar) != size:
                th_bar = [0.0] * size
            th_bar[:] = theta
            cache.th_bar = th_bar

            epi_bar = cache.epi_bar
            if epi_bar is None or len(epi_bar) != size:
                epi_bar = [0.0] * size
            epi_bar[:] = epi
            cache.epi_bar = epi_bar

            vf_bar = cache.vf_bar
            if vf_bar is None or len(vf_bar) != size:
                vf_bar = [0.0] * size
            vf_bar[:] = vf
            cache.vf_bar = vf_bar

            if w_topo != 0.0 and degs is not None:
                if isinstance(degs, dict):
                    deg_size = len(nodes)
                else:
                    deg_size = len(degs)
                deg_bar = cache.deg_bar
                if deg_bar is None or len(deg_bar) != deg_size:
                    deg_bar = [0.0] * deg_size
                if isinstance(degs, dict):
                    for i, node in enumerate(nodes):
                        deg_bar[i] = float(degs.get(node, 0.0))
                else:
                    for i, value in enumerate(degs):
                        deg_bar[i] = float(value)
                cache.deg_bar = deg_bar
            else:
                deg_bar = None
                cache.deg_bar = None
        else:
            th_bar = list(theta)
            epi_bar = list(epi)
            vf_bar = list(vf)
            deg_bar = list(degs) if w_topo != 0.0 and degs is not None else None
    return th_bar, epi_bar, vf_bar, deg_bar


def _compute_neighbor_means(
    G: TNFRGraph,
    data: MutableMapping[str, Any],
    *,
    x: Sequence[float],
    y: Sequence[float],
    epi_sum: Sequence[float],
    vf_sum: Sequence[float],
    count: Sequence[float] | Any,
    epi_count: Sequence[float] | None = None,
    deg_sum: Sequence[float] | None = None,
    degs: Mapping[Any, float] | Sequence[float] | None = None,
) -> tuple[Sequence[float], Sequence[float], Sequence[float], Sequence[float] | None]:
    """Return neighbour mean arrays for ΔNFR."""
    from ..mathematics.unified_numerical import np

    w_topo = data["w_topo"]
    theta = data["theta"]
    cache: DnfrCache | None = data.get("cache")
    is_numpy = np is not None and isinstance(count, np.ndarray)
    if epi_count is None:
        epi_count = count
    th_bar, epi_bar, vf_bar, deg_bar = _init_bar_arrays(data, degs=degs)

    if is_numpy:
        n = count.shape[0]
        mask = count > 0
        if not np.any(mask):
            return th_bar, epi_bar, vf_bar, deg_bar

        inv = _ensure_cached_array(cache, "neighbor_inv_count_np", (n,))
        inv.fill(0.0)
        np.divide(1.0, count, out=inv, where=mask)

        temp = _ensure_cached_array(cache, "neighbor_mean_tmp_np", (n,))
        np.arctan2(y, x, out=temp)

        theta_src = data.get("theta_np")
        if theta_src is None:
            theta_src = np.asarray(theta, dtype=float)
        # Computational extension at zero consumed floating sums. A tiny
        # nonzero resultant is not a derived physical zero or a stable target.
        zero_mask = (np.asarray(x) == 0.0) & (np.asarray(y) == 0.0)
        np.copyto(temp, theta_src, where=zero_mask)
        np.copyto(th_bar, temp, where=mask, casting="unsafe")

        if not data.get("stable_linear_gradients"):
            epi_denominator = np.asarray(epi_count, dtype=float)
            np.divide(
                epi_sum, epi_denominator, out=epi_bar, where=epi_denominator > 0.0
            )
            np.divide(vf_sum, count, out=vf_bar, where=mask)
        if w_topo != 0.0 and deg_bar is not None and deg_sum is not None:
            np.divide(deg_sum, count, out=deg_bar, where=mask)
        return th_bar, epi_bar, vf_bar, deg_bar

    n = len(theta)
    for i in range(n):
        c = count[i]
        if not c:
            continue
        inv = 1.0 / float(c)
        if x[i] == 0.0 and y[i] == 0.0:
            th_bar[i] = theta[i]
        else:
            th_bar[i] = math.atan2(y[i], x[i])
        if not data.get("stable_linear_gradients"):
            if epi_count[i] > 0.0:
                epi_bar[i] = epi_sum[i] / epi_count[i]
            vf_bar[i] = vf_sum[i] * inv
        if w_topo != 0.0 and deg_bar is not None and deg_sum is not None:
            deg_bar[i] = deg_sum[i] * inv
    return th_bar, epi_bar, vf_bar, deg_bar


def _linear_neighbor_gradients(
    G: TNFRGraph, data: Mapping[str, Any]
) -> tuple[Sequence[float], Sequence[float]]:
    """Compute linear differences using live weights and unique neighbors.

    Frequency and topology use unweighted arithmetic neighborhood means;
    the separate phase channel uses an unweighted circular mean. Only
    EPI implements the weighted diffusion ``D^-1 W EPI - EPI``; a zero
    weighted degree leaves that channel inactive, including at sink nodes.
    """
    nodes = data["nodes"]
    indices = data["idx"]
    epi = data["epi"]
    vf = data["vf"]
    if np is not None:
        source, target = data.get("edge_src"), data.get("edge_dst")
        if source is None or target is None:
            source, target = _build_edge_index_arrays(G, nodes, indices)
        weights = (
            _build_edge_weight_array(G, nodes, source, target)
            if data["w_epi"] != 0.0
            else None
        )
        return (
            edge_mean_differences(
                epi, source, target, weights, coefficient=data["w_epi"]
            ),
            edge_mean_differences(vf, source, target, coefficient=data["w_vf"]),
        )
    epi_gradient: list[float] = []
    vf_gradient: list[float] = []
    multigraph = G.is_multigraph()
    for node in nodes:
        neighbors = list(G.neighbors(node))
        if data["w_epi"] != 0.0:
            weights = []
            for neighbor in neighbors:
                edge_data = G[node][neighbor]
                weights.append(
                    float(
                        sum(attrs.get("weight", 1.0) for attrs in edge_data.values())
                        if multigraph
                        else edge_data.get("weight", 1.0)
                    )
                )
            epi_gradient.append(
                mean_neighbor_difference(
                    float(epi[indices[node]]),
                    [float(epi[indices[n]]) for n in neighbors],
                    weights,
                    coefficient=data["w_epi"],
                )
            )
        else:
            epi_gradient.append(0.0)
        vf_gradient.append(
            mean_neighbor_difference(
                float(vf[indices[node]]),
                [float(vf[indices[n]]) for n in neighbors],
                coefficient=data["w_vf"],
            )
            if data["w_vf"] != 0.0
            else 0.0
        )
    return epi_gradient, vf_gradient


def _compute_dnfr_common(
    G: TNFRGraph,
    data: MutableMapping[str, Any],
    *,
    x: Sequence[float],
    y: Sequence[float],
    epi_sum: Sequence[float],
    vf_sum: Sequence[float],
    count: Sequence[float] | None,
    deg_sum: Sequence[float] | None = None,
    degs: Sequence[float] | None = None,
    n_jobs: int | None = None,
    profile: MutableMapping[str, Any] | None = None,
) -> None:
    """Compute neighbour means and apply ΔNFR gradients.

    Parameters
    ----------
    profile : MutableMapping[str, float] or None, optional
        Mutable mapping that records wall-clock durations for the neighbour
        mean computation (``"dnfr_neighbor_means"``), the gradient assembly
        (``"dnfr_gradient_assembly"``) and the final in-place writes to the
        graph (``"dnfr_inplace_write"``).
    """
    start_timer, stop_timer = _profile_start_stop(
        profile,
        keys=("dnfr_neighbor_means", "dnfr_gradient_assembly", "dnfr_inplace_write"),
    )

    neighbor_timer = start_timer()
    data["stable_linear_gradients"] = True
    data["epi_gradient"], data["vf_gradient"] = _linear_neighbor_gradients(G, data)
    th_bar, epi_bar, vf_bar, deg_bar = _compute_neighbor_means(
        G,
        data,
        x=x,
        y=y,
        epi_sum=epi_sum,
        vf_sum=vf_sum,
        count=count,
        deg_sum=deg_sum,
        degs=degs,
    )
    stop_timer("dnfr_neighbor_means", neighbor_timer)
    _apply_dnfr_gradients(
        G,
        data,
        th_bar,
        epi_bar,
        vf_bar,
        deg_bar,
        degs,
        n_jobs=n_jobs,
        profile=profile,
    )


def _compute_delta_nfr_from_relative_sources(
    G, *, nodes, phase, neighbors, phase_sources, profile=None
) -> None:
    """Mix admitted relative-resultant phase sources with native EPI pressure.

    This private relational adapter consumes the same captured Arg(C+iS)/pi
    as each explicitly admitted relational phase law. It never reconstructs
    that branch by subtracting a rounded global circular mean. The caller owns
    the selected phase-domain admission and any resultant certification;
    matching live node order, phase and support are checked here, but supplied
    source values are not independently authenticated.

    Only the configured EPI/phase mixture is supported. All values and final
    pressures are validated before graph writes; existing default, fused and
    phase-only pressure paths retain their numerical contracts.
    """
    nodes = tuple(nodes)
    if nodes != tuple(G):
        raise ValueError("relative phase sources require the live node order")
    if (
        G.is_directed()
        or G.is_multigraph()
        or any(
            left == right
            or finite_represented_real(attrs.get("weight", 1.0), "conductance")[0]
            != 1.0
            for left, right, attrs in G.edges(data=True)
        )
    ):
        raise ValueError("relative phase sources require simple unit support")
    indices = {node: index for index, node in enumerate(nodes)}
    captured_phase = tuple(
        finite_represented_real(value, "captured pressure phase")[0] for value in phase
    )
    live_phase = tuple(_read_phase(G.nodes[node]) for node in nodes)
    if captured_phase != live_phase:
        raise ValueError("relative phase sources require the live phase state")
    captured_neighbors = tuple(tuple(row) for row in neighbors)
    live_neighbors = tuple(
        tuple(indices[neighbor] for neighbor in G.neighbors(node)) for node in nodes
    )
    if (
        any(type(index) is not int for row in captured_neighbors for index in row)
        or captured_neighbors != live_neighbors
    ):
        raise ValueError("relative phase sources require the live support order")
    sources = tuple(
        finite_represented_real(value, "relative phase source")[1]
        for value in phase_sources
    )
    if len(sources) != len(nodes) or any(abs(value) >= 1 for value in sources):
        raise ValueError("relative phase sources must match nodes and lie in (-1,1)")
    weights = _resolve_dnfr_weights(G)
    if (
        weights.get("vf", 0.0) != 0.0
        or weights.get("topo", 0.0) != 0.0
        or weights.get("phase", 0.0) <= 0.0
    ):
        raise ValueError("relative sources require the EPI/positive-phase mixture")
    data = {
        "nodes": nodes,
        "idx": indices,
        "epi": tuple(_read_scalar_epi(G.nodes[node]) for node in nodes),
        "vf": tuple(_read_capacity(G.nodes[node]) for node in nodes),
        "w_epi": weights.get("epi", 0.0),
        "w_vf": 0.0,
    }
    epi_gradient, _ = _linear_neighbor_gradients(G, data)
    for index, epi_term in enumerate(epi_gradient):
        if epi_term == 0.0 and data["w_epi"] != 0.0 and live_neighbors[index]:
            exact_difference = sum(
                (
                    Fraction(data["epi"][neighbor]) - Fraction(data["epi"][index])
                    for neighbor in live_neighbors[index]
                ),
                Fraction(0),
            )
            if exact_difference:
                raise ValueError("nonzero EPI pressure underflows to represented zero")
    phase_weight = Fraction(weights["phase"])
    pressures = tuple(
        finite_represented_real(
            Fraction(float(epi_term)) + phase_weight * source,
            "relative-source pressure",
        )[0]
        for epi_term, source in zip(epi_gradient, sources, strict=True)
    )
    _require_finite_pressure(pressures)
    for node, pressure in zip(nodes, pressures, strict=True):
        set_dnfr(G, node, pressure)
    _write_dnfr_metadata(
        G,
        weights=weights,
        hook_name="relative_resultant_canonical",
        note="relative_resultant_canonical: shared captured relational phase source",
    )
    if profile is not None:
        profile["dnfr_path"] = "relative_resultant_canonical"


def _reset_numpy_buffer(
    buffer: Any | None,
    size: int,
) -> Any:
    from ..mathematics.unified_numerical import np

    if (
        buffer is None
        or getattr(buffer, "shape", None) is None
        or buffer.shape[0] != size
    ):
        return np.zeros(size, dtype=float)
    buffer.fill(0.0)
    return buffer


def _init_neighbor_sums(
    data: MutableMapping[str, Any],
) -> NeighborStats:
    """Initialise containers for neighbour sums."""
    from ..mathematics.unified_numerical import np

    nodes = data["nodes"]
    n = len(nodes)
    w_topo = data["w_topo"]
    cache: DnfrCache | None = data.get("cache")

    def _reset_list(buffer: list[float] | None, value: float = 0.0) -> list[float]:
        if buffer is None or len(buffer) != n:
            return [value] * n
        for i in range(n):
            buffer[i] = value
        return buffer

    if np is not None:
        if cache is not None:
            x = cache.neighbor_x_np
            y = cache.neighbor_y_np
            epi_sum = cache.neighbor_epi_sum_np
            vf_sum = cache.neighbor_vf_sum_np
            count = cache.neighbor_count_np
            x = _reset_numpy_buffer(x, n)
            y = _reset_numpy_buffer(y, n)
            epi_sum = _reset_numpy_buffer(epi_sum, n)
            vf_sum = _reset_numpy_buffer(vf_sum, n)
            count = _reset_numpy_buffer(count, n)
            cache.neighbor_x_np = x
            cache.neighbor_y_np = y
            cache.neighbor_epi_sum_np = epi_sum
            cache.neighbor_vf_sum_np = vf_sum
            cache.neighbor_count_np = count
            cache.neighbor_x = _reset_list(cache.neighbor_x)
            cache.neighbor_y = _reset_list(cache.neighbor_y)
            cache.neighbor_epi_sum = _reset_list(cache.neighbor_epi_sum)
            cache.neighbor_vf_sum = _reset_list(cache.neighbor_vf_sum)
            cache.neighbor_count = _reset_list(cache.neighbor_count)
            if w_topo != 0.0:
                deg_sum = _reset_numpy_buffer(cache.neighbor_deg_sum_np, n)
                cache.neighbor_deg_sum_np = deg_sum
                cache.neighbor_deg_sum = _reset_list(cache.neighbor_deg_sum)
            else:
                cache.neighbor_deg_sum_np = None
                cache.neighbor_deg_sum = None
                deg_sum = None
        else:
            x = np.zeros(n, dtype=float)
            y = np.zeros(n, dtype=float)
            epi_sum = np.zeros(n, dtype=float)
            vf_sum = np.zeros(n, dtype=float)
            count = np.zeros(n, dtype=float)
            deg_sum = np.zeros(n, dtype=float) if w_topo != 0.0 else None
        degs = None
    else:
        if cache is not None:
            x = _reset_list(cache.neighbor_x)
            y = _reset_list(cache.neighbor_y)
            epi_sum = _reset_list(cache.neighbor_epi_sum)
            vf_sum = _reset_list(cache.neighbor_vf_sum)
            count = _reset_list(cache.neighbor_count)
            cache.neighbor_x = x
            cache.neighbor_y = y
            cache.neighbor_epi_sum = epi_sum
            cache.neighbor_vf_sum = vf_sum
            cache.neighbor_count = count
            if w_topo != 0.0:
                deg_sum = _reset_list(cache.neighbor_deg_sum)
                cache.neighbor_deg_sum = deg_sum
            else:
                cache.neighbor_deg_sum = None
                deg_sum = None
        else:
            x = [0.0] * n
            y = [0.0] * n
            epi_sum = [0.0] * n
            vf_sum = [0.0] * n
            count = [0.0] * n
            deg_sum = [0.0] * n if w_topo != 0.0 else None
        deg_list = data.get("deg_list")
        if w_topo != 0.0 and deg_list is not None:
            degs = deg_list
        else:
            degs = None
    return x, y, epi_sum, vf_sum, count, deg_sum, degs


def _prefer_sparse_accumulation(n: int, edge_count: int | None) -> bool:
    """Return ``True`` when neighbour sums should use edge accumulation."""

    if n <= 1 or not edge_count:
        return False
    possible_edges = n * (n - 1)
    if possible_edges <= 0:
        return False
    density = edge_count / possible_edges
    return density <= _SPARSE_DENSITY_THRESHOLD


def _accumulate_neighbors_dense(
    G: TNFRGraph,
    data: MutableMapping[str, Any],
    *,
    x: Any,
    y: Any,
    epi_sum: Any,
    vf_sum: Any,
    count: Any,
    deg_sum: Any | None,
) -> NeighborStats:
    """Vectorised neighbour accumulation using a dense adjacency matrix."""

    nodes = data["nodes"]
    if not nodes:
        return x, y, epi_sum, vf_sum, count, deg_sum, None

    A = data.get("A")
    if A is None:
        return _accumulate_neighbors_numpy(
            G,
            data,
            x=x,
            y=y,
            epi_sum=epi_sum,
            vf_sum=vf_sum,
            count=count,
            deg_sum=deg_sum,
        )

    cache: DnfrCache | None = data.get("cache")
    n = len(nodes)

    state = _ensure_numpy_state_vectors(data)
    linear_zero = (
        np.zeros(n, dtype=float) if data.get("stable_linear_gradients") else None
    )
    vectors = [
        state["cos"],
        state["sin"],
        linear_zero if linear_zero is not None else state["epi"],
        linear_zero if linear_zero is not None else state["vf"],
    ]

    components = _ensure_cached_array(cache, "dense_components_np", (n, 4))
    accum = _ensure_cached_array(cache, "dense_accum_np", (n, 4))

    # ``components`` retains the last source copies so callers relying on
    # cached buffers (e.g. diagnostics) still observe meaningful values.
    np.copyto(components, np.column_stack(vectors), casting="unsafe")

    np.matmul(A, components, out=accum)

    np.copyto(x, accum[:, 0], casting="unsafe")
    np.copyto(y, accum[:, 1], casting="unsafe")
    np.copyto(epi_sum, accum[:, 2], casting="unsafe")
    np.copyto(vf_sum, accum[:, 3], casting="unsafe")

    degree_counts = data.get("dense_degree_np")
    if degree_counts is None or getattr(degree_counts, "shape", (0,))[0] != n:
        degree_counts = None
    if degree_counts is None and cache is not None:
        cached_counts = cache.dense_degree_np
        if cached_counts is not None and getattr(cached_counts, "shape", (0,))[0] == n:
            degree_counts = cached_counts
    if degree_counts is None:
        degree_counts = A.sum(axis=1)
        if cache is not None:
            cache.dense_degree_np = degree_counts
    data["dense_degree_np"] = degree_counts
    np.copyto(count, degree_counts, casting="unsafe")

    degs = None
    if deg_sum is not None:
        deg_array = data.get("deg_array")
        if deg_array is None:
            deg_array = _resolve_numpy_degree_array(
                data,
                count,
                cache=cache,
            )
        if deg_array is None:
            deg_sum.fill(0.0)
        else:
            np.matmul(A, deg_array, out=deg_sum)
            degs = deg_array

    return x, y, epi_sum, vf_sum, count, deg_sum, degs


def _accumulate_neighbors_broadcasted(
    *,
    edge_src: Any,
    edge_dst: Any,
    cos: Any,
    sin: Any,
    epi: Any,
    vf: Any,
    x: Any,
    y: Any,
    epi_sum: Any,
    vf_sum: Any,
    count: Any | None,
    deg_sum: Any | None,
    deg_array: Any | None,
    cache: DnfrCache | None,
    chunk_size: int | None = None,
) -> dict[str, Any]:
    """Accumulate neighbour contributions using direct indexed reductions.

    Array reuse strategy for non-chunked blocks:
    --------------------------------------------
    This function optimizes memory usage by reusing cached destination arrays:

    1. **Accumulator reuse**: The `accum` matrix (component_rows × n) is cached
       across invocations when signature remains stable. For non-chunked paths,
       it's zero-filled (accum.fill(0.0)) rather than reallocated.

    2. **Workspace reuse**: The `workspace` buffer (component_rows × edge_count)
       stores intermediate edge values. In non-chunked mode with sufficient
       workspace size, edge values are extracted once into workspace rows
       via np.take(..., out=workspace[row, :]) to avoid repeated allocations.

    3. **Destination array writes**: np.bincount results are written to accum
       rows via np.copyto(..., casting="unsafe"), reusing the same memory
       across all components (cos, sin, epi, vf, count, deg).

    4. **Deg_array optimization**: When deg_array is provided and topology
       weight is active, degree values are extracted into workspace and
       accumulated via bincount, maintaining the reuse pattern.

    The non-chunked path achieves minimal temporary allocations by:
    - Reusing cached accum and workspace buffers
    - Extracting all edge values into workspace in a single pass
    - Writing bincount results directly to destination rows

    Note: np.bincount does not support an `out` parameter, so its return
    value must be copied to the destination. The workspace pattern minimizes
    the number of temporary arrays created during edge value extraction.

    This approach maintains ΔNFR computational accuracy (Invariant #1:
    Nodal Equation Integrity) while
    reducing memory footprint for repeated accumulations with stable topology.
    """
    from ..mathematics.unified_numerical import np

    n = x.shape[0]
    edge_count = int(edge_src.size)

    include_count = count is not None
    use_topology = deg_sum is not None and deg_array is not None

    component_rows = 4 + (1 if include_count else 0) + (1 if use_topology else 0)

    if edge_count:
        if chunk_size is None:
            resolved_chunk = edge_count
        else:
            try:
                resolved_chunk = int(chunk_size)
            except (TypeError, ValueError):
                resolved_chunk = edge_count
            else:
                if resolved_chunk <= 0:
                    resolved_chunk = edge_count
        resolved_chunk = max(1, min(edge_count, resolved_chunk))
    else:
        resolved_chunk = 0

    use_chunks = bool(edge_count and resolved_chunk < edge_count)

    if cache is not None:
        base_signature = (id(edge_src), id(edge_dst), n, edge_count)
        cache.edge_signature = base_signature
        signature = (base_signature, component_rows)
        previous_signature = cache.neighbor_accum_signature

        accum = cache.neighbor_accum_np
        if (
            accum is None
            or getattr(accum, "shape", None) != (component_rows, n)
            or previous_signature != signature
        ):
            accum = np.zeros((component_rows, n), dtype=float)
            cache.neighbor_accum_np = accum
        else:
            accum.fill(0.0)

        workspace = cache.neighbor_edge_values_np
        if use_chunks:
            workspace_length = resolved_chunk
        else:
            # For non-chunked path, allocate workspace to hold edge_count values
            # so we can extract edge values without temporary allocations
            workspace_length = edge_count if edge_count else component_rows
        if workspace_length:
            expected_shape = (component_rows, workspace_length)
            if workspace is None or getattr(workspace, "shape", None) != expected_shape:
                workspace = np.empty(expected_shape, dtype=float)
        else:
            workspace = None
        cache.neighbor_edge_values_np = workspace

        cache.neighbor_accum_signature = signature
    else:
        accum = np.zeros((component_rows, n), dtype=float)
        # For non-chunked path without cache, allocate workspace for edge values
        workspace_length = (
            edge_count
            if (not use_chunks and edge_count)
            else (resolved_chunk if use_chunks else component_rows)
        )
        workspace = (
            np.empty((component_rows, workspace_length), dtype=float)
            if workspace_length
            else None
        )

    if edge_count:
        row = 0
        cos_row = row
        row += 1
        sin_row = row
        row += 1
        epi_row = row
        row += 1
        vf_row = row
        row += 1
        count_row = row if include_count and count is not None else None
        if count_row is not None:
            row += 1
        deg_row = row if use_topology and deg_array is not None else None

        edge_src_int = edge_src.astype(np.intp, copy=False)
        edge_dst_int = edge_dst.astype(np.intp, copy=False)

        if use_chunks:
            chunk_step = resolved_chunk if resolved_chunk else edge_count
            chunk_indices = range(0, edge_count, chunk_step)

            for start in chunk_indices:
                end = min(start + chunk_step, edge_count)
                if start >= end:
                    continue
                src_slice = edge_src_int[start:end]
                dst_slice = edge_dst_int[start:end]
                slice_len = end - start
                if slice_len <= 0:
                    continue

                if workspace is not None:
                    chunk_matrix = workspace[:, :slice_len]
                else:
                    chunk_matrix = np.empty((component_rows, slice_len), dtype=float)

                np.take(cos, dst_slice, out=chunk_matrix[cos_row, :slice_len])
                np.take(sin, dst_slice, out=chunk_matrix[sin_row, :slice_len])
                np.take(epi, dst_slice, out=chunk_matrix[epi_row, :slice_len])
                np.take(vf, dst_slice, out=chunk_matrix[vf_row, :slice_len])

                if count_row is not None:
                    chunk_matrix[count_row, :slice_len].fill(1.0)
                if deg_row is not None and deg_array is not None:
                    np.take(deg_array, dst_slice, out=chunk_matrix[deg_row, :slice_len])

                def _accumulate_into(
                    target_row: int | None,
                    values: np.ndarray | None = None,
                    *,
                    unit_weight: bool = False,
                ) -> None:
                    if target_row is None:
                        return
                    row_view = accum[target_row]
                    if unit_weight:
                        np.add.at(row_view, src_slice, 1.0)
                    else:
                        if values is None:
                            return
                        np.add.at(row_view, src_slice, values)

                _accumulate_into(cos_row, chunk_matrix[cos_row, :slice_len])
                _accumulate_into(sin_row, chunk_matrix[sin_row, :slice_len])
                _accumulate_into(epi_row, chunk_matrix[epi_row, :slice_len])
                _accumulate_into(vf_row, chunk_matrix[vf_row, :slice_len])

                if count_row is not None:
                    _accumulate_into(count_row, unit_weight=True)

                if deg_row is not None and deg_array is not None:
                    _accumulate_into(deg_row, chunk_matrix[deg_row, :slice_len])
        else:
            # Non-chunked path: reuse workspace to minimize temporary allocations.
            # When workspace is available with sufficient size, extract edge values
            # into workspace rows before passing to bincount.
            if workspace is not None and workspace.shape[1] >= edge_count:
                # Verify workspace has enough rows for all components
                # workspace has shape (component_rows, edge_count)
                required_rows = max(
                    cos_row + 1,
                    sin_row + 1,
                    epi_row + 1,
                    vf_row + 1,
                    (count_row + 1) if count_row is not None else 0,
                    (deg_row + 1) if deg_row is not None else 0,
                )
                if workspace.shape[0] >= required_rows:
                    # Reuse workspace rows for edge value extraction
                    np.take(cos, edge_dst_int, out=workspace[cos_row, :edge_count])
                    np.take(sin, edge_dst_int, out=workspace[sin_row, :edge_count])
                    np.take(epi, edge_dst_int, out=workspace[epi_row, :edge_count])
                    np.take(vf, edge_dst_int, out=workspace[vf_row, :edge_count])

                    def _apply_full_bincount(
                        target_row: int | None,
                        values: np.ndarray | None = None,
                        *,
                        unit_weight: bool = False,
                    ) -> None:
                        if target_row is None:
                            return
                        if values is None and not unit_weight:
                            return
                        if unit_weight:
                            component_accum = np.bincount(
                                edge_src_int,
                                minlength=n,
                            )
                        else:
                            component_accum = np.bincount(
                                edge_src_int,
                                weights=values,
                                minlength=n,
                            )
                        np.copyto(
                            accum[target_row, :n],
                            component_accum[:n],
                            casting="unsafe",
                        )

                    _apply_full_bincount(cos_row, workspace[cos_row, :edge_count])
                    _apply_full_bincount(sin_row, workspace[sin_row, :edge_count])
                    _apply_full_bincount(epi_row, workspace[epi_row, :edge_count])
                    _apply_full_bincount(vf_row, workspace[vf_row, :edge_count])

                    if count_row is not None:
                        _apply_full_bincount(count_row, unit_weight=True)

                    if deg_row is not None and deg_array is not None:
                        np.take(
                            deg_array, edge_dst_int, out=workspace[deg_row, :edge_count]
                        )
                        _apply_full_bincount(deg_row, workspace[deg_row, :edge_count])
                else:
                    # Workspace doesn't have enough rows, fall back to temporary arrays
                    def _apply_full_bincount(
                        target_row: int | None,
                        values: np.ndarray | None = None,
                        *,
                        unit_weight: bool = False,
                    ) -> None:
                        if target_row is None:
                            return
                        if values is None and not unit_weight:
                            return
                        if unit_weight:
                            component_accum = np.bincount(
                                edge_src_int,
                                minlength=n,
                            )
                        else:
                            component_accum = np.bincount(
                                edge_src_int,
                                weights=values,
                                minlength=n,
                            )
                        np.copyto(
                            accum[target_row, :n],
                            component_accum[:n],
                            casting="unsafe",
                        )

                    _apply_full_bincount(cos_row, np.take(cos, edge_dst_int))
                    _apply_full_bincount(sin_row, np.take(sin, edge_dst_int))
                    _apply_full_bincount(epi_row, np.take(epi, edge_dst_int))
                    _apply_full_bincount(vf_row, np.take(vf, edge_dst_int))

                    if count_row is not None:
                        _apply_full_bincount(count_row, unit_weight=True)

                    if deg_row is not None and deg_array is not None:
                        _apply_full_bincount(deg_row, np.take(deg_array, edge_dst_int))
            else:
                # Fallback: no workspace or insufficient width, use temporary arrays
                def _apply_full_bincount(
                    target_row: int | None,
                    values: np.ndarray | None = None,
                    *,
                    unit_weight: bool = False,
                ) -> None:
                    if target_row is None:
                        return
                    if values is None and not unit_weight:
                        return
                    if unit_weight:
                        component_accum = np.bincount(
                            edge_src_int,
                            minlength=n,
                        )
                    else:
                        component_accum = np.bincount(
                            edge_src_int,
                            weights=values,
                            minlength=n,
                        )
                    np.copyto(
                        accum[target_row, :n],
                        component_accum[:n],
                        casting="unsafe",
                    )

                _apply_full_bincount(cos_row, np.take(cos, edge_dst_int))
                _apply_full_bincount(sin_row, np.take(sin, edge_dst_int))
                _apply_full_bincount(epi_row, np.take(epi, edge_dst_int))
                _apply_full_bincount(vf_row, np.take(vf, edge_dst_int))

                if count_row is not None:
                    _apply_full_bincount(count_row, unit_weight=True)

                if deg_row is not None and deg_array is not None:
                    _apply_full_bincount(deg_row, np.take(deg_array, edge_dst_int))
    else:
        accum.fill(0.0)
        if workspace is not None:
            workspace.fill(0.0)

    row = 0
    np.copyto(x, accum[row], casting="unsafe")
    row += 1
    np.copyto(y, accum[row], casting="unsafe")
    row += 1
    np.copyto(epi_sum, accum[row], casting="unsafe")
    row += 1
    np.copyto(vf_sum, accum[row], casting="unsafe")
    row += 1

    if include_count and count is not None:
        np.copyto(count, accum[row], casting="unsafe")
        row += 1

    if use_topology and deg_sum is not None:
        np.copyto(deg_sum, accum[row], casting="unsafe")

    return {
        "accumulator": accum,
        "edge_values": workspace,
    }


def _build_neighbor_sums_common(
    G: TNFRGraph,
    data: MutableMapping[str, Any],
    *,
    use_numpy: bool,
    n_jobs: int | None = None,
) -> NeighborStats:
    """Build neighbour accumulators honouring cached NumPy buffers when possible."""

    nodes = data["nodes"]
    cache: DnfrCache | None = data.get("cache")
    np_module = np
    has_numpy_buffers = _has_cached_numpy_buffers(data, cache)

    # Fallback: when np is None but we have cached NumPy buffers,
    # attempt to retrieve NumPy from sys.modules to avoid losing vectorization.
    # This preserves ΔNFR semantics (Invariant #1: Nodal Equation Integrity)
    # and maintains performance.
    if use_numpy and np_module is None and has_numpy_buffers:
        candidate = sys.modules.get("numpy")
        # Validate the candidate module has required NumPy attributes
        if (
            candidate is not None
            and hasattr(candidate, "ndarray")
            and hasattr(candidate, "empty")
        ):
            np_module = candidate

    if np_module is not None:
        if not nodes:
            return _init_neighbor_sums(data)

        x, y, epi_sum, vf_sum, count, deg_sum, degs = _init_neighbor_sums(data)

        # Reuse centralized sparse/dense decision from _prepare_dnfr_data.
        # The decision logic at lines 785-807 already computed prefer_sparse
        # and dense_override based on graph density and user flags.
        prefer_sparse = data.get("prefer_sparse")
        if prefer_sparse is None:
            # Fallback: recompute if not set (defensive, should be rare)
            prefer_sparse = _prefer_sparse_accumulation(
                len(nodes), data.get("edge_count")
            )
            data["prefer_sparse"] = prefer_sparse

        use_dense = False
        A = data.get("A")
        dense_override = data.get("dense_override", False)

        # Apply centralized decision: dense path requires adjacency matrix
        # and either high graph density or explicit dense_override flag.
        if use_numpy and A is not None:
            shape = getattr(A, "shape", (0, 0))
            matrix_valid = shape[0] == len(nodes) and shape[1] == len(nodes)
            if matrix_valid and (dense_override or not prefer_sparse):
                use_dense = True

        if use_dense:
            accumulator = _accumulate_neighbors_dense
        else:
            _ensure_numpy_state_vectors(data)
            accumulator = _accumulate_neighbors_numpy
        return accumulator(
            G,
            data,
            x=x,
            y=y,
            epi_sum=epi_sum,
            vf_sum=vf_sum,
            count=count,
            deg_sum=deg_sum,
        )

    if not nodes:
        return _init_neighbor_sums(data)

    x, y, epi_sum, vf_sum, count, deg_sum, degs_list = _init_neighbor_sums(data)
    idx = data["idx"]
    epi = data["epi"]
    vf = data["vf"]
    if data.get("stable_linear_gradients"):
        # Linear channels are reduced separately from pairwise differences.
        epi = vf = [0.0] * len(nodes)
    cos_th = data["cos_theta"]
    sin_th = data["sin_theta"]
    deg_list = data.get("deg_list")

    effective_jobs = _resolve_parallel_jobs(n_jobs, len(nodes))
    if effective_jobs:
        neighbor_indices: list[list[int]] = []
        for node in nodes:
            indices: list[int] = []
            for v in G.neighbors(node):
                indices.append(idx[v])
            neighbor_indices.append(indices)

        chunk_results = []
        with ProcessPoolExecutor(max_workers=effective_jobs) as executor:
            futures = []
            for start, end in _iter_chunk_offsets(len(nodes), effective_jobs):
                if start == end:
                    continue
                futures.append(
                    executor.submit(
                        _neighbor_sums_worker,
                        start,
                        end,
                        neighbor_indices,
                        cos_th,
                        sin_th,
                        epi,
                        vf,
                        x[start:end],
                        y[start:end],
                        epi_sum[start:end],
                        vf_sum[start:end],
                        count[start:end],
                        deg_sum[start:end] if deg_sum is not None else None,
                        deg_list,
                        degs_list,
                    )
                )
            for future in futures:
                chunk_results.append(future.result())

        for (
            start,
            chunk_x,
            chunk_y,
            chunk_epi,
            chunk_vf,
            chunk_count,
            chunk_deg,
        ) in sorted(chunk_results, key=lambda item: item[0]):
            end = start + len(chunk_x)
            x[start:end] = chunk_x
            y[start:end] = chunk_y
            epi_sum[start:end] = chunk_epi
            vf_sum[start:end] = chunk_vf
            count[start:end] = chunk_count
            if deg_sum is not None and chunk_deg is not None:
                deg_sum[start:end] = chunk_deg
        return x, y, epi_sum, vf_sum, count, deg_sum, degs_list

    for i, node in enumerate(nodes):
        deg_i = degs_list[i] if degs_list is not None else 0.0
        x_i = x[i]
        y_i = y[i]
        epi_i = epi_sum[i]
        vf_i = vf_sum[i]
        count_i = count[i]
        deg_acc = deg_sum[i] if deg_sum is not None else 0.0
        for v in G.neighbors(node):
            j = idx[v]
            cos_j = cos_th[j]
            sin_j = sin_th[j]
            epi_j = epi[j]
            vf_j = vf[j]
            x_i += cos_j
            y_i += sin_j
            epi_i += epi_j
            vf_i += vf_j
            count_i += 1
            if deg_sum is not None:
                deg_acc += deg_list[j] if deg_list is not None else deg_i
        x[i] = x_i
        y[i] = y_i
        epi_sum[i] = epi_i
        vf_sum[i] = vf_i
        count[i] = count_i
        if deg_sum is not None:
            deg_sum[i] = deg_acc
    return x, y, epi_sum, vf_sum, count, deg_sum, degs_list


def _accumulate_neighbors_numpy(
    G: TNFRGraph,
    data: MutableMapping[str, Any],
    *,
    x: Any,
    y: Any,
    epi_sum: Any,
    vf_sum: Any,
    count: Any | None,
    deg_sum: Any | None,
) -> NeighborStats:
    """Vectorised neighbour accumulation reusing cached NumPy buffers."""
    from ..mathematics.unified_numerical import np

    nodes = data["nodes"]
    if not nodes:
        return x, y, epi_sum, vf_sum, count, deg_sum, None

    cache: DnfrCache | None = data.get("cache")

    state = _ensure_numpy_state_vectors(data)
    cos_th = state["cos"]
    sin_th = state["sin"]
    epi = state["epi"]
    vf = state["vf"]
    if data.get("stable_linear_gradients"):
        epi = vf = np.zeros(len(nodes), dtype=float)

    edge_src = data.get("edge_src")
    edge_dst = data.get("edge_dst")
    if edge_src is None or edge_dst is None:
        edge_src, edge_dst = _build_edge_index_arrays(G, nodes, data["idx"])
        data["edge_src"] = edge_src
        data["edge_dst"] = edge_dst
        if cache is not None:
            cache.edge_src = edge_src
            cache.edge_dst = edge_dst
    if edge_src is not None:
        data["edge_count"] = int(edge_src.size)

    cached_deg_array = data.get("deg_array")

    # Memory optimization: When we have a cached degree array and need a count
    # buffer, we can reuse the degree array buffer as the destination for counts.
    # This works because:
    #   1. For undirected graphs, node degree equals neighbor count
    #   2. The degree array is already allocated and the right size
    #   3. We avoid allocating an extra row in the accumulator matrix
    # When reuse_count_from_deg is True:
    #   - We copy cached_deg_array into the count buffer before accumulation
    #   - We pass count_for_accum=None to _accumulate_neighbors_broadcasted
    #   - After accumulation, we restore count = cached_deg_array (line 2121)
    reuse_count_from_deg = bool(count is not None and cached_deg_array is not None)
    count_for_accum = count
    if count is not None:
        if reuse_count_from_deg:
            # Pre-fill count with degrees (will be returned as-is since accumulator
            # skips the count row when count_for_accum=None)
            np.copyto(count, cached_deg_array, casting="unsafe")
            count_for_accum = None
        else:
            count.fill(0.0)

    deg_array = None
    if deg_sum is not None:
        deg_sum.fill(0.0)
        deg_array = _resolve_numpy_degree_array(
            data, count if count is not None else None, cache=cache
        )
    elif cached_deg_array is not None:
        deg_array = cached_deg_array

    edge_count = int(edge_src.size) if edge_src is not None else 0
    chunk_hint = data.get("neighbor_chunk_hint")
    if chunk_hint is None:
        chunk_hint = G.graph.get("DNFR_CHUNK_SIZE")
    resolved_neighbor_chunk = (
        resolve_chunk_size(
            chunk_hint,
            edge_count,
            minimum=1,
            approx_bytes_per_item=_DNFR_APPROX_BYTES_PER_EDGE,
            clamp_to=None,
        )
        if edge_count
        else 0
    )
    data["neighbor_chunk_hint"] = chunk_hint
    data["neighbor_chunk_size"] = resolved_neighbor_chunk

    accum = _accumulate_neighbors_broadcasted(
        edge_src=edge_src,
        edge_dst=edge_dst,
        cos=cos_th,
        sin=sin_th,
        epi=epi,
        vf=vf,
        x=x,
        y=y,
        epi_sum=epi_sum,
        vf_sum=vf_sum,
        count=count_for_accum,
        deg_sum=deg_sum,
        deg_array=deg_array,
        cache=cache,
        chunk_size=resolved_neighbor_chunk,
    )

    data["neighbor_accum_np"] = accum.get("accumulator")
    edge_values = accum.get("edge_values")
    data["neighbor_edge_values_np"] = edge_values
    if edge_values is not None:
        width = getattr(edge_values, "shape", (0, 0))[1]
        data["neighbor_chunk_size"] = int(width)
    else:
        data["neighbor_chunk_size"] = resolved_neighbor_chunk
    if cache is not None:
        data["neighbor_accum_signature"] = cache.neighbor_accum_signature
    if reuse_count_from_deg and cached_deg_array is not None:
        count = cached_deg_array
    degs = deg_array if deg_sum is not None and deg_array is not None else None
    return x, y, epi_sum, vf_sum, count, deg_sum, degs


def _compute_dnfr(
    G: TNFRGraph,
    data: MutableMapping[str, Any],
    *,
    use_numpy: bool | None = None,
    n_jobs: int | None = None,
    profile: MutableMapping[str, Any] | None = None,
) -> None:
    """Compute ΔNFR using neighbour sums.

    Parameters
    ----------
    G : nx.Graph
        Graph on which the computation is performed.
    data : dict
        Precomputed ΔNFR data as returned by :func:`_prepare_dnfr_data`.
    use_numpy : bool | None, optional
        Backwards compatibility flag. When ``True`` the function eagerly
        prepares NumPy buffers (if available). When ``False`` the engine still
        prefers the vectorised path whenever NumPy is available
        and the graph does not set ``vectorized_dnfr`` to ``False``.
    profile : MutableMapping[str, float] or None, optional
        Mutable mapping that aggregates wall-clock durations for neighbour
        accumulation and records which execution branch was used. The
        ``"dnfr_neighbor_accumulation"`` bucket gathers the time spent inside
        :func:`_build_neighbor_sums_common`, while ``"dnfr_path"`` stores the
        string ``"vectorized"`` or ``"fallback"`` describing the active
        implementation.
    """
    start_timer, stop_timer = _profile_start_stop(
        profile,
        keys=("dnfr_neighbor_accumulation",),
    )

    np_module = np
    data["dnfr_numpy_available"] = bool(np_module)
    vector_disabled = not _should_vectorize(G, np_module)
    prefer_dense = np_module is not None and not vector_disabled
    if use_numpy is True and np_module is not None:
        prefer_dense = True
    if use_numpy is False or vector_disabled:
        prefer_dense = False
    data["dnfr_used_numpy"] = bool(prefer_dense and np_module is not None)
    if profile is not None:
        profile["dnfr_path"] = "vectorized" if data["dnfr_used_numpy"] else "fallback"

    # Optimization: Use fused kernel if available
    # This bypasses the separate accumulation/mean/gradient steps
    if data["dnfr_used_numpy"]:
        nodes = data["nodes"]
        cache = data.get("cache")
        edge_src = data.get("edge_src")
        edge_dst = data.get("edge_dst")

        if edge_src is None or edge_dst is None:
            edge_src, edge_dst = _build_edge_index_arrays(G, nodes, data["idx"])
            data["edge_src"] = edge_src
            data["edge_dst"] = edge_dst
            if cache is not None:
                cache.edge_src = edge_src
                cache.edge_dst = edge_dst

        state = _ensure_numpy_state_vectors(data)

        # Inspect individual live weights: sum(weights) == edge_count does
        # not imply unit weights (e.g. 0.5 and 1.5). The shared builder also
        # sums parallel edges and retains the exact unit-weight kernel path.
        edge_weight = (
            _build_edge_weight_array(G, nodes, edge_src, edge_dst)
            if data["w_epi"] != 0.0
            else None
        )

        # Note: accumulate_both_directions=False because _build_edge_index_arrays
        # already generates bidirectional edges for undirected graphs (via G.neighbors).
        # compute_fused_gradients_symmetric expects w_-prefixed weight keys
        # (the convention used by the optimized_numpy backend); data carries
        # them as w_phase/w_epi/w_vf/w_topo from _prepare_dnfr_data. Passing the
        # bare _dnfr_weights dict ({phase,epi,vf,topo}) silently zeroed ΔNFR.
        dnfr_values = compute_fused_gradients_symmetric(
            edge_src=edge_src,
            edge_dst=edge_dst,
            phase=state["theta"],
            epi=state["epi"],
            vf=state["vf"],
            weights={
                "w_phase": data["w_phase"],
                "w_epi": data["w_epi"],
                "w_vf": data["w_vf"],
                "w_topo": data["w_topo"],
            },
            edge_weight=edge_weight,
            accumulate_both_directions=False,
        )

        if profile is not None:
            profile["dnfr_path"] = "fused_canonical"

        write_timer = start_timer()
        for i, n in enumerate(nodes):
            set_dnfr(G, n, float(dnfr_values[i]))
        stop_timer("dnfr_inplace_write", write_timer)

        return

    data["n_jobs"] = n_jobs
    data["stable_linear_gradients"] = True
    neighbor_timer = start_timer()
    res = _build_neighbor_sums_common(
        G,
        data,
        use_numpy=prefer_dense,
        n_jobs=n_jobs,
    )
    stop_timer("dnfr_neighbor_accumulation", neighbor_timer)
    if res is None:
        return
    x, y, epi_sum, vf_sum, count, deg_sum, degs = res
    _compute_dnfr_common(
        G,
        data,
        x=x,
        y=y,
        epi_sum=epi_sum,
        vf_sum=vf_sum,
        count=count,
        deg_sum=deg_sum,
        degs=degs,
        n_jobs=n_jobs,
        profile=profile,
    )


def default_compute_delta_nfr(
    G: TNFRGraph,
    *,
    cache_size: int | None = 1,
    n_jobs: int | None = None,
    profile: MutableMapping[str, Any] | None = None,
) -> None:
    """Compute ΔNFR by mixing phase, EPI, νf and a topological term.

    Parameters
    ----------
    G : nx.Graph
        Graph on which the computation is performed.
    cache_size : int | None, optional
        Maximum number of edge configurations cached in ``G.graph``. ``None``
        is unlimited; values <= 0 disable reusable accumulation caches.
        Defaults to ``1``.
    n_jobs : int | None, optional
        Parallel worker count for the pure-Python accumulation path. ``None``
        or values <= 1 preserve the serial behaviour. The vectorised NumPy
        branch ignores this parameter as it already operates in bulk.
    profile : MutableMapping[str, float] or None, optional
        Mutable mapping that aggregates the wall-clock timings captured during
        the ΔNFR computation. The mapping receives the buckets documented in
        :func:`_prepare_dnfr_data` and :func:`_compute_dnfr`, plus
        ``"dnfr_neighbor_means"``, ``"dnfr_gradient_assembly"`` and
        ``"dnfr_inplace_write"`` describing the internal stages of
        :func:`_compute_dnfr_common`. ``"dnfr_path"`` reflects whether the
        vectorised or fallback implementation executed the call.
    """
    if profile is not None:
        for key in (
            "dnfr_cache_rebuild",
            "dnfr_neighbor_accumulation",
            "dnfr_neighbor_means",
            "dnfr_gradient_assembly",
            "dnfr_inplace_write",
        ):
            profile.setdefault(key, 0.0)

    data = _prepare_dnfr_data(G, cache_size=cache_size, profile=profile)
    _write_dnfr_metadata(
        G,
        weights=data["weights"],
        hook_name="default_compute_delta_nfr",
    )
    _compute_dnfr(G, data, n_jobs=n_jobs, profile=profile)
    if not data.get("dnfr_numpy_available"):
        cache = data.get("cache")
        cache_size = data.get("cache_size")
        caching_enabled = isinstance(cache, DnfrCache) and (
            cache_size is None or int(cache_size) > 0
        )
        if isinstance(cache, DnfrCache) and not caching_enabled:
            for attr in (
                "neighbor_x_np",
                "neighbor_y_np",
                "neighbor_epi_sum_np",
                "neighbor_vf_sum_np",
                "neighbor_count_np",
                "neighbor_deg_sum_np",
                "neighbor_inv_count_np",
                "neighbor_cos_avg_np",
                "neighbor_sin_avg_np",
                "neighbor_mean_tmp_np",
                "neighbor_mean_length_np",
                "neighbor_accum_np",
                "neighbor_edge_values_np",
            ):
                setattr(cache, attr, None)
            cache.neighbor_accum_signature = None


def _invoke_dnfr_hook(func: Callable[..., Any], *args: Any, **kwargs: Any) -> Any:
    """Call a pressure callback without retrying an internal ``TypeError``.

    Only the optional worker hint may be omitted for legacy signatures. For
    uninspectable callables, the compatibility retry requires an argument-binding
    failure at this call site, before any Python callback body was entered.
    """
    if "n_jobs" not in kwargs:
        return func(*args, **kwargs)
    try:
        signature = inspect.signature(func)
    except (TypeError, ValueError):
        signature = None
    if signature is not None:
        parameters = signature.parameters
        worker_parameter = parameters.get("n_jobs")
        supports_n_jobs = (
            worker_parameter is not None
            and worker_parameter.kind
            in (inspect.Parameter.POSITIONAL_OR_KEYWORD, inspect.Parameter.KEYWORD_ONLY)
        ) or any(
            parameter.kind == inspect.Parameter.VAR_KEYWORD
            for parameter in parameters.values()
        )
        if not supports_n_jobs:
            kwargs = {key: value for key, value in kwargs.items() if key != "n_jobs"}
        return func(*args, **kwargs)
    try:
        return func(*args, **kwargs)
    except TypeError as exc:
        traceback = exc.__traceback__
        message = str(exc)
        binding_failure = (
            traceback is not None
            and traceback.tb_next is None
            and ("n_jobs" in message or "takes no keyword arguments" in message)
        )
        if not binding_failure:
            raise
        kwargs = {key: value for key, value in kwargs.items() if key != "n_jobs"}
        return func(*args, **kwargs)


def set_delta_nfr_hook(
    G: TNFRGraph,
    func: DeltaNFRHook,
    *,
    name: str | None = None,
    note: str | None = None,
) -> None:
    """set a stable hook to compute ΔNFR.

    The callable should accept ``(G, *[, n_jobs])`` and is responsible for
    writing ``ALIAS_DNFR`` in each node. ``n_jobs`` is optional and ignored by
    hooks that do not support parallel execution. Basic metadata in
    ``G.graph`` is updated accordingly. Internal Python callback errors propagate
    without a second invocation; already performed callback writes are not
    rolled back by this adapter.
    """

    def _wrapped(graph: TNFRGraph, *args: Any, **kwargs: Any) -> None:
        _invoke_dnfr_hook(func, graph, *args, **kwargs)

    _wrapped.__name__ = getattr(func, "__name__", "custom_dnfr")
    _wrapped.__doc__ = getattr(func, "__doc__", _wrapped.__doc__)

    G.graph["compute_delta_nfr"] = _wrapped
    G.graph["_dnfr_hook_name"] = str(name or getattr(func, "__name__", "custom_dnfr"))
    if "_dnfr_weights" not in G.graph:
        _configure_dnfr_weights(G)
    if note:
        meta = G.graph.get("_DNFR_META", {})
        meta["note"] = str(note)
        G.graph["_DNFR_META"] = meta


def _weighted_hook_gradient(
    func: Callable[[TNFRGraph, NodeId, Mapping[str, Any]], float],
    G: TNFRGraph,
    node: NodeId,
    nd: Mapping[str, Any],
    coefficient: float,
) -> float:
    """Apply linear coefficients before the final difference is materialized."""
    if isinstance(func, _NeighborAverageGradient):
        return func(G, node, nd, coefficient=coefficient)
    return coefficient * float(func(G, node, nd))


def _dnfr_hook_chunk_worker(
    G: TNFRGraph,
    node_ids: Sequence[NodeId],
    grad_items: tuple[
        tuple[str, Callable[[TNFRGraph, NodeId, Mapping[str, Any]], float]],
        ...,
    ],
    weights: Mapping[str, float],
) -> list[tuple[NodeId, float]]:
    """Compute weighted gradients for ``node_ids``.

    The helper is defined at module level so it can be pickled by
    :class:`concurrent.futures.ProcessPoolExecutor`.
    """

    results: list[tuple[NodeId, float]] = []
    for node in node_ids:
        nd = G.nodes[node]
        total = 0.0
        for name, func in grad_items:
            w = float(weights.get(name, 0.0))
            if w:
                total += _weighted_hook_gradient(func, G, node, nd, w)
        results.append((node, total))
    return results


def _apply_dnfr_hook(
    G: TNFRGraph,
    grads: Mapping[str, Callable[[TNFRGraph, NodeId, Mapping[str, Any]], float]],
    *,
    weights: Mapping[str, float],
    hook_name: str,
    note: str | None = None,
    n_jobs: int | None = None,
) -> None:
    """Compute and store ΔNFR using ``grads``.

    Parameters
    ----------
    G : nx.Graph
        Graph whose nodes will receive the ΔNFR update.
    grads : dict
        Mapping from component names to callables with signature
        ``(G, node, data) -> float`` returning the gradient contribution.
    weights : dict
        Weight per component; missing entries default to ``0``.
    hook_name : str
        Friendly identifier stored in ``G.graph`` metadata.
    note : str | None, optional
        Additional documentation recorded next to the hook metadata.
    n_jobs : int | None, optional
        Optional worker count for the pure-Python execution path. When NumPy
        is available the helper always prefers the vectorised implementation
        and ignores ``n_jobs`` because the computation already happens in
        bulk.
    """

    nodes_data: list[tuple[NodeId, Mapping[str, Any]]] = list(G.nodes(data=True))
    if not nodes_data:
        _write_dnfr_metadata(G, weights=weights, hook_name=hook_name, note=note)
        return

    np_module = np
    if np_module is not None:
        totals = np_module.zeros(len(nodes_data), dtype=float)
        for name, func in grads.items():
            w = float(weights.get(name, 0.0))
            if w == 0.0:
                continue
            values = np_module.fromiter(
                (_weighted_hook_gradient(func, G, n, nd, w) for n, nd in nodes_data),
                dtype=float,
                count=len(nodes_data),
            )
            with np_module.errstate(over="ignore", invalid="ignore"):
                np_module.add(totals, values, out=totals)
        _require_finite_pressure(totals)
        for idx, (n, _) in enumerate(nodes_data):
            set_dnfr(G, n, float(totals[idx]))
        _write_dnfr_metadata(G, weights=weights, hook_name=hook_name, note=note)
        return

    effective_jobs = _resolve_parallel_jobs(n_jobs, len(nodes_data))
    results: list[tuple[NodeId, float]] | None = None
    if effective_jobs:
        grad_items = tuple(grads.items())
        # ProcessPoolExecutor requires picklable arguments. Instead of explicitly
        # testing with pickle.dumps (which poses security risks), we attempt
        # parallelization and gracefully fall back to serial on any failure.
        try:
            chunk_results: list[tuple[NodeId, float]] = []
            with ProcessPoolExecutor(max_workers=effective_jobs) as executor:
                futures = []
                node_ids: list[NodeId] = [n for n, _ in nodes_data]
                for start, end in _iter_chunk_offsets(len(node_ids), effective_jobs):
                    if start == end:
                        continue
                    futures.append(
                        executor.submit(
                            _dnfr_hook_chunk_worker,
                            G,
                            node_ids[start:end],
                            grad_items,
                            weights,
                        )
                    )
                for future in futures:
                    chunk_results.extend(future.result())
            results = chunk_results
        except Exception:
            # Parallel execution failed (pickle, executor, or worker error)
            # Fall back to serial processing
            results = None

    if results is None:
        results = []
        for n, nd in nodes_data:
            total = 0.0
            for name, func in grads.items():
                w = float(weights.get(name, 0.0))
                if w:
                    total += _weighted_hook_gradient(func, G, n, nd, w)
            results.append((n, total))

    _require_finite_pressure([value for _, value in results])
    for node, value in results:
        set_dnfr(G, node, float(value))

    _write_dnfr_metadata(G, weights=weights, hook_name=hook_name, note=note)


# --- Example hooks (optional) ---


class _PhaseGradient:
    """Callable computing the phase contribution using cached trig values."""

    __slots__ = ("cos", "sin", "theta")

    def __init__(
        self,
        cos_map: Mapping[NodeId, float],
        sin_map: Mapping[NodeId, float],
        theta_map: Mapping[NodeId, float],
    ) -> None:
        self.cos: Mapping[NodeId, float] = cos_map
        self.sin: Mapping[NodeId, float] = sin_map
        self.theta: Mapping[NodeId, float] = theta_map

    def __call__(
        self,
        G: TNFRGraph,
        n: NodeId,
        nd: Mapping[str, Any],
    ) -> float:
        th_i = self.theta[n]
        neighbors = list(G.neighbors(n))
        if len(neighbors) == 2:
            certified = certified_two_neighbor_phase(
                th_i,
                self.theta[neighbors[0]],
                self.theta[neighbors[1]],
            )
            if certified is not None:
                return certified.delta / math.pi
        x = sum(self.cos[v] for v in neighbors)
        y = sum(self.sin[v] for v in neighbors)
        if x == 0.0 and y == 0.0:
            return 0.0
        return angle_diff(math.atan2(y, x), _pressure_phase_coordinate(th_i)) / math.pi


class _NeighborAverageGradient:
    """Callable computing neighbour averages for scalar attributes."""

    __slots__ = ("alias", "values")

    def __init__(
        self,
        alias: tuple[str, ...],
        values: MutableMapping[NodeId, float],
    ) -> None:
        self.alias: tuple[str, ...] = alias
        self.values: MutableMapping[NodeId, float] = values

    def __call__(
        self,
        G: TNFRGraph,
        n: NodeId,
        nd: Mapping[str, Any],
        *,
        coefficient: float = 1.0,
    ) -> float:
        val = self.values.get(n)
        if val is None:
            val = (
                _read_scalar_epi(nd) if self.alias == ALIAS_EPI else _read_capacity(nd)
            )
            self.values[n] = val
        neighbors = list(G.neighbors(n))
        if not neighbors:
            return 0.0
        values = []
        for neigh in neighbors:
            neigh_val = self.values.get(neigh)
            if neigh_val is None:
                neigh_val = (
                    _read_scalar_epi(G.nodes[neigh])
                    if self.alias == ALIAS_EPI
                    else _read_capacity(G.nodes[neigh], val)
                )
                self.values[neigh] = neigh_val
            values.append(neigh_val)
        return mean_neighbor_difference(val, values, coefficient=coefficient)


def dnfr_phase_only(G: TNFRGraph, *, n_jobs: int | None = None) -> None:
    """Compute the circular neighbor-argument phase pressure alone.

    This is not the pairwise sine-coupling law. It uses the default pressure's
    local midpoint certificate and zero represented-resultant extension.

    Parameters
    ----------
    G : nx.Graph
        Graph whose nodes receive the ΔNFR assignment.
    n_jobs : int | None, optional
        Parallel worker hint used when NumPy is unavailable. Defaults to
        serial execution.
    """

    trig = compute_theta_trig((n, _read_phase(nd)) for n, nd in G.nodes(data=True))
    g_phase = _PhaseGradient(trig.cos, trig.sin, trig.theta)
    _apply_dnfr_hook(
        G,
        {"phase": g_phase},
        weights={"phase": 1.0},
        hook_name="dnfr_phase_only",
        note="Example hook.",
        n_jobs=n_jobs,
    )


def dnfr_epi_vf_mixed(G: TNFRGraph, *, n_jobs: int | None = None) -> None:
    """Mix unweighted support EPI and capacity gradients with fixed halves.

    This optional model ignores edge conductances; it is not the default
    weighted EPI channel with phase disabled.

    Parameters
    ----------
    G : nx.Graph
        Graph whose nodes receive the ΔNFR assignment.
    n_jobs : int | None, optional
        Parallel worker hint used when NumPy is unavailable. Defaults to
        serial execution.
    """

    epi_values = {n: _read_scalar_epi(nd) for n, nd in G.nodes(data=True)}
    vf_values = {n: _read_capacity(nd) for n, nd in G.nodes(data=True)}
    grads = {
        "epi": _NeighborAverageGradient(ALIAS_EPI, epi_values),
        "vf": _NeighborAverageGradient(ALIAS_VF, vf_values),
    }
    _apply_dnfr_hook(
        G,
        grads,
        weights={"phase": 0.0, "epi": 0.5, "vf": 0.5},
        hook_name="dnfr_epi_vf_mixed",
        note="Unweighted unique-support gradients; fixed half coefficients.",
        n_jobs=n_jobs,
    )


def dnfr_laplacian(G: TNFRGraph, *, n_jobs: int | None = None) -> None:
    """Unweighted support Laplacian over EPI and capacity with raw coefficients.

    The public EPI/capacity coefficients are not normalized here. Unlike the
    default EPI channel, this optional model ignores edge conductances.

    Parameters
    ----------
    G : nx.Graph
        Graph whose nodes receive the ΔNFR assignment.
    n_jobs : int | None, optional
        Parallel worker hint used when NumPy is unavailable. Defaults to
        serial execution.
    """

    weights_cfg = get_param(G, "DNFR_WEIGHTS")
    wE = float(weights_cfg.get("epi", DEFAULTS["DNFR_WEIGHTS"]["epi"]))
    wV = float(weights_cfg.get("vf", DEFAULTS["DNFR_WEIGHTS"]["vf"]))

    epi_values = {n: _read_scalar_epi(nd) for n, nd in G.nodes(data=True)}
    vf_values = {n: _read_capacity(nd) for n, nd in G.nodes(data=True)}
    grads = {
        "epi": _NeighborAverageGradient(ALIAS_EPI, epi_values),
        "vf": _NeighborAverageGradient(ALIAS_VF, vf_values),
    }
    _apply_dnfr_hook(
        G,
        grads,
        weights={"epi": wE, "vf": wV},
        hook_name="dnfr_laplacian",
        note="Unweighted unique-support gradients; raw EPI/capacity coefficients.",
        n_jobs=n_jobs,
    )


def compute_delta_nfr_hamiltonian(
    G: TNFRGraph,
    *,
    hbar_str: float | None = None,
    cache_hamiltonian: bool = True,
    profile: MutableMapping[str, Any] | None = None,
) -> None:
    r"""Write the auxiliary Hamiltonian's identically zero node read-out.

    This opt-in compatibility hook constructs the auxiliary matrix
    H_int = H_coh + H_freq + H_coupling and reads:

    .. math::
        \Delta\text{NFR}_n = \frac{i}{\hbar_{str}} \langle n | [\hat{H}_{int}, \rho_n] | n \rangle

    For the localized projector rho_n, the diagonal commutator is exactly
    H[n,n] - H[n,n] = 0. Selecting this hook therefore replaces every stored
    pressure by zero. It does not derive the default pressure, evolve a density
    matrix, or establish equilibrium of the full engine with other writers.

    Parameters
    ----------
    G : TNFRGraph
        Graph with nodes containing 'nu_f', 'phase', 'epi', 'si' attributes
    hbar_str : float, optional
        Finite nonzero real scale of the auxiliary matrix model. If None,
        uses ``G.graph.get('HBAR_STR', 1.0)``.
    cache_hamiltonian : bool, default=True
        Retain the latest construction in ``G.graph['_hamiltonian_cache']``.
        Every invocation rebuilds from current attributes, support and
        configuration; the retained object is not reused. False removes any
        previously retained snapshot after successful construction.
    profile : MutableMapping[str, float] or None, optional
        Mutable mapping that accumulates wall-clock timings:

        - ``"hamiltonian_construction"``: Time to build H_int
        - ``"hamiltonian_computation"``: Time to compute all ΔNFR values
        - ``"hamiltonian_write"``: Time to write results to nodes

    Notes
    -----

    The auxiliary affinity is distinct from canonical C(t). Matrix construction
    is dense and uses quadratic storage. The null read-out requires no matrix
    commutator products. Node-count/checksum cache reuse is insufficient because
    the matrix also consumes graph attributes and configuration.

    Examples
    --------

    **Basic usage**:

    >>> import networkx as nx
    >>> from tnfr.dynamics.dnfr import compute_delta_nfr_hamiltonian
    >>> G = nx.cycle_graph(10)
    >>> for node in G.nodes:
    ...     G.nodes[node].update({
    ...         'nu_f': 1.0, 'phase': 0.0, 'epi': 1.0, 'si': 0.8
    ...     })
    >>> compute_delta_nfr_hamiltonian(G)
    >>> # Every stored node pressure is now zero.

    **With profiling**:

    >>> profile = {}
    >>> compute_delta_nfr_hamiltonian(G, profile=profile)
    >>> print(f"Construction: {profile['hamiltonian_construction']:.3f}s")
    >>> print(f"Computation: {profile['hamiltonian_computation']:.3f}s")

    **Integration with dynamics**:

    >>> from tnfr.dynamics import set_delta_nfr_hook
    >>> set_delta_nfr_hook(G, compute_delta_nfr_hamiltonian, name="hamiltonian")
    >>> # Pressure refreshes now select the auxiliary null read-out.

    See Also
    --------
    tnfr.operators.hamiltonian.InternalHamiltonian : Core Hamiltonian class
    default_compute_delta_nfr : Configured multichannel graph pressure
    set_delta_nfr_hook : Register custom ΔNFR computation

    References
    ----------

    ``theory/STRUCTURAL_STABILITY_AND_DYNAMICS.md``, section 5
    """
    from ..operators.hamiltonian import InternalHamiltonian

    # Initialize profiling
    start_timer, stop_timer = _profile_start_stop(
        profile,
        keys=(
            "hamiltonian_construction",
            "hamiltonian_computation",
            "hamiltonian_write",
        ),
    )

    # Resolve the auxiliary model's scale.
    if hbar_str is None:
        hbar_str = G.graph.get("HBAR_STR", 1.0)

    # Rebuild: a node count/checksum cannot certify all consumed state/config.
    cache_key = "_hamiltonian_cache"
    timer = start_timer()
    ham = InternalHamiltonian(G, hbar_str=hbar_str)
    if cache_hamiltonian:
        G.graph[cache_key] = ham
    else:
        G.graph.pop(cache_key, None)
    stop_timer("hamiltonian_construction", timer)

    # Compute ΔNFR for all nodes
    timer = start_timer()

    delta_nfr_values = {}
    for node in ham.nodes:
        delta_nfr = ham.compute_node_delta_nfr(node)
        delta_nfr_values[node] = delta_nfr

    stop_timer("hamiltonian_computation", timer)

    # Write results to graph nodes
    timer = start_timer()

    for node, delta_val in delta_nfr_values.items():
        set_dnfr(G, node, delta_val)

    stop_timer("hamiltonian_write", timer)

    # Write metadata
    _write_dnfr_metadata(
        G,
        weights={"hamiltonian": 1.0},
        hook_name="compute_delta_nfr_hamiltonian",
        note="Auxiliary localized-projector commutator: identically zero pressure",
    )