Reference

Base classes

The shared machinery (ParamSet, MetaStore, Trajectory, SystemBase), the four family bases users subclass — ContinuousSystem (ODEs), DelaySystem (DDEs), DiscreteMap (maps) and StochasticSystem (diagonal-Itô SDEs) — and the System protocol the analysis toolkit is written against.

ParamSet

ParamSet(data: dict[str, Any])

Bases: MutableMapping[str, Any]

Ordered, fixed-key parameter container.

Keys are frozen at construction time — you can change values but not add or remove keys. Supports both dict-style (p["sigma"]) and attribute-style (p.sigma) read/write.

PARAMETER	DESCRIPTION
data	Initial key→value mapping. All future writes must use existing keys. TYPE: dict

RAISES	DESCRIPTION
AttributeError	On attribute-style read or write of an undeclared key (p.unknown / p.unknown = ...).
KeyError	On item-style write of an undeclared key (p["unknown"] = ...).
InvalidInputError	On any attempt to delete a key (the key set is frozen).

Examples:

>>> p = ParamSet({"sigma": 10.0, "rho": 28.0})
>>> p.sigma
10.0
>>> p.sigma = 15.0
>>> p["sigma"]
15.0
>>> p.unknown = 5.0            # raises AttributeError

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

def __init__(self, data: dict[str, Any]) -> None:
    object.__setattr__(self, "_data", dict(data))

as_tuple

as_tuple() -> tuple[Any, ...]

Return parameter values as a tuple (insertion order).

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

def as_tuple(self) -> tuple[Any, ...]:
    """Return parameter values as a tuple (insertion order)."""
    return tuple(self._data.values())

as_dict

as_dict() -> dict[str, Any]

Return a shallow copy as a plain dict.

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

def as_dict(self) -> dict[str, Any]:
    """Return a shallow copy as a plain dict."""
    return dict(self._data)

param_hash

param_hash() -> int

Return a process-stable 64-bit integer hash of the current parameter values.

Uses MD5 over a JSON-serialised representation so the result is reproducible across Python process restarts (unlike hash()).

The hash backs cache keys for per-system lowering / lambdify caches. At 64 bits the birthday-paradox collision probability reaches 50 % only around 2^32 ≈ 4·10⁹ distinct parameter sets, which is well beyond any realistic parameter sweep. The previous 32-bit width hit the same threshold at only 2^16 ≈ 65 000 sets, which a sufficiently large sweep could plausibly reach — and a collision there would silently return a compiled artifact built for a different parameter set.

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

def param_hash(self) -> int:
    """
    Return a process-stable 64-bit integer hash of the current parameter values.

    Uses MD5 over a JSON-serialised representation so the result is
    reproducible across Python process restarts (unlike ``hash()``).

    The hash backs cache keys for per-system lowering / lambdify caches.
    At 64 bits the birthday-paradox
    collision probability reaches 50 % only around ``2^32 ≈ 4·10⁹``
    distinct parameter sets, which is well beyond any realistic
    parameter sweep.  The previous 32-bit width hit the same threshold
    at only ``2^16 ≈ 65 000`` sets, which a sufficiently large sweep
    could plausibly reach — and a collision there would silently
    return a compiled artifact built for a different parameter set.
    """
    import hashlib
    import json

    s = json.dumps(list(self._data.items()), default=str)
    return int(hashlib.md5(s.encode()).hexdigest()[:16], 16)

MetaStore

MetaStore()

Bases: MutableMapping[str, Any]

Append-with-history metadata store for computed results.

Behaves like a dict for everyday use (meta["lyapunov_spectrum"] reads/writes the latest value), but every write is appended rather than overwritten, so earlier results survive::

sys.meta.record("lyapunov_spectrum", spec, dt=0.1, final_time=200.0)
sys.meta["lyapunov_spectrum"]            # latest value
sys.meta.history("lyapunov_spectrum")    # every record, with context

Equality compares the latest values against a plain dict (or another MetaStore), preserving sys.meta == {} style assertions.

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

def __init__(self) -> None:
    self._records: dict[str, list[dict[str, Any]]] = {}

record

record(key: str, value: Any, **context: Any) -> Any

Append value under key with optional context kwargs.

Each call stores a new record {"value", "context", "timestamp"} — earlier records under the same key are preserved (retrievable with :meth:history), and meta[key] returns the most recent value.

PARAMETER	DESCRIPTION
key	The result name (e.g. "lyapunov_spectrum"). TYPE: str
value	The computed value to store. TYPE: Any
**context	Free-form context recorded alongside the value (e.g. dt, final_time), surfaced by :meth:history. TYPE: Any DEFAULT: {}

RETURNS	DESCRIPTION
Any	value unchanged, so a result can be recorded and returned in one expression (return self.meta.record("k", v)).

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

def record(self, key: str, value: Any, **context: Any) -> Any:
    """Append ``value`` under ``key`` with optional context kwargs.

    Each call stores a new record ``{"value", "context", "timestamp"}`` —
    earlier records under the same key are preserved (retrievable with
    :meth:`history`), and ``meta[key]`` returns the most recent value.

    Parameters
    ----------
    key : str
        The result name (e.g. ``"lyapunov_spectrum"``).
    value : Any
        The computed value to store.
    **context
        Free-form context recorded alongside the value (e.g. ``dt``,
        ``final_time``), surfaced by :meth:`history`.

    Returns
    -------
    Any
        ``value`` unchanged, so a result can be recorded and returned in one
        expression (``return self.meta.record("k", v)``).
    """
    import time

    self._records.setdefault(key, []).append(
        {"value": value, "context": context, "timestamp": time.time()}
    )
    return value

history

history(key: str) -> list[dict[str, Any]]

Return every record for key (oldest first), with context.

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

def history(self, key: str) -> list[dict[str, Any]]:
    """Return every record for ``key`` (oldest first), with context."""
    return list(self._records.get(key, []))

latest

latest() -> dict[str, Any]

Return a plain dict of the latest value per key.

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

def latest(self) -> dict[str, Any]:
    """Return a plain dict of the latest value per key."""
    return {k: recs[-1]["value"] for k, recs in self._records.items()}

Trajectory

Trajectory(
    t: ndarray,
    y: ndarray,
    system: Any,
    meta: dict[str, Any] | None = None,
)

The result of integrating or iterating a dynamical system.

Supports tuple-unpacking for backward compatibility::

t, y = system.integrate(final_time=100)

ATTRIBUTE	DESCRIPTION
t	Time points (or step indices for discrete maps). TYPE: (ndarray, shape(T))
y	State at each time point. TYPE: (ndarray, shape(T, dim))
system	Back-reference to the system that produced this trajectory. TYPE: SystemBase
meta	Provenance: system name, params snapshot, solver, tolerances, ic. TYPE: dict

Examples:

>>> traj = lor.integrate(final_time=100)
>>> traj.dim
3
>>> traj["x"]            # named component (via the class's ``variables``)
array([...])
>>> traj.after(20.0)     # drop transient
Trajectory(n_steps=..., dim=3, t=[20.0, 100.0])
>>> t, y = traj          # tuple-unpack still works

Source code in src/tsdynamics/data/trajectory.py

def __init__(
    self,
    t: np.ndarray,
    y: np.ndarray,
    system: Any,
    meta: dict[str, Any] | None = None,
) -> None:
    self.t = np.asarray(t)
    self.y = np.asarray(y)
    self.system = system
    self.meta = dict(meta) if meta else {}
    self._kdtree: cKDTree | None = None

variables property

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

Component names declared by the system (instance attr or class ClassVar).

dim property

dim: int

State-space dimension.

n_steps property

n_steps: int

Number of time steps.

component

component(i: int | str) -> ndarray

Return a single state component.

PARAMETER	DESCRIPTION
i	Component index, or component name when the system declares variables. TYPE: int or str

RETURNS	DESCRIPTION
(ndarray, shape(T))

Source code in src/tsdynamics/data/trajectory.py

def component(self, i: int | str) -> np.ndarray:
    """
    Return a single state component.

    Parameters
    ----------
    i : int or str
        Component index, or component name when the system declares
        ``variables``.

    Returns
    -------
    ndarray, shape (T,)
    """
    if isinstance(i, str):
        i = self._component_index(i)
    return self.y[:, i]

after

after(t0: float) -> Trajectory

Drop the initial transient.

PARAMETER	DESCRIPTION
t0	Keep only time points t >= t0. TYPE: float

RETURNS	DESCRIPTION
Trajectory

Source code in src/tsdynamics/data/trajectory.py

def after(self, t0: float) -> Trajectory:
    """
    Drop the initial transient.

    Parameters
    ----------
    t0 : float
        Keep only time points ``t >= t0``.

    Returns
    -------
    Trajectory
    """
    mask = self.t >= t0
    return Trajectory(self.t[mask], self.y[mask], self.system, meta=self.meta)

to_plot_spec

to_plot_spec(
    kind: str | None = None,
    *,
    components: int
    | str
    | Sequence[int | str]
    | None = None,
    animate: bool | dict[str, Any] | Animation = False,
    **kind_kw: Any,
) -> PlotSpec

Describe this trajectory as a backend-agnostic :class:PlotSpec.

This is the one front door for trajectory plotting — every common view goes through here, so the parameterised viz.producers builders stay an internal detail.

Auto-dispatch With kind=None the semantic kind follows the number of selected components (after applying components=): 1 → TIME_SERIES, 2 → PHASE_PORTRAIT_2D, 3 → PHASE_PORTRAIT_3D, and 4+ → SPACETIME (a Lorenz-96-style field image, not a misleading 3-D portrait of the first three coordinates). A discrete-map orbit draws with a SCATTER mark (a point sequence) rather than a line. A Poincaré-section trajectory (carrying meta["plot_kind"] of "poincare_section") is recognised and drawn as its in-plane scatter.

Selecting components components= picks what to draw — a name, an index, or a sequence of them: components="x" (a single time series), components= ["y0", "y1", "y2"] (a 3-D portrait of three chosen channels). The auto-dispatch then keys off how many you selected.

Overriding the kind kind= forces any member of the closed :class:~tsdynamics.viz.spec.PlotKind vocabulary — e.g. kind="time_series" to overlay component-vs-time on a 3-D trajectory, or kind="spacetime" to image it. The recipe kind="delay" builds a delay-coordinate embedding x(t) vs x(t - tau); pass tau (in time units) via **kind_kw.

Per-kind options (**kind_kw) Options valid for one kind only are accepted as keywords rather than cluttering the signature — tau (required for kind="delay", converted from time units to samples via meta["dt"]), color_by (time series / phase portraits — a named field "time"/"speed"/"sagitta"/"curvature"/"acceleration"/ "arclength"/"index", a per-point array, or a callable f(trajectory) -> array), transpose (spacetime). Passing one to the wrong kind raises :class:~tsdynamics.errors.InvalidParameterError.

The :mod:tsdynamics.viz import is local to this method (lazy), and the spec carries no rendering code, so building a spec (or importing :mod:tsdynamics) never imports matplotlib / Plotly.

PARAMETER	DESCRIPTION
kind	Override the auto-dispatched kind (a PlotKind value, or the "delay" recipe). None auto-dispatches. TYPE: str DEFAULT: None
components	Which state components to draw (names or indices). None uses all. TYPE: int or str or sequence of int/str DEFAULT: None
animate	Turn the spec into a reveal animation. True uses sensible per-kind defaults (a moving head on portraits / spacetime, off for a plain time series); a dict overrides individual :class:~tsdynamics.viz.spec.Animation fields; an :class:~tsdynamics.viz.spec.Animation is used as-is. Tweak further with the chainable spec.animate() / .trail() / .head() / .camera() / .clock() methods. TYPE: bool or dict or Animation DEFAULT: False
**kind_kw	Per-kind options (see above). TYPE: Any DEFAULT: {}

RETURNS	DESCRIPTION
PlotSpec

Source code in src/tsdynamics/data/trajectory.py

def to_plot_spec(
    self,
    kind: str | None = None,
    *,
    components: int | str | Sequence[int | str] | None = None,
    animate: bool | dict[str, Any] | Animation = False,
    **kind_kw: Any,
) -> PlotSpec:
    """
    Describe this trajectory as a backend-agnostic :class:`PlotSpec`.

    This is the **one front door** for trajectory plotting — every common
    view goes through here, so the parameterised ``viz.producers`` builders
    stay an internal detail.

    Auto-dispatch
        With ``kind=None`` the semantic kind follows the number of selected
        components (after applying ``components=``): 1 → ``TIME_SERIES``,
        2 → ``PHASE_PORTRAIT_2D``, 3 → ``PHASE_PORTRAIT_3D``, and **4+ →
        ``SPACETIME``** (a Lorenz-96-style field image, *not* a misleading
        3-D portrait of the first three coordinates).  A discrete-map orbit
        draws with a ``SCATTER`` mark (a point sequence) rather than a line.
        A Poincaré-section trajectory (carrying ``meta["plot_kind"]`` of
        ``"poincare_section"``) is recognised and drawn as its in-plane
        scatter.

    Selecting components
        ``components=`` picks what to draw — a name, an index, or a sequence
        of them: ``components="x"`` (a single time series), ``components=
        ["y0", "y1", "y2"]`` (a 3-D portrait of three chosen channels).  The
        auto-dispatch then keys off how many you selected.

    Overriding the kind
        ``kind=`` forces any member of the closed
        :class:`~tsdynamics.viz.spec.PlotKind` vocabulary — e.g.
        ``kind="time_series"`` to overlay component-vs-time on a 3-D
        trajectory, or ``kind="spacetime"`` to image it.  The recipe
        ``kind="delay"`` builds a delay-coordinate embedding ``x(t)`` vs
        ``x(t - tau)``; pass ``tau`` (in **time units**) via ``**kind_kw``.

    Per-kind options (``**kind_kw``)
        Options valid for one kind only are accepted as keywords rather than
        cluttering the signature — ``tau`` (required for ``kind="delay"``,
        converted from time units to samples via ``meta["dt"]``),
        ``color_by`` (time series / phase portraits — a named field
        ``"time"``/``"speed"``/``"sagitta"``/``"curvature"``/``"acceleration"``/
        ``"arclength"``/``"index"``, a per-point array, or a callable
        ``f(trajectory) -> array``), ``transpose`` (spacetime).  Passing one to
        the wrong kind raises
        :class:`~tsdynamics.errors.InvalidParameterError`.

    The :mod:`tsdynamics.viz` import is local to this method (lazy), and the
    spec carries no rendering code, so building a spec (or importing
    :mod:`tsdynamics`) never imports matplotlib / Plotly.

    Parameters
    ----------
    kind : str, optional
        Override the auto-dispatched kind (a ``PlotKind`` value, or the
        ``"delay"`` recipe).  ``None`` auto-dispatches.
    components : int or str or sequence of int/str, optional
        Which state components to draw (names or indices).  ``None`` uses all.
    animate : bool or dict or Animation, optional
        Turn the spec into a reveal animation.  ``True`` uses sensible per-kind
        defaults (a moving head on portraits / spacetime, off for a plain time
        series); a dict overrides individual
        :class:`~tsdynamics.viz.spec.Animation` fields; an
        :class:`~tsdynamics.viz.spec.Animation` is used as-is.  Tweak further
        with the chainable ``spec.animate()`` / ``.trail()`` / ``.head()`` /
        ``.camera()`` / ``.clock()`` methods.
    **kind_kw
        Per-kind options (see above).

    Returns
    -------
    PlotSpec
    """
    from tsdynamics.viz.spec import PlotKind

    all_names = self.variables or tuple(f"y{i}" for i in range(self.dim))

    # A Poincaré section carries its intent in meta; honour it before the
    # dimensionality dispatch (only for the unmodified default view).
    if (
        kind is None
        and components is None
        and not kind_kw
        and str(self.meta.get("plot_kind", "")) == PlotKind.POINCARE_SECTION
    ):
        return self._with_animation(self._poincare_section_spec(all_names), animate)

    # The ``"field"`` / ``"spatial_field"`` recipe routes before component
    # resolution: here ``components=`` selects a *field block* (e.g. Gray–
    # Scott's "u"/"v"), not a state component, so it must not be resolved
    # against the per-cell labels.
    if kind is not None and _KIND_ALIASES.get(kind, kind) == "spatial_field":
        self._validate_kind_kw("spatial_field", kind_kw)
        field = self._spatial_field_spec(components)
        return self._with_animation(field, animate)

    sel = self._resolve_components(components, all_names)
    sel_names = tuple(all_names[i] for i in sel)
    n_sel = len(sel)
    discrete = self._is_discrete()

    # Resolve the routing key: a friendly alias (``"delay"``) first, else the
    # auto kind from the number of selected components.
    route = _KIND_ALIASES.get(kind, kind) if kind is not None else _auto_route(n_sel)
    self._validate_kind_kw(route, kind_kw)

    if route == "delay_embedding":
        delay = self._delay_spec(sel, sel_names, kind_kw, explicit=components is not None)
        return self._with_animation(delay, animate)

    spec_kind = PlotKind(route)  # "delay" never reaches here (aliased above)
    ys = self.y[:, sel]

    if spec_kind == PlotKind.SPACETIME:
        image = self._spacetime_spec(ys, sel_names, transpose=bool(kind_kw.get("transpose")))
        return self._with_animation(image, animate)

    color_by = kind_kw.get("color_by")
    if spec_kind == PlotKind.TIME_SERIES:
        series = self._time_series_spec(sel, ys, sel_names, discrete, color_by)
        return self._with_animation(series, animate)

    portrait = self._phase_portrait_spec(spec_kind, sel, ys, sel_names, discrete, color_by)
    return self._with_animation(portrait, animate)

