Advanced Dynamics

Beyond basic Kuramoto, SPO provides several additional dynamical modules.

Variational Free Energy Predictor

Implementation of Friston's Free Energy Principle mapped to Kuramoto dynamics. The system minimises its own prediction error as a coupling objective.

• Precision-weighted prediction error (maps to coupling K_ij)
• KL complexity term on precision
• Online precision estimation from error variance
• Forward prediction with error injection into the UPDE right-hand side

from scpn_phase_orchestrator.upde.prediction import VariationalPredictor

predictor = VariationalPredictor(n=16, dt=0.01)
# prediction_error = predictor.step(phases, predicted_phases)
# error_coupling() returns ε_gain·ε_i for injection into UPDE

prediction

Forward and variational prediction models for validated UPDE phase states.

The module supplies a linear prediction-error model and a variational free-energy predictor over one-dimensional oscillator phase vectors. Public constructors and update methods validate oscillator counts, positive time steps, finite phase/frequency arrays, and precision vectors before mutating internal weights, sufficient statistics, or error histories. The implementation is a concrete numerical mechanism and does not claim to formalize phenomenological time-consciousness.

Classes

PredictionState dataclass

PredictionState(
    predicted_phases: FloatArray,
    prediction_error: FloatArray,
    mean_error: float,
    weights: FloatArray,
)

Snapshot of the forward prediction model after one update step.

PredictionModel

PredictionModel(
    n_oscillators: int,
    learning_rate: float = 0.01,
    error_gain: float = 0.1,
)

Linear forward model for phase prediction.

Predicts θ̂(t+dt) from θ(t) using learned weights W: θ̂(t+dt) = θ(t) + dt · (ω + W · sin(Δθ))

Prediction error ε = θ_actual - θ̂ (wrapped to [-π, π]). Weights updated via gradient descent on ε²: W += η · ε ⊗ sin(Δθ)

The prediction error signal can be injected into the UPDE as an additional coupling term, implementing a form of predictive coding where the system minimizes its own prediction error.

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

def __init__(
    self,
    n_oscillators: int,
    learning_rate: float = 0.01,
    error_gain: float = 0.1,
):
    self._n = _validate_positive_int("n_oscillators", n_oscillators)
    self._lr = _validate_nonnegative_float("learning_rate", learning_rate)
    self._error_gain = _validate_nonnegative_float("error_gain", error_gain)
    self._W = np.zeros((self._n, self._n), dtype=np.float64)
    self._prev_phases: FloatArray | None = None
    self._prev_predicted: FloatArray | None = None

Attributes

weights property

weights: FloatArray

Copy of the current learned weight matrix W.

Returns

FloatArray Copy of the current learned weight matrix W.

error_gain property

error_gain: float

Scaling factor applied to prediction error before injection.

Returns

float Scaling factor applied to prediction error before injection.

Methods:

predict

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

Predict phases at next timestep.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). dt : float Integration step size.

Returns

FloatArray The predicted phases at the next timestep.

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

def predict(self, phases: FloatArray, omegas: FloatArray, dt: float) -> FloatArray:
    """Predict phases at next timestep.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    dt : float
        Integration step size.

    Returns
    -------
    FloatArray
        The predicted phases at the next timestep.
    """
    phases = _validate_vector("phases", phases, self._n)
    omegas = _validate_vector("omegas", omegas, self._n)
    dt = _validate_positive_float("dt", dt)
    diff = phases[np.newaxis, :] - phases[:, np.newaxis]
    coupling_pred = np.sum(self._W * np.sin(diff), axis=1)
    predicted: FloatArray = (phases + dt * (omegas + coupling_pred)) % TWO_PI
    return predicted

update

update(
    phases: FloatArray, omegas: FloatArray, dt: float
) -> PredictionState

Compute prediction error and update weights.

Call once per timestep AFTER the solver step.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). dt : float Integration step size.

Returns

PredictionState The updated prediction state after one learning step.

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

def update(
    self, phases: FloatArray, omegas: FloatArray, dt: float
) -> PredictionState:
    """Compute prediction error and update weights.

    Call once per timestep AFTER the solver step.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    dt : float
        Integration step size.

    Returns
    -------
    PredictionState
        The updated prediction state after one learning step.
    """
    phases = _validate_vector("phases", phases, self._n)
    omegas = _validate_vector("omegas", omegas, self._n)
    dt = _validate_positive_float("dt", dt)
    if self._prev_phases is None or self._prev_predicted is None:
        # First call — no prediction to compare
        predicted = self.predict(phases, omegas, dt)
        self._prev_phases = phases.copy()
        self._prev_predicted = predicted
        return PredictionState(
            predicted_phases=predicted,
            prediction_error=np.zeros(self._n),
            mean_error=0.0,
            weights=self._W.copy(),
        )

    # Prediction error: actual - predicted (wrapped to [-π, π])
    error = phases - self._prev_predicted
    error = (error + np.pi) % TWO_PI - np.pi

    # Weight update: gradient descent on Σ ε_i²
    diff = self._prev_phases[np.newaxis, :] - self._prev_phases[:, np.newaxis]
    sin_diff = np.sin(diff)
    self._W += self._lr * np.outer(error, np.ones(self._n)) * sin_diff

    # Predict next step
    predicted = self.predict(phases, omegas, dt)

    self._prev_phases = phases.copy()
    self._prev_predicted = predicted

    return PredictionState(
        predicted_phases=predicted,
        prediction_error=error,
        mean_error=float(np.mean(np.abs(error))),
        weights=self._W.copy(),
    )

error_coupling

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

Prediction-error signal for injection into UPDE.

Returns ε_gain · ε_i, where ε_i = θ_actual - θ̂_predicted. Add this to the UPDE derivative to implement predictive coding: dθ/dt = ω + K·sin(Δθ) + gain·ε

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). dt : float Integration step size.

Returns

FloatArray The prediction-error coupling signal for UPDE injection.

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

def error_coupling(
    self, phases: FloatArray, omegas: FloatArray, dt: float
) -> FloatArray:
    """Prediction-error signal for injection into UPDE.

    Returns ε_gain · ε_i, where ε_i = θ_actual - θ̂_predicted.
    Add this to the UPDE derivative to implement predictive coding:
      dθ/dt = ω + K·sin(Δθ) + gain·ε

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    dt : float
        Integration step size.

    Returns
    -------
    FloatArray
        The prediction-error coupling signal for UPDE injection.
    """
    phases = _validate_vector("phases", phases, self._n)
    _validate_vector("omegas", omegas, self._n)
    _validate_positive_float("dt", dt)
    if self._prev_predicted is None:
        return np.zeros(self._n)
    error = phases - self._prev_predicted
    error = (error + np.pi) % TWO_PI - np.pi
    out: FloatArray = self._error_gain * error
    return out

reset

reset() -> None

Zero the weight matrix and clear phase history.

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

def reset(self) -> None:
    """Zero the weight matrix and clear phase history."""
    self._W[:] = 0.0
    self._prev_phases = None
    self._prev_predicted = None

VariationalState dataclass

VariationalState(
    predicted_phases: FloatArray,
    error: FloatArray,
    free_energy: float,
    precision: FloatArray,
    complexity: float,
)

Snapshot of the variational predictor after one update step.

VariationalPredictor

VariationalPredictor(
    n_oscillators: int,
    prior_precision: float = 1.0,
    learning_rate: float = 0.01,
)

Variational free energy minimization for phase prediction.

Implements the formal mapping between SCPN phase dynamics and Friston's Free Energy Principle:

F = E_q[log q(theta) - log p(theta, y)] ~ prediction_error^2 / (2 * precision) + complexity

where

theta = phase states (sufficient statistics mu in FEP) y = observed phases q(theta) = recognition density (Gaussian, parameterized by mu, Sigma) prediction_error = y - f(mu) (sensory prediction error) precision = 1/sigma^2 (inverse variance, maps to coupling K) complexity = KL[q||p] (prior deviation cost)

The UPDE coupling term K_ij * sin(theta_j - theta_i) maps to precision-weighted prediction error under Laplace approximation (Friston 2010, Eq. 4).

This is NOT a claim to formalize Husserl's protention. It is a concrete numerical implementation of the mathematical correspondence between Kuramoto coupling and variational inference.

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

def __init__(
    self,
    n_oscillators: int,
    prior_precision: float = 1.0,
    learning_rate: float = 0.01,
):
    self._n = _validate_positive_int("n_oscillators", n_oscillators)
    self._lr = _validate_nonnegative_float("learning_rate", learning_rate)
    self._prior_precision = _validate_positive_float(
        "prior_precision", prior_precision
    )
    # Precision matrix (diagonal): initialized to prior_precision.
    # Under the FEP-Kuramoto correspondence, precision_ij ~ K_ij.
    self._precision = np.full(self._n, self._prior_precision, dtype=np.float64)
    self._mu = np.zeros(self._n, dtype=np.float64)
    self._omegas: FloatArray | None = None
    self._error_history: list[FloatArray] = []
    # Exponential moving average decay for precision updates
    self._ema_alpha = 0.1

