Delay Embedding — Phase-Space Reconstruction

Why this module is exposed

Delay embedding is how SPO can recover state-space structure when only scalar observations are available. It lets monitoring and anomaly workflows operate from raw traces without requiring a separate external embedding stack.

The module is intentionally strict at boundaries so a downstream controller uses phase-space features with deterministic semantics.

The monitor.embedding module reconstructs phase-space trajectories from scalar oscillator traces. It provides delay-coordinate embedding, Fraser-Swinney average mutual information, nearest-neighbor distances, and Python-side wrappers for optimal delay and embedding dimension.

Operational use

Use this module when a monitor needs geometry and short-term predictability from non-vector observations (for example, single-channel physiology traces mapped into state-space structure).

A practical sequence is:

1. Build a validated embedded matrix (delay_embed or auto_embed);
2. Select delay/dimension from information-theoretic and neighborhood diagnostics;
3. Feed the reconstructed state to downstream monitor logic;
4. Validate resulting geometrical signals against replay baselines before using them in policy experiments.

Keep this boundary explicit: embedding is a feature-construction layer, not a replacement for raw trace quality checks.

API

from scpn_phase_orchestrator.monitor.embedding import (
    delay_embed,
    mutual_information,
    nearest_neighbor_distances,
    optimal_delay,
    optimal_dimension,
    auto_embed,
)

delay_embed(signal, delay, dimension) returns the standard delay-coordinate matrix:

v(t) = [x(t), x(t + tau), x(t + 2 tau), ...]

mutual_information(signal, lag, n_bins) estimates average mutual information for delay selection. nearest_neighbor_distances(embedded) supports false-nearest-neighbor dimension selection.

Direct backend boundary

Go, Julia, and Mojo direct bridge calls share the same typed pre-dispatch contract before optional runtime loading. Signal payloads must be finite real one-dimensional float64 arrays and must reject numeric-string aliases before float coercion. Delay, dimension, lag, bin-count, row-count, and embedding-dimension controls must be integer values in the public API domain. Embedded nearest-neighbor payloads must be finite real flat float64 arrays whose length matches T*m, with numeric-string aliases rejected at the same pre-dispatch boundary.

Boolean aliases, object-dtype complex aliases, complex samples, non-finite samples, numeric-string aliases, non-vector signal payloads, malformed flattened embedding lengths, invalid delay/dimension requests, invalid lags, and invalid bin counts are rejected before shared-library, Julia, or subprocess execution.

Direct backend return payloads are validated before they are handed back to the public monitor boundary. The public dispatcher repeats the same physics-facing checks after backend fallback resolution so a shape-correct optional backend cannot silently return the wrong phase-space reconstruction. Delay-embedding outputs must have the exact (T_effective, dimension) shape and match the mathematical indexing x[t + k*tau]; numeric-string aliases and object-dtype complex aliases are rejected before float coercion; mutual-information outputs must be finite non-negative real scalars; nearest-neighbor outputs must contain finite non-negative distances and integral in-range neighbor indices, with self-neighbors rejected for non-trivial embeddings. Malformed Mojo text output is normalised to deterministic ValueError failures rather than leaking parser exceptions.

Invariants

delay_embed is exact indexing and should match across backends without tolerance. Mutual information is non-negative. Nearest-neighbor distances are finite non-negative values with integer neighbor indices in range and no self-neighbor for non-trivial inputs. Optional backend outputs that violate those invariants are rejected and the dispatcher falls back to the next available backend rather than returning a corrupted embedding or a truncated neighbor index.

Practical usage profile

Teams typically use this module during inspection and replay pipelines:

• validate embedding settings with optimal_delay/optimal_dimension,
• extract trajectory geometry with delay coordinates,
• use downstream monitors on the reconstructed state space instead of direct raw samples.

That pattern keeps signal reconstruction and control logic in one audited path.

embedding

Delay-embedding analysis with a 5-backend fallback chain.

Three compute primitives on the multi-language chain:

• :func:delay_embed — time-delay embedding matrix.
• :func:mutual_information — Fraser-Swinney 1986 average mutual information.
• :func:nearest_neighbor_distances — brute-force k=1 kNN in the embedded space (consumed by FNN).

Two wrappers stay Python-side (they are control flow over the primitives):

• :func:optimal_delay — first local minimum of MI (Fraser-Swinney).
• :func:optimal_dimension — Kennel-Brown-Abarbanel 1992 FNN.
• :func:auto_embed — convenience that chains optimal_delay, optimal_dimension, and :func:delay_embed.