plot

plot(backend: str | None = None, **kwargs: Any) -> Any

Render this trajectory via a visualization backend.

Sugar over :meth:to_plot_spec: the spec-shaping keywords (kind, components, and the per-kind options tau / color_by / transpose) are peeled off and passed to :meth:to_plot_spec; the remaining keywords are inline spec tweaks (xlabel / yscale / title / …) or backend keyword arguments (see :meth:tsdynamics.viz.spec.Plottable.plot).

The viz package is imported lazily here (not at module scope) so plain import tsdynamics never pulls it in — honouring the no-backend-on-import contract. Raises :class:~tsdynamics.viz.spec.VisualizationNotInstalled until a backend is registered.

Source code in src/tsdynamics/data/trajectory.py

def plot(self, backend: str | None = None, **kwargs: Any) -> Any:
    """Render this trajectory via a visualization backend.

    Sugar over :meth:`to_plot_spec`: the spec-shaping keywords (``kind``,
    ``components``, and the per-kind options ``tau`` / ``color_by`` /
    ``transpose``) are peeled off and passed to :meth:`to_plot_spec`; the
    remaining keywords are inline spec tweaks (``xlabel`` / ``yscale`` /
    ``title`` / …) or backend keyword arguments (see
    :meth:`tsdynamics.viz.spec.Plottable.plot`).

    The viz package is imported lazily here (not at module scope) so plain
    ``import tsdynamics`` never pulls it in — honouring the
    no-backend-on-import contract.  Raises
    :class:`~tsdynamics.viz.spec.VisualizationNotInstalled` until a backend
    is registered.
    """
    from tsdynamics.viz.spec import _apply_inline_tweaks

    spec_kw = {k: kwargs.pop(k) for k in list(kwargs) if k in _PLOT_SPEC_KEYS}
    spec = self.to_plot_spec(**spec_kw)
    backend_kw = _apply_inline_tweaks(spec, kwargs)
    return spec.render(backend, **backend_kw)

minmax

minmax() -> tuple[ndarray, ndarray]

Return per-component (minima, maxima), each of shape (dim,).

Source code in src/tsdynamics/data/trajectory.py

def minmax(self) -> tuple[np.ndarray, np.ndarray]:
    """Return per-component ``(minima, maxima)``, each of shape ``(dim,)``."""
    return self.y.min(axis=0), self.y.max(axis=0)

standardize

standardize() -> Trajectory

Return a copy with zero mean and unit standard deviation per component.

The applied transform is recorded in meta["standardized"].

Source code in src/tsdynamics/data/trajectory.py

def standardize(self) -> Trajectory:
    """
    Return a copy with zero mean and unit standard deviation per component.

    The applied transform is recorded in ``meta["standardized"]``.
    """
    mean = self.y.mean(axis=0)
    std = self.y.std(axis=0)
    std = np.where(std < np.finfo(float).tiny, 1.0, std)
    return Trajectory(
        self.t,
        (self.y - mean) / std,
        self.system,
        meta={**self.meta, "standardized": {"mean": mean, "std": std}},
    )

neighbors

neighbors(q: Any, k: int = 1) -> tuple[ndarray, ndarray]

Nearest trajectory points to query point(s) q.

Builds a KD-tree lazily on first call and caches it; subsequent queries are O(log T).

PARAMETER	DESCRIPTION
q	Query point(s). TYPE: (array - like, shape(dim) or (m, dim))
k	Number of neighbours per query point. TYPE: int DEFAULT: 1

RETURNS	DESCRIPTION
(distances, indices)	As returned by :meth:scipy.spatial.cKDTree.query.

Source code in src/tsdynamics/data/trajectory.py

def neighbors(self, q: Any, k: int = 1) -> tuple[np.ndarray, np.ndarray]:
    """
    Nearest trajectory points to query point(s) ``q``.

    Builds a KD-tree lazily on first call and caches it; subsequent
    queries are O(log T).

    Parameters
    ----------
    q : array-like, shape (dim,) or (m, dim)
        Query point(s).
    k : int
        Number of neighbours per query point.

    Returns
    -------
    (distances, indices)
        As returned by :meth:`scipy.spatial.cKDTree.query`.
    """
    from scipy.spatial import cKDTree

    if self._kdtree is None:
        self._kdtree = cKDTree(self.y)
    return cast(
        "tuple[np.ndarray, np.ndarray]",
        self._kdtree.query(np.asarray(q, dtype=float), k=k),
    )

set_distance

set_distance(
    other: Any, *, method: str = "centroid"
) -> float

Distance to another point set (Trajectory or array), as a set.

method is "centroid" (default), "hausdorff", or "minimum" — see :func:tsdynamics.data.set_distance. The matching primitive behind attractor deduplication and continuation.

Source code in src/tsdynamics/data/trajectory.py

def set_distance(self, other: Any, *, method: str = "centroid") -> float:
    """
    Distance to another point set (Trajectory or array), as a set.

    ``method`` is ``"centroid"`` (default), ``"hausdorff"``, or
    ``"minimum"`` — see :func:`tsdynamics.data.set_distance`.  The
    matching primitive behind attractor deduplication and continuation.
    """
    from tsdynamics.data import set_distance

    return set_distance(
        self, other, method=cast('Literal["centroid", "hausdorff", "minimum"]', method)
    )

SystemBase

SystemBase(
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
)

Bases: SystemPlottable

Abstract base class for all dynamical systems.

Provides: - params — a :class:ParamSet holding the system's parameter values. Attribute access on the system is transparently forwarded to params. - dim — integer state-space dimension. - ic — optional initial conditions array. - meta — dict for storing computed metadata (Lyapunov spectra, etc.). - copy() / with_params() for safe cloning. - resolve_ic() for uniform IC resolution across subclasses.

Class-level declarations

Subclasses should declare at class level::

class Lorenz(ContinuousSystem):
    params = {"sigma": 10.0, "rho": 28.0, "beta": 8/3}
    dim = 3

Constructor overrides

Individual instances can override params and/or ic::

lor = Lorenz(params={"rho": 30.0}, ic=[1.0, 0.0, 0.0])

The constructor raises :class:~tsdynamics.errors.InvalidParameterError (a ValueError subclass) for any unknown parameter key, so a typo such as params={"rhoo": 30.0} fails loudly instead of being silently ignored.

See Also

ParamSet : the fixed-key parameter container behind params. MetaStore : the append-with-history store behind meta. resolve_ic : the uniform initial-condition resolution helper.

Initialise a system from its class defaults plus instance overrides.

PARAMETER	DESCRIPTION
params	Per-instance parameter overrides. Every key must already exist in the class-level :attr:params defaults; unknown keys raise. TYPE: dict DEFAULT: None
ic	Initial conditions. Stored on self.ic (as a float array) and used by :meth:resolve_ic when no explicit ic is later supplied. TYPE: array - like DEFAULT: None
dim	State-space dimension override for variable-dimension systems. Falls back to the class-level :attr:dim when omitted. TYPE: int DEFAULT: None
field_shape	Spatial grid shape override for a spatially-extended system (see :attr:_field_shape). Falls back to the class-level value. TYPE: tuple of int DEFAULT: None

RAISES	DESCRIPTION
InvalidParameterError	If params contains a key that is not a declared parameter.

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

def __init__(
    self,
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
) -> None:
    """Initialise a system from its class defaults plus instance overrides.

    Parameters
    ----------
    params : dict, optional
        Per-instance parameter overrides.  Every key must already exist in
        the class-level :attr:`params` defaults; unknown keys raise.
    ic : array-like, optional
        Initial conditions.  Stored on ``self.ic`` (as a ``float`` array) and
        used by :meth:`resolve_ic` when no explicit ``ic`` is later supplied.
    dim : int, optional
        State-space dimension override for variable-dimension systems.  Falls
        back to the class-level :attr:`dim` when omitted.
    field_shape : tuple of int, optional
        Spatial grid shape override for a spatially-extended system (see
        :attr:`_field_shape`).  Falls back to the class-level value.

    Raises
    ------
    InvalidParameterError
        If ``params`` contains a key that is not a declared parameter.
    """
    # Build ParamSet from class defaults + constructor overrides
    defaults = dict(type(self).params)
    if params:
        unknown = set(params) - set(defaults)
        if unknown:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: unknown parameter(s) "
                f"{sorted(unknown)}. Declared: {sorted(defaults)}"
            )
        defaults.update(params)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    object.__setattr__(self, "dim", resolved_dim)

    # field_shape (spatially-extended systems): constructor arg > class
    # attribute.  Set via object.__setattr__ so a custom-N instance overrides
    # the class default without tripping the ClassVar assignment rule (mirrors
    # ``dim``).  ``_provenance`` reads the instance value onto ``traj.meta``.
    resolved_field_shape = field_shape if field_shape is not None else type(self)._field_shape
    object.__setattr__(self, "_field_shape", resolved_field_shape)

    # Initial conditions
    ic_arr = np.asarray(ic, dtype=float) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)

    # Metadata store: computed properties (Lyapunov, etc.) accumulate here
    # with history — repeated runs append instead of overwriting.
    object.__setattr__(self, "meta", MetaStore())

lyap property

lyap: Any

Lyapunov-exponent estimators bound to this system.

A cached :class:~tsdynamics.families._accessors.LyapunovAccessor exposing .spectrum() / .maximal() / .from_data() — each delegating to :func:tsdynamics.analysis.lyapunov_spectrum, :func:~tsdynamics.analysis.max_lyapunov and :func:~tsdynamics.analysis.lyapunov_from_data with this system bound.

chaos property

chaos: Any

Chaos indicators bound to this system.

A cached :class:~tsdynamics.families._accessors.ChaosAccessor exposing .gali() / .expansion_entropy() / .zero_one() — delegating to :func:tsdynamics.analysis.gali, :func:~tsdynamics.analysis.expansion_entropy and :func:~tsdynamics.analysis.zero_one_test.

dims property

dims: Any

Fractal-dimension estimators bound to this system.

A cached :class:~tsdynamics.families._accessors.DimensionsAccessor (.correlation() / .generalized() / …) delegating to the *_dimension free functions. These consume a point set; omitting the data argument runs the system first (an implicit integration).

recurrence property

recurrence: Any

Recurrence-quantification estimators bound to this system.

A cached :class:~tsdynamics.families._accessors.RecurrenceAccessor (.matrix() / .rqa() / .windowed()) delegating to :func:tsdynamics.analysis.recurrence_matrix, :func:~tsdynamics.analysis.rqa and :func:~tsdynamics.analysis.windowed_rqa.

entropy property

entropy: Any

Entropy / complexity estimators bound to this system.

A cached :class:~tsdynamics.families._accessors.EntropyAccessor (.permutation() / .sample() / …) delegating to the entropy free functions. These consume a scalar series; omitting data runs the system first.

surrogate property

surrogate: Any

Surrogate generators + nonlinearity tests bound to this system.

A cached :class:~tsdynamics.families._accessors.SurrogateAccessor (.test() / .generate() / …) delegating to :func:tsdynamics.analysis.surrogate_test, :func:~tsdynamics.analysis.surrogates and the surrogate statistics.

copy

copy() -> SystemBase

Return a deep copy with the same class, params, and ic.

The copy has its own independent params and meta stores, so mutating the clone's parameters or recording metadata on it never affects the original.

RETURNS	DESCRIPTION
SystemBase	A fresh instance of the same subclass with copied params and ic and an empty meta store.

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

def copy(self) -> SystemBase:
    """
    Return a deep copy with the same class, params, and ic.

    The copy has its own independent ``params`` and ``meta`` stores, so
    mutating the clone's parameters or recording metadata on it never
    affects the original.

    Returns
    -------
    SystemBase
        A fresh instance of the same subclass with copied ``params`` and
        ``ic`` and an empty ``meta`` store.
    """
    return type(self)(
        params=cast(ParamSet, self.params).as_dict(),
        ic=self.ic.copy() if self.ic is not None else None,
    )

with_params

with_params(**overrides: Any) -> SystemBase

Return a new system with some parameters overridden.

Does not mutate self. Designed for parameter sweeps::

for rho in np.linspace(0, 50, 200):
    traj = base_system.with_params(rho=rho).integrate(final_time=50)

PARAMETER	DESCRIPTION
**overrides	New parameter values. Keys must exist in params. TYPE: Any DEFAULT: {}

RETURNS	DESCRIPTION
SystemBase	New instance of the same subclass.

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

def with_params(self, **overrides: Any) -> SystemBase:
    """
    Return a **new** system with some parameters overridden.

    Does not mutate ``self``.  Designed for parameter sweeps::

        for rho in np.linspace(0, 50, 200):
            traj = base_system.with_params(rho=rho).integrate(final_time=50)

    Parameters
    ----------
    **overrides
        New parameter values.  Keys must exist in ``params``.

    Returns
    -------
    SystemBase
        New instance of the same subclass.
    """
    new_p = {**cast(ParamSet, self.params).as_dict(), **overrides}
    return type(self)(params=new_p, ic=self.ic)

resolve_ic

resolve_ic(ic: Any | None = None) -> ndarray

Resolve initial conditions consistently.

Priority:

1. ic argument (if provided)
2. self.ic (set by a previous integration / iteration)
3. type(self).default_ic (class-level default, if declared)
4. Random U[0, 1)^dim

The resolved IC is stored in self.ic so subsequent calls without an explicit ic reproduce the same initial state.

PARAMETER	DESCRIPTION
ic	TYPE: array - like or None DEFAULT: None

RETURNS	DESCRIPTION
(ndarray, shape(dim))

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

def resolve_ic(self, ic: Any | None = None) -> np.ndarray:
    """
    Resolve initial conditions consistently.

    Priority:

    1. ``ic`` argument (if provided)
    2. ``self.ic`` (set by a previous integration / iteration)
    3. ``type(self).default_ic`` (class-level default, if declared)
    4. Random ``U[0, 1)^dim``

    The resolved IC is stored in ``self.ic`` so subsequent calls without
    an explicit ``ic`` reproduce the same initial state.

    Parameters
    ----------
    ic : array-like or None

    Returns
    -------
    ndarray, shape (dim,)
    """
    if ic is not None:
        arr = np.asarray(ic, dtype=float).reshape(self.dim)
    elif self.ic is not None:
        arr = np.asarray(self.ic, dtype=float).reshape(self.dim)
    elif type(self).default_ic is not None:
        arr = np.asarray(type(self).default_ic, dtype=float).reshape(self.dim)
    else:
        arr = np.random.rand(cast(int, self.dim))
    object.__setattr__(self, "ic", arr.copy())
    return arr

fixed_points

fixed_points(**kwargs: Any) -> Any

Find fixed points / equilibria of this system.

Delegates to :func:tsdynamics.analysis.fixed_points with this system bound — returns the same list of :class:~tsdynamics.analysis.FixedPoint.

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

def fixed_points(self, **kwargs: Any) -> Any:
    """Find fixed points / equilibria of this system.

    Delegates to :func:`tsdynamics.analysis.fixed_points` with this system
    bound — returns the same list of
    :class:`~tsdynamics.analysis.FixedPoint`.
    """
    from tsdynamics.analysis import fixed_points

    return fixed_points(self, **kwargs)

poincare

poincare(
    section: Any = None,
    at: float = 0.0,
    *,
    plane: tuple[Any, ...] | None = None,
    direction: int = +1,
    **kwargs: Any,
) -> Any

Build a :class:~tsdynamics.derived.PoincareMap of this flow.

The friendly section= (a component index or name) + at= (the crossing value) spelling is sugar over the wrapper's plane tuple; an explicit plane=(normal, offset) may be passed instead for an arbitrary-normal plane. Calling .run(...) (or .trajectory(...)) on the returned map collects crossings — the returned object is exactly PoincareMap(self, plane, direction=...).

PARAMETER	DESCRIPTION
section	State component whose level set defines the section. A string is resolved against the system's variables. Ignored when an explicit plane is given. TYPE: int or str DEFAULT: None
at	The crossing value for section (the plane offset). TYPE: float DEFAULT: 0.0
plane	The raw (component_index, value) or (normal, offset) tuple passed straight to :class:~tsdynamics.derived.PoincareMap. Takes precedence over section / at. TYPE: tuple DEFAULT: None
direction	Crossing direction (sign). TYPE: int DEFAULT: +1
**kwargs	Forwarded to :class:~tsdynamics.derived.PoincareMap (dt, max_time). TYPE: Any DEFAULT: {}

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