Attributes

precision property

precision: FloatArray

Copy of the current per-oscillator precision vector.

Returns

FloatArray Copy of the current per-oscillator precision vector.

Methods:

free_energy

free_energy(
    predicted: FloatArray,
    observed: FloatArray,
    precision: FloatArray,
) -> float

Variational free energy F.

F = sum_i [ (y_i - f(mu_i))^2 * pi_i / 2 ] + sum_i [ log(pi_i) ]

First term: precision-weighted prediction error (accuracy). Second term: log-precision (complexity under Gaussian q). The sign convention follows Friston (2010): F is minimized.

Parameters

predicted : FloatArray Predicted phases in radians, shape (N,). observed : FloatArray Observed phases in radians, shape (N,). precision : FloatArray Per-oscillator precision vector, shape (N,).

Returns

float The variational free energy F.

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

def free_energy(
    self,
    predicted: FloatArray,
    observed: FloatArray,
    precision: FloatArray,
) -> float:
    """Variational free energy F.

    F = sum_i [ (y_i - f(mu_i))^2 * pi_i / 2 ] + sum_i [ log(pi_i) ]

    First term: precision-weighted prediction error (accuracy).
    Second term: log-precision (complexity under Gaussian q).
    The sign convention follows Friston (2010): F is minimized.

    Parameters
    ----------
    predicted : FloatArray
        Predicted phases in radians, shape ``(N,)``.
    observed : FloatArray
        Observed phases in radians, shape ``(N,)``.
    precision : FloatArray
        Per-oscillator precision vector, shape ``(N,)``.

    Returns
    -------
    float
        The variational free energy ``F``.
    """
    predicted = _validate_vector("predicted", predicted, self._n)
    observed = _validate_vector("observed", observed, self._n)
    precision = _validate_positive_vector("precision", precision, self._n)
    error = observed - predicted
    error = (error + np.pi) % TWO_PI - np.pi
    accuracy = float(np.sum(error**2 * precision / 2.0))
    # KL complexity: log-precision acts as a regularizer pulling
    # precision toward values where q(theta) stays close to prior p(theta).
    complexity = float(np.sum(np.log(np.maximum(precision, 1e-12))))
    return accuracy + complexity

update

update(
    phases: FloatArray, omegas: FloatArray, dt: float
) -> VariationalState

One variational update step.

1. Predict phases from current sufficient statistics mu.
2. Compute precision-weighted prediction error.
3. Update mu (gradient descent on F).
4. Update precision from error statistics.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). dt : float Integration step size.

Returns

VariationalState The updated variational state after one step.

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

def update(
    self, phases: FloatArray, omegas: FloatArray, dt: float
) -> VariationalState:
    """One variational update step.

    1. Predict phases from current sufficient statistics mu.
    2. Compute precision-weighted prediction error.
    3. Update mu (gradient descent on F).
    4. Update precision from error statistics.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    dt : float
        Integration step size.

    Returns
    -------
    VariationalState
        The updated variational state after one step.
    """
    phases = _validate_vector("phases", phases, self._n)
    omegas = _validate_vector("omegas", omegas, self._n)
    dt = _validate_positive_float("dt", dt)
    self._omegas = omegas

    # Forward model: f(mu) = mu + dt * omega (simplest generative model)
    predicted = (self._mu + dt * omegas) % TWO_PI

    # Prediction error wrapped to [-pi, pi]
    error = phases - predicted
    error = (error + np.pi) % TWO_PI - np.pi

    # Free energy before update
    fe = self.free_energy(predicted, phases, self._precision)

    # Complexity: KL divergence between current and prior precision.
    # For diagonal Gaussian: KL = 0.5 * sum(pi/pi_0 - 1 - log(pi/pi_0))
    ratio = self._precision / self._prior_precision
    complexity = float(0.5 * np.sum(ratio - 1.0 - np.log(np.maximum(ratio, 1e-12))))

    # Gradient descent on F w.r.t. mu:
    # dF/dmu = -precision * error  (since F ~ precision * error^2 / 2)
    # mu_new = mu - lr * dF/dmu = mu + lr * precision * error
    # type ignore: modulo preserves ndarray shape, but mypy narrows to scalar.
    self._mu = (self._mu + self._lr * self._precision * error) % TWO_PI  # type: ignore[assignment]

    # Update precision from error variance (online).
    # Precision = 1/variance. Use EMA of squared error as variance estimate.
    self._error_history.append(error.copy())
    if len(self._error_history) > 1:
        recent = np.array(self._error_history[-min(50, len(self._error_history)) :])
        var_estimate = np.mean(recent**2, axis=0)
        # EMA blend: pi_new = (1-a)*pi_old + a*(1/var)
        new_prec = 1.0 / np.maximum(var_estimate, 1e-8)
        self._precision = (
            1.0 - self._ema_alpha
        ) * self._precision + self._ema_alpha * new_prec

    return VariationalState(
        predicted_phases=predicted,
        error=error,
        free_energy=fe,
        precision=self._precision.copy(),
        complexity=complexity,
    )

precision_weighted_coupling

precision_weighted_coupling() -> FloatArray

Precision matrix interpretable as K_ij.

Under the FEP-Kuramoto correspondence (Friston 2010, Laplace approximation), the coupling matrix K_ij maps to the off-diagonal elements of the precision matrix of the generative model.

This returns diag(precision) as the simplest such mapping. For a full N x N coupling matrix, use np.diag(result).

Returns

FloatArray Precision matrix interpretable as K_ij.

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

def precision_weighted_coupling(self) -> FloatArray:
    """Precision matrix interpretable as K_ij.

    Under the FEP-Kuramoto correspondence (Friston 2010, Laplace
    approximation), the coupling matrix K_ij maps to the off-diagonal
    elements of the precision matrix of the generative model.

    This returns diag(precision) as the simplest such mapping.
    For a full N x N coupling matrix, use np.diag(result).

    Returns
    -------
    FloatArray
        Precision matrix interpretable as K_ij.
    """
    return np.diag(self._precision)

reset

reset() -> None

Reset precision to prior, zero sufficient statistics, clear history.

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

def reset(self) -> None:
    """Reset precision to prior, zero sufficient statistics, clear history."""
    self._precision[:] = self._prior_precision
    self._mu[:] = 0.0
    self._omegas = None
    self._error_history.clear()

Combinatorial Hodge Decomposition of Coupling

Decomposes the Kuramoto coupling current f_ij = ½(K_ij + K_ji)·sin(θ_j − θ_i) into three L²-orthogonal edge flows via combinatorial Hodge theory on the simplicial complex of oscillators (Jiang, Lim, Yao & Ye 2011):

• Gradient: conservative, curl-free flow grad(s) from a node potential.
• Curl: divergence-free rotational flow bounded by triangles.
• Harmonic: topological residual in the kernel of the Hodge 1-Laplacian; its dimension is the first Betti number β₁ — non-zero only when the graph carries cycles that no triangle fills.

Answers: "Is this synchronisation conservative, rotational, or carried by an irreducible topological cycle?"

from scpn_phase_orchestrator.coupling.hodge import hodge_decomposition

result = hodge_decomposition(K, phases)
result.gradient   # (N, N) antisymmetric conservative flow
result.curl       # (N, N) rotational flow
result.harmonic   # (N, N) topological flow
result.betti_one  # number of independent unfilled cycles

See the Coupling API reference for the full hodge_decomposition and HodgeResult signatures.

Simplicial (3-Body) Coupling

Higher-order interactions beyond pairwise: the 3-body term induces explosive (first-order) synchronization transitions.

dθ_i/dt = ω_i + (σ₁/N) Σ_j K_ij sin(θ_j - θ_i)
                + (σ₂/N²) Σ_{j,k} sin(θ_j + θ_k - 2θ_i)

Gambuzza et al. 2023, Nature Physics; Tang et al. 2025.

simplicial

Pairwise + all-to-all 3-body (simplicial) Kuramoto with a 5-backend chain.

Model

dθ_i/dt = ω_i
          + (σ₁/N) · Σ_j A_ij · sin(θ_j − θ_i)
          + (σ₂/N²) · Σ_{j,k} sin(θ_j + θ_k − 2θ_i)
          + ζ · sin(ψ − θ_i)

σ₂ > 0 drives explosive (first-order) transitions and shrinks basins of attraction while improving the locking stability of already-synchronous states (Gambuzza et al. 2023; Tang et al. 2025).

Closed form for the 3-body sum

Expanding sin(θ_j + θ_k − 2θ_i) = sin((θ_j − θ_i) + (θ_k − θ_i)) and separating the cross terms gives

Σ_{j,k} sin(θ_j + θ_k − 2θ_i) = 2 · S_i · C_i

with

S_i = Σ_j sin(θ_j − θ_i) = (Σ sin θ)·cos θ_i − (Σ cos θ)·sin θ_i
C_i = Σ_j cos(θ_j − θ_i) = (Σ cos θ)·cos θ_i + (Σ sin θ)·sin θ_i

