Analysis Toolkit

Purpose and scope

The toolkit is structured for operational readability: each monitor returns a different failure or regime signal before the scalar order parameter alone would show anything unusual. In practice, teams use this page as a first-pass selection guide, then tune thresholds against their domain trajectories.

How operators should use this layer

Treat the toolkit as a multi-signal diagnostic funnel rather than a single alarm source:

• start with fast, broad monitors (order_parameter, lyapunov);
• add structural monitors (plv, chimera, winding) when patterns localise;
• then confirm causal or thermodynamic interpretations (coupling_est, entropy_prod, itpc, pid).

This sequencing reduces false positives and gives policy teams a reproducible rationale for each escalation before any actuation change is promoted.

The sections below are ordered from global coherence to higher-order coupling relationships because that mirrors a typical diagnostic flow: global stability -> phase alignment structure -> causality and stability risk -> topological and thermodynamic drift.

SPO provides 12 dynamical monitors — most oscillator simulators have 1-2. Each monitor detects a different aspect of the dynamics that scalar R misses.

Selecting a minimal monitor set

For a first production run, teams usually start with:

• order_parameter for synchronization trend and collapse detection,
• lyapunov for local stability margin,
• plv for pairwise synchrony topology,
• one supervisory metric (evs or pid) for domain-facing interpretability.

Expanding to all monitors is recommended only after a baseline is stable; this keeps false alarm fatigue manageable while keeping observability depth.

Order Parameter & PLV

Standard Kuramoto order parameter R = |⟨exp(iθ)⟩| and Phase-Locking Value matrix PLV_ij = |⟨exp(i(θ_i - θ_j))⟩_t|.

order_params

Kuramoto order parameter family with 5-backend fallback chain.

Follows the AttnRes-level module standard (feedback_module_standard_attnres.md):

• compute_order_parameter — R and mean phase ψ.
• compute_plv — phase-locking value between two equal-length phase series.
• compute_layer_coherence — R restricted to a layer.

Each kernel is available in five languages — Rust, Mojo, Julia, Go, Python. AVAILABLE_BACKENDS reports detected backends in canonical fallback order, while ACTIVE_BACKEND is selected by a small import-time hot-path probe so slow external wrappers do not displace the faster local path.

Functions:

compute_order_parameter

compute_order_parameter(
    phases: FloatArray,
) -> tuple[float, float]

Kuramoto global order parameter (R, ψ).

R = |mean(exp(i · θ))|; ψ = arg(mean(exp(i · θ))) mod 2π.

Parameters

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

Returns

tuple[float, float] The (R, ψ) Kuramoto order parameter and mean phase.

Notes

R carries a positive small-sample bias: N uniformly random phases give E[R²] = 1/N, so R reads about 1/√N (≈ 0.35 at N = 8) even with no coherence. A monitor comparing coherence across small populations should use :func:debiased_squared_order_parameter, whose expectation is 0 under uniformity, rather than reading the raw R as if it were unbiased.

Source code in src/scpn_phase_orchestrator/upde/order_params.py

def compute_order_parameter(phases: FloatArray) -> tuple[float, float]:
    """Kuramoto global order parameter ``(R, ψ)``.

    ``R = |mean(exp(i · θ))|``;
    ``ψ = arg(mean(exp(i · θ))) mod 2π``.

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

    Returns
    -------
    tuple[float, float]
        The ``(R, ψ)`` Kuramoto order parameter and mean phase.

    Notes
    -----
    ``R`` carries a positive small-sample bias: ``N`` uniformly random phases give
    ``E[R²] = 1/N``, so ``R`` reads about ``1/√N`` (≈ 0.35 at ``N = 8``) even with
    no coherence. A monitor comparing coherence across small populations should use
    :func:`debiased_squared_order_parameter`, whose expectation is 0 under
    uniformity, rather than reading the raw ``R`` as if it were unbiased.
    """
    phases = _validate_phases("phases", phases)
    if phases.size == 0:
        return (0.0, 0.0)
    backend_fn = _dispatch("order_parameter")
    if backend_fn is not None:
        fn = cast("Callable[[FloatArray], tuple[float, float]]", backend_fn)
        p = np.ascontiguousarray(phases.ravel(), dtype=np.float64)
        r, psi = fn(p)
        return (
            validate_unit_interval_output(r, name="coherence magnitude"),
            validate_mean_phase_output(psi),
        )

    return _python_order_parameter(phases)

debiased_squared_order_parameter

debiased_squared_order_parameter(
    phases: FloatArray,
) -> float

Return the small-N-bias-corrected squared Kuramoto order parameter.

The raw magnitude R = |mean(exp(iθ))| has a positive small-sample bias: N uniformly random phases give E[R²] = 1/N, so R reads about 1/√N (≈ 0.35 at N = 8) even with no coherence. This estimator removes that floor. It is the pairwise phase consistency of Vinck et al. (2010),

(N · R² − 1) / (N − 1),

whose expectation is 0 under uniform phases and 1 at perfect synchrony. It can be slightly negative on a finite anti-aligned sample — that is honest, not an error. Reuses the accelerated :func:compute_order_parameter for R; the debiasing itself is a scalar correction.

Parameters

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

Returns

float The debiased squared order parameter, in [-1 / (N − 1), 1].

Raises

ValueError If fewer than two phases are supplied — the correction is undefined for a single oscillator (N − 1 = 0).

References

Vinck, M., van Wingerden, M., Womelsdorf, T., Fries, P., & Pennartz, C. M. A. (2010). The pairwise phase consistency: a bias-free measure of rhythmic neuronal synchronization. NeuroImage, 51(1), 112–122.

Source code in src/scpn_phase_orchestrator/upde/order_params.py

def debiased_squared_order_parameter(phases: FloatArray) -> float:
    """Return the small-N-bias-corrected squared Kuramoto order parameter.

    The raw magnitude ``R = |mean(exp(iθ))|`` has a positive small-sample bias:
    ``N`` uniformly random phases give ``E[R²] = 1/N``, so ``R`` reads about
    ``1/√N`` (≈ 0.35 at ``N = 8``) even with no coherence. This estimator removes
    that floor. It is the pairwise phase consistency of Vinck et al. (2010),

    ``(N · R² − 1) / (N − 1)``,

    whose expectation is 0 under uniform phases and 1 at perfect synchrony. It can
    be slightly negative on a finite anti-aligned sample — that is honest, not an
    error. Reuses the accelerated :func:`compute_order_parameter` for ``R``; the
    debiasing itself is a scalar correction.

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

    Returns
    -------
    float
        The debiased squared order parameter, in ``[-1 / (N − 1), 1]``.

    Raises
    ------
    ValueError
        If fewer than two phases are supplied — the correction is undefined for a
        single oscillator (``N − 1 = 0``).

    References
    ----------
    Vinck, M., van Wingerden, M., Womelsdorf, T., Fries, P., & Pennartz, C. M. A.
    (2010). The pairwise phase consistency: a bias-free measure of rhythmic
    neuronal synchronization. *NeuroImage*, 51(1), 112–122.
    """
    validated = _validate_phases("phases", phases)
    n = int(validated.size)
    if n < 2:
        raise ValueError("debiased order parameter requires at least two phases")
    coherence, _ = compute_order_parameter(validated)
    return float((n * coherence * coherence - 1.0) / (n - 1))