def poincare(
    self,
    section: Any = None,
    at: float = 0.0,
    *,
    plane: tuple[Any, ...] | None = None,
    direction: int = +1,
    **kwargs: Any,
) -> Any:
    """Build a :class:`~tsdynamics.derived.PoincareMap` of this flow.

    The friendly ``section=`` (a component index or name) + ``at=`` (the
    crossing value) spelling is sugar over the wrapper's ``plane`` tuple; an
    explicit ``plane=(normal, offset)`` may be passed instead for an
    arbitrary-normal plane.  Calling ``.run(...)`` (or ``.trajectory(...)``)
    on the returned map collects crossings — the returned object is exactly
    ``PoincareMap(self, plane, direction=...)``.

    Parameters
    ----------
    section : int or str, optional
        State component whose level set defines the section.  A string is
        resolved against the system's ``variables``.  Ignored when an
        explicit ``plane`` is given.
    at : float, default 0.0
        The crossing value for ``section`` (the plane offset).
    plane : tuple, optional
        The raw ``(component_index, value)`` or ``(normal, offset)`` tuple
        passed straight to :class:`~tsdynamics.derived.PoincareMap`.  Takes
        precedence over ``section`` / ``at``.
    direction : int, default +1
        Crossing direction (sign).
    **kwargs
        Forwarded to :class:`~tsdynamics.derived.PoincareMap` (``dt``,
        ``max_time``).
    """
    from tsdynamics.derived import PoincareMap

    if plane is None:
        from tsdynamics.errors import InvalidParameterError

        if section is None:
            raise InvalidParameterError(
                "poincare() needs either `section=` (with `at=`) or an explicit `plane=`."
            )
        comp = section
        if isinstance(comp, str):
            names = getattr(type(self), "variables", None)
            if names is None:
                raise InvalidParameterError(
                    f"{type(self).__name__} declares no `variables`; "
                    f"pass an integer `section=` (component index)."
                )
            comp = names.index(comp)
        plane = (int(comp), float(at))
    return PoincareMap(self, plane, direction=direction, **kwargs)

stroboscope

stroboscope(
    period: float | None = None, **kwargs: Any
) -> Any

Build a :class:~tsdynamics.derived.StroboscopicMap of this forced flow.

When period is omitted the forcing period is inferred from the system — from a forcing_period / drive_period hook (used verbatim) or a drive_frequency / omega hook (taken as the angular drive frequency, so the period is 2*pi / omega). The catalogue's forced systems (e.g. :class:~tsdynamics.systems.Duffing, whose autonomising phase obeys zdot = omega) follow the omega convention, so ForcedDuffing().stroboscope() just works. Pass period= to override the inference (or when no drive hook exists). Equivalent to StroboscopicMap(self, period).

PARAMETER	DESCRIPTION
period	The forcing period. When None (the default) it is inferred from the system; a system with no recognised drive hook raises, asking for an explicit period=. TYPE: float DEFAULT: None
**kwargs	Forwarded to :class:~tsdynamics.derived.StroboscopicMap. TYPE: Any DEFAULT: {}

RAISES	DESCRIPTION
InvalidParameterError	When period is omitted and the system exposes no drive hook to infer it from.

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

def stroboscope(self, period: float | None = None, **kwargs: Any) -> Any:
    """Build a :class:`~tsdynamics.derived.StroboscopicMap` of this forced flow.

    When ``period`` is omitted the forcing period is **inferred from the
    system** — from a ``forcing_period`` / ``drive_period`` hook (used
    verbatim) or a ``drive_frequency`` / ``omega`` hook (taken as the
    angular drive frequency, so the period is ``2*pi / omega``).  The
    catalogue's forced systems (e.g. :class:`~tsdynamics.systems.Duffing`,
    whose autonomising phase obeys ``zdot = omega``) follow the ``omega``
    convention, so ``ForcedDuffing().stroboscope()`` just works.  Pass
    ``period=`` to override the inference (or when no drive hook exists).
    Equivalent to ``StroboscopicMap(self, period)``.

    Parameters
    ----------
    period : float, optional
        The forcing period.  When ``None`` (the default) it is inferred from
        the system; a system with no recognised drive hook raises, asking for
        an explicit ``period=``.
    **kwargs
        Forwarded to :class:`~tsdynamics.derived.StroboscopicMap`.

    Raises
    ------
    InvalidParameterError
        When ``period`` is omitted and the system exposes no drive hook to
        infer it from.
    """
    from tsdynamics.derived import StroboscopicMap

    if period is None:
        from tsdynamics.errors import invalid_value
        from tsdynamics.families._accessors import infer_forcing_period

        try:
            period = infer_forcing_period(self)
        except KeyError as err:
            raise invalid_value(
                f"period for {type(self).__name__}.stroboscope()",
                value=None,
                rule="could not be inferred from the system",
                hint=(
                    "pass an explicit `period=` (e.g. `2 * np.pi / omega`), or give "
                    "the system a `drive_frequency` / `forcing_period` attribute."
                ),
            ) from err
    return StroboscopicMap(self, period, **kwargs)

tangent

tangent(k: int | None = None, **kwargs: Any) -> Any

Build a :class:~tsdynamics.derived.TangentSystem (state plus k deviation vectors).

Equivalent to TangentSystem(self, k, ...) — the Lyapunov engine.

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

def tangent(self, k: int | None = None, **kwargs: Any) -> Any:
    """Build a :class:`~tsdynamics.derived.TangentSystem` (state plus ``k`` deviation vectors).

    Equivalent to ``TangentSystem(self, k, ...)`` — the Lyapunov engine.
    """
    from tsdynamics.derived import TangentSystem

    return TangentSystem(self, k, **kwargs)

project

project(*components: Any, **kwargs: Any) -> Any

Build a :class:~tsdynamics.derived.ProjectedSystem onto components.

Accepts component indices or names (resolved against variables), as positional arguments (self.project("x", "z")) or a single sequence (self.project(["x", "z"])). Equivalent to ProjectedSystem(self, components).

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

def project(self, *components: Any, **kwargs: Any) -> Any:
    """Build a :class:`~tsdynamics.derived.ProjectedSystem` onto ``components``.

    Accepts component indices or names (resolved against ``variables``), as
    positional arguments (``self.project("x", "z")``) or a single sequence
    (``self.project(["x", "z"])``).  Equivalent to
    ``ProjectedSystem(self, components)``.
    """
    from tsdynamics.derived import ProjectedSystem

    if len(components) == 1 and not isinstance(components[0], (str, bytes)):
        first = components[0]
        try:
            comps = list(first)
        except TypeError:
            comps = [first]
    else:
        comps = list(components)
    return ProjectedSystem(self, comps, **kwargs)

ensemble

ensemble(states: Any) -> Any

Build an :class:~tsdynamics.derived.EnsembleSystem over states.

Equivalent to EnsembleSystem(self, states) — many copies stepped in lockstep.

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

def ensemble(self, states: Any) -> Any:
    """Build an :class:`~tsdynamics.derived.EnsembleSystem` over ``states``.

    Equivalent to ``EnsembleSystem(self, states)`` — many copies stepped in
    lockstep.
    """
    from tsdynamics.derived import EnsembleSystem

    return EnsembleSystem(self, states)

ContinuousSystem

ContinuousSystem(
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
)

Bases: SystemBase, ABC

Base class for ODE-based dynamical systems, integrated on the engine.

Subclass contract

1. Declare params = {...} and dim = N at class level.
2. Implement _equations as a @staticmethod returning a length-dim sequence of SymEngine symbolic expressions.
3. Optionally mark integer or loop-structural parameters in _structural_params — these are baked into the lowered tape rather than exposed as runtime control parameters.

Lowering

Each system is lowered once to an in-process IR tape with no warmup; the engine reads non-structural parameters live from the system on every run, so a parameter change never triggers a re-lowering.

Class-level attributes

_structural_params : frozenset[str] Parameter names that appear as integer loop bounds or affect the symbolic structure of _equations. These are baked in at compile time. For most systems this is empty (the default).

Example — Lorenz96 uses ``N`` to build the list comprehension::

    _structural_params = frozenset({"N"})

_default_method : str Default integrator name (default "RK45").

Examples:

>>> lor = Lorenz()
>>> traj = lor.integrate(final_time=100, dt=0.01)
>>> t, y = traj          # tuple-unpack
>>> lor.sigma = 15.0     # change param — zero recompile cost
>>> traj2 = lor.integrate(final_time=100)

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

def __init__(
    self,
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
) -> None:
    """Initialise a system from its class defaults plus instance overrides.

    Parameters
    ----------
    params : dict, optional
        Per-instance parameter overrides.  Every key must already exist in
        the class-level :attr:`params` defaults; unknown keys raise.
    ic : array-like, optional
        Initial conditions.  Stored on ``self.ic`` (as a ``float`` array) and
        used by :meth:`resolve_ic` when no explicit ``ic`` is later supplied.
    dim : int, optional
        State-space dimension override for variable-dimension systems.  Falls
        back to the class-level :attr:`dim` when omitted.
    field_shape : tuple of int, optional
        Spatial grid shape override for a spatially-extended system (see
        :attr:`_field_shape`).  Falls back to the class-level value.

    Raises
    ------
    InvalidParameterError
        If ``params`` contains a key that is not a declared parameter.
    """
    # Build ParamSet from class defaults + constructor overrides
    defaults = dict(type(self).params)
    if params:
        unknown = set(params) - set(defaults)
        if unknown:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: unknown parameter(s) "
                f"{sorted(unknown)}. Declared: {sorted(defaults)}"
            )
        defaults.update(params)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    object.__setattr__(self, "dim", resolved_dim)

    # field_shape (spatially-extended systems): constructor arg > class
    # attribute.  Set via object.__setattr__ so a custom-N instance overrides
    # the class default without tripping the ClassVar assignment rule (mirrors
    # ``dim``).  ``_provenance`` reads the instance value onto ``traj.meta``.
    resolved_field_shape = field_shape if field_shape is not None else type(self)._field_shape
    object.__setattr__(self, "_field_shape", resolved_field_shape)

    # Initial conditions
    ic_arr = np.asarray(ic, dtype=float) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)

    # Metadata store: computed properties (Lyapunov, etc.) accumulate here
    # with history — repeated runs append instead of overwriting.
    object.__setattr__(self, "meta", MetaStore())

is_discrete property

is_discrete: bool

ODEs are continuous-time systems.

reinit

reinit(
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    method: str | None = None,
    rtol: float = 1e-06,
    atol: float = 1e-09,
    backend: str | None = None,
) -> None

(Re)start the incremental stepper from state u at time t.

PARAMETER	DESCRIPTION
u	Initial state (falls back to self.ic, then random). TYPE: array - like DEFAULT: None
t	Start time (default 0.0). TYPE: float DEFAULT: None
params	Parameter overrides applied (in place) before restarting. TYPE: dict DEFAULT: None
method	Stepper configuration, as in :meth:integrate. TYPE: str | None DEFAULT: None
rtol	Stepper configuration, as in :meth:integrate. TYPE: str | None DEFAULT: None
atol	Stepper configuration, as in :meth:integrate. TYPE: str | None DEFAULT: None
backend	Stepper configuration, as in :meth:integrate. TYPE: str | None DEFAULT: None

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

def reinit(
    self,
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    method: str | None = None,
    rtol: float = 1e-6,
    atol: float = 1e-9,
    backend: str | None = None,
) -> None:
    """
    (Re)start the incremental stepper from state ``u`` at time ``t``.

    Parameters
    ----------
    u : array-like, optional
        Initial state (falls back to ``self.ic``, then random).
    t : float, optional
        Start time (default 0.0).
    params : dict, optional
        Parameter overrides applied (in place) before restarting.
    method, rtol, atol, backend
        Stepper configuration, as in :meth:`integrate`.
    """
    if params:
        for k, v in params.items():
            self.params[k] = v
    t0 = float(t) if t is not None else 0.0
    ic_arr = self.resolve_ic(u)
    from tsdynamics.engine.run import resolve_backend

    # Honour the requested backend through the stepping protocol: ``reference``
    # is the wheel-free pure-Python oracle and a real, supported stepping path
    # (one ``dt`` chunk per ``step`` through the reference ODE integrator), not
    # something to silently coerce to ``"interp"`` (the old behaviour, which
    # made the oracle unreachable via ``reinit``/``step``/``state`` — diagnosis
    # #5).  ``resolve_backend`` raises ``InvalidParameterError`` for an unknown
    # name (it subclasses ``ValueError``).
    self._step_backend = resolve_backend(
        backend if backend is not None else self._default_backend
    )

    from tsdynamics import solvers
    from tsdynamics.engine.problem import ode_problem

    # Lower the tape ONCE here and reuse it for every step() — a sweep reinits
    # thousands of times, so re-lowering per step would dominate the cost.
    # Resolve the step method first so an implicit kernel (bdf/rosenbrock/
    # trbdf2) gets a Jacobian-carrying tape — step() reuses this exact tape,
    # so the Jacobian must be baked in here or the engine refuses the step.
    #
    # ``method="auto"`` is the same a-priori auto-stiffness contract
    # ``integrate``/``ensemble`` honour (FIX-AUTOSTIFF): probe the Jacobian
    # spectrum at the start state and let the solver registry pick the implicit
    # ``bdf`` on a stiff RHS or the explicit ``rk45`` otherwise
    # (:func:`tsdynamics.solvers.recommend`).  Without this branch the stepping
    # path crashed with an opaque "unknown solver method 'auto'" — the one
    # entry point that rejected the advertised value (diagnosis P2-1), exactly
    # where auto-stiffness matters for a stiff system stepped incrementally
    # (Poincaré / basins / streaming).
    self._step_method = method or self._default_method
    if solvers.normalize(self._step_method) == "auto":
        resolution = solvers.recommend(self, family="ode", ic=ic_arr, t=t0)
    else:
        resolution = solvers.resolve(self._step_method)
    # Cache the canonical kernel name (``"RK45"`` → ``"rk45"``) so the per-step
    # loop hits the engine core directly without re-resolving every call.
    self._step_method_canonical = resolution.name
    prob = ode_problem(self, ic=ic_arr, t0=t0, **resolution.build_kwargs)
    # Marshal the (loop-invariant) tape to its engine wire arrays once here so the
    # per-``dt`` ``step`` loop reuses them instead of rebuilding the tuple on every
    # call (stream WS-INVHOIST).  Re-derived on each ``reinit`` (a sweep reinits
    # per parameter value), so a re-lowered tape is always picked up.  Marshal before
    # publishing the problem so a (theoretical) marshalling failure can't leave a new
    # ``_engine_problem`` paired with a stale ``_step_tape_arrays``.
    step_arrays = prob.tape.to_arrays()
    self._engine_problem = prob
    self._step_tape_arrays = step_arrays
    self._step_rtol = float(rtol)
    self._step_atol = float(atol)
    self._state_now = ic_arr.copy()
    self._t_now = t0
    # Drop any prior durable stepper handle: the next ``step`` rebuilds it from
    # the freshly lowered tape + new live state (WS-STEPPER).  Built lazily so a
    # cold ``state()``/``time()`` (or a ``reference`` flow) never forces the
    # compiled engine here — exactly as the pre-handle path reached the engine
    # only inside ``step``.
    self._ode_stepper = None

step

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

Advance the system by dt (default 0.01) and return the new state.

The first call performs an implicit :meth:reinit. Parameter changes made after reinit take effect on the next reinit, not on a live stepper.

Notes

Each call advances exactly one dt from the live (state, t), returning byte-for-byte the trajectory the released per-dt path produced — there is no batching, so the numbers are unchanged for every method, adaptive or fixed-step (streams WS-STEPBUF, WS-INVHOIST, WS-STEPPER). The amortisation is durable: the first step after a :meth:reinit builds an opaque resumable engine handle (:class:tsdynamics._rust.OdeStepper, via :func:~tsdynamics.engine.run.make_ode_stepper) that owns the built tape evaluator + solver once and carries the live (u, t) across calls; every later step is one :func:~tsdynamics.engine.run.step_advance on that handle — the tape is never re-marshalled into the engine again. So a constant-dt stepping loop (Poincaré refinement, basins over flows) skips not only the solver-registry resolve, the implicit-Jacobian decision, the output-grid build, provenance assembly and the :class:Trajectory wrap that the full :meth:integrate entry point pays, but also the per-step tape re-marshalling and tape rebuild the pre-handle stepping core paid. The control-parameter vector is still read live each step, so the live-stepper semantics are unchanged.

Why this stays answer-exact: the engine handle's advance(dt) re-seeds a fresh solver and state for each dt segment (the adaptive controller is re-seeded each step, exactly as the released per-dt integrate_dense did), so the numbers are bit-for-bit identical to that path — verified in the engine's own test. A batch-ahead variant that integrated a whole chunk in one engine call was rejected (WS-STEPBUF): a chunked adaptive integration is not equal to N single-dt integrations (the controller would carry its step/error state across output nodes), which silently corrupted sensitive consumers such as max_lyapunov. The durable handle amortises the build/marshalling, never the numerics.

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

