UPDE Engine

The Unified Phase Dynamics Engine (UPDE) is SPO's core integrator subsystem. It provides 19 ODE engine variants covering standard Kuramoto, Bayesian uncertainty propagation, amplitude dynamics (Stuart-Landau), higher-order interactions (simplicial), inertial systems (power grids), stochastic resonance, geometric integration, time delays, financial markets, spatial-phase coupling (swarmalators), hypergraph k-body coupling, mean-field reduction, variational prediction, adjoint gradients, and bifurcation continuation.

Pipeline position

CouplingBuilder.build() ──→ K_nm, α
                                │
Oscillators.extract() ──→ θ, ω │
                                │
Drivers.compute() ──→ Ψ        │
                                ↓
              ┌─────── UPDEEngine.step(θ, ω, K, ζ, Ψ, α) ───────┐
              │                                                    │
              │    Euler / RK4 / RK45 (adaptive)                   │
              │    Optional: Rust FFI via spo_kernel                │
              │                                                    │
              └────────────── θ_new ∈ [0, 2π)^N ─────────────────┘
                                │
                                ↓
                     compute_order_parameter(θ) → R, ψ
                                │
             BayesianUPDE → posterior predictive R ± sigma
                                │
                                ↓
                     RegimeManager.evaluate() → Regime

The engine is the computational core of SPO. Every subsystem feeds into it (coupling, oscillators, drivers) or consumes its output (order parameters, monitors, supervisor).

Engine variants

Engine	State	ODE	Use case
UPDEEngine	θ ∈ [0,2π)^N	Kuramoto	General synchronisation
BayesianUPDE	θ plus sampled K,ω	Monte Carlo UPDE	Safety-tier uncertainty quantification
SparseUPDEEngine	θ ∈ [0,2π)^N	Sparse Kuramoto	High-N scalability (\(O(N \log N)\))
SheafUPDEEngine	\(\vec{\theta} \in \mathbb{R}^{N \times D}\)	Cellular Sheaf	Multi-dimensional block coupling
StuartLandauEngine	[θ,r] ∈ R^{2N}	Stuart-Landau	Amplitude dynamics
SimplicialEngine	θ ∈ [0,2π)^N	3-body Kuramoto	Triadic/group synchronization
InertialEngine	[θ,ω̇] ∈ R^{2N}	Swing equation	Power grids
SwarmalatorEngine	[x,θ] ∈ R^{(D+1)N}	Position + phase	Swarm robotics
StochasticInjector	θ ∈ [0,2π)^N	Euler-Maruyama	Noise resonance
GeometricEngine	z ∈ C^N	SO(2) exponential	Long simulations
DelayedEngine	θ + buffer	Delayed Kuramoto	Transport delays
MarketEngine	θ from Hilbert	Price → phase	Financial markets
SplittingEngine	θ ∈ [0,2π)^N	Symplectic split	Energy-preserving
HypergraphEngine	θ ∈ [0,2π)^N	k-body coupling	Mixed-order
OttAntosenReduction	z ∈ C	Mean-field ODE	Fast prediction
PredictionModel	θ ∈ [0,2π)^N	Error injection	FEP-Kuramoto
AdjointGradient	∂R/∂K	Finite diff / JAX	Optimisation

Performance budgets

Operation	N	Budget	Rust path
UPDEEngine.step()	8	< 50 μs	~ 30 μs
UPDEEngine.step()	64	< 1 ms	~ 0.3 ms
UPDEEngine.step()	128	< 5 ms	~ 1 ms
compute_order_parameter()	256	< 100 μs	~ 2 μs
StuartLandauEngine.step()	32	< 2 ms	—
SplittingEngine.step()	64	< 1 ms	—
DelayedEngine.step()	32	< 1 ms	—

Core Kuramoto Engine

First-order Kuramoto ODE: dθ_i/dt = ω_i + Σ_j K_ij sin(θ_j - θ_i - α_ij) + ζ sin(Ψ - θ_i). Supports Euler, RK4, and RK45 (adaptive) integration. Optional Rust FFI acceleration via spo_kernel.PyUPDEStepper.

Direct Go, Julia, and Mojo accelerator entrypoints share the same boundary contract before optional runtime loading: phase and frequency vectors must be finite real one-dimensional float64 arrays with matching length; coupling and phase-lag matrices must be finite real square matrices (or flattened square matrices) matching oscillator count; the coupling diagonal must be exactly zero to exclude self-coupling; dt, atol, and rtol must be positive finite scalars; n_steps must be a non-negative integer; and n_substeps must be a positive integer. A zero-step direct call returns a copy of the initial phase vector without requiring the optional backend binary or runtime. Mojo subprocess output must contain exactly one raw stdout line per oscillator phase; blank, truncated, or overlong output is rejected before final phase validation.

The public stateless upde_run() and upde_run_omega_schedule() entrypoints use that same core-owned contract before backend selection. Boolean, complex, and numeric-string aliases are rejected before conversion for phase, frequency, coupling, phase-lag, schedule, scalar-control, and count inputs. Python and optional-backend results then pass through one finite real-vector cardinality check before publication; Julia returns retain their source dtype until this check, while Mojo stdout remains an explicitly parsed text protocol.

engine

Stateful :class:UPDEEngine.

The batched integrator is stateless — see :mod:scpn_phase_orchestrator.upde._run and the re-exported :func:upde_run. This module keeps the state-heavy observer: the class pre-allocates scratch buffers for the chosen method, holds a reentrant lock for thread-safety, and retains _last_dt across RK45 step calls.

Bayesian UPDE Uncertainty Propagation

Samples natural frequencies and coupling matrices from explicit distributions, runs the existing UPDE kernel for each draw, and reports posterior-predictive R ± sigma with credible intervals and audit diagnostics.

Public phase, frequency, coupling, phase-lag, posterior-fit, and Gaussian distribution arrays reject boolean, complex, and numeric-string aliases before conversion while preserving real numeric-object arrays. Custom distribution samples replay the same source-type, shape, and finiteness checks before Monte Carlo execution, and drive controls must be finite real scalars. Reserved NumPyro and BlackJAX names remain explicitly fail-closed.

bayesian

Uncertainty propagation for Kuramoto UPDE rollouts.

The shipped backend is deterministic NumPy Monte Carlo over explicit distributions for omega and K_nm. Probabilistic-programming backends are reserved as fail-closed names until their samplers are implemented and benchmarked against this reproducible baseline.

Classes

GaussianArrayDistribution dataclass

GaussianArrayDistribution(
    mean: object,
    std: object,
    non_negative: bool = False,
    zero_diagonal: bool = False,
)

Independent Gaussian array distribution with optional matrix guards.

Attributes

shape property

shape: tuple[int, ...]

Return the event shape sampled by this Gaussian distribution.

Returns

tuple[int, ...] Return the event shape sampled by this Gaussian distribution.

Methods:

sample

sample(rng: Generator, n_samples: int) -> FloatArray

Draw finite Gaussian samples with optional support guards applied.

Parameters

rng : np.random.Generator NumPy random generator used for sampling. n_samples : int Number of samples to draw.

Returns

FloatArray Finite Gaussian samples, shape (n_samples, *event_shape).

Raises

TypeError If rng is not a NumPy random generator.

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

def sample(self, rng: np.random.Generator, n_samples: int) -> FloatArray:
    """Draw finite Gaussian samples with optional support guards applied.

    Parameters
    ----------
    rng : np.random.Generator
        NumPy random generator used for sampling.
    n_samples : int
        Number of samples to draw.

    Returns
    -------
    FloatArray
        Finite Gaussian samples, shape ``(n_samples, *event_shape)``.

    Raises
    ------
    TypeError
        If ``rng`` is not a NumPy random generator.
    """
    if not isinstance(rng, np.random.Generator):
        raise TypeError("rng must be a numpy.random.Generator")
    sample_count = _validate_positive_integer(
        n_samples,
        name="n_samples",
        minimum=1,
    )
    draws = rng.normal(
        loc=np.asarray(self.mean, dtype=np.float64),
        scale=np.asarray(self.std, dtype=np.float64),
        size=(sample_count, *self.shape),
    )
    samples = np.asarray(draws, dtype=np.float64)
    if self.non_negative:
        samples = np.maximum(samples, 0.0)
    if self.zero_diagonal:
        diag = np.arange(self.shape[0])
        samples[:, diag, diag] = 0.0
    return np.ascontiguousarray(samples, dtype=np.float64)

BayesianUPDEConfig dataclass

BayesianUPDEConfig(
    n_samples: int = 128,
    seed: int | None = None,
    dt: float = 0.01,
    n_steps: int = 1,
    method: MethodName = "rk4",
    credible_interval: float = 0.95,
    backend: BackendName = "numpy",
    n_substeps: int = 1,
    atol: float = 1e-06,
    rtol: float = 0.001,
)

Configuration for Bayesian UPDE uncertainty propagation.

BayesianUPDEResult dataclass

BayesianUPDEResult(
    r_samples: FloatArray,
    final_phase_samples: FloatArray,
    omega_mean: FloatArray,
    knm_mean: FloatArray,
    r_mean: float,
    r_sigma: float,
    r_lower: float,
    r_upper: float,
    psi_mean: float,
    sample_count: int,
    credible_interval: float,
    backend: str,
    method: str,
)

Posterior predictive order-parameter summary.

Attributes

r_plus_minus property

r_plus_minus: tuple[float, float]

Return the compact R ± sigma pair.

Returns

tuple[float, float] Return the compact R ± sigma pair.

Methods:

to_audit_record

to_audit_record() -> dict[str, object]

Return JSON-safe uncertainty diagnostics.

Returns

dict[str, object] Return JSON-safe uncertainty diagnostics.

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

def to_audit_record(self) -> dict[str, object]:
    """Return JSON-safe uncertainty diagnostics.

    Returns
    -------
    dict[str, object]
        Return JSON-safe uncertainty diagnostics.
    """
    return {
        "kind": "bayesian_upde",
        "backend": self.backend,
        "method": self.method,
        "sample_count": self.sample_count,
        "credible_interval": self.credible_interval,
        "r_summary": {
            "mean": self.r_mean,
            "sigma": self.r_sigma,
            "lower": self.r_lower,
            "upper": self.r_upper,
            "plus_minus": [self.r_mean, self.r_sigma],
        },
        "psi_mean": self.psi_mean,
        "omega_mean": self.omega_mean.tolist(),
        "knm_mean": self.knm_mean.tolist(),
        "final_phase_mean": np.mean(self.final_phase_samples, axis=0).tolist(),
        "diagnostics": {
            "finite_samples": bool(
                np.all(np.isfinite(self.r_samples))
                and np.all(np.isfinite(self.final_phase_samples))
            ),
            "r_min": float(np.min(self.r_samples)),
            "r_max": float(np.max(self.r_samples)),
        },
    }

BayesianBackendStatus dataclass

BayesianBackendStatus(
    backend: str,
    available: bool,
    fail_closed: bool,
    reason: str,
    sample_count: int,
)

Execution status for one Bayesian UPDE backend name.

Methods:

to_audit_record

to_audit_record() -> dict[str, object]

Return JSON-safe backend availability diagnostics.

Returns

dict[str, object] Return JSON-safe backend availability diagnostics.

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

def to_audit_record(self) -> dict[str, object]:
    """Return JSON-safe backend availability diagnostics.

    Returns
    -------
    dict[str, object]
        Return JSON-safe backend availability diagnostics.
    """
    return {
        "kind": "bayesian_backend_status",
        "backend": self.backend,
        "available": self.available,
        "fail_closed": self.fail_closed,
        "reason": self.reason,
        "sample_count": self.sample_count,
    }

GaussianUPDEPosteriorFit dataclass

GaussianUPDEPosteriorFit(
    omega: GaussianArrayDistribution,
    knm: GaussianArrayDistribution,
    residual_rmse: float,
    sample_count: int,
    dt: float,
    ridge: float,
    backend: str = "numpy_lstsq",
)

Gaussian posterior approximation fitted from observed phase trajectories.

Methods:

to_audit_record

to_audit_record() -> dict[str, object]

Return JSON-safe posterior-fit diagnostics.

Returns

dict[str, object] Return JSON-safe posterior-fit diagnostics.

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

def to_audit_record(self) -> dict[str, object]:
    """Return JSON-safe posterior-fit diagnostics.

    Returns
    -------
    dict[str, object]
        Return JSON-safe posterior-fit diagnostics.
    """
    return {
        "kind": "gaussian_upde_posterior_fit",
        "backend": self.backend,
        "sample_count": self.sample_count,
        "dt": self.dt,
        "ridge": self.ridge,
        "residual_rmse": self.residual_rmse,
        "omega": {
            "mean": np.asarray(self.omega.mean, dtype=np.float64).tolist(),
            "std": np.asarray(self.omega.std, dtype=np.float64).tolist(),
        },
        "knm": {
            "mean": np.asarray(self.knm.mean, dtype=np.float64).tolist(),
            "std": np.asarray(self.knm.std, dtype=np.float64).tolist(),
        },
        "diagnostics": {
            "finite": bool(
                np.all(np.isfinite(np.asarray(self.omega.mean)))
                and np.all(np.isfinite(np.asarray(self.omega.std)))
                and np.all(np.isfinite(np.asarray(self.knm.mean)))
                and np.all(np.isfinite(np.asarray(self.knm.std)))
                and np.isfinite(self.residual_rmse)
            ),
            "zero_diagonal": bool(
                np.allclose(np.diag(np.asarray(self.knm.mean)), 0.0)
                and np.allclose(np.diag(np.asarray(self.knm.std)), 0.0)
            ),
        },
    }

Functions:

fit_gaussian_upde_posterior

fit_gaussian_upde_posterior(
    phase_trajectory: object,
    *,
    dt: float,
    alpha: object | None = None,
    ridge: float = 1e-06,
    coupling_std_floor: float = 1e-06,
    omega_std_floor: float = 1e-06,
) -> GaussianUPDEPosteriorFit

Fit Gaussian omega and K_nm priors from observed phases.

The estimator is a deterministic NumPy ridge least-squares baseline. It fits the Kuramoto right-hand side independently per target oscillator:

d theta_i / dt = omega_i + sum_j K_ij sin(theta_j - theta_i - alpha_ij).

The result is intentionally review-only: it produces distributions that can feed :func:bayesian_upde_run, but it does not apply control actions.

Parameters

phase_trajectory : object Observed phase trajectory, shape (T, N). dt : float Integration step size. alpha : object | None Phase-lag matrix in radians, shape (N, N), or None for no lag. ridge : float Ridge-regularisation strength for the prior fit. coupling_std_floor : float Lower bound on the fitted coupling standard deviation. omega_std_floor : float Lower bound on the fitted natural-frequency standard deviation.

Returns

GaussianUPDEPosteriorFit The fitted Gaussian omega and K_nm prior.

Raises

ValueError If the phase trajectory is empty or non-finite.

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

def fit_gaussian_upde_posterior(
    phase_trajectory: object,
    *,
    dt: float,
    alpha: object | None = None,
    ridge: float = 1e-6,
    coupling_std_floor: float = 1e-6,
    omega_std_floor: float = 1e-6,
) -> GaussianUPDEPosteriorFit:
    """Fit Gaussian ``omega`` and ``K_nm`` priors from observed phases.

    The estimator is a deterministic NumPy ridge least-squares baseline. It
    fits the Kuramoto right-hand side independently per target oscillator:

    ``d theta_i / dt = omega_i + sum_j K_ij sin(theta_j - theta_i - alpha_ij)``.

    The result is intentionally review-only: it produces distributions that can
    feed :func:`bayesian_upde_run`, but it does not apply control actions.

    Parameters
    ----------
    phase_trajectory : object
        Observed phase trajectory, shape ``(T, N)``.
    dt : float
        Integration step size.
    alpha : object | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    ridge : float
        Ridge-regularisation strength for the prior fit.
    coupling_std_floor : float
        Lower bound on the fitted coupling standard deviation.
    omega_std_floor : float
        Lower bound on the fitted natural-frequency standard deviation.

    Returns
    -------
    GaussianUPDEPosteriorFit
        The fitted Gaussian ``omega`` and ``K_nm`` prior.

    Raises
    ------
    ValueError
        If the phase trajectory is empty or non-finite.
    """
    trajectory = _as_finite_array(phase_trajectory, name="phase_trajectory")
    if trajectory.ndim != 2:
        raise ValueError(
            f"phase_trajectory must be a 2-D array, got shape {trajectory.shape}"
        )
    if trajectory.shape[0] < 3:
        raise ValueError("phase_trajectory must contain at least three samples")
    dt_value = _validate_positive_finite(dt, name="dt")
    ridge_value = _validate_non_negative_finite(ridge, name="ridge")
    coupling_floor = _validate_non_negative_finite(
        coupling_std_floor,
        name="coupling_std_floor",
    )
    omega_floor = _validate_non_negative_finite(
        omega_std_floor,
        name="omega_std_floor",
    )
    n_samples, n_oscillators = trajectory.shape
    alpha_array = (
        np.zeros((n_oscillators, n_oscillators), dtype=np.float64)
        if alpha is None
        else _as_finite_array(alpha, name="alpha")
    )
    if alpha_array.shape != (n_oscillators, n_oscillators):
        raise ValueError(
            f"alpha must have shape {(n_oscillators, n_oscillators)}, "
            f"got {alpha_array.shape}"
        )

    unwrapped = np.unwrap(trajectory, axis=0)
    derivatives = np.diff(unwrapped, axis=0) / dt_value
    theta = trajectory[:-1]
    omega_mean = np.empty(n_oscillators, dtype=np.float64)
    omega_std = np.empty(n_oscillators, dtype=np.float64)
    knm_mean = np.zeros((n_oscillators, n_oscillators), dtype=np.float64)
    knm_std = np.zeros((n_oscillators, n_oscillators), dtype=np.float64)
    residuals: list[FloatArray] = []

    for target in range(n_oscillators):
        source_indices = [source for source in range(n_oscillators) if source != target]
        features = np.column_stack(
            [
                np.ones(theta.shape[0], dtype=np.float64),
                *[
                    np.sin(
                        theta[:, source]
                        - theta[:, target]
                        - alpha_array[target, source]
                    )
                    for source in source_indices
                ],
            ]
        )
        target_derivative = derivatives[:, target]
        gram = features.T @ features
        penalty = ridge_value * np.eye(gram.shape[0], dtype=np.float64)
        penalty[0, 0] = 0.0
        coeffs = np.linalg.solve(gram + penalty, features.T @ target_derivative)
        predicted = features @ coeffs
        residual = target_derivative - predicted
        residuals.append(residual)
        dof = max(1, features.shape[0] - features.shape[1])
        residual_sigma = float(np.sqrt(float(residual @ residual) / dof))
        covariance = residual_sigma**2 * np.linalg.pinv(gram + penalty)
        coefficient_std = np.sqrt(np.maximum(np.diag(covariance), 0.0))
        omega_mean[target] = coeffs[0]
        omega_std[target] = max(float(coefficient_std[0]), omega_floor)
        for offset, source in enumerate(source_indices, start=1):
            value = max(float(coeffs[offset]), 0.0)
            knm_mean[target, source] = value
            knm_std[target, source] = (
                max(float(coefficient_std[offset]), coupling_floor)
                if value > 0.0
                else 0.0
            )

    residual_vector = np.concatenate(residuals)
    residual_rmse = float(np.sqrt(float(np.mean(residual_vector**2))))
    return GaussianUPDEPosteriorFit(
        omega=GaussianArrayDistribution(omega_mean, omega_std),
        knm=GaussianArrayDistribution(
            knm_mean,
            knm_std,
            non_negative=True,
            zero_diagonal=True,
        ),
        residual_rmse=residual_rmse,
        sample_count=n_samples,
        dt=dt_value,
        ridge=ridge_value,
    )

audit_bayesian_backend_status

audit_bayesian_backend_status(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
    backends: tuple[BackendName, ...] = (
        "numpy",
        "numpyro",
        "blackjax",
    ),
) -> tuple[BayesianBackendStatus, ...]

Probe Bayesian backend names without silently accepting unsupported ones.

Parameters

phases : object Oscillator phases in radians, shape (N,). omega : object Natural-frequency distribution or array. knm : object Coupling matrix K_nm, shape (N, N). alpha : object Phase-lag matrix in radians, shape (N, N). Use a zero matrix for no lag. zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. config : BayesianUPDEConfig | None Optional configuration object, or None for defaults. backends : tuple[BackendName, ...] Backend names to probe, in priority order.

Returns

tuple[BayesianBackendStatus, ...] Per-backend availability diagnostics.

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

