Coupling

The coupling subsystem builds, adapts, and analyses the inter-oscillator coupling matrix K_nm — the central object in Kuramoto dynamics. K_ij determines how strongly oscillator j pulls oscillator i toward synchrony.

The subsystem spans 27 source files: public API modules for construction (knm), geometry constraints, phase lag estimation, template management, Hodge decomposition, spectral analysis, plasticity, transfer-entropy adaptation, causal inference, connectome generation, E/I balance, attention residuals, spatial modulation, and a universal Bayesian prior, plus validated backend bridge files.

Pipeline position

CouplingBuilder.build() ──→ K_nm, α ──→ UPDEEngine.step()
       ↑                                       │
  UniversalPrior                                ↓
  LagModel.estimate ────→ α            compute_order_parameter()
  connectome loader ─────→ K_nm                 │
  auto-coupling-estimation ← raw phase time series
  plasticity/TE ←────────────────── phase history

CouplingBuilder is the entry point of the SPO pipeline. Every engine variant consumes (phases, omegas, knm, zeta, psi, alpha), so the coupling matrix and phase-lag matrix are required for any simulation. For data-first onboarding, auto_coupling_estimation() infers an initial directed coupling graph from phase time series before review, projection, or engine execution. The inference boundary requires finite real phase samples and enforces the transfer-entropy invariant that directed scores are non-negative with no self-edge diagonal. Across the coupling public boundary, boolean aliases mean Python bool, NumPy boolean scalars, and object arrays containing either form; those inputs are rejected before any float coercion.

---

K_nm Construction

CouplingBuilder

Builds coupling matrices from parameters.

Methods:

Method	Signature	Description
build	(n_layers, base_strength, decay_alpha) → CouplingState	Exponential-decay K_nm
build_scpn_physics	(k_base=0.45, alpha_decay=0.3) → CouplingState	16-layer SCPN physics
build_with_amplitude	(n, base, decay, amp_str, amp_dec) → CouplingState	Phase + amplitude K
apply_handshakes	(state, path) → CouplingState	Overlay from JSON spec
switch_template	(state, name, templates) → CouplingState	Runtime topology switch

apply_handshakes() parses the JSON specification fail-closed: non-finite constants, duplicate object keys, non-list matrix payloads, self-coupled entries, and out-of-range layer indices are rejected before any K_nm entries are modified.

CouplingState (frozen dataclass)

Field	Type	Description
knm	NDArray	Phase coupling matrix K_ij
alpha	NDArray	Phase-lag matrix α_ij
active_template	str	Name of active template
knm_r	NDArray \| None	Amplitude coupling (Stuart-Landau)

Coupling equation

For the standard Kuramoto model, the coupling enters as:

dθ_i/dt = ω_i + Σ_j K_ij sin(θ_j - θ_i - α_ij) + ζ sin(Ψ - θ_i)

K_ij is the (i,j) entry of the coupling matrix. The matrix must satisfy:

1. Square: K ∈ R^{N×N}
2. Symmetric: K_ij = K_ji (undirected coupling; directed via asymmetric K)
3. Non-negative: K_ij ≥ 0
4. Zero diagonal: K_ii = 0 (no self-coupling)

Exponential-decay construction

CouplingBuilder.build(n, base_strength, decay_alpha) produces:

K_ij = base_strength × exp(-decay_alpha × |i - j|),  K_ii = 0

This generates nearest-neighbour-dominant coupling with exponential fall-off — appropriate for layered systems where adjacent layers interact more strongly than distant ones.

SCPN physics construction

build_scpn_physics(k_base=0.45, alpha_decay=0.3) produces a 16×16 matrix using three coupling mechanisms:

1. Adjacent layers (|i-j| = 1): timescale matching via SCPN_LAYER_TIMESCALES (Quantum: 1e-15s to Social: 3.15e7s)
2. Near-neighbour (|i-j| ≤ 3): geometric mean of adjacent couplings
3. Distant (|i-j| > 3): exponential decay from k_base

The 16 SCPN layers span 22 orders of magnitude in timescale:

Layer	Name	Timescale
L1	Quantum	1e-15 s
L2	Sub-nuclear	1e-12 s
L3	Atomic	1e-10 s
L4	Molecular	1e-9 s
L5	Cellular	1e-3 s
L6	Neural	1e-2 s
L7	Synaptic	1e-1 s
L8	Circuit	1 s
L9	Regional	10 s
L10	Behavioural	60 s
L11	Cognitive	600 s
L12	Social	3600 s
L13	Cultural	86400 s
L14	Evolutionary	3.15e6 s
L15	Cosmological	3.15e7 s
L16	Director (meta)	—

Performance: build(100) < 10 ms, build_scpn_physics() < 5 ms.

knm

Coupling-matrix builders for generic and SCPN-layer topologies.

CouplingBuilder constructs deterministic K_nm, alpha, and optional amplitude-coupling snapshots from validated scalar parameters. The Rust-backed path and NumPy fallback share the same public contract: finite non-boolean inputs, positive layer counts, zero diagonals, and explicit template labels for runtime/audit reporting.

Classes

CouplingState dataclass

CouplingState(
    knm: FloatArray,
    alpha: FloatArray,
    active_template: str,
    knm_r: FloatArray | None = None,
)

Immutable snapshot of phase/amplitude coupling matrices and template.

CouplingBuilder

Builds Knm coupling matrices.

Methods:

build

build(
    n_layers: int, base_strength: float, decay_alpha: float
) -> CouplingState

Build an exponentially decayed phase-coupling matrix.

Parameters

n_layers Number of hierarchy layers or oscillators represented in the square coupling matrix. base_strength Coupling strength before distance decay is applied. decay_alpha Non-negative exponential decay coefficient in exp(-decay_alpha * |i - j|).

Returns

CouplingState Coupling snapshot with knm and alpha matrices of shape (n_layers, n_layers). The diagonal of knm is zero and alpha is initialised to zeros.

Notes

When the Rust extension is available, construction dispatches to spo_kernel.PyCouplingBuilder and preserves the same output contract as the NumPy fallback.

Raises

ValueError If the layer count or coupling parameters are invalid.

Source code in src/scpn_phase_orchestrator/coupling/knm.py

def build(
    self, n_layers: int, base_strength: float, decay_alpha: float
) -> CouplingState:
    """Build an exponentially decayed phase-coupling matrix.

    Parameters
    ----------
    n_layers
        Number of hierarchy layers or oscillators represented in the
        square coupling matrix.
    base_strength
        Coupling strength before distance decay is applied.
    decay_alpha
        Non-negative exponential decay coefficient in
        ``exp(-decay_alpha * |i - j|)``.

    Returns
    -------
    CouplingState
        Coupling snapshot with ``knm`` and ``alpha`` matrices of shape
        ``(n_layers, n_layers)``. The diagonal of ``knm`` is zero and
        ``alpha`` is initialised to zeros.

    Notes
    -----
    When the Rust extension is available, construction dispatches to
    ``spo_kernel.PyCouplingBuilder`` and preserves the same output
    contract as the NumPy fallback.

    Raises
    ------
    ValueError
        If the layer count or coupling parameters are invalid.
    """
    n_layers = _validate_positive_int(n_layers, name="n_layers")
    base_strength = _validate_finite_float(
        base_strength,
        name="base_strength",
        lower_bound=0.0,
    )
    decay_alpha = _validate_finite_float(
        decay_alpha,
        name="decay_alpha",
        lower_bound=0.0,
    )
    if _HAS_RUST:  # pragma: no cover
        from spo_kernel import PyCouplingBuilder

        try:
            d = PyCouplingBuilder().build(n_layers, base_strength, decay_alpha)
            n = _validate_positive_int(d["n"], name="rust n_layers")
            if n != n_layers:
                raise ValueError("Rust coupling output layer count mismatch")
            rust_knm, rust_alpha = _validate_coupling_output(
                d["knm"], d["alpha"], n_layers=n_layers
            )
            return CouplingState(
                knm=rust_knm, alpha=rust_alpha, active_template="default"
            )
        except Exception as exc:
            _fallback_reason = exc
    idx = np.arange(n_layers)
    dist = np.abs(idx[:, np.newaxis] - idx[np.newaxis, :])
    knm = base_strength * np.exp(-decay_alpha * dist)
    np.fill_diagonal(knm, 0.0)
    alpha = np.zeros((n_layers, n_layers), dtype=np.float64)
    return CouplingState(knm=knm, alpha=alpha, active_template="default")

build_scpn_physics

build_scpn_physics(
    k_base: float = 0.45, alpha_decay: float = 0.3
) -> CouplingState

Build 16×16 K_nm using SCPN layer physics.

Three coupling mechanisms (Paper 0, HolonomicAtlas v2.4.0): - Adjacent: timescale matching with calibration anchors - Near-neighbor (|n-m|=2): geometric mean of intermediate path - Distant (|n-m|>=3): exponential decay with cross-hierarchy boosts

Returns CouplingState with 16×16 matrix.

Parameters

k_base : float Base coupling strength. alpha_decay : float Exponential decay rate of the coupling.

Returns

CouplingState The 16×16 coupling state from SCPN layer physics.

Source code in src/scpn_phase_orchestrator/coupling/knm.py

def build_scpn_physics(
    self,
    k_base: float = 0.45,
    alpha_decay: float = 0.3,
) -> CouplingState:
    """Build 16×16 K_nm using SCPN layer physics.

    Three coupling mechanisms (Paper 0, HolonomicAtlas v2.4.0):
    - Adjacent: timescale matching with calibration anchors
    - Near-neighbor (|n-m|=2): geometric mean of intermediate path
    - Distant (|n-m|>=3): exponential decay with cross-hierarchy boosts

    Returns CouplingState with 16×16 matrix.

    Parameters
    ----------
    k_base : float
        Base coupling strength.
    alpha_decay : float
        Exponential decay rate of the coupling.

    Returns
    -------
    CouplingState
        The 16×16 coupling state from SCPN layer physics.
    """
    k_base = _validate_finite_float(
        k_base,
        name="k_base",
        lower_bound=0.0,
        inclusive=False,
    )
    alpha_decay = _validate_finite_float(
        alpha_decay,
        name="alpha_decay",
        lower_bound=0.0,
    )
    K = np.zeros((16, 16))

    # Pass 1: Adjacent layers. Use anchors where available.
    for n in range(1, 16):
        m = n + 1
        if (n, m) in SCPN_CALIBRATION_ANCHORS:
            val = SCPN_CALIBRATION_ANCHORS[(n, m)]
        else:
            val = self._adjacent_coupling(n, m, k_base)
        K[n - 1, m - 1] = val
        K[m - 1, n - 1] = val

    # Pass 2: Near-neighbor (|n-m|=2), geometric mean of intermediate path
    for n in range(1, 15):
        m = n + 2
        mid = (n + m) // 2
        k1 = K[n - 1, mid - 1]
        k2 = K[mid - 1, m - 1]
        val = float(np.sqrt(k1 * k2))
        # Frequency penalty for large timescale mismatch
        tau_n = SCPN_LAYER_TIMESCALES.get(n, 1.0)
        tau_m = SCPN_LAYER_TIMESCALES.get(m, 1.0)
        if n != 16 and m != 16 and tau_n > 0 and tau_m > 0:
            omega_n = 2.0 * np.pi / tau_n
            omega_m = 2.0 * np.pi / tau_m
            omega_avg = (omega_n + omega_m) / 2.0
            penalty = 1.0 + abs(omega_n - omega_m) / omega_avg * 0.1
            val /= penalty
        val = float(np.clip(val, 0.01, 0.4))
        K[n - 1, m - 1] = val
        K[m - 1, n - 1] = val

    # Pass 3: Distant (|n-m|>=3), exponential decay
    for n in range(1, 17):
        for m in range(n + 3, 17):
            dist = abs(n - m)
            val = k_base * np.exp(-alpha_decay * dist)
            val = float(np.clip(val, 0.001, 0.2))
            K[n - 1, m - 1] = val
            K[m - 1, n - 1] = val

    # Cross-hierarchy boosts (Paper 0, Section 5)
    _set_symmetric(K, 1, 16, max(K[0, 15], 0.05))  # Quantum-Meta
    _set_symmetric(K, 5, 7, max(K[4, 6], 0.15))  # Psycho-Symbolic

    alpha = np.zeros((16, 16), dtype=np.float64)
    return CouplingState(knm=K, alpha=alpha, active_template="scpn_physics")