The Rust backend exposes native optimal_delay_rust and optimal_dimension_rust entry points; when Rust is active those wrappers use the native path for maximum throughput. The Python fallback composes the primitives through the dispatcher.

MI and NN are exposed by Julia / Go / Mojo / Python only — Rust does not expose standalone MI or kNN FFI; those slots dispatch to the next available backend in the chain.

Classes

EmbeddingResult dataclass

EmbeddingResult(
    trajectory: FloatArray,
    delay: int,
    dimension: int,
    T_effective: int,
)

Delay-embedding output.

Methods:

__post_init__

__post_init__() -> None

Validate and normalise the embedded trajectory record.

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

def __post_init__(self) -> None:
    """Validate and normalise the embedded trajectory record."""
    trajectory = _validate_embedded(self.trajectory)
    delay = _validate_int_at_least(self.delay, name="delay", minimum=1)
    dimension = _validate_int_at_least(self.dimension, name="dimension", minimum=1)
    t_effective = _validate_int_at_least(
        self.T_effective,
        name="T_effective",
        minimum=0,
    )
    if trajectory.shape != (t_effective, dimension):
        raise ValueError(
            f"trajectory shape {trajectory.shape} does not match "
            f"(T_effective={t_effective}, dimension={dimension})"
        )
    self.trajectory = trajectory
    self.delay = delay
    self.dimension = dimension
    self.T_effective = t_effective

Functions:

delay_embed

delay_embed(
    signal: object, delay: object, dimension: object
) -> FloatArray

Time-delay embedding: v(t) = [x(t), x(t+τ), x(t+2τ), …].

Parameters

signal : object Real-valued time series, shape (T,). delay : object Embedding delay τ in samples. dimension : object Embedding dimension.

Returns

FloatArray The time-delay embedding, shape (M, dimension).

Raises

ValueError If delay or dimension is non-positive or too large for the signal.

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