So the 3-body contribution is evaluated in O(N²) (not O(N³)) using two global sums plus the per-node sincos expansion. All five backends use this identity; the pairwise path matches the Rust kernel's sincos expansion on the alpha-zero branch and the direct sin(diff) form otherwise, giving bit-exact parity.

Classes

SimplicialEngine

SimplicialEngine(
    n_oscillators: int, dt: float, sigma2: float = 0.0
)

Pairwise + simplicial (3-body, all-to-all) Kuramoto stepper.

The engine's geometry is (n, dt, σ₂); the step itself is stateless: (phases, omegas, K, α, ζ, ψ) → new_phases.

Initialise the simplicial Kuramoto stepper.

Parameters

n_oscillators : int Number of oscillators in the fixed engine geometry. dt : float Positive Euler timestep in seconds. sigma2 : float, default=0.0 Non-negative all-to-all triadic coupling strength.

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

def __init__(self, n_oscillators: int, dt: float, sigma2: float = 0.0):
    """Initialise the simplicial Kuramoto stepper.

    Parameters
    ----------
    n_oscillators : int
        Number of oscillators in the fixed engine geometry.
    dt : float
        Positive Euler timestep in seconds.
    sigma2 : float, default=0.0
        Non-negative all-to-all triadic coupling strength.
    """
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._dt = _validate_positive_float(dt, name="dt")
    self._sigma2 = _validate_nonnegative_float(sigma2, name="sigma2")

Attributes

sigma2 property writable

sigma2: float

Return the configured all-to-all triadic coupling strength.

Returns

float Return the configured all-to-all triadic coupling strength.

Methods:

step

step(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
) -> FloatArray

Advance one pairwise-plus-simplicial Kuramoto timestep.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray Phase-lag matrix in radians, shape (N, N), or None for no lag.

Returns

FloatArray The phases after one pairwise-plus-simplicial step.

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
) -> FloatArray:
    """Advance one pairwise-plus-simplicial Kuramoto timestep.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.

    Returns
    -------
    FloatArray
        The phases after one pairwise-plus-simplicial step.
    """
    return self.run(phases, omegas, knm, zeta, psi, alpha, n_steps=1)

run

run(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    n_steps: int,
) -> FloatArray

Integrate pairwise-plus-simplicial Kuramoto dynamics.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray Phase-lag matrix in radians, shape (N, N), or None for no lag. n_steps : int Number of integration steps to run.

Returns

FloatArray The final phases after n_steps simplicial steps.

Raises