apply_handshakes

apply_handshakes(
    state: CouplingState, handshakes_path: str | Path
) -> CouplingState

Overlay documented inter-layer couplings from JSON spec.

Reads the KNM_MATRIX_COMPLETE_SPECIFICATION.json format: each entry has from_layer, to_layer, coupling_strength. Negative values are preserved (inhibitory coupling).

Parameters

state : CouplingState The coupling state to transform. handshakes_path : str | Path Path to the JSON inter-layer handshake spec.

Returns

CouplingState The coupling state with the documented inter-layer couplings overlaid.

Raises

ValueError If the handshake spec is malformed or out of range.

Source code in src/scpn_phase_orchestrator/coupling/knm.py

def apply_handshakes(
    self, state: CouplingState, handshakes_path: str | Path
) -> CouplingState:
    """Overlay documented inter-layer couplings from JSON spec.

    Reads the KNM_MATRIX_COMPLETE_SPECIFICATION.json format:
    each entry has from_layer, to_layer, coupling_strength.
    Negative values are preserved (inhibitory coupling).

    Parameters
    ----------
    state : CouplingState
        The coupling state to transform.
    handshakes_path : str | Path
        Path to the JSON inter-layer handshake spec.

    Returns
    -------
    CouplingState
        The coupling state with the documented inter-layer couplings overlaid.

    Raises
    ------
    ValueError
        If the handshake spec is malformed or out of range.
    """
    path = Path(handshakes_path)
    data = _loads_knm_json(path.read_text(encoding="utf-8"))
    matrix = data.get("matrix")
    if not isinstance(matrix, list):
        raise ValueError("handshake matrix must be a list")
    knm = state.knm.copy()
    n_layers = knm.shape[0]
    for idx, entry in enumerate(matrix):
        if not isinstance(entry, dict):
            raise ValueError(f"matrix[{idx}] must be a mapping")
        fr = _validate_layer_index(
            entry.get("from_layer"),
            name="from_layer",
            n_layers=n_layers,
        )
        to = _validate_layer_index(
            entry.get("to_layer"),
            name="to_layer",
            n_layers=n_layers,
        )
        if fr == to:
            raise ValueError("handshake self-coupling entries are not physical")
        strength = _validate_finite_float(
            entry.get("coupling_strength"),
            name="coupling_strength",
        )
        knm[fr - 1, to - 1] = strength
        # Symmetric unless negative (directional inhibition)
        if strength >= 0:
            knm[to - 1, fr - 1] = strength
    return CouplingState(
        knm=knm,
        alpha=state.alpha.copy(),
        active_template="scpn_handshakes",
        knm_r=state.knm_r,
    )

build_with_amplitude

build_with_amplitude(
    n_layers: int,
    base_strength: float,
    decay_alpha: float,
    amp_strength: float,
    amp_decay: float,
) -> CouplingState

Build phase + amplitude coupling matrices together.

Parameters

n_layers : int Number of SCPN hierarchy layers. base_strength : float Base coupling strength at zero layer separation. decay_alpha : float Exponential decay rate of the coupling across layer separation. amp_strength : float Base amplitude-coupling strength. amp_decay : float Exponential decay rate of the amplitude coupling.

Returns

CouplingState The coupling state with both phase and amplitude coupling matrices.

Source code in src/scpn_phase_orchestrator/coupling/knm.py

def build_with_amplitude(
    self,
    n_layers: int,
    base_strength: float,
    decay_alpha: float,
    amp_strength: float,
    amp_decay: float,
) -> CouplingState:
    """Build phase + amplitude coupling matrices together.

    Parameters
    ----------
    n_layers : int
        Number of SCPN hierarchy layers.
    base_strength : float
        Base coupling strength at zero layer separation.
    decay_alpha : float
        Exponential decay rate of the coupling across layer separation.
    amp_strength : float
        Base amplitude-coupling strength.
    amp_decay : float
        Exponential decay rate of the amplitude coupling.

    Returns
    -------
    CouplingState
        The coupling state with both phase and amplitude coupling matrices.
    """
    amp_strength = _validate_finite_float(
        amp_strength,
        name="amp_strength",
        lower_bound=0.0,
    )
    amp_decay = _validate_finite_float(
        amp_decay,
        name="amp_decay",
        lower_bound=0.0,
    )
    phase = self.build(n_layers, base_strength, decay_alpha)
    idx = np.arange(n_layers)
    dist = np.abs(idx[:, np.newaxis] - idx[np.newaxis, :])
    knm_r = amp_strength * np.exp(-amp_decay * dist)
    np.fill_diagonal(knm_r, 0.0)
    return CouplingState(
        knm=phase.knm,
        alpha=phase.alpha,
        active_template=phase.active_template,
        knm_r=knm_r,
    )

switch_template

switch_template(
    state: CouplingState,
    template_name: str,
    templates: dict[str, FloatArray],
) -> CouplingState

Replace the active K_nm with a named template matrix.

Parameters

state : CouplingState The coupling state to transform. template_name : str Name of the template to activate. templates : dict[str, FloatArray] Mapping of template name to coupling matrix.

Returns

CouplingState The coupling state with the named template installed as K_nm.

Raises

KeyError If template_name is not in templates. ValueError If the template matrix is invalid.

Source code in src/scpn_phase_orchestrator/coupling/knm.py

def switch_template(
    self,
    state: CouplingState,
    template_name: str,
    templates: dict[str, FloatArray],
) -> CouplingState:
    """Replace the active K_nm with a named template matrix.

    Parameters
    ----------
    state : CouplingState
        The coupling state to transform.
    template_name : str
        Name of the template to activate.
    templates : dict[str, FloatArray]
        Mapping of template name to coupling matrix.

    Returns
    -------
    CouplingState
        The coupling state with the named template installed as ``K_nm``.

    Raises
    ------
    KeyError
        If ``template_name`` is not in ``templates``.
    ValueError
        If the template matrix is invalid.
    """
    if template_name not in templates:
        raise KeyError(f"Template {template_name!r} not found")
    template = np.asarray(templates[template_name], dtype=np.float64)
    if template.shape != state.knm.shape:
        raise ValueError(
            f"template shape {template.shape}, expected {state.knm.shape}"
        )
    if not np.all(np.isfinite(template)):
        raise ValueError("template values must be finite")
    if not np.allclose(np.diag(template), 0.0, rtol=0.0, atol=1e-15):
        raise ValueError("template self-coupling diagonal must be zero")
    return CouplingState(
        knm=template.copy(),
        alpha=state.alpha.copy(),
        active_template=template_name,
        knm_r=state.knm_r,
    )

---

Geometry Constraints

Enforces structural invariants on K_nm.

Constraint classes

Class	project(knm) behaviour
SymmetryConstraint	Returns (K + K^T) / 2
NonNegativeConstraint	Clamps negative entries to 0

Validation

validate_knm(knm, atol=1e-12) accepts only finite real square matrices and checks all four invariants: symmetric, non-negative, zero diagonal, and boolean/complex aliases rejected before numeric projection. Raises ValueError on violation.

project_knm(knm, constraints) applies constraints sequentially, then zeros the diagonal. Built-in and custom constraints are fail-closed: each constraint must be a GeometryConstraint, preserve the matrix shape, and return finite real square K_nm values before the next projection step.

geometry_constraints

Projection and validation helpers for coupling-matrix geometry.

Geometry constraints project candidate K_nm matrices onto simple feasible sets such as symmetry and non-negativity. validate_knm enforces the runtime matrix contract used by domainpacks and UPDE handoff: square, symmetric, non-negative, and zero diagonal within tolerance.

Classes

GeometryConstraint

Bases: ABC

Base class for K_nm matrix geometry constraints.

Methods:

project abstractmethod

project(knm: FloatArray) -> FloatArray

Project knm onto the feasible set defined by this constraint.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

FloatArray The projection of knm onto the constraint's feasible set.

Source code in src/scpn_phase_orchestrator/coupling/geometry_constraints.py

@abstractmethod
def project(self, knm: FloatArray) -> FloatArray:
    """Project *knm* onto the feasible set defined by this constraint.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    FloatArray
        The projection of *knm* onto the constraint's feasible set.
    """
    ...

SymmetryConstraint

Bases: GeometryConstraint

Enforce K_nm symmetry: K -> (K + K^T) / 2.

Methods:

project

project(knm: FloatArray) -> FloatArray

Return the symmetric part of knm.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

FloatArray The symmetric part of knm.

Source code in src/scpn_phase_orchestrator/coupling/geometry_constraints.py

def project(self, knm: FloatArray) -> FloatArray:
    """Return the symmetric part of *knm*.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    FloatArray
        The symmetric part of *knm*.
    """
    knm = _validate_knm_matrix(knm)
    result: FloatArray = 0.5 * (knm + knm.T)
    return result

NonNegativeConstraint

Bases: GeometryConstraint

Clamp negative entries to zero.

Methods:

project

project(knm: FloatArray) -> FloatArray

Return knm with all negative entries replaced by 0.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

FloatArray knm with negative entries clipped to zero.

Source code in src/scpn_phase_orchestrator/coupling/geometry_constraints.py

def project(self, knm: FloatArray) -> FloatArray:
    """Return *knm* with all negative entries replaced by 0.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    FloatArray
        *knm* with negative entries clipped to zero.
    """
    knm = _validate_knm_matrix(knm)
    result: FloatArray = np.maximum(knm, 0.0)
    return result

Functions:

validate_knm

validate_knm(
    knm: FloatArray, *, atol: float = 1e-12
) -> None

Check that a coupling matrix is square, symmetric, non-negative, zero-diagonal.

Raises ValueError with a specific message on the first violation found.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N). atol : float Absolute tolerance for the validity checks.

Raises

ValueError If knm is not square, symmetric, non-negative, and zero-diagonal.

Source code in src/scpn_phase_orchestrator/coupling/geometry_constraints.py

def validate_knm(knm: FloatArray, *, atol: float = 1e-12) -> None:
    """Check that a coupling matrix is square, symmetric, non-negative, zero-diagonal.

    Raises ValueError with a specific message on the first violation found.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    atol : float
        Absolute tolerance for the validity checks.

    Raises
    ------
    ValueError
        If ``knm`` is not square, symmetric, non-negative, and zero-diagonal.
    """
    knm = _validate_knm_matrix(knm)
    if not np.allclose(knm, knm.T, atol=atol):
        raise ValueError("Knm is not symmetric")
    if np.any(knm < -atol):
        raise ValueError("Knm contains negative entries")
    diag_max = float(np.max(np.abs(np.diag(knm))))
    if diag_max > atol:
        raise ValueError(f"Knm diagonal is non-zero (max |diag| = {diag_max:.2e})")

project_knm

project_knm(
    knm: FloatArray, constraints: list[GeometryConstraint]
) -> FloatArray

Apply all geometry constraints sequentially, then zero the diagonal.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N). constraints : list[GeometryConstraint] Geometry constraints applied in sequence.

Returns

FloatArray The coupling matrix after applying every constraint and zeroing the diagonal.

Raises

ValueError If a constraint produces an invalid coupling matrix.

Source code in src/scpn_phase_orchestrator/coupling/geometry_constraints.py

def project_knm(knm: FloatArray, constraints: list[GeometryConstraint]) -> FloatArray:
    """Apply all geometry constraints sequentially, then zero the diagonal.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    constraints : list[GeometryConstraint]
        Geometry constraints applied in sequence.

    Returns
    -------
    FloatArray
        The coupling matrix after applying every constraint and zeroing the diagonal.

    Raises
    ------
    ValueError
        If a constraint produces an invalid coupling matrix.
    """
    result = _validate_knm_matrix(knm).copy()
    for c in constraints:
        if not isinstance(c, GeometryConstraint):
            raise ValueError("geometry constraint must be a GeometryConstraint")
        projected = _validate_knm_matrix(
            c.project(result), name="geometry constraint output"
        )
        if projected.shape != result.shape:
            raise ValueError("geometry constraint output shape must match Knm shape")
        result = projected
    np.fill_diagonal(result, 0.0)
    return result

---

Phase Lag Estimation

Estimates inter-oscillator phase lags α_ij from observed time series or known physical distances.

From distances

LagModel.estimate_from_distances(distances, speed) computes:

α_ij = 2π × distances[i,j] / speed

Inputs must be a finite real square physical-distance matrix with non-negative entries, a zero diagonal, and symmetric pair distances, plus a finite positive propagation speed. Boolean aliases and complex/object-complex distance payloads are rejected before numeric coercion because transport delays are ordered real quantities. Returns an antisymmetric matrix: α_ij = -α_ji. This encodes the fact that if signal from i reaches j with positive lag, then j reaches i with negative lag. Directed or asymmetric empirical delays belong in build_alpha_matrix, not in the physical-distance constructor.

From cross-correlation

LagModel().estimate_lag(signal_a, signal_b, sample_rate) finds the cross-correlation peak lag in seconds between two signals. Signals must be finite real one-dimensional arrays with equal non-zero length and non-zero variance. The sample-rate must be a finite positive real value. Constant, boolean, complex/object-complex, non-finite, or length-mismatched signals are rejected before cross-correlation because they do not define a reliable phase-lag estimate.

Matrix construction

build_alpha_matrix(lag_estimates, n_layers, carrier_freq_hz=1.0) converts pairwise lag estimates (in seconds) to a phase-offset matrix (in radians):

α_ij = 2π × carrier_freq_hz × lag_seconds_ij

Performance: estimate_from_distances(64×64) < 5 ms.

lags

Phase-lag estimation and alpha-matrix construction.

LagModel converts physical distances or observed signal offsets into antisymmetric phase-lag matrices consumed by the UPDE engine. The module is kept dependency-light and deterministic while rejecting non-physical distance, sample, carrier, and speed inputs before the resulting lags enter closed-loop runs.

Classes

LagModel

Phase-lag estimation and alpha matrix construction.

Methods:

estimate_from_distances staticmethod

estimate_from_distances(
    distances: FloatArray, speed: float
) -> FloatArray

Build antisymmetric alpha matrix from pairwise distances and speed.

alpha[i,j] = 2*pi * distances[i,j] / speed. Matches the Rust LagModel::estimate_from_distances algorithm.

Parameters

distances : FloatArray Pairwise distance matrix, shape (N, N). speed : float Signal propagation speed.

Returns

FloatArray The antisymmetric phase-lag matrix, shape (N, N).

Source code in src/scpn_phase_orchestrator/coupling/lags.py

@staticmethod
def estimate_from_distances(distances: FloatArray, speed: float) -> FloatArray:
    """Build antisymmetric alpha matrix from pairwise distances and speed.

    alpha[i,j] = 2*pi * distances[i,j] / speed.
    Matches the Rust ``LagModel::estimate_from_distances`` algorithm.

    Parameters
    ----------
    distances : FloatArray
        Pairwise distance matrix, shape ``(N, N)``.
    speed : float
        Signal propagation speed.

    Returns
    -------
    FloatArray
        The antisymmetric phase-lag matrix, shape ``(N, N)``.
    """
    distances = _validate_distances(distances)
    speed = _validate_positive_real(speed, name="speed")
    # Antisymmetric propagation lag: alpha[i, j] = 2*pi*distance[i, j]/speed
    # for i < j and its negative below the diagonal. Vectorised over the strict
    # upper triangle (distances are symmetric with a zero diagonal), which is
    # bit-identical to the per-element form and removes the O(N^2) Python loop.
    scaled = (2.0 * np.pi * distances) / speed
    upper: FloatArray = np.triu(scaled, k=1)
    alpha: FloatArray = upper - upper.T
    return np.ascontiguousarray(alpha, dtype=np.float64)

estimate_lag

estimate_lag(
    signal_a: FloatArray,
    signal_b: FloatArray,
    sample_rate: float,
) -> float

Cross-correlation peak lag in seconds between two signals.

Parameters

signal_a : FloatArray First signal, shape (T,). signal_b : FloatArray Second signal, shape (T,). sample_rate : float Sampling rate in Hz.

Returns

float The cross-correlation peak lag in seconds.

Raises

ValueError If the signals differ in length or the sample rate is invalid.

Source code in src/scpn_phase_orchestrator/coupling/lags.py

def estimate_lag(
    self, signal_a: FloatArray, signal_b: FloatArray, sample_rate: float
) -> float:
    """Cross-correlation peak lag in seconds between two signals.

    Parameters
    ----------
    signal_a : FloatArray
        First signal, shape ``(T,)``.
    signal_b : FloatArray
        Second signal, shape ``(T,)``.
    sample_rate : float
        Sampling rate in Hz.

    Returns
    -------
    float
        The cross-correlation peak lag in seconds.

    Raises
    ------
    ValueError
        If the signals differ in length or the sample rate is invalid.
    """
    signal_a = _validate_signal(signal_a, name="signal_a")
    signal_b = _validate_signal(signal_b, name="signal_b")
    if signal_a.shape != signal_b.shape:
        raise ValueError("signals must be the same finite one-dimensional arrays")
    sample_rate = _validate_positive_real(sample_rate, name="sample_rate")
    corr = np.correlate(
        signal_a - signal_a.mean(), signal_b - signal_b.mean(), mode="full"
    )
    peak_idx = np.argmax(corr)
    lag_samples = peak_idx - (len(signal_a) - 1)
    return float(lag_samples / sample_rate)

build_alpha_matrix

build_alpha_matrix(
    lag_estimates: dict[tuple[int, int], float],
    n_layers: int,
    carrier_freq_hz: float = 1.0,
) -> FloatArray

Pairwise lag estimates (seconds) to phase-offset matrix (radians).

alpha[i,j] = 2picarrier_freq_hz*lag[i,j], antisymmetric. carrier_freq_hz defaults to 1.0 for backward compatibility.

Parameters

lag_estimates : dict[tuple[int, int], float] Mapping of oscillator pair to estimated lag in seconds. n_layers : int Number of SCPN hierarchy layers. carrier_freq_hz : float Carrier frequency in Hz used to convert lag to phase.

Returns

FloatArray The phase-offset matrix in radians, shape (N, N).

Source code in src/scpn_phase_orchestrator/coupling/lags.py

def build_alpha_matrix(
    self,
    lag_estimates: dict[tuple[int, int], float],
    n_layers: int,
    carrier_freq_hz: float = 1.0,
) -> FloatArray:
    """Pairwise lag estimates (seconds) to phase-offset matrix (radians).

    alpha[i,j] = 2*pi*carrier_freq_hz*lag[i,j], antisymmetric.
    ``carrier_freq_hz`` defaults to 1.0 for backward compatibility.

    Parameters
    ----------
    lag_estimates : dict[tuple[int, int], float]
        Mapping of oscillator pair to estimated lag in seconds.
    n_layers : int
        Number of SCPN hierarchy layers.
    carrier_freq_hz : float
        Carrier frequency in Hz used to convert lag to phase.

    Returns
    -------
    FloatArray
        The phase-offset matrix in radians, shape ``(N, N)``.
    """
    n_layers = _validate_n_layers(n_layers)
    carrier_freq_hz = _validate_positive_real(
        carrier_freq_hz, name="carrier_freq_hz"
    )
    alpha: FloatArray = np.zeros((n_layers, n_layers), dtype=np.float64)
    for indices, lag in lag_estimates.items():
        i, j, lag = _validate_lag_entry(indices, lag, n_layers=n_layers)
        offset = 2.0 * np.pi * carrier_freq_hz * lag
        alpha[i, j] = offset
        alpha[j, i] = -offset
    return alpha

---

Coupling Templates

Pre-configured coupling topologies for regime-dependent switching.

KnmTemplate (frozen dataclass)

Field	Type	Description
name	str	Template identifier
knm	NDArray	Coupling matrix
alpha	NDArray	Phase-lag matrix
description	str	Human-readable description

KnmTemplateSet

Registry for named templates:

• add(template) — register (overwrites existing with same name)
• get(name) → KnmTemplate — retrieve (raises KeyError if missing, error message lists available names)
• list_names() → list[str] — all registered names

Usage: The supervisor can switch coupling topology at runtime by calling CouplingBuilder.switch_template(state, name, templates) when a regime transition occurs (e.g., switching from all-to-all to nearest-neighbour when entering DEGRADED regime).

templates

Named coupling-template registry for runtime K/alpha switching.

Templates bundle a phase-coupling matrix, phase-lag matrix, description, and stable name. The registry is intentionally in-memory and deterministic: retrieval fails explicitly with available names when a requested template is not registered, leaving persistence and validation to the binding/template owner.

Classes

KnmTemplate dataclass

KnmTemplate(
    name: str,
    knm: FloatArray,
    alpha: FloatArray,
    description: str,
)

Named K_nm coupling matrix with associated phase-lag matrix.

KnmTemplateSet

KnmTemplateSet()

Registry of named K_nm templates for runtime switching.

Source code in src/scpn_phase_orchestrator/coupling/templates.py

def __init__(self) -> None:
    self._templates: dict[str, KnmTemplate] = {}

Methods:

add

add(template: KnmTemplate) -> None

Register a template, overwriting any existing one with the same name.

Parameters

template : KnmTemplate The K_nm template to register.

Raises

TypeError If template is not a KnmTemplate. ValueError If the template is invalid.

Source code in src/scpn_phase_orchestrator/coupling/templates.py

def add(self, template: KnmTemplate) -> None:
    """Register a template, overwriting any existing one with the same name.

    Parameters
    ----------
    template : KnmTemplate
        The ``K_nm`` template to register.

    Raises
    ------
    TypeError
        If ``template`` is not a ``KnmTemplate``.
    ValueError
        If the template is invalid.
    """
    if not isinstance(template, KnmTemplate):
        raise TypeError(f"template must be KnmTemplate, got {template!r}")
    normalized_name = template.name.strip()
    if not normalized_name:
        raise ValueError("template name must be a non-empty string")
    if (
        not isinstance(template.description, str)
        or not template.description.strip()
    ):
        raise ValueError("template description must be a non-empty string")
    if template.knm.ndim != 2 or template.alpha.ndim != 2:
        raise ValueError("template knm/alpha must be 2D matrices")
    if not np.issubdtype(template.knm.dtype, np.floating) or not np.issubdtype(
        template.alpha.dtype, np.floating
    ):
        raise ValueError("template knm/alpha must use floating-point dtypes")
    if template.knm.shape != template.alpha.shape:
        raise ValueError("template knm and alpha must have identical shapes")
    if template.knm.shape[0] != template.knm.shape[1]:
        raise ValueError("template knm must be square")
    if not np.isfinite(template.knm).all() or not np.isfinite(template.alpha).all():
        raise ValueError("template knm/alpha must contain only finite values")
    if any(
        isinstance(v, bool) or not isinstance(v, Real)
        for v in (float(np.min(template.knm)), float(np.max(template.knm)))
    ):
        raise ValueError("template knm must contain numeric real values")
    if not isfinite(float(np.min(template.alpha))) or not isfinite(
        float(np.max(template.alpha))
    ):
        raise ValueError("template alpha must contain finite real values")
    self._templates[normalized_name] = KnmTemplate(
        name=normalized_name,
        knm=_matrix_copy(template.knm),
        alpha=_matrix_copy(template.alpha),
        description=template.description,
    )

get

get(name: str) -> KnmTemplate

Retrieve a template by name. Raises KeyError if not found.

Parameters

name : str Name to look up.

Returns

KnmTemplate The registered template with the given name.

Raises

KeyError If no template with that name is registered.

Source code in src/scpn_phase_orchestrator/coupling/templates.py

def get(self, name: str) -> KnmTemplate:
    """Retrieve a template by name. Raises KeyError if not found.

    Parameters
    ----------
    name : str
        Name to look up.

    Returns
    -------
    KnmTemplate
        The registered template with the given name.

    Raises
    ------
    KeyError
        If no template with that name is registered.
    """
    if not isinstance(name, str) or not name.strip():
        raise KeyError(f"template name must be a non-empty string, got {name!r}")
    normalized = name.strip()
    try:
        template = self._templates[normalized]
        return KnmTemplate(
            name=template.name,
            knm=_matrix_copy(template.knm),
            alpha=_matrix_copy(template.alpha),
            description=template.description,
        )
    except KeyError:
        available = ", ".join(sorted(self._templates)) or "(none)"
        msg = f"Unknown template {name!r}; available: {available}"
        raise KeyError(msg) from None

list_names

list_names() -> list[str]

Return all registered template names.

Returns

list[str] Return all registered template names.

Source code in src/scpn_phase_orchestrator/coupling/templates.py

def list_names(self) -> list[str]:
    """Return all registered template names.

    Returns
    -------
    list[str]
        Return all registered template names.
    """
    return list(self._templates.keys())

---

Combinatorial Hodge Decomposition

Decomposes the Kuramoto coupling current into three L²-orthogonal edge-flow components via combinatorial Hodge theory (Jiang, Lim, Yao & Ye 2011, Statistical ranking and combinatorial Hodge theory, Math. Program. 127 (1):203–244):

coupling current  f = gradient ⊕ curl ⊕ harmonic

The oscillator network is treated as a simplicial complex (V, E, T): vertices are oscillators, edges are the pairs {i, j} with non-zero symmetric coupling, and triangles are the 3-cliques of that graph (or an explicit user-supplied set). The decomposed object is the alternating edge flow

f_ij = ½(K_ij + K_ji) · sin(θ_j − θ_i)

— the canonical coupling current, built from the symmetric coupling part so it satisfies f_ji = −f_ij. With node–edge incidence B1 and edge–triangle incidence B2:

gradient = B1ᵀ · L0⁺ · (B1 f)     # curl-free conservative flow
curl     = B2  · L2⁺ · (B2ᵀ f)    # divergence-free rotational flow
harmonic = f − gradient − curl    # ker of the Hodge 1-Laplacian

where L0 = B1 B1ᵀ and L2 = B2ᵀ B2. Because B1 B2 = 0, the three components are mutually L²-orthogonal.

HodgeResult (dataclass)

Field	Type	Physical meaning
gradient	NDArray (N, N)	Conservative (curl-free) flow grad(s)
curl	NDArray (N, N)	Rotational (divergence-free) flow bounded by triangles
harmonic	NDArray (N, N)	Topological residual in ker(L1) (non-zero only on cycles not filled by triangles)
flow	NDArray (N, N)	The input alternating coupling current
potential	NDArray (N,)	Minimum-norm node potential s with gradient = grad(s)
betti_one	int	First Betti number β₁ — dimension of the harmonic subspace

Each flow matrix is antisymmetric (M[i, j] is the flow on the oriented edge i → j, M[j, i] = −M[i, j]).

Interpretation

• Gradient-dominated: the current is a node-potential difference and the system relaxes towards a fixed phase configuration.
• Curl-dominated: circulation around filled triangles — local cyclic frustration with no global potential.
• Harmonic component: flows around topological cycles that no triangle bounds; its dimension equals the first Betti number β₁. On a triangle-free graph carrying a cycle (for example, a 4-cycle), a circulating current is purely harmonic — the topological content that a plain symmetric/antisymmetric matrix split cannot represent. In the SCPN identity-coherence model this is the identity invariant that persists across regime changes.

hodge_decomposition(knm, phases, triangles=None) computes all three components; pass an explicit triangles list of node triples to override the default 3-clique fill.

Because the decomposition relies on two least-squares pseudoinverse solves, exact cross-language parity is not attainable; the dispatcher validates each accelerated backend against the NumPy reference within rtol = 1e-10 / atol = 1e-12 (matching the spectral solver) and falls back to NumPy only after the backend has returned a valid Hodge payload.

Direct accelerator boundary contract: the public Python dispatcher, public Rust wrapper, and the Go, Julia, and Mojo Hodge adapters reject numeric-string aliases before Python, NumPy, shared-library, Julia, or subprocess coercion. The public surface applies the boundary to knm, phases, and explicit triangle nodes; the direct adapters apply it to counts, flattened coupling, phase, edge, triangle, backend-output, and Julia raw-return payloads. The shared typed float64 path also rejects boolean aliases, complex or non-finite payloads, malformed flattened n*n coupling buffers, phase vectors whose length does not match n, and invalid oscillator counts before optional runtime loading. After backend execution, the same output validator checks that gradient, curl, and harmonic are finite real non-boolean (N, N) or flattened N*N antisymmetric matrices before publication or parity fallback. Malformed backend outputs raise immediately; fallback is reserved for validated numerical parity mismatches. Empty Hodge systems return empty components without requiring optional runtimes, matching the public Python special case.

hodge

Combinatorial Hodge (Helmholtz–Hodge) decomposition of the Kuramoto current.

Exposes a 5-backend fallback chain.

Model

The oscillator network is treated as a simplicial complex (V, E, T): vertices V are the oscillators, edges E are the unordered pairs {i, j} (i < j) carrying non-zero symmetric coupling, and triangles T are the 2-simplices (3-cliques of the coupling graph, or an explicit user-supplied set).

The decomposed object is the alternating edge flow — the Kuramoto coupling current on the reference orientation i → j (i < j):

f_{ij} = K^{sym}_{ij} · sin(θ_j − θ_i),   K^{sym} = ½(K + Kᵀ)

which satisfies f_{ji} = −f_{ij} exactly, the defining property of a 1-cochain. Using the symmetric part of K keeps the current alternating even when K encodes directed coupling; for the standard symmetric Kuramoto model K^{sym} = K.

Boundary operators

• B1 (|V| × |E|) is the node–edge incidence ∂₁: for edge e = (i, j) with i < j, B1[i, e] = −1 and B1[j, e] = +1. The discrete gradient is grad(s) = B1ᵀ s with (B1ᵀ s)_{ij} = s_j − s_i; the divergence is its adjoint B1.
• B2 (|E| × |T|) is the edge–triangle incidence ∂₂: for triangle t = {i, j, k} with i < j < k the simplicial boundary ∂[i, j, k] = [j, k] − [i, k] + [i, j] gives B2[(i,j), t] = +1, B2[(j,k), t] = +1, B2[(i,k), t] = −1. The discrete curl is curl(f) = B2ᵀ f.

Because ∂₁ ∂₂ = B1 B2 = 0 (boundary of a boundary is empty), the gradient image im(B1ᵀ) and the curl image im(B2) are L²-orthogonal.

Decomposition

With graph Laplacian L0 = B1 B1ᵀ and triangle Laplacian L2 = B2ᵀ B2:

f_grad = B1ᵀ · L0⁺ · (B1 f)     (curl-free, conservative)
f_curl = B2  · L2⁺ · (B2ᵀ f)    (divergence-free, rotational)
f_harm = f − f_grad − f_curl    (harmonic: ker of the Hodge
                                 1-Laplacian L1 = B1ᵀB1 + B2 B2ᵀ)

The three components are mutually L²-orthogonal (Jiang, Lim, Yao & Ye 2011, Theorem 2.4). The harmonic part is both divergence-free and curl-free; its dimension equals the first Betti number

β₁ = |E| − rank(B1) − rank(B2)

i.e. the number of independent cycles not bounded by triangles. On a triangle-free graph with a cycle (e.g. a square) a circulating current is purely harmonic — the topological content that a plain symmetric/antisymmetric matrix split cannot represent.

Output

:class:HodgeResult returns the three components and the input current as antisymmetric (N, N) flow matrices (M[i, j] is the flow on edge i → j), the minimum-norm node potential s such that f_grad = grad(s), and the integer betti_one (β₁).

Numerics

The decomposition needs two least-squares solves (L0⁺, L2⁺), so exact cross-language parity is not attainable; the dispatcher validates each accelerated backend against the NumPy reference within rtol = 1e-10 / atol = 1e-12 (matching the spectral solver) and falls back to NumPy on any valid numerical mismatch. Malformed backend payloads fail closed before parity fallback.

Reference

Jiang, Lim, Yao & Ye 2011, Statistical ranking and combinatorial Hodge theory, Math. Program. 127 (1):203–244.

Classes

HodgeResult dataclass

HodgeResult(
    gradient: FloatArray,
    curl: FloatArray,
    harmonic: FloatArray,
    flow: FloatArray,
    potential: FloatArray,
    betti_one: int,
)

Decompose the Kuramoto coupling current into three L²-orthogonal flows.

Each flow matrix is antisymmetric: M[i, j] is the flow on the oriented edge i → j and M[j, i] = −M[i, j].

Attributes

gradient: Conservative (curl-free) component ``grad(s)``.
curl: Rotational (divergence-free) component bounded by triangles.
harmonic: Topological residual in ``ker(L1)``; non-zero exactly
    when the graph carries cycles not filled by triangles.
flow: The input alternating coupling current
    ``K^{sym}_{ij} · sin(θ_j − θ_i)``.
potential: Minimum-norm node potential ``s`` with
    ``gradient = grad(s)``.
betti_one: First Betti number ``β₁`` — the dimension of the
    harmonic subspace.

Functions:

hodge_decomposition

hodge_decomposition(
    knm: FloatArray,
    phases: FloatArray,
    triangles: Sequence[Sequence[int]] | None = None,
) -> HodgeResult

Decompose the Kuramoto coupling current into orthogonal edge flows.

Parameters

knm : FloatArray Square (N, N) coupling matrix; the symmetric part defines the edge support and the current magnitude. phases : FloatArray (N,) oscillator phases. triangles : Sequence[Sequence[int]] | None Optional explicit 2-simplices as node triples; each must reference existing edges. When omitted, all 3-cliques of the coupling graph are used.

Returns

HodgeResult :class:HodgeResult with the three flow components as antisymmetric (N, N) matrices, the input current, the node potential, and the first Betti number.

Source code in src/scpn_phase_orchestrator/coupling/hodge.py

def hodge_decomposition(
    knm: FloatArray,
    phases: FloatArray,
    triangles: Sequence[Sequence[int]] | None = None,
) -> HodgeResult:
    """Decompose the Kuramoto coupling current into orthogonal edge flows.

    Parameters
    ----------
    knm : FloatArray
        Square ``(N, N)`` coupling matrix; the symmetric part defines the edge support
        and the current magnitude.
    phases : FloatArray
        ``(N,)`` oscillator phases.
    triangles : Sequence[Sequence[int]] | None
        Optional explicit 2-simplices as node triples; each must reference existing
        edges. When omitted, all 3-cliques of the coupling graph are used.

    Returns
    -------
    HodgeResult
        :class:`HodgeResult` with the three flow components as antisymmetric ``(N, N)``
        matrices, the input current, the node potential, and the first Betti number.
    """
    phases = _validate_phase_vector(phases, name="phases")
    n = int(phases.size)
    if n == 0:
        empty = np.zeros((0, 0), dtype=np.float64)
        return HodgeResult(
            gradient=empty,
            curl=empty.copy(),
            harmonic=empty.copy(),
            flow=empty.copy(),
            potential=np.array([], dtype=np.float64),
            betti_one=0,
        )

    k = _validate_coupling_matrix(knm, expected_n=n)
    k_sym = 0.5 * (k + k.T)
    edges, tri = _simplicial_complex(k_sym, n, triangles)

    gradient, curl, harmonic, flow_matrix, potential, betti = _decompose(
        k, phases, edges, tri
    )
    reference: HodgeTuple = (gradient, curl, harmonic)

    n_edges = int(edges.shape[0])
    n_tris = int(tri.shape[0])
    edges_flat = np.ascontiguousarray(edges.ravel(), dtype=np.int64)
    tris_flat = np.ascontiguousarray(tri.ravel(), dtype=np.int64)
    k_flat = np.ascontiguousarray(k.ravel(), dtype=np.float64)

    backend_fn = _dispatch()
    if backend_fn is not None:
        backend_output = _normalise_backend_output(
            backend_fn(
                k_flat,
                phases,
                n,
                edges_flat,
                n_edges,
                tris_flat,
                n_tris,
            ),
            expected_n=n,
        )
        if _backend_matches_reference(backend_output, reference):
            g, c, h = backend_output
            return HodgeResult(
                gradient=g,
                curl=c,
                harmonic=h,
                flow=flow_matrix,
                potential=potential,
                betti_one=betti,
            )

    return HodgeResult(
        gradient=gradient,
        curl=curl,
        harmonic=harmonic,
        flow=flow_matrix,
        potential=potential,
        betti_one=betti,
    )

---

Spectral Analysis

Algebraic graph-theoretic properties of the coupling network.

Functions

Function	Returns	Description
graph_laplacian(knm)	NDArray	L = D - W (combinatorial Laplacian)
fiedler_value(knm)	float	λ₂(L) — algebraic connectivity
fiedler_vector(knm)	NDArray	Eigenvector of λ₂
critical_coupling(omegas, knm)	float	K_c = max|Δω| / λ₂
fiedler_partition(knm)	(list, list)	Network bisection via Fiedler sign
spectral_gap(knm)	float	λ₃ - λ₂ (cluster clarity)
sync_convergence_rate(knm, omegas, γ_max)	float	μ = K·λ₂·cos(γ)/N

Critical coupling estimate

The Dörfler-Bullo bound gives the minimum coupling strength for synchronisation:

K_c = max_{i,j} |ω_i - ω_j| / λ₂(L)

where λ₂ is the Fiedler eigenvalue (algebraic connectivity). Networks with higher λ₂ synchronise more easily.

Direct accelerator boundary contract: Go, Julia, and Mojo spectral adapters use one shared typed float64 validation path before loading shared-library, Julia, or subprocess runtimes. The contract rejects boolean aliases, numeric-string aliases, complex or non-finite flattened coupling payloads, non-vector inputs, malformed n*n buffer lengths, and invalid oscillator counts. Empty spectral problems return empty eigenvalue and Fiedler vectors without optional runtime loading. After backend execution, the same shared output validator is replayed for the direct Go, Julia, and Mojo adapters and for the public optional primitive path: returned eigenvalues and the Fiedler vector must be finite real non-boolean, non-numeric-string vectors of length N, eigenvalues must be non-negative and sorted ascending, and the Fiedler vector must be non-zero for N > 1. Malformed backend physics payloads raise immediately; fallback remains reserved for loader or runtime unavailability. Public spectral helpers enforce the same real-valued boundary on coupling matrices, frequency vectors, gamma_max, optional primitive eigensystem outputs, and Rust fast-path scalar/vector returns. Boolean aliases are not coerced into weights or frequencies, and complex-valued aliases are rejected before NumPy can discard imaginary components. Numeric-string aliases are rejected before Python, NumPy, Rust, Julia, Go, or Mojo can widen them into ordinary floating-point weights, frequencies, scalar controls, or eigensystem payloads.

spectral

Symmetric eigendecomposition of the combinatorial graph Laplacian L = D − A.

Exposes a 5-backend fallback chain.

For asymmetric measured coupling, the undirected adjacency is the reciprocal magnitude average A = (|W| + |Wᵀ|) / 2 with zeroed diagonal. Degrees are then computed from A. This preserves the combinatorial Laplacian contract: symmetric positive-semidefinite L, zero row sums, and 1 ∈ ker L.

Primitive

spectral_eig(W_flat, n) → (eigvals, fiedler) — eigenvalues ascending + Fiedler eigenvector (column 1 of the sorted decomposition).

Backend chain

• Rust: pre-existing fiedler_value_rust, fiedler_vector_rust, spectral_gap_rust, critical_coupling_rust, sync_convergence_rate_rust FFI fast paths are wired individually (each exposes a direct entry, no round-trip through the primitive).
• Julia: LinearAlgebra.eigen(Symmetric(L)) — LAPACK dsyev underneath, same numerics as NumPy.
• Go: gonum.org/v1/gonum/mat :: EigenSym — pure-Go symmetric solver, sub-1e-12 drift vs LAPACK on well-conditioned Laplacians.
• Mojo: LAPACK dsyev_ via the std.ffi.OwnedDLHandle pattern (same as _lapack_test.mojo).
• Python: np.linalg.eigh — LAPACK-backed reference.

Derived functions (fiedler_value, fiedler_vector, spectral_gap) route through the primitive on non-Rust backends. critical_coupling and sync_convergence_rate are composites that reuse fiedler_value.

References: Dörfler & Bullo 2014, Automatica 50(6):1539-1564; Dörfler & Bullo 2013, IEEE Proc. 102(10):1539-1564.

Functions:

graph_laplacian

graph_laplacian(knm: FloatArray) -> FloatArray

Combinatorial graph Laplacian L = D − A.

A is the reciprocal undirected magnitude adjacency (|W| + |Wᵀ|) / 2 with zero diagonal, so asymmetric measured couplings produce one symmetric edge weight before node degrees are computed.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

FloatArray The combinatorial graph Laplacian L = D − A.

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def graph_laplacian(knm: FloatArray) -> FloatArray:
    """Combinatorial graph Laplacian ``L = D − A``.

    ``A`` is the reciprocal undirected magnitude adjacency
    ``(|W| + |Wᵀ|) / 2`` with zero diagonal, so asymmetric measured
    couplings produce one symmetric edge weight before node degrees
    are computed.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    FloatArray
        The combinatorial graph Laplacian ``L = D − A``.
    """
    knm = _validate_coupling_matrix(knm)
    w = np.abs(knm)
    np.fill_diagonal(w, 0.0)
    adjacency = 0.5 * (w + w.T)
    np.fill_diagonal(adjacency, 0.0)
    degrees = adjacency.sum(axis=1)
    return cast("FloatArray", np.diag(degrees) - adjacency)

spectral_eig

spectral_eig(
    knm: FloatArray,
) -> tuple[FloatArray, FloatArray]

Symmetric eigendecomposition of L = D − |W|.

Returns (eigvals ascending, fiedler vector). Thin wrapper over the dispatched backend primitive; python reference is a direct np.linalg.eigh.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

tuple[FloatArray, FloatArray] The eigenvalues and eigenvectors of the symmetric Laplacian.

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def spectral_eig(knm: FloatArray) -> tuple[FloatArray, FloatArray]:
    """Symmetric eigendecomposition of ``L = D − |W|``.

    Returns ``(eigvals ascending, fiedler vector)``. Thin wrapper
    over the dispatched backend primitive; ``python`` reference
    is a direct ``np.linalg.eigh``.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    tuple[FloatArray, FloatArray]
        The eigenvalues and eigenvectors of the symmetric Laplacian.
    """
    knm = _validate_coupling_matrix(knm)
    n = knm.shape[0]
    flat = np.ascontiguousarray(knm.ravel(), dtype=np.float64)
    return _spectral_eig_checked(flat, n)

fiedler_value

fiedler_value(knm: FloatArray) -> float

Return the algebraic connectivity λ₂(L) (Dörfler-Bullo 2014).

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

float The algebraic connectivity λ₂(L).

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def fiedler_value(knm: FloatArray) -> float:
    """Return the algebraic connectivity ``λ₂(L)`` (Dörfler-Bullo 2014).

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    float
        The algebraic connectivity ``λ₂(L)``.
    """
    knm = _validate_coupling_matrix(knm)
    n = knm.shape[0]
    if n < 2:
        return 0.0
    flat = np.ascontiguousarray(knm.ravel(), dtype=np.float64)
    if ACTIVE_BACKEND == "rust":
        return _validate_non_negative_scalar(
            _rust_bundle()["fv"](flat, n), name="Fiedler value"
        )
    eigvals, _ = _spectral_eig_checked(flat, n)
    return float(eigvals[1])

fiedler_vector

fiedler_vector(knm: FloatArray) -> FloatArray

Return the λ₂ eigenvector partitioning the graph into clusters.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

FloatArray The λ₂ eigenvector partitioning the graph.

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def fiedler_vector(knm: FloatArray) -> FloatArray:
    """Return the ``λ₂`` eigenvector partitioning the graph into clusters.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    FloatArray
        The ``λ₂`` eigenvector partitioning the graph.
    """
    knm = _validate_coupling_matrix(knm)
    n = knm.shape[0]
    flat = np.ascontiguousarray(knm.ravel(), dtype=np.float64)
    if ACTIVE_BACKEND == "rust":
        return _validate_rust_fiedler_vector(_rust_bundle()["fvec"](flat, n), n=n)
    _, fiedler = _spectral_eig_checked(flat, n)
    return fiedler

critical_coupling

critical_coupling(
    omegas: FloatArray, knm: FloatArray
) -> float

Dörfler-Bullo critical coupling K_c = Δω / λ₂.

Returns +inf if the graph is disconnected (λ₂ ≈ 0).

Parameters

omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

float The Dörfler-Bullo critical coupling K_c.

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def critical_coupling(omegas: FloatArray, knm: FloatArray) -> float:
    """Dörfler-Bullo critical coupling ``K_c = Δω / λ₂``.

    Returns ``+inf`` if the graph is disconnected
    (``λ₂ ≈ 0``).

    Parameters
    ----------
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    float
        The Dörfler-Bullo critical coupling ``K_c``.
    """
    knm = _validate_coupling_matrix(knm)
    n = knm.shape[0]
    omegas = _validate_omegas(omegas, expected_n=n)
    if ACTIVE_BACKEND == "rust":
        flat = np.ascontiguousarray(knm.ravel(), dtype=np.float64)
        o = np.ascontiguousarray(omegas, dtype=np.float64)
        return _validate_non_negative_scalar(
            _rust_bundle()["kc"](o, flat, n),
            name="critical coupling",
            allow_infinite=True,
        )
    lambda2 = fiedler_value(knm)
    if lambda2 < 1e-12:
        return float("inf")
    omega_spread = float(np.max(omegas) - np.min(omegas))
    return omega_spread / lambda2

fiedler_partition

fiedler_partition(
    knm: FloatArray,
) -> tuple[list[int], list[int]]

Bisect the network using sign(v₂).

Returns (group_positive, group_negative) — indices of oscillators in each partition.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

tuple[list[int], list[int]] The two index lists of the sign(v₂) bisection.

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def fiedler_partition(knm: FloatArray) -> tuple[list[int], list[int]]:
    """Bisect the network using ``sign(v₂)``.

    Returns ``(group_positive, group_negative)`` — indices
    of oscillators in each partition.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    tuple[list[int], list[int]]
        The two index lists of the ``sign(v₂)`` bisection.
    """
    v2 = fiedler_vector(knm)
    pos = [i for i, val in enumerate(v2) if val >= 0]
    neg = [i for i, val in enumerate(v2) if val < 0]
    return pos, neg

spectral_gap

spectral_gap(knm: FloatArray) -> float

Return the gap between λ₂ and λ₃ (two-cluster cleanliness).

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N).

Returns

float The gap between λ₂ and λ₃.

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def spectral_gap(knm: FloatArray) -> float:
    """Return the gap between ``λ₂`` and ``λ₃`` (two-cluster cleanliness).

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    float
        The gap between ``λ₂`` and ``λ₃``.
    """
    knm = _validate_coupling_matrix(knm)
    n = knm.shape[0]
    if n < 3:
        return 0.0
    off_diag = np.abs(knm[~np.eye(n, dtype=bool)])
    if off_diag.size and np.allclose(off_diag, off_diag[0], rtol=1e-12, atol=1e-12):
        return 0.0
    flat = np.ascontiguousarray(knm.ravel(), dtype=np.float64)
    if ACTIVE_BACKEND == "rust":
        return _validate_non_negative_scalar(
            _rust_bundle()["sg"](flat, n), name="spectral gap"
        )
    eigvals, _ = _spectral_eig_checked(flat, n)
    return float(eigvals[2] - eigvals[1])

sync_convergence_rate

sync_convergence_rate(
    knm: FloatArray,
    omegas: FloatArray,
    gamma_max: float = 0.0,
) -> float

Estimate the convergence rate from λ₂ (Dörfler-Bullo 2014 §III.B).

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N). omegas : FloatArray Natural frequencies in rad/s, shape (N,). gamma_max : float Maximum phase-lag γ across edges.

Returns

float The estimated synchronisation convergence rate.

Source code in src/scpn_phase_orchestrator/coupling/spectral.py

