Reference

Derived systems¶

Wrappers that re-present an existing system through a new lens while keeping the System protocol intact, so analysis functions compose with them transparently. Prose introduction: the mental model.

PoincareMap ¶

PoincareMap(
    system: Any,
    plane: tuple[Any, ...],
    *,
    direction: int | str = +1,
    dt: float = 0.01,
    max_time: float = 10000.0,
)

Bases: DerivedSystem

Present a flow as the discrete map of its crossings through a hyperplane.

One step() advances the underlying system until the trajectory crosses the section plane in the chosen direction, refines the crossing point by cubic Hermite interpolation of the bracketing samples (using the system's numeric RHS for endpoint derivatives — O(dt⁴) accuracy), and returns the full-dimensional crossing state.

Because a PoincareMap is a discrete system, everything written for maps applies to flows through it — e.g. an orbit diagram over a PoincareMap is a bifurcation diagram of the flow.

PARAMETER	DESCRIPTION
`system`	A continuous-time system (ODE or DDE). TYPE: `System`
`plane`	The section, in any of three spellings: `(axis, c)` where `axis` is a component name (resolved against the system's `variables`) or an integer index, for the section `y_axis = c`; `(axis, c, direction)`, the same with the crossing direction (`"up"` / `"down"` / `"both"`) as a third element; or `(normal, offset)` with an arbitrary normal vector, for the section `normal · y = offset`. Examples: `plane=("y", 0.0)`, `plane=("y", 0.0, "up")`, `plane=(1, 0.0)`, `plane=([1, 0, 0], 0.0)`. TYPE: `tuple`
`direction`	Count only crossings with `d(normal·y)/dt > 0` (`+1` / `"up"`, default), `< 0` (`-1` / `"down"`), or both (`0` / `"both"`). A direction given inside `plane` (its third element) overrides this. TYPE: `+1, -1, 0} or {"up", "down", "both"` DEFAULT: `+1`
`dt`	March step used for crossing detection. The refinement makes the crossing itself far more accurate than `dt`; this only needs to be small enough not to skip crossings. TYPE: `float` DEFAULT: `0.01`
`max_time`	Raise if no crossing is found within this much time (e.g. the plane misses the attractor). TYPE: `float` DEFAULT: `10000.0`

Examples:

>>> pmap = PoincareMap(Rossler(), plane=("x", 0.0, "up"))
>>> section = pmap.trajectory(500)         # 500 crossings → PoincareSection
>>> section.y.shape
(500, 3)

Source code in src/tsdynamics/derived/poincare.py

def __init__(
    self,
    system: Any,
    plane: tuple[Any, ...],
    *,
    direction: int | str = +1,
    dt: float = 0.01,
    max_time: float = 1e4,
) -> None:
    super().__init__(system)
    plane, direction = _resolve_section_plane(system, plane, direction)
    normal, offset = self._parse_plane(system.dim, plane)
    self.plane = plane
    self._normal = normal
    self._offset = offset
    self.direction = direction
    self.dt = float(dt)
    self.max_time = float(max_time)

    # Numeric RHS for Hermite endpoint derivatives; falls back to linear
    # interpolation (O(dt²)) for systems without one (e.g. DDEs).
    self._rhs = system._rhs_numeric() if hasattr(system, "_rhs_numeric") else None

    self._u_cross: np.ndarray | None = None
    self._t_cross: float | None = None
    self._n_cross = 0

is_discrete `property` ¶

is_discrete: bool

A Poincaré map is a discrete view of the flow.

crossing_count `property` ¶

crossing_count: int

Return the number of crossings recorded so far.

step ¶

step(n_or_dt: int | None = None) -> ndarray

Advance to the n-th next crossing and return it (full-dim coords).

Source code in src/tsdynamics/derived/poincare.py

def step(self, n_or_dt: int | None = None) -> np.ndarray:
    """Advance to the ``n``-th next crossing and return it (full-dim coords)."""
    n = int(n_or_dt) if n_or_dt is not None else 1
    for _ in range(n):
        self._advance_to_crossing()
    assert self._u_cross is not None  # set by _advance_to_crossing (n >= 1)
    return self._u_cross.copy()

state ¶

state() -> ndarray

Return the last crossing point (or the inner state before any crossing).

Source code in src/tsdynamics/derived/poincare.py

def state(self) -> np.ndarray:
    """Return the last crossing point (or the inner state before any crossing)."""
    if self._u_cross is not None:
        return self._u_cross.copy()
    return cast(np.ndarray, self.system.state())

set_state ¶

set_state(u: Any) -> None

Overwrite the inner flow state and reset crossing bookkeeping.

Source code in src/tsdynamics/derived/poincare.py

def set_state(self, u: Any) -> None:
    """Overwrite the inner flow state and reset crossing bookkeeping."""
    self.system.set_state(u)
    self._u_cross = None
    self._t_cross = None

time ¶

time() -> float

Return the continuous time of the last crossing (inner time before any).

Source code in src/tsdynamics/derived/poincare.py

def time(self) -> float:
    """Return the continuous time of the last crossing (inner time before any)."""
    return self._t_cross if self._t_cross is not None else self.system.time()

as_events ¶

as_events() -> list[Event]

Return the section as a one-element [Event] for system.run(events=...).

A Poincaré section is an event: the crossing of g(u) = normal·u − offset in the map's :attr:direction. This exposes it as an :class:~tsdynamics.engine.run.Event so the general events= API reproduces the section — PoincareMap is one consumer of the same wired engine seam (stream WS-EVENTSAPI / WS-CROSSKERNEL). Driven at the same fixed-step march (method="rk4" at this map's dt) from the same initial condition, the crossings of inner.run(events=pmap.as_events(), ...) match :meth:trajectory to the engine's refinement accuracy.

Examples:

>>> pmap = PoincareMap(Rossler(), plane=("y", 0.0, "up"), dt=0.01)
>>> sol = Rossler().run(final_time=400, dt=0.01, method="rk4",
...                     events=pmap.as_events())
>>> sol.meta["y_events"][0][:5].shape      # the same crossing states
(5, 3)

Source code in src/tsdynamics/derived/poincare.py

def as_events(self) -> list[Event]:
    """Return the section as a one-element ``[Event]`` for ``system.run(events=...)``.

    A Poincaré section *is* an event: the crossing of ``g(u) = normal·u −
    offset`` in the map's :attr:`direction`.  This exposes it as an
    :class:`~tsdynamics.engine.run.Event` so the general ``events=`` API
    reproduces the section — ``PoincareMap`` is one consumer of the same
    wired engine seam (stream WS-EVENTSAPI / WS-CROSSKERNEL).  Driven at the
    same fixed-step march (``method="rk4"`` at this map's ``dt``) from the
    same initial condition, the crossings of
    ``inner.run(events=pmap.as_events(), ...)`` match
    :meth:`trajectory` to the engine's refinement accuracy.

    Examples
    --------
    >>> pmap = PoincareMap(Rossler(), plane=("y", 0.0, "up"), dt=0.01)
    >>> sol = Rossler().run(final_time=400, dt=0.01, method="rk4",
    ...                     events=pmap.as_events())
    >>> sol.meta["y_events"][0][:5].shape      # the same crossing states
    (5, 3)
    """
    from tsdynamics.engine.run import Event

    return [Event((self._normal.copy(), self._offset), direction=self.direction)]

reinit ¶

reinit(u: Any | None = None, **kwargs: Any) -> None

Restart the inner flow and clear crossing bookkeeping.

Source code in src/tsdynamics/derived/poincare.py

def reinit(self, u: Any | None = None, **kwargs: Any) -> None:
    """Restart the inner flow and clear crossing bookkeeping."""
    self.system.reinit(u, **kwargs)
    # Parameter values are baked into the numeric RHS — rebuild it so the
    # Hermite refinement matches the (possibly re-parametrized) dynamics.
    if hasattr(self.system, "_rhs_numeric"):
        self._rhs = self.system._rhs_numeric()
    self._u_cross = None
    self._t_cross = None
    self._n_cross = 0

trajectory ¶

trajectory(
    steps: int = 100,
    *,
    transient: int = 0,
    backend: str | None = None,
    **kwargs: Any,
) -> PoincareSection

Collect crossings as a :class:PoincareSection.

t holds the continuous crossing times; y the full-dimensional crossing states. transient crossings are discarded first. The returned :class:PoincareSection is a :class:~tsdynamics.data.Trajectory carrying section intent (so a renderer draws the in-plane scatter) plus a .summary() / .to_dict() readout.

For an ordinary (non-stiff) ODE on the compiled engine this marches the whole attractor and refines every crossing in one engine call (the wired Rust integrate_events, stream WS-CROSSKERNEL) — ~100× faster than the per-dt Python loop it replaces. DDEs, systems without a numeric RHS, stiff defaults, and backend="reference" keep the Python loop. The engine path is answer-identical to that loop's fixed-step (rk4) refinement; see :mod:tsdynamics.derived._crossings.

PARAMETER	DESCRIPTION
`steps`	Number of crossings to collect. TYPE: `int` DEFAULT: `100`
`transient`	Number of leading crossings to discard. TYPE: `int` DEFAULT: `0`
`backend`	Engine evaluator for the fast path; defaults to the inner system's backend. `"reference"` forces the pure-Python loop. TYPE: `('interp', 'jit', 'reference')` DEFAULT: `"interp"`

RETURNS	DESCRIPTION
`PoincareSection`	The collected crossings (continuous times in `t`, full-dimensional states in `y`), carrying section plot intent and a `.summary()` / `.to_dict()` readout.

RAISES	DESCRIPTION
`ConvergenceError`	If the inner flow diverges (non-finite state) or no crossing is found within `max_time` of marching — the plane misses the attractor, or the crossing :attr:`direction` is wrong.

Notes

Live-cursor semantics. trajectory advances the inner system as a side effect, so a subsequent :meth:step continues forward rather than re-yielding the crossings just collected. The two collection paths leave the inner cursor at slightly different places: the Python loop stops just past the last collected crossing, while the engine path stops at the span end it marched to (which can be a little beyond the last crossing). Both invariants — that step() resumes after the collected crossings — hold; do not rely on the exact cursor offset. Call :meth:reinit first if you need a deterministic restart point.

Source code in src/tsdynamics/derived/poincare.py

def trajectory(
    self,
    steps: int = 100,
    *,
    transient: int = 0,
    backend: str | None = None,
    **kwargs: Any,
) -> PoincareSection:
    """
    Collect crossings as a :class:`PoincareSection`.

    ``t`` holds the continuous crossing times; ``y`` the full-dimensional
    crossing states.  ``transient`` crossings are discarded first.  The
    returned :class:`PoincareSection` is a :class:`~tsdynamics.data.Trajectory`
    carrying section intent (so a renderer draws the in-plane scatter) plus a
    ``.summary()`` / ``.to_dict()`` readout.

    For an ordinary (non-stiff) ODE on the compiled engine this marches the
    whole attractor and refines every crossing in **one engine call** (the
    wired Rust ``integrate_events``, stream WS-CROSSKERNEL) — ~100× faster than
    the per-``dt`` Python loop it replaces.  DDEs, systems without a numeric
    RHS, stiff defaults, and ``backend="reference"`` keep the Python loop.  The
    engine path is answer-identical to that loop's fixed-step (``rk4``)
    refinement; see :mod:`tsdynamics.derived._crossings`.

    Parameters
    ----------
    steps : int
        Number of crossings to collect.
    transient : int
        Number of leading crossings to discard.
    backend : {"interp", "jit", "reference"}, optional
        Engine evaluator for the fast path; defaults to the inner system's
        backend.  ``"reference"`` forces the pure-Python loop.

    Returns
    -------
    PoincareSection
        The collected crossings (continuous times in ``t``, full-dimensional
        states in ``y``), carrying section plot intent and a
        ``.summary()`` / ``.to_dict()`` readout.

    Raises
    ------
    ConvergenceError
        If the inner flow diverges (non-finite state) or no crossing is found
        within ``max_time`` of marching — the plane misses the attractor, or
        the crossing :attr:`direction` is wrong.

    Notes
    -----
    **Live-cursor semantics.**  ``trajectory`` advances the *inner* system as
    a side effect, so a subsequent :meth:`step` continues forward rather than
    re-yielding the crossings just collected.  The two collection paths leave
    the inner cursor at slightly different places: the Python loop stops just
    past the **last collected** crossing, while the engine path stops at the
    **span end** it marched to (which can be a little beyond the last
    crossing).  Both invariants — that ``step()`` resumes *after* the
    collected crossings — hold; do not rely on the exact cursor offset.  Call
    :meth:`reinit` first if you need a deterministic restart point.
    """
    if kwargs:
        self.reinit(kwargs.pop("ic", None), **kwargs)

    if _crossings.engine_eligible(self.system, backend):
        from tsdynamics.engine.run import EngineNotAvailableError

        try:
            times, points = self._engine_trajectory(steps, transient, backend)
        except EngineNotAvailableError:
            times, points = self._python_trajectory(steps, transient)
    else:
        times, points = self._python_trajectory(steps, transient)

    meta = {
        "derived": "PoincareMap",
        # Section intent, so a renderer draws the 2-D in-plane scatter rather
        # than mistaking the full-dimensional crossing states for a flow line
        # (the string value of viz.PlotKind.POINCARE_SECTION).
        "plot_kind": "poincare_section",
        "plane": self.plane,
        "direction": self.direction,
        "dt": self.dt,
        "system": type(self.system).__name__,
        "params": self.params.as_dict(),
    }
    return PoincareSection(t=times, y=points, system=self.system, meta=meta)

StroboscopicMap ¶

StroboscopicMap(system: Any, period: float)

Bases: DerivedSystem

Present a forced flow as the discrete map of once-per-period samples.

One step() advances the underlying continuous system by exactly one forcing period and returns the new state. Orbit diagrams over a StroboscopicMap are the standard way to study forced oscillators (Duffing, forced van der Pol, ...).

PARAMETER	DESCRIPTION
`system`	A continuous-time system. TYPE: `System`
`period`	Sampling period (the forcing period). TYPE: `float`

Examples:

>>> smap = StroboscopicMap(ForcedVanDerPol(), period=2 * np.pi / 0.63)
>>> samples = smap.trajectory(300, transient=100)

Source code in src/tsdynamics/derived/stroboscopic.py

def __init__(self, system: Any, period: float) -> None:
    super().__init__(system)
    if period <= 0:
        raise ValueError(f"period must be positive, got {period}")
    self.period = float(period)

is_discrete `property` ¶

is_discrete: bool

A stroboscopic map is a discrete view of the flow.

step ¶

step(n_or_dt: int | None = None) -> ndarray

Advance n forcing periods (default 1) and return the new state.

n_or_dt is a period count, not a time increment — the wrapper presents a discrete map, so the argument is coerced to an integer with int() (a float is truncated toward zero, matching the discrete-view convention; pass a whole number to be explicit).

PARAMETER	DESCRIPTION
`n_or_dt`	Number of forcing periods to advance. `None` advances one period. TYPE: `int` DEFAULT: `None`

RETURNS	DESCRIPTION
`ndarray`	The full-dimensional state after `n` periods.

Source code in src/tsdynamics/derived/stroboscopic.py

def step(self, n_or_dt: int | None = None) -> np.ndarray:
    """Advance ``n`` forcing periods (default 1) and return the new state.

    ``n_or_dt`` is a **period count**, not a time increment — the wrapper
    presents a discrete map, so the argument is coerced to an integer with
    ``int()`` (a float is *truncated toward zero*, matching the discrete-view
    convention; pass a whole number to be explicit).

    Parameters
    ----------
    n_or_dt : int, optional
        Number of forcing periods to advance.  ``None`` advances one period.

    Returns
    -------
    numpy.ndarray
        The full-dimensional state after ``n`` periods.
    """
    n = int(n_or_dt) if n_or_dt is not None else 1
    return cast(np.ndarray, self.system.step(n * self.period))

time ¶

time() -> float

Return the inner flow time.

Source code in src/tsdynamics/derived/stroboscopic.py

def time(self) -> float:
    """Return the inner flow time."""
    return cast(float, self.system.time())

trajectory ¶

trajectory(
    steps: int = 100, *, transient: int = 0, **kwargs: Any
) -> Trajectory

Collect steps once-per-period samples (after transient periods).

Sampling starts from the inner flow's live cursor and advances by exactly one forcing :attr:period per sample; transient leading periods are stepped through and discarded first. Any keyword (ic=, solver options) triggers a :meth:reinit before sampling.

PARAMETER	DESCRIPTION
`steps`	Number of once-per-period samples to collect. TYPE: `int` DEFAULT: `100`
`transient`	Number of leading periods to step through and discard. TYPE: `int` DEFAULT: `0`
`**kwargs`	Forwarded to :meth:`reinit` (`ic` is popped) when non-empty. TYPE: `Any` DEFAULT: `{}`

RETURNS	DESCRIPTION
`Trajectory`	The strobed samples — continuous sample times in `t` (one period apart), full-dimensional states in `y`.

Source code in src/tsdynamics/derived/stroboscopic.py

def trajectory(self, steps: int = 100, *, transient: int = 0, **kwargs: Any) -> Trajectory:
    """Collect ``steps`` once-per-period samples (after ``transient`` periods).

    Sampling starts from the inner flow's **live cursor** and advances by
    exactly one forcing :attr:`period` per sample; ``transient`` leading
    periods are stepped through and discarded first.  Any keyword (``ic=``,
    solver options) triggers a :meth:`reinit` before sampling.

    Parameters
    ----------
    steps : int
        Number of once-per-period samples to collect.
    transient : int
        Number of leading periods to step through and discard.
    **kwargs
        Forwarded to :meth:`reinit` (``ic`` is popped) when non-empty.

    Returns
    -------
    Trajectory
        The strobed samples — continuous sample times in ``t`` (one period
        apart), full-dimensional states in ``y``.
    """
    if kwargs:
        self.reinit(kwargs.pop("ic", None), **kwargs)
    if transient:
        self.system.step(transient * self.period)
    times = np.empty(steps)
    points = np.empty((steps, self.system.dim))
    for k in range(steps):
        points[k] = self.step()
        times[k] = self.system.time()
    meta = {
        "derived": "StroboscopicMap",
        "period": self.period,
        "system": type(self.system).__name__,
        "params": self.params.as_dict(),
    }
    return Trajectory(t=times, y=points, system=self.system, meta=meta)

to_plot_spec ¶

to_plot_spec(
    kind: str | None = None, *, steps: int = 300
) -> PlotSpec

Describe the strobe sampling as a scatter of sampled states.

A stroboscopic map is a discrete sampling — once per forcing period — so the natural picture is a cloud of sampled points (the strobed orbit / attractor), not a connected flow line. This collects steps samples and builds a 2-D / 3-D SCATTER spec over the first two / three components (a 1-D system is a sample-index time series of dots).

The :mod:tsdynamics.viz.spec import is lazy, so building a spec never pulls in a plotting backend.

PARAMETER	DESCRIPTION
`kind`	Override the auto-dispatched semantic kind (e.g. `"phase_portrait_2d"`). `None` (the default) dispatches on the sampled dimensionality. TYPE: `str` DEFAULT: `None`
`steps`	Number of once-per-period samples to collect. Default `300`. TYPE: `int` DEFAULT: `300`

RETURNS	DESCRIPTION
`PlotSpec`

Notes

Sampling starts from the inner flow's live cursor (this calls :meth:trajectory, which steps the wrapped system as a side effect with no transient discarded), so the picture reflects wherever the system currently sits — :meth:reinit first for a deterministic start state, or burn the transient in beforehand to image the attractor rather than the approach to it.

Source code in src/tsdynamics/derived/stroboscopic.py

def to_plot_spec(self, kind: str | None = None, *, steps: int = 300) -> PlotSpec:
    """Describe the strobe sampling as a **scatter** of sampled states.

    A stroboscopic map is a *discrete* sampling — once per forcing period —
    so the natural picture is a cloud of sampled points (the strobed orbit /
    attractor), **not** a connected flow line.  This collects ``steps``
    samples and builds a 2-D / 3-D ``SCATTER`` spec over the first two / three
    components (a 1-D system is a sample-index time series of dots).

    The :mod:`tsdynamics.viz.spec` import is lazy, so building a spec never
    pulls in a plotting backend.

    Parameters
    ----------
    kind : str, optional
        Override the auto-dispatched semantic kind (e.g.
        ``"phase_portrait_2d"``).  ``None`` (the default) dispatches on the
        sampled dimensionality.
    steps : int, optional
        Number of once-per-period samples to collect.  Default ``300``.

    Returns
    -------
    PlotSpec

    Notes
    -----
    Sampling starts from the inner flow's **live cursor** (this calls
    :meth:`trajectory`, which steps the wrapped system as a side effect with
    no transient discarded), so the picture reflects wherever the system
    currently sits — :meth:`reinit` first for a deterministic start state, or
    burn the transient in beforehand to image the attractor rather than the
    approach to it.
    """
    from tsdynamics.viz.spec import Axis, Layer, PlotKind, PlotSpec

    section = self.trajectory(steps)
    names = self.variables or tuple(f"y{i}" for i in range(self.system.dim))
    title = f"Stroboscopic map — {type(self.system).__name__}"

    if self.system.dim == 1:
        spec_kind = PlotKind(kind) if kind is not None else PlotKind.TIME_SERIES
        return PlotSpec(
            kind=spec_kind,
            ndim=1,
            title=title,
            x=Axis(label="sample"),
            y=Axis(label=names[0]),
            layers=[
                Layer(
                    PlotKind.SCATTER,
                    {"x": np.arange(section.y.shape[0], dtype=float), "y": section.y[:, 0]},
                )
            ],
        )

    if kind is None:
        want_3d = self.system.dim >= 3
    else:
        want_3d = PlotKind(kind) == PlotKind.PHASE_PORTRAIT_3D
    spec_kind = (
        PlotKind(kind)
        if kind is not None
        else (PlotKind.PHASE_PORTRAIT_3D if want_3d else PlotKind.PHASE_PORTRAIT_2D)
    )
    cols: dict[str, np.ndarray] = {"x": section.y[:, 0], "y": section.y[:, 1]}
    z = None
    if want_3d:
        cols["z"] = section.y[:, 2]
        z = Axis(label=names[2])
    return PlotSpec(
        kind=spec_kind,
        ndim=3 if want_3d else 2,
        aspect="equal",
        title=title,
        x=Axis(label=names[0]),
        y=Axis(label=names[1]),
        z=z,
        layers=[Layer(PlotKind.SCATTER, cols)],
    )

TangentSystem ¶

TangentSystem(
    system: Any,
    k: int | None = None,
    *,
    backend: str | None = None,
)

Bases: DerivedSystem

Evolve a system together with k deviation (tangent) vectors.

Each step() advances the state and the deviation vectors, then QR-reorthonormalises; :meth:growths exposes the per-step logarithmic stretch factors log |diag R| and :meth:exponents their running time-average — the Lyapunov spectrum estimate. :meth:lyapunov_spectrum wraps that into the standard burn-in + time-averaged estimate, and is the single implementation every family's lyapunov_spectrum delegates to.

Implementation per family

Maps: pure NumPy — W ← J(x)·W with _jacobian evaluated at the pre-image (the correct tangent-map convention), then QR.
ODEs: the extended ODE (state ⊕ k tangent vectors, see :mod:tsdynamics.derived._variational) is lowered to an engine tape and integrated per step on the Rust engine (or the pure-Python reference oracle), then QR-reorthonormalised here. Select the variational backend with backend=: "interp" (default), "jit", or "reference".
DDEs: not supported — tangent dynamics of a DDE lives in an infinite-dimensional history space; use DelaySystem.lyapunov_spectrum (the engine DDE Lyapunov estimator).

PARAMETER	DESCRIPTION
`system`	The base system whose tangent dynamics to evolve. TYPE: `DiscreteMap or ContinuousSystem`
`k`	Number of deviation vectors (`1 ≤ k ≤ system.dim`). Defaults to the full state dimension. TYPE: `int` DEFAULT: `None`
`backend`	ODE variational backend (ignored for maps): `"interp"` (default), `"jit"`, or `"reference"`. TYPE: `str` DEFAULT: `None`

Examples:

>>> tang = TangentSystem(Henon(), k=2)
>>> tang.reinit([0.1, 0.1])
>>> for _ in range(5000):
...     tang.step()
>>> tang.exponents()        # ≈ [0.42, -1.62]

Source code in src/tsdynamics/derived/tangent.py

def __init__(self, system: Any, k: int | None = None, *, backend: str | None = None) -> None:
    if isinstance(system, DelaySystem):
        raise NotImplementedError(
            "TangentSystem does not support delay systems — their tangent space is the "
            "infinite-dimensional history space. Use DelaySystem.lyapunov_spectrum."
        )
    if isinstance(system, DiscreteMap):
        self._mode = "map"
    elif isinstance(system, ContinuousSystem):
        self._mode = "ode"
    else:
        raise TypeError(
            f"TangentSystem needs a DiscreteMap or ContinuousSystem, "
            f"got {type(system).__name__}."
        )
    super().__init__(system)
    self.k = int(k) if k is not None else system.dim
    if not 1 <= self.k <= system.dim:
        raise ValueError(f"k must be in [1, {system.dim}], got {self.k}")

    # Both maps and ODEs select among the same three backends: the compiled
    # engine (``interp``/``jit``) or the pure-Python ``reference`` oracle.  For
    # maps ``interp``/``jit`` run the Rust QR tangent-map kernel (stream
    # perf/map-lyapunov-kernel) and ``reference`` the pure-NumPy QR loop (also
    # the transparent fallback when the engine declines — a non-lowering
    # ``_step`` or an absent wheel); for ODEs they select the variational
    # integration backend as before.
    self._backend = (backend or "interp").lower()
    if self._backend not in _ENGINE_BACKENDS:
        raise ValueError(
            f"unknown {'map' if self._mode == 'map' else 'ODE'} tangent backend "
            f"{backend!r}; choose from {sorted(_ENGINE_BACKENDS)}."
        )

    self._W: np.ndarray | None = None  # (dim, k) deviation vectors (map mode)
    self._ext_tape: Any = None  # extended variational tape (ode engine mode)
    self._ext_tape_key: Any = None  # structural-param key the tape was built for
    self._ext_tape_arrays: Any = None  # the extended tape's engine wire arrays (cached)
    self._z: np.ndarray | None = None  # extended state (ode engine mode)
    self._t = 0.0
    # ODE integration options (engine mode), captured at reinit.
    self._method: str | None = None
    self._rtol = 1e-6
    self._atol = 1e-9
    self._integrator_kwargs: dict[str, Any] = {}
    # Resumable engine stepper amortisation (stream perf/lyapunov-stepper): when
    # the resolved kernel is *explicit* and the backend is the compiled engine
    # (``interp``/``jit``), the per-dt-chunk loop reuses one durable
    # ``OdeStepper`` handle (built once over the extended variational tape,
    # re-seeded after each QR via ``set_state``) instead of constructing a fresh
    # ``ODEProblem`` + calling ``run.integrate`` every chunk.  ``advance(dt)`` is
    # bit-for-bit identical to the per-dt ``integrate_dense`` it supersedes, so
    # the spectrum is unchanged.  ``make_ode_stepper`` rejects an implicit kernel
    # at construction, so a stiff base flow (an implicit ``_default_method``) and
    # the ``reference`` oracle keep the per-chunk ``run.integrate`` path.
    self._ode_stepper: Any = None  # the durable OdeStepper handle (explicit only)
    self._step_kernel: str | None = None  # canonical kernel name (e.g. "rk45")
    self._step_explicit_engine = False  # eligible for the stepper fast path?
    self._last_growths = np.zeros(self.k)
    self._sum_growths = np.zeros(self.k)
    self._elapsed = 0.0
    # The batch engine map kernel (:meth:`_map_spectrum_engine`) returns only
    # the time-averaged spectrum and the interval count, never the final
    # orthonormal frame nor the last interval's per-step stretch.  When that
    # path runs, the streaming accessors ``deviations()`` / ``growths()`` cannot
    # be reproduced and are marked stale (set True there, cleared by ``reinit``)
    # so they raise an honest error instead of returning the long-run average.
    self._map_engine_stale = False

is_discrete `property` ¶

is_discrete: bool

Match the wrapped system's time semantics.

reinit ¶

reinit(
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    **kwargs: Any,
) -> None

Restart state, deviation vectors, and accumulated growth sums.

Source code in src/tsdynamics/derived/tangent.py

def reinit(
    self,
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    **kwargs: Any,
) -> None:
    """Restart state, deviation vectors, and accumulated growth sums."""
    self._reset_accumulators()

    if self._mode == "map":
        self.system.reinit(u, t=t, params=params)
        self._W = np.eye(self.system.dim)[:, : self.k]
        return

    # ode mode — apply any parameter overrides to the inner system first.
    if params:
        for key, value in params.items():
            self.params[key] = value
    ic_arr = self.system.resolve_ic(u)
    t0 = float(t) if t is not None else 0.0
    self._t = t0
    self._method = kwargs.pop("method", None)
    self._rtol = kwargs.pop("rtol", 1e-6)
    self._atol = kwargs.pop("atol", 1e-9)
    self._integrator_kwargs = dict(kwargs)

    self._reinit_ode_engine(ic_arr)

step ¶

step(n_or_dt: float | None = None) -> ndarray

Advance state + deviation vectors and reorthonormalise.

For maps n_or_dt is the number of iterations (QR after each); for ODEs it is the time increment (default 0.1). Returns the new state.

Source code in src/tsdynamics/derived/tangent.py

def step(self, n_or_dt: float | None = None) -> np.ndarray:
    """
    Advance state + deviation vectors and reorthonormalise.

    For maps ``n_or_dt`` is the number of iterations (QR after each);
    for ODEs it is the time increment (default 0.1).
    Returns the new state.
    """
    if self._mode == "map":
        return self._step_map(int(n_or_dt) if n_or_dt is not None else 1)
    dt = float(n_or_dt) if n_or_dt is not None else self._ODE_DEFAULT_DT
    return self._step_ode_engine(dt)

state ¶

state() -> ndarray

Return a copy of the current base-system state.

Source code in src/tsdynamics/derived/tangent.py

def state(self) -> np.ndarray:
    """Return a copy of the current base-system state."""
    if self._mode == "map":
        return cast(np.ndarray, self.system.state())
    # ODE mode: the constructor guarantees an engine backend, so the extended
    # state ``self._z`` always carries the base state in its leading slots.
    if self._z is None:
        self.reinit()
    assert self._z is not None  # reinit() seeds the extended state in ODE mode
    return self._z[: self.system.dim].copy()

set_state ¶

set_state(u: Any) -> None

Overwrite the base state (map mode only — ODE tangent vectors would desync).

Source code in src/tsdynamics/derived/tangent.py

def set_state(self, u: Any) -> None:
    """Overwrite the base state (map mode only — ODE tangent vectors would desync)."""
    if self._mode == "map":
        self.system.set_state(u)
    else:
        raise NotImplementedError(
            "TangentSystem(ode).set_state is not supported — the tangent vectors "
            "would desynchronise. Use reinit(u)."
        )

time ¶

time() -> float

Return the current time / iteration count.

Source code in src/tsdynamics/derived/tangent.py

def time(self) -> float:
    """Return the current time / iteration count."""
    return self.system.time() if self._mode == "map" else self._t

deviations ¶

deviations() -> ndarray

Return the current orthonormal deviation vectors, shape (dim, k).

Available on every supported backend — maps and the ODE engine backends both carry the deviation matrix explicitly.

Source code in src/tsdynamics/derived/tangent.py

def deviations(self) -> np.ndarray:
    """Return the current orthonormal deviation vectors, shape ``(dim, k)``.

    Available on every supported backend — maps and the ODE engine
    backends both carry the deviation matrix explicitly.
    """
    if self._mode == "map":
        if self._map_engine_stale:
            raise RuntimeError(
                "deviations() is unavailable after a batch engine map "
                "lyapunov_spectrum: the Rust kernel returns only the averaged "
                "spectrum, not the final orthonormal frame. Call reinit() and "
                "step() to drive the streaming tangent frame explicitly, or use "
                "backend='reference' for the streaming pure-Python QR loop."
            )
        if self._W is None:
            raise RuntimeError("deviations() is available after reinit()")
        return self._W.copy()
    # ODE engine backend (interp/jit/reference): the constructor rejects any
    # other backend, so the deviation vectors are always carried explicitly.
    if self._z is None:
        raise RuntimeError("deviations() is available after reinit()")
    return split_extended(self._z, self.system.dim, self.k)[1]

growths ¶

growths() -> ndarray

Return the log stretch factors log|diag R| from the most recent step.

RAISES	DESCRIPTION
`RuntimeError`	In map mode, after a batch engine `lyapunov_spectrum` call: the Rust kernel returns only the time-averaged spectrum, not the last interval's per-step stretch, so the most-recent-step value cannot be reproduced. Drive the frame with `reinit()` + `step()` (or use `backend='reference'`) for streaming per-step growths.

Source code in src/tsdynamics/derived/tangent.py

def growths(self) -> np.ndarray:
    """Return the log stretch factors ``log|diag R|`` from the most recent step.

    Raises
    ------
    RuntimeError
        In map mode, after a *batch* engine ``lyapunov_spectrum`` call: the
        Rust kernel returns only the time-averaged spectrum, not the last
        interval's per-step stretch, so the most-recent-step value cannot be
        reproduced. Drive the frame with ``reinit()`` + ``step()`` (or use
        ``backend='reference'``) for streaming per-step growths.
    """
    if self._mode == "map" and self._map_engine_stale:
        raise RuntimeError(
            "growths() is unavailable after a batch engine map "
            "lyapunov_spectrum: the Rust kernel returns only the averaged "
            "spectrum, not the most recent step's log|diag R|. Call reinit() and "
            "step() to drive the streaming tangent frame explicitly, or use "
            "backend='reference' for the streaming pure-Python QR loop."
        )
    return self._last_growths.copy()

exponents ¶

exponents() -> ndarray

Return the running Lyapunov-spectrum estimate (accumulated growths / elapsed).

Source code in src/tsdynamics/derived/tangent.py

def exponents(self) -> np.ndarray:
    """Return the running Lyapunov-spectrum estimate (accumulated growths / elapsed)."""
    if self._elapsed == 0.0:
        return np.zeros(self.k)
    return self._sum_growths / self._elapsed

convergence ¶

convergence(
    steps: int = 2000,
    n_or_dt: float | None = None,
    *,
    ic: Any | None = None,
) -> tuple[ndarray, ndarray]

Record the running Lyapunov estimates as they converge.

Reinitialises the tangent frame and steps it steps times, capturing the running :meth:exponents estimate after every step. The estimates settle as the time-average accumulates — the curve a user inspects to judge whether a Lyapunov run has converged.

PARAMETER	DESCRIPTION
`steps`	Number of tangent steps to record. Default `2000`. TYPE: `int` DEFAULT: `2000`
`n_or_dt`	Per-step increment (iterations for a map, `dt` for an ODE). `None` uses the family default. TYPE: `float` DEFAULT: `None`
`ic`	Initial condition; `None` resolves the system's default. TYPE: `array - like` DEFAULT: `None`

RETURNS	DESCRIPTION
`(times, estimates)`	`times` shape `(steps,)`; `estimates` shape `(steps, k)` — the running estimate of each of the `k` exponents after each step.

Source code in src/tsdynamics/derived/tangent.py

def convergence(
    self,
    steps: int = 2000,
    n_or_dt: float | None = None,
    *,
    ic: Any | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Record the running Lyapunov estimates as they converge.

    Reinitialises the tangent frame and steps it ``steps`` times, capturing
    the running :meth:`exponents` estimate after every step.  The estimates
    settle as the time-average accumulates — the curve a user inspects to
    judge whether a Lyapunov run has converged.

    Parameters
    ----------
    steps : int, optional
        Number of tangent steps to record.  Default ``2000``.
    n_or_dt : float, optional
        Per-step increment (iterations for a map, ``dt`` for an ODE).
        ``None`` uses the family default.
    ic : array-like, optional
        Initial condition; ``None`` resolves the system's default.

    Returns
    -------
    (times, estimates)
        ``times`` shape ``(steps,)``; ``estimates`` shape ``(steps, k)`` —
        the running estimate of each of the ``k`` exponents after each step.
    """
    if steps <= 0:
        raise ValueError(f"steps must be positive, got {steps}")
    # Random-IC retry on divergence, mirroring the map Lyapunov path
    # (:meth:`_map_spectrum`): when called without an explicit ``ic`` the
    # system resolves a random one, which for a small-basin map (Hénon) can
    # land off-attractor and diverge.  Re-draw a fresh IC and retry; an
    # explicitly supplied ``ic`` that diverges is the caller's, so re-raise.
    max_retries = 10
    for attempt in range(max_retries):
        use_ic = ic if attempt == 0 else None
        times = np.empty(steps)
        estimates = np.empty((steps, self.k))
        try:
            self.reinit(use_ic)
            for s in range(steps):
                self.step(n_or_dt)
                times[s] = self.time()
                estimates[s] = self.exponents()
        except RuntimeError:
            if ic is not None or attempt == max_retries - 1:
                raise
            # Force a fresh random IC on the next attempt.
            object.__setattr__(self.system, "ic", None)
            continue
        return times, estimates
    # Unreachable: the loop returns on success or re-raises on the last
    # attempt; this satisfies the type checker that all paths are covered.
    raise RuntimeError(  # pragma: no cover
        f"{type(self.system).__name__}: tangent convergence diverged from every IC."
    )

to_plot_spec ¶

to_plot_spec(
    kind: str | None = None,
    *,
    steps: int = 2000,
    n_or_dt: float | None = None,
    ic: Any | None = None,
) -> PlotSpec

Describe the Lyapunov-estimate convergence as a :class:PlotSpec.

Builds a :data:~tsdynamics.viz.spec.PlotKind.DIAGNOSTIC_CURVE of each exponent's running estimate against time — a labelled family of lines (one LINE layer per exponent, legended), the standard read-out for "has the Lyapunov spectrum settled?". The estimates are collected via :meth:convergence.

The :mod:tsdynamics.viz.spec import is lazy, so building a spec never pulls in a plotting backend.

PARAMETER	DESCRIPTION
`kind`	Override the semantic kind. `None` (the default) uses `DIAGNOSTIC_CURVE`. TYPE: `str` DEFAULT: `None`
`steps`	Number of tangent steps to record. Default `2000`. TYPE: `int` DEFAULT: `2000`
`n_or_dt`	Per-step increment (iterations for a map, `dt` for an ODE). TYPE: `float` DEFAULT: `None`
`ic`	Initial condition; `None` resolves the system's default. TYPE: `array - like` DEFAULT: `None`

RETURNS	DESCRIPTION
`PlotSpec`

Source code in src/tsdynamics/derived/tangent.py

def to_plot_spec(
    self,
    kind: str | None = None,
    *,
    steps: int = 2000,
    n_or_dt: float | None = None,
    ic: Any | None = None,
) -> PlotSpec:
    """Describe the Lyapunov-estimate convergence as a :class:`PlotSpec`.

    Builds a :data:`~tsdynamics.viz.spec.PlotKind.DIAGNOSTIC_CURVE` of each
    exponent's running estimate against time — a labelled family of lines
    (one ``LINE`` layer per exponent, legended), the standard read-out for
    "has the Lyapunov spectrum settled?".  The estimates are collected via
    :meth:`convergence`.

    The :mod:`tsdynamics.viz.spec` import is lazy, so building a spec never
    pulls in a plotting backend.

    Parameters
    ----------
    kind : str, optional
        Override the semantic kind.  ``None`` (the default) uses
        ``DIAGNOSTIC_CURVE``.
    steps : int, optional
        Number of tangent steps to record.  Default ``2000``.
    n_or_dt : float, optional
        Per-step increment (iterations for a map, ``dt`` for an ODE).
    ic : array-like, optional
        Initial condition; ``None`` resolves the system's default.

    Returns
    -------
    PlotSpec
    """
    from tsdynamics.viz.spec import Axis, Layer, Legend, PlotKind, PlotSpec

    times, estimates = self.convergence(steps, n_or_dt, ic=ic)
    spec_kind = PlotKind(kind) if kind is not None else PlotKind.DIAGNOSTIC_CURVE
    layers = [
        Layer(
            PlotKind.LINE,
            {"x": times, "y": estimates[:, i]},
            label=f"$\\lambda_{{{i + 1}}}$",
        )
        for i in range(self.k)
    ]
    return PlotSpec(
        kind=spec_kind,
        ndim=2,
        title=f"Lyapunov convergence — {type(self.system).__name__}",
        x=Axis(label="iteration" if self.is_discrete else "time"),
        y=Axis(label="Lyapunov estimate"),
        layers=layers,
        legend=Legend(),
    )

lyapunov_spectrum ¶

lyapunov_spectrum(**kwargs: Any) -> ndarray

Estimate the Lyapunov spectrum — the unified engine for every family.

Map and ODE families both delegate their lyapunov_spectrum here, so the QR/variational machinery lives in exactly one place. Mode-specific keywords:

maps: steps (default 5000), ic, reortho_interval (1).
ODEs: final_time (200.0), dt (0.1), ic, burn_in (50.0), method, rtol (1e-6), atol (1e-9), and any extra integrator keywords.

The estimate is recorded in self.meta['lyapunov_spectrum'] (the inner system's :class:~tsdynamics.families.base.MetaStore).

RETURNS	DESCRIPTION
`(ndarray, shape(k))`	Lyapunov exponents, largest first (QR order).

Source code in src/tsdynamics/derived/tangent.py

def lyapunov_spectrum(self, **kwargs: Any) -> np.ndarray:
    """
    Estimate the Lyapunov spectrum — the unified engine for every family.

    Map and ODE families both delegate their ``lyapunov_spectrum`` here, so
    the QR/variational machinery lives in exactly one place.  Mode-specific
    keywords:

    - **maps**: ``steps`` (default 5000), ``ic``, ``reortho_interval`` (1).
    - **ODEs**: ``final_time`` (200.0), ``dt`` (0.1), ``ic``, ``burn_in``
      (50.0), ``method``, ``rtol`` (1e-6), ``atol`` (1e-9), and any extra
      integrator keywords.

    The estimate is recorded in ``self.meta['lyapunov_spectrum']`` (the inner
    system's :class:`~tsdynamics.families.base.MetaStore`).

    Returns
    -------
    ndarray, shape (k,)
        Lyapunov exponents, largest first (QR order).
    """
    if self._mode == "map":
        return self._lyapunov_spectrum_map(**kwargs)
    return self._lyapunov_spectrum_ode(**kwargs)

EnsembleSystem ¶

EnsembleSystem(system: Any, states: Any)

Many copies of one system, advanced synchronously from different states.

Used for two-trajectory Lyapunov estimates, basin sampling, and ensemble statistics. Members are independent copies — parameters are shared at construction, states are per-member.

PARAMETER	DESCRIPTION
`system`	The template system (copied per member; the original is untouched). TYPE: `System`
`states`	One initial state per member. TYPE: `(array - like, shape(m, dim))`

Examples:

>>> ens = EnsembleSystem(Lorenz(), [[1, 1, 1], [1.001, 1, 1]])
>>> ens.step(0.01)
array([[...], [...]])

Source code in src/tsdynamics/derived/ensemble.py

def __init__(self, system: Any, states: Any) -> None:
    states_arr = np.atleast_2d(np.asarray(states, dtype=float))
    if states_arr.shape[1] != system.dim:
        raise ValueError(f"states must have shape (m, {system.dim}), got {states_arr.shape}")
    self.template = system
    self.members = []
    for s in states_arr:
        member = system.copy()
        member.reinit(s)
        self.members.append(member)

size `property` ¶

size: int

Number of ensemble members.

dim `property` ¶

dim: int

State-space dimension of each member.

is_discrete `property` ¶

is_discrete: bool

Match the template system's time semantics.

step ¶

step(n_or_dt: float | int | None = None) -> ndarray

Advance every member synchronously and return the stacked states.

PARAMETER	DESCRIPTION
`n_or_dt`	The per-member increment (a `dt` for a flow, an iteration count for a map). `None` uses each member's default. TYPE: `float or int` DEFAULT: `None`

RETURNS	DESCRIPTION
`ndarray`	The new states, shape `(size, dim)` — one row per member.

Source code in src/tsdynamics/derived/ensemble.py

def step(self, n_or_dt: float | int | None = None) -> np.ndarray:
    """Advance every member synchronously and return the stacked states.

    Parameters
    ----------
    n_or_dt : float or int, optional
        The per-member increment (a ``dt`` for a flow, an iteration count for
        a map).  ``None`` uses each member's default.

    Returns
    -------
    numpy.ndarray
        The new states, shape ``(size, dim)`` — one row per member.
    """
    return np.array([m.step(n_or_dt) for m in self.members])

states ¶

states() -> ndarray

Return the current member states, shape (size, dim).

Source code in src/tsdynamics/derived/ensemble.py

def states(self) -> np.ndarray:
    """Return the current member states, shape ``(size, dim)``."""
    return np.array([m.state() for m in self.members])

set_states ¶

set_states(states: Any) -> None

Overwrite every member's state.

PARAMETER	DESCRIPTION
`states`	One new state per member, in member order. TYPE: `(array - like, shape(size, dim))`

RAISES	DESCRIPTION
`ValueError`	If `states` is not exactly `(size, dim)`.

Source code in src/tsdynamics/derived/ensemble.py

def set_states(self, states: Any) -> None:
    """Overwrite every member's state.

    Parameters
    ----------
    states : array-like, shape (size, dim)
        One new state per member, in member order.

    Raises
    ------
    ValueError
        If ``states`` is not exactly ``(size, dim)``.
    """
    states_arr = np.atleast_2d(np.asarray(states, dtype=float))
    if states_arr.shape != (self.size, self.dim):
        raise ValueError(f"expected shape {(self.size, self.dim)}, got {states_arr.shape}")
    for member, s in zip(self.members, states_arr, strict=True):
        member.set_state(s)

time ¶

time() -> float

Return the common member time.

Source code in src/tsdynamics/derived/ensemble.py

def time(self) -> float:
    """Return the common member time."""
    return self.members[0].time() if self.members else 0.0

collect ¶

collect(
    steps: int, n_or_dt: float | int | None = None
) -> tuple[ndarray, ndarray]

Step every member steps times and stack the sampled states.

Advances the whole ensemble synchronously, recording each member's state after every step. This is the trajectory collector the static fan chart (:meth:to_plot_spec) summarises into a median line + percentile band.

PARAMETER	DESCRIPTION
`steps`	Number of samples to collect (one per step). TYPE: `int`
`n_or_dt`	The per-step increment forwarded to each member's `step` (a `dt` for a flow, an iteration count for a map). `None` uses the member default. TYPE: `float or int` DEFAULT: `None`

RETURNS	DESCRIPTION
`(times, states)`	`times` shape `(steps,)` (the common member time after each step); `states` shape `(steps, size, dim)` — sample, member, component.

Source code in src/tsdynamics/derived/ensemble.py

def collect(
    self, steps: int, n_or_dt: float | int | None = None
) -> tuple[np.ndarray, np.ndarray]:
    """Step every member ``steps`` times and stack the sampled states.

    Advances the whole ensemble synchronously, recording each member's state
    after every step.  This is the trajectory collector the static fan chart
    (:meth:`to_plot_spec`) summarises into a median line + percentile band.

    Parameters
    ----------
    steps : int
        Number of samples to collect (one per step).
    n_or_dt : float or int, optional
        The per-step increment forwarded to each member's ``step`` (a ``dt``
        for a flow, an iteration count for a map).  ``None`` uses the member
        default.

    Returns
    -------
    (times, states)
        ``times`` shape ``(steps,)`` (the common member time after each step);
        ``states`` shape ``(steps, size, dim)`` — sample, member, component.
    """
    if steps <= 0:
        raise ValueError(f"steps must be positive, got {steps}")
    times = np.empty(steps)
    states = np.empty((steps, self.size, self.dim))
    for k in range(steps):
        states[k] = self.step(n_or_dt)
        times[k] = self.time()
    return times, states

to_plot_spec ¶

to_plot_spec(
    kind: str | None = None,
    *,
    steps: int = 200,
    component: int = 0,
    band: float = 90.0,
) -> PlotSpec

Describe the ensemble as a static fan chart (median + percentile band).

Collects the ensemble's evolution of one component and summarises the spread across members at each time as a shaded percentile band (an AREA layer carrying "lo" / "hi" band edges, with lo <= hi) under the across-member median line — the standard, animation-free way to read an ensemble's dispersion. This is not an animation: it is one :data:~tsdynamics.viz.spec.PlotKind.ENSEMBLE_FAN static spec.

The :mod:tsdynamics.viz.spec import is lazy, so building a spec never pulls in a plotting backend.

PARAMETER	DESCRIPTION
`kind`	Override the semantic kind. `None` (the default) uses `ENSEMBLE_FAN`. TYPE: `str` DEFAULT: `None`
`steps`	Number of samples to collect across the ensemble. Default `200`. TYPE: `int` DEFAULT: `200`
`component`	Which state component to chart. Default `0`. TYPE: `int` DEFAULT: `0`
`band`	Central percentile mass to shade (`90` → the 5th–95th percentile band). Default `90.0`; clamped to `(0, 100]`. TYPE: `float` DEFAULT: `90.0`

RETURNS	DESCRIPTION
`PlotSpec`

Source code in src/tsdynamics/derived/ensemble.py

def to_plot_spec(
    self, kind: str | None = None, *, steps: int = 200, component: int = 0, band: float = 90.0
) -> PlotSpec:
    """Describe the ensemble as a **static fan chart** (median + percentile band).

    Collects the ensemble's evolution of one component and summarises the
    spread across members at each time as a shaded percentile band (an
    ``AREA`` layer carrying ``"lo"`` / ``"hi"`` band edges, with ``lo <= hi``)
    under the across-member **median** line — the standard, animation-free way
    to read an ensemble's dispersion.  This is **not** an animation: it is one
    :data:`~tsdynamics.viz.spec.PlotKind.ENSEMBLE_FAN` static spec.

    The :mod:`tsdynamics.viz.spec` import is lazy, so building a spec never
    pulls in a plotting backend.

    Parameters
    ----------
    kind : str, optional
        Override the semantic kind.  ``None`` (the default) uses
        ``ENSEMBLE_FAN``.
    steps : int, optional
        Number of samples to collect across the ensemble.  Default ``200``.
    component : int, optional
        Which state component to chart.  Default ``0``.
    band : float, optional
        Central percentile mass to shade (``90`` → the 5th–95th percentile
        band).  Default ``90.0``; clamped to ``(0, 100]``.

    Returns
    -------
    PlotSpec
    """
    from tsdynamics.viz.spec import Axis, Layer, PlotKind, PlotSpec

    if not 0 <= component < self.dim:
        raise ValueError(f"component must be in [0, {self.dim}), got {component}")
    if not 0.0 < band <= 100.0:
        raise ValueError(f"band must be in (0, 100], got {band}")

    times, states = self.collect(steps)
    comp = states[:, :, component]  # (steps, size)
    lo_pct = (100.0 - band) / 2.0
    hi_pct = 100.0 - lo_pct
    lo = np.percentile(comp, lo_pct, axis=1)
    hi = np.percentile(comp, hi_pct, axis=1)
    median = np.median(comp, axis=1)
    # The band edges are percentiles of the same sample, so lo <= hi holds by
    # construction; enforce it defensively against any float ordering quirk.
    lo = np.minimum(lo, hi)

    names = getattr(type(self.template), "variables", None)
    ylabel = names[component] if names is not None else f"y{component}"
    spec_kind = PlotKind(kind) if kind is not None else PlotKind.ENSEMBLE_FAN
    return PlotSpec(
        kind=spec_kind,
        ndim=2,
        title=f"Ensemble fan — {type(self.template).__name__} (n={self.size})",
        x=Axis(label="iteration" if self.is_discrete else "time"),
        y=Axis(label=ylabel),
        layers=[
            Layer(
                PlotKind.AREA,
                {"x": times, "y": median, "lo": lo, "hi": hi},
                label=f"{int(round(band))}% band",
                style={"alpha": 0.3},
            ),
            Layer(
                PlotKind.LINE,
                {"x": times, "y": median},
                label="median",
            ),
        ],
    )

ProjectedSystem ¶

ProjectedSystem(
    system: Any,
    components: Any,
    *,
    complete: Callable[[ndarray], Any] | None = None,
)

Bases: DerivedSystem

View a system through a subset of its components.

The full system is stepped underneath; only state()/step() outputs are projected. set_state needs the inverse direction and therefore requires a complete callable mapping a projected state back to a full state.

PARAMETER	DESCRIPTION
`system`	The full system. TYPE: `System`
`components`	Component indices (or names, when the system declares `variables`). TYPE: `sequence of int or str`
`complete`	`complete(u_projected) -> u_full` used by `set_state`/`reinit` when given projected-dimensional inputs. TYPE: `callable` DEFAULT: `None`

Examples:

>>> proj = ProjectedSystem(Lorenz(), ["x", "z"])
>>> proj.step(0.01).shape
(2,)

Source code in src/tsdynamics/derived/projected.py

def __init__(
    self,
    system: Any,
    components: Any,
    *,
    complete: Callable[[np.ndarray], Any] | None = None,
) -> None:
    super().__init__(system)
    names = getattr(type(system), "variables", None)
    idx = []
    for c in components:
        if isinstance(c, str):
            if names is None:
                raise ValueError(
                    f"{type(system).__name__} declares no `variables`; "
                    f"use integer component indices."
                )
            idx.append(names.index(c))
        else:
            idx.append(int(c))
    if not idx:
        raise ValueError("components must be non-empty")
    self.components = tuple(idx)
    self.complete = complete

dim `property` ¶

dim: int

Dimension of the projected view.

variables `property` ¶

variables: tuple[str, ...] | None

Component names of the projected view (the inner names, subset).

Overrides :class:DerivedSystem's pass-through (which would return the inner system's full names and mislabel the projected columns). Returns None when the inner system declares no variables.

step ¶

step(n_or_dt: float | int | None = None) -> ndarray

Advance the full system; return the projected new state.

Source code in src/tsdynamics/derived/projected.py

def step(self, n_or_dt: float | int | None = None) -> np.ndarray:
    """Advance the full system; return the projected new state."""
    return cast(np.ndarray, self.system.step(n_or_dt)[list(self.components)])

state ¶

state() -> ndarray

Return the projected current state.

Source code in src/tsdynamics/derived/projected.py

def state(self) -> np.ndarray:
    """Return the projected current state."""
    return cast(np.ndarray, self.system.state()[list(self.components)])

set_state ¶

set_state(u: Any) -> None

Overwrite the state (projected inputs need a complete callable).

Source code in src/tsdynamics/derived/projected.py

def set_state(self, u: Any) -> None:
    """Overwrite the state (projected inputs need a ``complete`` callable)."""
    u_arr = np.asarray(u, dtype=float)
    self.system.set_state(self._to_full_state(u_arr, verb="set_state"))

reinit ¶

reinit(u: Any | None = None, **kwargs: Any) -> None

Restart the full system (projected inputs need a complete callable).

Source code in src/tsdynamics/derived/projected.py

def reinit(self, u: Any | None = None, **kwargs: Any) -> None:
    """Restart the full system (projected inputs need a ``complete`` callable)."""
    if u is not None:
        u = self._to_full_state(np.asarray(u, dtype=float), verb="reinit")
    self.system.reinit(u, **kwargs)

trajectory ¶

trajectory(*args: Any, **kwargs: Any) -> Trajectory

Full-system trajectory with projected columns.

Source code in src/tsdynamics/derived/projected.py

def trajectory(self, *args: Any, **kwargs: Any) -> Trajectory:
    """Full-system trajectory with projected columns."""
    traj = self.system.trajectory(*args, **kwargs)
    meta = {**traj.meta, "projected": self.components}
    # Back-reference ``self`` (not the inner system): the returned ``y`` holds
    # only the projected columns, and ``self.variables`` names exactly those —
    # so ``traj["x"]`` resolves to the right column and an unknown name raises
    # KeyError rather than silently mislabelling or IndexError-ing.
    return Trajectory(traj.t, traj.y[:, list(self.components)], self, meta=meta)

WrappedSystem ¶

WrappedSystem(
    step_fn: Callable[[ndarray, float], Any],
    *,
    dim: int,
    is_discrete: bool = True,
    initial: Any | None = None,
    default_dt: float = 1.0,
    variables: tuple[str, ...] | None = None,
)

Wrap any external stepping rule as a first-class :class:System.

Give it a step_fn(state, n_or_dt) -> new_state and a dimension, and the whole analysis toolkit (orbit diagrams, Lyapunov-from-rescaling, Poincaré sections, ensembles, basins) applies to your own simulation code — a foreign ODE solver, an agent-based model, a hardware-in-the-loop rig, anything that advances a state vector.

PARAMETER	DESCRIPTION
`step_fn`	`step_fn(state, n_or_dt) -> new_state`. `new_state` is array-like of length `dim`. `n_or_dt` is the iteration count (discrete) or the time increment (continuous); pass it through to your stepper. TYPE: `callable`
`dim`	State-space dimension. TYPE: `int`
`is_discrete`	Whether `n_or_dt` counts iterations (True) or measures time (False). TYPE: `bool` DEFAULT: `True`
`initial`	Default initial state used when `reinit` gets no explicit state. TYPE: `array - like` DEFAULT: `None`
`default_dt`	Step taken by `step()` when called with no argument. TYPE: `float` DEFAULT: `1.0`
`variables`	Component names, enabling `traj["x"]` on produced trajectories. TYPE: `tuple of str` DEFAULT: `None`

Examples:

>>> # a plain logistic map written by hand
>>> import numpy as np
>>> def step(u, n):
...     x = u[0]
...     for _ in range(int(n)):
...         x = 3.9 * x * (1 - x)
...     return [x]
>>> sysm = WrappedSystem(step, dim=1, is_discrete=True, initial=[0.5])
>>> traj = sysm.trajectory(500)
>>> import tsdynamics as ts
>>> ts.max_lyapunov(sysm, ic=[0.3]) > 0          # chaotic
True

Source code in src/tsdynamics/families/wrapped.py

def __init__(
    self,
    step_fn: Callable[[np.ndarray, float], Any],
    *,
    dim: int,
    is_discrete: bool = True,
    initial: Any | None = None,
    default_dt: float = 1.0,
    variables: tuple[str, ...] | None = None,
) -> None:
    self._step_fn = step_fn
    self.dim = int(dim)
    self._is_discrete = bool(is_discrete)
    self._initial = None if initial is None else np.asarray(initial, dtype=float).reshape(dim)
    self._default_dt = float(default_dt)
    self.variables = variables

    self._state: np.ndarray | None = None
    self._t: float = 0.0

is_discrete `property` ¶

is_discrete: bool

Whether stepping counts iterations rather than time.

reinit ¶

reinit(
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
) -> None

Restart from state u (falls back to initial, then zeros).

Source code in src/tsdynamics/families/wrapped.py

def reinit(
    self,
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
) -> None:
    """Restart from state ``u`` (falls back to ``initial``, then zeros)."""
    if u is not None:
        self._state = np.asarray(u, dtype=float).reshape(self.dim)
    elif self._initial is not None:
        self._state = self._initial.copy()
    else:
        self._state = np.zeros(self.dim)
    self._t = float(t) if t is not None else 0.0

step ¶

step(n_or_dt: float | None = None) -> ndarray

Advance by n_or_dt (default default_dt) and return the new state.

Source code in src/tsdynamics/families/wrapped.py

def step(self, n_or_dt: float | None = None) -> np.ndarray:
    """Advance by ``n_or_dt`` (default ``default_dt``) and return the new state."""
    if self._state is None:
        self.reinit()
    assert self._state is not None
    amount = self._default_dt if n_or_dt is None else n_or_dt
    new = np.asarray(self._step_fn(self._state, amount), dtype=np.float64).reshape(self.dim)
    if not np.all(np.isfinite(new)):
        raise ConvergenceError("WrappedSystem.step produced non-finite state")
    self._state = new
    self._t += amount
    return self._state.copy()

state ¶

state() -> ndarray

Return a copy of the current state.

Source code in src/tsdynamics/families/wrapped.py

def state(self) -> np.ndarray:
    """Return a copy of the current state."""
    if self._state is None:
        self.reinit()
    assert self._state is not None
    return self._state.copy()

set_state ¶

set_state(u: Any) -> None

Overwrite the current state.

Source code in src/tsdynamics/families/wrapped.py

def set_state(self, u: Any) -> None:
    """Overwrite the current state."""
    self._state = np.asarray(u, dtype=float).reshape(self.dim)

time ¶

time() -> float

Return the current time (continuous) or iteration count (discrete).

Source code in src/tsdynamics/families/wrapped.py

def time(self) -> float:
    """Return the current time (continuous) or iteration count (discrete)."""
    return self._t

copy ¶

copy() -> WrappedSystem

Return a fresh wrapper sharing the same step rule (independent state).

Source code in src/tsdynamics/families/wrapped.py

def copy(self) -> WrappedSystem:
    """Return a fresh wrapper sharing the same step rule (independent state)."""
    return WrappedSystem(
        self._step_fn,
        dim=self.dim,
        is_discrete=self._is_discrete,
        initial=self._initial,
        default_dt=self._default_dt,
        variables=self.variables,
    )

trajectory ¶

trajectory(
    n: int | None = None,
    *,
    transient: int = 0,
    ic: Any | None = None,
    final_time: float | None = None,
    dt: float | None = None,
) -> Trajectory

Step the wrapper repeatedly (after transient) and collect a Trajectory.

The wrapper accepts both the map-style positional sample count n and — for a continuous wrapper — the family-uniform final_time / dt spelling, so a generic protocol caller (which calls trajectory(final_time=…, dt=…)) works without raising TypeError.

Exactly one count source must be supplied: either n (the number of samples to collect), or final_time. When final_time is given the sample count is round(final_time / dt) and the per-step increment is dt, which defaults to default_dt for both continuous and discrete wrappers. (A discrete wrapper that wants the "one time unit = one iteration" convention should construct itself with default_dt=1.0, the constructor default.)

PARAMETER	DESCRIPTION
`n`	Number of samples to collect. Mutually exclusive with `final_time`. TYPE: `int` DEFAULT: `None`
`transient`	Number of leading samples to discard (steps taken but not recorded). TYPE: `int` DEFAULT: `0`
`ic`	Initial state for the run (falls back to `initial`, then zeros). TYPE: `array - like` DEFAULT: `None`
`final_time`	Integration horizon, an alternative to `n` for continuous wrappers (and accepted for discrete ones). The sample count is `round(final_time / dt)`. TYPE: `float` DEFAULT: `None`
`dt`	Per-step increment used with `final_time` (and as the `step` argument). Defaults to `default_dt`. TYPE: `float` DEFAULT: `None`

RETURNS	DESCRIPTION
`Trajectory`	`n` recorded samples with their times.

RAISES	DESCRIPTION
`InvalidInputError`	If neither `n` nor `final_time` is given, if both are, or if the resolved sample count is not positive.

Examples:

>>> import numpy as np
>>> flow = lambda u, dt: [u[0] * np.exp(0.5 * dt)]
>>> w = WrappedSystem(flow, dim=1, is_discrete=False, default_dt=0.1)
>>> traj = w.trajectory(final_time=1.0, dt=0.1)   # family-uniform spelling
>>> traj.y.shape
(10, 1)

Source code in src/tsdynamics/families/wrapped.py

def trajectory(
    self,
    n: int | None = None,
    *,
    transient: int = 0,
    ic: Any | None = None,
    final_time: float | None = None,
    dt: float | None = None,
) -> Trajectory:
    """Step the wrapper repeatedly (after ``transient``) and collect a Trajectory.

    The wrapper accepts **both** the map-style positional sample count ``n``
    and — for a *continuous* wrapper — the family-uniform ``final_time`` /
    ``dt`` spelling, so a generic protocol caller (which calls
    ``trajectory(final_time=…, dt=…)``) works without raising ``TypeError``.

    Exactly one count source must be supplied: either ``n`` (the number of
    samples to collect), or ``final_time``.  When ``final_time`` is given the
    sample count is ``round(final_time / dt)`` and the per-step increment is
    ``dt``, which defaults to ``default_dt`` for **both** continuous and
    discrete wrappers.  (A discrete wrapper that wants the "one time unit =
    one iteration" convention should construct itself with ``default_dt=1.0``,
    the constructor default.)

    Parameters
    ----------
    n : int, optional
        Number of samples to collect.  Mutually exclusive with ``final_time``.
    transient : int, default 0
        Number of leading samples to discard (steps taken but not recorded).
    ic : array-like, optional
        Initial state for the run (falls back to ``initial``, then zeros).
    final_time : float, optional
        Integration horizon, an alternative to ``n`` for continuous wrappers
        (and accepted for discrete ones).  The sample count is
        ``round(final_time / dt)``.
    dt : float, optional
        Per-step increment used with ``final_time`` (and as the ``step``
        argument).  Defaults to ``default_dt``.

    Returns
    -------
    Trajectory
        ``n`` recorded samples with their times.

    Raises
    ------
    InvalidInputError
        If neither ``n`` nor ``final_time`` is given, if both are, or if the
        resolved sample count is not positive.

    Examples
    --------
    >>> import numpy as np
    >>> flow = lambda u, dt: [u[0] * np.exp(0.5 * dt)]
    >>> w = WrappedSystem(flow, dim=1, is_discrete=False, default_dt=0.1)
    >>> traj = w.trajectory(final_time=1.0, dt=0.1)   # family-uniform spelling
    >>> traj.y.shape
    (10, 1)
    """
    step_dt = self._default_dt if dt is None else float(dt)
    if final_time is not None:
        if n is not None:
            raise InvalidInputError(
                "WrappedSystem.trajectory: pass either n or final_time, not both."
            )
        n = int(round(float(final_time) / step_dt))
    elif n is None:
        raise InvalidInputError(
            "WrappedSystem.trajectory requires n (sample count) or final_time."
        )
    n = int(n)
    if n <= 0:
        raise InvalidInputError(f"WrappedSystem.trajectory needs a positive count, got {n}.")

    self.reinit(ic)
    for _ in range(transient):
        self.step(step_dt)
    ts = np.empty(n)
    ys = np.empty((n, self.dim))
    for i in range(n):
        ys[i] = self.step(step_dt)
        ts[i] = self._t
    meta = {"system": "WrappedSystem", "is_discrete": self._is_discrete}
    return Trajectory(t=ts, y=ys, system=self, meta=meta)

DerivedSystem ¶

DerivedSystem(system: Any)

Base for wrappers that present an existing system through a new lens.

A derived system implements the :class:~tsdynamics.families.System protocol by delegating to a wrapped system, transforming what "one step" or "the state" means (Poincaré crossings, stroboscopic samples, projections...).

Parameters and metadata are forwarded to the wrapped system, and with_params re-parametrizes the inner system and rebuilds the wrapper, so parameter sweeps compose: an orbit diagram over a PoincareMap is a bifurcation diagram of the underlying flow.

Source code in src/tsdynamics/derived/_base.py

def __init__(self, system: Any) -> None:
    self.system = system

with_params ¶

with_params(**overrides: Any) -> DerivedSystem

Return a new wrapper of the same kind around a re-parametrized copy.

Source code in src/tsdynamics/derived/_base.py

def with_params(self, **overrides: Any) -> DerivedSystem:
    """Return a new wrapper of the same kind around a re-parametrized copy."""
    return self._rebuild(self.system.with_params(**overrides))

copy ¶

copy() -> DerivedSystem

Return a new wrapper of the same kind around a copy of the inner system.

Source code in src/tsdynamics/derived/_base.py

def copy(self) -> DerivedSystem:
    """Return a new wrapper of the same kind around a copy of the inner system."""
    return self._rebuild(self.system.copy())

trajectory ¶

trajectory(*args: Any, **kwargs: Any) -> Trajectory

Produce the wrapper's trajectory — subclasses implement the lens-specific collection.

Source code in src/tsdynamics/derived/_base.py

def trajectory(self, *args: Any, **kwargs: Any) -> Trajectory:
    """Produce the wrapper's trajectory — subclasses implement the lens-specific collection."""
    raise NotImplementedError

run ¶

run(*args: Any, **kwargs: Any) -> Trajectory

Produce the wrapper's trajectory — the alias of :meth:trajectory.

run is the library's canonical trajectory-producer verb (a flow's Lorenz().run(...), a map's Henon().run(...)), so a fluent derived view reads left-to-right with the same verb at the end::

section = Rossler().poincare(section="y", at=0.0).run(steps=500)

It forwards verbatim to this wrapper's :meth:trajectory, so the two are byte-identical and every wrapper-specific keyword (transient, ...) is honoured. trajectory remains the member the structural System protocol requires; run is the discoverable spelling.

Source code in src/tsdynamics/derived/_base.py

def run(self, *args: Any, **kwargs: Any) -> Trajectory:
    """Produce the wrapper's trajectory — the alias of :meth:`trajectory`.

    ``run`` is the library's canonical trajectory-producer verb (a flow's
    ``Lorenz().run(...)``, a map's ``Henon().run(...)``), so a fluent
    derived view reads left-to-right with the same verb at the end::

        section = Rossler().poincare(section="y", at=0.0).run(steps=500)

    It forwards verbatim to this wrapper's :meth:`trajectory`, so the two are
    byte-identical and every wrapper-specific keyword (``transient``, ...) is
    honoured.  ``trajectory`` remains the member the structural ``System``
    protocol requires; ``run`` is the discoverable spelling.
    """
    return self.trajectory(*args, **kwargs)

to_plot_spec ¶

to_plot_spec(kind: str | None = None) -> PlotSpec

Describe this derived view as a backend-agnostic :class:PlotSpec.

The default delegates to the wrapper's own :meth:trajectory: it collects the lens-specific trajectory (Poincaré crossings, projected columns, ...) and forwards to that trajectory's :meth:~tsdynamics.data.Trajectory.to_plot_spec. Subclasses whose natural picture is not a single trajectory line — a stroboscopic scatter, an ensemble fan, a Lyapunov convergence curve — override this with their own spec builder.

The :mod:tsdynamics.viz.spec import stays inside the trajectory method (lazy), so building a spec never pulls in a plotting backend.

PARAMETER	DESCRIPTION
`kind`	Override the auto-dispatched semantic kind with any member of the closed :class:`~tsdynamics.viz.spec.PlotKind` vocabulary. `None` (the default) lets the underlying trajectory auto-dispatch. TYPE: `str` DEFAULT: `None`

RETURNS	DESCRIPTION
`PlotSpec`

Source code in src/tsdynamics/derived/_base.py

def to_plot_spec(self, kind: str | None = None) -> PlotSpec:
    """Describe this derived view as a backend-agnostic :class:`PlotSpec`.

    The default delegates to the wrapper's own :meth:`trajectory`: it
    collects the lens-specific trajectory (Poincaré crossings, projected
    columns, ...) and forwards to that trajectory's
    :meth:`~tsdynamics.data.Trajectory.to_plot_spec`.  Subclasses whose
    natural picture is *not* a single trajectory line — a stroboscopic
    scatter, an ensemble fan, a Lyapunov convergence curve — override this
    with their own spec builder.

    The :mod:`tsdynamics.viz.spec` import stays inside the trajectory method
    (lazy), so building a spec never pulls in a plotting backend.

    Parameters
    ----------
    kind : str, optional
        Override the auto-dispatched semantic kind with any member of the
        closed :class:`~tsdynamics.viz.spec.PlotKind` vocabulary.  ``None``
        (the default) lets the underlying trajectory auto-dispatch.

    Returns
    -------
    PlotSpec
    """
    return self.trajectory().to_plot_spec(kind=kind)

Derived systems¶

PoincareMap ¶

is_discrete property ¶

crossing_count property ¶

step ¶

state ¶

set_state ¶

time ¶

as_events ¶

reinit ¶

trajectory ¶

StroboscopicMap ¶

is_discrete property ¶

step ¶

time ¶

trajectory ¶

to_plot_spec ¶

TangentSystem ¶

is_discrete property ¶

reinit ¶

step ¶

state ¶

set_state ¶

time ¶

deviations ¶

growths ¶

exponents ¶

convergence ¶

to_plot_spec ¶

lyapunov_spectrum ¶

EnsembleSystem ¶

size property ¶

dim property ¶

is_discrete property ¶

step ¶

states ¶

set_states ¶

time ¶

collect ¶

to_plot_spec ¶

ProjectedSystem ¶

dim property ¶

variables property ¶

step ¶

state ¶

set_state ¶

reinit ¶

trajectory ¶

WrappedSystem ¶

is_discrete property ¶

reinit ¶

step ¶

state ¶

set_state ¶

time ¶

copy ¶

trajectory ¶

DerivedSystem ¶

with_params ¶

copy ¶

trajectory ¶

run ¶

to_plot_spec ¶

is_discrete `property` ¶

crossing_count `property` ¶

is_discrete `property` ¶

is_discrete `property` ¶

size `property` ¶

dim `property` ¶

is_discrete `property` ¶

dim `property` ¶

variables `property` ¶

is_discrete `property` ¶