ValueError If n_steps is negative or the state arrays are invalid.

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    n_steps: int,
) -> FloatArray:
    """Integrate pairwise-plus-simplicial Kuramoto dynamics.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    FloatArray
        The final phases after ``n_steps`` simplicial steps.

    Raises
    ------
    ValueError
        If ``n_steps`` is negative or the state arrays are invalid.
    """
    n_steps = _validate_nonnegative_int(n_steps, name="n_steps")
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    omegas64 = _validate_state_array(omegas, name="omegas", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    if np.any(np.diag(knm64) != 0.0):
        raise ValueError("knm diagonal must be exactly zero")
    alpha64 = _validate_state_array(alpha, name="alpha", shape=(self._n, self._n))
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    if n_steps == 0:
        return np.asarray(phases64, dtype=np.float64).copy()
    knm_flat = knm64.ravel()
    alpha_flat = alpha64.ravel()
    backend_fn = _dispatch()
    if backend_fn is not None:
        try:
            backend_out = backend_fn(
                phases64,
                omegas64,
                knm_flat,
                alpha_flat,
                self._n,
                zeta,
                psi,
                float(self._sigma2),
                float(self._dt),
                int(n_steps),
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            return _validate_backend_result(
                _python_run(
                    phases64,
                    omegas64,
                    knm_flat,
                    alpha_flat,
                    self._n,
                    zeta,
                    psi,
                    float(self._sigma2),
                    float(self._dt),
                    int(n_steps),
                ),
                name="backend output",
                n=self._n,
            )
        return _validate_backend_result(
            backend_out,
            name="backend output",
            n=self._n,
        )
    return _validate_backend_result(
        _python_run(
            phases64,
            omegas64,
            knm_flat,
            alpha_flat,
            self._n,
            zeta,
            psi,
            float(self._sigma2),
            float(self._dt),
            int(n_steps),
        ),
        name="backend output",
        n=self._n,
    )

order_parameter

order_parameter(phases: FloatArray) -> float

Compute the standard Kuramoto R = ||.

Parameters

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

Returns

float The Kuramoto order parameter R.

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

def order_parameter(self, phases: FloatArray) -> float:
    """Compute the standard Kuramoto R = |<exp(iθ)>|.

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

    Returns
    -------
    float
        The Kuramoto order parameter ``R``.
    """
    phases64 = _validate_state_array(
        phases,
        name="phases",
        shape=(self._n,),
    )
    return float(np.abs(np.mean(np.exp(1j * phases64))))

Functions:

Time-Delayed Coupling

Circular buffer supports arbitrary time delays with automatic fallback to instantaneous coupling. Time delays generate "effective higher-order interactions for free" (Ciszak et al. 2025).

delay

Time-delayed Kuramoto buffer and engine with validated phase history.

DelayBuffer stores copied finite phase snapshots in a bounded deque, and DelayedEngine advances phases with delayed coupling, optional external forcing, and Rust acceleration when available. Constructors and step inputs reject non-positive dimensions, non-finite scalars, shape-mismatched arrays, and boolean or numeric-string aliases before integration so delayed history never aliases invalid caller state.

Classes

DelayBuffer

DelayBuffer(n_oscillators: int, max_delay_steps: int)

Circular buffer storing phase history for delayed coupling.

Stores last max_delay_steps snapshots. Retrieves phases from delay_steps steps ago.

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

def __init__(self, n_oscillators: int, max_delay_steps: int):
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._max = _validate_positive_int(max_delay_steps, name="max_delay_steps")
    self._buffer: deque[FloatArray] = deque(maxlen=self._max)

Attributes

length property

length: int

Number of snapshots currently stored.

Returns

int Number of snapshots currently stored.

Methods:

push

push(phases: FloatArray) -> None

Append a phase snapshot to the buffer.

Parameters

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

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

def push(self, phases: FloatArray) -> None:
    """Append a phase snapshot to the buffer.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    """
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    self._buffer.append(phases64.copy())

get_delayed

get_delayed(delay_steps: int) -> FloatArray | None

Return phases from delay_steps ago, or None if not enough history.

Parameters

delay_steps : int Number of steps in the past to retrieve from the delay buffer.

Returns

FloatArray | None The phase snapshot from delay_steps ago, or None if history is short.

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

def get_delayed(self, delay_steps: int) -> FloatArray | None:
    """Return phases from `delay_steps` ago, or None if not enough history.

    Parameters
    ----------
    delay_steps : int
        Number of steps in the past to retrieve from the delay buffer.

    Returns
    -------
    FloatArray | None
        The phase snapshot from ``delay_steps`` ago, or ``None`` if history is
        short.
    """
    delay = _validate_positive_int(delay_steps, name="delay_steps")
    if delay > len(self._buffer):
        return None
    return self._buffer[-delay]

clear

clear() -> None

Discard all stored phase snapshots.

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

def clear(self) -> None:
    """Discard all stored phase snapshots."""
    self._buffer.clear()

DelayedEngine

DelayedEngine(
    n_oscillators: int, dt: float, delay_steps: int = 1
)

Kuramoto with time-delayed coupling.

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

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

def __init__(self, n_oscillators: int, dt: float, delay_steps: int = 1):
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._dt = _validate_positive_float(dt, name="dt")
    self._delay_steps = _validate_positive_int(delay_steps, name="delay_steps")
    self._buffer: deque[FloatArray] = deque(maxlen=self._delay_steps + 1)

Attributes

delay_steps property

delay_steps: int

Return the configured discrete coupling delay.

Returns

int Return the configured discrete coupling delay.

Methods:

step

step(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    step_idx: int = 0,
) -> FloatArray

Advance one delayed Kuramoto timestep from validated state arrays.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray | None Phase-lag matrix in radians, shape (N, N), or None for no lag. step_idx : int Zero-based index of the current step, used to address delayed coupling history.

Returns

FloatArray The phases after one delayed Kuramoto step, in [0, 2π).

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    step_idx: int = 0,
) -> FloatArray:
    """Advance one delayed Kuramoto timestep from validated state arrays.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    step_idx : int
        Zero-based index of the current step, used to address delayed coupling
        history.

    Returns
    -------
    FloatArray
        The phases after one delayed Kuramoto step, in ``[0, 2π)``.
    """
    del step_idx
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    omegas64 = _validate_state_array(omegas, name="omegas", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    alpha64: FloatArray
    if alpha is None:
        alpha64 = np.zeros((self._n, self._n), dtype=np.float64)
    else:
        alpha64 = _validate_state_array(
            alpha,
            name="alpha",
            shape=(self._n, self._n),
        )
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    self._buffer.append(phases64.copy())
    delayed = self._buffer[0] if len(self._buffer) > self._delay_steps else phases64
    diff = delayed[np.newaxis, :] - phases64[:, np.newaxis] - alpha64
    coupling = np.sum(knm64 * np.sin(diff), axis=1)
    dtheta = omegas64 + coupling
    if zeta != 0.0:
        dtheta += zeta * np.sin(psi - phases64)
    step_out: FloatArray = (phases64 + self._dt * dtheta) % TWO_PI
    return step_out

run

run(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    n_steps: int = 100,
) -> FloatArray

Run delayed Kuramoto integration for n_steps validated steps.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray | None Phase-lag matrix in radians, shape (N, N), or None for no lag. n_steps : int Number of integration steps to run.

Returns

FloatArray The final phases after n_steps delayed Kuramoto steps.

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    n_steps: int = 100,
) -> FloatArray:
    """Run delayed Kuramoto integration for ``n_steps`` validated steps.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    FloatArray
        The final phases after ``n_steps`` delayed Kuramoto steps.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    omegas64 = _validate_state_array(omegas, name="omegas", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    alpha64: FloatArray
    if alpha is None:
        alpha64 = np.zeros((self._n, self._n), dtype=np.float64)
    else:
        alpha64 = _validate_state_array(
            alpha,
            name="alpha",
            shape=(self._n, self._n),
        )
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    knm_flat = np.ascontiguousarray(knm64.ravel(), dtype=np.float64)
    alpha_flat = np.ascontiguousarray(alpha64.ravel(), dtype=np.float64)
    backend_fn = _dispatch()
    if backend_fn is not None:
        try:
            return _validate_phase_output(
                backend_fn(
                    phases64,
                    omegas64,
                    knm_flat,
                    alpha_flat,
                    self._n,
                    zeta,
                    psi,
                    self._dt,
                    self._delay_steps,
                    n_steps,
                ),
                n_oscillators=self._n,
            )
        except (ImportError, RuntimeError, OSError, KeyError, ValueError):
            pass
    return _python_run(
        phases64,
        omegas64,
        knm_flat,
        alpha_flat,
        self._n,
        zeta,
        psi,
        self._dt,
        self._delay_steps,
        n_steps,
    )

Functions:

Stochastic Resonance with Optimal Noise

Euler-Maruyama integration with automatic optimal noise tuning. Counter-intuitive: noise at D* INCREASES synchronization.

Optimal noise: D* ≈ K·R_det/2 (Tselios et al. 2025). Self-consistency via modified Bessel transcendental equation (Acebrón et al. 2005).

stochastic

Stochastic noise injection and noise-level sweeps for UPDE phase dynamics.

StochasticInjector owns a local random generator and applies Euler-Maruyama phase noise under validated non-negative diffusion and positive time-step parameters. find_optimal_noise sweeps finite non-negative candidate noise levels against a supplied UPDE engine and reports the best coherence profile without changing the engine configuration or caller-provided input arrays outside normal engine stepping.

Classes

NoiseProfile dataclass

NoiseProfile(
    D: float, R_achieved: float, R_deterministic: float
)

Validated noise-sweep result linking diffusion to bounded order.

StochasticInjector

StochasticInjector(D: float, seed: int | None = None)

Add calibrated noise to phase dynamics.

Euler-Maruyama: θ_i(t+dt) = θ_i(t) + f(θ)dt + √(2Ddt) * ξ_i where ξ_i ~ N(0,1) i.i.d.

Tselios et al. 2025 — stochastic resonance in Kuramoto networks.

Create an injector with finite D and an optional valid seed.

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

def __init__(self, D: float, seed: int | None = None):
    """Create an injector with finite ``D`` and an optional valid seed."""
    self._D = _validate_finite_non_negative(D, name="D")
    self._rng = np.random.default_rng(_validate_optional_seed(seed))

Attributes

D property writable

D: float

Return the configured non-negative diffusion coefficient.

Returns

float Return the configured non-negative diffusion coefficient.

Methods:

inject

inject(phases: FloatArray, dt: float) -> FloatArray

Add Wiener noise to phases: θ += √(2D*dt) * N(0,1).

Parameters

phases : FloatArray Finite real numeric oscillator phases in radians, shape (N,). Boolean, complex, and numeric-string aliases are rejected. dt : float Integration step size.

Returns

FloatArray The phases with added Wiener noise.

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

def inject(self, phases: FloatArray, dt: float) -> FloatArray:
    """Add Wiener noise to phases: θ += √(2D*dt) * N(0,1).

    Parameters
    ----------
    phases : FloatArray
        Finite real numeric oscillator phases in radians, shape ``(N,)``.
        Boolean, complex, and numeric-string aliases are rejected.
    dt : float
        Integration step size.

    Returns
    -------
    FloatArray
        The phases with added Wiener noise.
    """
    dt = _validate_finite_positive(dt, name="dt")
    phases = _validate_phases(phases)
    if self._D == 0.0:
        return phases
    noise = self._rng.standard_normal(len(phases))
    result: FloatArray = (phases + np.sqrt(2.0 * self._D * dt) * noise) % TWO_PI
    return result

Functions:

optimal_D

optimal_D(K: float, R_det: float) -> float

Estimate optimal noise for stochastic resonance.

D* ≈ K·R_det/2 (common noise case). Tselios et al. 2025.

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

def optimal_D(K: float, R_det: float) -> float:
    """Estimate optimal noise for stochastic resonance.

    D* ≈ K·R_det/2 (common noise case).
    Tselios et al. 2025.
    """
    return K * R_det / 2.0

find_optimal_noise

find_optimal_noise(
    engine: UPDEEngine,
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray,
    D_range: FloatArray | None = None,
    n_steps: int = 500,
    seed: int = 42,
) -> NoiseProfile

Sweep noise levels, return D that maximizes R.

Uses the engine to simulate n_steps at each D value.

Parameters

engine : UPDEEngine The UPDE engine used to integrate each trial. phases_init : FloatArray Initial oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). alpha : FloatArray Phase-lag matrix in radians, shape (N, N), or None for no lag. D_range : FloatArray | None Finite non-negative real numeric diffusion coefficients to sweep, or None for the default range. Coercive aliases are rejected. n_steps : int Number of integration steps to run. seed : int Non-negative non-boolean seed for the deterministic RNG.

Returns

NoiseProfile The noise profile whose diffusion D maximises R.

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

def find_optimal_noise(
    engine: UPDEEngine,
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray,
    D_range: FloatArray | None = None,
    n_steps: int = 500,
    seed: int = 42,
) -> NoiseProfile:
    """Sweep noise levels, return D that maximizes R.

    Uses the engine to simulate n_steps at each D value.

    Parameters
    ----------
    engine : UPDEEngine
        The UPDE engine used to integrate each trial.
    phases_init : FloatArray
        Initial oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    alpha : FloatArray
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    D_range : FloatArray | None
        Finite non-negative real numeric diffusion coefficients to sweep, or
        ``None`` for the default range. Coercive aliases are rejected.
    n_steps : int
        Number of integration steps to run.
    seed : int
        Non-negative non-boolean seed for the deterministic RNG.

    Returns
    -------
    NoiseProfile
        The noise profile whose diffusion ``D`` maximises ``R``.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    seed = _validate_seed(seed)
    D_range = _validate_noise_range(D_range)
    if D_range is None:
        K_mean = float(np.mean(knm[knm > 0])) if np.any(knm > 0) else 1.0
        D_range = np.linspace(0.0, K_mean, 11, dtype=np.float64)

    best_D = 0.0
    best_R = 0.0
    R_det = 0.0

    for i, D in enumerate(D_range):
        phases = phases_init.copy()
        injector = StochasticInjector(D, seed=seed + i)
        for _ in range(n_steps):
            phases = engine.step(phases, omegas, knm, 0.0, 0.0, alpha)
            if D > 0:
                phases = injector.inject(phases, engine._dt)
        R, _ = compute_order_parameter(phases)
        if i == 0:
            R_det = R
        if best_R < R:
            best_R = R
            best_D = float(D)

    return NoiseProfile(D=best_D, R_achieved=best_R, R_deterministic=R_det)

Geometric (Torus-Preserving) Integrator

Symplectic Euler on T^N via SO(2) exponential map. Works in unit complex representation z_i = exp(iθ_i), avoiding mod 2π discontinuity errors that cause subtle numerical drift in standard integrators.

Essential for long timescale simulations where standard Euler accumulates phase wrapping errors.

geometric

Torus-preserving symplectic Euler integrator on T^N = (S¹)^N.

Exposes a 5-backend fallback chain.

Scheme

Each phase is lifted to the unit circle z_i = exp(iθ_i); the Kuramoto derivative ω_eff_i is computed in the tangent space, and z_i is advanced by the exponential map

z_i(t + dt) = z_i(t) · exp(i · ω_eff_i · dt)

followed by renormalisation to the unit circle. This avoids the mod-2π discontinuity that introduces subtle truncation errors in standard integrators when trajectories cross θ = 0.

Across the five backends the (z_re, z_im) state is carried in between steps (no atan2 round-trip per step), matching the Rust kernel spo-engine/src/geometric.rs bit-for-bit. The pairwise derivative uses the sincos expansion on the alpha == 0 branch and the direct atan2 + sin(diff) form otherwise.

Classes

TorusEngine

TorusEngine(n_oscillators: int, dt: float)

Symplectic Euler on T^N with 5-backend dispatch.

Store (n, dt); step / run are stateless in (θ, ω, K, α, ζ, ψ).

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

def __init__(self, n_oscillators: int, dt: float):
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._dt = _validate_positive_float(dt, name="dt")

Methods:

step

step(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
) -> FloatArray

One torus step; returns phases in [0, 2π).

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray Phase-lag matrix in radians, shape (N, N), or None for no lag.

Returns

FloatArray The phases after one torus step, in [0, 2π).

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
) -> FloatArray:
    """One torus step; returns phases in ``[0, 2π)``.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.

    Returns
    -------
    FloatArray
        The phases after one torus step, in ``[0, 2π)``.
    """
    return self.run(phases, omegas, knm, zeta, psi, alpha, n_steps=1)