def delay_embed(
    signal: object,
    delay: object,
    dimension: object,
) -> FloatArray:
    """Time-delay embedding: ``v(t) = [x(t), x(t+τ), x(t+2τ), …]``.

    Parameters
    ----------
    signal : object
        Real-valued time series, shape ``(T,)``.
    delay : object
        Embedding delay ``τ`` in samples.
    dimension : object
        Embedding dimension.

    Returns
    -------
    FloatArray
        The time-delay embedding, shape ``(M, dimension)``.

    Raises
    ------
    ValueError
        If ``delay`` or ``dimension`` is non-positive or too large for the signal.
    """
    s = _validate_signal(signal)
    delay = _validate_int_at_least(delay, name="delay", minimum=1)
    dimension = _validate_int_at_least(dimension, name="dimension", minimum=1)
    t_eff = int(s.size) - (dimension - 1) * delay
    if t_eff <= 0:
        msg = (
            f"Signal too short (T={s.size}) for delay={delay}, "
            f"dimension={dimension}: need T > {(dimension - 1) * delay}"
        )
        raise ValueError(msg)

    backend_fn = _dispatch("de")
    if backend_fn is not None:
        fn = cast("Callable[[FloatArray, int, int], FloatArray]", backend_fn)
        try:
            return _validate_delay_embedding_output(
                fn(s, delay, dimension),
                signal=s,
                delay=delay,
                t_effective=t_eff,
                dimension=dimension,
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            backend_fn = None

    indices = np.arange(dimension) * delay
    rows = np.arange(t_eff)[:, np.newaxis] + indices[np.newaxis, :]
    trajectory: FloatArray = np.asarray(s[rows], dtype=np.float64)
    return _validate_delay_embedding_output(
        trajectory,
        signal=s,
        delay=delay,
        t_effective=t_eff,
        dimension=dimension,
    )

mutual_information

mutual_information(
    signal: object, lag: object, n_bins: object = 32
) -> float

Fraser-Swinney 1986 average mutual information at lag.

Parameters

signal : object Real-valued time series, shape (T,). lag : object Lag in samples. n_bins : object Number of histogram bins.

Returns

float The average mutual information at the given lag.

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

def mutual_information(
    signal: object,
    lag: object,
    n_bins: object = 32,
) -> float:
    """Fraser-Swinney 1986 average mutual information at ``lag``.

    Parameters
    ----------
    signal : object
        Real-valued time series, shape ``(T,)``.
    lag : object
        Lag in samples.
    n_bins : object
        Number of histogram bins.

    Returns
    -------
    float
        The average mutual information at the given lag.
    """
    s = _validate_signal(signal)
    lag = _validate_int_at_least(lag, name="lag", minimum=0)
    n_bins = _validate_int_at_least(n_bins, name="n_bins", minimum=2)
    if s.size - lag <= 0:
        return 0.0

    backend_fn = _dispatch("mi")
    if backend_fn is not None:
        fn = cast("Callable[[FloatArray, int, int], float]", backend_fn)
        try:
            return _validate_non_negative_scalar(
                fn(s, lag, n_bins),
                name="mutual_information",
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            backend_fn = None

    t_total = s.size - lag
    x = s[:t_total]
    y = s[lag : lag + t_total]
    hist_xy, _, _ = np.histogram2d(x, y, bins=n_bins)
    hist_x = hist_xy.sum(axis=1)
    hist_y = hist_xy.sum(axis=0)
    total = hist_xy.sum()
    if total <= 0:
        return 0.0
    p_xy = hist_xy / total
    p_x = hist_x / total
    p_y = hist_y / total
    mi = 0.0
    for i in range(n_bins):
        for j in range(n_bins):
            if p_xy[i, j] > 0 and p_x[i] > 0 and p_y[j] > 0:
                mi += p_xy[i, j] * np.log(p_xy[i, j] / (p_x[i] * p_y[j]))
    return _validate_non_negative_scalar(mi, name="mutual_information")

nearest_neighbor_distances

nearest_neighbor_distances(
    embedded: object,
) -> tuple[FloatArray, IntArray]

Brute-force k = 1 kNN on the rows of embedded.

Parameters

embedded : object Delay-embedded trajectory, shape (M, dimension).

Returns

tuple[FloatArray, IntArray] The nearest-neighbour distances and their indices.

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

def nearest_neighbor_distances(
    embedded: object,
) -> tuple[FloatArray, IntArray]:
    """Brute-force ``k = 1`` kNN on the rows of ``embedded``.

    Parameters
    ----------
    embedded : object
        Delay-embedded trajectory, shape ``(M, dimension)``.

    Returns
    -------
    tuple[FloatArray, IntArray]
        The nearest-neighbour distances and their indices.
    """
    e = _validate_embedded(embedded)
    t, m = int(e.shape[0]), int(e.shape[1])
    if t == 0:
        return np.zeros(0, dtype=np.float64), np.zeros(0, dtype=np.int64)

    backend_fn = _dispatch("nn")
    if backend_fn is not None:
        fn = cast(
            "Callable[[FloatArray, int, int], tuple[FloatArray, IntArray]]",
            backend_fn,
        )
        try:
            dist, idx = fn(np.ascontiguousarray(e.ravel(), dtype=np.float64), t, m)
            return _validate_nn_output(dist, idx, n_points=t)
        except (ImportError, RuntimeError, OSError, KeyError):
            backend_fn = None

    nn_dist = np.full(t, np.inf)
    nn_idx = np.zeros(t, dtype=np.int64)
    for i in range(t):
        diffs = e - e[i]
        dists = np.sqrt(np.sum(diffs**2, axis=1))
        dists[i] = np.inf
        j = int(np.argmin(dists))
        nn_dist[i] = dists[j]
        nn_idx[i] = j
    return _validate_nn_output(nn_dist, nn_idx, n_points=t)

optimal_delay

optimal_delay(
    signal: object,
    max_lag: object = 100,
    n_bins: object = 32,
) -> int

First local minimum of :func:mutual_information vs lag.

Parameters

signal : object Real-valued time series, shape (T,). max_lag : object Largest lag to search. n_bins : object Number of histogram bins.

Returns

int The first mutual-information minimum, as a lag in samples.

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

def optimal_delay(
    signal: object,
    max_lag: object = 100,
    n_bins: object = 32,
) -> int:
    """First local minimum of :func:`mutual_information` vs ``lag``.

    Parameters
    ----------
    signal : object
        Real-valued time series, shape ``(T,)``.
    max_lag : object
        Largest lag to search.
    n_bins : object
        Number of histogram bins.

    Returns
    -------
    int
        The first mutual-information minimum, as a lag in samples.
    """
    s = _validate_signal(signal)
    max_lag = _validate_int_at_least(max_lag, name="max_lag", minimum=1)
    n_bins = _validate_int_at_least(n_bins, name="n_bins", minimum=2)

    if ACTIVE_BACKEND == "rust":
        try:
            fn = cast(
                "Callable[[FloatArray, int, int], int]",
                _load_backend("rust")["optimal_delay"],
            )
            return int(fn(s, max_lag, n_bins))
        except (ImportError, RuntimeError, OSError, KeyError):
            max_lag = int(max_lag)

    max_lag = min(max_lag, s.size // 2)
    mi_values = np.array([mutual_information(s, lag, n_bins) for lag in range(max_lag)])
    for i in range(1, len(mi_values) - 1):
        if mi_values[i] < mi_values[i - 1] and mi_values[i] < mi_values[i + 1]:
            return i
    return 1

optimal_dimension

optimal_dimension(
    signal: object,
    delay: object,
    max_dim: object = 10,
    rtol: object = 15.0,
    atol: object = 2.0,
) -> int

Kennel-Brown-Abarbanel 1992 FNN to select embedding dimension.

Parameters

signal : object Real-valued time series, shape (T,). delay : object Embedding delay τ in samples. max_dim : object Largest embedding dimension to test. rtol : object Relative tolerance for the false-nearest-neighbour test. atol : object Absolute tolerance for the false-nearest-neighbour test.

Returns

int The selected embedding dimension.

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

def optimal_dimension(
    signal: object,
    delay: object,
    max_dim: object = 10,
    rtol: object = 15.0,
    atol: object = 2.0,
) -> int:
    """Kennel-Brown-Abarbanel 1992 FNN to select embedding dimension.

    Parameters
    ----------
    signal : object
        Real-valued time series, shape ``(T,)``.
    delay : object
        Embedding delay ``τ`` in samples.
    max_dim : object
        Largest embedding dimension to test.
    rtol : object
        Relative tolerance for the false-nearest-neighbour test.
    atol : object
        Absolute tolerance for the false-nearest-neighbour test.

    Returns
    -------
    int
        The selected embedding dimension.
    """
    s = _validate_signal(signal)
    delay = _validate_int_at_least(delay, name="delay", minimum=1)
    max_dim = _validate_int_at_least(max_dim, name="max_dim", minimum=1)
    rtol = _validate_non_negative_real(rtol, name="rtol")
    atol = _validate_non_negative_real(atol, name="atol")
    sigma = float(np.std(s))
    if sigma == 0:
        return 1

    if ACTIVE_BACKEND == "rust":
        try:
            fn = cast(
                "Callable[..., int]",
                _load_backend("rust")["optimal_dimension"],
            )
            return int(
                fn(s, delay, max_dim, rtol, atol),
            )
        except (ImportError, RuntimeError, OSError, KeyError):
            max_dim = int(max_dim)

    for m in range(1, max_dim + 1):
        t_next = s.size - m * delay
        if t_next <= 1:
            return m
        emb_m = delay_embed(s, delay, m)
        t_m = emb_m.shape[0]
        nn_dist, nn_idx = nearest_neighbor_distances(emb_m)

        n_false = 0
        n_valid = 0
        for i in range(t_m):
            j = int(nn_idx[i])
            d = nn_dist[i]
            if d == 0 or not np.isfinite(d):
                continue
            i_next = i + m * delay
            j_next = j + m * delay
            if i_next >= s.size or j_next >= s.size:
                continue
            n_valid += 1
            extra = abs(s[i_next] - s[j_next])
            if extra / d > rtol:
                n_false += 1
                continue
            new_dist = (d * d + extra * extra) ** 0.5
            if new_dist / sigma > atol:
                n_false += 1
        fnn_frac = n_false / n_valid if n_valid > 0 else 0.0
        if fnn_frac < 0.01:
            return m
    return max_dim

auto_embed

auto_embed(
    signal: object,
    max_lag: object = 100,
    max_dim: object = 10,
) -> EmbeddingResult

optimal_delay ∘ optimal_dimension ∘ delay_embed.

Parameters

signal : object Real-valued time series, shape (T,). max_lag : object Largest lag to search. max_dim : object Largest embedding dimension to test.

Returns

EmbeddingResult The auto-selected delay/dimension embedding result.

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

def auto_embed(
    signal: object,
    max_lag: object = 100,
    max_dim: object = 10,
) -> EmbeddingResult:
    """``optimal_delay`` ∘ ``optimal_dimension`` ∘ ``delay_embed``.

    Parameters
    ----------
    signal : object
        Real-valued time series, shape ``(T,)``.
    max_lag : object
        Largest lag to search.
    max_dim : object
        Largest embedding dimension to test.

    Returns
    -------
    EmbeddingResult
        The auto-selected delay/dimension embedding result.
    """
    tau = optimal_delay(signal, max_lag)
    m = optimal_dimension(signal, tau, max_dim)
    traj = delay_embed(signal, tau, m)
    return EmbeddingResult(
        trajectory=traj,
        delay=tau,
        dimension=m,
        T_effective=int(traj.shape[0]),
    )