compute_plv

compute_plv(
    phases_a: FloatArray, phases_b: FloatArray
) -> float

Phase-locking value between two equal-length phase series.

PLV = |mean(exp(i · (φ_a − φ_b)))| over samples.

Parameters

phases_a : FloatArray First phase series in radians, shape (T,). phases_b : FloatArray Second phase series in radians, shape (T,).

Returns

float The phase-locking value between the two series.

Raises

ValueError If the two phase series have different lengths.

Source code in src/scpn_phase_orchestrator/upde/order_params.py

def compute_plv(phases_a: FloatArray, phases_b: FloatArray) -> float:
    """Phase-locking value between two equal-length phase series.

    PLV = ``|mean(exp(i · (φ_a − φ_b)))|`` over samples.

    Parameters
    ----------
    phases_a : FloatArray
        First phase series in radians, shape ``(T,)``.
    phases_b : FloatArray
        Second phase series in radians, shape ``(T,)``.

    Returns
    -------
    float
        The phase-locking value between the two series.

    Raises
    ------
    ValueError
        If the two phase series have different lengths.
    """
    phases_a = _validate_phases("phases_a", phases_a)
    phases_b = _validate_phases("phases_b", phases_b)
    if phases_a.size != phases_b.size:
        raise ValueError(
            f"PLV requires equal-length arrays, got {phases_a.size} vs {phases_b.size}"
        )
    if phases_a.size == 0:
        return 0.0
    backend_fn = _dispatch("plv")
    if backend_fn is not None:
        fn = cast(
            "Callable[[FloatArray, FloatArray], float]",
            backend_fn,
        )
        a = np.ascontiguousarray(phases_a.ravel(), dtype=np.float64)
        b = np.ascontiguousarray(phases_b.ravel(), dtype=np.float64)
        return validate_unit_interval_output(fn(a, b), name="coherence magnitude")

    return validate_unit_interval_output(
        float(np.abs(np.mean(np.exp(1j * (phases_a - phases_b))))),
        name="coherence magnitude",
    )

compute_layer_coherence

compute_layer_coherence(
    phases: FloatArray, layer_mask: BoolArray | IntArray
) -> float

Return the order parameter R for the oscillators in layer_mask.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). layer_mask : BoolArray | IntArray Boolean mask or integer index array selecting the layer's oscillators.

Returns

float The Kuramoto order parameter R for the selected oscillators.

Source code in src/scpn_phase_orchestrator/upde/order_params.py

def compute_layer_coherence(
    phases: FloatArray, layer_mask: BoolArray | IntArray
) -> float:
    """Return the order parameter R for the oscillators in ``layer_mask``.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    layer_mask : BoolArray | IntArray
        Boolean mask or integer index array selecting the layer's oscillators.

    Returns
    -------
    float
        The Kuramoto order parameter ``R`` for the selected oscillators.
    """
    phases = _validate_phases("phases", phases)
    if phases.size == 0:
        return 0.0
    indices = _layer_indices(layer_mask, phases.size)
    if indices.size == 0:
        return 0.0
    backend_fn = _dispatch("layer_coherence")
    if backend_fn is not None:
        fn = cast("Callable[[FloatArray, IntArray], float]", backend_fn)
        p = np.ascontiguousarray(phases.ravel(), dtype=np.float64)
        return validate_unit_interval_output(
            fn(p, indices),
            name="coherence magnitude",
        )

    # ``indices`` is non-empty (checked above) and every entry is a valid
    # in-range oscillator, so the gathered sub-population is never empty.
    sub = phases[indices]
    z = np.mean(np.exp(1j * sub))
    return validate_unit_interval_output(
        float(np.abs(z)),
        name="coherence magnitude",
    )

Phase-Amplitude Coupling (PAC)

Modulation index (MI) via Tort et al. 2010. Bins low-frequency phase, computes mean amplitude per bin, KL divergence from uniform. N×N PAC matrix: entry [i,j] = MI(phase_i, amplitude_j).

Central to neuroscience — cross-frequency coupling between brain oscillation bands (theta-gamma, alpha-beta).

pac

Tort 2010 phase-amplitude coupling with 5-backend fallback chain.

Follows feedback_module_standard_attnres.md:

• modulation_index — scalar Tort 2010 MI on a single (θ_low, a_high) pair of time series.
• pac_matrix — (N, N) pairwise MI matrix over N oscillator phase / amplitude channels.
• pac_gate — pure-Python boolean gate on an MI value (no backend dispatch needed — trivial comparison).

All compute kernels are available in Rust, Mojo, Julia, Go, Python. AVAILABLE_BACKENDS reports detected backends in canonical fallback order, while ACTIVE_BACKEND is selected by a small import-time hot-path probe so slow external wrappers do not displace the faster local path.

Functions:

modulation_index

modulation_index(
    theta_low: FloatArray,
    amp_high: FloatArray,
    n_bins: int = 18,
) -> float

Phase-amplitude coupling via Tort et al. 2010, J. Neurophysiol.

Bins amplitude by phase, computes KL divergence from uniform, returns the modulation index normalised to [0, 1] by log(n_bins).

Parameters

theta_low : FloatArray Low-frequency driver phase in radians, shape (T,). amp_high : FloatArray High-frequency amplitude envelope, shape (T,). n_bins : int Number of phase bins used for the modulation-index histogram.

Returns

float The Tort modulation index of phase-amplitude coupling.

Raises

ValueError If n_bins is not a positive integer or inputs mismatch.

Source code in src/scpn_phase_orchestrator/upde/pac.py

def modulation_index(
    theta_low: FloatArray, amp_high: FloatArray, n_bins: int = 18
) -> float:
    """Phase-amplitude coupling via Tort et al. 2010, J. Neurophysiol.

    Bins amplitude by phase, computes KL divergence from uniform,
    returns the modulation index normalised to ``[0, 1]`` by
    ``log(n_bins)``.

    Parameters
    ----------
    theta_low : FloatArray
        Low-frequency driver phase in radians, shape ``(T,)``.
    amp_high : FloatArray
        High-frequency amplitude envelope, shape ``(T,)``.
    n_bins : int
        Number of phase bins used for the modulation-index histogram.

    Returns
    -------
    float
        The Tort modulation index of phase-amplitude coupling.

    Raises
    ------
    ValueError
        If ``n_bins`` is not a positive integer or inputs mismatch.
    """
    n_bins = _validate_n_bins(n_bins)
    theta_low = _validate_signal("theta_low", theta_low)
    amp_high = _validate_signal("amp_high", amp_high)
    if np.any(amp_high < 0.0):
        raise ValueError("amp_high must contain non-negative amplitudes")
    if theta_low.size == 0 or amp_high.size == 0:
        return 0.0

    backend_fn = _dispatch("modulation_index")
    if backend_fn is not None:
        fn = cast(
            "Callable[[FloatArray, FloatArray, int], float]",
            backend_fn,
        )
        return _validate_mi_value(
            "modulation_index backend",
            fn(
                np.ascontiguousarray(theta_low, dtype=np.float64),
                np.ascontiguousarray(amp_high, dtype=np.float64),
                n_bins,
            ),
        )

    return _validate_mi_value(
        "modulation_index",
        _modulation_index_python(theta_low, amp_high, n_bins),
    )