run

run(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    n_steps: int,
) -> FloatArray

Integrate torus phase dynamics for the requested number of steps.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray Phase-lag matrix in radians, shape (N, N), or None for no lag. n_steps : int Number of integration steps to run.

Returns

FloatArray The final finite torus phases after n_steps torus steps, in [0, 2π).

Raises

ValueError If the submitted state is malformed or an optional backend returns a phase vector outside the public torus contract.

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    n_steps: int,
) -> FloatArray:
    """Integrate torus phase dynamics for the requested number of steps.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    FloatArray
        The final finite torus phases after ``n_steps`` torus steps, in
        ``[0, 2π)``.

    Raises
    ------
    ValueError
        If the submitted state is malformed or an optional backend returns
        a phase vector outside the public torus contract.
    """
    n_steps = _validate_nonnegative_int(n_steps, name="n_steps")
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    omegas64 = _validate_state_array(omegas, name="omegas", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    alpha64 = _validate_state_array(
        alpha,
        name="alpha",
        shape=(self._n, self._n),
    )
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    knm_flat = knm64.ravel()
    alpha_flat = alpha64.ravel()
    backend_fn = _dispatch()
    if backend_fn is not None:
        return validate_torus_output(
            backend_fn(
                phases64,
                omegas64,
                knm_flat,
                alpha_flat,
                self._n,
                float(zeta),
                float(psi),
                float(self._dt),
                n_steps,
            ),
            n=self._n,
        )
    return _python_torus_run(
        phases64,
        omegas64,
        knm_flat,
        alpha_flat,
        self._n,
        float(zeta),
        float(psi),
        float(self._dt),
        n_steps,
    )

order_parameter

order_parameter(phases: FloatArray) -> float

Compute the standard Kuramoto R = ||.

Parameters

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

Returns

float The Kuramoto order parameter R.

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

def order_parameter(self, phases: FloatArray) -> float:
    """Compute the standard Kuramoto R = |<exp(iθ)>|.

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

    Returns
    -------
    float
        The Kuramoto order parameter ``R``.
    """
    phases64 = _validate_state_array(
        phases,
        name="phases",
        shape=(self._n,),
    )
    return float(np.abs(np.mean(np.exp(1j * phases64))))

Functions:

Ott-Antonsen Mean-Field Reduction

Exact analytical reduction for Lorentzian frequency distributions:

dz/dt = -(Δ + iω₀)z + (K/2)(z - |z|²z)

Critical coupling K_c = 2Δ. Steady-state: R_ss = √(1 - 2Δ/K). Used by the PredictiveSupervisor as a fast forward model for MPC (O(1) vs O(N) for full simulation).

reduction

Exact mean-field (Ott-Antonsen) reduction for globally-coupled Kuramoto.

Uses a Lorentzian g(ω) and a 5-backend fallback chain.

Dynamics

On the Ott-Antonsen manifold the full N-oscillator Kuramoto system reduces to a single complex-scalar ODE:

dz/dt = −(Δ + iω₀)·z + (K/2)·(z − |z|²·z)

with z = R·e^{iψ} the mean-field order parameter, Δ the half-width of the Lorentzian g(ω), ω₀ its centre, and K the coupling strength.

Steady-state R_ss = √(1 − 2Δ/K) for K > K_c = 2Δ and R_ss = 0 below. Reference: Ott & Antonsen 2008, Chaos 18(3):037113.

Numerics

run(z0, n_steps) is the compute-kernel path: a tight RK4 loop on the real/imaginary components of z. This is dispatched across Rust / Mojo / Julia / Go / Python with bit-exact parity (scalar ODE, no reduction identities, no global sums — the only differences between backends are the rounding order of the k1..k4 accumulation, which matches exactly).

The scalar-output helpers K_c, steady_state_R and predict_from_oscillators stay native Python + optional Rust — they are O(1) arithmetic or O(N) percentile work and do not benefit from multi-language chains.

Classes

OAState dataclass

OAState(z: complex, R: float, psi: float, K_c: float)

Ott-Antonsen mean-field state: order parameter and critical coupling.

OttAntonsenReduction

OttAntonsenReduction(
    omega_0: float, delta: float, K: float, dt: float = 0.01
)

Ott-Antonsen mean-field reduction for globally-coupled Kuramoto.

The class stores (ω₀, Δ, K, dt) and exposes K_c, steady_state_R(), step(z), run(z0, n_steps) and predict_from_oscillators(omegas, K). run is dispatched across the 5-backend chain; the scalar helpers stay Python + optional Rust.

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

def __init__(
    self,
    omega_0: float,
    delta: float,
    K: float,
    dt: float = 0.01,
):
    omega_0 = _validate_finite_real(omega_0, name="omega_0")
    delta = _validate_finite_real(delta, name="delta")
    K = _validate_finite_real(K, name="K")
    dt = _validate_finite_real(dt, name="dt")
    if delta < 0:
        raise ValueError(f"delta (half-width) must be non-negative, got {delta}")
    if dt <= 0.0:
        raise ValueError(f"dt must be positive, got {dt}")
    self._omega_0 = omega_0
    self._delta = delta
    self._K = K
    self._dt = dt

Attributes

K_c property

K_c: float

Critical coupling K_c = 2Δ.

Returns

float Critical coupling K_c = 2Δ.

Methods:

steady_state_R

steady_state_R() -> float

Return the analytical steady-state R_ss = √(1 − 2Δ/K) for K > K_c.

Returns

float Return the analytical steady-state R_ss = √(1 − 2Δ/K) for K > K_c.

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

def steady_state_R(self) -> float:
    """Return the analytical steady-state ``R_ss = √(1 − 2Δ/K)`` for ``K > K_c``.

    Returns
    -------
    float
        Return the analytical steady-state ``R_ss = √(1 − 2Δ/K)`` for ``K > K_c``.
    """
    if _HAS_RUST_SCALAR:
        value = _rust_steady_state_r(self._delta, self._K)
    elif self.K_c >= self._K:
        value = 0.0
    else:
        value = (1.0 - 2.0 * self._delta / self._K) ** 0.5
    return _reduction_validation.validate_oa_steady_state_output(value)

step

step(z: complex) -> complex

Single RK4 step on the OA ODE.

Parameters

z : complex Complex Ott-Antonsen order parameter.

Returns

complex The complex order parameter after one RK4 step.

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

def step(self, z: complex) -> complex:
    """Single RK4 step on the OA ODE.

    Parameters
    ----------
    z : complex
        Complex Ott-Antonsen order parameter.

    Returns
    -------
    complex
        The complex order parameter after one RK4 step.
    """
    z = _validate_finite_complex(z, name="z")
    re, im, _, _ = self._run_scalar(z.real, z.imag, n_steps=1)
    return complex(re, im)

run

run(z0: complex, n_steps: int) -> OAState

Integrate n_steps RK4 steps; return the final OAState.

Parameters

z0 : complex Initial complex Ott-Antonsen order parameter. n_steps : int Number of integration steps to run.

Returns

OAState The final OAState after n_steps RK4 steps.

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

def run(self, z0: complex, n_steps: int) -> OAState:
    """Integrate ``n_steps`` RK4 steps; return the final ``OAState``.

    Parameters
    ----------
    z0 : complex
        Initial complex Ott-Antonsen order parameter.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    OAState
        The final ``OAState`` after ``n_steps`` RK4 steps.
    """
    z0 = _validate_finite_complex(z0, name="z0")
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    re, im, r, psi = self._run_scalar(z0.real, z0.imag, n_steps)
    return OAState(z=complex(re, im), R=r, psi=psi, K_c=self.K_c)

predict_from_oscillators

predict_from_oscillators(
    omegas: FloatArray, K: float
) -> OAState

Fit a Lorentzian to omegas and return the relaxed OAState.

Parameters

omegas : FloatArray Natural frequencies in rad/s, shape (N,). K : float Global coupling strength.

Returns

OAState The relaxed OAState for the fitted Lorentzian.

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

def predict_from_oscillators(
    self,
    omegas: FloatArray,
    K: float,
) -> OAState:
    """Fit a Lorentzian to ``omegas`` and return the relaxed ``OAState``.

    Parameters
    ----------
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    K : float
        Global coupling strength.

    Returns
    -------
    OAState
        The relaxed ``OAState`` for the fitted Lorentzian.
    """
    omegas64 = _reduction_validation.validate_oa_frequency_sample(
        omegas,
        name="omegas",
    )
    K = _validate_finite_real(K, name="K")
    if _HAS_RUST_SCALAR:
        omega_0, delta = _rust_fit_lorentzian(omegas64)
    else:
        omega_0 = float(np.median(omegas64))
        q75, q25 = np.percentile(omegas64, [75, 25])
        delta = (q75 - q25) / 2.0 if q75 > q25 else 0.01
    reducer = OttAntonsenReduction(omega_0, delta, K, dt=self._dt)
    return reducer.run(complex(0.01, 0.0), n_steps=int(10.0 / self._dt))

Functions:

Second-Order Inertial Kuramoto (Power Grids)

The swing equation models power grid transient stability:

m_i θ̈_i + d_i θ̇_i = P_i + Σ_j K_ij sin(θ_j - θ_i)

where m_i is rotor inertia, d_i is damping, P_i is power injection (positive = generator, negative = load), K_ij is line susceptance.

Desynchronization = cascading blackout (Iberian Peninsula, April 2025).

from scpn_phase_orchestrator.upde.inertial import InertialKuramotoEngine

engine = InertialKuramotoEngine(n=100, dt=0.01)
theta, omega, theta_traj, omega_traj = engine.run(
    theta0, omega0, power, knm, inertia, damping, n_steps=10000
)
freq_dev = engine.frequency_deviation(omega)  # Hz deviation
R = engine.coherence(theta)  # phase coherence

inertial

Second-order (swing-equation) Kuramoto with a 5-backend fallback chain.

Model

Each oscillator has a phase θ_i and a "frequency-deviation" ω_i ≡ dθ_i/dt. The swing equation is

M_i · d²θ_i/dt² + D_i · dθ_i/dt = P_i + Σ_j K_ij · sin(θ_j − θ_i)

and is advanced with classical explicit RK4 on the (θ, ω) pair. This is the power-grid form used in Filatrella-Nielsen-Mallick 2008.

Numerics

The derivative uses the sin(θ_j − θ_i) = sin(θ_j)·cos(θ_i) − cos(θ_j)·sin(θ_i) expansion so that floating-point rounding matches the Rust kernel (spo-engine/src/inertial.rs) bit-for-bit. All five backends (Rust, Mojo, Julia, Go, Python) agree within ~1e-14 on the canonical all-to-all test problem; the dispatcher selects the fastest available path.

Classes

InertialKuramotoEngine

InertialKuramotoEngine(n: int, dt: float = 0.01)

Second-order swing-equation Kuramoto stepper with 5-backend dispatch.

The engine's geometry is (n, dt); the step itself is stateless: (θ, ω, P, K, M, D) → (θ', ω').

Initialise the stateless inertial Kuramoto stepper geometry.

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

def __init__(self, n: int, dt: float = 0.01) -> None:
    """Initialise the stateless inertial Kuramoto stepper geometry."""
    self._n = _validate_positive_int(n, name="n")
    self._dt = _validate_positive_float(dt, name="dt")

Methods:

step

step(
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
) -> tuple[FloatArray, FloatArray]

Advance one second-order inertial Kuramoto timestep.

Parameters

theta : FloatArray Oscillator phases in radians, shape (N,). omega_dot : FloatArray Instantaneous frequency deviations in rad/s, shape (N,). power : FloatArray Per-oscillator power injection in the swing equation, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). inertia : FloatArray Per-oscillator inertia coefficients, shape (N,). damping : FloatArray Per-oscillator damping coefficients, shape (N,).

