Oscillators

Phase extraction from raw signals via canonical Physical (P), Informational (I), and Symbolic (S) channels, plus named extension channels. The P/I/S decomposition is the default abstraction that makes SPO domain-agnostic, but deployments are not limited to three channels: any signal that exhibits periodic or quasi-periodic behaviour can map onto one or more named channels.

Pipeline position

Raw signals ──→ PhysicalExtractor  ──→ PhaseState(θ, ω, quality)
Event streams ──→ InformationalExtractor ──→ PhaseState(θ, ω, quality)
Sequences ──→ SymbolicExtractor ──→ PhaseState(θ, ω, quality)
                                              │
                                              ↓
                                    PhaseQualityScorer
                                              │
                                    ┌─────────┼──────────┐
                                    ↓         ↓          ↓
                              θ array    ω array    quality mask
                                    │         │          │
                                    ↓         ↓          ↓
                           UPDEEngine.step(phases, omegas, knm * mask, ...)

Oscillators are the input adapters of the SPO pipeline. They convert raw domain signals into the (θ, ω) vectors that the engine requires. Quality scores gate which oscillators participate in coupling.

---

The PIS Model

Every domain signal decomposes into one or more oscillator channels. The canonical channels are:

Channel	Signal type	Extraction method	Example domains
P (Physical)	Continuous waveforms	Hilbert transform	EEG, ECG, vibration, voltage, plasma
I (Informational)	Event streams, rates	Inter-event interval	Network traffic, API calls, manufacturing
S (Symbolic)	Categorical sequences	Ring mapping	Protocols, language, music, genetics

Not every domain uses all three canonical channels. A pure physics domain (tokamak plasma) might use only P. A pure IT domain (microservices) might use only I. Larger deployments can add named extension channels such as thermal, market_sentiment, or operator_intent while preserving the same PhaseState contract. The binding specification declares which channels are active.

---

Phase State

PhaseState (dataclass)

Field	Type	Range	Description
theta	float	[0, 2π)	Phase angle
omega	float	R	Instantaneous frequency (rad/s)
amplitude	float	≥ 0	Signal strength (SNR proxy)
quality	float	[0, 1]	Extraction confidence
channel	str	Identifier	Binding channel (P, I, S, or named extension)
node_id	str	—	Unique oscillator identifier

Quality scores gate downstream processing: low-quality oscillators are downweighted in coupling and excluded from regime classification.

---

Extractor Interface

All channel extractors implement the PhaseExtractor abstract base class:

from numpy.typing import NDArray
import numpy as np

FloatArray = NDArray[np.float64]

class PhaseExtractor(ABC):
    @abstractmethod
    def extract(self, signal: FloatArray, sample_rate: float) -> list[PhaseState]: ...

    @abstractmethod
    def quality_score(self, phase_states: list[PhaseState]) -> float: ...

The extract method receives a raw signal window and sample rate, and returns one or more PhaseState objects. The quality_score method computes an aggregate quality for the extraction.

---

Physical Extraction (P)

PhysicalExtractor

PhysicalExtractor(node_id: str = "phys_0")

Uses the analytic signal (Hilbert transform) to decompose a real-valued waveform into instantaneous phase and amplitude:

z(t) = x(t) + i H[x(t)]
θ(t) = arg(z(t)),  A(t) = |z(t)|
ω = 2π × median(instantaneous frequency)

Quality metric

_envelope_quality(signal, analytic) returns quality based on the coefficient of variation (CV) of the analytic signal envelope:

quality = clip(1.0 - CV(|z(t)|), 0, 1)

Clean sinusoids have near-constant envelope (CV ≈ 0, quality ≈ 1.0). Noisy signals have variable envelope (high CV, low quality).

Validation

• Rejects empty signals, single-sample signals, and 2-D arrays with ValueError("1-D with >= 2 samples")
• Returns channel = "P", node_id from constructor

Rust acceleration

When spo_kernel is importable, uses spo_kernel.physical_extract() for the core computation. Python fallback uses scipy Hilbert transform. Parity verified in tests/test_oscillator_physical.py::test_rust_python_parity.

Performance: extract(1s @ 1kHz) < 5 ms.

physical

Physical-channel phase extraction from continuous numeric waveforms.

PhysicalExtractor validates finite one-dimensional real signals and positive sample rates, then derives instantaneous phase, angular frequency, amplitude, and envelope-quality metadata via the Hilbert transform. Optional Rust acceleration preserves the same PhaseState contract as the NumPy/SciPy path.

Classes

PhysicalExtractor

PhysicalExtractor(
    node_id: str = "phys_0",
    *,
    band: tuple[float, float]
    | Sequence[float]
    | None = None,
    filter_order: int = 4,
    edge_trim: int | None = None,
)

Bases: PhaseExtractor

Extracts instantaneous phase from continuous waveforms via Hilbert transform.

By default the extractor Hilbert-transforms the raw broadband signal and reports the trailing (endpoint) instantaneous phase — the historical behaviour, retained bit-for-bit so existing bindings and sealed evidence are unchanged. Two optional, opt-in refinements are available for callers who need a cleaner estimate:

• band=(low_hz, high_hz) applies a zero-phase Butterworth band-pass (scipy.signal.filtfilt) before the Hilbert transform, isolating the phase of interest instead of mixing all spectral content.
• edge_trim (or, when a band is set, an automatic filter-transient trim) discards edge samples where the FFT-Hilbert and filtfilt transients are worst, so the reported phase is taken from the last reliable interior sample rather than the artefact-prone endpoint.

Both refinements are applied identically before the NumPy and Rust paths, so the accelerated kernel needs no change and stays bit-parity with the reference.

Source code in src/scpn_phase_orchestrator/oscillators/physical.py

def __init__(
    self,
    node_id: str = "phys_0",
    *,
    band: tuple[float, float] | Sequence[float] | None = None,
    filter_order: int = 4,
    edge_trim: int | None = None,
):
    self._node_id = _validate_node_id(node_id)
    self._band = _validate_band(band)
    self._filter_order = _validate_filter_order(filter_order)
    self._edge_trim = _validate_edge_trim(edge_trim)

Methods:

extract

extract(
    signal: FloatArray, sample_rate: float
) -> list[PhaseState]

Extract instantaneous phase from a 1-D waveform via Hilbert transform.

Parameters

signal : FloatArray Input signal, shape (T,). sample_rate : float Sampling rate in Hz.

Returns

list[PhaseState] Instantaneous phase from a 1-D waveform via Hilbert transform.

Source code in src/scpn_phase_orchestrator/oscillators/physical.py