pac_matrix

pac_matrix(
    phases_history: FloatArray,
    amplitudes_history: FloatArray,
    n_bins: int = 18,
) -> FloatArray

Return the (N, N) PAC matrix [i, j] = MI(phase_i, amplitude_j).

Parameters

phases_history : FloatArray (T, N) phase time series. amplitudes_history : FloatArray (T, N) amplitude time series. n_bins : int number of phase bins.

Returns

FloatArray FloatArray The (N, N) phase-amplitude coupling matrix.

Raises

ValueError If n_bins is not positive or the histories have mismatched shapes.

Source code in src/scpn_phase_orchestrator/upde/pac.py

def pac_matrix(
    phases_history: FloatArray,
    amplitudes_history: FloatArray,
    n_bins: int = 18,
) -> FloatArray:
    """Return the ``(N, N)`` PAC matrix ``[i, j] = MI(phase_i, amplitude_j)``.

    Parameters
    ----------
    phases_history : FloatArray
        ``(T, N)`` phase time series.
    amplitudes_history : FloatArray
        ``(T, N)`` amplitude time series.
    n_bins : int
        number of phase bins.

    Returns
    -------
    FloatArray
        FloatArray The ``(N, N)`` phase-amplitude coupling matrix.

    Raises
    ------
    ValueError
        If ``n_bins`` is not positive or the histories have mismatched shapes.
    """
    n_bins = _validate_n_bins(n_bins)
    phases_history = _validate_history("phases_history", phases_history)
    amplitudes_history = _validate_history("amplitudes_history", amplitudes_history)
    if np.any(amplitudes_history < 0.0):
        raise ValueError("amplitudes_history must contain non-negative amplitudes")
    t, n = phases_history.shape
    if t == 0 or n == 0:
        return np.zeros((n, n), dtype=np.float64)
    if amplitudes_history.shape != (t, n):
        raise ValueError("phases and amplitudes must have the same shape")

    backend_fn = _dispatch("pac_matrix")
    if backend_fn is not None:
        fn = cast(
            ("Callable[[FloatArray, FloatArray, int, int, int], FloatArray]"),
            backend_fn,
        )
        flat = fn(
            np.ascontiguousarray(phases_history.ravel(order="C"), dtype=np.float64),
            np.ascontiguousarray(amplitudes_history.ravel(order="C"), dtype=np.float64),
            t,
            n,
            n_bins,
        )
        matrix = validate_pac_matrix_output(flat, n=n)
        return matrix.reshape((n, n), order="C")

    result = np.zeros((n, n), dtype=np.float64)
    for i in range(n):
        for j in range(n):
            result[i, j] = modulation_index(
                phases_history[:, i], amplitudes_history[:, j], n_bins
            )
    return result

pac_gate

pac_gate(pac_value: float, threshold: float = 0.3) -> bool

Binary gate: True when PAC exceeds threshold.

Pure-Python helper; no dispatcher — the comparison is trivial.

Parameters

pac_value : float A phase-amplitude coupling value. threshold : float Decision threshold.

Returns

bool True when the PAC value exceeds the threshold.

Source code in src/scpn_phase_orchestrator/upde/pac.py

def pac_gate(pac_value: float, threshold: float = 0.3) -> bool:
    """Binary gate: ``True`` when PAC exceeds ``threshold``.

    Pure-Python helper; no dispatcher — the comparison is trivial.

    Parameters
    ----------
    pac_value : float
        A phase-amplitude coupling value.
    threshold : float
        Decision threshold.

    Returns
    -------
    bool
        ``True`` when the PAC value exceeds the threshold.
    """
    pac_value = _validate_finite_real("pac_value", pac_value)
    threshold = _validate_finite_real("threshold", threshold)
    return pac_value >= threshold

Chimera State Detection

Detects chimera states: coexisting coherent and incoherent clusters within the same network. Uses local order parameter R_i based on neighborhood coupling.

• Coherent: R_i > 0.7
• Incoherent: R_i < 0.3
• Boundary: in-between
• Chimera index = boundary_count / N

Detects phase transitions that global R misses.

chimera

Chimera state detection with a 5-backend fallback chain.

Kuramoto & Battogtokh 2002, Nonlinear Phenomena in Complex Systems 5:380–385. An oscillator i is coherent when its local order parameter R_i = |⟨exp(i(θ_j − θ_i))⟩_{j ∈ N(i)}| exceeds the coherence threshold, incoherent when it falls below the incoherence threshold. The chimera index is the fraction of oscillators that sit in the boundary band in between.

Compute surface:

• :func:local_order_parameter — (N,) per-oscillator R_i vector; the coupling diagonal must be zero so self-coupling is never counted as a neighbour.
• :func:detect_chimera — classification wrapper returning :class:ChimeraState.

Entrainment Verification Score (EVS)

Three-criterion battery for rigorous entrainment validation:

1. ITPC (inter-trial phase coherence) persistence
2. Survival during stimulus pause
3. Frequency specificity (ratio at target vs control frequency)

Distinguishes true entrainment from broadband phase-locking artifacts.

EVS and phase-locking metrics for finite two-dimensional phase recordings.

The module implements ITPC, persistence across pauses, and frequency-specificity checks for Entrainment Verification Signals. A Rust extension is used when available while the Python fallback remains the reference-compatible path. Inputs are normalized to finite trials x time phase arrays, pause indices are bounds-checked, and candidate frequency vectors must match the trial axis before evidence is reported.

Partial Information Decomposition (PID)

Decomposes mutual information into:

• Redundancy: shared information from both oscillator groups
• Synergy: information present only in the joint group

Detects when groups carry synergistic (non-redundant) information about global phase (Williams & Beer 2010).

pid

Partial information decomposition (PID) about global synchronisation.

Decomposes two oscillator groups with a 5-backend fallback chain.

Model

Williams & Beer 2010 (Nonnegative Decomposition of Multivariate Information, arXiv:1004.2515) decompose the information two sources carry about a target into redundant, unique, and synergistic parts. Estimating it needs a distribution, so the input is a phase history (T, N) (T timesteps, N oscillators). Each timestep is reduced to three circular observables:

• target Y_t — the global order-parameter phase ∠⟨e^{iθ}⟩ over all oscillators,
• source A_t — the group-A order-parameter phase,
• source B_t — the group-B order-parameter phase.

The three series are binned into n_bins equal-width phase bins and the joint distribution is estimated over the T samples.

Decomposition

With the specific information I_spec(Y=y; S) = Σ_s p(s|y)·log[p(y|s)/p(y)]:

redundancy  I_red = Σ_y p(y)·min( I_spec(Y=y; A), I_spec(Y=y; B) )
synergy     I_syn = MI(A,B; Y) − MI(A; Y) − MI(B; Y) + I_red