Returns

tuple[FloatArray, FloatArray] The (θ, ω̇) state after one second-order step.

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

def step(
    self,
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
) -> tuple[FloatArray, FloatArray]:
    """Advance one second-order inertial Kuramoto timestep.

    Parameters
    ----------
    theta : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omega_dot : FloatArray
        Instantaneous frequency deviations in rad/s, shape ``(N,)``.
    power : FloatArray
        Per-oscillator power injection in the swing equation, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    inertia : FloatArray
        Per-oscillator inertia coefficients, shape ``(N,)``.
    damping : FloatArray
        Per-oscillator damping coefficients, shape ``(N,)``.

    Returns
    -------
    tuple[FloatArray, FloatArray]
        The ``(θ, ω̇)`` state after one second-order step.
    """
    theta64 = _validate_state_array(theta, name="theta", shape=(self._n,))
    omega_dot64 = _validate_state_array(
        omega_dot,
        name="omega_dot",
        shape=(self._n,),
    )
    power64 = _validate_state_array(power, name="power", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    inertia64 = _validate_positive_state_array(
        inertia,
        name="inertia",
        shape=(self._n,),
    )
    damping64 = _validate_nonnegative_state_array(
        damping,
        name="damping",
        shape=(self._n,),
    )
    knm_flat = knm64.ravel()
    backend_fn = _dispatch()
    if backend_fn is not None:
        new_theta, new_omega = backend_fn(
            theta64,
            omega_dot64,
            power64,
            knm_flat,
            inertia64,
            damping64,
            self._n,
            self._dt,
        )
        return _validate_backend_output(new_theta, new_omega, n=self._n)
    new_theta, new_omega = _python_step(
        theta64,
        omega_dot64,
        power64,
        knm_flat,
        inertia64,
        damping64,
        self._n,
        self._dt,
    )
    return _validate_backend_output(new_theta, new_omega, n=self._n)

run

run(
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
    n_steps: int,
) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray]

Integrate inertial Kuramoto dynamics and return final state plus traces.

Parameters

theta : FloatArray Oscillator phases in radians, shape (N,). omega_dot : FloatArray Instantaneous frequency deviations in rad/s, shape (N,). power : FloatArray Per-oscillator power injection in the swing equation, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). inertia : FloatArray Per-oscillator inertia coefficients, shape (N,). damping : FloatArray Per-oscillator damping coefficients, shape (N,). n_steps : int Number of integration steps to run.

Returns

tuple[FloatArray, FloatArray, FloatArray, FloatArray] The final (θ, ω̇) plus the θ and ω̇ traces.

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