def extract(self, signal: FloatArray, sample_rate: float) -> list[PhaseState]:
    """Extract instantaneous phase from a 1-D waveform via Hilbert transform.

    Parameters
    ----------
    signal : FloatArray
        Input signal, shape ``(T,)``.
    sample_rate : float
        Sampling rate in Hz.

    Returns
    -------
    list[PhaseState]
        Instantaneous phase from a 1-D waveform via Hilbert transform.
    """
    signal = _validate_signal(signal)
    sample_rate = _validate_sample_rate(sample_rate)
    from scipy.signal import hilbert  # noqa: PLC0415

    filtered = self._apply_bandpass(signal, sample_rate)
    analytic = hilbert(filtered)

    trim = self._resolve_edge_trim(analytic.shape[0])
    if trim:
        filtered = filtered[trim : filtered.shape[0] - trim]
        analytic = analytic[trim : analytic.shape[0] - trim]

    if _rust_physical_extract is not None:
        try:
            theta, omega, amplitude, quality = _rust_physical_extract(
                np.ascontiguousarray(np.real(analytic)),
                np.ascontiguousarray(np.imag(analytic)),
                sample_rate,
            )
        except Exception:
            theta, omega, amplitude, quality = self._python_extract(
                filtered, analytic, sample_rate
            )
    else:
        theta, omega, amplitude, quality = self._python_extract(
            filtered, analytic, sample_rate
        )

    return [
        PhaseState(
            theta=theta,
            omega=omega,
            amplitude=amplitude,
            quality=quality,
            channel="P",
            node_id=self._node_id,
        )
    ]

quality_score

quality_score(phase_states: list[PhaseState]) -> float

Mean extraction quality across phase states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states.

Returns

float Mean extraction quality across phase states.

Source code in src/scpn_phase_orchestrator/oscillators/physical.py

def quality_score(self, phase_states: list[PhaseState]) -> float:
    """Mean extraction quality across phase states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.

    Returns
    -------
    float
        Mean extraction quality across phase states.
    """
    if not phase_states:
        return 0.0
    return float(np.mean([ps.quality for ps in phase_states]))

---

Informational Extraction (I)

InformationalExtractor

InformationalExtractor(node_id: str = "info_0")

Converts event timestamps into phase oscillators:

1. Compute inter-event intervals: τ_k = t_k - t_{k-1}
2. Median frequency: f = 1 / median(τ)
3. Angular frequency: ω = 2πf
4. Phase: θ = (2πf × total_duration) mod 2π
5. Amplitude: mean instantaneous frequency
6. Quality: 1/(1 + CV(τ)) where CV = std(τ)/mean(τ)

Edge cases

Input	Result
Single timestamp	θ=0, ω=0, quality=0
Identical timestamps	θ=0, ω=0, quality=0
Two timestamps	Valid extraction from one interval
Regular events	quality ≈ 1.0
Irregular events	quality < 0.9

Performance: extract(100 timestamps) < 500 μs.

informational

Informational-channel phase extraction from event timestamp streams.

InformationalExtractor treats sorted event timestamps as a cadence signal, deriving phase from inter-event intervals and reporting interval-regularity quality. Non-numeric, boolean, complex, non-finite, or unsorted timestamp inputs fail before they can seed runtime phase state.

Classes

InformationalExtractor

InformationalExtractor(node_id: str = 'info_0')

Bases: PhaseExtractor

Extracts phase from event timestamps (spike trains, discrete events).

Converts inter-event intervals to instantaneous frequency, then derives phase via cumulative integral of frequency.

Source code in src/scpn_phase_orchestrator/oscillators/informational.py

def __init__(self, node_id: str = "info_0"):
    self._node_id = _validate_node_id(node_id)

Methods:

extract

extract(
    signal: FloatArray, sample_rate: float
) -> list[PhaseState]

Extract phase states from event timestamps.

Parameters

signal : FloatArray 1-D array of event timestamps in seconds (sorted ascending). sample_rate : float not used for timestamps but kept for interface consistency.

Returns

list[PhaseState] The result.

Raises

ValueError If the inputs are invalid or inconsistent.

Source code in src/scpn_phase_orchestrator/oscillators/informational.py

def extract(self, signal: FloatArray, sample_rate: float) -> list[PhaseState]:
    """Extract phase states from event timestamps.

    Parameters
    ----------
    signal : FloatArray
        1-D array of event timestamps in seconds (sorted ascending).
    sample_rate : float
        not used for timestamps but kept for interface consistency.

    Returns
    -------
    list[PhaseState]
        The result.

    Raises
    ------
    ValueError
        If the inputs are invalid or inconsistent.
    """
    signal = _validate_signal(signal)
    raw_intervals = np.diff(signal)
    if np.any(raw_intervals < 0.0):
        raise ValueError("signal timestamps must be sorted ascending")
    if len(signal) < 2:
        return [
            PhaseState(
                theta=0.0,
                omega=0.0,
                amplitude=0.0,
                quality=0.0,
                channel="I",
                node_id=self._node_id,
            )
        ]

    intervals = raw_intervals[raw_intervals > 0]
    if len(intervals) == 0:
        return [
            PhaseState(
                theta=0.0,
                omega=0.0,
                amplitude=0.0,
                quality=0.0,
                channel="I",
                node_id=self._node_id,
            )
        ]

    if _rust_event_phase is not None:
        try:
            theta, omega_median, quality = _rust_event_phase(signal)
            inst_freq = 1.0 / intervals  # Hz (amplitude is mean frequency)
            amplitude = float(np.mean(inst_freq))
            return [
                PhaseState(
                    theta=float(theta),
                    omega=float(omega_median),
                    amplitude=amplitude,
                    quality=float(quality),
                    channel="I",
                    node_id=self._node_id,
                )
            ]
        except Exception:
            # Preserve validated Python semantics on optional-kernel failures.
            quality = 0.0

    inst_freq = 1.0 / intervals  # Hz
    omega_median = float(np.median(inst_freq)) * TWO_PI  # rad/s

    total_time = float(signal[-1] - signal[0])
    omega_median_hz = float(np.median(inst_freq))
    cumulative_phase = TWO_PI * omega_median_hz * total_time
    theta = float(cumulative_phase % TWO_PI)

    # Quality: inverse coefficient of variation of intervals (regularity)
    cv = (
        float(np.std(intervals) / np.mean(intervals))
        if np.mean(intervals) > 0
        else 1.0
    )
    quality = float(np.clip(1.0 / (1.0 + cv), 0.0, 1.0))

    return [
        PhaseState(
            theta=theta,
            omega=omega_median,
            amplitude=float(np.mean(inst_freq)),
            quality=quality,
            channel="I",
            node_id=self._node_id,
        )
    ]

quality_score

quality_score(phase_states: list[PhaseState]) -> float

Mean interval-regularity quality across phase states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states.

Returns

float Mean interval-regularity quality across phase states.

Source code in src/scpn_phase_orchestrator/oscillators/informational.py

def quality_score(self, phase_states: list[PhaseState]) -> float:
    """Mean interval-regularity quality across phase states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.

    Returns
    -------
    float
        Mean interval-regularity quality across phase states.
    """
    if not phase_states:
        return 0.0
    return float(np.mean([ps.quality for ps in phase_states]))

---

Symbolic Extraction (S)

SymbolicExtractor

SymbolicExtractor(n_states: int, node_id: str = "sym", mode: str = "ring")

Parameter	Type	Description
n_states	int	Vocabulary size (≥ 2)
node_id	str	Oscillator identifier
mode	str	"ring" or "graph"

Ring mode

Maps state index s to phase: θ_s = 2πs / N (mod 2π). Equispaced phases with gap = 2π/N.

Graph mode

Cumulative transition distances normalised to [0, 2π).