I_red is the Williams & Beer I_min redundancy; the unique information of each source is MI(S; Y) − I_red and MI(A; Y) = I_red + U_A holds by construction. All terms are non-negative.

A single snapshot (T = 1) carries no distributional information, so every component is 0; meaningful decomposition needs T ≥ 2.

Functions:

redundancy

redundancy(
    phases: FloatArray,
    group_a: list[int] | IntArray,
    group_b: list[int] | IntArray,
    n_bins: int = _DEFAULT_BINS,
) -> float

Redundant information both groups share about the global phase.

I_red = Σ_y p(y)·min(I_spec(Y=y; A), I_spec(Y=y; B)) (Williams & Beer 2010 I_min). phases is a (T, N) phase history.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). group_a : list[int] | IntArray Indices of the first oscillator group. group_b : list[int] | IntArray Indices of the second oscillator group. n_bins : int Number of histogram bins.

Returns

float The redundant information the groups share about the global phase.

Source code in src/scpn_phase_orchestrator/monitor/pid.py

def redundancy(
    phases: FloatArray,
    group_a: list[int] | IntArray,
    group_b: list[int] | IntArray,
    n_bins: int = _DEFAULT_BINS,
) -> float:
    """Redundant information both groups share about the global phase.

    ``I_red = Σ_y p(y)·min(I_spec(Y=y; A), I_spec(Y=y; B))`` (Williams & Beer
    2010 ``I_min``). ``phases`` is a ``(T, N)`` phase history.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    group_a : list[int] | IntArray
        Indices of the first oscillator group.
    group_b : list[int] | IntArray
        Indices of the second oscillator group.
    n_bins : int
        Number of histogram bins.

    Returns
    -------
    float
        The redundant information the groups share about the global phase.
    """
    red, _ = _decompose(phases, group_a, group_b, n_bins)
    return _validate_pid_scalar(red, name="redundancy")

synergy

synergy(
    phases: FloatArray,
    group_a: list[int] | IntArray,
    group_b: list[int] | IntArray,
    n_bins: int = _DEFAULT_BINS,
) -> float

Synergistic information present only in the joint (A, B).

I_syn = MI(A,B; Y) − MI(A; Y) − MI(B; Y) + I_red. Positive synergy means the combined group carries information about the global state that neither subgroup carries alone. phases is a (T, N) phase history.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). group_a : list[int] | IntArray Indices of the first oscillator group. group_b : list[int] | IntArray Indices of the second oscillator group. n_bins : int Number of histogram bins.

Returns

float The synergistic information present only in the joint (A, B).

Source code in src/scpn_phase_orchestrator/monitor/pid.py

def synergy(
    phases: FloatArray,
    group_a: list[int] | IntArray,
    group_b: list[int] | IntArray,
    n_bins: int = _DEFAULT_BINS,
) -> float:
    """Synergistic information present only in the joint ``(A, B)``.

    ``I_syn = MI(A,B; Y) − MI(A; Y) − MI(B; Y) + I_red``. Positive synergy means
    the combined group carries information about the global state that neither
    subgroup carries alone. ``phases`` is a ``(T, N)`` phase history.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    group_a : list[int] | IntArray
        Indices of the first oscillator group.
    group_b : list[int] | IntArray
        Indices of the second oscillator group.
    n_bins : int
        Number of histogram bins.

    Returns
    -------
    float
        The synergistic information present only in the joint ``(A, B)``.
    """
    _, syn = _decompose(phases, group_a, group_b, n_bins)
    return _validate_pid_scalar(syn, name="synergy")

Lyapunov Exponent

Real-time estimation of the maximal Lyapunov exponent. Positive = chaos, zero = edge of chaos (critical), negative = stable attractor.

lyapunov

Lyapunov stability monitor with a 5-backend fallback chain.

Two public surfaces:

• :class:LyapunovGuard — stateful observer that tracks the Lyapunov function V(θ) = -(K/2N) Σ_ij A_ij cos(θ_i − θ_j), its numerical time derivative, and basin-of-attraction membership (van Hemmen & Wreszinski 1993). Single-backend NumPy; inexpensive per call.
• :func:lyapunov_spectrum — full Lyapunov spectrum via periodic QR reorthogonalisation (Benettin 1980 / Shimada-Nagashima 1979). Multi- backend; the heavy kernel is dispatched to Rust → Mojo → Julia → Go → Python in order of availability.

Classes

LyapunovState dataclass

LyapunovState(
    V: float,
    dV_dt: float,
    in_basin: bool,
    max_phase_diff: float,
)

Lyapunov function V, dV/dt, basin membership, and max phase diff.

Methods:

__post_init__

__post_init__() -> None

Normalize scalar aliases and reject invalid state fields.

Source code in src/scpn_phase_orchestrator/monitor/lyapunov.py

def __post_init__(self) -> None:
    """Normalize scalar aliases and reject invalid state fields."""
    v_value = _validate_finite_real(self.V, name="V")
    dv_dt_value = _validate_finite_real(self.dV_dt, name="dV_dt")
    if not isinstance(self.in_basin, bool):
        raise ValueError(f"in_basin must be a boolean flag, got {self.in_basin!r}")
    max_phase_diff = _validate_non_negative_real(
        self.max_phase_diff,
        name="max_phase_diff",
    )
    if max_phase_diff > np.pi:
        raise ValueError(
            f"max_phase_diff must be <= pi for geodesic phase distance, "
            f"got {self.max_phase_diff!r}"
        )

    self.V = v_value
    self.dV_dt = dv_dt_value
    self.max_phase_diff = max_phase_diff

LyapunovGuard

LyapunovGuard(basin_threshold: object = np.pi / 2)

Lyapunov stability monitor for Kuramoto networks.

V(θ) = -(K/2N) Σ_{i,j} A_ij cos(θ_i - θ_j)

dV/dt ≤ 0 for gradient flow (Kuramoto is gradient on V). Basin of attraction: max|θ_i - θ_j| < π/2 for connected pairs.

van Hemmen & Wreszinski 1993, J. Stat. Phys. 72:145-166.

Create a guard with a validated geodesic basin threshold.

Source code in src/scpn_phase_orchestrator/monitor/lyapunov.py

def __init__(self, basin_threshold: object = np.pi / 2):
    """Create a guard with a validated geodesic basin threshold."""
    basin_threshold = _validate_positive_real(
        basin_threshold,
        name="basin_threshold",
    )
    if basin_threshold > np.pi:
        raise ValueError(
            f"basin_threshold must be <= pi for geodesic phase distance, "
            f"got {basin_threshold!r}"
        )
    self._basin_threshold = basin_threshold
    self._prev_V: float | None = None

Methods:

evaluate

evaluate(phases: object, knm: object) -> LyapunovState

Compute Lyapunov function, its time derivative, and basin check.

Parameters

phases : object Oscillator phases in radians, shape (N,). knm : object Coupling matrix K_nm, shape (N, N).

Returns

LyapunovState The Lyapunov value, its derivative, and the basin-check result.

Source code in src/scpn_phase_orchestrator/monitor/lyapunov.py