def step(self, n_or_dt: float | None = None) -> np.ndarray:
    """
    Advance the system by ``dt`` (default 0.01) and return the new state.

    The first call performs an implicit :meth:`reinit`.  Parameter changes
    made after ``reinit`` take effect on the next ``reinit``, not on a
    live stepper.

    Notes
    -----
    Each call advances **exactly one** ``dt`` from the live ``(state, t)``,
    returning byte-for-byte the trajectory the released per-``dt`` path
    produced — there is no batching, so the numbers are unchanged for every
    method, adaptive or fixed-step (streams WS-STEPBUF, WS-INVHOIST,
    WS-STEPPER).  The amortisation is durable: the first ``step`` after a
    :meth:`reinit` builds an opaque resumable engine handle
    (:class:`tsdynamics._rust.OdeStepper`, via
    :func:`~tsdynamics.engine.run.make_ode_stepper`) that owns the built tape
    evaluator + solver once and carries the live ``(u, t)`` across calls; every
    later ``step`` is one :func:`~tsdynamics.engine.run.step_advance` on that
    handle — the tape is **never re-marshalled into the engine again**.  So a
    constant-``dt`` stepping loop (Poincaré refinement, basins over flows) skips
    not only the solver-registry resolve, the implicit-Jacobian decision, the
    output-grid build, provenance assembly and the :class:`Trajectory` wrap that
    the full :meth:`integrate` entry point pays, but also the per-step tape
    re-marshalling and tape *rebuild* the pre-handle stepping core paid.  The
    control-parameter vector is still read live each step, so the live-stepper
    semantics are unchanged.

    Why this stays answer-exact: the engine handle's ``advance(dt)`` re-seeds a
    *fresh* solver and state for each ``dt`` segment (the adaptive controller is
    re-seeded each step, exactly as the released per-``dt`` ``integrate_dense``
    did), so the numbers are bit-for-bit identical to that path — verified in
    the engine's own test.  A batch-ahead variant that integrated a whole chunk
    in one engine call was rejected (WS-STEPBUF): a chunked adaptive integration
    is *not* equal to N single-``dt`` integrations (the controller would carry
    its step/error state across output nodes), which silently corrupted
    sensitive consumers such as ``max_lyapunov``.  The durable handle amortises
    the *build/marshalling*, never the numerics.
    """
    from tsdynamics.engine.run import make_ode_stepper, step_advance

    if self._engine_problem is None:
        self.reinit()
    # After ``reinit`` (run above when cold) these stepping-state attributes are
    # always populated; narrow them for the typed engine call below.
    assert self._state_now is not None
    assert self._step_method_canonical is not None
    dt = float(n_or_dt) if n_or_dt is not None else self._default_step_dt

    # The wheel-free pure-Python oracle: advance one ``dt`` chunk through the
    # reference ODE integrator (the same path ``integrate(backend="reference")``
    # uses), so the protocol exposes the oracle honestly instead of secretly
    # running the compiled engine (diagnosis #5).  Off the durable engine-handle
    # fast path (reference owns no ``OdeStepper``); it is the validation backend,
    # not a hot loop.
    if self._step_backend == "reference":
        return self._step_reference(dt)

    t0 = self._t_now
    tf = t0 + dt
    # Preserve the released ``step`` span contract exactly: a non-positive ``dt``
    # or a non-forward window (``t0 + dt == t0`` at a large ``t0``) must raise the
    # canonical :class:`~tsdynamics.errors.InvalidParameterError` (a
    # ``ValueError``), not silently no-op on the engine handle.  When the span
    # clears ``1e-9`` the regime is unambiguously forward+positive, so the happy
    # path goes straight to the handle (never touching ``make_output_grid``);
    # only the sub-``1e-9`` remainder defers to the helper for that identical
    # loud-footgun error (and the byte-identical degenerate-grid behaviour).
    if not tf - t0 > 1e-9:
        from tsdynamics.utils.grids import make_output_grid

        # Raises InvalidParameterError for dt <= 0 / a non-forward window; for a
        # valid-but-tiny span it returns the (possibly single-node) grid, which
        # the per-``dt`` engine core integrated identically — reproduce that here
        # via the same lean core to stay byte-identical for the rare small step.
        t_eval = make_output_grid(t0, tf, dt)
        from tsdynamics.engine.run import _step_continuous

        y = _step_continuous(
            self._step_tape_arrays,
            self._state_now,
            self._engine_problem.params_vec(),
            t_eval,
            method=self._step_method_canonical,
            rtol=self._step_rtol,
            atol=self._step_atol,
            jit=self._step_backend == "jit",
            name=type(self).__name__,
        )
        state = np.asarray(y[-1], dtype=float)
        # The handle (if already built) no longer mirrors the live point after
        # this off-handle advance; drop it so the next ``step`` rebuilds from the
        # synced ``_state_now``/``_t_now``.
        object.__setattr__(self, "_ode_stepper", None)
        object.__setattr__(self, "_t_now", tf)
        object.__setattr__(self, "_state_now", state.copy())
        return state.copy()

    # Build the durable resumable handle lazily on the first ``step`` after a
    # ``reinit`` (so a cold ``state()``/``time()`` never forces the engine), from
    # the live state + the tape arrays + the solver config ``reinit`` cached.
    if self._ode_stepper is None:
        stepper = make_ode_stepper(
            self._step_tape_arrays,
            self._state_now,
            self._t_now,
            method=self._step_method_canonical,
            rtol=self._step_rtol,
            atol=self._step_atol,
            jit=self._step_backend == "jit",
        )
        object.__setattr__(self, "_ode_stepper", stepper)

    state = step_advance(
        self._ode_stepper,
        dt,
        self._engine_problem.params_vec(),
        name=type(self).__name__,
    )
    # The state/time advance writes private framework attributes that always pass
    # straight through ``SystemBase.__setattr__`` (underscore-prefixed) — go direct
    # to ``object.__setattr__`` so the hot loop skips the param-typo guard's
    # ``params`` membership check on every step (WS-INVHOIST).  The engine handle
    # is the live-state authority; ``_state_now``/``_t_now`` mirror it so
    # ``state()``/``time()`` and ``set_state`` stay consistent.
    object.__setattr__(self, "_t_now", self._t_now + dt)
    object.__setattr__(self, "_state_now", state.copy())
    return state.copy()

state

state() -> ndarray

Return a copy of the current state (implicit reinit if cold).

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

def state(self) -> np.ndarray:
    """Return a copy of the current state (implicit ``reinit`` if cold)."""
    if self._state_now is None:
        self.reinit()
    assert self._state_now is not None  # set by reinit above
    return self._state_now.copy()

set_state

set_state(u: Any) -> None

Overwrite the current state without changing the current time.

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

def set_state(self, u: Any) -> None:
    """Overwrite the current state without changing the current time."""
    u_arr = np.asarray(u, dtype=float).reshape(self.dim)
    if self._engine_problem is None:
        self.reinit(u_arr)
    else:
        self._state_now = u_arr.copy()
        # Drop the durable stepper handle: its live point no longer mirrors the
        # reseated state, so the next ``step`` rebuilds it from ``_state_now`` /
        # ``_t_now`` (WS-STEPPER).  Cheaper and less error-prone than reseating
        # the handle in place, and ``set_state`` is not on the hot stepping path.
        self._ode_stepper = None

time

time() -> float

Return the current stepper time.

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

def time(self) -> float:
    """Return the current stepper time."""
    return self._t_now

trajectory

trajectory(
    final_time: float = 100.0,
    *,
    dt: float = 0.02,
    transient: float = 0.0,
    **kwargs: Any,
) -> Trajectory

Protocol-uniform trajectory: integrate plus optional transient drop.

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

def trajectory(
    self,
    final_time: float = 100.0,
    *,
    dt: float = 0.02,
    transient: float = 0.0,
    **kwargs: Any,
) -> Trajectory:
    """Protocol-uniform trajectory: ``integrate`` plus optional transient drop."""
    traj = self.integrate(final_time=transient + final_time, dt=dt, **kwargs)
    return traj.after(transient) if transient > 0 else traj

jacobian_sym

jacobian_sym() -> list[list[Any]]

Return the symbolic Jacobian of _equations, differentiated by SymEngine.

Rows are d f_i / d y(j) for the current structural parameters; non-structural parameters appear as symbols. Hand-written _jacobian methods on system classes are never used at runtime — this autogenerated form is the single source of truth (the test suite cross-checks hand-written ones against it).

RETURNS	DESCRIPTION
list of ``dim`` rows, each a list of ``dim`` SymEngine expressions.

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

def jacobian_sym(self) -> list[list[Any]]:
    """
    Return the symbolic Jacobian of ``_equations``, differentiated by SymEngine.

    Rows are ``d f_i / d y(j)`` for the *current* structural parameters;
    non-structural parameters appear as symbols.  Hand-written
    ``_jacobian`` methods on system classes are never used at runtime —
    this autogenerated form is the single source of truth (the test suite
    cross-checks hand-written ones against it).

    Returns
    -------
    list of ``dim`` rows, each a list of ``dim`` SymEngine expressions.
    """
    import symengine

    from tsdynamics.engine.symbols import state_time_symbols

    y, t_sym = state_time_symbols()

    dim = cast(int, self.dim)
    struct_vals = self._structural_vals()
    control_syms = {k: symengine.Symbol(k) for k in self._control_params()}
    f_sym = list(type(self)._equations(y, t_sym, **{**struct_vals, **control_syms}))
    if len(f_sym) != dim:
        raise InvalidParameterError(
            f"_equations must return {dim} expressions, got {len(f_sym)}"
        )
    return [
        [_resolve_derivative_nodes(symengine.sympify(fi).diff(y(j))) for j in range(dim)]
        for fi in f_sym
    ]

jacobian

jacobian(u: Any, t: float = 0.0) -> ndarray

Evaluate the (autogenerated) Jacobian numerically at state u.

PARAMETER	DESCRIPTION
u	State at which to evaluate. TYPE: (array - like, shape(dim))
t	Time (matters only for non-autonomous systems). TYPE: float DEFAULT: 0.0

RETURNS	DESCRIPTION
(ndarray, shape(dim, dim))

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

def jacobian(self, u: Any, t: float = 0.0) -> np.ndarray:
    """
    Evaluate the (autogenerated) Jacobian numerically at state ``u``.

    Parameters
    ----------
    u : array-like, shape (dim,)
        State at which to evaluate.
    t : float
        Time (matters only for non-autonomous systems).

    Returns
    -------
    ndarray, shape (dim, dim)
    """
    dim = cast(int, self.dim)
    _, jac_fn, control_names = self._build_lambdified()
    vals = [float(self.params[k]) for k in control_names]
    arg = np.concatenate([np.asarray(u, dtype=float).ravel(), [t], vals])
    return np.asarray(jac_fn(arg), dtype=float).reshape(dim, dim)

run

run(
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    events: Any = None,
    **kwargs: Any,
) -> Trajectory

Produce a trajectory — the one canonical verb for every family.

run is the unified trajectory producer: it answers the same call for flows, maps, DDEs and SDEs, dispatching on :attr:is_discrete. For a continuous-time system (this family) it integrates the flow, so run is a thin alias of :meth:integrate and forwards every keyword to it unchanged.

PARAMETER	DESCRIPTION
final_time	End of the integration window. Default 100.0. TYPE: float DEFAULT: 100.0
dt	Output sampling interval. The internal stepper is adaptive. TYPE: float DEFAULT: 0.02
events	Detect events along the flow (a SciPy-shaped events= API). Each element is an :class:~tsdynamics.engine.run.Event, a bare g(y, t) callable carrying .direction / .terminal attributes (the SciPy convention), or a plane tuple (("y", 0.0, "up")). A terminal event stops the integration at its first crossing (arbitrary stopping). The returned trajectory carries each event's crossings in meta["t_events"] / meta["y_events"] (one array per event, aligned with events), plus meta["terminated"]. This wires the same compiled event engine :class:~tsdynamics.derived.poincare.PoincareMap uses; PoincareMap.as_events() shows the section as one such event. TYPE: sequence DEFAULT: None
**kwargs	Forwarded verbatim to :meth:integrate (t0, ic, method, rtol, atol, backend, …). TYPE: Any DEFAULT: {}

RETURNS	DESCRIPTION
Trajectory	Identical to :meth:integrate when events is None — run adds no behaviour. With events set, the dense trajectory (truncated at the first terminal crossing) plus the per-event crossings in meta.

See Also

integrate : The family-specific spelling (a permanent alias of run).

Examples:

>>> traj = Lorenz().run(final_time=100, dt=0.01)
>>> Henon().run(n=5000)            # the same verb iterates a map
>>> sol = Lorenz().run(final_time=50, events=[("z", 27.0, "up")])
>>> sol.meta["t_events"][0].shape    # times z=27 was crossed upward
(... ,)

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

def run(
    self,
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    events: Any = None,
    **kwargs: Any,
) -> Trajectory:
    """
    Produce a trajectory — the one canonical verb for every family.

    ``run`` is the unified trajectory producer: it answers the same call for
    flows, maps, DDEs and SDEs, dispatching on :attr:`is_discrete`.  For a
    continuous-time system (this family) it integrates the flow, so
    ``run`` is a thin alias of :meth:`integrate` and forwards every keyword
    to it unchanged.

    Parameters
    ----------
    final_time : float
        End of the integration window. Default 100.0.
    dt : float
        Output sampling interval. The internal stepper is adaptive.
    events : sequence, optional
        Detect events along the flow (a SciPy-shaped ``events=`` API).  Each
        element is an :class:`~tsdynamics.engine.run.Event`, a bare
        ``g(y, t)`` callable carrying ``.direction`` / ``.terminal``
        attributes (the SciPy convention), or a plane tuple
        (``("y", 0.0, "up")``).  A **terminal** event stops the integration at
        its first crossing (arbitrary stopping).  The returned trajectory
        carries each event's crossings in ``meta["t_events"]`` /
        ``meta["y_events"]`` (one array per event, aligned with ``events``),
        plus ``meta["terminated"]``.  This wires the same compiled event
        engine :class:`~tsdynamics.derived.poincare.PoincareMap` uses;
        ``PoincareMap.as_events()`` shows the section as one such event.
    **kwargs
        Forwarded verbatim to :meth:`integrate` (``t0``, ``ic``, ``method``,
        ``rtol``, ``atol``, ``backend``, …).

    Returns
    -------
    Trajectory
        Identical to :meth:`integrate` when ``events`` is ``None`` — ``run``
        adds no behaviour.  With ``events`` set, the dense trajectory
        (truncated at the first terminal crossing) plus the per-event
        crossings in ``meta``.

    See Also
    --------
    integrate : The family-specific spelling (a permanent alias of ``run``).

    Examples
    --------
    >>> traj = Lorenz().run(final_time=100, dt=0.01)
    >>> Henon().run(n=5000)            # the same verb iterates a map
    >>> sol = Lorenz().run(final_time=50, events=[("z", 27.0, "up")])
    >>> sol.meta["t_events"][0].shape    # times z=27 was crossed upward
    (... ,)
    """
    return self.integrate(final_time=final_time, dt=dt, events=events, **kwargs)

integrate

integrate(
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    method: str | None = None,
    rtol: float = 1e-06,
    atol: float = 1e-09,
    backend: str | None = None,
    events: Any = None,
    **integrator_kwargs: Any,
) -> Trajectory

Integrate the ODE and return a :class:~tsdynamics.families.Trajectory.

PARAMETER	DESCRIPTION
final_time	End of integration window. Default 100.0. TYPE: float DEFAULT: 100.0
dt	Output sampling interval. The internal stepper is adaptive. TYPE: float DEFAULT: 0.02
t0	Start time. Default 0.0. Allows warm restarts from a non-zero time (the IC is interpreted as the state at t0). TYPE: float DEFAULT: 0.0
ic	Initial state at t0. Falls back to self.ic, then U[0, 1)^dim. TYPE: array - like DEFAULT: None
method	Solver name, resolved by the solver registry (default "RK45"): explicit (RK45 / DOP853 / tsit5 / dop853) or implicit / stiff (bdf / rosenbrock / trbdf2). Pass "auto" to select a kernel by a-priori auto-stiffness — the Jacobian spectrum at the start state is probed and bdf chosen on a stiff RHS, rk45 otherwise (:func:tsdynamics.solvers.recommend; a one-point heuristic, so a reliably-stiff system should still declare _default_method). TYPE: str DEFAULT: None
rtol	Solver tolerances (default 1e-6 / 1e-9). TYPE: float DEFAULT: 1e-06
atol	Solver tolerances (default 1e-6 / 1e-9). TYPE: float DEFAULT: 1e-06
backend	Where the ODE is integrated. Defaults to _default_backend ("interp"). "interp" / "jit" — the Rust engine (the zero-warmup SSA-tape interpreter or the Cranelift JIT) via the shared engine seam (:func:tsdynamics.engine.run.integrate). "reference" — the dependency-light pure-Python oracle (the lowered tape integrated with SciPy); the engine's validation backend, usable without the compiled wheel. TYPE: ('interp', 'jit', 'reference') DEFAULT: "interp"
events	Detect events along the flow (the SciPy-shaped events= API; see :meth:run). Each element is an :class:~tsdynamics.engine.run.Event, a bare g(y, t) callable carrying .direction / .terminal attributes, or a plane tuple (("y", 0.0, "up")). A terminal event stops the integration at its first crossing; the returned trajectory carries each event's crossings in meta["t_events"] / meta["y_events"] (aligned with events) plus meta["terminated"]. TYPE: sequence DEFAULT: None

RETURNS	DESCRIPTION
Trajectory	Supports tuple-unpacking: t, y = sys.integrate(...).

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