Quality scoring

Transition type	Quality
Single step (	Δs
Stalled (Δs = 0)	0.2
Large jump (	Δs
First state (no prior)	0.5

Omega derivation

ω is derived from consecutive phase differences divided by dt (1/sample_rate). For ring mode with single steps: ω = 2π/(N·dt).

Performance: extract(1000 states) < 1 ms.

symbolic

Symbolic-channel phase extraction from discrete state sequences.

SymbolicExtractor maps integer state indices onto ring or graph-walk phases for semiotic and finite-state systems. It rejects invalid state counts, non-integer signals, boolean arrays, complex arrays, and invalid sample rates so symbolic phases remain explicit and deterministic.

Classes

SymbolicExtractor

SymbolicExtractor(
    n_states: int,
    node_id: str = "sym",
    mode: str = "ring",
    *,
    initial_transition_quality: float = SYMBOLIC_INITIAL_TRANSITION_QUALITY_BASELINE,
)

Bases: PhaseExtractor

Phase extraction from discrete symbolic state sequences.

Maps discrete state indices to phases on the unit circle via theta = 2pis/N (ring-phase) or via graph-walk position.

Configure the symbolic oscillator over n_states discrete states.

Parameters

n_states : int total number of discrete states N. node_id : str identifier for generated PhaseState objects. mode : str "ring" for ring-phase, "graph" for graph-walk phase. initial_transition_quality : float quality assigned when no transition evidence exists yet (first sample / insufficient history).

Source code in src/scpn_phase_orchestrator/oscillators/symbolic.py

def __init__(
    self,
    n_states: int,
    node_id: str = "sym",
    mode: str = "ring",
    *,
    initial_transition_quality: float = (
        SYMBOLIC_INITIAL_TRANSITION_QUALITY_BASELINE
    ),
):
    """Configure the symbolic oscillator over ``n_states`` discrete states.

    Parameters
    ----------
    n_states : int
        total number of discrete states N.
    node_id : str
        identifier for generated PhaseState objects.
    mode : str
        "ring" for ring-phase, "graph" for graph-walk phase.
    initial_transition_quality : float
        quality assigned when no transition evidence exists yet (first sample /
        insufficient history).
    """
    n_states = _validate_n_states(n_states)
    if mode not in ("ring", "graph"):
        raise ValueError(f"mode must be 'ring' or 'graph', got {mode!r}")
    self._n_states = n_states
    self._node_id = _validate_node_id(node_id)
    self._mode = mode
    self._initial_transition_quality = _validate_initial_transition_quality(
        initial_transition_quality
    )

Methods:

extract

extract(
    signal: FloatArray | IntArray, sample_rate: float
) -> list[PhaseState]

Map discrete state indices to phases on the unit circle.

Parameters

signal : FloatArray | IntArray Input signal, shape (T,). sample_rate : float Sampling rate in Hz.

Returns

list[PhaseState] Discrete state indices to phases on the unit circle.

Source code in src/scpn_phase_orchestrator/oscillators/symbolic.py

def extract(
    self, signal: FloatArray | IntArray, sample_rate: float
) -> list[PhaseState]:
    """Map discrete state indices to phases on the unit circle.

    Parameters
    ----------
    signal : FloatArray | IntArray
        Input signal, shape ``(T,)``.
    sample_rate : float
        Sampling rate in Hz.

    Returns
    -------
    list[PhaseState]
        Discrete state indices to phases on the unit circle.
    """
    indices = _validate_signal(signal)
    sample_rate = _validate_sample_rate(sample_rate)
    if self._mode == "ring":
        if _HAS_RUST_SYMBOLIC:
            thetas = np.asarray(
                _rust_ring_phases(indices, self._n_states),
                dtype=np.float64,
            )
        else:
            thetas = TWO_PI * indices / self._n_states
    else:
        # Graph-walk: cumulative phase from state transitions
        # Each step adds phase proportional to the transition distance
        if len(indices) < 2:
            if _HAS_RUST_SYMBOLIC:
                thetas = np.asarray(
                    _rust_ring_phases(indices, self._n_states),
                    dtype=np.float64,
                )
            else:
                thetas = TWO_PI * indices.astype(np.float64) / self._n_states
        else:
            steps = np.abs(np.diff(indices)).astype(np.float64)
            cumulative = np.concatenate([[0.0], np.cumsum(steps)])
            total = cumulative[-1] if cumulative[-1] > 0 else 1.0
            if _HAS_RUST_SYMBOLIC:
                thetas = np.asarray(
                    _rust_graph_walk_phases(indices, self._n_states),
                    dtype=np.float64,
                )
            else:
                thetas = TWO_PI * cumulative / total

    thetas = thetas % TWO_PI
    dt = 1.0 / sample_rate
    omegas = np.zeros_like(thetas)
    if len(thetas) > 1:
        dtheta = np.diff(thetas)
        # Unwrap jumps larger than pi
        dtheta = (dtheta + np.pi) % TWO_PI - np.pi
        omegas[1:] = dtheta / dt

    states = []
    rust_qualities: FloatArray | None = None
    if _HAS_RUST_SYMBOLIC:
        rust_qualities = np.asarray(
            _rust_transition_qualities(
                indices,
                self._n_states,
                self._initial_transition_quality,
            ),
            dtype=np.float64,
        )
    for i in range(len(thetas)):
        states.append(
            PhaseState(
                theta=float(thetas[i]),
                omega=float(omegas[i]),
                amplitude=1.0,
                quality=(
                    float(rust_qualities[i])
                    if rust_qualities is not None
                    else self._transition_quality(indices, i)
                ),
                channel="S",
                node_id=self._node_id,
            )
        )
    return states

quality_score

quality_score(phase_states: list[PhaseState]) -> float

Mean transition quality across phase states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states.

Returns

float Mean transition quality across phase states.

Source code in src/scpn_phase_orchestrator/oscillators/symbolic.py

def quality_score(self, phase_states: list[PhaseState]) -> float:
    """Mean transition quality across phase states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.

    Returns
    -------
    float
        Mean transition quality across phase states.
    """
    if not phase_states:
        return 0.0
    return float(np.mean([ps.quality for ps in phase_states]))

---

Wavelet-ridge extractor (physical channel)

WaveletExtractor is a band-adaptive alternative to the Hilbert extractor: it computes a complex Morlet continuous wavelet transform across a log-spaced frequency bank, selects the dominant energy ridge over a cone-of-influence-safe interior region, and reads the analytic phase along that ridge. The terminal phase is extrapolated from a COI-safe interior sample at the ridge frequency, so the corrupted signal edge is avoided. The Morlet wavelet is ψ(t) = π^(−1/4)·exp(i·ω₀·t/s)·exp(−(t/s)²/2)/√s with ω₀ = 6; a pure-NumPy path is used because SciPy 1.15 removed cwt/morlet2. Prefer it over Hilbert when a dominant oscillation sits in broadband noise or slow drift.

wavelet

Physical-channel phase extraction via a complex Morlet wavelet ridge.

WaveletExtractor computes a continuous wavelet transform with a bank of complex Morlet wavelets, selects the dominant scale (the energy ridge over a cone-of-influence-safe interior region), and reads the analytic phase along that ridge. Unlike the broadband Hilbert transform, the wavelet ridge is band-adaptive and therefore robust to broadband noise and slow drift around a dominant oscillation. The terminal phase is extrapolated from a COI-safe interior sample using the ridge frequency, avoiding the corrupted signal edge.

The complex Morlet is psi(t) = pi^(-1/4) * exp(i*w0*t/s) * exp(-(t/s)^2/2) / sqrt(s) with w0 = 6 (the standard admissibility-approximating choice). The peak frequency of scale s is f = w0 * fs / (2*pi*s). A pure-Python/NumPy path is used because SciPy 1.15 removed cwt/morlet2 and PyWavelets is not a required dependency; the extractor preserves the PhaseState contract.

Classes

WaveletExtractor

WaveletExtractor(node_id: str = 'wav_0')

Bases: PhaseExtractor

Extracts instantaneous phase from a waveform via a complex Morlet ridge.

Source code in src/scpn_phase_orchestrator/oscillators/wavelet.py

def __init__(self, node_id: str = "wav_0"):
    self._node_id = _validate_node_id(node_id)

Methods:

extract

extract(
    signal: FloatArray, sample_rate: float
) -> list[PhaseState]

Extract phase from a 1-D waveform via the dominant Morlet ridge.

Parameters

signal : FloatArray Input signal, shape (T,). sample_rate : float Sampling rate in Hz.

Returns

list[PhaseState] A single PhaseState on channel "P" carrying the terminal phase, the ridge angular frequency, amplitude, and ridge-regularity quality.

Source code in src/scpn_phase_orchestrator/oscillators/wavelet.py

def extract(self, signal: FloatArray, sample_rate: float) -> list[PhaseState]:
    """Extract phase from a 1-D waveform via the dominant Morlet ridge.

    Parameters
    ----------
    signal : FloatArray
        Input signal, shape ``(T,)``.
    sample_rate : float
        Sampling rate in Hz.

    Returns
    -------
    list[PhaseState]
        A single `PhaseState` on channel ``"P"`` carrying the terminal phase,
        the ridge angular frequency, amplitude, and ridge-regularity quality.
    """
    signal = _validate_signal(signal)
    sample_rate = _validate_sample_rate(sample_rate)

    centred = signal - float(np.mean(signal))
    rms = float(np.sqrt(np.mean(np.square(centred))))
    amplitude = float(np.sqrt(2.0) * rms)
    n = centred.size

    freqs = self._frequency_grid(n, sample_rate)
    if freqs.size == 0 or rms < 1e-15:
        theta = 0.0 if centred[-1] >= 0.0 else float(np.pi)
        return [self._state(theta % TWO_PI, 0.0, amplitude, 0.0)]

    scales = _MORLET_W0 * sample_rate / (TWO_PI * freqs)
    ridge_idx, coeffs, ridge_scale = self._ridge(centred, scales)
    f_ridge = float(freqs[ridge_idx])
    omega = float(TWO_PI * f_ridge)

    theta = self._terminal_phase(coeffs, ridge_scale, omega, sample_rate, n)
    quality = self._ridge_quality(coeffs, ridge_scale)

    return [self._state(theta, omega, amplitude, quality)]

quality_score

quality_score(phase_states: list[PhaseState]) -> float

Mean extraction quality across phase states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states.

Returns

float Mean extraction quality across phase states.

Source code in src/scpn_phase_orchestrator/oscillators/wavelet.py

def quality_score(self, phase_states: list[PhaseState]) -> float:
    """Mean extraction quality across phase states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.

    Returns
    -------
    float
        Mean extraction quality across phase states.
    """
    if not phase_states:
        return 0.0
    return float(np.mean([ps.quality for ps in phase_states]))

---

Zero-crossing extractor (physical channel)

ZeroCrossingExtractor recovers phase from interpolated zero crossings: each crossing is a half-cycle (π of phase advance), absolute phase is anchored to the crossing direction (a rising crossing ≡ 0, a falling crossing ≡ π, the sine convention), and a Schmitt-trigger deadband (a fraction of the RMS) suppresses spurious noise-induced crossings. Angular frequency comes from the mean half-period and quality from the regularity of the half-period intervals. Prefer it for sharply non-sinusoidal periodic signals where a single analytic phase is ill-defined.

zero_crossing

Physical-channel phase extraction from zero crossings.

ZeroCrossingExtractor recovers instantaneous phase from a real waveform by locating its zero crossings (with sub-sample linear interpolation), treating each crossing as a half-cycle (pi of phase advance), and anchoring absolute phase to the crossing direction (rising crossing ≡ phase 0, falling crossing ≡ phase pi, matching the sine convention). It is robust to a constant offset (the mean is removed) and reports an angular frequency from the mean half-period and a quality from the regularity of the half-period intervals. The extractor produces the same PhaseState contract as the Hilbert-based PhysicalExtractor but is preferable for sharply non-sinusoidal periodic signals where a single dominant analytic phase is ill-defined.

Classes

ZeroCrossingExtractor

ZeroCrossingExtractor(
    node_id: str = "zc_0", *, hysteresis: float = 0.5
)

Bases: PhaseExtractor

Extracts instantaneous phase from a waveform via interpolated zero crossings.

Source code in src/scpn_phase_orchestrator/oscillators/zero_crossing.py

def __init__(self, node_id: str = "zc_0", *, hysteresis: float = 0.5):
    self._node_id = _validate_node_id(node_id)
    if not isfinite(hysteresis) or not (0.0 <= hysteresis < 1.0):
        raise ValueError("hysteresis must be a fraction in [0, 1)")
    self._hysteresis = float(hysteresis)

Methods:

extract

extract(
    signal: FloatArray, sample_rate: float
) -> list[PhaseState]

Extract phase from a 1-D waveform via interpolated zero crossings.

Parameters

signal : FloatArray Input signal, shape (T,). sample_rate : float Sampling rate in Hz.

Returns

list[PhaseState] A single PhaseState on channel "P" carrying the terminal phase, angular frequency, amplitude, and crossing-regularity quality.

Source code in src/scpn_phase_orchestrator/oscillators/zero_crossing.py

def extract(self, signal: FloatArray, sample_rate: float) -> list[PhaseState]:
    """Extract phase from a 1-D waveform via interpolated zero crossings.

    Parameters
    ----------
    signal : FloatArray
        Input signal, shape ``(T,)``.
    sample_rate : float
        Sampling rate in Hz.

    Returns
    -------
    list[PhaseState]
        A single `PhaseState` on channel ``"P"`` carrying the terminal phase,
        angular frequency, amplitude, and crossing-regularity quality.
    """
    signal = _validate_signal(signal)
    sample_rate = _validate_sample_rate(sample_rate)

    centred = signal - float(np.mean(signal))
    rms = float(np.sqrt(np.mean(np.square(centred))))
    amplitude = float(np.sqrt(2.0) * rms)

    crossings, rising = self._confirmed_crossings(centred, self._hysteresis * rms)

    if crossings.size < 2 or rms < 1e-15:
        # Too few crossings (near-DC or sub-half-cycle window): no reliable
        # phase. Report a terminal phase consistent with the last sample sign.
        theta = 0.0 if centred[-1] >= 0.0 else float(np.pi)
        return [
            PhaseState(
                theta=theta % TWO_PI,
                omega=0.0,
                amplitude=amplitude,
                quality=0.0,
                channel="P",
                node_id=self._node_id,
            )
        ]

    half_periods = np.diff(crossings)  # samples per half-cycle
    mean_half = float(np.mean(half_periods))
    half_seconds = mean_half / sample_rate
    omega = float(np.pi / half_seconds)  # pi advance per half-period

    last_cross = float(crossings[-1])
    anchor = 0.0 if bool(rising[-1]) else float(np.pi)
    elapsed_seconds = (len(centred) - 1 - last_cross) / sample_rate
    theta = (anchor + omega * elapsed_seconds) % TWO_PI

    quality = self._interval_quality(half_periods)

    return [
        PhaseState(
            theta=float(theta),
            omega=omega,
            amplitude=amplitude,
            quality=quality,
            channel="P",
            node_id=self._node_id,
        )
    ]

quality_score

quality_score(phase_states: list[PhaseState]) -> float

Mean extraction quality across phase states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states.

Returns

float Mean extraction quality across phase states.

Source code in src/scpn_phase_orchestrator/oscillators/zero_crossing.py

def quality_score(self, phase_states: list[PhaseState]) -> float:
    """Mean extraction quality across phase states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.

    Returns
    -------
    float
        Mean extraction quality across phase states.
    """
    if not phase_states:
        return 0.0
    return float(np.mean([ps.quality for ps in phase_states]))

---

Quality Scoring

PhaseQualityScorer

Method	Signature	Description
score	(states) → float	Amplitude-weighted mean quality
detect_collapse	(states, threshold=0.1) → bool	True if >50% below threshold
downweight_mask	(states, min_quality=0.3) → NDArray[np.float64]	Weight array, zeros below min

Downweight mask in pipeline

The mask is applied to the coupling matrix before engine evaluation:

mask = scorer.downweight_mask(states, min_quality=0.3)
knm_gated = knm * mask[:, None] * mask[None, :]
# Low-quality oscillators decoupled from high-quality ones

This prevents noisy phase estimates from corrupting the synchronisation dynamics. Only oscillators with quality ≥ min_quality participate.

Performance: downweight_mask(100 states) < 50 μs.

quality

Quality aggregation and collapse detection for extracted phase states.

The scorer turns per-oscillator extraction quality into weighted aggregate signals for runtime gating and diagnostics. Empty state sets collapse to safe defaults, low-quality states can be masked, and amplitude weighting prevents near-zero signals from dominating quality summaries.

Classes

PhaseQualityScorer

PhaseQualityScorer(
    collapse_threshold: float = 0.1,
    min_quality: float = 0.3,
)

Aggregate quality scoring and collapse detection for phase state arrays.

Source code in src/scpn_phase_orchestrator/oscillators/quality.py

def __init__(self, collapse_threshold: float = 0.1, min_quality: float = 0.3):
    if not np.isfinite(collapse_threshold):
        raise ValueError("collapse_threshold must be finite")
    if not np.isfinite(min_quality):
        raise ValueError("min_quality must be finite")
    if not 0.0 <= collapse_threshold <= 1.0:
        raise ValueError("collapse_threshold must be in [0, 1]")
    if not 0.0 <= min_quality <= 1.0:
        raise ValueError("min_quality must be in [0, 1]")
    self._collapse_threshold = float(collapse_threshold)
    self._min_quality = float(min_quality)
    self._rust = (
        _RustPhaseQualityScorer(
            collapse_threshold=self._collapse_threshold,
            min_quality=self._min_quality,
        )
        if _RustPhaseQualityScorer is not None
        else None
    )

Methods:

score

score(phase_states: list[PhaseState]) -> float

Weighted average quality across all phase states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states.

Returns

float Weighted average quality across all phase states.

Source code in src/scpn_phase_orchestrator/oscillators/quality.py

def score(self, phase_states: list[PhaseState]) -> float:
    """Weighted average quality across all phase states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.

    Returns
    -------
    float
        Weighted average quality across all phase states.
    """
    if not phase_states:
        return 0.0
    qualities = np.array([ps.quality for ps in phase_states])
    amplitudes = np.array([ps.amplitude for ps in phase_states])
    if self._rust is not None:
        return float(self._rust.score(qualities.tolist(), amplitudes.tolist()))
    weights = np.maximum(amplitudes, 1e-12)
    return float(np.average(qualities, weights=weights))

detect_collapse

detect_collapse(
    phase_states: list[PhaseState], threshold: float = 0.1
) -> bool

Return True if quality is below threshold for the majority of states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states. threshold : float Decision threshold.

Returns

bool True if quality is below threshold for the majority of states.

Raises

ValueError If the inputs are invalid or inconsistent.

Source code in src/scpn_phase_orchestrator/oscillators/quality.py

def detect_collapse(
    self, phase_states: list[PhaseState], threshold: float = 0.1
) -> bool:
    """Return True if quality is below threshold for the majority of states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.
    threshold : float
        Decision threshold.

    Returns
    -------
    bool
        True if quality is below threshold for the majority of states.

    Raises
    ------
    ValueError
        If the inputs are invalid or inconsistent.
    """
    if not phase_states:
        return True
    if not np.isfinite(threshold):
        raise ValueError("threshold must be finite")
    threshold = float(threshold)
    if not 0.0 <= threshold <= 1.0:
        raise ValueError("threshold must be in [0, 1]")
    if self._rust is not None and threshold == self._collapse_threshold:
        qualities = [ps.quality for ps in phase_states]
        return bool(self._rust.is_collapsed(qualities))
    below = sum(1 for ps in phase_states if ps.quality < threshold)
    return below > len(phase_states) / 2

downweight_mask

downweight_mask(
    phase_states: list[PhaseState], min_quality: float = 0.3
) -> FloatArray

Weight array in [0,1], zeros below min_quality.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states. min_quality : float Minimum extraction quality.

Returns

FloatArray Weight array in [0,1], zeros below min_quality.

Raises

ValueError If the inputs are invalid or inconsistent.

Source code in src/scpn_phase_orchestrator/oscillators/quality.py

def downweight_mask(
    self, phase_states: list[PhaseState], min_quality: float = 0.3
) -> FloatArray:
    """Weight array in [0,1], zeros below min_quality.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.
    min_quality : float
        Minimum extraction quality.

    Returns
    -------
    FloatArray
        Weight array in [0,1], zeros below min_quality.

    Raises
    ------
    ValueError
        If the inputs are invalid or inconsistent.
    """
    if not phase_states:
        return np.array([], dtype=np.float64)
    if not np.isfinite(min_quality):
        raise ValueError("min_quality must be finite")
    min_quality = float(min_quality)
    if not 0.0 <= min_quality <= 1.0:
        raise ValueError("min_quality must be in [0, 1]")
    qualities = np.array([ps.quality for ps in phase_states])
    if self._rust is not None and min_quality == self._min_quality:
        return np.asarray(
            self._rust.downweight_mask(qualities.tolist()),
            dtype=np.float64,
        )
    mask = np.where(qualities >= min_quality, qualities, 0.0)
    return mask.astype(np.float64)

---

Base Types

base

Shared phase-state contract and extractor interface.

PhaseState is the typed handoff record from channel-specific extractors into binding, quality scoring, and UPDE initialisation. PhaseExtractor defines the minimal extraction/quality interface implemented by physical waveform, informational event, and symbolic state-sequence extractors.

Classes

PhaseState dataclass

PhaseState(
    theta: float,
    omega: float,
    amplitude: float,
    quality: float,
    channel: str,
    node_id: str,
)

Extracted phase, frequency, amplitude, and quality for one oscillator.

PhaseExtractor

Bases: ABC

Abstract base for signal-to-phase extraction algorithms.

Methods:

extract abstractmethod

extract(
    signal: FloatArray, sample_rate: float
) -> list[PhaseState]

Extract phase states from a raw signal at the given sample rate.

Parameters

signal : FloatArray Input signal, shape (T,). sample_rate : float Sampling rate in Hz.

Returns

list[PhaseState] Phase states from a raw signal at the given sample rate.

Source code in src/scpn_phase_orchestrator/oscillators/base.py

@abstractmethod
def extract(self, signal: FloatArray, sample_rate: float) -> list[PhaseState]:
    """Extract phase states from a raw signal at the given sample rate.

    Parameters
    ----------
    signal : FloatArray
        Input signal, shape ``(T,)``.
    sample_rate : float
        Sampling rate in Hz.

    Returns
    -------
    list[PhaseState]
        Phase states from a raw signal at the given sample rate.
    """
    ...

quality_score abstractmethod

quality_score(phase_states: list[PhaseState]) -> float

Aggregate quality metric (0..1) over a set of extracted phase states.

Parameters

phase_states : list[PhaseState] Extracted per-oscillator phase states.

Returns

float Aggregate quality metric (0..1) over a set of extracted phase states.

Source code in src/scpn_phase_orchestrator/oscillators/base.py

@abstractmethod
def quality_score(self, phase_states: list[PhaseState]) -> float:
    """Aggregate quality metric (0..1) over a set of extracted phase states.

    Parameters
    ----------
    phase_states : list[PhaseState]
        Extracted per-oscillator phase states.

    Returns
    -------
    float
        Aggregate quality metric (0..1) over a set of extracted phase states.
    """
    ...

Phase Initialisation

Utilities for deterministic and random initial phase generation used in simulation setup and reproducible experiment seeds.

init_phases

Synthetic binding-aware initial phase generation.

extract_initial_phases uses the binding's oscillator families to generate small deterministic synthetic signals per P/I/S channel, extract their phases, and produce a finite initial phase vector for UPDE startup. It validates omega length, seed, symbolic state counts, and extractor families before falling back to seeded random phases for unsupported channel families.

Classes

Functions:

extract_initial_phases

extract_initial_phases(
    spec: BindingSpec, omegas: FloatArray, seed: int = 42
) -> FloatArray

Extract initial phases from channels defined in binding_spec.

For each oscillator, generates a synthetic signal matching the family channel or extractor semantics and extracts the phase. Falls back to random phase if extraction fails.

Returns (n_osc,) array of initial phases in [0, 2*pi).

Parameters

spec : BindingSpec The binding specification. omegas : FloatArray Natural frequencies in rad/s, shape (N,). seed : int Seed for the deterministic RNG.

Returns

FloatArray Initial phases from channels defined in binding_spec.

Source code in src/scpn_phase_orchestrator/oscillators/init_phases.py

def extract_initial_phases(
    spec: BindingSpec,
    omegas: FloatArray,
    seed: int = 42,
) -> FloatArray:
    """Extract initial phases from channels defined in binding_spec.

    For each oscillator, generates a synthetic signal matching the family
    channel or extractor semantics and extracts the phase. Falls back to
    random phase if extraction fails.

    Returns (n_osc,) array of initial phases in [0, 2*pi).

    Parameters
    ----------
    spec : BindingSpec
        The binding specification.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    seed : int
        Seed for the deterministic RNG.

    Returns
    -------
    FloatArray
        Initial phases from channels defined in binding_spec.
    """
    omegas = _validate_omegas(omegas, expected_count=_oscillator_count(spec))
    seed = _validate_seed(seed)
    rng = np.random.default_rng(seed)
    n_osc = len(omegas)
    phases = np.zeros(n_osc)
    t = np.linspace(0, 1.0, 256)

    families = spec.oscillator_families
    physical_extractors: dict[str, PhaseExtractor] = {
        name: build_extractor(name, node_id="init_physical")
        for name in ("hilbert", "wavelet", "zero_crossing")
    }
    informational_extractor = InformationalExtractor(node_id="init_informational")
    symbolic_n_states = _get_n_states(families)
    symbolic_extractor = SymbolicExtractor(
        n_states=symbolic_n_states,
        node_id="init_symbolic",
    )
    symbolic_pending: list[tuple[int, int]] = []

    osc_idx = 0
    for layer in spec.layers:
        for _osc_id in layer.oscillator_ids:
            omega = omegas[osc_idx]
            family = _resolve_family(layer.family, families, osc_idx)
            channel = family.channel if family is not None else "P"
            extractor_type = family.extractor_type if family is not None else "hilbert"

            if channel == "P" or extractor_type in _PHYSICAL_EXTRACTORS:
                signal = np.sin(omega * TWO_PI * t) + rng.normal(0, 0.1, len(t))
                algorithm = resolve_extractor_type(extractor_type)
                extractor = physical_extractors.get(
                    algorithm, physical_extractors["hilbert"]
                )
                states = extractor.extract(signal, sample_rate=256.0)
                phases[osc_idx] = states[-1].theta if states else rng.uniform(0, TWO_PI)

            elif channel == "I" or extractor_type in _INFORMATIONAL_EXTRACTORS:
                n_events = max(3, int(omega * 10))
                timestamps = np.sort(rng.uniform(0, 1.0, n_events))
                states = informational_extractor.extract(timestamps, sample_rate=1.0)
                phases[osc_idx] = states[-1].theta if states else rng.uniform(0, TWO_PI)

            elif channel == "S" or extractor_type in _SYMBOLIC_EXTRACTORS:
                state_idx = int(rng.integers(0, symbolic_n_states))
                symbolic_pending.append((osc_idx, state_idx))

            else:
                phases[osc_idx] = rng.uniform(0, TWO_PI)

            osc_idx += 1

    if symbolic_pending:
        symbolic_indices = np.array(
            [state_idx for _, state_idx in symbolic_pending],
            dtype=np.int64,
        )
        symbolic_states = symbolic_extractor.extract(
            cast(FloatArray, symbolic_indices),
            sample_rate=1.0,
        )
        for (phase_idx, _), state in zip(
            symbolic_pending, symbolic_states, strict=True
        ):
            phases[phase_idx] = state.theta

    return phases

---

Phase reduction

Model-free phase reduction: a dependency-light evaluator of a trained phase autoencoder (see nn.phase_autoencoder) that recovers the asymptotic phase Θ(x) and the phase-sensitivity function Z(θ) — the phase response curve — from frozen NumPy weights, with no JAX on the control path.

phase_reduction

A pure-NumPy evaluator for a trained phase autoencoder.

The phase autoencoder (nn.phase_autoencoder) is trained with JAX, but the asymptotic phase and the phase-sensitivity function it learns are needed on the control path, which must stay dependency-light. This module evaluates the trained encoder/decoder — frozen to plain NumPy weights — without importing JAX.

Given the trained encoder g(x) = (Ỹ₁, Ỹ₂, Ỹ₃) (a ReLU multilayer perceptron), the asymptotic phase is Θ(x) = atan2(Ỹ₂, Ỹ₁) (the unit-circle normalisation cancels inside atan2). The phase-sensitivity function — the gradient of the phase with respect to the state, evaluated on the limit cycle — is

Z(θ) = ∇ₓ Θ |_{x = decode(cos θ, sin θ, 0)}
     = (∂Θ/∂Ỹ₁) ∇ₓ Ỹ₁ + (∂Θ/∂Ỹ₂) ∇ₓ Ỹ₂,

with ∂Θ/∂Ỹ₁ = −Ỹ₂/(Ỹ₁²+Ỹ₂²), ∂Θ/∂Ỹ₂ = Ỹ₁/(Ỹ₁²+Ỹ₂²) and the encoder Jacobian computed by exact reverse-mode through the ReLU network. This is the phase response curve (Nakao 2016) recovered model-free from data.

References

• Yawata, Fukami, Taira & Nakao 2024, Chaos 34, 063111 — phase autoencoder.
• Nakao 2016, Contemp. Phys. 57, 188 — phase reduction theory.

Classes

PhaseReductionWeights dataclass

PhaseReductionWeights(
    encoder_weights: tuple[FloatArray, ...],
    encoder_biases: tuple[FloatArray, ...],
    decoder_weights: tuple[FloatArray, ...],
    decoder_biases: tuple[FloatArray, ...],
    omega: float,
    decay: float,
    state_dim: int,
)

Frozen encoder/decoder weights and (ω, λ) of a phase autoencoder.

Parameters

encoder_weights, encoder_biases : tuple[numpy.ndarray, ...] Per-layer encoder weight matrices (out, in) and bias vectors. decoder_weights, decoder_biases : tuple[numpy.ndarray, ...] Per-layer decoder weight matrices and bias vectors. omega : float The learned angular frequency ω. decay : float The learned amplitude decay λ < 0. state_dim : int The oscillator state dimension n.

PhaseReducer dataclass

PhaseReducer(weights: PhaseReductionWeights)

A dependency-light evaluator of a trained phase autoencoder.

Parameters

weights : PhaseReductionWeights The frozen encoder/decoder weights and (ω, λ).

Attributes

omega property

omega: float

The learned angular frequency ω.

Returns

float The learned angular frequency ω.

decay property

decay: float

The learned amplitude decay λ < 0.

Returns

float The learned amplitude decay λ, strictly negative.

Methods:

asymptotic_phase

asymptotic_phase(state: FloatArray) -> float

Return the asymptotic phase Θ(x) = atan2(Ỹ₂, Ỹ₁) of a state.

Parameters

state : numpy.ndarray The oscillator state x of shape (n,).

Returns

float The asymptotic phase in (−π, π].

Raises

ValueError If state is not a finite vector of length state_dim.

Source code in src/scpn_phase_orchestrator/oscillators/phase_reduction.py

def asymptotic_phase(self, state: FloatArray) -> float:
    """Return the asymptotic phase ``Θ(x) = atan2(Ỹ₂, Ỹ₁)`` of a state.

    Parameters
    ----------
    state : numpy.ndarray
        The oscillator state ``x`` of shape ``(n,)``.

    Returns
    -------
    float
        The asymptotic phase in ``(−π, π]``.

    Raises
    ------
    ValueError
        If ``state`` is not a finite vector of length ``state_dim``.
    """
    vector = _validate_state(state, name="state", dim=self.weights.state_dim)
    raw = self._encode_raw(vector)
    return float(np.arctan2(raw[1], raw[0]))

encode_observables

encode_observables(states: FloatArray) -> FloatArray

Lift a batch of states to the unnormalised encoder latent (K, 3).

These are the model-free Koopman observables: the learned coordinate in which the phase autoencoder's dynamics are (approximately) linear, so a :class:~scpn_phase_orchestrator.monitor.koopman_edmd.KoopmanPredictor fitted in them captures nonlinear oscillator dynamics that the analytic dictionaries miss.

Parameters

states : numpy.ndarray A batch of states of shape (K, state_dim).

Returns

numpy.ndarray The unnormalised latent batch of shape (K, 3).

Raises

ValueError If states is not a finite (K, state_dim) array.

Source code in src/scpn_phase_orchestrator/oscillators/phase_reduction.py

def encode_observables(self, states: FloatArray) -> FloatArray:
    """Lift a batch of states to the unnormalised encoder latent ``(K, 3)``.

    These are the model-free Koopman observables: the learned coordinate in
    which the phase autoencoder's dynamics are (approximately) linear, so a
    :class:`~scpn_phase_orchestrator.monitor.koopman_edmd.KoopmanPredictor`
    fitted in them captures nonlinear oscillator dynamics that the analytic
    dictionaries miss.

    Parameters
    ----------
    states : numpy.ndarray
        A batch of states of shape ``(K, state_dim)``.

    Returns
    -------
    numpy.ndarray
        The unnormalised latent batch of shape ``(K, 3)``.

    Raises
    ------
    ValueError
        If ``states`` is not a finite ``(K, state_dim)`` array.
    """
    matrix = _validate_states(states, name="states", dim=self.weights.state_dim)
    activation = matrix
    weights = self.weights.encoder_weights
    biases = self.weights.encoder_biases
    last = len(weights) - 1
    for index, (weight, bias) in enumerate(zip(weights, biases, strict=True)):
        activation = activation @ weight.T + bias
        if index != last:
            activation = np.maximum(activation, 0.0)
    return np.ascontiguousarray(activation, dtype=np.float64)

reconstruct

reconstruct(phase: float) -> FloatArray

Reconstruct the on-cycle state at a phase via the decoder.

Parameters

phase : float The phase θ to reconstruct on the limit cycle.

Returns

numpy.ndarray The decoded state decode(cos θ, sin θ, 0) of shape (n,).

Source code in src/scpn_phase_orchestrator/oscillators/phase_reduction.py

def reconstruct(self, phase: float) -> FloatArray:
    """Reconstruct the on-cycle state at a phase via the decoder.

    Parameters
    ----------
    phase : float
        The phase ``θ`` to reconstruct on the limit cycle.

    Returns
    -------
    numpy.ndarray
        The decoded state ``decode(cos θ, sin θ, 0)`` of shape ``(n,)``.
    """
    latent = np.array(
        [np.cos(float(phase)), np.sin(float(phase)), 0.0], dtype=np.float64
    )
    return np.ascontiguousarray(self._decode(latent), dtype=np.float64)

phase_sensitivity

phase_sensitivity(phase: float) -> FloatArray

Return the phase-sensitivity function Z(θ) = ∇ₓ Θ on the cycle.

Parameters

phase : float The phase θ on the limit cycle at which to evaluate Z.

Returns

numpy.ndarray The phase response curve value Z(θ) of shape (n,).

Source code in src/scpn_phase_orchestrator/oscillators/phase_reduction.py

def phase_sensitivity(self, phase: float) -> FloatArray:
    """Return the phase-sensitivity function ``Z(θ) = ∇ₓ Θ`` on the cycle.

    Parameters
    ----------
    phase : float
        The phase ``θ`` on the limit cycle at which to evaluate ``Z``.

    Returns
    -------
    numpy.ndarray
        The phase response curve value ``Z(θ)`` of shape ``(n,)``.
    """
    state = self.reconstruct(phase)
    raw = self._encode_raw(state)
    denominator = raw[0] ** 2 + raw[1] ** 2 + 1.0e-12
    dphase_dy1 = -raw[1] / denominator
    dphase_dy2 = raw[0] / denominator
    jacobian = self._encoder_jacobian(state)
    sensitivity = dphase_dy1 * jacobian[0] + dphase_dy2 * jacobian[1]
    return np.ascontiguousarray(sensitivity, dtype=np.float64)

---

Extractor factory

build_extractor maps a binding extractor_type — a channel alias (physical/informational/symbolic) or a canonical algorithm name (hilbert/wavelet/zero_crossing/event/ring/graph) — to the concrete PhaseExtractor that implements it. Aliases resolve through resolve_extractor_type; an unknown type raises ValueError (fail-closed) rather than silently degrading to a default algorithm.

factory

Construct the phase extractor named by a binding extractor_type.

build_extractor maps a domainpack extractor_type (a channel alias such as physical/informational/symbolic or a canonical algorithm name such as hilbert/wavelet/zero_crossing/event/ring/graph) to the concrete PhaseExtractor that implements it. Aliases are resolved through resolve_extractor_type; an unknown type raises ValueError (fail-closed) rather than silently degrading to a default algorithm.

Classes

Functions:

build_extractor

build_extractor(
    extractor_type: str,
    *,
    node_id: str = "extractor",
    n_states: int = 2,
    config: Mapping[str, object] | None = None,
) -> PhaseExtractor

Build the PhaseExtractor for a binding extractor_type.

Parameters

extractor_type : str A channel alias (physical/informational/symbolic) or a canonical algorithm name (hilbert/wavelet/zero_crossing/ event/ring/graph). node_id : str Identifier stamped onto the extractor's emitted PhaseState records. n_states : int Number of discrete states for symbolic (ring/graph) extractors; ignored by the continuous and event extractors. config : Mapping[str, object] | None Optional oscillator-family config from the binding spec. For the physical/hilbert extractor its band/filter_order/edge_trim keys select the opt-in zero-phase band-pass and edge-trim; other extractors ignore it.

Returns

PhaseExtractor The extractor implementing extractor_type.

Raises

ValueError If extractor_type does not resolve to a known algorithm.

Source code in src/scpn_phase_orchestrator/oscillators/factory.py

def build_extractor(
    extractor_type: str,
    *,
    node_id: str = "extractor",
    n_states: int = 2,
    config: Mapping[str, object] | None = None,
) -> PhaseExtractor:
    """Build the `PhaseExtractor` for a binding ``extractor_type``.

    Parameters
    ----------
    extractor_type : str
        A channel alias (``physical``/``informational``/``symbolic``) or a
        canonical algorithm name (``hilbert``/``wavelet``/``zero_crossing``/
        ``event``/``ring``/``graph``).
    node_id : str
        Identifier stamped onto the extractor's emitted `PhaseState` records.
    n_states : int
        Number of discrete states for symbolic (``ring``/``graph``) extractors;
        ignored by the continuous and event extractors.
    config : Mapping[str, object] | None
        Optional oscillator-family ``config`` from the binding spec. For the
        physical/``hilbert`` extractor its ``band``/``filter_order``/``edge_trim``
        keys select the opt-in zero-phase band-pass and edge-trim; other extractors
        ignore it.

    Returns
    -------
    PhaseExtractor
        The extractor implementing ``extractor_type``.

    Raises
    ------
    ValueError
        If ``extractor_type`` does not resolve to a known algorithm.
    """
    algorithm = resolve_extractor_type(extractor_type)
    if algorithm == "hilbert":
        return PhysicalExtractor(node_id=node_id, **_physical_kwargs(config))
    if algorithm == "wavelet":
        return WaveletExtractor(node_id=node_id)
    if algorithm == "zero_crossing":
        return ZeroCrossingExtractor(node_id=node_id)
    if algorithm == "event":
        return InformationalExtractor(node_id=node_id)
    if algorithm in ("ring", "graph"):
        return SymbolicExtractor(n_states=n_states, node_id=node_id, mode=algorithm)
    raise ValueError(
        f"unknown extractor_type {extractor_type!r} (resolved {algorithm!r}); "
        "expected one of hilbert, wavelet, zero_crossing, event, ring, graph "
        "or an alias physical/informational/symbolic"
    )

---

Cross-channel composition

A domain can use multiple channels simultaneously. The binding spec declares which channels are active and how they map to oscillator indices:

layers:
  - name: voltage
    channel: P
    indices: [0, 1, 2, 3]
  - name: event_rate
    channel: I
    indices: [4, 5]
  - name: protocol_state
    channel: S
    indices: [6, 7]

All channels produce PhaseState with the same fields, so the engine treats them uniformly. The channel field enables channel-aware analysis (e.g., computing R separately for P and I oscillators).

Rust FFI acceleration

PhysicalExtractor uses spo_kernel.physical_extract() when the Rust extension is installed. The Rust path computes the Hilbert transform and phase extraction in a single pass, avoiding Python/NumPy overhead for large signals.

Parity is verified in tests/test_oscillator_physical.py::test_rust_python_parity with tolerance atol=1e-10 for phase, rtol=0.01 for frequency.

---

Performance summary

Operation	Budget	Rust	Notes
PhysicalExtractor.extract(1s @ 1kHz)	< 5 ms	< 1 ms	Hilbert transform
InformationalExtractor.extract(100 ts)	< 500 μs	—	numpy operations
SymbolicExtractor.extract(1000 states)	< 1 ms	—	ring mapping
PhaseQualityScorer.downweight_mask(100)	< 50 μs	—	array comparison

Domain examples

Neuroscience (EEG)

# 64-channel EEG → 64 P-channel oscillators
extractor = PhysicalExtractor(node_id="eeg")
for ch in range(64):
    states = extractor.extract(eeg_data[ch], fs=256.0)
    phases[ch] = states[0].theta
    omegas[ch] = states[0].omega

Microservices (queue depths)

# 12 services → 12 I-channel oscillators
extractor = InformationalExtractor(node_id="svc")
for svc in services:
    timestamps = svc.request_timestamps()
    states = extractor.extract(timestamps, sample_rate=0.0)
    phases[svc.id] = states[0].theta

Genomic sequences

# DNA codons → S-channel oscillators
extractor = SymbolicExtractor(n_states=64, mode="ring")
codon_indices = encode_codons(sequence)
states = extractor.extract(codon_indices, sample_rate=1.0)