def evaluate(self, phases: object, knm: object) -> LyapunovState:
    """Compute Lyapunov function, its time derivative, and basin check.

    Parameters
    ----------
    phases : object
        Oscillator phases in radians, shape ``(N,)``.
    knm : object
        Coupling matrix ``K_nm``, shape ``(N, N)``.

    Returns
    -------
    LyapunovState
        The Lyapunov value, its derivative, and the basin-check result.
    """
    phase_values = _validate_vector(phases, name="phases")
    n = len(phase_values)
    coupling_matrix = _validate_matrix(knm, name="knm", expected_shape=(n, n))
    _validate_zero_diagonal(coupling_matrix, name="knm")
    if n == 0:
        return LyapunovState(V=0.0, dV_dt=0.0, in_basin=True, max_phase_diff=0.0)

    diff = phase_values[:, np.newaxis] - phase_values[np.newaxis, :]
    cos_diff = np.cos(diff)

    # Lyapunov fn for Kuramoto gradient system
    # V(θ) = -(1/2N) Σ K_ij cos(θ_i - θ_j)
    # van Hemmen & Wreszinski 1993, Eq. 2.3
    V = -0.5 * float(np.sum(coupling_matrix * cos_diff)) / n

    # Numerical dV/dt from consecutive calls
    dV_dt = 0.0
    if self._prev_V is not None:
        dV_dt = V - self._prev_V
    self._prev_V = V

    # Basin of attraction: all connected pairs within π/2 of each other
    # (sufficient condition for gradient convergence)
    connected = coupling_matrix > 0
    if np.any(connected):
        abs_diff = np.abs(diff)
        # Geodesic distance on S¹: min(|Δ|, 2π-|Δ|)
        abs_diff = np.minimum(abs_diff, 2 * np.pi - abs_diff)
        max_diff = float(np.max(abs_diff[connected]))
    else:
        max_diff = 0.0

    in_basin = max_diff < self._basin_threshold

    return LyapunovState(
        V=V,
        dV_dt=dV_dt,
        in_basin=in_basin,
        max_phase_diff=max_diff,
    )

reset

reset() -> None

Clear cached previous V, so next evaluate() reports dV/dt = 0.

Source code in src/scpn_phase_orchestrator/monitor/lyapunov.py

def reset(self) -> None:
    """Clear cached previous V, so next evaluate() reports dV/dt = 0."""
    self._prev_V = None

Functions:

lyapunov_spectrum

lyapunov_spectrum(
    phases_init: object,
    omegas: object,
    knm: object,
    alpha: object,
    dt: object = 0.01,
    n_steps: object = 1000,
    qr_interval: object = 10,
    zeta: object = 0.0,
    psi: object = 0.0,
) -> FloatArray

Full Lyapunov spectrum (all N exponents) via QR decomposition.

Evolves N perturbation vectors alongside the Kuramoto ODE. Every qr_interval steps, QR-reorthogonalises and accumulates growth rates from the diagonal of R.

Benettin et al. 1980, Meccanica 15:9-20. Shimada & Nagashima 1979, Prog. Theor. Phys. 61:1605-1616.

Dispatches to the first available backend per the SPO fallback chain (Rust → Mojo → Julia → Go → Python). All five produce the same exponents up to floating-point rounding; the dispatcher's choice only affects wall-clock cost.

Parameters

phases_init : object (N,) initial phases. omegas : object (N,) natural frequencies. knm : object (N, N) coupling matrix. alpha : object (N, N) phase-lag matrix. dt : object integration timestep. n_steps : object total integration steps. qr_interval : object steps between QR reorthogonalisations. zeta : object driver strength. psi : object target driver phase.

Returns

FloatArray (N,) array of Lyapunov exponents, sorted descending.

Raises

ValueError If the integration parameters are invalid.

Source code in src/scpn_phase_orchestrator/monitor/lyapunov.py