def audit_bayesian_backend_status(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
    backends: tuple[BackendName, ...] = ("numpy", "numpyro", "blackjax"),
) -> tuple[BayesianBackendStatus, ...]:
    """Probe Bayesian backend names without silently accepting unsupported ones.

    Parameters
    ----------
    phases : object
        Oscillator phases in radians, shape ``(N,)``.
    omega : object
        Natural-frequency distribution or array.
    knm : object
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    alpha : object
        Phase-lag matrix in radians, shape ``(N, N)``. Use a zero matrix for no
        lag.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    config : BayesianUPDEConfig | None
        Optional configuration object, or ``None`` for defaults.
    backends : tuple[BackendName, ...]
        Backend names to probe, in priority order.

    Returns
    -------
    tuple[BayesianBackendStatus, ...]
        Per-backend availability diagnostics.
    """
    base_config = config or BayesianUPDEConfig(n_samples=8, seed=0, n_steps=1)
    statuses: list[BayesianBackendStatus] = []
    for backend in backends:
        backend_config = replace(base_config, backend=backend)
        try:
            result = bayesian_upde_run(
                phases,
                omega=omega,
                knm=knm,
                alpha=alpha,
                zeta=zeta,
                psi=psi,
                config=backend_config,
            )
        except NotImplementedError as exc:
            statuses.append(
                BayesianBackendStatus(
                    backend=backend,
                    available=False,
                    fail_closed=True,
                    reason=str(exc),
                    sample_count=0,
                )
            )
        else:
            statuses.append(
                BayesianBackendStatus(
                    backend=backend,
                    available=True,
                    fail_closed=False,
                    reason="executed",
                    sample_count=result.sample_count,
                )
            )
    return tuple(statuses)

bayesian_upde_run

bayesian_upde_run(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
) -> BayesianUPDEResult

Run UPDE over sampled omega and K_nm distributions.

Parameters

phases : object Oscillator phases in radians, shape (N,). omega : object Natural-frequency distribution or array. knm : object Coupling matrix K_nm, shape (N, N). alpha : object Phase-lag matrix in radians, shape (N, N). Use a zero matrix for no lag. zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. config : BayesianUPDEConfig | None Optional configuration object, or None for defaults.

Returns

BayesianUPDEResult The Bayesian UPDE result with uncertainty diagnostics.

Raises

NotImplementedError If the requested backend is not implemented. ValueError If the sampled inputs are invalid.

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

def bayesian_upde_run(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
) -> BayesianUPDEResult:
    """Run UPDE over sampled ``omega`` and ``K_nm`` distributions.

    Parameters
    ----------
    phases : object
        Oscillator phases in radians, shape ``(N,)``.
    omega : object
        Natural-frequency distribution or array.
    knm : object
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    alpha : object
        Phase-lag matrix in radians, shape ``(N, N)``. Use a zero matrix for no
        lag.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    config : BayesianUPDEConfig | None
        Optional configuration object, or ``None`` for defaults.

    Returns
    -------
    BayesianUPDEResult
        The Bayesian UPDE result with uncertainty diagnostics.

    Raises
    ------
    NotImplementedError
        If the requested backend is not implemented.
    ValueError
        If the sampled inputs are invalid.
    """
    resolved = config or BayesianUPDEConfig()
    if resolved.backend != "numpy":
        raise NotImplementedError(
            f"{resolved.backend} Bayesian UPDE backend is not implemented; "
            "use backend='numpy' for reproducible Monte Carlo propagation"
        )

    phase_array = _as_finite_array(phases, name="phases")
    if phase_array.ndim != 1:
        raise ValueError(f"phases must be 1-D, got shape {phase_array.shape}")
    n_oscillators = int(phase_array.shape[0])
    alpha_array = _as_finite_array(alpha, name="alpha")
    if alpha_array.shape != (n_oscillators, n_oscillators):
        raise ValueError(
            f"alpha must have shape {(n_oscillators, n_oscillators)}, "
            f"got {alpha_array.shape}"
        )
    zeta_value = _validate_finite_real(zeta, name="zeta")
    psi_value = _validate_finite_real(psi, name="psi")

    rng = np.random.default_rng(resolved.seed)
    omega_samples = _sample_array(
        omega,
        expected_shape=(n_oscillators,),
        n_samples=resolved.n_samples,
        rng=rng,
        name="omega",
    )
    knm_samples = _sample_array(
        knm,
        expected_shape=(n_oscillators, n_oscillators),
        n_samples=resolved.n_samples,
        rng=rng,
        name="knm",
    )

    final_phase_samples: FloatArray = np.empty(
        (resolved.n_samples, n_oscillators), dtype=np.float64
    )
    r_samples: FloatArray = np.empty(resolved.n_samples, dtype=np.float64)
    psi_samples: FloatArray = np.empty(resolved.n_samples, dtype=np.float64)
    for idx in range(resolved.n_samples):
        final_phases = upde_run(
            phase_array,
            omega_samples[idx],
            knm_samples[idx],
            alpha_array,
            zeta_value,
            psi_value,
            resolved.dt,
            resolved.n_steps,
            method=resolved.method,
            n_substeps=resolved.n_substeps,
            atol=resolved.atol,
            rtol=resolved.rtol,
        )
        final_phase_samples[idx] = final_phases
        order_r, order_psi = compute_order_parameter(final_phases)
        r_samples[idx] = order_r
        psi_samples[idx] = order_psi

    tail = (1.0 - resolved.credible_interval) / 2.0
    lower, upper = np.quantile(r_samples, [tail, 1.0 - tail])
    return BayesianUPDEResult(
        r_samples=np.ascontiguousarray(r_samples, dtype=np.float64),
        final_phase_samples=np.ascontiguousarray(
            final_phase_samples,
            dtype=np.float64,
        ),
        omega_mean=np.mean(omega_samples, axis=0),
        knm_mean=np.mean(knm_samples, axis=0),
        r_mean=float(np.mean(r_samples)),
        r_sigma=float(np.std(r_samples, ddof=1)),
        r_lower=float(lower),
        r_upper=float(upper),
        psi_mean=float(np.mean(psi_samples)),
        sample_count=resolved.n_samples,
        credible_interval=resolved.credible_interval,
        backend=resolved.backend,
        method=resolved.method,
    )

JAX-Accelerated Kuramoto Engine

Optional JAX implementation for GPU-oriented Kuramoto rollouts. It preserves the same validated inputs and phase wrapping semantics as the NumPy engine. Kuramoto and Stuart-Landau state, frequency, growth, coupling, amplitude- coupling, and phase-lag arrays reject boolean, complex, and numeric-string aliases before host conversion or device dispatch. Finite real numeric-object arrays remain supported.

jax_engine

GPU-accelerated Kuramoto solver via JAX JIT compilation.

Raises ImportError if JAX is not installed. Check HAS_JAX before use. Usage: from scpn_phase_orchestrator.upde.jax_engine import HAS_JAX if HAS_JAX: from scpn_phase_orchestrator.upde.jax_engine import JaxUPDEEngine engine = JaxUPDEEngine(n, dt=0.01)

Classes

JaxUPDEEngine

JaxUPDEEngine(
    n: int, dt: float = 0.01, method: str = "rk4"
)

JAX-accelerated Kuramoto/UPDE integrator.

GPU-compiled via jax.jit. First call triggers XLA compilation (~1-3s), subsequent calls run at native speed.

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

def __init__(self, n: int, dt: float = 0.01, method: str = "rk4") -> None:
    if not HAS_JAX:
        msg = "JAX not installed. Install with: pip install jax jaxlib"
        raise ImportError(msg)
    self._n = _validate_positive_int(n, name="n")
    self._dt = _validate_positive_float(dt, name="dt")
    self._method = _validate_method(method)
    euler_fn, rk4_fn = _build_jax_step()
    self._euler = euler_fn
    self._rk4 = rk4_fn

Methods:

step

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

Advance phases by one Kuramoto step on GPU via JIT-compiled JAX.

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 JIT-compiled Kuramoto step.

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
) -> FloatArray:
    """Advance phases by one Kuramoto step on GPU via JIT-compiled JAX.

    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 JIT-compiled Kuramoto step.
    """
    phases = _validate_array(phases, name="phases", shape=(self._n,))
    omegas = _validate_array(omegas, name="omegas", shape=(self._n,))
    knm = _validate_array(knm, name="knm", shape=(self._n, self._n))
    alpha = _validate_array(alpha, name="alpha", shape=(self._n, self._n))
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")

    jp = jnp.asarray(phases)
    jo = jnp.asarray(omegas)
    jk = jnp.asarray(knm)
    ja = jnp.asarray(alpha)

    if self._method == "rk4":
        result = self._rk4(jp, jo, jk, zeta, psi, ja, self._dt)
    else:
        result = self._euler(jp, jo, jk, zeta, psi, ja, self._dt)

    return np.asarray(result)

JaxStuartLandauEngine

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

JAX-accelerated Stuart-Landau integrator (RK4 only).

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

def __init__(self, n: int, dt: float = 0.01) -> None:
    if not HAS_JAX:
        msg = "JAX not installed. Install with: pip install jax jaxlib"
        raise ImportError(msg)
    self._n = _validate_positive_int(n, name="n")
    self._dt = _validate_positive_float(dt, name="dt")
    self._sl_rk4 = _build_jax_sl_step()

Methods:

step

step(
    state: FloatArray,
    omegas: FloatArray,
    mu: FloatArray,
    knm: FloatArray,
    knm_r: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    epsilon: float = 1.0,
) -> FloatArray

Advance Stuart-Landau state by one RK4 step via JIT-compiled JAX.

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

def step(
    self,
    state: FloatArray,
    omegas: FloatArray,
    mu: FloatArray,
    knm: FloatArray,
    knm_r: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    epsilon: float = 1.0,
) -> FloatArray:
    """Advance Stuart-Landau state by one RK4 step via JIT-compiled JAX."""
    state = _validate_array(state, name="state", shape=(2 * self._n,))
    omegas = _validate_array(omegas, name="omegas", shape=(self._n,))
    mu = _validate_array(mu, name="mu", shape=(self._n,))
    knm = _validate_array(knm, name="knm", shape=(self._n, self._n))
    knm_r = _validate_array(knm_r, name="knm_r", shape=(self._n, self._n))
    alpha = _validate_array(alpha, name="alpha", shape=(self._n, self._n))
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    epsilon = _validate_finite_float(epsilon, name="epsilon")

    js = jnp.asarray(state)
    result = self._sl_rk4(
        js,
        jnp.asarray(omegas),
        jnp.asarray(mu),
        jnp.asarray(knm),
        jnp.asarray(knm_r),
        zeta,
        psi,
        jnp.asarray(alpha),
        epsilon,
        self._dt,
    )
    return np.asarray(result)

Stuart-Landau Amplitude Engine

Phase + amplitude dynamics with supercritical/subcritical Hopf bifurcation. State vector: [θ₁...θₙ | r₁...rₙ]. Amplitude coupling via K_r matrix. Amplitudes clamped non-negative after each integration step.

stuart_landau

Thread-safe Stuart-Landau phase-amplitude integrator with backend parity.

StuartLandauEngine advances paired phase and amplitude state vectors using Euler, RK4, or adaptive RK45 methods, with optional Rust acceleration when the kernel is installed. Constructor and step validation reject invalid dimensions, methods, non-finite arrays, and non-finite forcing before solver state changes. The instance lock protects reusable scratch buffers and adaptive timestep state for concurrent callers.

Classes

StuartLandauEngine

StuartLandauEngine(
    n_oscillators: int,
    dt: float,
    method: str = "euler",
    atol: float = 1e-06,
    rtol: float = 0.001,
)

Coupled Stuart-Landau (phase-amplitude) integrator.

State vector layout: state[:n] = phases θ, state[n:] = amplitudes r.

Phase ODE (Acebrón et al. 2005, Rev. Mod. Phys. 77(1)): dθ_i/dt = ω_i + Σ_j K_ij sin(θ_j - θ_i - α_ij) + ζ sin(Ψ - θ_i)

Amplitude ODE

dr_i/dt = (μ_i - r_i²)·r_i + ε Σ_j K^r_ij · r_j · cos(θ_j - θ_i - α_ij)

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

def __init__(
    self,
    n_oscillators: int,
    dt: float,
    method: str = "euler",
    atol: float = 1e-6,
    rtol: float = 1e-3,
):
    n_oscillators = _validate_positive_int(
        n_oscillators,
        name="n_oscillators",
    )
    dt = _validate_positive_float(dt, name="dt")
    atol = _validate_positive_float(atol, name="atol")
    rtol = _validate_positive_float(rtol, name="rtol")
    if method not in ("euler", "rk4", "rk45"):
        msg = f"Unknown method {method!r}, expected 'euler', 'rk4', or 'rk45'"
        raise ValueError(msg)
    self._n = n_oscillators
    self._dt = dt
    self._method = method
    self._atol = atol
    self._rtol = rtol
    self._last_dt = dt

    self._use_rust = False
    try:  # pragma: no cover
        import spo_kernel  # noqa: PLC0415

        self._rust = spo_kernel.PyStuartLandauStepper(
            n_oscillators, dt=dt, method=method, n_substeps=1, atol=atol, rtol=rtol
        )
        self._use_rust = True
    except ImportError:  # pragma: no cover — Rust FFI optional
        pass

    n = n_oscillators
    self._phase_diff = np.empty((n, n), dtype=np.float64)
    self._sin_diff = np.empty((n, n), dtype=np.float64)
    self._cos_diff = np.empty((n, n), dtype=np.float64)
    self._scratch_dtheta = np.empty(n, dtype=np.float64)
    self._scratch_dr = np.empty(n, dtype=np.float64)
    self._scratch_deriv = np.empty(2 * n, dtype=np.float64)

    if method == "rk45":
        self._ks = [np.empty(2 * n, dtype=np.float64) for _ in range(7)]
        self._err_buf = np.empty(2 * n, dtype=np.float64)

    # Serialise concurrent step() callers on this instance so the
    # pre-allocated scratch arrays above are not shared across threads.
    self._lock = threading.RLock()

Attributes

last_dt property

last_dt: float

Last accepted timestep (adapts with RK45).

Returns

float Last accepted timestep (adapts with RK45).

Methods:

step

step(
    state: FloatArray,
    omegas: FloatArray,
    mu: FloatArray,
    knm: FloatArray,
    knm_r: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    epsilon: float = 1.0,
) -> FloatArray

Advance (θ, r) by one timestep. Returns new state (2N,).

Parameters

state : FloatArray Finite real numeric Stuart-Landau state [θ; r], shape (2N,). Boolean, complex, and numeric-string aliases are rejected. omegas : FloatArray Natural frequencies in rad/s, shape (N,). mu : FloatArray Per-oscillator linear growth parameters μ, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). knm_r : FloatArray Amplitude coupling matrix, 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. epsilon : float Finite-difference perturbation size.

Returns

FloatArray The finite real numeric [θ; r] state, shape (2N,).

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

def step(
    self,
    state: FloatArray,
    omegas: FloatArray,
    mu: FloatArray,
    knm: FloatArray,
    knm_r: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    epsilon: float = 1.0,
) -> FloatArray:
    """Advance (θ, r) by one timestep. Returns new state (2N,).

    Parameters
    ----------
    state : FloatArray
        Finite real numeric Stuart-Landau state ``[θ; r]``, shape
        ``(2N,)``. Boolean, complex, and numeric-string aliases are
        rejected.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    mu : FloatArray
        Per-oscillator linear growth parameters ``μ``, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    knm_r : FloatArray
        Amplitude coupling matrix, 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.
    epsilon : float
        Finite-difference perturbation size.

    Returns
    -------
    FloatArray
        The finite real numeric ``[θ; r]`` state, shape ``(2N,)``.
    """
    (
        state,
        omegas,
        mu,
        knm,
        knm_r,
        zeta,
        psi,
        alpha,
        epsilon,
    ) = self._validate(state, omegas, mu, knm, knm_r, zeta, psi, alpha, epsilon)
    with self._lock:
        if self._use_rust:  # pragma: no cover
            result = _validate_state_array(
                self._rust.step(
                    state,
                    omegas,
                    mu,
                    knm.ravel(),
                    knm_r.ravel(),
                    zeta,
                    psi,
                    alpha.ravel(),
                    epsilon,
                ),
                name="Rust state output",
                shape=(2 * self._n,),
                finite_message="Rust state output contains NaN or Inf",
            )
            self._last_dt = _validate_positive_float(
                self._rust.last_dt,
                name="Rust last_dt",
            )
            return result
        p: _Params = (omegas, mu, knm, knm_r, zeta, psi, alpha, epsilon)
        if self._method == "euler":
            return self._euler_step(state, p)
        if self._method == "rk45":
            return self._rk45_step(state, p)
        return self._rk4_step(state, p)

compute_order_parameter

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

Amplitude-weighted Kuramoto: Z = mean(r_i · exp(i·θ_i)).

Parameters

state : FloatArray Finite real numeric Stuart-Landau state [θ; r], shape (2N,). Coercive aliases are rejected.

Returns

tuple[float, float] The amplitude-weighted (R, ψ) order parameter and phase.

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

def compute_order_parameter(self, state: FloatArray) -> tuple[float, float]:
    """Amplitude-weighted Kuramoto: Z = mean(r_i · exp(i·θ_i)).

    Parameters
    ----------
    state : FloatArray
        Finite real numeric Stuart-Landau state ``[θ; r]``, shape
        ``(2N,)``. Coercive aliases are rejected.

    Returns
    -------
    tuple[float, float]
        The amplitude-weighted ``(R, ψ)`` order parameter and phase.
    """
    n = self._n
    state = _validate_state_array(
        state,
        name="state",
        shape=(2 * n,),
        finite_message="state contains NaN or Inf",
    )
    z = np.mean(state[n:] * np.exp(1j * state[:n]))
    return float(np.abs(z)), float(np.angle(z) % TWO_PI)

compute_mean_amplitude

compute_mean_amplitude(state: FloatArray) -> float

Mean amplitude across all oscillators.

Parameters

state : FloatArray Finite real numeric Stuart-Landau state [θ; r], shape (2N,). Coercive aliases are rejected.

Returns

float The mean amplitude across all oscillators.

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

def compute_mean_amplitude(self, state: FloatArray) -> float:
    """Mean amplitude across all oscillators.

    Parameters
    ----------
    state : FloatArray
        Finite real numeric Stuart-Landau state ``[θ; r]``, shape
        ``(2N,)``. Coercive aliases are rejected.

    Returns
    -------
    float
        The mean amplitude across all oscillators.
    """
    state = _validate_state_array(
        state,
        name="state",
        shape=(2 * self._n,),
        finite_message="state contains NaN or Inf",
    )
    return float(np.mean(state[self._n :]))

Simplicial (3-Body) Engine

Higher-order interactions beyond pairwise coupling. The 3-body term σ₂/N² Σ_{j,k} sin(θ_j + θ_k - 2θ_i) produces explosive (first-order) synchronization transitions not achievable with pairwise coupling alone. Vectorized via trig identity: 2·S_i·C_i where S = Σsin(Δθ), C = Σcos(Δθ).

Use this engine when the physical or learned topology contains group effects that cannot be decomposed into independent pairwise edges. Practical examples include neural assemblies with co-active triplets, multi-agent coordination under triangular constraints, reaction loops, power-network group modes, and simplicial-complex or hypergraph-derived topology where synchronization thresholds depend on 3-body motifs.

Direct Go, Julia, and Mojo simplicial accelerator entrypoints share the same validated torus boundary before optional runtime loading: phase and frequency vectors must be finite real one-dimensional float64 arrays matching the oscillator count; flattened pairwise coupling and phase-lag buffers must have exactly N*N values; pairwise self-coupling K_ii must be zero because the pairwise graph represents interactions between distinct oscillators; zeta, psi, sigma2, dt, and n_steps must be finite non-boolean controls with non-negative triadic strength, positive timestep, and non-negative step count. Zero-step direct calls return a copy of the input phases without loading the optional runtime. Direct input arrays, shared/public backend outputs, and Julia raw returns reject numeric-string aliases before float coercion. Backend outputs must be finite torus phases in [0, 2*pi). The public dispatcher and Rust wrapper apply that same output contract to optional backend returns before exposing SimplicialEngine.run() results, so backend physics-contract faults raise instead of falling through as trusted higher-order synchronization evidence.

Gambuzza et al. 2023, Nature Physics; Tang et al. 2025. Detailed documentation: Simplicial (3-body) — detailed reference

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:

Second-Order Inertial Engine (Power Grids)

Swing equation: m_i θ̈_i + d_i θ̇_i = P_i + Σ_j K_ij sin(θ_j - θ_i). Models power grid transient stability where m_i is generator inertia, d_i is damping, P_i is power injection (positive = generation, negative = load), and K_ij is transmission line susceptance. RK4 integration.

Includes frequency_deviation() (Hz from nominal — >0.5 Hz triggers load shedding in real grids) and coherence() (phase-lock measure).

Filatrella et al. 2008; Dörfler & Bullo 2014.

Optional inertial backend outputs are validated before public return: theta and omega_dot must keep oscillator cardinality and finite values, with returned phases inside [0, 2*pi). Backend loader/runtime unavailability may fall through to Python; malformed backend physics payloads raise instead of becoming swing-equation state. Public and direct inertial state vectors, coupling buffers, scalar controls, step counts, frequency_deviation() inputs, coherence() inputs, optional backend outputs, and direct Julia raw returns reject numeric-string aliases before Python, NumPy, or accelerator coercion.

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 Regime Detection

Extracts instantaneous phase from price/return time series via Hilbert transform, computes Kuramoto order parameter R(t) across assets, classifies synchronization regimes (desync/transition/synchronised), and detects crash early warning signals (R crossing threshold from below).

Direct Go, Julia, and Mojo market accelerator entrypoints share a validated float64 boundary before optional runtime loading: flattened phase payloads must be finite real vectors with exactly T*N values; T, N, and PLV window controls must be positive non-boolean integers; the PLV window must not exceed T; backend R(t) outputs must have length T and lie in [0, 1]; rolling PLV outputs must have the expected (T-window+1)*N*N cardinality, lie in [0, 1], preserve unit diagonals, and remain symmetric. The public market dispatcher applies the same output contract to optional backend returns before exposing market_order_parameter() or market_plv() results, so backend physics-contract faults propagate instead of falling through as trusted market evidence.

R(t) → 1 preceding market crashes documented for Black Monday 1987 and the 2008 financial crisis (arXiv:1109.1167).

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 with both spatial position x_i ∈ R^D and oscillator phase θ_i ∈ S¹. Phase modulates spatial attraction (J parameter); spatial proximity modulates phase coupling (K/|x_ij|). D-dimensional (2D, 3D supported).

Five collective states: static sync, static async, static phase wave, splintered phase wave, active phase wave — depending on J and K signs.

O'Keeffe, Hong, Strogatz, Nature Communications 2017.

Direct Go, Julia, and Mojo swarmalator accelerator entrypoints share a validated position-phase boundary before optional runtime loading: positions must be finite real float64 values with shape (N, D) or exactly N*D flattened values without numeric-string aliases; phase and frequency vectors must be finite real one-dimensional float64 arrays of length N without numeric-string aliases; N, D, and dt must be positive and non-string typed; and attraction, repulsion, phase-attraction modulation, and phase-coupling coefficients must be finite real controls. Backend outputs must return finite positions and torus phases in [0, 2*pi) without numeric-string aliases, and Mojo stdout must contain exactly N*D + N scalar lines. The public SwarmalatorEngine.step() dispatcher and Rust wrapper apply the same output contract before publication, including object-dtype boolean-alias rejection, while public constructor controls, state arrays, scalar controls, step counts, order-parameter inputs, optional backend outputs, and direct Julia raw returns reject numeric-string aliases before Python, NumPy, or accelerator coercion. Loader and runtime unavailability still fall back to Python.

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:

Stochastic Engine

Euler-Maruyama integration with Gaussian noise injection. Includes automatic optimal noise tuning: D ≈ K·R_det/2 (Tselios et al. 2025). Counter-intuitive: noise at D INCREASES synchronization (stochastic resonance). Self-consistency solved via modified Bessel equation (Acebrón et al. 2005).

Public injection validates phase arrays even when D=0: inputs must be finite, one-dimensional, real numeric payloads, with boolean, complex, and numeric-string aliases rejected before the no-op return or noise arithmetic. Noise sweeps apply the same source-type contract to D_range, validate non-negative integer seeds before range arithmetic, and publish only physical NoiseProfile records (D >= 0, both order parameters in [0, 1]).

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

Symplectic Euler on T^N using SO(2) exponential map: z_i = exp(iθ_i). Avoids mod 2π discontinuity errors that accumulate in standard integrators over long simulations. Essential for multi-hour or multi-day simulations where phase wrapping drift becomes significant. The public dispatcher validates optional backend outputs before publication: selected Rust, Go, Julia, or Mojo returns must be finite phase vectors with the same oscillator cardinality, values in [0, 2*pi), and no numeric-string aliases. Public constructor, state, scalar-control, and order-parameter phase inputs plus direct Go/Julia/Mojo phase, frequency, coupling, phase-lag, scalar, count, and backend-output boundaries reject numeric-string aliases before float coercion or optional native runtime loading. Detailed documentation: Geometric (SO(2)) — detailed reference

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:

Time-Delayed Coupling Engine

The circular buffer supports one fixed discrete delay for every coupling pair, with current-phase history used while the buffer fills. Time delays generate "effective higher-order interactions for free" (Ciszak et al. 2025) because the delayed coupling mixes information across multiple timescales. Public and direct phase, frequency, coupling, and phase-lag arrays reject numeric-string aliases before float coercion. The same pre-coercion contract guards optional-backend outputs and direct Julia raw returns before phase cardinality, finiteness, and [0, 2*pi) validation.

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:

Ott-Antonsen Mean-Field Reduction

Exact analytical reduction for globally-coupled Kuramoto with Lorentzian frequency distribution. Reduces N-oscillator system to a single complex ODE: 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) computation vs O(N) for full simulation).

Direct Go, Julia, and Mojo Ott-Antonsen accelerator entrypoints share the same scalar boundary before optional runtime loading: the complex order parameter must lie inside the OA unit disk, Lorentzian width must be non-negative, timestep and step count must be positive, and all scalar controls must be finite non-boolean real values. Backend outputs are accepted only when the returned complex state remains in the OA unit disk, R matches |z|, and psi matches atan2(Im(z), Re(z)), preserving the physical mean-field state contract across the polyglot chain. The public dispatcher applies the same output contract to optional backend returns before publishing OAState, so inconsistent R, inconsistent psi, or boolean-alias scalar outputs fail closed instead of becoming mean-field evidence. Empirical frequency samples reject boolean, complex, and numeric-string aliases before Lorentzian fitting. The direct Rust runner uses the same pre-coercion scalar-input contract as Go, Julia, and Mojo, while its optional steady-state helper must return a finite real order parameter in [0, 1] before publication.

Detailed documentation: Ott-Antonsen Reduction — detailed reference

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:

Variational Free Energy Predictor

Implementation of Friston's Free Energy Principle mapped to Kuramoto dynamics. Precision-weighted prediction error drives coupling updates; KL divergence provides a complexity penalty. Online precision estimation from error variance.

Includes PredictionModel (forward prediction with error injection) and VariationalPredictor (FEP-Kuramoto correspondence). Their public phase, frequency, predicted-state, observed-state, and precision vectors reject boolean, complex, and numeric-string aliases before float64 conversion. Real numeric object arrays remain supported, and malformed array protocol or conversion payloads are normalized to field-specific ValueError failures before predictor state can mutate.

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

Adjoint Gradient Computation

Finite-difference and JAX-autodiff gradients of the synchronization cost (1 - R) with respect to the coupling matrix K_nm. Used for gradient-based coupling optimization without the overhead of forward-mode differentiation.

Both paths validate a non-empty finite real phase vector, matching frequency vector, square coupling and phase-lag matrices, and a zero self-coupling diagonal before simulation or optional-backend import. Step counts must be positive non-boolean integers. The finite-difference perturbation and JAX timestep must be positive finite reals; finite-difference drive scalars must be finite real values. Boolean, complex, and numeric-string array aliases fail closed instead of being coerced by NumPy or JAX.

adjoint

Adjoint and finite-difference sensitivities for UPDE coupling gradients.

The module defines the synchronization cost 1 - R and two gradient paths: a deterministic NumPy finite-difference estimator over coupling entries and a diffrax continuous-adjoint implementation when the optional JAX/diffrax stack is installed. The finite-difference path mutates only local coupling copies for each perturbation; the continuous-adjoint path integrates the same Kuramoto-Sakaguchi field and differentiates through the solver, agreeing with the finite-difference reference in direction and — in the fine-dt limit — in magnitude. It fails explicitly with ImportError when the dependency is absent instead of silently claiming accelerated gradients, and never mutates the process-global jax_enable_x64 flag.

Classes

Functions:

cost_R

cost_R(phases: FloatArray) -> float

Cost 1 − R (minimise to maximise synchronisation).

Parameters

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

Returns

float The cost 1 − R for the supplied phases.

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

def cost_R(phases: FloatArray) -> float:
    """Cost ``1 − R`` (minimise to maximise synchronisation).

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

    Returns
    -------
    float
        The cost ``1 − R`` for the supplied phases.
    """
    R, _ = compute_order_parameter(phases)
    return 1.0 - R

gradient_knm_fd

gradient_knm_fd(
    engine: UPDEEngine,
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray,
    n_steps: int = 100,
    epsilon: float = 0.0001,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray

Finite-difference gradient of cost_R w.r.t. knm entries.

For each K_ij (i≠j), perturbs by ±ε and measures the effect on R after n_steps. Returns gradient matrix ∂(1-R)/∂K_ij.

Complexity: O(N² · N² · n_steps) = O(N⁴ · n_steps) because each of the ~N² off-diagonal entries requires a full N-step simulation that is itself O(N²) per step. Use gradient_knm_jax() for anything beyond N≈16. The adjoint method via diffrax reduces this to O(n_steps) but requires JAX.

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). Use a zero matrix for no lag. n_steps : int Positive number of integration steps to run. epsilon : float Strictly positive finite-difference perturbation size. zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians.

Returns

FloatArray The finite-difference gradient of the cost with respect to knm.

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

def gradient_knm_fd(
    engine: UPDEEngine,
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray,
    n_steps: int = 100,
    epsilon: float = 1e-4,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray:
    """Finite-difference gradient of cost_R w.r.t. knm entries.

    For each K_ij (i≠j), perturbs by ±ε and measures the effect on R
    after n_steps. Returns gradient matrix ∂(1-R)/∂K_ij.

    Complexity: O(N² · N² · n_steps) = O(N⁴ · n_steps) because each of
    the ~N² off-diagonal entries requires a full N-step simulation that is
    itself O(N²) per step. Use gradient_knm_jax() for anything beyond N≈16.
    The adjoint method via diffrax reduces this to O(n_steps) but requires JAX.

    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)``. Use a zero matrix for no
        lag.
    n_steps : int
        Positive number of integration steps to run.
    epsilon : float
        Strictly positive finite-difference perturbation size.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.

    Returns
    -------
    FloatArray
        The finite-difference gradient of the cost with respect to ``knm``.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    epsilon = _validate_positive_float(epsilon, name="epsilon")
    zeta = _validate_finite_real(zeta, name="zeta")
    psi = _validate_finite_real(psi, name="psi")
    phases_init, omegas, knm, alpha = _validate_adjoint_arrays(
        phases_init,
        omegas,
        knm,
        alpha,
    )
    n = knm.shape[0]
    grad = np.zeros((n, n), dtype=np.float64)

    for i in range(n):
        for j in range(n):
            if i == j:
                continue

            knm_plus = knm.copy()
            knm_plus[i, j] += epsilon
            p_plus = engine.run(
                phases_init, omegas, knm_plus, zeta, psi, alpha, n_steps
            )
            c_plus = cost_R(p_plus)

            knm_minus = knm.copy()
            knm_minus[i, j] -= epsilon
            p_minus = engine.run(
                phases_init, omegas, knm_minus, zeta, psi, alpha, n_steps
            )
            c_minus = cost_R(p_minus)

            grad[i, j] = (c_plus - c_minus) / (2 * epsilon)

    return grad

gradient_knm_jax

gradient_knm_jax(
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray,
    n_steps: int = 100,
    dt: float = 0.01,
) -> FloatArray

Gradient of cost_R w.r.t. knm via a diffrax continuous adjoint.

Integrates the Kuramoto-Sakaguchi field

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

over [0, n_steps·dt] with an adaptive solver (diffrax.Tsit5) and a RecursiveCheckpointAdjoint, then differentiates cost_R of the final phases with respect to knm using reverse-mode autodiff. This is the O(1)-memory continuous-adjoint path the finite-difference estimator in :func:gradient_knm_fd approximates; the two agree in direction and, in the fine-dt limit, in magnitude (the finite-difference reference differentiates the discrete explicit-Euler map, so a fixed dt leaves an O(dt) discretisation gap).

The solver runs in JAX's active default precision — it does not mutate the process-global jax_enable_x64 flag, so it does not perturb the float32 default the rest of the differentiable stack relies on. Callers that need float64 gradients must enable x64 at process start-up themselves.

Parameters

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). Use a zero matrix for no lag; this is the regime in which the gradient matches the drive-free (ζ = ψ = 0) finite-difference reference. n_steps : int Positive number of nominal steps; the integration horizon is n_steps · dt. dt : float Positive finite nominal step size setting the integration horizon and the adaptive solver's initial step.

Returns

FloatArray The continuous-adjoint gradient of the cost with respect to knm.

Raises

ImportError If the optional JAX/diffrax stack is not installed.

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

def gradient_knm_jax(
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray,
    n_steps: int = 100,
    dt: float = 0.01,
) -> FloatArray:
    """Gradient of ``cost_R`` w.r.t. ``knm`` via a diffrax continuous adjoint.

    Integrates the Kuramoto-Sakaguchi field

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

    over ``[0, n_steps·dt]`` with an adaptive solver (``diffrax.Tsit5``) and a
    ``RecursiveCheckpointAdjoint``, then differentiates ``cost_R`` of the final
    phases with respect to ``knm`` using reverse-mode autodiff. This is the
    ``O(1)``-memory continuous-adjoint path the finite-difference estimator in
    :func:`gradient_knm_fd` approximates; the two agree in direction and, in the
    fine-``dt`` limit, in magnitude (the finite-difference reference
    differentiates the discrete explicit-Euler map, so a fixed ``dt`` leaves an
    ``O(dt)`` discretisation gap).

    The solver runs in JAX's active default precision — it does **not** mutate
    the process-global ``jax_enable_x64`` flag, so it does not perturb the
    float32 default the rest of the differentiable stack relies on. Callers that
    need float64 gradients must enable x64 at process start-up themselves.

    Parameters
    ----------
    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)``. Use a zero matrix for no
        lag; this is the regime in which the gradient matches the drive-free
        (``ζ = ψ = 0``) finite-difference reference.
    n_steps : int
        Positive number of nominal steps; the integration horizon is
        ``n_steps · dt``.
    dt : float
        Positive finite nominal step size setting the integration horizon and
        the adaptive solver's initial step.

    Returns
    -------
    FloatArray
        The continuous-adjoint gradient of the cost with respect to ``knm``.

    Raises
    ------
    ImportError
        If the optional JAX/diffrax stack is not installed.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    dt = _validate_positive_float(dt, name="dt")
    phases_init, omegas, knm, alpha = _validate_adjoint_arrays(
        phases_init,
        omegas,
        knm,
        alpha,
    )
    import diffrax
    import jax
    import jax.numpy as jnp

    theta0 = jnp.asarray(phases_init)
    om = jnp.asarray(omegas)
    al = jnp.asarray(alpha)
    t1 = float(n_steps) * dt

    def _field(_t: Any, theta: Any, coupling: Any) -> Any:
        """Kuramoto-Sakaguchi tangent-space derivative dθ/dt."""
        diff = theta[jnp.newaxis, :] - theta[:, jnp.newaxis] - al
        return om + jnp.sum(coupling * jnp.sin(diff), axis=1)

    def _cost(knm_j: Any) -> Any:
        """Continuous-adjoint roll-out to the sync cost ``1 − R``."""
        solution = diffrax.diffeqsolve(
            diffrax.ODETerm(_field),
            diffrax.Tsit5(),
            t0=0.0,
            t1=t1,
            dt0=dt,
            y0=theta0,
            args=knm_j,
            stepsize_controller=diffrax.PIDController(rtol=1e-6, atol=1e-6),
            adjoint=diffrax.RecursiveCheckpointAdjoint(),
            max_steps=16384,
        )
        theta_final = solution.ys[-1]
        z = jnp.mean(jnp.exp(1j * theta_final))
        return 1.0 - jnp.abs(z)

    grad = jax.grad(_cost)(jnp.asarray(knm))
    return np.asarray(grad, dtype=np.float64)

Order Parameters & Metrics

Kuramoto order parameter R (global coherence), PLV (pairwise phase-locking value), and layer coherence (R for oscillator subsets). Optional Rust acceleration.

Direct Go, Julia, and Mojo order-parameter entrypoints share the same typed boundary before optional runtime loading: phase payloads must be one-dimensional finite real vectors, PLV inputs must be equal-length phase vectors, and layer-coherence indices must be unique in-range oscillator indices. Empty zero-measure calls return the neutral value without loading an optional runtime: (0.0, 0.0) for global order, 0.0 for PLV, and 0.0 for layer coherence. Backend outputs are accepted only when physical: R, PLV, and layer coherence must be finite values in [0, 1], and mean phase must be finite before being canonicalised to the public [0, 2*pi) convention. Public phase-vector inputs and shared backend scalar outputs reject numeric-string aliases before float coercion, so text cannot acquire numeric provenance at either publication boundary. Mojo stdout remains an explicit text protocol and is parsed before the shared typed scalar validator runs.

The release benchmark gate for this surface is:

python benchmarks/order_params_benchmark.py --parity-gate --sizes 64 --calls 1

It records Rust/Mojo/Julia/Go/Python status, timing, unavailable-toolchain reasons, deterministic hashes, and tolerance-bounded parity against the forced Python reference for global R, mean phase, PLV, and layer coherence. The reference-suite snapshot exposes this gate as order_parameter_polyglot. The public dispatcher applies the same scalar output contract to optional backend returns before publication: order-parameter magnitudes, PLV, and layer coherence must be finite real values in [0, 1], with boolean aliases rejected instead of widened to synthetic 0.0 or 1.0 evidence and numeric strings rejected instead of converted into typed backend evidence.

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",
    )

metrics

Immutable diagnostic dataclasses for UPDE lock and layer state snapshots.

The module contains data-only records for pairwise lock signatures, per-layer coherence and amplitude summaries, and aggregate UPDE diagnostic state. These types intentionally perform no simulation, I/O, or mutation beyond dataclass construction; producers are responsible for numeric validation before emitting snapshots.

Classes

LockSignature dataclass

LockSignature(
    source_layer: int,
    target_layer: int,
    plv: float,
    mean_lag: float,
)

Phase-locking value and mean lag between two layers.

LayerState dataclass

LayerState(
    R: float,
    psi: float,
    lock_signatures: dict[str, LockSignature] = dict(),
    mean_amplitude: float = 0.0,
    amplitude_spread: float = 0.0,
)

Per-layer diagnostics: order parameter R, mean phase, amplitude stats.

UPDEState dataclass

UPDEState(
    layers: list[LayerState],
    cross_layer_alignment: FloatArray,
    stability_proxy: float,
    regime_id: str,
    mean_amplitude: float = 0.0,
    pac_max: float = 0.0,
    subcritical_fraction: float = 0.0,
    boundary_violation_count: int = 0,
    imprint_mean: float = 0.0,
)

Full UPDE diagnostic snapshot: per-layer states and aggregates.

Phase-Amplitude Coupling (PAC)

Modulation index (MI) via Tort et al. 2010. Bins low-frequency phase, computes mean high-frequency amplitude per bin, KL divergence from uniform. Produces N×N PAC matrix: entry [i,j] = MI(phase_i, amplitude_j). Direct Go, Julia, and Mojo PAC entrypoints now share the same typed boundary before optional runtime loading: phase and amplitude payloads must be finite real float64 vectors, amplitudes must be non-negative, n_bins must be a non-boolean integer of at least two, and matrix calls require exactly T*N flattened phase and amplitude samples. Empty common MI windows return zero without loading optional runtimes. Backend MI and pairwise PAC outputs are accepted only when finite, correctly sized, and inside the physical [0, 1] interval. Direct Mojo PAC output must contain exactly one scalar line for modulation-index calls or exactly N×N scalar lines for matrix calls; blank, non-finite, truncated, or overlong text output is rejected before public assembly. The public PAC dispatcher applies the same output contract to optional backend returns before exposing modulation_index() or pac_matrix() results, so backend physics-contract faults fail closed instead of being clipped into synthetic PAC evidence. Central to neuroscience cross-frequency coupling analysis.

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

Envelope & Numerics

Amplitude envelope extraction and numerical integration utilities (DP54 coefficients, error estimation, step size control). The public envelope dispatcher validates optional backend RMS-envelope outputs before publication: extracted envelopes must keep input cardinality, remain finite, stay non-negative, and reject numeric-string aliases before coercion; modulation-depth outputs must be finite scalars inside [0, 1] and reject numeric-string aliases as well. Public amplitude/envelope inputs and window share the same alias boundary. Loader/runtime failures still fall through to the Python floor. Detailed documentation: Envelope (RMS) — detailed reference

envelope

Sliding-window RMS envelope and modulation-depth statistic.

Exposes a 5-backend fallback chain. AVAILABLE_BACKENDS keeps the canonical fallback order; ACTIVE_BACKEND is chosen by a small hot-path probe so slow external wrappers do not displace the faster local path.

The sliding-window RMS uses the O(T) cumulative-sum form: compute cs[i] = Σ_{k < i} x_k², then rms[i] = sqrt((cs[i+w] − cs[i]) / w) for valid indices, with a front-pad of the first valid value. The 1-D path is on the 5-backend chain; the 2-D (T, N) batched path stays pure NumPy because the Rust FFI is 1-D-only and the vectorised NumPy form is already near-optimal at realistic N.

Classes

EnvelopeState dataclass

EnvelopeState(
    mean_amplitude: float,
    amplitude_spread: float,
    modulation_depth: float,
    subcritical_count: int,
)

Snapshot of amplitude envelope statistics.

Functions:

extract_envelope

extract_envelope(
    amplitudes_history: FloatArray, window: int = 10
) -> FloatArray

Sliding-window RMS envelope.

Parameters

amplitudes_history : FloatArray (T,) or (T, N) amplitude time series. window : int RMS window length in samples.

Returns

FloatArray Same shape as input; the first window − 1 entries are front-padded with the first valid RMS value.

Raises

ValueError If window is not a positive integer no larger than the history.

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

def extract_envelope(
    amplitudes_history: FloatArray,
    window: int = 10,
) -> FloatArray:
    """Sliding-window RMS envelope.

    Parameters
    ----------
    amplitudes_history : FloatArray
        ``(T,)`` or ``(T, N)`` amplitude time series.
    window : int
        RMS window length in samples.

    Returns
    -------
    FloatArray
        Same shape as input; the first ``window − 1`` entries are front-padded with the
        first valid RMS value.

    Raises
    ------
    ValueError
        If ``window`` is not a positive integer no larger than the history.
    """
    if _contains_numeric_string_alias(amplitudes_history):
        raise ValueError("amplitudes_history must not contain numeric-string aliases")
    if _is_numeric_string_alias(window):
        raise ValueError("window must not be a numeric-string alias")
    amplitudes = np.asarray(amplitudes_history, dtype=np.float64)
    if amplitudes.size == 0:
        return amplitudes.copy()
    if window < 1:
        raise ValueError(f"window must be >= 1, got {window}")
    if not np.all(np.isfinite(amplitudes)):
        raise ValueError("amplitudes_history must contain only finite values")

    if amplitudes.ndim == 1:
        if window >= amplitudes.size:
            return _extract_1d_python(amplitudes, int(window))
        backend_fn = _dispatch("extract")
        if backend_fn is not None:
            fn = cast("Callable[[FloatArray, int], FloatArray]", backend_fn)
            return validate_extract_envelope_output(
                fn(amplitudes, int(window)),
                n=int(amplitudes.size),
            )
        return _extract_1d_python(amplitudes, int(window))

    if amplitudes.ndim != 2:
        raise ValueError("amplitudes_history must be 1-D or 2-D")

    # 2-D path stays pure NumPy.
    sq = amplitudes**2
    if window >= sq.shape[0]:
        rms = np.sqrt(np.mean(sq, axis=0))
        return np.tile(rms, (sq.shape[0], 1))
    cs = np.cumsum(sq, axis=0)
    cs = np.vstack([np.zeros((1, sq.shape[1]), dtype=np.float64), cs])
    rms = np.sqrt((cs[window:] - cs[:-window]) / window)
    first = rms[0] if rms.shape[0] > 0 else np.zeros(sq.shape[1])
    return np.vstack([np.tile(first, (window - 1, 1)), rms])

envelope_modulation_depth

envelope_modulation_depth(envelope: FloatArray) -> float

Modulation depth (max − min) / (max + min) ∈ [0, 1].

Returns 0.0 for empty or non-positive envelopes.

Parameters

envelope : FloatArray An amplitude-envelope time series, shape (T,).

Returns

float The modulation depth (max − min) / (max + min) in [0, 1].

Raises

ValueError If envelope contains numeric-string aliases.

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

def envelope_modulation_depth(envelope: FloatArray) -> float:
    """Modulation depth ``(max − min) / (max + min) ∈ [0, 1]``.

    Returns ``0.0`` for empty or non-positive envelopes.

    Parameters
    ----------
    envelope : FloatArray
        An amplitude-envelope time series, shape ``(T,)``.

    Returns
    -------
    float
        The modulation depth ``(max − min) / (max + min)`` in ``[0, 1]``.

    Raises
    ------
    ValueError
        If ``envelope`` contains numeric-string aliases.
    """
    if _contains_numeric_string_alias(envelope):
        raise ValueError("envelope must not contain numeric-string aliases")
    if envelope.size == 0:
        return 0.0
    backend_fn = _dispatch("mod")
    if backend_fn is not None:
        fn = cast("Callable[[FloatArray], float]", backend_fn)
        return validate_envelope_modulation_output(fn(envelope))
    flat = envelope.ravel()
    vmax = float(np.max(flat))
    vmin = float(np.min(flat))
    denom = vmax + vmin
    if denom <= 0.0:
        return 0.0
    return float((vmax - vmin) / denom)

numerics

Numerical integration configuration and explicit-step stability checks.

IntegrationConfig records solver tolerances and method selection, while check_stability provides a CFL-like phase-step bound for explicit Kuramoto integration. The helper is deliberately conservative and side-effect free: it does not adapt solvers or clamp parameters, it only reports whether the supplied derivative bound keeps a single step below a half-cycle.

Classes

IntegrationConfig dataclass

IntegrationConfig(
    dt: float,
    substeps: int = 1,
    method: str = "euler",
    max_dt: float = 0.01,
    atol: float = 1e-06,
    rtol: float = 0.001,
)

Numerical integration parameters for the phase ODE solver.

Functions:

check_stability

check_stability(
    dt: float, max_omega: float, max_coupling: float
) -> bool

CFL-like stability bound for explicit Kuramoto integration.

Analogous to Courant–Friedrichs–Lewy (1928); see docs/specs/upde_numerics.md. dt * max_deriv < pi ensures phase change stays below half-cycle per step.

Parameters

dt : float Integration step size. max_omega : float Largest absolute natural frequency in the system. max_coupling : float Largest absolute coupling magnitude in the system.

Returns

bool True when the timestep satisfies the CFL-like stability bound.

Raises

ValueError If dt, max_omega, or max_coupling is not finite and positive.

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

def check_stability(dt: float, max_omega: float, max_coupling: float) -> bool:
    """CFL-like stability bound for explicit Kuramoto integration.

    Analogous to Courant–Friedrichs–Lewy (1928); see docs/specs/upde_numerics.md.
    dt * max_deriv < pi ensures phase change stays below half-cycle per step.

    Parameters
    ----------
    dt : float
        Integration step size.
    max_omega : float
        Largest absolute natural frequency in the system.
    max_coupling : float
        Largest absolute coupling magnitude in the system.

    Returns
    -------
    bool
        ``True`` when the timestep satisfies the CFL-like stability bound.

    Raises
    ------
    ValueError
        If ``dt``, ``max_omega``, or ``max_coupling`` is not finite and positive.
    """
    if (
        type(dt) in {float, int}
        and type(max_omega) in {float, int}
        and type(max_coupling) in {float, int}
    ):
        if not math.isfinite(dt) or dt <= 0.0:
            raise ValueError("dt must be a finite positive real")
        if not math.isfinite(max_omega) or max_omega < 0.0:
            raise ValueError("max_omega must be a finite non-negative real")
        if not math.isfinite(max_coupling) or max_coupling < 0.0:
            raise ValueError("max_coupling must be a finite non-negative real")
        max_deriv = max_omega + max_coupling
        if max_deriv == 0.0:
            return True
        return dt * max_deriv < math.pi
    dt_value = _validate_positive_finite(dt, name="dt")
    omega_bound = _validate_non_negative_finite(max_omega, name="max_omega")
    coupling_bound = _validate_non_negative_finite(max_coupling, name="max_coupling")
    max_deriv = omega_bound + coupling_bound
    if max_deriv == 0.0:
        return True
    # pi threshold: phase change per step must stay below half-cycle
    return dt_value * max_deriv < math.pi

Splitting Engine

Operator-splitting UPDE integrator for stiff regimes and deterministic phase update decomposition.

splitting

Strang second-order operator splitting for the Kuramoto ODE.

Exposes a 5-backend fallback chain.

Scheme

Split dθ/dt = ω + Σ_j K_ij · sin(θ_j − θ_i − α_ij) + ζ · sin(ψ − θ_i) into

A: dθ/dt = ω              (exact rotation)
B: dθ/dt = coupling       (RK4 on the nonlinear part)

and compose symmetrically as A(dt/2) → B(dt) → A(dt/2) (Strang scheme, second-order in dt).

Why split?

The ω flow is linear, so it has no truncation error; folding it into a monolithic RK45 burns integrator budget on a solvable direction while damping unrelated accuracy in the nonlinear direction. Reference: Hairer, Lubich & Wanner 2006, Geometric Numerical Integration §II.5.

Numerics

The B-stage RK4 uses the Rust kernel's sin(θ_j − θ_i) = sin(θ_j)·cos(θ_i) − cos(θ_j)·sin(θ_i) expansion on the alpha-zero branch so that floating-point rounding matches Rust (spo-engine/src/splitting.rs) bit-for-bit. Nonzero alpha falls back to the direct sin(diff) form in all five backends.

Classes

SplittingEngine

SplittingEngine(n_oscillators: int, dt: float)

Strang-split Kuramoto stepper with 5-backend dispatch.

The engine's geometry is (n, dt); the step is stateless.

Create a Strang-splitting engine for n_oscillators and dt.

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

def __init__(self, n_oscillators: int, dt: float):
    """Create a Strang-splitting engine for ``n_oscillators`` and ``dt``."""
    n_oscillators = _validate_positive_int(
        n_oscillators,
        name="n_oscillators",
    )
    dt = _validate_nonzero_finite_float(dt, name="dt")
    # Negative dt is intentional for symplectic-reversibility
    # checks — Strang is time-reversible, so stepping forward
    # then with dt → −dt returns to the starting state.
    self._n = n_oscillators
    self._dt = dt

Methods:

step

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

One Strang-split step: A(dt/2) → B(dt) → A(dt/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 Strang-split step.

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
) -> FloatArray:
    """One Strang-split step: A(dt/2) → B(dt) → A(dt/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 Strang-split 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

Apply repeated Strang-split phase integration 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 phases after n_steps Strang-split steps.

Raises

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

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    n_steps: int,
) -> FloatArray:
    """Apply repeated Strang-split phase integration 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 phases after ``n_steps`` Strang-split steps.

    Raises
    ------
    ValueError
        If ``n_steps`` is negative or the state arrays are invalid.
    """
    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")
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    knm_flat = knm64.ravel()
    alpha_flat = alpha64.ravel()
    backend_fn = _dispatch() if self._dt > 0.0 else None
    # Non-Rust backends also validate dt > 0 internally, so
    # negative dt (symplectic-reversibility checks) falls back
    # to the Python reference.
    if backend_fn is not None:
        return _validate_backend_output(
            backend_fn(
                phases64,
                omegas64,
                knm_flat,
                alpha_flat,
                self._n,
                float(zeta),
                float(psi),
                float(self._dt),
                int(n_steps),
            ),
            n=self._n,
        )
    return _python_run(
        phases64,
        omegas64,
        knm_flat,
        alpha_flat,
        self._n,
        float(zeta),
        float(psi),
        float(self._dt),
        int(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/splitting.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:

Bifurcation Continuation

Traces the synchronization transition R(K) as a function of coupling strength. The incoherent state (R≈0) bifurcates to partial synchronization (R>0) at the critical coupling K_c.

Two interfaces:

• trace_sync_transition(): sweep R(K) over a range of coupling strengths
• find_critical_coupling(): binary search for K_c with configurable precision

Analytical reference: K_c = 2/(π g(0)) for Lorentzian g(ω) with half-width Δ gives K_c = 2Δ (Kuramoto 1975, Strogatz 2000).

Usage:

from scpn_phase_orchestrator.upde.bifurcation import (
    trace_sync_transition, find_critical_coupling,
)

# Sweep R(K) curve
diagram = trace_sync_transition(omegas, K_range=(0, 5), n_points=50)
print(f"K_c ≈ {diagram.K_critical}")

# Precise K_c via binary search
Kc = find_critical_coupling(omegas, tol=0.05)

bifurcation

Bifurcation continuation for Kuramoto synchronisation transitions.

Traces steady-state order parameter R as a function of coupling strength K using pseudo-arclength continuation (Keller 1977). Detects critical coupling K_c where the incoherent state R ≈ 0 bifurcates to partial synchronisation R > 0.

Analytical reference: K_c = 2 / (π g(0)) for Lorentzian g(ω) with half-width Δ → K_c = 2Δ (Kuramoto 1975, Strogatz 2000).

5-backend chain via delegation

The single-trial kernel steady_state_r(phases, omegas, knm, alpha, k_scale, dt, n_transient, n_measure) → R is already dispatched across Rust / Mojo / Julia / Go / Python in :mod:scpn_phase_orchestrator.upde.basin_stability. This module delegates to it rather than re-implementing the Euler trial integrator, which means every trace_sync_transition / find_critical_coupling call in the Python-composite branch inherits the full fallback chain for free.

The two composite Rust kernels — trace_sync_transition_rust (batched K-sweep) and find_critical_coupling_bif_rust (binary search inside Rust) — are preserved as one-shot fast paths: a single FFI call amortises the per-K boundary overhead better than the N_points × dispatch-call path.

Classes

BifurcationPoint dataclass

BifurcationPoint(K: float, R: float, stable: bool)

One sampled point on a Kuramoto synchronisation branch.

BifurcationDiagram dataclass

BifurcationDiagram(
    points: list[BifurcationPoint] = list(),
    K_critical: float | None = None,
)

Ordered bifurcation samples plus optional critical coupling.

Attributes

K_values property

K_values: FloatArray

Return continuation coupling strengths in diagram order.

Returns

FloatArray Return continuation coupling strengths in diagram order.

R_values property

R_values: FloatArray

Return continuation order parameters in diagram order.

Returns

FloatArray Return continuation order parameters in diagram order.

Functions:

trace_sync_transition

trace_sync_transition(
    omegas: FloatArray,
    knm_template: FloatArray | None = None,
    alpha: FloatArray | None = None,
    K_range: tuple[float, float] = (0.0, 5.0),
    n_points: int = 50,
    dt: float = 0.01,
    n_transient: int = 2000,
    n_measure: int = 500,
    seed: int = 42,
) -> BifurcationDiagram

Trace R(K) for the Kuramoto synchronisation transition.

Sweeps coupling strength K from K_range[0] to K_range[1], running the ODE to steady state at each point, and returns a :class:BifurcationDiagram with the (K, R) pairs plus the estimated critical coupling K_c.

When the Rust composite kernel is available, the whole sweep is batched into a single FFI call. Otherwise the function loops in Python and each trial is dispatched through the 5-backend chain inherited from :func:basin_stability.steady_state_r.

Parameters

omegas : FloatArray Finite real numeric natural frequencies in rad/s, shape (N,). Boolean, complex, and numeric-string aliases are rejected. knm_template : FloatArray | None Unit coupling template scaled along the continuation, or None for all-to-all. Must be finite, real numeric, and zero-diagonal. alpha : FloatArray | None Finite real numeric phase-lag matrix in radians, shape (N, N), or None for no lag. K_range : tuple[float, float] Inclusive (min, max) coupling-strength range to scan. n_points : int Number of coupling points sampled across the range. dt : float Integration step size. n_transient : int Number of transient steps discarded before measurement. n_measure : int Number of steps averaged to measure the order parameter. seed : int Seed for the deterministic RNG.

Returns

BifurcationDiagram The traced R(K) bifurcation diagram.

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

def trace_sync_transition(
    omegas: FloatArray,
    knm_template: FloatArray | None = None,
    alpha: FloatArray | None = None,
    K_range: tuple[float, float] = (0.0, 5.0),
    n_points: int = 50,
    dt: float = 0.01,
    n_transient: int = 2000,
    n_measure: int = 500,
    seed: int = 42,
) -> BifurcationDiagram:
    """Trace R(K) for the Kuramoto synchronisation transition.

    Sweeps coupling strength ``K`` from ``K_range[0]`` to
    ``K_range[1]``, running the ODE to steady state at each point,
    and returns a :class:`BifurcationDiagram` with the ``(K, R)``
    pairs plus the estimated critical coupling ``K_c``.

    When the Rust composite kernel is available, the whole sweep
    is batched into a single FFI call. Otherwise the function
    loops in Python and each trial is dispatched through the
    5-backend chain inherited from
    :func:`basin_stability.steady_state_r`.

    Parameters
    ----------
    omegas : FloatArray
        Finite real numeric natural frequencies in rad/s, shape ``(N,)``.
        Boolean, complex, and numeric-string aliases are rejected.
    knm_template : FloatArray | None
        Unit coupling template scaled along the continuation, or ``None`` for
        all-to-all. Must be finite, real numeric, and zero-diagonal.
    alpha : FloatArray | None
        Finite real numeric phase-lag matrix in radians, shape ``(N, N)``, or
        ``None`` for no lag.
    K_range : tuple[float, float]
        Inclusive ``(min, max)`` coupling-strength range to scan.
    n_points : int
        Number of coupling points sampled across the range.
    dt : float
        Integration step size.
    n_transient : int
        Number of transient steps discarded before measurement.
    n_measure : int
        Number of steps averaged to measure the order parameter.
    seed : int
        Seed for the deterministic RNG.

    Returns
    -------
    BifurcationDiagram
        The traced ``R(K)`` bifurcation diagram.
    """
    omegas = _validate_omegas(omegas)
    n = int(omegas.shape[0])
    K_range = _validate_k_range(K_range)
    n_points = _validate_integral(n_points, name="n_points", minimum=2)
    dt = _validate_positive_float(dt, name="dt")
    n_transient = _validate_integral(n_transient, name="n_transient", minimum=0)
    n_measure = _validate_integral(n_measure, name="n_measure", minimum=0)
    seed = _validate_integral(seed, name="seed", minimum=0)

    if knm_template is None:
        knm_template = _default_coupling(n)
    else:
        knm_template = _validate_matrix(
            knm_template,
            name="knm_template",
            n=n,
            require_zero_diagonal=True,
        )
    if alpha is None:
        alpha = np.zeros((n, n), dtype=np.float64)
    else:
        alpha = _validate_matrix(alpha, name="alpha", n=n)

    rng = np.random.default_rng(seed)
    phases_init = rng.uniform(0, 2 * np.pi, n)
    diagram = BifurcationDiagram()

    if _HAS_COMPOSITE_RUST:
        o = np.ascontiguousarray(omegas, dtype=np.float64)
        k = np.ascontiguousarray(knm_template.ravel(), dtype=np.float64)
        a = np.ascontiguousarray(alpha.ravel(), dtype=np.float64)
        p = np.ascontiguousarray(phases_init, dtype=np.float64)
        kv, rv, kc = _rust_trace(
            o,
            k,
            a,
            n,
            p,
            K_range[0],
            K_range[1],
            n_points,
            dt,
            n_transient,
            n_measure,
        )
        kv, rv = _validate_rust_trace_result(
            kv,
            rv,
            n_points=n_points,
            K_range=K_range,
        )
        critical = _validate_optional_critical_coupling(kc)
        for i in range(len(kv)):
            diagram.points.append(
                BifurcationPoint(
                    K=float(kv[i]),
                    R=float(rv[i]),
                    stable=True,
                ),
            )
        if critical is not None:
            diagram.K_critical = critical
        return diagram

    # Composite Rust unavailable — loop in Python, each trial
    # dispatched through the basin_stability 5-backend chain.
    K_values = np.linspace(K_range[0], K_range[1], n_points)
    for K_val in K_values:
        R = _steady_state_R_dispatch(
            phases_init,
            omegas,
            K_val,
            knm_template,
            alpha,
            dt,
            n_transient,
            n_measure,
        )
        diagram.points.append(
            BifurcationPoint(K=float(K_val), R=R, stable=True),
        )

    R_arr = diagram.R_values
    threshold = 0.1
    crossings = np.where(
        (R_arr[:-1] < threshold) & (R_arr[1:] >= threshold),
    )[0]
    if len(crossings) > 0:
        idx = crossings[0]
        K_lo, K_hi = float(K_values[idx]), float(K_values[idx + 1])
        R_lo, R_hi = float(R_arr[idx]), float(R_arr[idx + 1])
        frac = (threshold - R_lo) / (R_hi - R_lo)
        diagram.K_critical = K_lo + frac * (K_hi - K_lo)
    return diagram

find_critical_coupling

find_critical_coupling(
    omegas: FloatArray,
    knm_template: FloatArray | None = None,
    dt: float = 0.01,
    n_transient: int = 3000,
    n_measure: int = 1000,
    tol: float = 0.05,
    seed: int = 42,
) -> float

Binary-search the critical coupling K_c where R crosses 0.1.

More precise than :func:trace_sync_transition when only K_c is needed. Returns nan if no transition is found in [0, 20].

Parameters

omegas : FloatArray Finite real numeric natural frequencies in rad/s, shape (N,). Boolean, complex, and numeric-string aliases are rejected. knm_template : FloatArray | None Unit coupling template scaled along the continuation, or None for all-to-all. Must be finite, real numeric, and zero-diagonal. dt : float Integration step size. n_transient : int Number of transient steps discarded before measurement. n_measure : int Number of steps averaged to measure the order parameter. tol : float Convergence tolerance for the binary search. seed : int Seed for the deterministic RNG.

Returns

float The critical coupling K_c where R first crosses 0.1.

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

def find_critical_coupling(
    omegas: FloatArray,
    knm_template: FloatArray | None = None,
    dt: float = 0.01,
    n_transient: int = 3000,
    n_measure: int = 1000,
    tol: float = 0.05,
    seed: int = 42,
) -> float:
    """Binary-search the critical coupling ``K_c`` where ``R`` crosses 0.1.

    More precise than :func:`trace_sync_transition` when only
    ``K_c`` is needed. Returns ``nan`` if no transition is found
    in ``[0, 20]``.

    Parameters
    ----------
    omegas : FloatArray
        Finite real numeric natural frequencies in rad/s, shape ``(N,)``.
        Boolean, complex, and numeric-string aliases are rejected.
    knm_template : FloatArray | None
        Unit coupling template scaled along the continuation, or ``None`` for
        all-to-all. Must be finite, real numeric, and zero-diagonal.
    dt : float
        Integration step size.
    n_transient : int
        Number of transient steps discarded before measurement.
    n_measure : int
        Number of steps averaged to measure the order parameter.
    tol : float
        Convergence tolerance for the binary search.
    seed : int
        Seed for the deterministic RNG.

    Returns
    -------
    float
        The critical coupling ``K_c`` where ``R`` first crosses 0.1.
    """
    omegas = _validate_omegas(omegas)
    n = int(omegas.shape[0])
    dt = _validate_positive_float(dt, name="dt")
    n_transient = _validate_integral(n_transient, name="n_transient", minimum=0)
    n_measure = _validate_integral(n_measure, name="n_measure", minimum=0)
    tol = _validate_positive_float(tol, name="tol")
    seed = _validate_integral(seed, name="seed", minimum=0)

    if knm_template is None:
        knm_template = _default_coupling(n)
    else:
        knm_template = _validate_matrix(
            knm_template,
            name="knm_template",
            n=n,
            require_zero_diagonal=True,
        )

    alpha = np.zeros((n, n))
    rng = np.random.default_rng(seed)
    phases_init = rng.uniform(0, 2 * np.pi, n)

    if _HAS_COMPOSITE_RUST:
        o = np.ascontiguousarray(omegas, dtype=np.float64)
        k = np.ascontiguousarray(knm_template.ravel(), dtype=np.float64)
        a = np.ascontiguousarray(alpha.ravel(), dtype=np.float64)
        p = np.ascontiguousarray(phases_init, dtype=np.float64)
        return _validate_find_critical_coupling_result(
            _rust_find_kc(
                o,
                k,
                a,
                n,
                p,
                dt,
                n_transient,
                n_measure,
                tol,
            ),
        )

    threshold = 0.1
    K_lo, K_hi = 0.0, 20.0

    R_hi = _steady_state_R_dispatch(
        phases_init,
        omegas,
        K_hi,
        knm_template,
        alpha,
        dt,
        n_transient,
        n_measure,
    )
    if R_hi < threshold:
        return float("nan")

    for _ in range(30):
        K_mid = (K_lo + K_hi) / 2
        R_mid = _steady_state_R_dispatch(
            phases_init,
            omegas,
            K_mid,
            knm_template,
            alpha,
            dt,
            n_transient,
            n_measure,
        )
        if R_mid < threshold:
            K_lo = K_mid
        else:
            K_hi = K_mid
        if K_hi - K_lo < tol:
            break

    return (K_lo + K_hi) / 2

Basin Stability

Monte Carlo estimation of the volume of the basin of attraction for the synchronised state. Basin stability S_B is the probability that a random initial condition converges to the synchronised attractor.

Procedure: Draw n_samples random phase configurations from [0, 2π)^N, integrate each to steady state, check if R_final > R_threshold. S_B = fraction that converge.

multi_basin_stability() classifies outcomes at multiple R thresholds to detect multi-stability (chimera states, partial synchronization).

Optional Rust, Go, Julia, and Mojo backend outputs are validated before public publication: each steady-state order-parameter scalar must be finite, non-boolean, non-numeric-string, and inside [0, 1]. Public and direct phase, frequency, flattened coupling, phase-lag, scalar-control, threshold, and count inputs reject numeric-string aliases before float coercion. Loader/runtime unavailability can still fall through to Python; malformed backend physics evidence fails closed.

Usage:

from scpn_phase_orchestrator.upde.basin_stability import (
    basin_stability, multi_basin_stability,
)

result = basin_stability(omegas, knm, n_samples=1000)
print(f"S_B = {result.S_B:.3f} ({result.n_converged}/{result.n_samples})")

# Multi-threshold detection
results = multi_basin_stability(omegas, knm, R_thresholds=(0.3, 0.6, 0.8))

References: Menck et al. 2013, Nature Physics 9:89-92.

basin_stability

Basin stability for Kuramoto synchronisation with a 5-backend fallback chain.

Monte Carlo estimation of the volume of the basin of attraction for the synchronised state. Basin stability S_B is the probability that a random initial condition converges to the synchronised attractor (Menck et al. 2013, Ji et al. 2014).

Kernel of the computation

The single-trial primitive is steady_state_r(phases_init, omegas, knm, alpha, dt, n_transient, n_measure) → R — explicit Euler integration of the Kuramoto ODE, transient discarded, time-averaged order parameter returned. The trial kernel has no RNG and is dispatched across Rust / Mojo / Julia / Go / Python (bit-exact parity on deterministic inputs).

RNG ownership

The Monte Carlo loop lives in Python: np.random.default_rng(seed) draws n_samples random phase vectors from [0, 2π)^N and calls the dispatched trial kernel once per IC. This is the dimension pattern — Python owns the randomness so the compute primitive stays deterministic and parity-testable. The original basin_stability_rust (seed-in → S_B-out) kernel is preserved as a one-shot fast path when all four arguments match, but regular use goes through the dispatched per-trial kernel.

Functions:

steady_state_r

steady_state_r(
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray | None = None,
    k_scale: float = 1.0,
    dt: float = 0.01,
    n_transient: int = 500,
    n_measure: int = 200,
) -> float

One-trial Kuramoto steady-state R (dispatched).

Integrates the Kuramoto ODE for n_transient + n_measure steps and returns the time-averaged order parameter over the latter window. Delegates to the fastest available backend.

Parameters

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 | None Phase-lag matrix in radians, shape (N, N), or None for no lag. k_scale : float Multiplicative scale applied to the coupling matrix. dt : float Integration step size. n_transient : int Number of transient steps discarded before measurement. n_measure : int Number of steps averaged to measure the order parameter.

Returns

float The steady-state Kuramoto order parameter R of the trial.

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

def steady_state_r(
    phases_init: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray | None = None,
    k_scale: float = 1.0,
    dt: float = 0.01,
    n_transient: int = 500,
    n_measure: int = 200,
) -> float:
    """One-trial Kuramoto steady-state R (dispatched).

    Integrates the Kuramoto ODE for ``n_transient + n_measure`` steps
    and returns the time-averaged order parameter over the latter
    window. Delegates to the fastest available backend.

    Parameters
    ----------
    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 | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    k_scale : float
        Multiplicative scale applied to the coupling matrix.
    dt : float
        Integration step size.
    n_transient : int
        Number of transient steps discarded before measurement.
    n_measure : int
        Number of steps averaged to measure the order parameter.

    Returns
    -------
    float
        The steady-state Kuramoto order parameter ``R`` of the trial.
    """
    phases_init = _validate_nonempty_vector(phases_init, name="phases_init")
    N = int(phases_init.shape[0])
    omegas = _validate_vector(omegas, name="omegas", shape=(N,))
    knm = _validate_vector(knm, name="knm", shape=(N, N))
    if alpha is None:
        alpha_flat = np.zeros(N * N, dtype=np.float64)
    else:
        alpha_flat = _validate_vector(alpha, name="alpha", shape=(N, N)).ravel()
    k_scale = _validate_finite_float(k_scale, name="k_scale")
    dt = _validate_positive_float(dt, name="dt")
    n_transient = _validate_integral(n_transient, name="n_transient", minimum=0)
    n_measure = _validate_integral(n_measure, name="n_measure", minimum=0)
    if n_measure == 0:
        return 0.0
    knm_flat = knm.ravel()
    backend_fn = _dispatch()
    if backend_fn is not None:
        return _basin_stability_validation.validate_basin_stability_output(
            backend_fn(
                phases_init,
                omegas,
                knm_flat,
                alpha_flat,
                N,
                k_scale,
                dt,
                n_transient,
                n_measure,
            )
        )
    return _python_steady_state_r(
        phases_init,
        omegas,
        knm_flat,
        alpha_flat,
        N,
        k_scale,
        dt,
        n_transient,
        n_measure,
    )

basin_stability

basin_stability(
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray | None = None,
    dt: float = 0.01,
    n_transient: int = 500,
    n_measure: int = 200,
    n_samples: int = 100,
    R_threshold: float = 0.8,
    seed: int = 42,
) -> BasinStabilityResult

Estimate basin stability of the synchronised state.

Draws n_samples random initial phase configurations from [0, 2π)^N, integrates each to steady state via the dispatched trial kernel, and classifies trials by R_final ≥ R_threshold.

Parameters

omegas : FloatArray (N,) natural frequencies. knm : FloatArray (N, N) coupling matrix. alpha : FloatArray | None (N, N) phase lags (default: zeros). dt : float Integration timestep. n_transient : int Transient steps to discard. n_measure : int Steps to average R over. n_samples : int Number of random initial conditions. R_threshold : float Threshold for classifying as "synchronised". seed : int RNG seed (owned by Python).

Returns

BasinStabilityResult BasinStabilityResult with S_B, R_final array, and counts.

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

def basin_stability(
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray | None = None,
    dt: float = 0.01,
    n_transient: int = 500,
    n_measure: int = 200,
    n_samples: int = 100,
    R_threshold: float = 0.8,
    seed: int = 42,
) -> BasinStabilityResult:
    """Estimate basin stability of the synchronised state.

    Draws ``n_samples`` random initial phase configurations from
    ``[0, 2π)^N``, integrates each to steady state via the dispatched
    trial kernel, and classifies trials by ``R_final ≥ R_threshold``.

    Parameters
    ----------
    omegas : FloatArray
        (N,) natural frequencies.
    knm : FloatArray
        (N, N) coupling matrix.
    alpha : FloatArray | None
        (N, N) phase lags (default: zeros).
    dt : float
        Integration timestep.
    n_transient : int
        Transient steps to discard.
    n_measure : int
        Steps to average R over.
    n_samples : int
        Number of random initial conditions.
    R_threshold : float
        Threshold for classifying as "synchronised".
    seed : int
        RNG seed (owned by Python).

    Returns
    -------
    BasinStabilityResult
        BasinStabilityResult with S_B, R_final array, and counts.
    """
    omegas = _validate_omegas(omegas)
    N = int(omegas.shape[0])
    knm = _validate_vector(knm, name="knm", shape=(N, N))
    if alpha is None:
        alpha_flat = np.zeros(N * N, dtype=np.float64)
    else:
        alpha_flat = _validate_vector(alpha, name="alpha", shape=(N, N)).ravel()
    dt = _validate_positive_float(dt, name="dt")
    n_transient = _validate_integral(n_transient, name="n_transient", minimum=0)
    n_measure = _validate_integral(n_measure, name="n_measure", minimum=0)
    n_samples = _validate_integral(n_samples, name="n_samples", minimum=0)
    R_threshold = _validate_unit_interval(R_threshold, name="R_threshold")
    seed = _validate_integral(seed, name="seed", minimum=0)

    R_finals = _monte_carlo_R_finals(
        omegas,
        knm.ravel(),
        alpha_flat,
        N,
        dt,
        n_transient,
        n_measure,
        n_samples,
        seed,
    )
    n_converged = int(np.sum(R_finals >= R_threshold))
    return BasinStabilityResult(
        S_B=n_converged / n_samples if n_samples > 0 else 0.0,
        n_samples=n_samples,
        n_converged=n_converged,
        R_final=R_finals,
        R_threshold=R_threshold,
    )

multi_basin_stability

multi_basin_stability(
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray | None = None,
    dt: float = 0.01,
    n_transient: int = 500,
    n_measure: int = 200,
    n_samples: int = 100,
    R_thresholds: tuple[float, ...] = (0.3, 0.6, 0.8),
    seed: int = 42,
) -> dict[str, BasinStabilityResult]

Basin stability at multiple synchronisation thresholds.

One Monte Carlo sweep; threshold classification repeated locally for each entry of R_thresholds.

Returns

Dict mapping ``"R>={thresh:.2f}"`` to BasinStabilityResult.

Parameters

omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). alpha : FloatArray | None Phase-lag matrix in radians, shape (N, N), or None for no lag. dt : float Integration step size. n_transient : int Number of transient steps discarded before measurement. n_measure : int Number of steps averaged to measure the order parameter. n_samples : int Number of random initial-condition samples. R_thresholds : tuple[float, ...] Order-parameter thresholds to evaluate basin stability at. seed : int Seed for the deterministic RNG.

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

def multi_basin_stability(
    omegas: FloatArray,
    knm: FloatArray,
    alpha: FloatArray | None = None,
    dt: float = 0.01,
    n_transient: int = 500,
    n_measure: int = 200,
    n_samples: int = 100,
    R_thresholds: tuple[float, ...] = (0.3, 0.6, 0.8),
    seed: int = 42,
) -> dict[str, BasinStabilityResult]:
    """Basin stability at multiple synchronisation thresholds.

    One Monte Carlo sweep; threshold classification repeated locally
    for each entry of ``R_thresholds``.

    Returns
    -------
        Dict mapping ``"R>={thresh:.2f}"`` to BasinStabilityResult.

    Parameters
    ----------
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    alpha : FloatArray | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    dt : float
        Integration step size.
    n_transient : int
        Number of transient steps discarded before measurement.
    n_measure : int
        Number of steps averaged to measure the order parameter.
    n_samples : int
        Number of random initial-condition samples.
    R_thresholds : tuple[float, ...]
        Order-parameter thresholds to evaluate basin stability at.
    seed : int
        Seed for the deterministic RNG.
    """
    omegas = _validate_omegas(omegas)
    N = int(omegas.shape[0])
    knm = _validate_vector(knm, name="knm", shape=(N, N))
    if alpha is None:
        alpha_flat = np.zeros(N * N, dtype=np.float64)
    else:
        alpha_flat = _validate_vector(alpha, name="alpha", shape=(N, N)).ravel()
    dt = _validate_positive_float(dt, name="dt")
    n_transient = _validate_integral(n_transient, name="n_transient", minimum=0)
    n_measure = _validate_integral(n_measure, name="n_measure", minimum=0)
    n_samples = _validate_integral(n_samples, name="n_samples", minimum=0)
    R_thresholds = _validate_thresholds(R_thresholds)
    seed = _validate_integral(seed, name="seed", minimum=0)

    R_finals = _monte_carlo_R_finals(
        omegas,
        knm.ravel(),
        alpha_flat,
        N,
        dt,
        n_transient,
        n_measure,
        n_samples,
        seed,
    )
    results: dict[str, BasinStabilityResult] = {}
    for thresh in R_thresholds:
        n_above = int(np.sum(R_finals >= thresh))
        results[f"R>={thresh:.2f}"] = BasinStabilityResult(
            S_B=n_above / n_samples if n_samples > 0 else 0.0,
            n_samples=n_samples,
            n_converged=n_above,
            R_final=R_finals,
            R_threshold=thresh,
        )
    return results

Hypergraph (k-Body) Coupling Engine

Generalized k-body Kuramoto interactions via explicit hyperedge lists. Extends beyond the simplicial engine's fixed 3-body coupling to arbitrary k-body interactions for any k ≥ 2.

For a k-hyperedge {i₁, ..., iₖ}, the coupling on oscillator iₘ is: σₖ · sin(Σ_{j≠m} θ_{iⱼ} - (k-1)·θ_{iₘ})

This generalizes: - k=2: sin(θ_j - θ_i) — standard Kuramoto - k=3: sin(θ_j + θ_k - 2θ_i) — simplicial - k=4: sin(θ_j + θ_k + θ_l - 3θ_i) — quartic interaction

Supports mixed-order interactions: some edges pairwise, some 3-body, some 4-body, in the same network.

Optional hypergraph backends are validated before their results publish: returned phase vectors must keep oscillator cardinality, contain finite real values, and remain in [0, 2*pi). Loader/runtime unavailability can still fall back to Python; malformed backend outputs raise instead of becoming simulation evidence. Public phase, frequency, optional matrix, scalar, count, and order-parameter inputs, direct backend vectors, index buffers, scalar controls, and backend outputs reject numeric-string aliases before Python, NumPy, or accelerator coercion.

Usage:

from scpn_phase_orchestrator.upde.hypergraph import HypergraphEngine

eng = HypergraphEngine(n_oscillators=8, dt=0.01)
eng.add_all_to_all(order=3, strength=0.5)  # all 3-body edges
eng.add_edge((0, 1, 2, 3), strength=0.2)   # one 4-body edge

phases = eng.run(phases_init, omegas, n_steps=1000,
                 pairwise_knm=knm)  # combine with standard coupling

References: Tanaka & Aoyagi 2011, Phys. Rev. Lett. 106:224101; Bick et al. 2023, Nat. Rev. Physics 5:307-317. Detailed documentation: Hypergraph (k-body) — detailed reference

hypergraph

Hypergraph Kuramoto with arbitrary k-body interactions beyond pairwise.

Exposes a 5-backend fallback chain.

Extends the standard Kuramoto model with k-body coupling terms for any k ≥ 2. The standard model (k=2) and simplicial model (k=3) are special cases.

For a k-hyperedge {i₁, …, iₖ}, the coupling on oscillator iₘ is

σₖ · sin( Σ_{j≠m} θ_{iⱼ} − (k−1)·θ_{iₘ} )

which generalises

k = 2:  sin(θ_j − θ_i)            — standard Kuramoto
k = 3:  sin(θ_j + θ_k − 2·θ_i)    — simplicial / triadic

The engine also accepts a dense pairwise coupling matrix pairwise_knm (optional) and an external-drive field (ζ, ψ).

Numerics

The pairwise-derivative loop uses the Rust kernel's sin(θ_j − θ_i) = sin(θ_j)·cos(θ_i) − cos(θ_j)·sin(θ_i) expansion in the alpha == 0 fast path so that floating-point rounding matches Rust (spo-engine/src/hypergraph.rs) bit-for-bit. Alpha ≠ 0 falls back to the direct sin(θ_j − θ_i − α) form in all five backends.

References

Tanaka & Aoyagi 2011, Phys. Rev. Lett. 106:224101. Skardal & Arenas 2019, Comm. Phys. 2:22. Bick, Gross, Harrington & Schaub 2023, Nat. Rev. Physics 5:307-317.

Classes

Hyperedge dataclass

Hyperedge(nodes: tuple[int, ...], strength: float = 1.0)

A k-body interaction among oscillators.

Attributes

nodes: Tuple of oscillator indices in this hyperedge.
strength: Coupling strength σₖ for this hyperedge.

Attributes

order property

order: int

Return the number of oscillators participating in the hyperedge.

Returns

int Return the number of oscillators participating in the hyperedge.

HypergraphEngine

HypergraphEngine(
    n_oscillators: int,
    dt: float,
    hyperedges: list[Hyperedge] | None = None,
)

Kuramoto engine with arbitrary k-body hypergraph coupling.

Supports mixed-order interactions: some edges can be pairwise, some 3-body, some 4-body, etc. Each Hyperedge specifies which oscillators participate and the coupling strength.

Initialise a validated hypergraph Kuramoto engine.

Parameters

n_oscillators : int Positive number of oscillators in the simulated system. dt : float Positive finite explicit-Euler step size. hyperedges : list[Hyperedge] | None Optional initial hyperedge definitions to validate and store.

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

def __init__(
    self,
    n_oscillators: int,
    dt: float,
    hyperedges: list[Hyperedge] | None = None,
) -> None:
    """Initialise a validated hypergraph Kuramoto engine.

    Parameters
    ----------
    n_oscillators : int
        Positive number of oscillators in the simulated system.
    dt : float
        Positive finite explicit-Euler step size.
    hyperedges : list[Hyperedge] | None
        Optional initial hyperedge definitions to validate and store.
    """
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._dt = _validate_positive_float(dt, name="dt")
    self._hyperedges = [
        _validate_hyperedge(edge, n_oscillators=self._n)
        for edge in (hyperedges or [])
    ]

Attributes

n_edges property

n_edges: int

Return the number of configured hyperedges.

Returns

int Return the number of configured hyperedges.

Methods:

add_edge

add_edge(
    nodes: tuple[int, ...], strength: float = 1.0
) -> None

Validate and append one explicit k-body hyperedge.

Parameters

nodes : tuple[int, ...] Indices of the oscillators participating in the hyperedge. strength : float Coupling strength assigned to the hyperedge(s).

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

def add_edge(self, nodes: tuple[int, ...], strength: float = 1.0) -> None:
    """Validate and append one explicit k-body hyperedge.

    Parameters
    ----------
    nodes : tuple[int, ...]
        Indices of the oscillators participating in the hyperedge.
    strength : float
        Coupling strength assigned to the hyperedge(s).
    """
    edge = _validate_hyperedge(
        Hyperedge(nodes=nodes, strength=strength),
        n_oscillators=self._n,
    )
    self._hyperedges.append(edge)

add_all_to_all

add_all_to_all(order: int, strength: float = 1.0) -> None

Add all C(N, order) hyperedges of given order.

Parameters

order : int Interaction order (number of oscillators per hyperedge). strength : float Coupling strength assigned to the hyperedge(s).

Raises

ValueError If order is outside 2..N.

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

def add_all_to_all(self, order: int, strength: float = 1.0) -> None:
    """Add all C(N, order) hyperedges of given order.

    Parameters
    ----------
    order : int
        Interaction order (number of oscillators per hyperedge).
    strength : float
        Coupling strength assigned to the hyperedge(s).

    Raises
    ------
    ValueError
        If ``order`` is outside ``2..N``.
    """
    from itertools import combinations

    order = _validate_positive_int(order, name="order")
    if order > self._n:
        raise ValueError(f"order must be <= n_oscillators, got {order!r}")
    validated_edge = _validate_hyperedge(
        Hyperedge(nodes=tuple(range(order)), strength=strength),
        n_oscillators=self._n,
    )
    for combo in combinations(range(self._n), order):
        self._hyperedges.append(
            Hyperedge(nodes=combo, strength=validated_edge.strength)
        )

step

step(
    phases: FloatArray,
    omegas: FloatArray,
    pairwise_knm: FloatArray | None = None,
    alpha: FloatArray | None = None,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray

One explicit-Euler step.

Parameters

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

Returns

FloatArray The phases after one explicit-Euler hypergraph step.

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    pairwise_knm: FloatArray | None = None,
    alpha: FloatArray | None = None,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray:
    """One explicit-Euler step.

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

    Returns
    -------
    FloatArray
        The phases after one explicit-Euler hypergraph step.
    """
    return self.run(
        phases,
        omegas,
        n_steps=1,
        pairwise_knm=pairwise_knm,
        alpha=alpha,
        zeta=zeta,
        psi=psi,
    )

run

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

Integrate n_steps Euler steps via the fastest backend; return phases.

Parameters

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

Returns

FloatArray The final phases after n_steps hypergraph steps.

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    n_steps: int,
    pairwise_knm: FloatArray | None = None,
    alpha: FloatArray | None = None,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray:
    """Integrate ``n_steps`` Euler steps via the fastest backend; return phases.

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

    Returns
    -------
    FloatArray
        The final phases after ``n_steps`` hypergraph 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,))
    knm_flat = _validate_optional_state_array(
        pairwise_knm,
        name="pairwise_knm",
        shape=(self._n, self._n),
    )
    alpha_flat = _validate_optional_state_array(
        alpha,
        name="alpha",
        shape=(self._n, self._n),
    )
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    en, eo, es = self._encode_edges()
    backend_fn = _dispatch()
    if backend_fn is not None:
        return _validate_backend_output(
            backend_fn(
                phases64,
                omegas64,
                self._n,
                en,
                eo,
                es,
                knm_flat,
                alpha_flat,
                float(zeta),
                float(psi),
                float(self._dt),
                int(n_steps),
            ),
            n=self._n,
        )
    return _validate_backend_output(
        _python_run(
            phases64,
            omegas64,
            self._n,
            en,
            eo,
            es,
            knm_flat,
            alpha_flat,
            float(zeta),
            float(psi),
            float(self._dt),
            int(n_steps),
        ),
        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/hypergraph.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:

Strang Splitting Engine

Symmetric operator splitting: A(dt/2) → B(dt) → A(dt/2) where A is exact rotation (ω·dt) and B is RK4 on coupling. Second-order accurate, time-reversible, preserves symplectic structure approximately.

Direct Go, Julia, and Mojo Strang-splitting accelerator entrypoints share a validated torus boundary before optional runtime loading: phase and frequency vectors must be finite real one-dimensional float64 arrays matching the oscillator count; flattened pairwise coupling and phase-lag buffers must have exactly N*N values; pairwise self-coupling K_ii must be zero; zeta and psi must be finite controls; and direct accelerator dt plus n_steps must be positive. Backend outputs must be finite torus phases in [0, 2*pi), and public state arrays, direct vectors/matrices, and backend outputs reject numeric-string aliases before float coercion. Mojo stdout must contain exactly one phase line per oscillator. The public SplittingEngine still supports negative dt for reversibility checks by using the Python reference path instead of direct optional accelerators. Detailed documentation: Strang Splitting — detailed reference

splitting

Strang second-order operator splitting for the Kuramoto ODE.

Exposes a 5-backend fallback chain.

Scheme

Split dθ/dt = ω + Σ_j K_ij · sin(θ_j − θ_i − α_ij) + ζ · sin(ψ − θ_i) into

A: dθ/dt = ω              (exact rotation)
B: dθ/dt = coupling       (RK4 on the nonlinear part)

and compose symmetrically as A(dt/2) → B(dt) → A(dt/2) (Strang scheme, second-order in dt).

Why split?

The ω flow is linear, so it has no truncation error; folding it into a monolithic RK45 burns integrator budget on a solvable direction while damping unrelated accuracy in the nonlinear direction. Reference: Hairer, Lubich & Wanner 2006, Geometric Numerical Integration §II.5.

Numerics

The B-stage RK4 uses the Rust kernel's sin(θ_j − θ_i) = sin(θ_j)·cos(θ_i) − cos(θ_j)·sin(θ_i) expansion on the alpha-zero branch so that floating-point rounding matches Rust (spo-engine/src/splitting.rs) bit-for-bit. Nonzero alpha falls back to the direct sin(diff) form in all five backends.

Classes

SplittingEngine

SplittingEngine(n_oscillators: int, dt: float)

Strang-split Kuramoto stepper with 5-backend dispatch.

The engine's geometry is (n, dt); the step is stateless.

Create a Strang-splitting engine for n_oscillators and dt.

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

def __init__(self, n_oscillators: int, dt: float):
    """Create a Strang-splitting engine for ``n_oscillators`` and ``dt``."""
    n_oscillators = _validate_positive_int(
        n_oscillators,
        name="n_oscillators",
    )
    dt = _validate_nonzero_finite_float(dt, name="dt")
    # Negative dt is intentional for symplectic-reversibility
    # checks — Strang is time-reversible, so stepping forward
    # then with dt → −dt returns to the starting state.
    self._n = n_oscillators
    self._dt = dt

Methods:

step

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

One Strang-split step: A(dt/2) → B(dt) → A(dt/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 Strang-split step.

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
) -> FloatArray:
    """One Strang-split step: A(dt/2) → B(dt) → A(dt/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 Strang-split 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

Apply repeated Strang-split phase integration 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 phases after n_steps Strang-split steps.

Raises

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

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float,
    psi: float,
    alpha: FloatArray,
    n_steps: int,
) -> FloatArray:
    """Apply repeated Strang-split phase integration 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 phases after ``n_steps`` Strang-split steps.

    Raises
    ------
    ValueError
        If ``n_steps`` is negative or the state arrays are invalid.
    """
    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")
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    knm_flat = knm64.ravel()
    alpha_flat = alpha64.ravel()
    backend_fn = _dispatch() if self._dt > 0.0 else None
    # Non-Rust backends also validate dt > 0 internally, so
    # negative dt (symplectic-reversibility checks) falls back
    # to the Python reference.
    if backend_fn is not None:
        return _validate_backend_output(
            backend_fn(
                phases64,
                omegas64,
                knm_flat,
                alpha_flat,
                self._n,
                float(zeta),
                float(psi),
                float(self._dt),
                int(n_steps),
            ),
            n=self._n,
        )
    return _python_run(
        phases64,
        omegas64,
        knm_flat,
        alpha_flat,
        self._n,
        float(zeta),
        float(psi),
        float(self._dt),
        int(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/splitting.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:

---

Sparse Engine

SparseUPDEEngine implements the Kuramoto model using a CSR (compressed sparse row) coupling matrix. This reduces memory from O(N^2) dense coupling storage to O(N + E), where E is the number of active directed edge connections.

It is designed for large-scale simulations (national power grids, social networks) where most oscillators are only coupled to local neighbours.

Features

• Scalability: Uses CSR row pointers, column indices, coupling values, and phase-lag values so sparse topologies avoid dense N x N allocation.
• FFI parity: Offloads sparse integration to the Rust backend when spo_kernel is available, with Python fallback preserving the same shape, finite-value, phase-bounds, and adaptive-timestep contracts. Optional-Rust step/run output must remain inside [0, 2*pi), and the backend last_dt must be positive and finite before the public diagnostic is updated.
• Input validation: Rejects malformed CSR row pointers, invalid oscillator indices, non-finite phase/frequency/coupling arrays, unsupported methods, and malformed optional-backend outputs before downstream workflows consume the result.

sparse_engine

Sparse CSR-style UPDE engine for validated oscillator coupling graphs.

The sparse engine advances phase vectors from row pointers, column indices, coupling values, and phase-lag values instead of dense N x N matrices. Inputs are checked for finite values, CSR monotonicity, edge-count consistency, valid oscillator indices, and method selection before stepping. Optional Rust execution and Python fallback preserve the same shape and bounds contracts.

Classes

SparseUPDEEngine

SparseUPDEEngine(
    n_oscillators: int,
    dt: float,
    method: str = "euler",
    atol: float = 1e-06,
    rtol: float = 0.001,
)

Kuramoto UPDE integrator with sparse coupling matrix support.

The SparseUPDEEngine solves the Universal Phase Dynamics Equation (UPDE) using a CSR (Compressed Sparse Row) representation for the coupling matrix K_nm and phase lags alpha_nm. This is critical for scaling to large-scale oscillator networks (e.g., N > 10,000) where the dense K_nm matrix would consume terabytes of RAM.

Mathematics: dtheta_i/dt = omega_i + sum_{j in neighbors(i)} K_ij sin(theta_j - theta_i - alpha_ij) + zeta sin(Psi - theta_i)

The integrator supports sub-microsecond in-place plasticity updates when running on the Rust FFI path, allowing the coupling topology to evolve concurrently with the phase dynamics.

Initialize the sparse integrator.

Parameters

n_oscillators : int Total number of oscillators N in the network. dt : float Integration timestep in seconds. method : str Numerical method ('euler', 'rk4', or 'rk45'). atol : float Absolute tolerance for adaptive RK45. rtol : float Relative tolerance for adaptive RK45.

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

def __init__(
    self,
    n_oscillators: int,
    dt: float,
    method: str = "euler",
    atol: float = 1e-6,
    rtol: float = 1e-3,
):
    """Initialize the sparse integrator.

    Parameters
    ----------
    n_oscillators : int
        Total number of oscillators N in the network.
    dt : float
        Integration timestep in seconds.
    method : str
        Numerical method ('euler', 'rk4', or 'rk45').
    atol : float
        Absolute tolerance for adaptive RK45.
    rtol : float
        Relative tolerance for adaptive RK45.
    """
    n_oscillators = _validate_positive_int(
        n_oscillators,
        name="n_oscillators",
    )
    dt = _validate_positive_float(dt, name="dt")
    atol = _validate_positive_float(atol, name="atol")
    rtol = _validate_positive_float(rtol, name="rtol")
    self._n = n_oscillators
    self._dt = dt
    if method not in ("euler", "rk4", "rk45"):
        msg = f"Unknown method {method!r}, expected 'euler', 'rk4', or 'rk45'"
        raise ValueError(msg)
    self._method = method
    self._atol = atol
    self._rtol = rtol
    self._last_dt = dt
    self._lock = threading.RLock()

    self._rust = None
    if _HAS_RUST:
        try:
            from spo_kernel import PySparseUPDEStepper

            self._rust = PySparseUPDEStepper(
                n_oscillators, dt, method, atol=atol, rtol=rtol
            )
        except ImportError:
            pass

Attributes

last_dt property

last_dt: float

Return the most recent accepted Python or Rust timestep.

Returns

float Positive finite timestep accepted by the sparse engine.

Methods:

step

step(
    phases: FloatArray,
    omegas: FloatArray,
    row_ptr: IntArray,
    col_indices: IntArray,
    knm_values: FloatArray,
    zeta: float,
    psi: float,
    alpha_values: FloatArray,
) -> FloatArray

Advance phases by one sparse timestep, return new phases in [0, 2*pi).

Parameters

phases : FloatArray Current phase vector [theta_1, ..., theta_N], shape (N,). omegas : FloatArray Natural frequency vector [omega_1, ..., omega_N], shape (N,). row_ptr : IntArray CSR row pointers, shape (N+1,). col_indices : IntArray CSR column indices, shape (E,). knm_values : FloatArray CSR coupling strengths, shape (E,). zeta : float External forcing strength (global scalar). psi : float Reference phase target (global scalar). alpha_values : FloatArray CSR phase lags, shape (E,).

Returns

FloatArray New phase vector [theta_1(t+dt), ..., theta_N(t+dt)], shape (N,).

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    row_ptr: IntArray,
    col_indices: IntArray,
    knm_values: FloatArray,
    zeta: float,
    psi: float,
    alpha_values: FloatArray,
) -> FloatArray:
    """Advance phases by one sparse timestep, return new phases in [0, 2*pi).

    Parameters
    ----------
    phases : FloatArray
        Current phase vector [theta_1, ..., theta_N], shape (N,).
    omegas : FloatArray
        Natural frequency vector [omega_1, ..., omega_N], shape (N,).
    row_ptr : IntArray
        CSR row pointers, shape (N+1,).
    col_indices : IntArray
        CSR column indices, shape (E,).
    knm_values : FloatArray
        CSR coupling strengths, shape (E,).
    zeta : float
        External forcing strength (global scalar).
    psi : float
        Reference phase target (global scalar).
    alpha_values : FloatArray
        CSR phase lags, shape (E,).

    Returns
    -------
    FloatArray
        New phase vector [theta_1(t+dt), ..., theta_N(t+dt)], shape (N,).
    """
    zeta = _validate_finite_real(zeta, name="zeta")
    psi = _validate_finite_real(psi, name="psi")
    self._validate_inputs(
        phases,
        omegas,
        row_ptr,
        col_indices,
        knm_values,
        alpha_values,
        zeta,
        psi,
    )
    with self._lock:
        if self._rust is not None:
            result = self._rust.step(
                np.ascontiguousarray(phases.ravel(), dtype=np.float64),
                np.ascontiguousarray(omegas.ravel(), dtype=np.float64),
                np.ascontiguousarray(row_ptr.ravel(), dtype=np.uint64),
                np.ascontiguousarray(col_indices.ravel(), dtype=np.uint64),
                np.ascontiguousarray(knm_values.ravel(), dtype=np.float64),
                zeta,
                psi,
                np.ascontiguousarray(alpha_values.ravel(), dtype=np.float64),
            )
            output = self._validate_rust_output(result)
            self._last_dt = _validate_positive_float(
                self._rust.last_dt,
                name="Rust last_dt",
            )
            return output

        if self._method == "euler":
            return self._euler_step(
                phases,
                omegas,
                row_ptr,
                col_indices,
                knm_values,
                zeta,
                psi,
                alpha_values,
            )
        if self._method == "rk45":
            return self._rk45_step(
                phases,
                omegas,
                row_ptr,
                col_indices,
                knm_values,
                zeta,
                psi,
                alpha_values,
            )
        return self._rk4_step(
            phases,
            omegas,
            row_ptr,
            col_indices,
            knm_values,
            zeta,
            psi,
            alpha_values,
        )

run

run(
    phases: FloatArray,
    omegas: FloatArray,
    row_ptr: IntArray,
    col_indices: IntArray,
    knm_values: FloatArray,
    zeta: float,
    psi: float,
    alpha_values: FloatArray,
    n_steps: int,
) -> FloatArray

Run multiple steps in a batch, return final phases.

Parameters

phases : FloatArray Initial phase vector. omegas : FloatArray Natural frequencies. row_ptr : IntArray CSR row pointers. col_indices : IntArray CSR column indices. knm_values : FloatArray CSR coupling strengths. zeta : float External forcing strength. psi : float Reference phase target. alpha_values : FloatArray CSR phase lags. n_steps : int Number of integration steps to perform.

Returns

FloatArray Final phase vector after n_steps.

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    row_ptr: IntArray,
    col_indices: IntArray,
    knm_values: FloatArray,
    zeta: float,
    psi: float,
    alpha_values: FloatArray,
    n_steps: int,
) -> FloatArray:
    """Run multiple steps in a batch, return final phases.

    Parameters
    ----------
    phases : FloatArray
        Initial phase vector.
    omegas : FloatArray
        Natural frequencies.
    row_ptr : IntArray
        CSR row pointers.
    col_indices : IntArray
        CSR column indices.
    knm_values : FloatArray
        CSR coupling strengths.
    zeta : float
        External forcing strength.
    psi : float
        Reference phase target.
    alpha_values : FloatArray
        CSR phase lags.
    n_steps : int
        Number of integration steps to perform.

    Returns
    -------
    FloatArray
        Final phase vector after n_steps.
    """
    n_steps = _validate_nonnegative_int(n_steps, name="n_steps")
    zeta = _validate_finite_real(zeta, name="zeta")
    psi = _validate_finite_real(psi, name="psi")
    self._validate_inputs(
        phases,
        omegas,
        row_ptr,
        col_indices,
        knm_values,
        alpha_values,
        zeta,
        psi,
    )
    if n_steps == 0:
        return phases.copy()
    with self._lock:
        if self._rust is not None:
            result = self._rust.run(
                np.ascontiguousarray(phases.ravel(), dtype=np.float64),
                np.ascontiguousarray(omegas.ravel(), dtype=np.float64),
                np.ascontiguousarray(row_ptr.ravel(), dtype=np.uint64),
                np.ascontiguousarray(col_indices.ravel(), dtype=np.uint64),
                np.ascontiguousarray(knm_values.ravel(), dtype=np.float64),
                zeta,
                psi,
                np.ascontiguousarray(alpha_values.ravel(), dtype=np.float64),
                n_steps,
            )
            output = self._validate_rust_output(result)
            self._last_dt = _validate_positive_float(
                self._rust.last_dt,
                name="Rust last_dt",
            )
            return output

        p = phases.copy()
        for _ in range(n_steps):
            p = self.step(
                p, omegas, row_ptr, col_indices, knm_values, zeta, psi, alpha_values
            )
        return p

---

Cellular Sheaf Engine

SheafUPDEEngine extends the Kuramoto model from scalar phases to multi-dimensional phase vectors. This implements a cellular-sheaf model of synchronization.

Instead of a single phase \(\theta_i\), each oscillator maintains a vector \(\vec{\theta}_i \in \mathbb{R}^D\). The scalar coupling \(K_{ij}\) is replaced by a restriction map—a block matrix \(B_{ij} \in \mathbb{R}^{D \times D}\) that maps the phase space of node \(j\) into the reference frame of node \(i\).

\[ \dot{\theta}_{i,d} = \omega_{i,d} + \sum_j \sum_k B_{ij}^{dk} \sin(\theta_{j,k} - \theta_{i,d}) + \zeta \sin(\Psi_d - \theta_{i,d}) \]

Features

• Cross-frequency coupling: Dimension \(k\) on node \(j\) can directly drive dimension \(d\) on node \(i\) through off-diagonal elements of \(B_{ij}\).
• Structured topology: Models opinion dynamics, multimodal synchronization, and anisotropic structural constraints natively.
• Rust parity: Uses PySheafUPDEStepper when available while preserving the same public state-shape, numeric-type, finiteness, and torus-domain contracts.
• Fail-closed arrays: Phase, frequency, restriction-map, and drive-target arrays reject boolean, complex, and numeric-string aliases before conversion.
• Fail-closed publication: Rust output must be a finite real flattened N * D torus state, and its adaptive timestep must be positive and finite. A zero-step run returns a validated independent copy without backend dispatch.

sheaf_engine

Cellular-sheaf UPDE integrator for multidimensional oscillator phases.

SheafUPDEEngine advances N x D phase matrices using restriction-map coupling blocks and optional Rust acceleration. It validates oscillator counts, dimensions, timestep/tolerances, solver method, forcing scalars, phase targets, and tensor shapes before integration. Instance-level locks protect reusable scratch buffers so concurrent callers cannot corrupt adaptive or fixed-step solver state.

Classes

SheafUPDEEngine

SheafUPDEEngine(
    n_oscillators: int,
    d_dimensions: int,
    dt: float,
    method: str = "euler",
    atol: float = 1e-06,
    rtol: float = 0.001,
)

Cellular Sheaf UPDE integrator for multi-dimensional phase vectors.

Phase per oscillator is a vector of dimension D. Restriction maps (coupling blocks) B_ij are D x D matrices mapping the phase space of oscillator j into the space of oscillator i.

Mathematics: d(theta_{i,d})/dt = omega_{i,d} + sum_j sum_k B_ij^{dk} sin(theta_{j,k} - theta_{i,d}) + zeta * sin(Psi_d - theta_{i,d})

This enables complex cross-frequency coupling and opinion dynamics over multidimensional belief spaces.

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

def __init__(
    self,
    n_oscillators: int,
    d_dimensions: int,
    dt: float,
    method: str = "euler",
    atol: float = 1e-6,
    rtol: float = 1e-3,
):
    n_oscillators = _validate_positive_int(
        n_oscillators,
        name="n_oscillators",
    )
    d_dimensions = _validate_positive_int(
        d_dimensions,
        name="d_dimensions",
    )
    dt = _validate_positive_float(dt, name="dt")
    atol = _validate_positive_float(atol, name="atol")
    rtol = _validate_positive_float(rtol, name="rtol")
    self._n = n_oscillators
    self._d = d_dimensions
    self._dt = dt
    if method not in ("euler", "rk4", "rk45"):
        msg = f"Unknown method {method!r}, expected 'euler', 'rk4', or 'rk45'"
        raise ValueError(msg)
    self._method = method
    self._atol = atol
    self._rtol = rtol
    self._last_dt = dt
    self._lock = threading.RLock()

    self._rust = None
    if _HAS_RUST:
        try:
            from spo_kernel import PySheafUPDEStepper

            self._rust = PySheafUPDEStepper(
                n_oscillators, d_dimensions, dt, method, atol=atol, rtol=rtol
            )
        except ImportError:
            pass

Attributes

last_dt property

last_dt: float

Return the most recent accepted Python or Rust timestep.

Returns

float Positive finite timestep accepted by the sheaf engine.

Methods:

step

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

Advance phases by one timestep.

Parameters

phases : FloatArray Current phase matrix [theta_i,d], shape (N, D). omegas : FloatArray Natural frequency matrix [omega_i,d], shape (N, D). restriction_maps : FloatArray Block matrix coupling [B_ij^{dk}], shape (N, N, D, D). zeta : float External forcing strength (global scalar). psi : FloatArray Reference phase target vector, shape (D,).

Returns

FloatArray New phase matrix, shape (N, D).

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    restriction_maps: FloatArray,
    zeta: float,
    psi: FloatArray,
) -> FloatArray:
    """Advance phases by one timestep.

    Parameters
    ----------
    phases : FloatArray
        Current phase matrix [theta_i,d], shape (N, D).
    omegas : FloatArray
        Natural frequency matrix [omega_i,d], shape (N, D).
    restriction_maps : FloatArray
        Block matrix coupling [B_ij^{dk}], shape (N, N, D, D).
    zeta : float
        External forcing strength (global scalar).
    psi : FloatArray
        Reference phase target vector, shape (D,).

    Returns
    -------
    FloatArray
        New phase matrix, shape (N, D).
    """
    phases, omegas, restriction_maps, zeta, psi = self._validate_inputs(
        phases,
        omegas,
        restriction_maps,
        zeta,
        psi,
    )
    with self._lock:
        if self._rust is not None:
            res = self._rust.step(
                np.ascontiguousarray(phases.ravel(), dtype=np.float64),
                np.ascontiguousarray(omegas.ravel(), dtype=np.float64),
                np.ascontiguousarray(restriction_maps.ravel(), dtype=np.float64),
                float(zeta),
                np.ascontiguousarray(psi.ravel(), dtype=np.float64),
            )
            output = _reshape_rust_result(
                res,
                name="step",
                shape=(self._n, self._d),
            )
            self._last_dt = _validate_positive_float(
                self._rust.last_dt,
                name="Rust last_dt",
            )
            return output

        if self._method == "euler":
            return self._euler_step(phases, omegas, restriction_maps, zeta, psi)
        if self._method == "rk45":
            return self._rk45_step(phases, omegas, restriction_maps, zeta, psi)
        return self._rk4_step(phases, omegas, restriction_maps, zeta, psi)

run

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

Run multiple steps in a batch, return final phases.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N, D). omegas : FloatArray Natural frequencies in rad/s, shape (N, D). restriction_maps : FloatArray Sheaf restriction maps, shape (N, N, D, D). zeta : float External drive strength ζ. psi : FloatArray External drive reference phase Ψ in radians, shape (D,). n_steps : int Number of integration steps to run. Zero returns an independent, validated copy without invoking the optional backend.

Returns

FloatArray The final phases after n_steps sheaf steps.

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    restriction_maps: FloatArray,
    zeta: float,
    psi: FloatArray,
    n_steps: int,
) -> FloatArray:
    """Run multiple steps in a batch, return final phases.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N, D)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N, D)``.
    restriction_maps : FloatArray
        Sheaf restriction maps, shape ``(N, N, D, D)``.
    zeta : float
        External drive strength ``ζ``.
    psi : FloatArray
        External drive reference phase ``Ψ`` in radians, shape ``(D,)``.
    n_steps : int
        Number of integration steps to run. Zero returns an independent,
        validated copy without invoking the optional backend.

    Returns
    -------
    FloatArray
        The final phases after ``n_steps`` sheaf steps.
    """
    n_steps = _validate_nonnegative_int(n_steps, name="n_steps")
    phases, omegas, restriction_maps, zeta, psi = self._validate_inputs(
        phases,
        omegas,
        restriction_maps,
        zeta,
        psi,
    )
    if n_steps == 0:
        return phases.copy()
    with self._lock:
        if self._rust is not None:
            res = self._rust.run(
                np.ascontiguousarray(phases.ravel(), dtype=np.float64),
                np.ascontiguousarray(omegas.ravel(), dtype=np.float64),
                np.ascontiguousarray(restriction_maps.ravel(), dtype=np.float64),
                float(zeta),
                np.ascontiguousarray(psi.ravel(), dtype=np.float64),
                n_steps,
            )
            output = _reshape_rust_result(
                res,
                name="run",
                shape=(self._n, self._d),
            )
            self._last_dt = _validate_positive_float(
                self._rust.last_dt,
                name="Rust last_dt",
            )
            return output

        p = phases.copy()
        for _ in range(n_steps):
            p = self.step(p, omegas, restriction_maps, zeta, psi)
        return p

Time-varying natural frequencies

UPDEEngine now accepts configured fixed or callable natural frequencies via omega=. If a call omits the omegas argument, the engine resolves the configured source at the current outer-step time and stores the resolved vector in omega_current. Callable schedules are materialised as finite (steps, n) matrices and dispatched through the Rust, Go, Julia, Mojo, or Python schedule runner when available.

Use this for drifting oscillators, moving-agent frequency shifts, chirps, thermal detuning, and Doppler preparation. The detailed contract is documented in UPDE — Time-varying omega.

PHA-C.2 Doppler-corrected UPDE

DopplerEngine adds graph-weighted relative-velocity detuning before each UPDE outer step. It consumes the PHA-C.5 omega(t) schedule contract and is used when phase locking depends on moving oscillators rather than fixed natural frequencies, such as counter-propagating plasmoids, mobile acoustic clocks, robot/sensor swarms, and moving-grid assets.

See UPDE — Doppler Engine for the mathematical contract, scalar/vector velocity handling, backend parity surface, and Mach-1 counter-propagating acceptance scenario.

PHA-C formal proof-obligation bridge

PHACKinematicProofObligation projects a verified end-to-end PHA-C acceptance record into fixed-point Lean obligations for SPOFormal.Kinematic and SPOFormal.Continuous. The manifest is review-only and non-actuating: it binds the accepted timeline hash, acceptance hash, spatial merge-window tolerance, phase tolerance, time-step units, horizon-time units, Gronwall budget trace, continuous horizon-drive budget, and certificate theorem names.

Predictive downstream consumers can now supply both relative_velocity_step_bound_m, coupling_residual_step_bound_m, and phase_drift_bound_rad when building the obligation. The residual bound is recorded separately as configured_coupling_residual_step_bound_units, then combined with the observed moving-frame kinematic residual before the sampled residual rate, discrete drive bound, and continuous drive bound are accepted. The phase drift bound is recorded separately as configured_phase_drift_bound_units, then added to observed phase dispersion before the phase margin is accepted. The manifest now also names the Lean PhaseBudgetBounds.budgetCertificate predicate and phase_budget_certificate_discharges_phase_lock theorem for that phase budget. The phase_budget_discharged field must replay that theorem condition exactly before the combined proof obligation can discharge. This keeps FRC or MIF specialisations from hiding residual uncertainty inside the relative-velocity term or phase uncertainty inside replay dispersion.

The same manifest now records the Lean KinematicBounds.acceptanceCertificate predicate and acceptance_certificate_discharges_runtime_preconditions theorem. The acceptance_certificate_discharged field is recomputed from the spatial Gronwall margin, phase-budget discharge, and moving-frame equation replay certificate before the manifest hash is accepted.

pha_c_formal_obligation

Deterministic Lean proof-obligation manifests for PHA-C acceptance records.

The PHA-C acceptance chain is runtime evidence. The Lean kinematic proofs are formal evidence. This module binds the two surfaces by projecting a verified PHACAcceptanceRecord into fixed-point natural-number obligations that match SPOFormal.Kinematic.KinematicBounds. The resulting manifest remains review-only and non-actuating; it is a reproducible bridge for release review, MIF/FRC specialisation, and benchmark gating.

Classes

PHACKinematicProofObligation dataclass

PHACKinematicProofObligation(
    schema_version: str,
    evidence_kind: str,
    claim_boundary: str,
    acceptance_claim_boundary: str,
    execution_disabled: bool,
    actuating: bool,
    lean_module: str,
    lean_certificate_predicate: str,
    lean_theorem: str,
    continuous_lean_module: str,
    continuous_certificate_predicate: str,
    continuous_theorem: str,
    phase_lean_module: str,
    phase_certificate_predicate: str,
    phase_theorem: str,
    acceptance_certificate_predicate: str,
    acceptance_certificate_theorem: str,
    fixed_point_scale_m: float,
    fixed_point_scale_rad: float,
    fixed_point_time_scale_s: float,
    time_step_s: float,
    time_scale_units_per_second: int,
    time_step_units: int,
    horizon_time_units: int,
    initial_tolerance_units: int,
    lipschitz_step_gain_units: int,
    relative_velocity_rate_bound_units_per_second: int,
    relative_velocity_step_bound_units: int,
    configured_coupling_residual_step_bound_units: int,
    coupling_residual_rate_bound_units_per_second: int,
    coupling_residual_step_bound_units: int,
    continuous_drive_rate_bound_units_per_second: int,
    continuous_horizon_drive_bound_units: int,
    continuous_linear_budget_units: int,
    continuous_margin_units: int,
    drive_bound_units: int,
    merge_window_tolerance_units: int,
    horizon_steps: int,
    linear_budget_units: int,
    gronwall_budget_units: int,
    gronwall_budget_margin_units: int,
    gronwall_budget_trace_sha256: str,
    window_budget_margin_units: int,
    phase_tolerance_units: int,
    max_phase_dispersion_units: int,
    configured_phase_drift_bound_units: int,
    phase_budget_units: int,
    phase_margin_units: int,
    phase_budget_discharged: bool,
    acceptance_kinematic_equations_validated: bool,
    acceptance_kinematic_summary_replay_tolerance: float,
    acceptance_kinematic_summary_replay_tolerance_units: int,
    acceptance_kinematic_summary_replay_tolerance_limit_units: int,
    acceptance_replay_certificate_discharged: bool,
    acceptance_certificate_discharged: bool,
    observed_velocity_step_units: int,
    kinematic_residual_units: int,
    path_length_units: int,
    max_spatial_dispersion_units: int,
    continuous_envelope_discharged: bool,
    proof_obligations_discharged: bool,
    acceptance_sha256: str,
    timeline_sha256: str,
    record_sha256: str,
)

Review-only fixed-point obligations linked to the Lean kinematic proof.

Methods:

to_dict

to_dict() -> dict[str, bool | float | int | str]

Return a verified JSON-safe canonical representation.

Returns

dict[str, bool | float | int | str] The verified JSON-safe canonical representation.

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

def to_dict(self) -> dict[str, bool | float | int | str]:
    """Return a verified JSON-safe canonical representation.

    Returns
    -------
    dict[str, bool | float | int | str]
        The verified JSON-safe canonical representation.
    """
    return pha_c_kinematic_proof_obligation_to_dict(self)

Functions:

pha_c_kinematic_proof_obligation_to_dict

pha_c_kinematic_proof_obligation_to_dict(
    obligation: PHACKinematicProofObligation,
) -> dict[str, bool | float | int | str]

Return a verified canonical JSON-safe proof-obligation manifest.

Parameters

obligation : PHACKinematicProofObligation The PHA-C kinematic proof obligation to verify and serialise.

Returns

dict[str, bool | float | int | str] The verified canonical JSON-safe proof-obligation manifest.

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

def pha_c_kinematic_proof_obligation_to_dict(
    obligation: PHACKinematicProofObligation,
) -> dict[str, bool | float | int | str]:
    """Return a verified canonical JSON-safe proof-obligation manifest.

    Parameters
    ----------
    obligation : PHACKinematicProofObligation
        The PHA-C kinematic proof obligation to verify and serialise.

    Returns
    -------
    dict[str, bool | float | int | str]
        The verified canonical JSON-safe proof-obligation manifest.
    """
    verified = verify_pha_c_kinematic_proof_obligation(obligation)
    payload = _dict_without_record_hash(verified)
    payload["record_sha256"] = verified.record_sha256
    return payload

build_pha_c_kinematic_proof_obligation

build_pha_c_kinematic_proof_obligation(
    record: PHACAcceptanceRecord,
    *,
    fixed_point_scale_m: float = PHA_C_FORMAL_DEFAULT_SCALE_M,
    fixed_point_scale_rad: float = PHA_C_FORMAL_DEFAULT_SCALE_RAD,
    fixed_point_time_scale_s: float = PHA_C_FORMAL_DEFAULT_TIME_SCALE_S,
    relative_velocity_step_bound_m: float = 0.0,
    coupling_residual_step_bound_m: float = 0.0,
    phase_drift_bound_rad: float = 0.0,
    lipschitz_step_gain_units: int = 0,
) -> PHACKinematicProofObligation

Project a verified PHA-C acceptance record into Lean proof obligations.

The default obligation is a replay certificate: the maximum observed spatial dispersion is already measured over the accepted trajectory, so the Lean drive term only includes explicitly supplied future relative-velocity slack and the signed moving-frame residual. MIF/FRC specialisations can provide non-zero relative_velocity_step_bound_m and coupling_residual_step_bound_m, phase_drift_bound_rad, and lipschitz_step_gain_units values when they want a predictive finite-horizon Gronwall certificate instead of a replay-only envelope.

Parameters

record : PHACAcceptanceRecord The PHA-C record to operate on. fixed_point_scale_m : float Spatial fixed-point scale in metres. fixed_point_scale_rad : float Phase fixed-point scale in radians. fixed_point_time_scale_s : float Temporal fixed-point scale in seconds. relative_velocity_step_bound_m : float Per-step relative-velocity bound in metres. coupling_residual_step_bound_m : float Per-step coupling-residual bound in metres. phase_drift_bound_rad : float Per-step phase-drift bound in radians. lipschitz_step_gain_units : int Lipschitz step-gain bound in dimensionless integer units.

Returns

PHACKinematicProofObligation The Lean proof-obligation projection of the acceptance record.

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

def build_pha_c_kinematic_proof_obligation(
    record: PHACAcceptanceRecord,
    *,
    fixed_point_scale_m: float = PHA_C_FORMAL_DEFAULT_SCALE_M,
    fixed_point_scale_rad: float = PHA_C_FORMAL_DEFAULT_SCALE_RAD,
    fixed_point_time_scale_s: float = PHA_C_FORMAL_DEFAULT_TIME_SCALE_S,
    relative_velocity_step_bound_m: float = 0.0,
    coupling_residual_step_bound_m: float = 0.0,
    phase_drift_bound_rad: float = 0.0,
    lipschitz_step_gain_units: int = 0,
) -> PHACKinematicProofObligation:
    """Project a verified PHA-C acceptance record into Lean proof obligations.

    The default obligation is a replay certificate: the maximum observed
    spatial dispersion is already measured over the accepted trajectory, so the
    Lean drive term only includes explicitly supplied future relative-velocity
    slack and the signed moving-frame residual. MIF/FRC specialisations can
    provide non-zero ``relative_velocity_step_bound_m`` and
    ``coupling_residual_step_bound_m``, ``phase_drift_bound_rad``, and
    ``lipschitz_step_gain_units`` values when they want a predictive
    finite-horizon Gronwall certificate instead of a replay-only envelope.

    Parameters
    ----------
    record : PHACAcceptanceRecord
        The PHA-C record to operate on.
    fixed_point_scale_m : float
        Spatial fixed-point scale in metres.
    fixed_point_scale_rad : float
        Phase fixed-point scale in radians.
    fixed_point_time_scale_s : float
        Temporal fixed-point scale in seconds.
    relative_velocity_step_bound_m : float
        Per-step relative-velocity bound in metres.
    coupling_residual_step_bound_m : float
        Per-step coupling-residual bound in metres.
    phase_drift_bound_rad : float
        Per-step phase-drift bound in radians.
    lipschitz_step_gain_units : int
        Lipschitz step-gain bound in dimensionless integer units.

    Returns
    -------
    PHACKinematicProofObligation
        The Lean proof-obligation projection of the acceptance record.
    """
    verified_record = verify_pha_c_acceptance_record(record)
    scale_m = _validate_positive_scale(fixed_point_scale_m, name="fixed_point_scale_m")
    scale_rad = _validate_positive_scale(
        fixed_point_scale_rad,
        name="fixed_point_scale_rad",
    )
    time_scale = _validate_positive_scale(
        fixed_point_time_scale_s,
        name="fixed_point_time_scale_s",
    )
    time_step_s = _validate_positive_scale(verified_record.dt, name="time_step_s")
    time_scale_units_per_second = _validate_int(
        _ceil_positive_ratio_units(
            1.0,
            time_scale,
            name="time_scale_units_per_second",
        ),
        name="time_scale_units_per_second",
        minimum=1,
    )
    time_step_units = _validate_int(
        _ceil_positive_ratio_units(
            time_step_s,
            time_scale,
            name="time_step_units",
        ),
        name="time_step_units",
        minimum=1,
    )
    horizon_steps = _validate_int(
        verified_record.step_count,
        name="step_count",
        minimum=1,
    )
    horizon_time_units = horizon_steps * time_step_units
    gain_units = _validate_int(
        lipschitz_step_gain_units,
        name="lipschitz_step_gain_units",
        minimum=0,
    )
    raw_relative_velocity_units = _nonnegative_units(
        relative_velocity_step_bound_m,
        scale=scale_m,
        name="relative_velocity_step_bound_m",
    )
    configured_residual_m = _validate_nonnegative_scalar(
        coupling_residual_step_bound_m,
        name="coupling_residual_step_bound_m",
    )
    observed_residual_m = _validate_nonnegative_scalar(
        verified_record.kinematic_residual_max_m,
        name="kinematic_residual_max_m",
    )
    configured_residual_units = _nonnegative_units(
        configured_residual_m,
        scale=scale_m,
        name="coupling_residual_step_bound_m",
    )
    observed_residual_units = _nonnegative_units(
        observed_residual_m,
        scale=scale_m,
        name="kinematic_residual_max_m",
    )
    raw_residual_units = max(configured_residual_units, observed_residual_units)
    residual_bound_m = max(configured_residual_m, observed_residual_m)
    relative_velocity_rate_units = max(
        _nonnegative_units(
            relative_velocity_step_bound_m / time_step_s,
            scale=scale_m,
            name="relative_velocity_rate_bound_m_per_s",
        ),
        _ceil_div_units(
            raw_relative_velocity_units * time_scale_units_per_second,
            time_step_units,
            name="relative_velocity_rate_bound_units_per_second",
        ),
    )
    residual_rate_units = max(
        _nonnegative_units(
            residual_bound_m / time_step_s,
            scale=scale_m,
            name="coupling_residual_rate_bound_m_per_s",
        ),
        _ceil_div_units(
            raw_residual_units * time_scale_units_per_second,
            time_step_units,
            name="coupling_residual_rate_bound_units_per_second",
        ),
    )
    relative_velocity_units = _ceil_div_units(
        relative_velocity_rate_units * time_step_units,
        time_scale_units_per_second,
        name="relative_velocity_step_bound_units",
    )
    residual_units = _ceil_div_units(
        residual_rate_units * time_step_units,
        time_scale_units_per_second,
        name="coupling_residual_step_bound_units",
    )
    drive_units = relative_velocity_units + residual_units
    initial_units = _nonnegative_units(
        verified_record.max_spatial_dispersion_m,
        scale=scale_m,
        name="max_spatial_dispersion_m",
    )
    merge_tolerance_units = _nonnegative_units(
        verified_record.spatial_tol_m,
        scale=scale_m,
        name="spatial_tol_m",
    )
    continuous_drive_rate_units = relative_velocity_rate_units + residual_rate_units
    continuous_horizon_drive_units = _ceil_div_units(
        continuous_drive_rate_units * horizon_time_units,
        time_scale_units_per_second,
        name="continuous_horizon_drive_bound_units",
    )
    continuous_linear_budget_units = initial_units + continuous_horizon_drive_units
    continuous_margin_units = merge_tolerance_units - continuous_linear_budget_units
    continuous_envelope_discharged = continuous_margin_units >= 0
    linear_budget_units = initial_units + horizon_steps * drive_units
    gronwall_trace_units = _gronwall_budget_trace(
        initial_tolerance_units=initial_units,
        lipschitz_step_gain_units=gain_units,
        drive_bound_units=drive_units,
        horizon_steps=horizon_steps,
    )
    gronwall_budget_units = gronwall_trace_units[-1]
    gronwall_budget_margin_units = merge_tolerance_units - gronwall_budget_units
    gronwall_trace_sha256 = _gronwall_budget_trace_sha256(
        trace_units=gronwall_trace_units,
        horizon_steps=horizon_steps,
    )
    window_margin_units = gronwall_budget_margin_units
    phase_tolerance_units = _nonnegative_units(
        verified_record.phase_tol_rad,
        scale=scale_rad,
        name="phase_tol_rad",
    )
    phase_dispersion_units = _nonnegative_units(
        verified_record.max_phase_dispersion_rad,
        scale=scale_rad,
        name="max_phase_dispersion_rad",
    )
    configured_phase_drift_units = _nonnegative_units(
        phase_drift_bound_rad,
        scale=scale_rad,
        name="phase_drift_bound_rad",
    )
    phase_budget_units = phase_dispersion_units + configured_phase_drift_units
    phase_margin_units = phase_tolerance_units - phase_budget_units
    phase_budget_discharged = phase_margin_units >= 0
    acceptance_tolerance_units = _nonnegative_units(
        verified_record.kinematic_summary_replay_tolerance,
        scale=scale_m,
        name="acceptance_kinematic_summary_replay_tolerance",
    )
    acceptance_tolerance_limit_units = _nonnegative_units(
        PHA_C_ACCEPTANCE_KINEMATIC_SUMMARY_REPLAY_TOLERANCE,
        scale=scale_m,
        name="acceptance_kinematic_summary_replay_tolerance_limit",
    )
    acceptance_replay_certificate_discharged = (
        verified_record.kinematic_equations_validated
        and acceptance_tolerance_units <= acceptance_tolerance_limit_units
    )
    acceptance_certificate_discharged = (
        gronwall_budget_margin_units >= 0
        and phase_budget_discharged
        and acceptance_replay_certificate_discharged
    )
    observed_velocity_step_units = _nonnegative_units(
        verified_record.max_abs_velocity_m_per_s * verified_record.dt,
        scale=scale_m,
        name="observed_velocity_step_m",
    )
    path_length_units = _nonnegative_units(
        verified_record.path_length_max_m,
        scale=scale_m,
        name="path_length_max_m",
    )
    discharged = (
        window_margin_units >= 0
        and phase_margin_units >= 0
        and phase_budget_discharged
        and acceptance_certificate_discharged
        and continuous_envelope_discharged
        and verified_record.execution_disabled
        and not verified_record.actuating
        and verified_record.claim_boundary == PHA_C_ACCEPTANCE_CLAIM_BOUNDARY
    )
    payload_without_hash: dict[str, Any] = {
        "schema_version": PHA_C_FORMAL_OBLIGATION_SCHEMA,
        "evidence_kind": PHA_C_FORMAL_OBLIGATION_EVIDENCE_KIND,
        "claim_boundary": PHA_C_FORMAL_OBLIGATION_CLAIM_BOUNDARY,
        "acceptance_claim_boundary": verified_record.claim_boundary,
        "execution_disabled": True,
        "actuating": False,
        "lean_module": PHA_C_FORMAL_LEAN_MODULE,
        "lean_certificate_predicate": PHA_C_FORMAL_CERTIFICATE_PREDICATE,
        "lean_theorem": PHA_C_FORMAL_CERTIFICATE_THEOREM,
        "continuous_lean_module": PHA_C_FORMAL_CONTINUOUS_LEAN_MODULE,
        "continuous_certificate_predicate": (
            PHA_C_FORMAL_CONTINUOUS_CERTIFICATE_PREDICATE
        ),
        "continuous_theorem": PHA_C_FORMAL_CONTINUOUS_CERTIFICATE_THEOREM,
        "phase_lean_module": PHA_C_FORMAL_PHASE_LEAN_MODULE,
        "phase_certificate_predicate": PHA_C_FORMAL_PHASE_CERTIFICATE_PREDICATE,
        "phase_theorem": PHA_C_FORMAL_PHASE_CERTIFICATE_THEOREM,
        "acceptance_certificate_predicate": (
            PHA_C_FORMAL_ACCEPTANCE_CERTIFICATE_PREDICATE
        ),
        "acceptance_certificate_theorem": (PHA_C_FORMAL_ACCEPTANCE_CERTIFICATE_THEOREM),
        "fixed_point_scale_m": scale_m,
        "fixed_point_scale_rad": scale_rad,
        "fixed_point_time_scale_s": time_scale,
        "time_step_s": time_step_s,
        "time_scale_units_per_second": time_scale_units_per_second,
        "time_step_units": time_step_units,
        "horizon_time_units": horizon_time_units,
        "initial_tolerance_units": initial_units,
        "lipschitz_step_gain_units": gain_units,
        "relative_velocity_rate_bound_units_per_second": relative_velocity_rate_units,
        "relative_velocity_step_bound_units": relative_velocity_units,
        "configured_coupling_residual_step_bound_units": configured_residual_units,
        "coupling_residual_rate_bound_units_per_second": residual_rate_units,
        "coupling_residual_step_bound_units": residual_units,
        "continuous_drive_rate_bound_units_per_second": continuous_drive_rate_units,
        "continuous_horizon_drive_bound_units": continuous_horizon_drive_units,
        "continuous_linear_budget_units": continuous_linear_budget_units,
        "continuous_margin_units": continuous_margin_units,
        "drive_bound_units": drive_units,
        "merge_window_tolerance_units": merge_tolerance_units,
        "horizon_steps": horizon_steps,
        "linear_budget_units": linear_budget_units,
        "gronwall_budget_units": gronwall_budget_units,
        "gronwall_budget_margin_units": gronwall_budget_margin_units,
        "gronwall_budget_trace_sha256": gronwall_trace_sha256,
        "window_budget_margin_units": window_margin_units,
        "phase_tolerance_units": phase_tolerance_units,
        "max_phase_dispersion_units": phase_dispersion_units,
        "configured_phase_drift_bound_units": configured_phase_drift_units,
        "phase_budget_units": phase_budget_units,
        "phase_margin_units": phase_margin_units,
        "phase_budget_discharged": phase_budget_discharged,
        "acceptance_kinematic_equations_validated": (
            verified_record.kinematic_equations_validated
        ),
        "acceptance_kinematic_summary_replay_tolerance": (
            verified_record.kinematic_summary_replay_tolerance
        ),
        "acceptance_kinematic_summary_replay_tolerance_units": (
            acceptance_tolerance_units
        ),
        "acceptance_kinematic_summary_replay_tolerance_limit_units": (
            acceptance_tolerance_limit_units
        ),
        "acceptance_replay_certificate_discharged": (
            acceptance_replay_certificate_discharged
        ),
        "acceptance_certificate_discharged": acceptance_certificate_discharged,
        "observed_velocity_step_units": observed_velocity_step_units,
        "kinematic_residual_units": observed_residual_units,
        "path_length_units": path_length_units,
        "max_spatial_dispersion_units": initial_units,
        "continuous_envelope_discharged": continuous_envelope_discharged,
        "proof_obligations_discharged": discharged,
        "acceptance_sha256": verified_record.acceptance_sha256,
        "timeline_sha256": verified_record.timeline_sha256,
    }
    return PHACKinematicProofObligation(
        **payload_without_hash,
        record_sha256=_sha256_json(payload_without_hash),
    )

verify_pha_c_kinematic_proof_obligation

verify_pha_c_kinematic_proof_obligation(
    obligation: PHACKinematicProofObligation,
) -> PHACKinematicProofObligation

Validate a PHA-C Lean proof-obligation manifest fail-closed.

Parameters

obligation : PHACKinematicProofObligation The PHA-C kinematic proof obligation to operate on.

Returns

PHACKinematicProofObligation The same obligation after fail-closed validation.

Raises

TypeError If the manifest has the wrong type. ValueError If the manifest fails validation.

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

def verify_pha_c_kinematic_proof_obligation(
    obligation: PHACKinematicProofObligation,
) -> PHACKinematicProofObligation:
    """Validate a PHA-C Lean proof-obligation manifest fail-closed.

    Parameters
    ----------
    obligation : PHACKinematicProofObligation
        The PHA-C kinematic proof obligation to operate on.

    Returns
    -------
    PHACKinematicProofObligation
        The same obligation after fail-closed validation.

    Raises
    ------
    TypeError
        If the manifest has the wrong type.
    ValueError
        If the manifest fails validation.
    """
    if not isinstance(obligation, PHACKinematicProofObligation):
        raise TypeError("obligation must be a PHACKinematicProofObligation")
    exact_strings = {
        "schema_version": PHA_C_FORMAL_OBLIGATION_SCHEMA,
        "evidence_kind": PHA_C_FORMAL_OBLIGATION_EVIDENCE_KIND,
        "claim_boundary": PHA_C_FORMAL_OBLIGATION_CLAIM_BOUNDARY,
        "acceptance_claim_boundary": PHA_C_ACCEPTANCE_CLAIM_BOUNDARY,
        "lean_module": PHA_C_FORMAL_LEAN_MODULE,
        "lean_certificate_predicate": PHA_C_FORMAL_CERTIFICATE_PREDICATE,
        "lean_theorem": PHA_C_FORMAL_CERTIFICATE_THEOREM,
        "continuous_lean_module": PHA_C_FORMAL_CONTINUOUS_LEAN_MODULE,
        "continuous_certificate_predicate": (
            PHA_C_FORMAL_CONTINUOUS_CERTIFICATE_PREDICATE
        ),
        "continuous_theorem": PHA_C_FORMAL_CONTINUOUS_CERTIFICATE_THEOREM,
        "phase_lean_module": PHA_C_FORMAL_PHASE_LEAN_MODULE,
        "phase_certificate_predicate": PHA_C_FORMAL_PHASE_CERTIFICATE_PREDICATE,
        "phase_theorem": PHA_C_FORMAL_PHASE_CERTIFICATE_THEOREM,
        "acceptance_certificate_predicate": (
            PHA_C_FORMAL_ACCEPTANCE_CERTIFICATE_PREDICATE
        ),
        "acceptance_certificate_theorem": (PHA_C_FORMAL_ACCEPTANCE_CERTIFICATE_THEOREM),
    }
    for field, expected in exact_strings.items():
        got = getattr(obligation, field)
        if got != expected:
            raise ValueError(f"{field} must be {expected!r}")

    if not _validate_bool(
        obligation.execution_disabled,
        name="execution_disabled",
    ):
        raise ValueError("execution_disabled must be true")
    if _validate_bool(obligation.actuating, name="actuating"):
        raise ValueError("actuating must be false")
    _validate_positive_scale(
        obligation.fixed_point_scale_m,
        name="fixed_point_scale_m",
    )
    _validate_positive_scale(
        obligation.fixed_point_scale_rad,
        name="fixed_point_scale_rad",
    )
    _validate_positive_scale(
        obligation.fixed_point_time_scale_s,
        name="fixed_point_time_scale_s",
    )
    _validate_positive_scale(obligation.time_step_s, name="time_step_s")
    nat_fields = (
        "time_scale_units_per_second",
        "time_step_units",
        "horizon_time_units",
        "initial_tolerance_units",
        "lipschitz_step_gain_units",
        "relative_velocity_rate_bound_units_per_second",
        "relative_velocity_step_bound_units",
        "configured_coupling_residual_step_bound_units",
        "coupling_residual_rate_bound_units_per_second",
        "coupling_residual_step_bound_units",
        "continuous_drive_rate_bound_units_per_second",
        "continuous_horizon_drive_bound_units",
        "continuous_linear_budget_units",
        "drive_bound_units",
        "merge_window_tolerance_units",
        "horizon_steps",
        "linear_budget_units",
        "gronwall_budget_units",
        "phase_tolerance_units",
        "max_phase_dispersion_units",
        "configured_phase_drift_bound_units",
        "phase_budget_units",
        "acceptance_kinematic_summary_replay_tolerance_units",
        "acceptance_kinematic_summary_replay_tolerance_limit_units",
        "observed_velocity_step_units",
        "kinematic_residual_units",
        "path_length_units",
        "max_spatial_dispersion_units",
    )
    for field in nat_fields:
        _validate_int(getattr(obligation, field), name=field, minimum=0)
    _validate_int(
        obligation.gronwall_budget_margin_units,
        name="gronwall_budget_margin_units",
        minimum=-(10**18),
    )
    _validate_int(
        obligation.window_budget_margin_units,
        name="window_budget_margin_units",
        minimum=-(10**18),
    )
    _validate_int(
        obligation.continuous_margin_units,
        name="continuous_margin_units",
        minimum=-(10**18),
    )
    _validate_sha256_hex(
        obligation.gronwall_budget_trace_sha256,
        name="gronwall_budget_trace_sha256",
    )
    _validate_int(
        obligation.phase_margin_units,
        name="phase_margin_units",
        minimum=-(10**18),
    )
    _validate_bool(
        obligation.phase_budget_discharged,
        name="phase_budget_discharged",
    )
    if not _validate_bool(
        obligation.acceptance_kinematic_equations_validated,
        name="acceptance_kinematic_equations_validated",
    ):
        raise ValueError("acceptance_kinematic_equations_validated must be true")
    kinematic_replay_tolerance = _validate_positive_scale(
        obligation.acceptance_kinematic_summary_replay_tolerance,
        name="acceptance_kinematic_summary_replay_tolerance",
    )
    if (
        kinematic_replay_tolerance
        != PHA_C_ACCEPTANCE_KINEMATIC_SUMMARY_REPLAY_TOLERANCE
    ):
        raise ValueError(
            "acceptance_kinematic_summary_replay_tolerance must match "
            "the acceptance constant"
        )
    expected_acceptance_tolerance_units = _nonnegative_units(
        kinematic_replay_tolerance,
        scale=obligation.fixed_point_scale_m,
        name="acceptance_kinematic_summary_replay_tolerance",
    )
    if (
        obligation.acceptance_kinematic_summary_replay_tolerance_units
        != expected_acceptance_tolerance_units
    ):
        raise ValueError(
            "acceptance_kinematic_summary_replay_tolerance_units must replay",
        )
    expected_acceptance_tolerance_limit_units = _nonnegative_units(
        PHA_C_ACCEPTANCE_KINEMATIC_SUMMARY_REPLAY_TOLERANCE,
        scale=obligation.fixed_point_scale_m,
        name="acceptance_kinematic_summary_replay_tolerance_limit",
    )
    if (
        obligation.acceptance_kinematic_summary_replay_tolerance_limit_units
        != expected_acceptance_tolerance_limit_units
    ):
        raise ValueError(
            "acceptance_kinematic_summary_replay_tolerance_limit_units must replay",
        )
    _validate_bool(
        obligation.acceptance_replay_certificate_discharged,
        name="acceptance_replay_certificate_discharged",
    )
    _validate_bool(
        obligation.acceptance_certificate_discharged,
        name="acceptance_certificate_discharged",
    )
    _validate_bool(
        obligation.continuous_envelope_discharged,
        name="continuous_envelope_discharged",
    )
    _validate_bool(
        obligation.proof_obligations_discharged,
        name="proof_obligations_discharged",
    )
    _validate_sha256_hex(obligation.acceptance_sha256, name="acceptance_sha256")
    _validate_sha256_hex(obligation.timeline_sha256, name="timeline_sha256")
    _validate_sha256_hex(obligation.record_sha256, name="record_sha256")

    expected_time_scale_units = _ceil_positive_ratio_units(
        1.0,
        obligation.fixed_point_time_scale_s,
        name="time_scale_units_per_second",
    )
    if obligation.time_scale_units_per_second != expected_time_scale_units:
        raise ValueError("time_scale_units_per_second must match time scale")
    expected_time_step_units = _ceil_positive_ratio_units(
        obligation.time_step_s,
        obligation.fixed_point_time_scale_s,
        name="time_step_units",
    )
    if obligation.time_step_units != expected_time_step_units:
        raise ValueError("time_step_units must match time step")
    expected_horizon_time_units = obligation.horizon_steps * obligation.time_step_units
    if obligation.horizon_time_units != expected_horizon_time_units:
        raise ValueError("horizon_time_units must match horizon and time step")

    expected_relative_velocity_units = _ceil_div_units(
        obligation.relative_velocity_rate_bound_units_per_second
        * obligation.time_step_units,
        obligation.time_scale_units_per_second,
        name="relative_velocity_step_bound_units",
    )
    if (
        obligation.relative_velocity_step_bound_units
        != expected_relative_velocity_units
    ):
        raise ValueError(
            "relative_velocity_step_bound_units must match sampled rate bound",
        )
    expected_residual_units = _ceil_div_units(
        obligation.coupling_residual_rate_bound_units_per_second
        * obligation.time_step_units,
        obligation.time_scale_units_per_second,
        name="coupling_residual_step_bound_units",
    )
    if obligation.coupling_residual_step_bound_units != expected_residual_units:
        raise ValueError(
            "coupling_residual_step_bound_units must match sampled rate bound",
        )

    expected_continuous_drive_rate = (
        obligation.relative_velocity_rate_bound_units_per_second
        + obligation.coupling_residual_rate_bound_units_per_second
    )
    if (
        obligation.continuous_drive_rate_bound_units_per_second
        != expected_continuous_drive_rate
    ):
        raise ValueError(
            "continuous_drive_rate_bound_units_per_second must equal rate sum",
        )
    expected_continuous_horizon_drive = _ceil_div_units(
        expected_continuous_drive_rate * obligation.horizon_time_units,
        obligation.time_scale_units_per_second,
        name="continuous_horizon_drive_bound_units",
    )
    if (
        obligation.continuous_horizon_drive_bound_units
        != expected_continuous_horizon_drive
    ):
        raise ValueError(
            "continuous_horizon_drive_bound_units must match sampled horizon rate",
        )
    expected_continuous_linear_budget = (
        obligation.initial_tolerance_units + expected_continuous_horizon_drive
    )
    if obligation.continuous_linear_budget_units != (expected_continuous_linear_budget):
        raise ValueError(
            "continuous_linear_budget_units must match continuous envelope",
        )
    expected_continuous_margin = (
        obligation.merge_window_tolerance_units - expected_continuous_linear_budget
    )
    if obligation.continuous_margin_units != expected_continuous_margin:
        raise ValueError("continuous_margin_units must match continuous envelope")
    expected_continuous_discharged = expected_continuous_margin >= 0
    if obligation.continuous_envelope_discharged != expected_continuous_discharged:
        raise ValueError(
            "continuous_envelope_discharged does not match certificate math",
        )

    expected_drive = (
        obligation.relative_velocity_step_bound_units
        + obligation.coupling_residual_step_bound_units
    )
    if obligation.drive_bound_units != expected_drive:
        raise ValueError("drive_bound_units must equal relative velocity plus residual")
    if obligation.kinematic_residual_units > (
        obligation.coupling_residual_step_bound_units
    ):
        raise ValueError("kinematic_residual_units must fit sampled residual bound")
    if obligation.configured_coupling_residual_step_bound_units > (
        obligation.coupling_residual_step_bound_units
    ):
        raise ValueError(
            "configured_coupling_residual_step_bound_units must fit residual bound",
        )
    if obligation.max_spatial_dispersion_units != obligation.initial_tolerance_units:
        raise ValueError("max_spatial_dispersion_units must mirror initial tolerance")
    expected_budget = (
        obligation.initial_tolerance_units
        + obligation.horizon_steps * obligation.drive_bound_units
    )
    if obligation.linear_budget_units != expected_budget:
        raise ValueError("linear_budget_units must match the Lean linear budget")
    expected_gronwall_trace = _gronwall_budget_trace(
        initial_tolerance_units=obligation.initial_tolerance_units,
        lipschitz_step_gain_units=obligation.lipschitz_step_gain_units,
        drive_bound_units=obligation.drive_bound_units,
        horizon_steps=obligation.horizon_steps,
    )
    expected_gronwall_budget = expected_gronwall_trace[-1]
    if obligation.gronwall_budget_units != expected_gronwall_budget:
        raise ValueError("gronwall_budget_units must match the Lean Gronwall budget")
    expected_gronwall_margin = (
        obligation.merge_window_tolerance_units - expected_gronwall_budget
    )
    if obligation.gronwall_budget_margin_units != expected_gronwall_margin:
        raise ValueError("gronwall_budget_margin_units must match the Lean margin")
    expected_trace_hash = _gronwall_budget_trace_sha256(
        trace_units=expected_gronwall_trace,
        horizon_steps=obligation.horizon_steps,
    )
    if obligation.gronwall_budget_trace_sha256 != expected_trace_hash:
        raise ValueError("gronwall_budget_trace_sha256 does not replay")
    expected_margin = expected_gronwall_margin
    if obligation.window_budget_margin_units != expected_margin:
        raise ValueError("window_budget_margin_units must match the Lean margin")
    expected_phase_budget = (
        obligation.max_phase_dispersion_units
        + obligation.configured_phase_drift_bound_units
    )
    if obligation.phase_budget_units != expected_phase_budget:
        raise ValueError(
            "phase_budget_units must match dispersion plus configured drift",
        )
    expected_phase_margin = obligation.phase_tolerance_units - expected_phase_budget
    if obligation.phase_margin_units != expected_phase_margin:
        raise ValueError(
            "phase_margin_units must match phase tolerance minus phase budget",
        )
    expected_phase_budget_discharged = expected_phase_margin >= 0
    if obligation.phase_budget_discharged != expected_phase_budget_discharged:
        raise ValueError(
            "phase_budget_discharged does not match phase certificate math",
        )
    expected_acceptance_replay_discharged = (
        obligation.acceptance_kinematic_equations_validated
        and obligation.acceptance_kinematic_summary_replay_tolerance_units
        <= obligation.acceptance_kinematic_summary_replay_tolerance_limit_units
    )
    if (
        obligation.acceptance_replay_certificate_discharged
        != expected_acceptance_replay_discharged
    ):
        raise ValueError(
            "acceptance_replay_certificate_discharged does not match Lean replay",
        )
    expected_acceptance_certificate_discharged = (
        expected_margin >= 0
        and expected_phase_budget_discharged
        and expected_acceptance_replay_discharged
    )
    if (
        obligation.acceptance_certificate_discharged
        != expected_acceptance_certificate_discharged
    ):
        raise ValueError(
            "acceptance_certificate_discharged does not match Lean certificate",
        )
    expected_discharged = (
        expected_margin >= 0
        and expected_phase_budget_discharged
        and expected_acceptance_certificate_discharged
        and expected_continuous_discharged
        and obligation.execution_disabled
        and not obligation.actuating
    )
    if obligation.proof_obligations_discharged != expected_discharged:
        raise ValueError("proof_obligations_discharged does not match certificate math")
    expected_hash = _sha256_json(_dict_without_record_hash(obligation))
    if obligation.record_sha256 != expected_hash:
        raise ValueError("record_sha256 does not match canonical obligation payload")
    return obligation