def integrate(
    self,
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    method: str | None = None,
    rtol: float = 1e-6,
    atol: float = 1e-9,
    backend: str | None = None,
    events: Any = None,
    **integrator_kwargs: Any,
) -> Trajectory:
    """
    Integrate the ODE and return a :class:`~tsdynamics.families.Trajectory`.

    Parameters
    ----------
    final_time : float
        End of integration window. Default 100.0.
    dt : float
        Output sampling interval. The internal stepper is adaptive.
    t0 : float
        Start time. Default 0.0. Allows warm restarts from a non-zero
        time (the IC is interpreted as the state at ``t0``).
    ic : array-like, optional
        Initial state at ``t0``. Falls back to ``self.ic``, then
        ``U[0, 1)^dim``.
    method : str, optional
        Solver name, resolved by the solver registry (default ``"RK45"``):
        explicit (``RK45`` / ``DOP853`` / ``tsit5`` / ``dop853``) or implicit
        / stiff (``bdf`` / ``rosenbrock`` / ``trbdf2``).  Pass ``"auto"`` to
        select a kernel by a-priori auto-stiffness — the Jacobian spectrum at
        the start state is probed and ``bdf`` chosen on a stiff RHS, ``rk45``
        otherwise (:func:`tsdynamics.solvers.recommend`; a one-point heuristic,
        so a reliably-stiff system should still declare ``_default_method``).
    rtol, atol : float
        Solver tolerances (default 1e-6 / 1e-9).
    backend : {"interp", "jit", "reference"}, optional
        Where the ODE is integrated.  Defaults to ``_default_backend``
        (``"interp"``).

        - ``"interp"`` / ``"jit"`` — the **Rust engine** (the zero-warmup
          SSA-tape interpreter or the Cranelift JIT) via the shared engine
          seam (:func:`tsdynamics.engine.run.integrate`).
        - ``"reference"`` — the dependency-light pure-Python oracle (the
          lowered tape integrated with SciPy); the engine's validation
          backend, usable without the compiled wheel.
    events : sequence, optional
        Detect events along the flow (the SciPy-shaped ``events=`` API; see
        :meth:`run`).  Each element is an
        :class:`~tsdynamics.engine.run.Event`, a bare ``g(y, t)`` callable
        carrying ``.direction`` / ``.terminal`` attributes, or a plane tuple
        (``("y", 0.0, "up")``).  A **terminal** event stops the integration at
        its first crossing; the returned trajectory carries each event's
        crossings in ``meta["t_events"]`` / ``meta["y_events"]`` (aligned with
        ``events``) plus ``meta["terminated"]``.

    Returns
    -------
    Trajectory
        Supports tuple-unpacking: ``t, y = sys.integrate(...)``.
    """
    if integrator_kwargs:
        # Anything left in **integrator_kwargs is an unrecognised keyword (every
        # valid argument is bound to an explicit parameter above). Reject it
        # instead of silently dropping a typo'd keyword (the WS-ERRADOPT footgun).
        from tsdynamics.errors import invalid_value

        bad = sorted(integrator_kwargs)[0]
        raise invalid_value(
            bad,
            integrator_kwargs[bad],
            rule="is not a valid integrate()/run() keyword",
            hint="check the keyword spelling (e.g. final_time, dt, t0, ic, method, rtol, atol).",
        )
    if events is not None:
        return self._run_events(
            final_time=final_time,
            dt=dt,
            events=events,
            t0=t0,
            ic=ic,
            method=method,
            rtol=rtol,
            atol=atol,
            backend=backend,
        )
    backend = backend if backend is not None else self._default_backend
    return self._dispatch(
        backend=backend,
        final_time=final_time,
        dt=dt,
        t0=t0,
        ic=ic,
        method=method or self._default_method,
        rtol=rtol,
        atol=atol,
    )

lyapunov_spectrum

lyapunov_spectrum(
    final_time: float = 200.0,
    dt: float = 0.1,
    *,
    ic: Any | None = None,
    n_exp: int | None = None,
    burn_in: float = 50.0,
    method: str | None = None,
    rtol: float = 1e-06,
    atol: float = 1e-09,
    backend: str = "interp",
    **integrator_kwargs: Any,
) -> ndarray

Estimate the Lyapunov spectrum of the flow.

Delegates to :class:~tsdynamics.derived.tangent.TangentSystem, the one backend-neutral variational/Lyapunov engine shared across families: the extended variational ODE (state ⊕ k tangent vectors) is integrated on the chosen backend per dt-chunk and QR-reorthonormalised. The Benettin time-averaging of the log-stretch rates follows the classical construction of Benettin et al. [1]_.

Results are stored in self.meta['lyapunov_spectrum'].

PARAMETER	DESCRIPTION
final_time	Averaging window length after burn-in. Default 200.0. TYPE: float DEFAULT: 200.0
dt	Sampling interval for local exponent accumulation. Default 0.1. TYPE: float DEFAULT: 0.1
ic	Initial state. Falls back to self.ic, then random. TYPE: array - like DEFAULT: None
n_exp	Number of exponents. Defaults to dim. TYPE: int DEFAULT: None
burn_in	Discard this much time before averaging. Default 50.0. TYPE: float DEFAULT: 50.0
method	Integrator (default "RK45"). TYPE: str DEFAULT: None
rtol	Tolerances. TYPE: float DEFAULT: 1e-06
atol	Tolerances. TYPE: float DEFAULT: 1e-06
backend	Backend on which the extended variational ODE is integrated. Defaults to "interp" (the zero-warmup Rust engine interpreter). "reference" is the dependency-light pure-Python oracle (usable without the compiled wheel); "jit" is the Cranelift JIT. Any other name is rejected by :class:~tsdynamics.derived.tangent.TangentSystem. TYPE: ('interp', 'jit', 'reference') DEFAULT: "interp"

RETURNS	DESCRIPTION
(ndarray, shape(n_exp))	Lyapunov exponents ordered from largest to smallest.

RAISES	DESCRIPTION
InvalidParameterError	If n_exp is given and not a positive integer.

References

.. [1] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, "Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them," Meccanica 15, 9-30 (1980).

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

def lyapunov_spectrum(
    self,
    final_time: float = 200.0,
    dt: float = 0.1,
    *,
    ic: Any | None = None,
    n_exp: int | None = None,
    burn_in: float = 50.0,
    method: str | None = None,
    rtol: float = 1e-6,
    atol: float = 1e-9,
    backend: str = "interp",
    **integrator_kwargs: Any,
) -> np.ndarray:
    """
    Estimate the Lyapunov spectrum of the flow.

    Delegates to :class:`~tsdynamics.derived.tangent.TangentSystem`, the one
    backend-neutral variational/Lyapunov engine shared across families: the
    *extended* variational ODE (state ⊕ ``k`` tangent vectors) is integrated
    on the chosen ``backend`` per dt-chunk and QR-reorthonormalised.  The
    Benettin time-averaging of the log-stretch rates follows the classical
    construction of Benettin et al. [1]_.

    Results are stored in ``self.meta['lyapunov_spectrum']``.

    Parameters
    ----------
    final_time : float
        Averaging window length after burn-in. Default 200.0.
    dt : float
        Sampling interval for local exponent accumulation. Default 0.1.
    ic : array-like, optional
        Initial state. Falls back to ``self.ic``, then random.
    n_exp : int, optional
        Number of exponents. Defaults to ``dim``.
    burn_in : float
        Discard this much time before averaging. Default 50.0.
    method : str, optional
        Integrator (default ``"RK45"``).
    rtol, atol : float
        Tolerances.
    backend : {"interp", "jit", "reference"}, optional
        Backend on which the extended variational ODE is integrated.
        Defaults to ``"interp"`` (the zero-warmup Rust engine interpreter).
        ``"reference"`` is the dependency-light pure-Python oracle (usable
        without the compiled wheel); ``"jit"`` is the Cranelift JIT.  Any
        other name is rejected by :class:`~tsdynamics.derived.tangent.TangentSystem`.

    Returns
    -------
    ndarray, shape (n_exp,)
        Lyapunov exponents ordered from largest to smallest.

    Raises
    ------
    InvalidParameterError
        If ``n_exp`` is given and not a positive integer.

    References
    ----------
    .. [1] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn,
       "Lyapunov characteristic exponents for smooth dynamical systems and
       for Hamiltonian systems; a method for computing all of them,"
       *Meccanica* 15, 9-30 (1980).
    """
    if n_exp is not None and n_exp <= 0:
        raise InvalidParameterError(f"n_exp must be a positive integer, got {n_exp!r}")
    from tsdynamics.derived.tangent import TangentSystem

    k = n_exp if n_exp is not None else self.dim
    return TangentSystem(self, k=k, backend=backend).lyapunov_spectrum(
        final_time=final_time,
        dt=dt,
        ic=ic,
        burn_in=burn_in,
        method=method,
        rtol=rtol,
        atol=atol,
        **integrator_kwargs,
    )

DelaySystem

DelaySystem(
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
)

Bases: SystemBase, ABC

Base class for delay differential systems (DDEs), integrated on the engine.

Subclass contract

1. Declare params = {...} and dim = N.
2. Implement _equations as a @staticmethod returning a length-dim sequence of SymEngine symbolic expressions. Use y(i, t - tau) for delayed state access.

Lowering

Each system is lowered once to an in-process IR tape with no warmup. Delay values directly affect the history-buffer structure, so they are baked into the tape rather than read live like the other parameters; a delay change re-lowers, while ordinary parameters are read live with no re-lowering.

DDEs typically need looser tolerances than ODEs (start with rtol=atol=1e-3).

History

Pass a history callable h(s) → sequence defining the past for s ≤ 0. If omitted, a constant past equal to ic is used.

.. note:: Provide a non-equilibrium history to avoid trivial Lyapunov exponents. lyapunov_spectrum starts from a constant past, so the workaround is to run integrate first with the desired history, then pass the end-state as ic to lyapunov_spectrum.

Examples:

>>> mg = MackeyGlass()
>>> hist = lambda s: [1.0 + 0.1 * np.sin(0.2 * s)]
>>> traj = mg.integrate(final_time=500, history=hist)
>>> exps = mg.lyapunov_spectrum(n_exp=2, ic=traj.y[-1])

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

def __init__(
    self,
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
) -> None:
    """Initialise a system from its class defaults plus instance overrides.

    Parameters
    ----------
    params : dict, optional
        Per-instance parameter overrides.  Every key must already exist in
        the class-level :attr:`params` defaults; unknown keys raise.
    ic : array-like, optional
        Initial conditions.  Stored on ``self.ic`` (as a ``float`` array) and
        used by :meth:`resolve_ic` when no explicit ``ic`` is later supplied.
    dim : int, optional
        State-space dimension override for variable-dimension systems.  Falls
        back to the class-level :attr:`dim` when omitted.
    field_shape : tuple of int, optional
        Spatial grid shape override for a spatially-extended system (see
        :attr:`_field_shape`).  Falls back to the class-level value.

    Raises
    ------
    InvalidParameterError
        If ``params`` contains a key that is not a declared parameter.
    """
    # Build ParamSet from class defaults + constructor overrides
    defaults = dict(type(self).params)
    if params:
        unknown = set(params) - set(defaults)
        if unknown:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: unknown parameter(s) "
                f"{sorted(unknown)}. Declared: {sorted(defaults)}"
            )
        defaults.update(params)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    object.__setattr__(self, "dim", resolved_dim)

    # field_shape (spatially-extended systems): constructor arg > class
    # attribute.  Set via object.__setattr__ so a custom-N instance overrides
    # the class default without tripping the ClassVar assignment rule (mirrors
    # ``dim``).  ``_provenance`` reads the instance value onto ``traj.meta``.
    resolved_field_shape = field_shape if field_shape is not None else type(self)._field_shape
    object.__setattr__(self, "_field_shape", resolved_field_shape)

    # Initial conditions
    ic_arr = np.asarray(ic, dtype=float) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)

    # Metadata store: computed properties (Lyapunov, etc.) accumulate here
    # with history — repeated runs append instead of overwriting.
    object.__setattr__(self, "meta", MetaStore())

is_discrete property

is_discrete: bool

DDEs are continuous-time systems.

reinit

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

(Re)start the incremental stepper from a constant past equal to u.

DDE state is a history function; the protocol restart uses a constant past (the same convention as lyapunov_spectrum). For a custom history, use :meth:integrate with history= and continue from traj.y[-1].

Stepping is forward-only and re-integrates from the constant past on the Rust DDE engine each call (the method of steps has no stateful one-step restart), so it is correct but O(steps²) — use :meth:integrate for a full trajectory.

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

def reinit(
    self,
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    rtol: float | None = None,
    atol: float | None = None,
    **kwargs: Any,
) -> None:
    """
    (Re)start the incremental stepper from a constant past equal to ``u``.

    DDE state is a history *function*; the protocol restart uses a
    constant past (the same convention as ``lyapunov_spectrum``).  For a
    custom history, use :meth:`integrate` with ``history=`` and continue
    from ``traj.y[-1]``.

    Stepping is forward-only and re-integrates from the constant past on the
    Rust DDE engine each call (the method of steps has no stateful one-step
    restart), so it is correct but ``O(steps²)`` — use :meth:`integrate` for
    a full trajectory.
    """
    if params:
        for k, v in params.items():
            self.params[k] = v
    if t is not None and float(t) != 0.0:
        raise NotImplementedError(
            "DelaySystem.reinit only supports t=0 (the past starts there)."
        )
    self._past_ic = self.resolve_ic(u)
    self._step_rtol = rtol
    self._step_atol = atol
    self._state_now = self._past_ic.copy()
    self._t_now = 0.0

step

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

Advance by dt (default 0.1, forward-only) and return the new state.

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

def step(self, n_or_dt: float | None = None) -> np.ndarray:
    """Advance by ``dt`` (default 0.1, forward-only) and return the new state."""
    if self._past_ic is None:
        self.reinit()
    dt = float(n_or_dt) if n_or_dt is not None else self._default_step_dt
    self._t_now = self._t_now + dt
    traj = self.integrate(
        final_time=self._t_now,
        dt=min(dt, self._t_now),
        ic=self._past_ic,
        rtol=self._step_rtol if self._step_rtol is not None else self._default_rtol,
        atol=self._step_atol if self._step_atol is not None else self._default_atol,
    )
    state = np.asarray(traj.y[-1], dtype=float)
    if not np.isfinite(state).all():
        raise ConvergenceError(
            f"{type(self).__name__}: DDE diverged at t={self._t_now:.6g} during step()."
        )
    self._state_now = state.copy()
    return state

state

state() -> ndarray

Return a copy of the current state (implicit reinit if cold).

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

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

set_state

set_state(u: Any) -> None

Not available for DDEs — their state is a whole history function.

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

def set_state(self, u: Any) -> None:
    """Not available for DDEs — their state is a whole history function."""
    raise NotImplementedError(
        f"{type(self).__name__}.set_state is impossible for delay systems: the "
        f"instantaneous state is a history function over [t - max_delay, t], not a "
        f"point.  Use reinit(u) to restart from a constant past, or integrate(...) "
        f"with a history callable."
    )

time

time() -> float

Return the current stepper time.

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

def time(self) -> float:
    """Return the current stepper time."""
    return self._t_now

trajectory

trajectory(
    final_time: float = 100.0,
    *,
    dt: float = 0.02,
    transient: float = 0.0,
    **kwargs: Any,
) -> Trajectory

Protocol-uniform trajectory: integrate plus optional transient drop.

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

def trajectory(
    self,
    final_time: float = 100.0,
    *,
    dt: float = 0.02,
    transient: float = 0.0,
    **kwargs: Any,
) -> Trajectory:
    """Protocol-uniform trajectory: ``integrate`` plus optional transient drop."""
    traj = self.integrate(final_time=transient + final_time, dt=dt, **kwargs)
    return traj.after(transient) if transient > 0 else traj

run

run(
    final_time: float = 100.0,
    dt: float = 0.02,
    **kwargs: Any,
) -> Trajectory

Produce a trajectory — the one canonical verb for every family.

run is the unified trajectory producer: it answers the same call for flows, maps, DDEs and SDEs, dispatching on :attr:is_discrete. For a delay system (this family) it integrates the DDE, so run is a thin alias of :meth:integrate and forwards every keyword to it unchanged.

PARAMETER	DESCRIPTION
final_time	Integration end time. Default 100.0. TYPE: float DEFAULT: 100.0
dt	Output sampling interval. TYPE: float DEFAULT: 0.02
**kwargs	Forwarded verbatim to :meth:integrate (ic, history, rtol, atol, backend, method). TYPE: Any DEFAULT: {}

RETURNS	DESCRIPTION
Trajectory	Identical to :meth:integrate — run adds no behaviour.

See Also

integrate : The family-specific spelling (a permanent alias of run).

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