def lyapunov_spectrum(
    phases_init: object,
    omegas: object,
    knm: object,
    alpha: object,
    dt: object = 0.01,
    n_steps: object = 1000,
    qr_interval: object = 10,
    zeta: object = 0.0,
    psi: object = 0.0,
) -> FloatArray:
    """Full Lyapunov spectrum (all N exponents) via QR decomposition.

    Evolves N perturbation vectors alongside the Kuramoto ODE. Every
    ``qr_interval`` steps, QR-reorthogonalises and accumulates growth
    rates from the diagonal of R.

    Benettin et al. 1980, Meccanica 15:9-20.
    Shimada & Nagashima 1979, Prog. Theor. Phys. 61:1605-1616.

    Dispatches to the first available backend per the SPO fallback
    chain (Rust → Mojo → Julia → Go → Python). All five produce the
    same exponents up to floating-point rounding; the dispatcher's
    choice only affects wall-clock cost.

    Parameters
    ----------
    phases_init : object
        (N,) initial phases.
    omegas : object
        (N,) natural frequencies.
    knm : object
        (N, N) coupling matrix.
    alpha : object
        (N, N) phase-lag matrix.
    dt : object
        integration timestep.
    n_steps : object
        total integration steps.
    qr_interval : object
        steps between QR reorthogonalisations.
    zeta : object
        driver strength.
    psi : object
        target driver phase.

    Returns
    -------
    FloatArray
        (N,) array of Lyapunov exponents, sorted descending.

    Raises
    ------
    ValueError
        If the integration parameters are invalid.
    """
    p = _validate_vector(phases_init, name="phases_init")
    n = int(p.size)
    if n < 1:
        raise ValueError("phases_init must contain at least one oscillator")
    o = _validate_vector(omegas, name="omegas")
    if o.shape != p.shape:
        raise ValueError(f"omegas shape {o.shape} does not match {p.shape}")
    k = _validate_matrix(knm, name="knm", expected_shape=(n, n))
    _validate_zero_diagonal(k, name="knm")
    a = _validate_matrix(alpha, name="alpha", expected_shape=(n, n))
    dt = _validate_positive_real(dt, name="dt")
    n_steps = _validate_int_at_least(n_steps, name="n_steps", minimum=0)
    qr_interval = _validate_int_at_least(
        qr_interval,
        name="qr_interval",
        minimum=1,
    )
    zeta = _validate_non_negative_real(zeta, name="zeta")
    psi = _validate_finite_real(psi, name="psi")
    backend_fn = _dispatch()
    if backend_fn is None:
        return _validate_spectrum_output(
            _lyapunov_spectrum_python(
                p,
                o,
                k,
                a,
                float(dt),
                int(n_steps),
                int(qr_interval),
                float(zeta),
                float(psi),
            ),
            n=n,
        )
    # Rust PyO3 binding takes flat (N*N,) row-major k/alpha; the other
    # backends accept the 2-D forms directly.
    if ACTIVE_BACKEND == "rust":
        try:
            return _validate_spectrum_output(
                backend_fn(
                    p,
                    o,
                    k.ravel(),
                    a.ravel(),
                    dt,
                    n_steps,
                    qr_interval,
                    zeta,
                    psi,
                ),
                n=n,
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            return _validate_spectrum_output(
                _lyapunov_spectrum_python(
                    p,
                    o,
                    k,
                    a,
                    float(dt),
                    int(n_steps),
                    int(qr_interval),
                    float(zeta),
                    float(psi),
                ),
                n=n,
            )
    try:
        return _validate_spectrum_output(
            backend_fn(
                p,
                o,
                k,
                a,
                float(dt),
                int(n_steps),
                int(qr_interval),
                float(zeta),
                float(psi),
            ),
            n=n,
        )
    except (ImportError, RuntimeError, OSError, KeyError):
        return _validate_spectrum_output(
            _lyapunov_spectrum_python(
                p,
                o,
                k,
                a,
                float(dt),
                int(n_steps),
                int(qr_interval),
                float(zeta),
                float(psi),
            ),
            n=n,
        )

Entropy Production

Measures thermodynamic irreversibility of the phase dynamics. Higher entropy production = system further from equilibrium.

entropy_prod

Overdamped-Kuramoto thermodynamic dissipation rate with a 5-backend chain.

Σ = Σ_i (dθ_i/dt)² · dt
dθ_i/dt = ω_i + (α / N) · Σ_j K_ij · sin(θ_j − θ_i)

Zero at frequency-locked fixed points; positive otherwise. Reference: Acebrón et al. 2005, Rev. Mod. Phys. 77:137–185.

Functions:

entropy_production_rate

entropy_production_rate(
    phases: object,
    omegas: object,
    knm: object,
    alpha: object,
    dt: object,
) -> float

Thermodynamic dissipation rate Σ (dθ/dt)² · dt.

dθ_i/dt = ω_i + (α / N) Σ_j K_ij sin(θ_j − θ_i). Zero at frequency-locked fixed points; positive otherwise.

Acebrón et al. 2005, Rev. Mod. Phys. 77:137–185.

Parameters

phases : object (N,) instantaneous phases in radians. omegas : object (N,) natural frequencies. knm : object (N, N) coupling matrix. alpha : object global coupling strength. dt : object integration timestep for the · dt factor.

Returns

float Non-negative dissipation scalar.

Raises

ValueError If the inputs are non-finite or mismatched.

Source code in src/scpn_phase_orchestrator/monitor/entropy_prod.py

def entropy_production_rate(
    phases: object,
    omegas: object,
    knm: object,
    alpha: object,
    dt: object,
) -> float:
    """Thermodynamic dissipation rate ``Σ (dθ/dt)² · dt``.

    ``dθ_i/dt = ω_i + (α / N) Σ_j K_ij sin(θ_j − θ_i)``. Zero at
    frequency-locked fixed points; positive otherwise.

    Acebrón et al. 2005, Rev. Mod. Phys. **77**:137–185.

    Parameters
    ----------
    phases : object
        ``(N,)`` instantaneous phases in radians.
    omegas : object
        ``(N,)`` natural frequencies.
    knm : object
        ``(N, N)`` coupling matrix.
    alpha : object
        global coupling strength.
    dt : object
        integration timestep for the ``· dt`` factor.

    Returns
    -------
    float
        Non-negative dissipation scalar.

    Raises
    ------
    ValueError
        If the inputs are non-finite or mismatched.
    """
    phases = _validate_vector(phases, name="phases")
    n = int(phases.size)
    omegas = _validate_vector(omegas, name="omegas")
    if omegas.shape != phases.shape:
        raise ValueError(f"omegas shape {omegas.shape} does not match {phases.shape}")
    knm = _validate_matrix(knm, name="knm", expected_shape=(n, n))
    alpha = _validate_finite_float(alpha, name="alpha")
    dt = _validate_finite_float(dt, name="dt")
    if dt < 0.0:
        raise ValueError(f"dt must be non-negative, got {dt!r}")
    if n == 0 or dt == 0.0:
        return 0.0
    backend_fn = _dispatch()
    if backend_fn is not None:
        try:
            backend_rate = backend_fn(phases, omegas, knm, alpha, dt)
        except (ImportError, RuntimeError, OSError, KeyError):
            backend_fn = None
        else:
            return _validate_entropy_rate(
                backend_rate,
                name="backend entropy rate",
            )

    diff = phases[np.newaxis, :] - phases[:, np.newaxis]
    coupling = np.sum(knm * np.sin(diff), axis=1)
    dtheta_dt = omegas + (alpha / n) * coupling
    return _validate_entropy_rate(np.sum(dtheta_dt**2) * dt)

Winding Number

Topological charge of phase trajectories. Counts how many times the phase wraps around the circle. Integer-valued topological invariant.

winding

Cumulative winding-number tracker with a 5-backend fallback chain.

w_i = floor(Σ_t wrap(Δθ_{i,t}) / 2π) where wrap(x) ∈ (−π, π]. Counts how many full 2π rotations each oscillator completes across a phase history; positive = counterclockwise, negative = clockwise.

Functions:

winding_numbers

winding_numbers(phases_history: FloatArray) -> IntArray

Cumulative winding number of each oscillator over a trajectory.

w_i = floor(Σ_t wrap(Δθ_{i,t}) / 2π) with wrap(x) ∈ (−π, π].

Parameters

phases_history : FloatArray (T, N) phases in radians.

Returns

IntArray (N,) int64 array of winding numbers.

Source code in src/scpn_phase_orchestrator/monitor/winding.py

def winding_numbers(phases_history: FloatArray) -> IntArray:
    """Cumulative winding number of each oscillator over a trajectory.

    ``w_i = floor(Σ_t wrap(Δθ_{i,t}) / 2π)`` with
    ``wrap(x) ∈ (−π, π]``.

    Parameters
    ----------
    phases_history : FloatArray
        ``(T, N)`` phases in radians.

    Returns
    -------
    IntArray
        ``(N,)`` int64 array of winding numbers.
    """
    phases_history = _validate_phase_history(phases_history)
    if phases_history.ndim != 2 or phases_history.shape[0] < 2:
        n = phases_history.shape[-1] if phases_history.ndim == 2 else 0
        return np.zeros(n, dtype=np.int64)

    t, n = int(phases_history.shape[0]), int(phases_history.shape[1])
    flat: FloatArray = np.ascontiguousarray(phases_history.ravel(), dtype=np.float64)
    expected = _winding_reference(phases_history)

    backend_fn = _dispatch()
    if backend_fn is not None:
        try:
            return _validate_backend_winding(
                backend_fn(flat, t, n),
                n=n,
                t=t,
                expected=expected,
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            n = int(n)

    return _validate_backend_winding(expected, n=n, t=t, expected=expected)

winding_vector

winding_vector(phases_history: FloatArray) -> IntArray

N-dimensional integer classification vector from winding numbers.

Alias for :func:winding_numbers; topologically distinct trajectories map to distinct integer-lattice points.

Parameters

phases_history : FloatArray Phase history, shape (T, N).

Returns

IntArray The integer winding classification vector.

Source code in src/scpn_phase_orchestrator/monitor/winding.py

def winding_vector(phases_history: FloatArray) -> IntArray:
    """N-dimensional integer classification vector from winding numbers.

    Alias for :func:`winding_numbers`; topologically distinct
    trajectories map to distinct integer-lattice points.

    Parameters
    ----------
    phases_history : FloatArray
        Phase history, shape ``(T, N)``.

    Returns
    -------
    IntArray
        The integer winding classification vector.
    """
    return winding_numbers(phases_history)

Inter-Trial Phase Coherence (ITPC)

Phase consistency across repeated trials or time windows. Standard neuroscience measure for event-related phase locking.

itpc

Lachaux 1999 inter-trial phase coherence with a 5-backend fallback chain.

Two kernels:

• :func:compute_itpc — ITPC across trials at each time point.
• :func:itpc_persistence — mean ITPC at stimulus-pause indices.

Functions:

compute_itpc

compute_itpc(phases_trials: object) -> FloatArray

Inter-Trial Phase Coherence at each time point.

ITPC = |mean(exp(i·θ))| across trials (Lachaux et al. 1999).

Parameters

phases_trials : object shape (n_trials, n_timepoints) — phases in radians. A 1-D input is treated as a single trial.

Returns

FloatArray (n_timepoints,) array of ITPC values in [0, 1].

Source code in src/scpn_phase_orchestrator/monitor/itpc.py

def compute_itpc(phases_trials: object) -> FloatArray:
    """Inter-Trial Phase Coherence at each time point.

    ``ITPC = |mean(exp(i·θ))|`` across trials (Lachaux et al. 1999).

    Parameters
    ----------
    phases_trials : object
        shape ``(n_trials, n_timepoints)`` — phases in radians. A 1-D input is treated
        as a single trial.

    Returns
    -------
    FloatArray
        ``(n_timepoints,)`` array of ITPC values in ``[0, 1]``.
    """
    phases = _validate_phases_trials(phases_trials)
    if phases.ndim == 1:
        return np.array([1.0], dtype=np.float64)
    if phases.shape[0] == 0:
        return np.array([], dtype=np.float64)
    n_trials, n_tp = phases.shape
    expected = _compute_itpc_reference(phases)

    backend_fn = _dispatch("itpc")
    if backend_fn is not None:
        try:
            if ACTIVE_BACKEND == "rust":
                fn_rust = cast(
                    "Callable[[FloatArray, int, int], FloatArray]",
                    backend_fn,
                )
                flat = np.ascontiguousarray(phases.ravel(), dtype=np.float64)
                return _validate_itpc_values(
                    fn_rust(flat, n_trials, n_tp),
                    n_timepoints=n_tp,
                    expected=expected,
                )
            fn = cast("Callable[[FloatArray, int, int], FloatArray]", backend_fn)
            return _validate_itpc_values(
                fn(phases.ravel(), int(n_trials), int(n_tp)),
                n_timepoints=n_tp,
                expected=expected,
                atol=1e-9 if ACTIVE_BACKEND == "mojo" else 1e-12,
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            backend_fn = None

    return _validate_itpc_values(expected, n_timepoints=n_tp, expected=expected)

itpc_persistence

itpc_persistence(
    phases_trials: object, pause_indices: object
) -> float

Mean ITPC at stimulus-pause indices.

Distinguishes true neural entrainment from evoked response: if ITPC remains high after the driving stimulus stops, oscillators have genuinely phase-locked. If it drops immediately, the response was merely evoked.

Parameters

phases_trials : object (n_trials, n_timepoints) phases in radians. pause_indices : object time-point indices falling within / after a pause.

Returns

float Mean ITPC across pause_indices. 0.0 if empty.

Source code in src/scpn_phase_orchestrator/monitor/itpc.py

def itpc_persistence(
    phases_trials: object,
    pause_indices: object,
) -> float:
    """Mean ITPC at stimulus-pause indices.

    Distinguishes true neural entrainment from evoked response: if ITPC
    remains high after the driving stimulus stops, oscillators have
    genuinely phase-locked. If it drops immediately, the response was
    merely evoked.

    Parameters
    ----------
    phases_trials : object
        ``(n_trials, n_timepoints)`` phases in radians.
    pause_indices : object
        time-point indices falling within / after a pause.

    Returns
    -------
    float
        Mean ITPC across ``pause_indices``. ``0.0`` if empty.
    """
    phases = _validate_phases_trials(phases_trials)
    pause_idx = _validate_pause_indices(pause_indices)
    if pause_idx.size == 0:
        return 0.0

    if phases.ndim == 1:
        phases = phases.reshape(1, -1)
    n_trials, n_tp = phases.shape
    itpc_full = _compute_itpc_reference(phases)
    valid = pause_idx[(pause_idx >= 0) & (pause_idx < itpc_full.size)]
    expected = 0.0 if valid.size == 0 else float(np.mean(itpc_full[valid]))

    backend_fn = _dispatch("persistence")
    if backend_fn is not None:
        try:
            if ACTIVE_BACKEND == "rust":
                fn_rust = cast(
                    "Callable[[FloatArray, int, int, IntArray], float]",
                    backend_fn,
                )
                return _validate_persistence_value(
                    fn_rust(
                        np.ascontiguousarray(phases.ravel(), dtype=np.float64),
                        n_trials,
                        n_tp,
                        np.ascontiguousarray(pause_idx, dtype=np.int64),
                    ),
                    expected=expected,
                )
            fn = cast("Callable[[FloatArray, int, int, IntArray], float]", backend_fn)
            return _validate_persistence_value(
                fn(phases.ravel(), int(n_trials), int(n_tp), pause_idx),
                expected=expected,
                atol=1e-9 if ACTIVE_BACKEND == "mojo" else 1e-12,
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            backend_fn = None

    return _validate_persistence_value(expected, expected=expected)

Coupling Estimation from Data

Two methods for inferring coupling from observed time series:

1. Basic: least-squares fit of dθ/dt - ω = Σ K_ij sin(θ_j - θ_i)
2. Harmonics: higher Fourier harmonics for non-sinusoidal coupling

The harmonics method captures real biological coupling shapes (Stankovski 2017).

coupling_est

Least-squares coupling estimators for observed phase trajectories.

The primary estimator fits pairwise sinusoidal Kuramoto coupling from phase-history derivatives and natural frequencies, returning a dense matrix with zero diagonal. The harmonics variant expands the regression library with higher Fourier sine and cosine terms. Both routines are offline inference helpers: they estimate parameters from caller-provided arrays and perform no runtime actuation or binding updates.

Functions:

estimate_coupling

estimate_coupling(
    phases: FloatArray, omegas: FloatArray, dt: float
) -> FloatArray

Estimate K_ij coupling matrix from observed phase trajectories.

Least-squares fit of dθ_i/dt - ω_i = Σ_j K_ij sin(θ_j - θ_i). Constructs the regression matrix from pairwise sin(Δθ) and solves for K_ij via pseudoinverse.

Parameters

phases : FloatArray (n_oscillators, n_timesteps) phase trajectories. omegas : FloatArray (n_oscillators,) natural frequencies. dt : float timestep between samples.

Returns

FloatArray (n_oscillators, n_oscillators) estimated coupling matrix K_ij.

Raises

ValueError If the inputs are invalid or inconsistent.

Source code in src/scpn_phase_orchestrator/autotune/coupling_est.py

def estimate_coupling(
    phases: FloatArray,
    omegas: FloatArray,
    dt: float,
) -> FloatArray:
    """Estimate K_ij coupling matrix from observed phase trajectories.

    Least-squares fit of dθ_i/dt - ω_i = Σ_j K_ij sin(θ_j - θ_i).
    Constructs the regression matrix from pairwise sin(Δθ) and solves
    for K_ij via pseudoinverse.

    Parameters
    ----------
    phases : FloatArray
        (n_oscillators, n_timesteps) phase trajectories.
    omegas : FloatArray
        (n_oscillators,) natural frequencies.
    dt : float
        timestep between samples.

    Returns
    -------
    FloatArray
        (n_oscillators, n_oscillators) estimated coupling matrix K_ij.

    Raises
    ------
    ValueError
        If the inputs are invalid or inconsistent.
    """
    phases, omegas, dt = _validate_inputs(phases, omegas, dt)
    n, T = phases.shape
    if T < 3:
        raise ValueError(f"Need >= 3 timesteps, got {T}")

    # Phase derivative: (dθ/dt)_i ≈ (θ_{t+1} - θ_{t-1}) / (2*dt)
    dphase = np.diff(np.unwrap(phases, axis=1), axis=1) / dt
    # Use interior points for derivative
    phases_mid = phases[:, :-1]
    T_eff = dphase.shape[1]

    knm = np.zeros((n, n), dtype=np.float64)

    for i in range(n):
        # Target: dθ_i/dt - ω_i at each timestep
        target = dphase[i, :] - omegas[i]

        # Regressor: sin(θ_j - θ_i) for each j, at each timestep
        regressors = np.sin(
            phases_mid[:, :T_eff] - phases_mid[i : i + 1, :T_eff]
        )  # (n, T_eff)

        # Least squares: target = K_i · regressors
        # K_i = target @ regressors^T @ (regressors @ regressors^T)^{-1}
        with contextlib.suppress(np.linalg.LinAlgError):
            coeffs = np.linalg.lstsq(regressors.T, target, rcond=None)[0]
            if np.all(np.isfinite(coeffs)):
                knm[i, :] = coeffs

    np.fill_diagonal(knm, 0.0)
    return knm

estimate_coupling_harmonics

estimate_coupling_harmonics(
    phases: FloatArray,
    omegas: FloatArray,
    dt: float,
    n_harmonics: int = 2,
) -> dict[str, FloatArray]

Estimate coupling with higher Fourier harmonics.

Fits: dθ_i/dt - ω_i = Σ_j Σ_k [a_jk sin(k·Δθ) + b_jk cos(k·Δθ)] for k = 1..n_harmonics.

Real biological oscillators have non-sinusoidal coupling (Stankovski 2017, Rev. Mod. Phys.).

Returns dict with keys 'sin_1', 'cos_1', 'sin_2', 'cos_2', ... each an (n, n) matrix of coefficients.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). dt : float Integration step size. n_harmonics : int Number of harmonics to fit.

Returns

dict[str, FloatArray] Coupling with higher Fourier harmonics.

Raises

ValueError If the inputs are invalid or inconsistent.

Source code in src/scpn_phase_orchestrator/autotune/coupling_est.py

def estimate_coupling_harmonics(
    phases: FloatArray,
    omegas: FloatArray,
    dt: float,
    n_harmonics: int = 2,
) -> dict[str, FloatArray]:
    """Estimate coupling with higher Fourier harmonics.

    Fits: dθ_i/dt - ω_i = Σ_j Σ_k [a_jk sin(k·Δθ) + b_jk cos(k·Δθ)]
    for k = 1..n_harmonics.

    Real biological oscillators have non-sinusoidal coupling
    (Stankovski 2017, Rev. Mod. Phys.).

    Returns dict with keys 'sin_1', 'cos_1', 'sin_2', 'cos_2', ...
    each an (n, n) matrix of coefficients.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    dt : float
        Integration step size.
    n_harmonics : int
        Number of harmonics to fit.

    Returns
    -------
    dict[str, FloatArray]
        Coupling with higher Fourier harmonics.

    Raises
    ------
    ValueError
        If the inputs are invalid or inconsistent.
    """
    phases, omegas, dt = _validate_inputs(phases, omegas, dt)
    n_harmonics = _validate_n_harmonics(n_harmonics)
    n, T = phases.shape
    if T < 3:
        raise ValueError(f"Need >= 3 timesteps, got {T}")

    dphase = np.diff(np.unwrap(phases, axis=1), axis=1) / dt
    phases_mid = phases[:, :-1]
    T_eff = dphase.shape[1]

    result: dict[str, FloatArray] = {}
    for k in range(1, n_harmonics + 1):
        result[f"sin_{k}"] = np.zeros((n, n), dtype=np.float64)
        result[f"cos_{k}"] = np.zeros((n, n), dtype=np.float64)

    for i in range(n):
        target = dphase[i, :] - omegas[i]
        diff = phases_mid[:, :T_eff] - phases_mid[i : i + 1, :T_eff]

        # Build regressor matrix: [sin(Δθ), cos(Δθ), sin(2Δθ), cos(2Δθ), ...]
        blocks = []
        for k in range(1, n_harmonics + 1):
            blocks.append(np.sin(k * diff))
            blocks.append(np.cos(k * diff))
        regressors = np.vstack(blocks)  # (2*n_harmonics*n, T_eff)

        try:
            coeffs = np.linalg.lstsq(regressors.T, target, rcond=None)[0]
        except np.linalg.LinAlgError:
            continue
        if not np.all(np.isfinite(coeffs)):
            continue

        # Unpack coefficients
        idx = 0
        for k in range(1, n_harmonics + 1):
            result[f"sin_{k}"][i, :] = coeffs[idx : idx + n]
            idx += n
            result[f"cos_{k}"][i, :] = coeffs[idx : idx + n]
            idx += n

    for key in result:
        np.fill_diagonal(result[key], 0.0)

    return result

Read with operational intent

Most monitors are most useful when compared over time and context, not as single point alarms. A practical dashboard should show short-term and rolling-window views side by side, then correlate alarms with known interventions.

The intended use is:

• detect onset conditions with one or two fast indicators,
• confirm with a slower structural monitor,
• only then trigger policy or supervisory changes.

Synthetic HCP Connectome Generation

Generates neuroscience-realistic coupling matrices inspired by the Human Connectome Project:

• Intra-hemispheric exponential distance decay
• Inter-hemispheric corpus callosum pattern
• Default Mode Network hub structure

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()

Monitoring stack as a decision chain

Treat this page as a decision chain for escalation, not a list of separate tools. The intended order is:

1. start with one global stability indicator,
2. confirm structural coherence with pairwise and topology-aware indicators,
3. apply causal or energetic checks before any bounded actuation proposal.

The sequence is designed to reduce false positives and preserve audit quality.

Minimal observability profile

A practical minimum profile for a first production pilot is:

• order_parameter for baseline synchrony,
• lyapunov for local stability trend,
• one causal or directional metric (coupling_est or itpc),
• one action governance metric (evs or pid).

This gives enough signal to decide whether a policy should stay static, reduce its scope, or escalate to broader review.