def sync_convergence_rate(
    knm: FloatArray,
    omegas: FloatArray,
    gamma_max: float = 0.0,
) -> float:
    """Estimate the convergence rate from ``λ₂`` (Dörfler-Bullo 2014 §III.B).

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    gamma_max : float
        Maximum phase-lag ``γ`` across edges.

    Returns
    -------
    float
        The estimated synchronisation convergence rate.
    """
    knm = _validate_coupling_matrix(knm)
    n = knm.shape[0]
    if n == 0:
        return 0.0
    omegas = _validate_omegas(omegas, expected_n=n)
    gamma_max = _validate_gamma_max(gamma_max)
    if ACTIVE_BACKEND == "rust":
        flat = np.ascontiguousarray(knm.ravel(), dtype=np.float64)
        o = np.ascontiguousarray(omegas, dtype=np.float64)
        return _validate_non_negative_scalar(
            _rust_bundle()["scr"](flat, o, n, gamma_max),
            name="sync convergence rate",
        )
    lambda2 = fiedler_value(knm)
    pos_vals = knm[knm > 0]
    k_eff = float(np.mean(pos_vals)) if pos_vals.size > 0 else 0.0
    return float(k_eff * lambda2 * np.cos(gamma_max) / n)

---

Three-Factor Hebbian Plasticity

Coupling adaptation rule inspired by biological synaptic plasticity:

ΔK_ij = lr × eligibility_ij × modulator × phase_gate

Functions

• compute_eligibility(phases) → NDArray(n,n): pairwise Hebbian trace cos(θ_j - θ_i) with zero diagonal. In-phase pairs → +1 (strengthen), anti-phase → -1 (weaken).

• three_factor_update(knm, eligibility, modulator, phase_gate, lr=0.01) → NDArray: applies the three-factor rule. Only modifies K when all three factors are active. The boundary enforces the same physical K_nm contract consumed by the UPDE engines: knm must be finite, real, non-negative, square, and zero-diagonal; eligibility must be finite, real, square, zero-diagonal, and bounded in [-1, 1]. Negative modulation can depress coupling but is clamped at zero, and the result always keeps a zero self-coupling diagonal.

Three factors

1. Eligibility (local): cos(Δθ) — pairwise Hebbian trace
2. Modulator (global): scalar from L16 director layer (dopamine analog)
3. Phase gate (global): Boolean from topological-integration gate

Reference: Friston 2005 on free energy and synaptic plasticity.

plasticity

Validated three-factor plasticity updates for coupling matrices.

The module computes pairwise phase eligibility traces and applies a modulator-gated Hebbian update to K_nm. Public functions reject boolean, non-numeric, non-finite, non-vector, non-square, and shape-mismatched inputs so plasticity cannot corrupt coupling state silently. The update preserves the Kuramoto coupling contract by requiring non-negative zero-diagonal K_nm, bounded zero-diagonal eligibility traces, and finite real scalar controls.

Functions:

compute_eligibility

compute_eligibility(phases: FloatArray) -> FloatArray

Pairwise Hebbian eligibility trace: cos(theta_j - theta_i).

Returns shape (n, n) with zero diagonal.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,).

Returns

FloatArray The pairwise Hebbian eligibility trace cos(θ_j − θ_i).

Source code in src/scpn_phase_orchestrator/coupling/plasticity.py

def compute_eligibility(phases: FloatArray) -> FloatArray:
    """Pairwise Hebbian eligibility trace: cos(theta_j - theta_i).

    Returns shape (n, n) with zero diagonal.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.

    Returns
    -------
    FloatArray
        The pairwise Hebbian eligibility trace ``cos(θ_j − θ_i)``.
    """
    phases = _validate_phase_vector(phases, name="phases")
    diffs = phases[np.newaxis, :] - phases[:, np.newaxis]
    elig = np.cos(diffs)
    np.fill_diagonal(elig, 0.0)
    result: FloatArray = elig
    return result

three_factor_update

three_factor_update(
    knm: FloatArray,
    eligibility: FloatArray,
    modulator: float,
    phase_gate: bool,
    lr: float = 0.01,
) -> FloatArray

Three-factor plasticity rule: K_ij += lr * eligibility_ij * M * gate.

Factors

1. eligibility — pairwise phase correlation (Hebbian trace)
2. modulator — scalar reward/error signal from L16 director
3. phase_gate — boolean from the topological-integration gate

Friston 2005, Philos. Trans. R. Soc. B 360:815-836 (free energy & synaptic plasticity).

Parameters

knm : FloatArray current coupling matrix, shape (n, n). eligibility : FloatArray Hebbian trace, shape (n, n). modulator : float scalar neuromodulatory signal. phase_gate : bool if False, no update occurs (integration gate below threshold). lr : float learning rate.

Returns

FloatArray Updated coupling matrix (new array, does not mutate input).

Raises

TypeError If an argument has the wrong type. ValueError If the eligibility or coupling shapes mismatch.

Source code in src/scpn_phase_orchestrator/coupling/plasticity.py

def three_factor_update(
    knm: FloatArray,
    eligibility: FloatArray,
    modulator: float,
    phase_gate: bool,
    lr: float = 0.01,
) -> FloatArray:
    """Three-factor plasticity rule: K_ij += lr * eligibility_ij * M * gate.

    Factors:
        1. eligibility — pairwise phase correlation (Hebbian trace)
        2. modulator — scalar reward/error signal from L16 director
        3. phase_gate — boolean from the topological-integration gate

    Friston 2005, Philos. Trans. R. Soc. B
    360:815-836 (free energy & synaptic plasticity).

    Parameters
    ----------
    knm : FloatArray
        current coupling matrix, shape (n, n).
    eligibility : FloatArray
        Hebbian trace, shape (n, n).
    modulator : float
        scalar neuromodulatory signal.
    phase_gate : bool
        if False, no update occurs (integration gate below threshold).
    lr : float
        learning rate.

    Returns
    -------
    FloatArray
        Updated coupling matrix (new array, does not mutate input).

    Raises
    ------
    TypeError
        If an argument has the wrong type.
    ValueError
        If the eligibility or coupling shapes mismatch.
    """
    knm = _validate_coupling_matrix(knm)
    eligibility = _validate_eligibility_matrix(eligibility)
    if eligibility.shape != knm.shape:
        raise ValueError(
            "eligibility shape "
            f"{eligibility.shape} does not match knm shape {knm.shape}"
        )
    modulator = _validate_finite_real(modulator, name="modulator")
    if not isinstance(phase_gate, bool):
        raise TypeError("phase_gate must be a bool")
    lr = _validate_learning_rate(lr)
    if not phase_gate:
        return knm.copy()
    delta = lr * eligibility * modulator
    updated = np.maximum(knm + delta, 0.0)
    np.fill_diagonal(updated, 0.0)
    result: FloatArray = updated
    return result

---

Transfer Entropy Adaptive Coupling

Directed causal adaptation that breaks symmetry:

K_ij(t+1) = (1 - decay) × K_ij(t) + lr × TE(i → j)

te_adapt_coupling(knm, phase_history, lr=0.01, decay=0.0, n_bins=8):

• Computes transfer entropy TE(i→j) for all pairs from phase history
• Updates coupling: pairs with causal influence get stronger
• Applies decay to forget old coupling structure
• Clamps K ≥ 0 and zeros diagonal
• Rejects boolean aliases in both knm and phase_history before numeric coercion

Unlike Hebbian plasticity (symmetric), TE captures directed information flow — oscillator i can influence j without j influencing i.

Reference: Lizier 2012, "Local Information Transfer as Spatiotemporal Filter." Detailed documentation: TE Adaptive — detailed reference

te_adaptive

Transfer-entropy-guided coupling adaptation for offline matrix updates.

te_adapt_coupling derives a directed transfer-entropy matrix from phase history and combines it with the current coupling matrix under learning-rate and decay parameters. The Python fallback clamps the returned coupling to non-negative values and clears self-coupling; the optional Rust path preserves the same dense N x N output contract. The helper returns a new matrix and does not mutate live solver state or apply actuation.

Functions:

te_adapt_coupling

te_adapt_coupling(
    knm: FloatArray,
    phase_history: FloatArray,
    lr: float = 0.01,
    decay: float = 0.0,
    n_bins: int = 8,
) -> FloatArray

Adapt coupling matrix using transfer entropy as learning signal.

K_ij(t+1) = (1-decay) * K_ij(t) + lr * TE(i→j)

Strengthens coupling along causal information flow channels. Weakens where there is no causal influence.

Lizier 2012, "Local Information Transfer as a Spatiotemporal Filter for Complex Systems," Physical Review E 77(2):026110.

Parameters

knm : FloatArray current (n, n) coupling matrix. phase_history : FloatArray (n, T) recent phase trajectories. lr : float learning rate for TE-based update. decay : float coupling decay rate per update (0 = no decay). n_bins : int histogram bins for TE estimation.

Returns

FloatArray FloatArray The coupling matrix adapted by the transfer-entropy learning signal.

Raises

RuntimeError If the transfer-entropy backend fails.

Source code in src/scpn_phase_orchestrator/coupling/te_adaptive.py

def te_adapt_coupling(
    knm: FloatArray,
    phase_history: FloatArray,
    lr: float = 0.01,
    decay: float = 0.0,
    n_bins: int = 8,
) -> FloatArray:
    """Adapt coupling matrix using transfer entropy as learning signal.

    K_ij(t+1) = (1-decay) * K_ij(t) + lr * TE(i→j)

    Strengthens coupling along causal information flow channels.
    Weakens where there is no causal influence.

    Lizier 2012, "Local Information Transfer as a Spatiotemporal Filter
    for Complex Systems," Physical Review E 77(2):026110.

    Parameters
    ----------
    knm : FloatArray
        current (n, n) coupling matrix.
    phase_history : FloatArray
        (n, T) recent phase trajectories.
    lr : float
        learning rate for TE-based update.
    decay : float
        coupling decay rate per update (0 = no decay).
    n_bins : int
        histogram bins for TE estimation.

    Returns
    -------
    FloatArray
        FloatArray The coupling matrix adapted by the transfer-entropy learning signal.

    Raises
    ------
    RuntimeError
        If the transfer-entropy backend fails.
    """
    knm = _validate_knm(knm)
    n = knm.shape[0]
    phase_history = _validate_phase_history(phase_history, n=n)
    lr = _validate_non_negative_real(lr, name="lr")
    decay = _validate_decay(decay)
    n_bins = _validate_n_bins(n_bins)
    te = _validate_transfer_entropy_scores(
        transfer_entropy_matrix(phase_history, n_bins=n_bins),
        n=n,
    )
    if _HAS_RUST:
        k_flat = np.ascontiguousarray(knm.ravel(), dtype=np.float64)
        t_flat = np.ascontiguousarray(te.ravel(), dtype=np.float64)
        result_flat = np.asarray(
            _rust_te_adapt(k_flat, t_flat, n, lr, decay),
            dtype=np.float64,
        )
        if result_flat.size != n * n:
            raise RuntimeError("TE adaptive backend returned wrong shape")
        return _validate_adapted_coupling(result_flat.reshape(n, n), n=n)
    knm_new = (1.0 - decay) * knm + lr * te
    np.fill_diagonal(knm_new, 0.0)
    result: FloatArray = np.maximum(knm_new, 0.0)
    return result

---

E/I Balance

Computes and adjusts excitatory/inhibitory coupling balance. The aggregate ratio summarises overall balance, while the four directed interaction-type means resolve it into the source→target block strengths that Kuroki & Mizuseki 2025 (Neural Computation 37 (7):1353–1372) identify as the control parameters of the EI-Kuramoto synchronised / bistable / desynchronised regimes.

EIBalance (dataclass)

Field	Type	Description
ratio	float	E/I balance ratio (excitatory_strength / inhibitory_strength)
excitatory_strength	float	Mean coupling from excitatory sources over all targets
inhibitory_strength	float	Mean coupling from inhibitory sources over all targets
is_balanced	bool	True if 0.8 ≤ ratio ≤ 1.2
e_to_e	float	Mean E→E interaction-type coupling
e_to_i	float	Mean E→I interaction-type coupling
i_to_e	float	Mean I→E interaction-type coupling
i_to_i	float	Mean I→I interaction-type coupling

Each aggregate strength is the count-weighted blend of its two outgoing interaction-type blocks (e.g. excitatory_strength blends e_to_e and e_to_i over the target-group sizes).

Functions

• compute_ei_balance(knm, excitatory_indices, inhibitory_indices) → EIBalance
• adjust_ei_ratio(knm, excitatory_indices, inhibitory_indices, target_ratio=1.0) → NDArray — scales inhibitory coupling to achieve target ratio

Both helpers reject boolean aliases in knm before computing row means or scaling inhibitory rows.

ei_balance

Excitatory/inhibitory balance summaries and adjustment helpers.

The module measures mean outgoing coupling from caller-specified excitatory and inhibitory index sets, then optionally rescales inhibitory rows toward a target ratio. Rust acceleration is used when available; the NumPy fallback preserves the same shape and summary contract for examples and deterministic tests.

Classes

EIBalance dataclass

EIBalance(
    ratio: float,
    excitatory_strength: float,
    inhibitory_strength: float,
    is_balanced: bool,
    e_to_e: float,
    e_to_i: float,
    i_to_e: float,
    i_to_i: float,
)

Summary of excitatory and inhibitory coupling balance.

excitatory_strength / inhibitory_strength aggregate the mean coupling from each source group over all targets, and ratio is their quotient. The four *_to_* block means resolve this into the directed interaction-type strengths (source group → target group) that Kuroki & Mizuseki 2025 identify as the control parameters of the synchronised, bistable, and desynchronised regimes of the EI-Kuramoto model.

Functions:

compute_ei_balance

compute_ei_balance(
    knm: FloatArray,
    excitatory_indices: list[int],
    inhibitory_indices: list[int],
) -> EIBalance

Compute E/I balance from coupling matrix and layer typing.

Kuroki & Mizuseki 2025, Neural Computation — E/I balance is the critical parameter for synchronization, not K or D.

ratio > 1: excitation-dominated (hypersynchrony risk) ratio < 1: inhibition-dominated (desynchronization risk) ratio ≈ 1: balanced (optimal for metastability)

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N). excitatory_indices : list[int] Indices of the excitatory oscillators. inhibitory_indices : list[int] Indices of the inhibitory oscillators.

Returns

EIBalance The E/I balance summary derived from the coupling typing.

Source code in src/scpn_phase_orchestrator/coupling/ei_balance.py

def compute_ei_balance(
    knm: FloatArray,
    excitatory_indices: list[int],
    inhibitory_indices: list[int],
) -> EIBalance:
    """Compute E/I balance from coupling matrix and layer typing.

    Kuroki & Mizuseki 2025, Neural Computation — E/I balance is the
    critical parameter for synchronization, not K or D.

    ratio > 1: excitation-dominated (hypersynchrony risk)
    ratio < 1: inhibition-dominated (desynchronization risk)
    ratio ≈ 1: balanced (optimal for metastability)

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    excitatory_indices : list[int]
        Indices of the excitatory oscillators.
    inhibitory_indices : list[int]
        Indices of the inhibitory oscillators.

    Returns
    -------
    EIBalance
        The E/I balance summary derived from the coupling typing.
    """
    knm = _validate_knm(knm)
    n = knm.shape[0]
    excitatory_indices = _validate_indices(excitatory_indices, n, "excitatory")
    inhibitory_indices = _validate_indices(inhibitory_indices, n, "inhibitory")

    if _HAS_RUST:
        k_flat = np.ascontiguousarray(knm.ravel())
        e_arr = np.array(excitatory_indices, dtype=np.int64)
        i_arr = np.array(inhibitory_indices, dtype=np.int64)
        ratio, e_str, i_str, balanced, e_to_e, e_to_i, i_to_e, i_to_i = _rust_ei(
            k_flat, n, e_arr, i_arr
        )
        return EIBalance(
            ratio=float(ratio),
            excitatory_strength=float(e_str),
            inhibitory_strength=float(i_str),
            is_balanced=bool(balanced),
            e_to_e=float(e_to_e),
            e_to_i=float(e_to_i),
            i_to_e=float(i_to_e),
            i_to_i=float(i_to_i),
        )

    e_mask = np.zeros(n, dtype=bool)
    i_mask = np.zeros(n, dtype=bool)
    for idx in excitatory_indices:
        e_mask[idx] = True
    for idx in inhibitory_indices:
        i_mask[idx] = True

    # Excitatory strength: mean coupling FROM excitatory oscillators
    e_strength = float(np.mean(knm[e_mask, :])) if np.any(e_mask) else 0.0
    # Inhibitory strength: mean coupling FROM inhibitory oscillators
    i_strength = float(np.mean(knm[i_mask, :])) if np.any(i_mask) else 0.0

    if i_strength < 1e-15:
        ratio = float("inf") if e_strength > 0 else 1.0
    else:
        ratio = e_strength / i_strength

    return EIBalance(
        ratio=ratio,
        excitatory_strength=e_strength,
        inhibitory_strength=i_strength,
        is_balanced=0.8 <= ratio <= 1.2,
        e_to_e=_block_mean(knm, e_mask, e_mask),
        e_to_i=_block_mean(knm, e_mask, i_mask),
        i_to_e=_block_mean(knm, i_mask, e_mask),
        i_to_i=_block_mean(knm, i_mask, i_mask),
    )

adjust_ei_ratio

adjust_ei_ratio(
    knm: FloatArray,
    excitatory_indices: list[int],
    inhibitory_indices: list[int],
    target_ratio: float = 1.0,
) -> FloatArray

Scale inhibitory coupling to achieve target E/I ratio.

Returns modified knm with inhibitory rows scaled so that E_strength / I_strength ≈ target_ratio.

Parameters

knm : FloatArray Coupling matrix K_nm, shape (N, N). excitatory_indices : list[int] Indices of the excitatory oscillators. inhibitory_indices : list[int] Indices of the inhibitory oscillators. target_ratio : float Target excitatory/inhibitory coupling ratio.

Returns

FloatArray The coupling matrix with inhibitory weights scaled to the target ratio.

Source code in src/scpn_phase_orchestrator/coupling/ei_balance.py

def adjust_ei_ratio(
    knm: FloatArray,
    excitatory_indices: list[int],
    inhibitory_indices: list[int],
    target_ratio: float = 1.0,
) -> FloatArray:
    """Scale inhibitory coupling to achieve target E/I ratio.

    Returns modified knm with inhibitory rows scaled so that
    E_strength / I_strength ≈ target_ratio.

    Parameters
    ----------
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    excitatory_indices : list[int]
        Indices of the excitatory oscillators.
    inhibitory_indices : list[int]
        Indices of the inhibitory oscillators.
    target_ratio : float
        Target excitatory/inhibitory coupling ratio.

    Returns
    -------
    FloatArray
        The coupling matrix with inhibitory weights scaled to the target ratio.
    """
    knm = _validate_knm(knm)
    target_ratio = _validate_target_ratio(target_ratio)
    n = knm.shape[0]
    excitatory_indices = _validate_indices(excitatory_indices, n, "excitatory")
    inhibitory_indices = _validate_indices(inhibitory_indices, n, "inhibitory")

    if _HAS_RUST:
        k_flat = np.ascontiguousarray(knm.ravel())
        e_arr = np.array(excitatory_indices, dtype=np.int64)
        i_arr = np.array(inhibitory_indices, dtype=np.int64)
        result_flat: FloatArray = np.asarray(
            _rust_adjust(k_flat, n, e_arr, i_arr, target_ratio),
        )
        return result_flat.reshape(n, n)

    balance = compute_ei_balance(knm, excitatory_indices, inhibitory_indices)
    if balance.inhibitory_strength < 1e-15 or balance.excitatory_strength < 1e-15:
        return knm.copy()

    current_ratio = balance.ratio
    if abs(current_ratio - target_ratio) < 1e-10:
        return knm.copy()

    # Scale inhibitory rows: I_new = I_old * (current_ratio / target_ratio)
    scale = current_ratio / target_ratio
    result: FloatArray = knm.copy()
    for idx in inhibitory_indices:
        result[idx, :] *= scale
    return result

---

Universal Bayesian Prior

Gaussian prior over coupling parameters, calibrated from the SCPN experimental programme.

CouplingPrior (dataclass)

Field	Type	Default
K_base	float	0.47
decay_alpha	float	0.25
K_c_estimate	float	0.0

UniversalPrior

• default() → CouplingPrior — MAP estimate (K_base=0.47, α=0.25)
• sample(rng=None, seed=None) → CouplingPrior — random draw from prior; seed must be an integer in the unsigned 64-bit range when provided
• estimate_Kc(omegas, n_layers) → CouplingPrior — combines prior with Dörfler-Bullo K_c for a finite one-dimensional frequency vector
• log_probability(K_base, decay_alpha) → float — unnormalised log-probability under Gaussian prior

estimate_Kc rejects boolean aliases in omegas, including NumPy boolean scalars carried inside object arrays, before constructing the prior graph.

Detailed documentation: Universal Prior — detailed reference

prior

Empirical domain-agnostic prior for coupling hyperparameters.

UniversalPrior provides default/sample/log-probability helpers for K_base and decay_alpha, plus a Dörfler-Bullo-style critical-coupling estimate over the spectral module. When the optional Rust kernel is importable the log-probability path dispatches there; otherwise the NumPy scalar fallback preserves the same Gaussian-prior contract.

Classes

CouplingPrior dataclass

CouplingPrior(
    K_base: float, decay_alpha: float, K_c_estimate: float
)

Coupling configuration: base strength, decay, and K_c estimate.

UniversalPrior

UniversalPrior(
    K_base_mean: float = _K_BASE_MEAN,
    K_base_std: float = _K_BASE_STD,
    decay_alpha_mean: float = _DECAY_ALPHA_MEAN,
    decay_alpha_std: float = _DECAY_ALPHA_STD,
)

Domain-agnostic coupling prior from 25-domainpack empirical distribution.

K_base ~ N(0.47, 0.09), decay_alpha ~ N(0.25, 0.07). Any new domain starts from this prior. Combined with Dörfler-Bullo K_c, collapses auto-tune from 5D optimization to 2D.

Source: R4-A3 cross-domain transfer analysis (Stankovski 2017, Rev. Mod. Phys.).

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def __init__(
    self,
    K_base_mean: float = _K_BASE_MEAN,
    K_base_std: float = _K_BASE_STD,
    decay_alpha_mean: float = _DECAY_ALPHA_MEAN,
    decay_alpha_std: float = _DECAY_ALPHA_STD,
):
    self._K_base_mean = _validate_finite_real(K_base_mean, name="K_base_mean")
    self._K_base_std = _validate_positive_real(K_base_std, name="K_base_std")
    self._decay_alpha_mean = _validate_finite_real(
        decay_alpha_mean, name="decay_alpha_mean"
    )
    self._decay_alpha_std = _validate_positive_real(
        decay_alpha_std, name="decay_alpha_std"
    )

Methods:

sample

sample(
    rng: Generator | None = None, seed: int | None = None
) -> CouplingPrior

Draw a random coupling configuration from the prior.

Pass rng for an explicit generator, or seed to create a seeded one. If neither is given, a fresh unseeded generator is used (NOT reproducible across sessions).

Parameters

rng : np.random.Generator | None NumPy random generator, or None to seed from seed. seed : int | None Seed for the deterministic RNG.

Returns

CouplingPrior A coupling configuration sampled from the prior.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def sample(
    self,
    rng: np.random.Generator | None = None,
    seed: int | None = None,
) -> CouplingPrior:
    """Draw a random coupling configuration from the prior.

    Pass ``rng`` for an explicit generator, or ``seed`` to create a
    seeded one. If neither is given, a fresh unseeded generator is used
    (NOT reproducible across sessions).

    Parameters
    ----------
    rng : np.random.Generator | None
        NumPy random generator, or ``None`` to seed from ``seed``.
    seed : int | None
        Seed for the deterministic RNG.

    Returns
    -------
    CouplingPrior
        A coupling configuration sampled from the prior.
    """
    if rng is None:
        seed = _validate_seed(seed)
        rng = np.random.default_rng(seed)
    K = max(0.01, rng.normal(self._K_base_mean, self._K_base_std))
    alpha = max(0.01, rng.normal(self._decay_alpha_mean, self._decay_alpha_std))
    return CouplingPrior(K_base=K, decay_alpha=alpha, K_c_estimate=0.0)

default

default() -> CouplingPrior

Return the MAP (maximum a posteriori) estimate = the means.

Returns

CouplingPrior Return the MAP (maximum a posteriori) estimate = the means.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def default(self) -> CouplingPrior:
    """Return the MAP (maximum a posteriori) estimate = the means.

    Returns
    -------
    CouplingPrior
        Return the MAP (maximum a posteriori) estimate = the means.
    """
    return CouplingPrior(
        K_base=self._K_base_mean,
        decay_alpha=self._decay_alpha_mean,
        K_c_estimate=0.0,
    )

estimate_Kc

estimate_Kc(
    omegas: FloatArray, n_layers: int
) -> CouplingPrior

Combine prior with Dörfler-Bullo K_c for given omegas.

K_c = max|ω_i - ω_j| / λ₂(L) where L is built from the prior's decay_alpha on a chain graph of n_layers.

Parameters

omegas : FloatArray Natural frequencies in rad/s, shape (N,). n_layers : int Number of SCPN hierarchy layers.

Returns

CouplingPrior The prior combined with the Dörfler-Bullo K_c for the given frequencies.

Raises

TypeError If an argument has the wrong type. ValueError If omegas or the layer count is invalid.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def estimate_Kc(self, omegas: FloatArray, n_layers: int) -> CouplingPrior:
    """Combine prior with Dörfler-Bullo K_c for given omegas.

    K_c = max|ω_i - ω_j| / λ₂(L) where L is built from the prior's
    decay_alpha on a chain graph of n_layers.

    Parameters
    ----------
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    n_layers : int
        Number of SCPN hierarchy layers.

    Returns
    -------
    CouplingPrior
        The prior combined with the Dörfler-Bullo ``K_c`` for the given frequencies.

    Raises
    ------
    TypeError
        If an argument has the wrong type.
    ValueError
        If ``omegas`` or the layer count is invalid.
    """
    if isinstance(n_layers, bool) or not isinstance(n_layers, Integral):
        raise TypeError(f"n_layers must be an integer, got {n_layers!r}")
    n_layers = int(n_layers)
    if n_layers <= 0:
        raise ValueError("n_layers must be a positive integer")
    from scpn_phase_orchestrator.coupling.spectral import critical_coupling

    prior = self.default()
    omega_values = _validate_frequency_vector(omegas)
    if omega_values.size != n_layers:
        raise ValueError("n_layers must match the length of omegas")
    idx = np.arange(n_layers)
    dist = np.abs(idx[:, np.newaxis] - idx[np.newaxis, :])
    knm: FloatArray = prior.K_base * np.exp(-prior.decay_alpha * dist)
    np.fill_diagonal(knm, 0.0)
    K_c = critical_coupling(omega_values, knm)
    return CouplingPrior(
        K_base=prior.K_base,
        decay_alpha=prior.decay_alpha,
        K_c_estimate=K_c,
    )

log_probability

log_probability(K_base: float, decay_alpha: float) -> float

Log-probability under the Gaussian prior (unnormalised).

Parameters

K_base : float Base coupling strength before spatial modulation. decay_alpha : float Exponential decay rate of the coupling across layer separation.

Returns

float The unnormalised log-probability of the configuration under the prior.

Raises

TypeError If an argument has the wrong type. ValueError If K_base or decay_alpha is out of range.

Source code in src/scpn_phase_orchestrator/coupling/prior.py

def log_probability(self, K_base: float, decay_alpha: float) -> float:
    """Log-probability under the Gaussian prior (unnormalised).

    Parameters
    ----------
    K_base : float
        Base coupling strength before spatial modulation.
    decay_alpha : float
        Exponential decay rate of the coupling across layer separation.

    Returns
    -------
    float
        The unnormalised log-probability of the configuration under the prior.

    Raises
    ------
    TypeError
        If an argument has the wrong type.
    ValueError
        If ``K_base`` or ``decay_alpha`` is out of range.
    """
    if isinstance(K_base, bool) or not isinstance(K_base, Real):
        raise TypeError("K_base must be a finite real value")
    if isinstance(decay_alpha, bool) or not isinstance(decay_alpha, Real):
        raise TypeError("decay_alpha must be a finite real value")
    K_base = float(K_base)
    decay_alpha = float(decay_alpha)
    if not np.isfinite(K_base) or not np.isfinite(decay_alpha):
        raise ValueError("K_base and decay_alpha must be finite real values")
    if _HAS_RUST:
        log_prob = float(
            _rust_log_prob(
                K_base,
                decay_alpha,
                self._K_base_mean,
                self._K_base_std,
                self._decay_alpha_mean,
                self._decay_alpha_std,
            )
        )
        if not np.isfinite(log_prob):
            raise ValueError(
                "Rust prior log-probability must be a finite real value"
            )
        return log_prob
    lp_K = -0.5 * ((K_base - self._K_base_mean) / self._K_base_std) ** 2
    lp_a = (
        -0.5 * ((decay_alpha - self._decay_alpha_mean) / self._decay_alpha_std) ** 2
    )
    return lp_K + lp_a

---

HCP Connectome Generator

Neuroscience-realistic coupling matrices.

Synthetic generator

load_hcp_connectome(n_regions, seed=42) generates a matrix with:

• Intra-hemispheric: exponential distance decay
• Inter-hemispheric: corpus callosum pattern (homotopic connections)
• Default Mode Network: hub structure with elevated coupling

Real data bridge

load_neurolib_hcp(n_regions=80) loads real HCP structural connectivity from the neurolib library. Supports n_regions from 2 to 80.

Both loaders validate optional-backend matrices before publication. Boolean, complex/object-complex, and numeric-string aliases are rejected before float64 conversion; finite real numeric-object matrices remain compatible. Shape, finiteness, non-negativity, symmetry, and zero-diagonal constraints are then replayed at the public Python boundary.

Performance: load_hcp_connectome(80) < 10 ms (Python), ~48 µs (Rust, 17.6x speedup). Detailed documentation: HCP Connectome — detailed reference

connectome

Synthetic and optional neurolib HCP coupling loaders.

load_hcp_connectome generates a deterministic HCP-inspired synthetic matrix with explicit non-real-data provenance. load_neurolib_hcp is the optional real HCP path and fails with an import error when neurolib is unavailable. Both paths return non-negative zero-diagonal structural coupling matrices suitable for examples, validation, and explicit downstream review.

Functions:

load_neurolib_hcp

load_neurolib_hcp(n_regions: int = 80) -> FloatArray

Load real HCP structural connectivity from neurolib.

Parameters

n_regions : int number of regions to return (max 80). If < 80, returns the top-left (n_regions, n_regions) submatrix.

Returns

FloatArray Symmetric non-negative coupling matrix, shape (n_regions, n_regions).

Raises

ImportError If neurolib is not installed. ValueError If n_regions < 2 or > 80.

Source code in src/scpn_phase_orchestrator/coupling/connectome.py

def load_neurolib_hcp(n_regions: int = 80) -> FloatArray:
    """Load real HCP structural connectivity from neurolib.

    Parameters
    ----------
    n_regions : int
        number of regions to return (max 80). If < 80, returns the top-left (n_regions,
        n_regions) submatrix.

    Returns
    -------
    FloatArray
        Symmetric non-negative coupling matrix, shape (n_regions, n_regions).

    Raises
    ------
    ImportError
        If neurolib is not installed.
    ValueError
        If n_regions < 2 or > 80.
    """
    try:
        # type ignore: neurolib is optional and currently lacks complete type metadata;
        # runtime availability is handled by the ModuleNotFoundError branch.
        from neurolib.utils.loadData import (  # type: ignore[import-untyped,import-not-found]
            Dataset,
        )
    except ModuleNotFoundError:
        raise ImportError(
            "neurolib is required for real HCP data: pip install neurolib"
        ) from None

    n_regions = _validate_n_regions(n_regions, max_regions=_NEUROLIB_HCP_SIZE)

    ds = Dataset("hcp")
    sc = _coerce_connectome_matrix(
        ds.Cmat,
        n_regions=_NEUROLIB_HCP_SIZE,
        source="neurolib HCP",
    )[:n_regions, :n_regions].copy()
    np.fill_diagonal(sc, 0.0)
    return _validate_connectome_matrix(
        sc,
        n_regions=n_regions,
        source="neurolib HCP",
    )

load_hcp_connectome

load_hcp_connectome(
    n_regions: int, seed: int = 42
) -> FloatArray

Generate a synthetic HCP-inspired coupling matrix.

Parameters

n_regions : int number of cortical regions (must be >= 2, even recommended).

Returns

FloatArray Symmetric coupling matrix, shape (n_regions, n_regions), zero diagonal.

Source code in src/scpn_phase_orchestrator/coupling/connectome.py

def load_hcp_connectome(n_regions: int, seed: int = 42) -> FloatArray:
    """Generate a synthetic HCP-inspired coupling matrix.

    Parameters
    ----------
    n_regions : int
        number of cortical regions (must be >= 2, even recommended).

    Returns
    -------
    FloatArray
        Symmetric coupling matrix, shape (n_regions, n_regions), zero diagonal.
    """
    n_regions = _validate_n_regions(n_regions)
    seed = _validate_seed(seed)
    backend = _rust_load_hcp if _HAS_RUST else None
    return _load_hcp_connectome_cached(n_regions, seed, _HAS_RUST, backend).copy()

---

Rust FFI acceleration

spo_kernel.PyCouplingBuilder provides Rust-accelerated K_nm construction. The Python implementation is the reference; the Rust path is selected automatically when spo_kernel is importable. Parity is verified in tests/test_rust_python_parity_performance.py.

Rust builder returns are inspected before numeric conversion: boolean, complex/object-complex, and numeric-string K_nm or alpha aliases are rejected and trigger the documented NumPy fallback. Finite real numeric-object matrices remain compatible; shape, finiteness, non-negativity, symmetry, and zero- diagonal checks still run before publication.

Performance summary

Operation	Budget	Measured
CouplingBuilder.build(100)	< 10 ms	~2 ms
build_scpn_physics()	< 5 ms	~1 ms
estimate_from_distances(64)	< 5 ms	~0.5 ms
load_hcp_connectome(80)	< 10 ms	~3 ms
validate_knm(64)	< 1 ms	~0.1 ms
graph_laplacian(64)	< 1 ms	~0.007 ms
fiedler_value(64)	< 1 ms	~0.12 ms

Spatial coupling modulation

SpatialCouplingModulator is the public PHA-C.1 coupling surface for systems where the effective phase coupling must depend on moving geometry instead of static oscillator labels. It turns a zero-diagonal base K_nm matrix and a position matrix into a physically constrained modulated coupling matrix.

Use it when spatial proximity, mobile agents, tissue geometry, sensor placement, or edge-node distance changes the strength of phase transfer. The default kernel is 1 / (1 + distance), which is bounded, finite at zero separation, symmetric for Euclidean positions, and preserves the zero self-coupling diagonal required by the oscillator engines.

The module also exposes exponential, power-law, and inverse-distance kernels. The inverse-distance form is reserved for Swarmalator compatibility and uses an epsilon-regularised denominator so the historical kernel remains bit-true without introducing singularities.

The reference implementation is NumPy. Rust, Go, Julia, and Mojo adapters are validated as optional accelerators and must reproduce the same invariants before their output is accepted: finite real-valued matrices, exact shape or flat cardinality, non-boolean and non-complex values, non-negative entries, zero diagonal, and symmetry preservation for symmetric inputs. Public positions, base coupling matrices, scalar decay controls, direct accelerator counts/forms/flat buffers, optional backend outputs, and raw Julia returns reject numeric-string aliases before float coercion. The public dispatcher preserves matrix-shaped output for callers after replaying the shared direct output validator; optional backend fallback remains limited to loader or runtime unavailability.

See Coupling - Spatial Modulator for examples, backend notes, and the benchmark contract.