def run(
    self,
    final_time: float = 100.0,
    dt: float = 0.02,
    **kwargs: Any,
) -> Trajectory:
    """
    Produce a trajectory — the one canonical verb for every family.

    ``run`` is the unified trajectory producer: it answers the same call for
    flows, maps, DDEs and SDEs, dispatching on :attr:`is_discrete`.  For a
    delay system (this family) it integrates the DDE, so ``run`` is a thin
    alias of :meth:`integrate` and forwards every keyword to it unchanged.

    Parameters
    ----------
    final_time : float
        Integration end time. Default 100.0.
    dt : float
        Output sampling interval.
    **kwargs
        Forwarded verbatim to :meth:`integrate` (``ic``, ``history``,
        ``rtol``, ``atol``, ``backend``, ``method``).

    Returns
    -------
    Trajectory
        Identical to :meth:`integrate` — ``run`` adds no behaviour.

    See Also
    --------
    integrate : The family-specific spelling (a permanent alias of ``run``).
    """
    return self.integrate(final_time=final_time, dt=dt, **kwargs)

integrate

integrate(
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    ic: Any | None = None,
    history: History = None,
    rtol: float | None = None,
    atol: float | None = None,
    backend: str | None = None,
    method: str = "rk45",
    **kwargs: Any,
) -> Trajectory

Integrate the DDE and return a :class:~tsdynamics.families.Trajectory.

PARAMETER	DESCRIPTION
final_time	Integration end time. Default 100.0. TYPE: float DEFAULT: 100.0
dt	Output sampling interval. TYPE: float DEFAULT: 0.02
ic	Used for constant past when history is None. Falls back to self.ic, then random. TYPE: array - like DEFAULT: None
history	h(s) → sequence of length dim for s ≤ 0. If None, a constant past equal to ic is used. TYPE: callable DEFAULT: None
rtol	Integration tolerances. DDEs typically need 1e-3; very tight tolerances can stall the solver. TYPE: float DEFAULT: None
atol	Integration tolerances. DDEs typically need 1e-3; very tight tolerances can stall the solver. TYPE: float DEFAULT: None
backend	Which engine integrates the DDE. Defaults to _default_backend ("interp"). "interp" / "jit" route — through the shared engine seam (:func:tsdynamics.engine.run.integrate) — to the Rust method-of-steps engine (history ring buffer + cubic-Hermite dense interpolation; stream E-DDE), reusing the explicit solver kernels. Only constant delays lower; a state-dependent delay raises. backend="reference" is unsupported for DDEs (there is no pure-Python delay integrator). TYPE: str DEFAULT: None
method	The explicit kernel ("rk45", "tsit5", "dop853", "rk4"); the method of steps drives explicit kernels only. "auto" is an explicit no-op here — it resolves to the DDE default (rk45) without an auto-stiffness probe. Auto-stiffness is an ODE-only feature: the one-point heuristic reads only the instantaneous Jacobian and would ignore the delay terms that shape a DDE's spectrum (it could even select an implicit kernel the method-of-steps engine cannot drive), so pass another explicit method= directly if you need one. TYPE: str DEFAULT: "rk45"

RETURNS	DESCRIPTION
Trajectory

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

def integrate(
    self,
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    ic: Any | None = None,
    history: History = None,
    rtol: float | None = None,
    atol: float | None = None,
    backend: str | None = None,
    method: str = "rk45",
    **kwargs: Any,
) -> Trajectory:
    """
    Integrate the DDE and return a :class:`~tsdynamics.families.Trajectory`.

    Parameters
    ----------
    final_time : float
        Integration end time. Default 100.0.
    dt : float
        Output sampling interval.
    ic : array-like, optional
        Used for constant past when ``history`` is ``None``.
        Falls back to ``self.ic``, then random.
    history : callable, optional
        ``h(s) → sequence`` of length ``dim`` for ``s ≤ 0``.
        If ``None``, a constant past equal to ``ic`` is used.
    rtol, atol : float
        Integration tolerances.  DDEs typically need 1e-3; very tight
        tolerances can stall the solver.
    backend : str, optional
        Which engine integrates the DDE.  Defaults to ``_default_backend``
        (``"interp"``).  ``"interp"`` / ``"jit"`` route — through the shared
        engine seam (:func:`tsdynamics.engine.run.integrate`) — to the Rust
        method-of-steps engine (history ring buffer + cubic-Hermite dense
        interpolation; stream E-DDE), reusing the explicit solver kernels.
        Only **constant** delays lower; a state-dependent delay raises.
        ``backend="reference"`` is unsupported for DDEs (there is no
        pure-Python delay integrator).
    method : str, default "rk45"
        The explicit kernel (``"rk45"``, ``"tsit5"``, ``"dop853"``,
        ``"rk4"``); the method of steps drives explicit kernels only.
        ``"auto"`` is an explicit **no-op** here — it resolves to the DDE
        default (``rk45``) without an auto-stiffness probe.  Auto-stiffness is
        an ODE-only feature: the one-point heuristic reads only the
        instantaneous Jacobian and would ignore the delay terms that shape a
        DDE's spectrum (it could even select an implicit kernel the
        method-of-steps engine cannot drive), so pass another explicit
        ``method=`` directly if you need one.

    Returns
    -------
    Trajectory
    """
    backend = backend if backend is not None else self._default_backend
    return self._integrate_engine(
        final_time,
        dt,
        ic=ic,
        history=history,
        rtol=rtol,
        atol=atol,
        backend=backend,
        method=method,
    )

lyapunov_spectrum

lyapunov_spectrum(
    final_time: float = 200.0,
    dt: float = 0.1,
    *,
    ic: Any | None = None,
    n_exp: int = 1,
    burn_in: float = 50.0,
    rtol: float | None = None,
    atol: float | None = None,
    backend: str | None = None,
    **kwargs: Any,
) -> ndarray

Estimate the n_exp leading Lyapunov exponents of the delay system.

The engine estimator (stream E-DDE-LYAP, result stored in self.meta['lyapunov_spectrum']) integrates the extended variational DDE on the Rust engine with a function-space Benettin renormalisation (:func:tsdynamics.families._dde_lyapunov.dde_lyapunov_spectrum): backend="interp" / "jit". "reference" is rejected (the engine has no pure-Python DDE integrator).

PARAMETER	DESCRIPTION
final_time	Averaging window after burn-in. Default 200.0. TYPE: float DEFAULT: 200.0
dt	Sampling interval (should divide the maximum delay). TYPE: float DEFAULT: 0.1
ic	Initial state. Provide the end-state of a prior integrate call so the trajectory starts on the attractor (recommended). TYPE: array - like DEFAULT: None
n_exp	Number of leading exponents to estimate. DDEs have infinitely many; choose consciously. Default 1. TYPE: int DEFAULT: 1
burn_in	Discard interval. Default 50.0. TYPE: float DEFAULT: 50.0
rtol	Integration tolerances. The engine path renormalises every delay window and uses defaults of 1e-7 / 1e-9. TYPE: float DEFAULT: None
atol	Integration tolerances. The engine path renormalises every delay window and uses defaults of 1e-7 / 1e-9. TYPE: float DEFAULT: None
backend	"interp" or "jit". Defaults to :attr:_default_backend ("interp"). TYPE: str DEFAULT: None

Notes

For best results, pass ic=traj.y[-1] from a prior integrate run — this places the trajectory on the attractor and avoids trivial exponents from equilibrium pasts.

RETURNS	DESCRIPTION
(ndarray, shape(n_exp))

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

def lyapunov_spectrum(
    self,
    final_time: float = 200.0,
    dt: float = 0.1,
    *,
    ic: Any | None = None,
    n_exp: int = 1,
    burn_in: float = 50.0,
    rtol: float | None = None,
    atol: float | None = None,
    backend: str | None = None,
    **kwargs: Any,
) -> np.ndarray:
    """
    Estimate the ``n_exp`` leading Lyapunov exponents of the delay system.

    The **engine** estimator (stream E-DDE-LYAP, result stored in
    ``self.meta['lyapunov_spectrum']``) integrates the extended variational
    DDE on the Rust engine with a function-space Benettin renormalisation
    (:func:`tsdynamics.families._dde_lyapunov.dde_lyapunov_spectrum`):
    ``backend="interp"`` / ``"jit"``.  ``"reference"`` is rejected (the engine
    has no pure-Python DDE integrator).

    Parameters
    ----------
    final_time : float
        Averaging window after burn-in. Default 200.0.
    dt : float
        Sampling interval (should divide the maximum delay).
    ic : array-like, optional
        Initial state. Provide the end-state of a prior ``integrate``
        call so the trajectory starts on the attractor (recommended).
    n_exp : int
        Number of leading exponents to estimate. DDEs have infinitely
        many; choose consciously. Default 1.
    burn_in : float
        Discard interval. Default 50.0.
    rtol, atol : float, optional
        Integration tolerances.  The engine path renormalises every delay
        window and uses defaults of ``1e-7`` / ``1e-9``.
    backend : str, optional
        ``"interp"`` or ``"jit"``.  Defaults to :attr:`_default_backend`
        (``"interp"``).

    Notes
    -----
    For best results, pass ``ic=traj.y[-1]`` from a prior ``integrate`` run —
    this places the trajectory on the attractor and avoids trivial exponents
    from equilibrium pasts.

    Returns
    -------
    ndarray, shape (n_exp,)
    """
    backend = backend if backend is not None else self._default_backend
    from tsdynamics.families._dde_lyapunov import dde_lyapunov_spectrum

    if kwargs:
        raise InvalidInputError(
            f"lyapunov_spectrum(backend={backend!r}) does not accept the "
            f"extra integration keyword(s) {sorted(kwargs)}."
        )
    exps = dde_lyapunov_spectrum(
        self,
        n_exp=n_exp,
        final_time=final_time,
        dt=dt,
        burn_in=burn_in,
        ic=ic,
        backend=backend,
        rtol=rtol if rtol is not None else 1e-7,
        atol=atol if atol is not None else 1e-9,
    )
    self.meta.record(
        "lyapunov_spectrum",
        exps,
        backend=backend,
        n_exp=n_exp,
        final_time=final_time,
        dt=dt,
        burn_in=burn_in,
    )
    return exps

DiscreteMap

DiscreteMap(
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
)

Bases: SystemBase

Base class for discrete maps iterated on the engine.

Subclass contract

1. Declare params = {...} and dim = N.
2. Implement _step and _jacobian as @staticmethod static methods. Parameters arrive as positional arguments in the order they appear in the class-level params dict.

Iteration

iterate lowers _step to an in-process IR tape and runs the engine's native map loop, with no warmup. The engine reads the current parameter values live on every run, so a parameter change never triggers a re-lowering.

Lyapunov spectrum

Computed in a single forward pass via QR decomposition of the Jacobian product — no redundant second iteration over the trajectory.

Examples:

>>> h = Henon()
>>> traj = h.iterate(steps=10_000)
>>> t_idx, X = traj          # tuple-unpack
>>> exps = h.lyapunov_spectrum(steps=5_000)
>>> h_variant = h.with_params(a=1.2)
>>> traj2 = h_variant.iterate(steps=10_000)

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

def __init__(
    self,
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
) -> None:
    """Initialise a system from its class defaults plus instance overrides.

    Parameters
    ----------
    params : dict, optional
        Per-instance parameter overrides.  Every key must already exist in
        the class-level :attr:`params` defaults; unknown keys raise.
    ic : array-like, optional
        Initial conditions.  Stored on ``self.ic`` (as a ``float`` array) and
        used by :meth:`resolve_ic` when no explicit ``ic`` is later supplied.
    dim : int, optional
        State-space dimension override for variable-dimension systems.  Falls
        back to the class-level :attr:`dim` when omitted.
    field_shape : tuple of int, optional
        Spatial grid shape override for a spatially-extended system (see
        :attr:`_field_shape`).  Falls back to the class-level value.

    Raises
    ------
    InvalidParameterError
        If ``params`` contains a key that is not a declared parameter.
    """
    # Build ParamSet from class defaults + constructor overrides
    defaults = dict(type(self).params)
    if params:
        unknown = set(params) - set(defaults)
        if unknown:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: unknown parameter(s) "
                f"{sorted(unknown)}. Declared: {sorted(defaults)}"
            )
        defaults.update(params)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    object.__setattr__(self, "dim", resolved_dim)

    # field_shape (spatially-extended systems): constructor arg > class
    # attribute.  Set via object.__setattr__ so a custom-N instance overrides
    # the class default without tripping the ClassVar assignment rule (mirrors
    # ``dim``).  ``_provenance`` reads the instance value onto ``traj.meta``.
    resolved_field_shape = field_shape if field_shape is not None else type(self)._field_shape
    object.__setattr__(self, "_field_shape", resolved_field_shape)

    # Initial conditions
    ic_arr = np.asarray(ic, dtype=float) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)

    # Metadata store: computed properties (Lyapunov, etc.) accumulate here
    # with history — repeated runs append instead of overwriting.
    object.__setattr__(self, "meta", MetaStore())

is_discrete property

is_discrete: bool

Maps are discrete-time systems.

reinit

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

(Re)start stepping from state u at iteration count t.

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

def reinit(
    self,
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
) -> None:
    """(Re)start stepping from state ``u`` at iteration count ``t``."""
    if params:
        for k, v in params.items():
            self.params[k] = v
    self._state_now = self.resolve_ic(u)
    self._n_now = int(t) if t is not None else 0

step

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

Advance n iterations and return the new state.

The first call performs an implicit :meth:reinit.

PARAMETER	DESCRIPTION
n_or_dt	Number of iterations to advance (default 1). Must be a positive whole number — a discrete map has no notion of a fractional step. TYPE: int DEFAULT: None

RETURNS	DESCRIPTION
(ndarray, shape(dim))	A copy of the state after n iterations.

RAISES	DESCRIPTION
InvalidParameterError	If n_or_dt is non-positive or not a whole number.
ConvergenceError	If the orbit diverges to a non-finite state within the n steps.

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

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

    The first call performs an implicit :meth:`reinit`.

    Parameters
    ----------
    n_or_dt : int, optional
        Number of iterations to advance (default 1).  Must be a positive
        whole number — a discrete map has no notion of a fractional step.

    Returns
    -------
    ndarray, shape (dim,)
        A copy of the state after ``n`` iterations.

    Raises
    ------
    InvalidParameterError
        If ``n_or_dt`` is non-positive or not a whole number.
    ConvergenceError
        If the orbit diverges to a non-finite state within the ``n`` steps.
    """
    if self._state_now is None:
        self.reinit()
    if n_or_dt is None:
        n = 1
    else:
        nf = float(n_or_dt)
        if not nf.is_integer() or nf < 1:
            raise InvalidParameterError(
                f"{type(self).__name__}.step takes a positive whole number of "
                f"iterations, got {n_or_dt!r} (fractional time steps have no "
                f"meaning for discrete maps)."
            )
        n = int(nf)
    assert self._state_now is not None
    x = self._state_now
    params = cast("ParamSet", self.params).as_tuple()

    step_fn = type(self)._step
    # A diverging orbit overflows to ``inf`` in pure-Python float64 arithmetic;
    # that is *expected* and is caught by the finite check below, so silence the
    # spurious NumPy over/under/invalid FP warnings the loop would otherwise emit
    # (under ``filterwarnings=error`` they would mask the real divergence signal).
    with np.errstate(all="ignore"):
        for _ in range(n):
            x = np.asarray(step_fn(x, *params), dtype=np.float64).ravel()
    if not np.isfinite(x).all():
        raise ConvergenceError(
            f"{type(self).__name__}: map diverged at iteration {self._n_now + n}."
        )
    self._state_now = np.asarray(x, dtype=float).reshape(self.dim)
    self._n_now += n
    return self._state_now.copy()

state

state() -> ndarray

Return a copy of the current state (implicit reinit if cold).

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

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

set_state

set_state(u: Any) -> None

Overwrite the current state.

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

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

time

time() -> float

Return the current iteration count.

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

def time(self) -> float:
    """Return the current iteration count."""
    return float(self._n_now)

trajectory

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

Protocol-uniform trajectory: iterate plus optional transient drop.

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

def trajectory(
    self,
    steps: int = 1000,
    *,
    transient: int = 0,
    **kwargs: Any,
) -> Trajectory:
    """Protocol-uniform trajectory: ``iterate`` plus optional transient drop."""
    traj = self.iterate(steps=transient + steps, **kwargs)
    return traj[transient:] if transient > 0 else traj

run

run(n: int = 1000, **kwargs: Any) -> Trajectory

Produce a trajectory — the one canonical verb for every family.

run is the unified trajectory producer: it answers the same call for flows, maps, DDEs and SDEs, dispatching on :attr:is_discrete. For a discrete map (this family) it iterates the map, so run is a thin alias of :meth:iterate. The number of iterations is named n (the canonical step-count keyword), forwarded to :meth:iterate as its steps argument; every other keyword is forwarded unchanged.

PARAMETER	DESCRIPTION
n	Number of iterations. Default 1000. TYPE: int DEFAULT: 1000
**kwargs	Forwarded verbatim to :meth:iterate (ic, max_retries, backend). TYPE: Any DEFAULT: {}

RETURNS	DESCRIPTION
Trajectory	Identical to :meth:iterate — run adds no behaviour.

See Also

iterate : The family-specific spelling (a permanent alias of run).

Examples:

>>> Henon().run(n=5000)
>>> Lorenz().run(final_time=100, dt=0.01)   # the same verb integrates a flow

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

def run(
    self,
    n: int = 1000,
    **kwargs: Any,
) -> Trajectory:
    """
    Produce a trajectory — the one canonical verb for every family.

    ``run`` is the unified trajectory producer: it answers the same call for
    flows, maps, DDEs and SDEs, dispatching on :attr:`is_discrete`.  For a
    discrete map (this family) it iterates the map, so ``run`` is a thin
    alias of :meth:`iterate`.  The number of iterations is named ``n`` (the
    canonical step-count keyword), forwarded to :meth:`iterate` as its
    ``steps`` argument; every other keyword is forwarded unchanged.

    Parameters
    ----------
    n : int
        Number of iterations. Default 1000.
    **kwargs
        Forwarded verbatim to :meth:`iterate` (``ic``, ``max_retries``,
        ``backend``).

    Returns
    -------
    Trajectory
        Identical to :meth:`iterate` — ``run`` adds no behaviour.

    See Also
    --------
    iterate : The family-specific spelling (a permanent alias of ``run``).

    Examples
    --------
    >>> Henon().run(n=5000)
    >>> Lorenz().run(final_time=100, dt=0.01)   # the same verb integrates a flow
    """
    return self.iterate(steps=n, **kwargs)

iterate

iterate(
    steps: int = 1000,
    ic: Any | None = None,
    max_retries: int = 10,
    *,
    backend: str | None = None,
) -> Trajectory

Iterate the map for steps steps on the engine.

PARAMETER	DESCRIPTION
steps	Number of iterations. Default 1000. TYPE: int DEFAULT: 1000
ic	Initial state. Falls back to self.ic, then random. TYPE: array - like DEFAULT: None
max_retries	Retry with a new random IC if divergence is detected (only when no explicit ic was given; an explicit ic that diverges raises). TYPE: int DEFAULT: 10
backend	Where the iteration runs. Defaults to _default_backend ("interp"). "interp" (default) / "jit" — the Rust engine's native map loop (interpreter or Cranelift JIT). Requires the compiled extension (:mod:tsdynamics._rust); until it is built these raise :class:~tsdynamics.engine.run.EngineNotAvailableError. "reference" — the lowered next-state tape, iterated in pure Python (the dependency-light oracle the engine is validated against). Every backend lowers _step to the engine IR, so it requires a map whose step traces symbolically (see :func:tsdynamics.engine.compile.lower_map); piecewise or numpy-ufunc steps raise :class:~tsdynamics.engine.compile.TapeCompileError. TYPE: str DEFAULT: None

RETURNS	DESCRIPTION
Trajectory	t is arange(steps) (integer step indices, not float times).

RAISES	DESCRIPTION
ConvergenceError	If an explicit ic diverges, or every random-IC retry diverges.
EngineNotAvailableError	If a Rust-engine backend ("interp" / "jit") is requested but the compiled extension is not built. This propagates immediately (it is not divergence) rather than consuming the retry budget.
TapeCompileError	If _step cannot be lowered to the engine IR (piecewise / ufunc).

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

def iterate(
    self,
    steps: int = 1000,
    ic: Any | None = None,
    max_retries: int = 10,
    *,
    backend: str | None = None,
) -> Trajectory:
    """
    Iterate the map for ``steps`` steps on the engine.

    Parameters
    ----------
    steps : int
        Number of iterations. Default 1000.
    ic : array-like, optional
        Initial state. Falls back to ``self.ic``, then random.
    max_retries : int
        Retry with a new random IC if divergence is detected (only when no
        explicit ``ic`` was given; an explicit ic that diverges raises).
    backend : str, optional
        Where the iteration runs.  Defaults to ``_default_backend``
        (``"interp"``).

        - ``"interp"`` (default) / ``"jit"`` — the Rust engine's native
          map loop (interpreter or Cranelift JIT).  Requires the compiled
          extension (:mod:`tsdynamics._rust`); until it is built these
          raise :class:`~tsdynamics.engine.run.EngineNotAvailableError`.
        - ``"reference"`` — the lowered next-state tape, iterated in pure
          Python (the dependency-light oracle the engine is validated
          against).

        Every backend lowers ``_step`` to the engine IR, so it requires a
        map whose step traces symbolically (see
        :func:`tsdynamics.engine.compile.lower_map`); piecewise or
        ``numpy``-ufunc steps raise
        :class:`~tsdynamics.engine.compile.TapeCompileError`.

    Returns
    -------
    Trajectory
        ``t`` is ``arange(steps)`` (integer step indices, not float times).

    Raises
    ------
    ConvergenceError
        If an explicit ``ic`` diverges, or every random-IC retry diverges.
    EngineNotAvailableError
        If a Rust-engine backend (``"interp"`` / ``"jit"``) is requested but
        the compiled extension is not built.  This propagates immediately
        (it is not divergence) rather than consuming the retry budget.
    TapeCompileError
        If ``_step`` cannot be lowered to the engine IR (piecewise / ufunc).
    """
    backend = backend if backend is not None else self._default_backend

    # Iterate on the Rust engine.  Preserve the random-IC retry only when no
    # explicit ``ic`` was given (a random draw can land off-basin); an
    # explicit ic that diverges raises loudly, the engine's contract.
    ic_explicit = ic is not None
    ic_arr = self.resolve_ic(ic)
    for attempt in range(max_retries):
        try:
            return self._iterate_engine(steps=steps, ic=ic_arr, backend=backend)
        except (ConvergenceError, ArithmeticError) as exc:
            # Catch ONLY divergence — :class:`ConvergenceError` (the engine /
            # reference "diverge loudly" signal, also a ``RuntimeError``) and the
            # arithmetic blow-ups (``OverflowError`` / ``FloatingPointError`` /
            # ``ZeroDivisionError``, all :class:`ArithmeticError`).  A missing /
            # broken engine surfaces as
            # :class:`~tsdynamics.engine.run.EngineNotAvailableError` (a
            # :class:`~tsdynamics.errors.BackendError`, hence a ``RuntimeError`` but
            # NOT a ``ConvergenceError``); narrowing the catch lets it — and any
            # other genuine fault, e.g. a ``backend="jit"`` compile failure —
            # propagate loudly instead of being mistaken for divergence and
            # silently burning the whole retry budget.
            if ic_explicit or attempt == max_retries - 1:
                raise
            # Off-basin random draw diverged; warn (not stdout) and retry from
            # a fresh random IC. Final exhaustion raises loudly below.
            warnings.warn(
                f"{type(self).__name__}.iterate: {exc} "
                "Retrying from a new random initial condition.",
                RuntimeWarning,
                stacklevel=2,
            )
            ic_arr = np.random.rand(cast(int, self.dim))
            object.__setattr__(self, "ic", ic_arr.copy())
    raise ConvergenceError(
        f"{type(self).__name__}.iterate exhausted {max_retries} "
        f"retries without a finite trajectory."
    )

lyapunov_spectrum

lyapunov_spectrum(
    steps: int = 5000,
    ic: Any | None = None,
    n_exp: int | None = None,
    reortho_interval: int = 1,
    *,
    backend: str | None = None,
) -> ndarray

QR-based Lyapunov spectrum.

Delegates to :class:~tsdynamics.derived.tangent.TangentSystem, the one backend-neutral variational/Lyapunov engine shared across families — a single forward pass evaluating the Jacobian alongside the trajectory, QR-reorthonormalising every reortho_interval steps, with a random-IC retry on divergence.

On the compiled-engine backends ("interp" default / "jit") the whole QR tangent-map iteration runs in one Rust kernel call (:func:tsdynamics.engine.run.map_lyapunov) — no per-step Python→FFI round-trip, so it is dramatically faster than the per-step NumPy loop. backend="reference" (and any map whose _step will not lower to the engine IR, or a wheel-free environment) runs the pure-Python QR loop — the oracle the engine is validated against.

Results are stored in self.meta['lyapunov_spectrum'].

PARAMETER	DESCRIPTION
steps	Number of iterations. Default 5000. TYPE: int DEFAULT: 5000
ic	Initial state. Falls back to self.ic, then random. TYPE: array - like DEFAULT: None
n_exp	Number of exponents. Defaults to dim. TYPE: int DEFAULT: None
reortho_interval	Reorthonormalise every this many steps. Default 1. TYPE: int DEFAULT: 1
backend	"interp" (default, the Rust kernel) / "jit" (Cranelift) / "reference" (the pure-Python QR loop). TYPE: str DEFAULT: None

RETURNS	DESCRIPTION
(ndarray, shape(n_exp))	Lyapunov exponents ordered from largest to smallest.

References

.. [1] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, "Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them," Meccanica 15, 9-30 (1980).

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

def lyapunov_spectrum(
    self,
    steps: int = 5000,
    ic: Any | None = None,
    n_exp: int | None = None,
    reortho_interval: int = 1,
    *,
    backend: str | None = None,
) -> np.ndarray:
    """
    QR-based Lyapunov spectrum.

    Delegates to :class:`~tsdynamics.derived.tangent.TangentSystem`, the one
    backend-neutral variational/Lyapunov engine shared across families — a
    single forward pass evaluating the Jacobian alongside the trajectory,
    QR-reorthonormalising every ``reortho_interval`` steps, with a random-IC
    retry on divergence.

    On the compiled-engine backends (``"interp"`` default / ``"jit"``) the whole
    QR tangent-map iteration runs in one Rust kernel call
    (:func:`tsdynamics.engine.run.map_lyapunov`) — no per-step Python→FFI
    round-trip, so it is dramatically faster than the per-step NumPy loop.
    ``backend="reference"`` (and any map whose ``_step`` will not lower to the
    engine IR, or a wheel-free environment) runs the pure-Python QR loop — the
    oracle the engine is validated against.

    Results are stored in ``self.meta['lyapunov_spectrum']``.

    Parameters
    ----------
    steps : int
        Number of iterations. Default 5000.
    ic : array-like, optional
        Initial state. Falls back to ``self.ic``, then random.
    n_exp : int, optional
        Number of exponents. Defaults to ``dim``.
    reortho_interval : int
        Reorthonormalise every this many steps. Default 1.
    backend : str, optional
        ``"interp"`` (default, the Rust kernel) / ``"jit"`` (Cranelift) /
        ``"reference"`` (the pure-Python QR loop).

    Returns
    -------
    ndarray, shape (n_exp,)
        Lyapunov exponents ordered from largest to smallest.

    References
    ----------
    .. [1] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn,
       "Lyapunov characteristic exponents for smooth dynamical systems and
       for Hamiltonian systems; a method for computing all of them,"
       *Meccanica* 15, 9-30 (1980).
    """
    from tsdynamics.derived.tangent import TangentSystem

    k = n_exp or self.dim
    return TangentSystem(self, k=k, backend=backend).lyapunov_spectrum(
        steps=steps, ic=ic, reortho_interval=reortho_interval
    )

StochasticSystem

StochasticSystem(
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
)

Bases: SystemBase, ABC

Base class for diagonal-Itô stochastic differential equations.

Subclass contract

1. Declare params = {...} and dim = N at class level.

2. Implement _drift and _diffusion as @staticmethods, each returning a length-dim sequence of SymEngine symbolic expressions (use y(i) for component i and t for time — no NumPy, no math, no Python if):

3. _drift(y, t, **params) is the deterministic part f;

4. _diffusion(y, t, **params) is the per-component diagonal noise coefficient g (so dX_k = f_k dt + g_k dW_k).

5. Optionally mark integer / loop-structural parameters in _structural_params (baked in at lowering time, like the ODE family).

Example

Geometric Brownian motion dX = μX dt + σX dW::

class GeometricBrownianMotion(StochasticSystem):
    params = {"mu": 0.1, "sigma": 0.3}
    dim = 1
    variables = ("x",)

    @staticmethod
    def _drift(y, t, mu, sigma):
        return [mu * y(0)]

    @staticmethod
    def _diffusion(y, t, mu, sigma):
        return [sigma * y(0)]

gbm = GeometricBrownianMotion()
traj = gbm.integrate(final_time=1.0, dt=0.01, ic=[1.0], seed=0)

Notes

Integration uses a fixed step (the step is the noise scale √dt). method selects "euler_maruyama" (default, order 0.5) or "milstein" (order 1.0). Pass seed for a reproducible noise realisation; the resolved seed is recorded in the trajectory's meta.

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

def __init__(
    self,
    params: dict[str, Any] | None = None,
    ic: Any | None = None,
    dim: int | None = None,
    field_shape: tuple[int, ...] | None = None,
) -> None:
    """Initialise a system from its class defaults plus instance overrides.

    Parameters
    ----------
    params : dict, optional
        Per-instance parameter overrides.  Every key must already exist in
        the class-level :attr:`params` defaults; unknown keys raise.
    ic : array-like, optional
        Initial conditions.  Stored on ``self.ic`` (as a ``float`` array) and
        used by :meth:`resolve_ic` when no explicit ``ic`` is later supplied.
    dim : int, optional
        State-space dimension override for variable-dimension systems.  Falls
        back to the class-level :attr:`dim` when omitted.
    field_shape : tuple of int, optional
        Spatial grid shape override for a spatially-extended system (see
        :attr:`_field_shape`).  Falls back to the class-level value.

    Raises
    ------
    InvalidParameterError
        If ``params`` contains a key that is not a declared parameter.
    """
    # Build ParamSet from class defaults + constructor overrides
    defaults = dict(type(self).params)
    if params:
        unknown = set(params) - set(defaults)
        if unknown:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: unknown parameter(s) "
                f"{sorted(unknown)}. Declared: {sorted(defaults)}"
            )
        defaults.update(params)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    object.__setattr__(self, "dim", resolved_dim)

    # field_shape (spatially-extended systems): constructor arg > class
    # attribute.  Set via object.__setattr__ so a custom-N instance overrides
    # the class default without tripping the ClassVar assignment rule (mirrors
    # ``dim``).  ``_provenance`` reads the instance value onto ``traj.meta``.
    resolved_field_shape = field_shape if field_shape is not None else type(self)._field_shape
    object.__setattr__(self, "_field_shape", resolved_field_shape)

    # Initial conditions
    ic_arr = np.asarray(ic, dtype=float) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)

    # Metadata store: computed properties (Lyapunov, etc.) accumulate here
    # with history — repeated runs append instead of overwriting.
    object.__setattr__(self, "meta", MetaStore())

is_discrete property

is_discrete: bool

SDEs are continuous-time systems.

reinit

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

(Re)start the incremental stepper from state u at time t.

PARAMETER	DESCRIPTION
u	Initial state (falls back to self.ic, then random). TYPE: array - like DEFAULT: None
t	Start time (default 0.0). TYPE: float DEFAULT: None
params	Parameter overrides applied (in place) before restarting. TYPE: dict DEFAULT: None
method	"euler_maruyama" (default) or "milstein". TYPE: str DEFAULT: None
seed	Seed for the noise stream (random if omitted) — set it for a reproducible path. TYPE: int DEFAULT: None
dt	Default step size for :meth:step (default 0.01). TYPE: float DEFAULT: None

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

def reinit(
    self,
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    method: str | None = None,
    seed: int | None = None,
    dt: float | None = None,
) -> None:
    """
    (Re)start the incremental stepper from state ``u`` at time ``t``.

    Parameters
    ----------
    u : array-like, optional
        Initial state (falls back to ``self.ic``, then random).
    t : float, optional
        Start time (default 0.0).
    params : dict, optional
        Parameter overrides applied (in place) before restarting.
    method : str, optional
        ``"euler_maruyama"`` (default) or ``"milstein"``.
    seed : int, optional
        Seed for the noise stream (random if omitted) — set it for a
        reproducible path.
    dt : float, optional
        Default step size for :meth:`step` (default ``0.01``).
    """
    if params:
        for k, v in params.items():
            self.params[k] = v
    t0 = float(t) if t is not None else 0.0
    ic_arr = self.resolve_ic(u)
    canon = self._resolve_method(method)
    base_seed = _resolve_seed(seed)
    step_dt = float(dt) if dt is not None else type(self)._default_step_dt

    problem = self._problem(ic=ic_arr, t0=t0, method=canon)
    self._stepper = {
        "method": canon,
        "drift": problem.drift,
        "diffusion": problem.diffusion,
        "p": problem.params_vec(),
        "rng": _SplitMix64(base_seed),
        "seed": base_seed,
        "dt": step_dt,
    }
    self._state_now = ic_arr.copy()
    self._t_now = t0

step

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

Advance by dt (default 0.01) and return the new state.

The first call performs an implicit :meth:reinit. Each call draws a fresh diagonal Wiener increment from the stepper's seeded stream, so repeated stepping traces one reproducible sample path.

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

def step(self, n_or_dt: float | None = None) -> np.ndarray:
    """
    Advance by ``dt`` (default ``0.01``) and return the new state.

    The first call performs an implicit :meth:`reinit`. Each call draws a
    fresh diagonal Wiener increment from the stepper's seeded stream, so
    repeated stepping traces one reproducible sample path.
    """
    if self._stepper is None:
        self.reinit()
    assert self._stepper is not None
    s = self._stepper
    dt = float(n_or_dt) if n_or_dt is not None else s["dt"]
    u = np.asarray(self._state_now, dtype=float)
    dw = _wiener(s["rng"], dt, u.size)
    u = _sde_step(s["method"], s["drift"], s["diffusion"], u, s["p"], self._t_now, dw, dt)
    if not np.all(np.isfinite(u)):
        raise ConvergenceError(
            f"{type(self).__name__}: SDE diverged at t={self._t_now + dt:.6g} during step()."
        )
    self._t_now += dt
    self._state_now = u.copy()
    return u

state

state() -> ndarray

Return a copy of the current state (implicit reinit if cold).

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

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

set_state

set_state(u: Any) -> None

Overwrite the current state without changing the current time.

Unlike a DDE (whose state is a whole history function), an SDE's state is a single Markovian point, so set_state is well-defined.

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

def set_state(self, u: Any) -> None:
    """Overwrite the current state without changing the current time.

    Unlike a DDE (whose state is a whole history function), an SDE's state is
    a single Markovian point, so ``set_state`` is well-defined.
    """
    u_arr = np.asarray(u, dtype=float).reshape(self.dim)
    if self._stepper is None:
        self.reinit(u_arr)
    else:
        self._state_now = u_arr.copy()

time

time() -> float

Return the current stepper time.

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

def time(self) -> float:
    """Return the current stepper time."""
    return self._t_now

trajectory

trajectory(
    final_time: float = 100.0,
    *,
    dt: float = 0.02,
    transient: float = 0.0,
    **kwargs: Any,
) -> Trajectory

Protocol-uniform trajectory: integrate plus optional transient drop.

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

def trajectory(
    self,
    final_time: float = 100.0,
    *,
    dt: float = 0.02,
    transient: float = 0.0,
    **kwargs: Any,
) -> Trajectory:
    """Protocol-uniform trajectory: ``integrate`` plus optional transient drop."""
    traj = self.integrate(final_time=transient + final_time, dt=dt, **kwargs)
    return traj.after(transient) if transient > 0 else traj

run

run(
    final_time: float = 100.0,
    dt: float = 0.02,
    **kwargs: Any,
) -> Trajectory

Produce a trajectory — the one canonical verb for every family.

run is the unified trajectory producer: it answers the same call for flows, maps, DDEs and SDEs, dispatching on :attr:is_discrete. For a stochastic system (this family) it integrates the SDE, so run is a thin alias of :meth:integrate and forwards every keyword to it unchanged.

PARAMETER	DESCRIPTION
final_time	End of the integration window. Default 100.0. TYPE: float DEFAULT: 100.0
dt	Fixed step size and output sampling interval (for an SDE the step is the noise scale). TYPE: float DEFAULT: 0.02
**kwargs	Forwarded verbatim to :meth:integrate (t0, ic, method, seed, backend). TYPE: Any DEFAULT: {}

RETURNS	DESCRIPTION
Trajectory	Identical to :meth:integrate — run adds no behaviour.

See Also

integrate : The family-specific spelling (a permanent alias of run).

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

def run(
    self,
    final_time: float = 100.0,
    dt: float = 0.02,
    **kwargs: Any,
) -> Trajectory:
    """
    Produce a trajectory — the one canonical verb for every family.

    ``run`` is the unified trajectory producer: it answers the same call for
    flows, maps, DDEs and SDEs, dispatching on :attr:`is_discrete`.  For a
    stochastic system (this family) it integrates the SDE, so ``run`` is a
    thin alias of :meth:`integrate` and forwards every keyword to it
    unchanged.

    Parameters
    ----------
    final_time : float
        End of the integration window. Default 100.0.
    dt : float
        Fixed step size *and* output sampling interval (for an SDE the step
        is the noise scale).
    **kwargs
        Forwarded verbatim to :meth:`integrate` (``t0``, ``ic``, ``method``,
        ``seed``, ``backend``).

    Returns
    -------
    Trajectory
        Identical to :meth:`integrate` — ``run`` adds no behaviour.

    See Also
    --------
    integrate : The family-specific spelling (a permanent alias of ``run``).
    """
    return self.integrate(final_time=final_time, dt=dt, **kwargs)

integrate

integrate(
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    method: str | None = None,
    seed: int | None = None,
    backend: str | None = None,
) -> Trajectory

Integrate the SDE and return a :class:~tsdynamics.families.Trajectory.

PARAMETER	DESCRIPTION
final_time	End of the integration window. TYPE: float DEFAULT: 100.0
dt	Fixed step size and output sampling interval — for an SDE the step is the noise scale (each increment is drawn ~ N(0, dt)), so the output grid and the discretisation share one dt. TYPE: float DEFAULT: 0.02
t0	Start time (the IC is the state at t0). TYPE: float DEFAULT: 0.0
ic	Initial state. Falls back to self.ic, then U[0, 1)^dim. TYPE: array - like DEFAULT: None
method	"euler_maruyama" (default, order 0.5) or "milstein" (order 1.0). TYPE: str DEFAULT: None
seed	Seed for the noise realisation (random if omitted). The resolved seed is recorded in traj.meta["seed"] so a run can be reproduced. .. note:: integrate(seed=s) draws its noise from the raw seed s, whereas a single-row batch ensemble([ic], seed=s)[0] draws from the per-index stream seed_for(s, 0) (the parallel-equals-serial contract of :meth:ensemble). These two derived seeds differ, so the two calls trace different sample paths for the same s — by design, not a bug. To reproduce one ensemble trajectory standalone, integrate with that index's derived seed. TYPE: int DEFAULT: None
backend	Defaults to _default_backend ("interp"). "interp" / "jit" dispatch the two-tape SDE call to the compiled Rust engine (:mod:tsdynamics._rust) via the dedicated run.sde_integrate_dense seam. The engine path reproduces the reference to floating-point tolerance under a fixed seed (see the module docstring on the Box–Muller ULP) and raises :class:~tsdynamics.engine.run.EngineNotAvailableError if the extension is not built. TYPE: str DEFAULT: None

RETURNS	DESCRIPTION
Trajectory	Supports tuple-unpacking: t, y = sys.integrate(...).

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

def integrate(
    self,
    final_time: float = 100.0,
    dt: float = 0.02,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    method: str | None = None,
    seed: int | None = None,
    backend: str | None = None,
) -> Trajectory:
    """
    Integrate the SDE and return a :class:`~tsdynamics.families.Trajectory`.

    Parameters
    ----------
    final_time : float, default 100.0
        End of the integration window.
    dt : float, default 0.02
        Fixed step size *and* output sampling interval — for an SDE the step
        is the noise scale (each increment is drawn ``~ N(0, dt)``), so the
        output grid and the discretisation share one ``dt``.
    t0 : float, default 0.0
        Start time (the IC is the state at ``t0``).
    ic : array-like, optional
        Initial state. Falls back to ``self.ic``, then ``U[0, 1)^dim``.
    method : str, optional
        ``"euler_maruyama"`` (default, order 0.5) or ``"milstein"``
        (order 1.0).
    seed : int, optional
        Seed for the noise realisation (random if omitted). The resolved seed
        is recorded in ``traj.meta["seed"]`` so a run can be reproduced.

        .. note::
           ``integrate(seed=s)`` draws its noise from the **raw** seed ``s``,
           whereas a single-row batch ``ensemble([ic], seed=s)[0]`` draws from
           the **per-index** stream ``seed_for(s, 0)`` (the
           parallel-equals-serial contract of :meth:`ensemble`).  These two
           derived seeds differ, so the two calls trace **different** sample
           paths for the same ``s`` — by design, not a bug.  To reproduce one
           ``ensemble`` trajectory standalone, integrate with that index's
           derived seed.
    backend : str, optional
        Defaults to ``_default_backend`` (``"interp"``).  ``"interp"`` /
        ``"jit"`` dispatch the two-tape
        SDE call to the compiled Rust engine (:mod:`tsdynamics._rust`) via the
        dedicated ``run.sde_integrate_dense`` seam.  The engine path reproduces
        the reference to floating-point tolerance under a fixed seed (see the
        module docstring on the Box–Muller ULP) and raises
        :class:`~tsdynamics.engine.run.EngineNotAvailableError` if the
        extension is not built.

    Returns
    -------
    Trajectory
        Supports tuple-unpacking: ``t, y = sys.integrate(...)``.
    """
    from tsdynamics.engine import run

    backend = backend if backend is not None else self._default_backend
    canon = self._resolve_method(method)
    base_seed = _resolve_seed(seed)
    ic_arr = self.resolve_ic(ic)
    problem = self._problem(ic=ic_arr, t0=t0, method=canon)
    t_eval = make_output_grid(t0, final_time, dt)

    backend_canon = run.resolve_backend(backend)
    if backend_canon == "reference":
        rng = _SplitMix64(base_seed)
        y = self._run_reference(problem, t_eval, dt, canon, rng)
        engine = "reference"
    else:
        y = run.sde_integrate_dense(
            problem, t_eval, dt=dt, method=canon, seed=base_seed, backend=backend_canon
        )
        engine = "rust"

    return Trajectory(
        t=t_eval,
        y=y,
        system=self,
        meta=self._provenance(
            family="sde",
            engine=engine,
            backend=backend_canon,
            method=canon,
            dt=dt,
            t0=t0,
            seed=base_seed,
            ic=ic_arr.copy(),
        ),
    )

ensemble

ensemble(
    ics: Any,
    *,
    final_time: float = 100.0,
    dt: float = 0.02,
    t0: float = 0.0,
    method: str | None = None,
    seed: int | None = None,
    backend: str | None = None,
) -> ndarray

Integrate a batch of initial conditions and return their final states.

ics is (n, dim); each row is integrated from t0 to final_time and its final state returned as a row of the (n, dim) result. Trajectory i draws its noise from a stream seeded with seed_for(seed, i) — depending only on the index — so the batch is reproducible and matches the compiled engine's per-index seeding (the parallel-equals-serial contract). A diverging trajectory yields a row of NaN rather than aborting the batch.

.. note:: Because every trajectory is seeded by index, ensemble([ic], seed=s) draws from seed_for(s, 0), which differs from the raw seed s used by integrate(seed=s). The two calls therefore trace different sample paths for the same s — by design (the index seeding is what makes a batch reproducible and parallel-safe).

PARAMETER	DESCRIPTION
ics	The batch of initial conditions. TYPE: (array - like, shape(n, dim))
final_time	As in :meth:integrate (seed is the ensemble's base seed). TYPE: float DEFAULT: 100.0
dt	As in :meth:integrate (seed is the ensemble's base seed). TYPE: float DEFAULT: 100.0
t0	As in :meth:integrate (seed is the ensemble's base seed). TYPE: float DEFAULT: 100.0
method	As in :meth:integrate (seed is the ensemble's base seed). TYPE: float DEFAULT: 100.0
seed	As in :meth:integrate (seed is the ensemble's base seed). TYPE: float DEFAULT: 100.0
backend	Defaults to _default_backend ("interp"). "interp" / "jit" fan the batch out on the compiled engine's rayon pool; "reference" is a pure-Python loop. All seed each trajectory by index, so the final states match across backends to floating-point tolerance. TYPE: str DEFAULT: None

RETURNS	DESCRIPTION
(ndarray, shape(n, dim))	Final states (rows of NaN for diverged trajectories).

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

def ensemble(
    self,
    ics: Any,
    *,
    final_time: float = 100.0,
    dt: float = 0.02,
    t0: float = 0.0,
    method: str | None = None,
    seed: int | None = None,
    backend: str | None = None,
) -> np.ndarray:
    """
    Integrate a batch of initial conditions and return their final states.

    ``ics`` is ``(n, dim)``; each row is integrated from ``t0`` to
    ``final_time`` and its final state returned as a row of the ``(n, dim)``
    result. Trajectory ``i`` draws its noise from a stream seeded with
    ``seed_for(seed, i)`` — depending only on the index — so the batch is
    **reproducible** and matches the compiled engine's per-index seeding
    (the parallel-equals-serial contract). A diverging trajectory yields a
    row of ``NaN`` rather than aborting the batch.

    .. note::
       Because every trajectory is seeded by index, ``ensemble([ic], seed=s)``
       draws from ``seed_for(s, 0)``, which differs from the **raw** seed ``s``
       used by ``integrate(seed=s)``.  The two calls therefore trace different
       sample paths for the same ``s`` — by design (the index seeding is what
       makes a batch reproducible and parallel-safe).

    Parameters
    ----------
    ics : array-like, shape (n, dim)
        The batch of initial conditions.
    final_time, dt, t0, method, seed
        As in :meth:`integrate` (``seed`` is the ensemble's base seed).
    backend : str, optional
        Defaults to ``_default_backend`` (``"interp"``).  ``"interp"`` /
        ``"jit"`` fan the batch out on the compiled engine's rayon pool;
        ``"reference"`` is a pure-Python loop.  All seed each trajectory by index, so the final
        states match across backends to floating-point tolerance.

    Returns
    -------
    ndarray, shape (n, dim)
        Final states (rows of ``NaN`` for diverged trajectories).
    """
    from tsdynamics.engine import run

    backend = backend if backend is not None else self._default_backend
    canon = self._resolve_method(method)
    base_seed = _resolve_seed(seed)
    ics = np.ascontiguousarray(ics, dtype=np.float64)
    if ics.ndim != 2 or ics.shape[1] != self.dim:
        raise InvalidParameterError(f"ics must be (n, {self.dim}); got shape {ics.shape}")

    backend_canon = run.resolve_backend(backend)
    if backend_canon != "reference":
        problem = self._problem(ic=ics[0], t0=t0, method=canon)
        return run.sde_ensemble_final(
            problem,
            ics,
            t0=t0,
            t1=float(final_time),
            dt=dt,
            method=canon,
            seed=base_seed,
            backend=backend_canon,
        )

    problem = self._problem(ic=ics[0], t0=t0, method=canon)
    drift, diffusion = problem.drift, problem.diffusion
    p = problem.params_vec()
    out = np.empty_like(ics)
    for i, ic in enumerate(ics):
        rng = _SplitMix64(_seed_for(base_seed, i))
        try:
            u, _ = _advance_to(
                canon,
                drift,
                diffusion,
                np.asarray(ic, dtype=float),
                p,
                t0,
                float(final_time),
                dt,
                rng,
                type(self).__name__,
            )
            out[i] = u
        except RuntimeError:
            out[i] = np.nan
    return cast(np.ndarray, out)

System

Bases: Protocol

Structural type for steppable dynamical systems.

is_discrete property

is_discrete: bool

True for iterated maps (and map-like wrappers such as Poincaré maps).

step

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

Advance the system and return the new state.

The argument is the number of iterations for a discrete map and the time increment dt for a continuous flow (each family supplies a sensible default when None). Calling step on a fresh system performs an implicit :meth:reinit first.

RETURNS	DESCRIPTION
ndarray	The state after advancing.

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

def step(self, n_or_dt: float | int | None = None) -> np.ndarray:
    """Advance the system and return the new state.

    The argument is the number of iterations for a discrete map and the time
    increment ``dt`` for a continuous flow (each family supplies a sensible
    default when ``None``).  Calling ``step`` on a fresh system performs an
    implicit :meth:`reinit` first.

    Returns
    -------
    ndarray
        The state after advancing.
    """
    ...

state

state() -> ndarray

Return a copy of the current state vector.

Calling state on a fresh system performs an implicit :meth:reinit first. The returned array is a copy, so mutating it never disturbs the live stepper.

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

def state(self) -> np.ndarray:
    """Return a copy of the current state vector.

    Calling ``state`` on a fresh system performs an implicit :meth:`reinit`
    first.  The returned array is a copy, so mutating it never disturbs the
    live stepper.
    """
    ...

set_state

set_state(u: Any) -> None

Overwrite the current state.

Not available for delay systems — a DDE's state is a whole history function, not a point, so its implementation raises NotImplementedError; use :meth:reinit to restart from a constant past instead.

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

def set_state(self, u: Any) -> None:
    """Overwrite the current state.

    Not available for delay systems — a DDE's state is a whole history
    function, not a point, so its implementation raises
    ``NotImplementedError``; use :meth:`reinit` to restart from a constant
    past instead.
    """
    ...

time

time() -> float

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

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

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

reinit

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

Restart the stepper from state u at time t.

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

def reinit(
    self,
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
) -> None:
    """Restart the stepper from state ``u`` at time ``t``."""
    ...

trajectory

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

Produce a trajectory on a uniform output grid (the alias of :meth:run).

.. note:: run is the canonical trajectory-producer verb (see the module docstring), but trajectory is the member the structural protocol requires: every family and wrapper implements it, whereas a few (WrappedSystem and the derived wrappers) expose only trajectory — not run — so requiring run here would make them fail isinstance(obj, System). Code written against System should call trajectory; code holding a concrete flow or map should prefer run.

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

def trajectory(self, *args: Any, **kwargs: Any) -> Trajectory:
    """Produce a trajectory on a uniform output grid (the alias of :meth:`run`).

    .. note::
        ``run`` is the canonical trajectory-producer verb (see the module
        docstring), but ``trajectory`` is the member the *structural*
        protocol requires: every family and wrapper implements it, whereas a
        few (``WrappedSystem`` and the derived wrappers) expose only
        ``trajectory`` — not ``run`` — so requiring ``run`` here would make
        them fail ``isinstance(obj, System)``.  Code written against
        ``System`` should call ``trajectory``; code holding a concrete flow
        or map should prefer ``run``.
    """
    ...