def run(
    self,
    theta: FloatArray,
    omega_dot: FloatArray,
    power: FloatArray,
    knm: FloatArray,
    inertia: FloatArray,
    damping: FloatArray,
    n_steps: int,
) -> tuple[
    FloatArray,
    FloatArray,
    FloatArray,
    FloatArray,
]:
    """Integrate inertial Kuramoto dynamics and return final state plus traces.

    Parameters
    ----------
    theta : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omega_dot : FloatArray
        Instantaneous frequency deviations in rad/s, shape ``(N,)``.
    power : FloatArray
        Per-oscillator power injection in the swing equation, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    inertia : FloatArray
        Per-oscillator inertia coefficients, shape ``(N,)``.
    damping : FloatArray
        Per-oscillator damping coefficients, shape ``(N,)``.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    tuple[FloatArray, FloatArray, FloatArray, FloatArray]
        The final ``(θ, ω̇)`` plus the ``θ`` and ``ω̇`` traces.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    theta_traj = np.empty((n_steps, self._n))
    omega_traj = np.empty((n_steps, self._n))
    th = _validate_state_array(theta, name="theta", shape=(self._n,)).copy()
    od = _validate_state_array(
        omega_dot,
        name="omega_dot",
        shape=(self._n,),
    ).copy()
    for i in range(n_steps):
        th, od = self.step(th, od, power, knm, inertia, damping)
        theta_traj[i] = th
        omega_traj[i] = od
    return th, od, theta_traj, omega_traj

frequency_deviation

frequency_deviation(omega_dot: FloatArray) -> float

Return maximum absolute frequency deviation in cycles per unit time.

Parameters

omega_dot : FloatArray Instantaneous frequency deviations in rad/s, shape (N,).

Returns

float The maximum absolute frequency deviation in cycles per unit time.

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

def frequency_deviation(self, omega_dot: FloatArray) -> float:
    """Return maximum absolute frequency deviation in cycles per unit time.

    Parameters
    ----------
    omega_dot : FloatArray
        Instantaneous frequency deviations in rad/s, shape ``(N,)``.

    Returns
    -------
    float
        The maximum absolute frequency deviation in cycles per unit time.
    """
    omega_dot64 = _validate_state_array(
        omega_dot,
        name="omega_dot",
        shape=(self._n,),
    )
    return float(np.max(np.abs(omega_dot64)) / TWO_PI)

coherence

coherence(theta: FloatArray) -> float

Return the Kuramoto order parameter for the supplied phases.

Parameters

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

Returns

float The Kuramoto order parameter R.

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

def coherence(self, theta: FloatArray) -> float:
    """Return the Kuramoto order parameter for the supplied phases.

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

    Returns
    -------
    float
        The Kuramoto order parameter ``R``.
    """
    theta64 = _validate_state_array(theta, name="theta", shape=(self._n,))
    return float(np.abs(np.mean(np.exp(1j * theta64))))

Functions:

Financial Market Synchronization

Detect market regimes via Kuramoto order parameter on asset price phases. R(t) → 1 precedes market crashes (Black Monday 1987, 2008 crisis).

from scpn_phase_orchestrator.upde.market import (
    extract_phase, market_order_parameter, detect_regimes, sync_warning,
)

phases = extract_phase(returns_matrix)  # Hilbert transform
R = market_order_parameter(phases)       # R(t) across assets
regimes = detect_regimes(R)              # 0=desync, 1=transition, 2=sync
warnings = sync_warning(R, threshold=0.7)  # crash early warning

market

Kuramoto-based financial market synchronisation analysis.

Exposes a 5-backend fallback chain.

Extracts instantaneous phase from price / return time series via the Hilbert transform (scipy.signal.hilbert — FFT-based, stays Python-side because the Rust/Go/Mojo backends do not ship an FFT), then dispatches the two post-processing compute kernels:

• market_order_parameter(phases) — R(t) = |⟨exp(iθ)⟩_N| at every timestep. O(T · N).
• market_plv(phases, window) — rolling phase-locking-value matrix between assets, O((T − W + 1) · N² · W) with a sincos precompute that eliminates trig from the inner loop.

The detect_regimes classifier and sync_warning crossing detector are O(T) masking / comparison operations; they stay pure NumPy. R(t) → 1 preceded Black Monday 1987 and the 2008 crash (arXiv:1109.1167; CEUR-WS Vol-915).

Functions:

extract_phase

extract_phase(series: FloatArray) -> FloatArray

Extract instantaneous phase from a time series via the Hilbert transform.

Stays Python-side because the transform is FFT-based (scipy.signal.hilbert) and the compiled backends do not ship an FFT library.

Parameters

series : FloatArray Real-valued time series, shape (T,).

Returns

FloatArray The instantaneous phase of the series in [0, 2π).

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

def extract_phase(series: FloatArray) -> FloatArray:
    """Extract instantaneous phase from a time series via the Hilbert transform.

    Stays Python-side because the transform is FFT-based
    (``scipy.signal.hilbert``) and the compiled backends do not
    ship an FFT library.

    Parameters
    ----------
    series : FloatArray
        Real-valued time series, shape ``(T,)``.

    Returns
    -------
    FloatArray
        The instantaneous phase of the series in ``[0, 2π)``.
    """
    series = _validate_series(series)
    analytic = hilbert(series, axis=0)
    phase: FloatArray = np.angle(analytic) % (2.0 * np.pi)
    return phase

market_order_parameter

market_order_parameter(phases: FloatArray) -> FloatArray

Return the Kuramoto order parameter R(t) across N assets.

Parameters

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

Returns

FloatArray The Kuramoto order parameter time series R(t).

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

def market_order_parameter(phases: FloatArray) -> FloatArray:
    """Return the Kuramoto order parameter ``R(t)`` across ``N`` assets.

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

    Returns
    -------
    FloatArray
        The Kuramoto order parameter time series ``R(t)``.
    """
    phases = _validate_phase_matrix(phases)
    T, N = phases.shape
    if T == 0:
        return np.empty(0, dtype=np.float64)
    flat = np.ascontiguousarray(phases.ravel(), dtype=np.float64)
    dispatched = _dispatch()
    if dispatched is not None:
        op_fn, _ = dispatched
        return validate_market_order_output(op_fn(flat, T, N), t=T)
    return _python_market_order_parameter(flat, T, N)

market_plv

market_plv(
    phases: FloatArray, window: int = 50
) -> FloatArray

Compute the rolling phase-locking-value matrix between assets.

Returns shape (T − window + 1, N, N).

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). window : int Sliding-window length in samples.

Returns

FloatArray The rolling phase-locking-value matrices, shape (T − window + 1, N, N).

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

def market_plv(phases: FloatArray, window: int = 50) -> FloatArray:
    """Compute the rolling phase-locking-value matrix between assets.

    Returns shape ``(T − window + 1, N, N)``.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    window : int
        Sliding-window length in samples.

    Returns
    -------
    FloatArray
        The rolling phase-locking-value matrices, shape ``(T − window + 1, N, N)``.
    """
    phases = _validate_phase_matrix(phases)
    window = _validate_positive_int(window, name="window")
    T, N = phases.shape
    if window > T or N == 0:
        return np.empty((0, N, N), dtype=np.float64)
    flat = np.ascontiguousarray(phases.ravel(), dtype=np.float64)
    dispatched = _dispatch()
    if dispatched is not None:
        _, plv_fn = dispatched
        result_flat = validate_market_plv_output(
            plv_fn(flat, T, N, window),
            t=T,
            n=N,
            window=window,
        )
    else:
        result_flat = _python_market_plv(flat, T, N, window)
    n_windows = T - window + 1
    return result_flat.reshape(n_windows, N, N)

detect_regimes

detect_regimes(
    R: FloatArray,
    sync_threshold: float = 0.7,
    desync_threshold: float = 0.3,
) -> IntArray

Classify market synchronisation regimes from R(t).

Returns int32 labels: 0 = desynchronised, 1 = transition, 2 = synchronised. O(T) masking; no multi-language port needed.

Parameters

R : FloatArray Order-parameter time series R(t), shape (T,). sync_threshold : float Order parameter above which the market is classed as synchronised. desync_threshold : float Order parameter below which the market is classed as desynchronised.

Returns

IntArray The per-timestep market regime labels.

Raises

ValueError If the thresholds are inconsistent or R is not 1-D.

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

def detect_regimes(
    R: FloatArray,
    sync_threshold: float = 0.7,
    desync_threshold: float = 0.3,
) -> IntArray:
    """Classify market synchronisation regimes from ``R(t)``.

    Returns ``int32`` labels: 0 = desynchronised, 1 = transition,
    2 = synchronised. O(T) masking; no multi-language port needed.

    Parameters
    ----------
    R : FloatArray
        Order-parameter time series ``R(t)``, shape ``(T,)``.
    sync_threshold : float
        Order parameter above which the market is classed as synchronised.
    desync_threshold : float
        Order parameter below which the market is classed as desynchronised.

    Returns
    -------
    IntArray
        The per-timestep market regime labels.

    Raises
    ------
    ValueError
        If the thresholds are inconsistent or ``R`` is not 1-D.
    """
    R = _validate_signal_vector(R, name="R")
    sync_threshold = _validate_finite_float(
        sync_threshold,
        name="sync_threshold",
    )
    desync_threshold = _validate_finite_float(
        desync_threshold,
        name="desync_threshold",
    )
    if sync_threshold < desync_threshold:
        raise ValueError(
            "sync_threshold must be greater than or equal to desync_threshold",
        )
    try:
        from spo_kernel import detect_regimes_rust as _rust_regimes

        flat = np.ascontiguousarray(R.ravel())
        return np.asarray(_rust_regimes(flat, sync_threshold, desync_threshold))
    except ImportError:
        pass
    regimes = np.ones(len(R), dtype=np.int32)
    mask_sync = sync_threshold <= R
    mask_desync = desync_threshold >= R
    regimes[mask_sync] = 2
    regimes[mask_desync] = 0
    return regimes

sync_warning

sync_warning(
    R: FloatArray,
    threshold: float = 0.7,
    lookback: int = 10,
) -> BoolArray

Detect synchronisation warnings where smoothed R crosses up.

Parameters

R : FloatArray Order-parameter time series R(t), shape (T,). threshold : float Decision threshold. lookback : int Number of past samples smoothed over before the crossing test.

Returns

BoolArray A per-timestep boolean mask of synchronisation warnings.

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

def sync_warning(
    R: FloatArray,
    threshold: float = 0.7,
    lookback: int = 10,
) -> BoolArray:
    """Detect synchronisation warnings where smoothed ``R`` crosses up.

    Parameters
    ----------
    R : FloatArray
        Order-parameter time series ``R(t)``, shape ``(T,)``.
    threshold : float
        Decision threshold.
    lookback : int
        Number of past samples smoothed over before the crossing test.

    Returns
    -------
    BoolArray
        A per-timestep boolean mask of synchronisation warnings.
    """
    R = _validate_signal_vector(R, name="R")
    threshold = _validate_finite_float(threshold, name="threshold")
    lookback = _validate_positive_int(lookback, name="lookback")
    if lookback > 1:
        kernel = np.ones(lookback) / lookback
        R_smooth = np.convolve(R, kernel, mode="same")
    else:
        R_smooth = R
    warnings = np.zeros(len(R), dtype=bool)
    for t in range(1, len(R)):
        if R_smooth[t] >= threshold and R_smooth[t - 1] < threshold:
            warnings[t] = True
    return warnings

Swarmalator Dynamics

Agents that are simultaneously self-propelled particles AND phase oscillators. Phase modulates spatial attraction; proximity modulates phase coupling.

Five collective states emerge depending on J and K: - J > 0, K > 0: static sync (clustered, phase-locked) - J > 0, K < 0: static async (clustered, anti-phase) - J < 0, K > 0: static phase wave (spatially ordered by phase) - J < 0, K < 0: splintered phase wave - |J| ≈ 0: active phase wave (rotating)

from scpn_phase_orchestrator.upde.swarmalator import SwarmalatorEngine

engine = SwarmalatorEngine(n=50, dim=2, dt=0.01, A=1.0, B=1.0, J=0.5, K=1.0)
pos, phases, pos_traj, phase_traj = engine.run(
    positions0, phases0, omegas, n_steps=5000
)

# Metrics
R = engine.phase_coherence(phases)
compactness = engine.spatial_coherence(pos)
corr = engine.phase_spatial_correlation(pos, phases)

O'Keeffe, Hong, Strogatz, Nature Communications 2017. Experimental: Nature Communications Dec 2025 (colloidal system).

swarmalator

Swarmalator step (position + phase) with a 5-backend fallback chain.

Swarmalators combine spatial attraction / repulsion with phase oscillator dynamics (O'Keeffe, Hong & Strogatz, Nat. Commun. 8:1504, 2017). Each agent has a position x_i ∈ ℝ^d and a phase θ_i; they co-evolve through attract/repulse + phase-coupling terms:

ẋ_i = (1/N) Σ_j (x_j − x_i) [(a + j·cos(θ_j − θ_i)) / |x_j − x_i|
                             − b / |x_j − x_i|²]
θ̇_i = ω_i + (k / N) Σ_j sin(θ_j − θ_i) / |x_j − x_i|

The repulsion b·(x_j − x_i) / |x_j − x_i|² is the canonical inverse-distance hard core of O'Keeffe-Hong-Strogatz (magnitude b / |x_j − x_i|), with a = A = 1, b = B = 1, j = J, k = K recovering the original model. A single regularisation constant ε = 1e-6 is added to |x_j − x_i|² (and inside the sqrt for the attraction/phase |x_j − x_i|) so the kernel is finite at coincident agents; it vanishes in the ε → 0 limit.

Classes

SwarmalatorEngine

SwarmalatorEngine(
    n_agents: int, dim: int = 2, dt: float = 0.01
)

Swarmalator stepper with 5-backend dispatch.

The engine is stateful in its (n_agents, dim, dt) geometry but the step contract is stateless: (pos, phases, omegas) → (new_pos, new_phases).

Initialise the stateless swarmalator stepper geometry.

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

def __init__(self, n_agents: int, dim: int = 2, dt: float = 0.01) -> None:
    """Initialise the stateless swarmalator stepper geometry."""
    self._n = _validate_positive_int(n_agents, name="n_agents")
    self._dim = _validate_positive_int(dim, name="dim")
    self._dt = _validate_positive_float(dt, name="dt")

Methods:

step

step(
    pos: FloatArray,
    phases: FloatArray,
    omegas: FloatArray,
    a: float = 1.0,
    b: float = 1.0,
    j: float = 1.0,
    k: float = 1.0,
) -> tuple[FloatArray, FloatArray]

Advance coupled swarmalator positions and phases by one step.

Parameters

pos Agent positions with shape (n_agents, dim). phases Agent phases in radians, shape (n_agents,). omegas Natural angular frequencies, shape (n_agents,). a Baseline spatial attraction coefficient. b Spatial repulsion coefficient. j Phase-dependent attraction modulation. k Phase-coupling coefficient.

Returns

tuple[FloatArray, FloatArray] Updated positions with shape (n_agents, dim) and updated phases wrapped into [0, 2*pi).

Notes

The dispatcher selects the first available accelerated backend and falls back to the NumPy reference path with the same state contract.

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

def step(
    self,
    pos: FloatArray,
    phases: FloatArray,
    omegas: FloatArray,
    a: float = 1.0,
    b: float = 1.0,
    j: float = 1.0,
    k: float = 1.0,
) -> tuple[FloatArray, FloatArray]:
    """Advance coupled swarmalator positions and phases by one step.

    Parameters
    ----------
    pos
        Agent positions with shape ``(n_agents, dim)``.
    phases
        Agent phases in radians, shape ``(n_agents,)``.
    omegas
        Natural angular frequencies, shape ``(n_agents,)``.
    a
        Baseline spatial attraction coefficient.
    b
        Spatial repulsion coefficient.
    j
        Phase-dependent attraction modulation.
    k
        Phase-coupling coefficient.

    Returns
    -------
    tuple[FloatArray, FloatArray]
        Updated positions with shape ``(n_agents, dim)`` and updated
        phases wrapped into ``[0, 2*pi)``.

    Notes
    -----
    The dispatcher selects the first available accelerated backend and
    falls back to the NumPy reference path with the same state contract.
    """
    pos64 = _validate_state_array(
        pos,
        name="pos",
        shape=(self._n, self._dim),
    )
    phases64 = _validate_state_array(
        phases,
        name="phases",
        shape=(self._n,),
    )
    omegas64 = _validate_state_array(
        omegas,
        name="omegas",
        shape=(self._n,),
    )
    a = _validate_finite_float(a, name="a")
    b = _validate_finite_float(b, name="b")
    j = _validate_finite_float(j, name="j")
    k = _validate_finite_float(k, name="k")

    backend_fn = _dispatch()
    if backend_fn is not None:
        return _validate_backend_output(
            *backend_fn(
                pos64,
                phases64,
                omegas64,
                self._n,
                self._dim,
                a,
                b,
                j,
                k,
                self._dt,
            ),
            n=self._n,
            dim=self._dim,
        )
    new_pos, new_phases = _python_step(
        pos64,
        phases64,
        omegas64,
        self._n,
        self._dim,
        a,
        b,
        j,
        k,
        self._dt,
    )
    return _validate_backend_output(new_pos, new_phases, n=self._n, dim=self._dim)

run

run(
    pos: FloatArray,
    phases: FloatArray,
    omegas: FloatArray,
    a: float = 1.0,
    b: float = 1.0,
    j: float = 1.0,
    k: float = 1.0,
    n_steps: int = 100,
) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray]

Integrate swarmalator positions and phases with trajectory capture.

Parameters

pos : FloatArray Swarmalator positions, shape (N, 2). phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). a : float Spatial attraction strength. b : float Spatial repulsion strength. j : float Phase-to-space coupling strength. k : float Space-to-phase coupling strength. n_steps : int Number of integration steps to run.

Returns

tuple[FloatArray, FloatArray, FloatArray, FloatArray] The final positions and phases plus their trajectory traces.

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

def run(
    self,
    pos: FloatArray,
    phases: FloatArray,
    omegas: FloatArray,
    a: float = 1.0,
    b: float = 1.0,
    j: float = 1.0,
    k: float = 1.0,
    n_steps: int = 100,
) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray]:
    """Integrate swarmalator positions and phases with trajectory capture.

    Parameters
    ----------
    pos : FloatArray
        Swarmalator positions, shape ``(N, 2)``.
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    a : float
        Spatial attraction strength.
    b : float
        Spatial repulsion strength.
    j : float
        Phase-to-space coupling strength.
    k : float
        Space-to-phase coupling strength.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    tuple[FloatArray, FloatArray, FloatArray, FloatArray]
        The final positions and phases plus their trajectory traces.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    curr_pos = _validate_state_array(
        pos,
        name="pos",
        shape=(self._n, self._dim),
    ).copy()
    curr_phases = _validate_state_array(
        phases,
        name="phases",
        shape=(self._n,),
    ).copy()
    omegas64 = _validate_state_array(
        omegas,
        name="omegas",
        shape=(self._n,),
    )
    pos_traj = np.empty((n_steps, self._n, self._dim))
    phase_traj = np.empty((n_steps, self._n))
    for i in range(n_steps):
        curr_pos, curr_phases = self.step(
            curr_pos,
            curr_phases,
            omegas64,
            a,
            b,
            j,
            k,
        )
        pos_traj[i] = curr_pos
        phase_traj[i] = curr_phases
    return curr_pos, curr_phases, pos_traj, phase_traj

order_parameter

order_parameter(phases: FloatArray) -> float

Return the Kuramoto order parameter for swarmalator phases.

Parameters

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

Returns

float The Kuramoto order parameter R.

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

def order_parameter(self, phases: FloatArray) -> float:
    """Return the Kuramoto order parameter for swarmalator phases.

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

    Returns
    -------
    float
        The Kuramoto order parameter ``R``.
    """
    phases64 = _validate_state_array(
        phases,
        name="phases",
        shape=(self._n,),
    )
    return float(np.abs(np.mean(np.exp(1j * phases64))))

Functions: