Skip to content

Reference

Base classes

The shared machinery (ParamSet, 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.

run is the one verb that produces data on every one of them; reinit / step / state / time are the stepping protocol underneath it.

ParamSet

ParamSet(data: dict[str, Any])

Bases: dict[str, Any]

Ordered, fixed-key parameter container (a dict you cannot re-key).

Keys are frozen at construction — 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 structural mutation (del, pop, clear, ...).

Examples:

>>> p = ParamSet({"sigma": 10.0, "rho": 28.0})
>>> p
{'sigma': 10.0, 'rho': 28.0}
>>> isinstance(p, dict)
True
>>> p.as_tuple()
(10.0, 28.0)
>>> p["sigma"] = 12.0
>>> p.sigma
12.0
>>> p["sigmaa"] = 1.0
Traceback (most recent call last):
    ...
KeyError: "Unknown parameter 'sigmaa'. Declared params: ['sigma', 'rho']"
Source code in src/tsdynamics/families/_params.py
def __init__(self, data: dict[str, Any]) -> None:
    super().__init__(data)

pop

pop(*args: Any, **kwargs: Any) -> Any

Refuse — a parameter set has fixed keys.

Source code in src/tsdynamics/families/_params.py
def pop(self, *args: Any, **kwargs: Any) -> Any:
    """Refuse — a parameter set has fixed keys."""
    from tsdynamics.errors import InvalidInputError

    raise InvalidInputError("Parameters are fixed-key — cannot pop.")

popitem

popitem() -> tuple[str, Any]

Refuse — a parameter set has fixed keys.

Source code in src/tsdynamics/families/_params.py
def popitem(self) -> tuple[str, Any]:
    """Refuse — a parameter set has fixed keys."""
    from tsdynamics.errors import InvalidInputError

    raise InvalidInputError("Parameters are fixed-key — cannot popitem.")

clear

clear() -> None

Refuse — a parameter set has fixed keys.

Source code in src/tsdynamics/families/_params.py
def clear(self) -> None:
    """Refuse — a parameter set has fixed keys."""
    from tsdynamics.errors import InvalidInputError

    raise InvalidInputError("Parameters are fixed-key — cannot clear.")

fromkeys classmethod

fromkeys(*args: Any, **kwargs: Any) -> Any

Refuse — a parameter set's keys come from the system that declares them.

The inherited dict.fromkeys was already uncallable in every spelling: it builds the new mapping by calling cls(), and ParamSet.__init__ requires data, so the class method died on a raw TypeError: ParamSet.__init__() missing 1 required positional argument: 'data' — an error about a constructor the caller never wrote. Its three fixed-key siblings (pop / popitem / clear) all answer by name, so this one does too.

Source code in src/tsdynamics/families/_params.py
@classmethod
def fromkeys(cls, *args: Any, **kwargs: Any) -> Any:
    """Refuse — a parameter set's keys come from the system that declares them.

    The inherited ``dict.fromkeys`` was already **uncallable in every
    spelling**: it builds the new mapping by calling ``cls()``, and
    ``ParamSet.__init__`` requires ``data``, so the class method died on a
    raw ``TypeError: ParamSet.__init__() missing 1 required positional
    argument: 'data'`` — an error about a constructor the caller never
    wrote.  Its three fixed-key siblings (``pop`` / ``popitem`` /
    ``clear``) all answer by name, so this one does too.
    """
    from tsdynamics.errors import InvalidInputError, remedy

    raise InvalidInputError(
        "Parameters are fixed-key — cannot build a ParamSet from bare keys. "
        "The keys are the ones the system declares, and the values are the "
        "point."
        + remedy(
            'ParamSet({"sigma": 10.0, "rho": 28.0})',
            "system.with_params(sigma=12.0)",
        )
    )

setdefault

setdefault(key: str, default: Any = None) -> Any

Return self[key]; refuse to insert an undeclared key.

Source code in src/tsdynamics/families/_params.py
def setdefault(self, key: str, default: Any = None, /) -> Any:
    """Return ``self[key]``; refuse to *insert* an undeclared key."""
    if key not in self:
        self._reject(key)
    return self[key]

update

update(*args: Any, **kwargs: Any) -> None

Update values in place; refuse any key that is not declared.

Source code in src/tsdynamics/families/_params.py
def update(self, *args: Any, **kwargs: Any) -> None:
    """Update values in place; refuse any key that is not declared."""
    incoming: dict[str, Any] = dict(*args, **kwargs)
    unknown = [k for k in incoming if k not in self]
    if unknown:
        self._reject(unknown[0])
    for k, v in incoming.items():
        dict.__setitem__(self, k, v)

as_tuple

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

Return parameter values as a tuple (declaration order).

This ordering is the tape contract: the engine's lowering hands the values to the kernel positionally, so re-ordering them silently swaps parameters. dict preserves insertion order, which is why a plain dict subclass is a safe home for it.

Source code in src/tsdynamics/families/_params.py
def as_tuple(self) -> tuple[Any, ...]:
    """Return parameter values as a tuple (declaration order).

    This ordering is the **tape contract**: the engine's lowering hands the
    values to the kernel positionally, so re-ordering them silently swaps
    parameters.  ``dict`` preserves insertion order, which is why a plain
    ``dict`` subclass is a safe home for it.
    """
    return tuple(self.values())

as_dict

as_dict() -> dict[str, Any]

Return a shallow copy as a plain dict.

Source code in src/tsdynamics/families/_params.py
def as_dict(self) -> dict[str, Any]:
    """Return a shallow copy as a plain ``dict``."""
    return dict(self)

copy

copy() -> ParamSet

Return an independent :class:ParamSet with the same keys and values.

Source code in src/tsdynamics/families/_params.py
def copy(self) -> ParamSet:
    """Return an independent :class:`ParamSet` with the same keys and values."""
    return ParamSet(dict(self))

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.

Source code in src/tsdynamics/families/_params.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.
    """
    import hashlib
    import json

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

Trajectory

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

The result of integrating or iterating a dynamical system.

A trajectory is a container of samples: len(traj) is the number of time samples and iterating it yields that many (t_i, y_i) pairs.

.. versionchanged:: 6.0 __iter__ used to be a tuple-unpacking convenience yielding the two columns (t, y), which put it in direct contradiction with __len__ (len(traj) said n_steps, len(list(traj)) said 2) and broke the container contract for any code that sizes an iterable and then walks it. Iteration now yields samples; the column unpacking moved to the explicit :meth:unpack::

    t, y = traj.unpack()          # was: t, y = traj
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.run(final_time=100)
>>> traj.dim
3
>>> traj["x"]            # named component (via the class's ``variables``)
array([...])
>>> traj.after(20.0)     # drop transient
Trajectory of Lorenz  (n: ..., dim: 3, var: x, y, z)   t ∈ [20, 100]
    measure it:  ts.analysis.find(traj)   ·   draw it:  ts.plot(traj)
>>> t, y = traj.unpack()          # the two columns
>>> for t_i, y_i in traj:         # ... or walk the samples
...     pass

Build a trajectory from time and state arrays.

PARAMETER DESCRIPTION
t

Time points (or step indices).

TYPE: (array_like, shape(T))

y

State at each time point. A 1-D array is read as a single component and stored as (T, 1).

TYPE: (array_like, shape(T) or (T, dim))

system

The system that produced this trajectory. Measured data has no system, so it defaults to None — every consumer in the library tolerates that, and it is what lets a user's own arrays reach the plotting and analysis layers::

traj = ts.Trajectory(t, y)
ts.plot(traj)

TYPE: SystemBase DEFAULT: None

meta

Provenance (system name, params, solver, dt, tolerances, ic).

TYPE: dict DEFAULT: None

variables

One name per state component — what :meth:component, traj["x"], :meth:to_frame's columns and every plot axis label use. This is the door measured data comes in by, which has no system to read declared names off::

traj = ts.Trajectory(t, x, variables=("voltage",))
traj["voltage"]

.. versionadded:: 6.0 It used to be reachable only as the undocumented meta={"variables": (...)}, and naming it raised a bare TypeError — on the one constructor the package docstring points data users at, so everything a measured signal plotted was labelled y0.

TYPE: sequence of str DEFAULT: None

RAISES DESCRIPTION
InvalidInputError

If variables does not name every component exactly once.

Source code in src/tsdynamics/data/trajectory.py
def __init__(
    self,
    t: np.ndarray,
    y: np.ndarray,
    system: Any = None,
    meta: dict[str, Any] | None = None,
    *,
    variables: Sequence[str] | None = None,
) -> None:
    """Build a trajectory from time and state arrays.

    Parameters
    ----------
    t : array_like, shape (T,)
        Time points (or step indices).
    y : array_like, shape (T,) or (T, dim)
        State at each time point.  A 1-D array is read as a single
        component and stored as ``(T, 1)``.
    system : SystemBase, optional
        The system that produced this trajectory.  **Measured data has no
        system**, so it defaults to ``None`` — every consumer in the
        library tolerates that, and it is what lets a user's own arrays
        reach the plotting and analysis layers::

            traj = ts.Trajectory(t, y)
            ts.plot(traj)

    meta : dict, optional
        Provenance (system name, params, solver, ``dt``, tolerances, ic).
    variables : sequence of str, optional
        One name per state component — what :meth:`component`, ``traj["x"]``,
        :meth:`to_frame`'s columns and every plot axis label use.  This is
        the door **measured** data comes in by, which has no system to read
        declared names off::

            traj = ts.Trajectory(t, x, variables=("voltage",))
            traj["voltage"]

        .. versionadded:: 6.0
            It used to be reachable only as the undocumented
            ``meta={"variables": (...)}``, and naming it raised a bare
            ``TypeError`` — on the one constructor the package docstring
            points data users at, so everything a measured signal plotted
            was labelled ``y0``.

    Raises
    ------
    InvalidInputError
        If ``variables`` does not name every component exactly once.
    """
    self.t = np.asarray(t)
    y_arr = np.asarray(y)
    self.y = y_arr[:, None] if y_arr.ndim == 1 else y_arr
    self.system = system
    self.meta = dict(meta) if meta else {}
    if variables is not None:
        self.meta["variables"] = _validated_variables(variables, int(self.y.shape[1]))
    self._kdtree: cKDTree | None = None
    self._unbounded: Unbounded | None = None
    self._unbounded_checked = False

unbounded property

unbounded: Unbounded | None

What this orbit escaped to, or None if it stayed bounded.

A run that blows up without reaching the engine's hard 1e150 guard comes back as an ordinary, finite, entirely meaningless trajectory — measured, Chua().run(ic=[500, 0, 0]) returns max|y| = 1.6e9, no exception, isnan all False. Every downstream answer taken from it then looks exactly like an honest one. This is the flag that tells the two apart: it is the line the repr prints, and what :class:~tsdynamics.derived.poincare.PoincareSection carries through.

Computed on first read and cached, so a trajectory nobody asks about pays nothing.

RETURNS DESCRIPTION
Unbounded or None

None when the orbit is bounded. Otherwise the record whose str is the warning line — see :class:tsdynamics._utils.escape.Unbounded.

Examples:

>>> import tsdynamics as ts
>>> lor = ts.systems.Lorenz()
>>> lor.run(final_time=5.0, ic=[1.0, 1.0, 1.0]).unbounded is None
True

shape property

shape: tuple[int, int]

(n_samples, dim) — the shape of :attr:y.

dt property

dt: float | None

The sampling interval, derived from the t axis.

The axis is the one truth, and None when there is no single answer — a non-uniform axis (a Poincaré section) or an axis too short to have a step. meta["dt"] used to record what the run asked for and slicing carried it verbatim, so a decimated trajectory reported the undecimated step and every per-unit-time estimator taken from meta was off by the decimation factor with no exception. Since v6 round 8 selecting rows re-derives meta["dt"] from the new axis (:data:_ROW_DERIVED_META), so the two agree — but this property is still the reading the library itself uses, because it cannot be absent.

Examples:

>>> import tsdynamics as ts
>>> tr = ts.systems.Lorenz().run(final_time=1.0, dt=0.01, ic=[1.0, 1.0, 1.0])
>>> round(tr.dt, 10)
0.01
>>> round(tr[::5].dt, 10)          # the decimated step, not the run's
0.05

variables property

variables: tuple[str, ...]

The name of every component — never None.

Resolution order: the names this trajectory carries in meta["variables"] (written by a column selection, so a sub-trajectory names its own columns), then the producing system's variables, then the generated y0 … y{dim-1}.

.. versionchanged:: 6.0 This used to be a pass-through to system.variables returning None for measured data — so a bare-array trajectory printed y0/y1 column headings in :meth:to_frame and on its plot, then refused traj["y0"]; and a column selection reported the inner system's full tuple, which made traj["x", "z"]["y"] silently return z. The names now come from one place and are always as long as the state block is wide.

dim property

dim: int

State-space dimension.

n_steps property

n_steps: int

Number of time steps.

unpack

unpack() -> tuple[ndarray, ndarray]

Return the two column arrays (t, y).

The explicit spelling of what t, y = traj used to do implicitly (see :meth:__iter__)::

t, y = system.run(final_time=100).unpack()
RETURNS DESCRIPTION
tuple of (ndarray, ndarray)

t of shape (T,) and y of shape (T, dim) — the live arrays, not copies.

Source code in src/tsdynamics/data/trajectory.py
def unpack(self) -> tuple[np.ndarray, np.ndarray]:
    """Return the two column arrays ``(t, y)``.

    The explicit spelling of what ``t, y = traj`` used to do implicitly (see
    :meth:`__iter__`)::

        t, y = system.run(final_time=100).unpack()

    Returns
    -------
    tuple of (ndarray, ndarray)
        ``t`` of shape ``(T,)`` and ``y`` of shape ``(T, dim)`` — the live
        arrays, not copies.
    """
    return self.t, self.y

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 — keep the samples at or after t0.

PARAMETER DESCRIPTION
t0

A cut on the recorded time axis, read in that axis's own unit: time units for a flow, the iteration index for a map, the crossing time for a section. (It is a cut, not a horizon: run(t0=...) says where the integration starts.)

TYPE: float

RETURNS DESCRIPTION
Trajectory

The tail, sharing this trajectory's system and provenance.

Source code in src/tsdynamics/data/trajectory.py
def after(self, t0: float) -> Trajectory:
    """
    Drop the initial transient — keep the samples at or after ``t0``.

    Parameters
    ----------
    t0 : float
        A **cut on the recorded time axis**, read in that axis's own unit:
        time units for a flow, the iteration index for a map, the crossing
        time for a section.  (It is a cut, not a horizon: ``run(t0=...)``
        says where the integration *starts*.)

    Returns
    -------
    Trajectory
        The tail, sharing this trajectory's system and provenance.
    """
    mask = self.t >= t0
    return self._rows(self.t[mask], self.y[mask])

to_frame

to_frame() -> Any

Return a :class:pandas.DataFrame view of the trajectory.

One row per time sample: the time vector becomes the index (named "t") and each state component a column, named from the system's variables when it declares them and y0…y{dim-1} otherwise. The provenance travels along on frame.attrs["meta"].

pandas is a soft dependency imported lazily here, so importing tsdynamics never pulls it in; a missing install raises an :class:ImportError pointing at pip install pandas (and at the traj.t / traj.y arrays, which need nothing extra).

RETURNS DESCRIPTION
DataFrame

Shape (n_steps, dim).

RAISES DESCRIPTION
ImportError

If :mod:pandas is not installed.

Source code in src/tsdynamics/data/trajectory.py
def to_frame(self) -> Any:
    """Return a :class:`pandas.DataFrame` view of the trajectory.

    One row per time sample: the time vector becomes the index (named
    ``"t"``) and each state component a column, named from the system's
    ``variables`` when it declares them and ``y0…y{dim-1}`` otherwise.  The
    provenance travels along on ``frame.attrs["meta"]``.

    ``pandas`` is a **soft** dependency imported lazily here, so importing
    tsdynamics never pulls it in; a missing install raises an
    :class:`ImportError` pointing at ``pip install pandas`` (and at the
    ``traj.t`` / ``traj.y`` arrays, which need nothing extra).

    Returns
    -------
    pandas.DataFrame
        Shape ``(n_steps, dim)``.

    Raises
    ------
    ImportError
        If :mod:`pandas` is not installed.
    """
    pd = self._require_pandas()
    names = self.variables
    frame = pd.DataFrame(
        self.y,
        index=pd.Index(self.t, name="t"),
        columns=list(names),
    )
    frame.attrs["meta"] = dict(self.meta)
    return frame

minmax

minmax() -> MinMax

Return the per-component extent — (lo, hi), each of shape (dim,).

Unpacks like the bare tuple it used to be (lo, hi = traj.minmax()) and names its halves (traj.minmax().hi), so a reader of the call site does not have to remember which one comes first.

Source code in src/tsdynamics/data/trajectory.py
def minmax(self) -> MinMax:
    """Return the per-component extent — ``(lo, hi)``, each of shape ``(dim,)``.

    Unpacks like the bare tuple it used to be (``lo, hi = traj.minmax()``)
    and names its halves (``traj.minmax().hi``), so a reader of the call
    site does not have to remember which one comes first.
    """
    return MinMax(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) -> Neighbors

Return the trajectory states nearest 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). A single point gives a result without the leading query axis; m points keep it.

TYPE: (array - like, shape(dim) or (m, dim))

k

Number of neighbours per query point. The k axis is always present, whatever k is.

TYPE: int DEFAULT: 1

RETURNS DESCRIPTION
Neighbors

.state — the neighbouring states, shape (k, dim) (or (m, k, dim)); .distance and .index — shape (k,) (or (m, k)). It unpacks as (distance, index), so d, i = traj.neighbors(q) reads as it always did.

.. versionchanged:: 6.0

This handed back scipy.spatial.cKDTree.query's raw (distances, indices) — two unlabelled arrays whose rank changed with k (k=1 gave two scalars, k=2 two arrays, so i[0] worked for one and raised for the other), and which never contained the thing a caller is usually after: the neighbouring states.

Examples:

>>> import numpy as np, tsdynamics as ts
>>> tr = ts.systems.Lorenz().run(final_time=1.0, dt=0.1, ic=[1.0, 1.0, 1.0])
>>> near = tr.neighbors([1.0, 1.0, 1.0], k=2)
>>> near.state.shape, near.distance.shape, near.index.shape
((2, 3), (2,), (2,))
Source code in src/tsdynamics/data/trajectory.py
def neighbors(self, q: Any, k: int = 1) -> Neighbors:
    """Return the trajectory **states** nearest 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).  A single point gives a result without the leading
        query axis; ``m`` points keep it.
    k : int, default 1
        Number of neighbours per query point.  The ``k`` axis is **always**
        present, whatever ``k`` is.

    Returns
    -------
    Neighbors
        ``.state`` — the neighbouring states, shape ``(k, dim)`` (or
        ``(m, k, dim)``); ``.distance`` and ``.index`` — shape ``(k,)`` (or
        ``(m, k)``).  It unpacks as ``(distance, index)``, so
        ``d, i = traj.neighbors(q)`` reads as it always did.

    .. versionchanged:: 6.0
        This handed back ``scipy.spatial.cKDTree.query``'s raw
        ``(distances, indices)`` — two unlabelled arrays whose *rank* changed
        with ``k`` (``k=1`` gave two scalars, ``k=2`` two arrays, so
        ``i[0]`` worked for one and raised for the other), and which never
        contained the thing a caller is usually after: the neighbouring
        states.

    Examples
    --------
    >>> import numpy as np, tsdynamics as ts
    >>> tr = ts.systems.Lorenz().run(final_time=1.0, dt=0.1, ic=[1.0, 1.0, 1.0])
    >>> near = tr.neighbors([1.0, 1.0, 1.0], k=2)
    >>> near.state.shape, near.distance.shape, near.index.shape
    ((2, 3), (2,), (2,))
    """
    from scipy.spatial import cKDTree

    if self._kdtree is None:
        self._kdtree = cKDTree(self.y)
    query = np.asarray(q, dtype=float)
    raw_d, raw_i = self._kdtree.query(query, k=k)
    # scipy squeezes the ``k`` axis when ``k == 1``; restore it so the answer
    # has one shape for every ``k``.  The *query* axis follows the input's
    # rank, which is what a caller who passed one point means.
    distance = np.asarray(raw_d, dtype=float)
    index = np.asarray(raw_i)
    if k == 1:
        distance, index = distance[..., None], index[..., None]
    # A query for more neighbours than there are points: scipy reports the
    # one-past-the-end index with an infinite distance.  Report the state as
    # NaN rather than indexing out of the array.
    missing = index >= self.y.shape[0]
    state = np.asarray(self.y, dtype=float)[np.where(missing, 0, index)]
    if bool(missing.any()):
        state = state.copy()
        state[missing] = np.nan
    return Neighbors(distance=distance, index=index, state=state)

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,
    seed: int | None = None,
    **param_kwargs: Any,
)

Bases: DeriveMixin, 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

Every declared parameter is a plain constructor keyword — the most natural line a user will type::

lor = Lorenz(rho=30.0, ic=[1.0, 0.0, 0.0])

The params= dict spelling stays supported and means exactly the same thing (it is the way to pass a name computed at runtime)::

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 — whether it arrives as params={"rhoo": 30.0} or as rhoo=30.0 — so a typo fails loudly, naming the declared parameters, instead of being silently ignored. Giving the same parameter through both channels is also an error: there is deliberately no precedence rule to memorise.

The five keywords :meth:__init__ reserves for itself (params, ic, dim, field_shape, seed) cannot also be parameter names — a class declaring one is refused at definition time by :meth:__init_subclass__ rather than silently shadowing it.

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. An ic given here is explicit: it is never silently replaced by a random draw (a diverging one raises instead).

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

seed

Seed for the initial-condition draw used when neither ic nor system.info.default_ic supplies one. The draw runs on a private :class:numpy.random.Generator, so it is reproducible and never touches the global numpy.random stream. When omitted a fresh OS-entropy seed is drawn on first use and recorded on traj.meta["ic_seed"], so any run can be reproduced exactly by passing that value back as seed=.

TYPE: int DEFAULT: None

**param_kwargs

Per-parameter overrides given as plain keywords — Lorenz(sigma=12.0) is exactly Lorenz(params={"sigma": 12.0}). Each name must be a declared parameter of this class; an unknown one raises and names the valid options. A name given both here and inside params= is an error rather than a silent precedence rule.

TYPE: Any DEFAULT: {}

RAISES DESCRIPTION
InvalidParameterError

If params or **param_kwargs contains a key that is not a declared parameter, or if the same parameter is given twice (once in params= and once as a keyword).

Examples:

>>> Lorenz(sigma=12.0).sigma            # the natural spelling
12.0
>>> Lorenz(params={"sigma": 12.0}).sigma   # equivalent, still supported
12.0
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,
    seed: int | None = None,
    **param_kwargs: Any,
) -> 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.
        An ``ic`` given here is **explicit**: it is never silently replaced by
        a random draw (a diverging one raises instead).
    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.
    seed : int, optional
        Seed for the **initial-condition draw** used when neither ``ic`` nor
        ``system.info.default_ic`` supplies one.  The draw runs on a private
        :class:`numpy.random.Generator`, so it is reproducible *and* never
        touches the global ``numpy.random`` stream.  When omitted a fresh
        OS-entropy seed is drawn on first use and recorded on
        ``traj.meta["ic_seed"]``, so any run can be reproduced exactly by
        passing that value back as ``seed=``.
    **param_kwargs
        Per-parameter overrides given as plain keywords — ``Lorenz(sigma=12.0)``
        is exactly ``Lorenz(params={"sigma": 12.0})``.  Each name must be a
        declared parameter of this class; an unknown one raises and names the
        valid options.  A name given **both** here and inside ``params=`` is an
        error rather than a silent precedence rule.

    Raises
    ------
    InvalidParameterError
        If ``params`` or ``**param_kwargs`` contains a key that is not a
        declared parameter, or if the same parameter is given twice (once in
        ``params=`` and once as a keyword).

    Examples
    --------
    >>> Lorenz(sigma=12.0).sigma            # the natural spelling
    12.0
    >>> Lorenz(params={"sigma": 12.0}).sigma   # equivalent, still supported
    12.0
    """
    # Build ParamSet from class defaults + constructor overrides.  The two
    # override channels (``params=`` and free keywords) are merged here, and
    # both are validated against the *declared* parameter names.
    defaults = dict(type(self).params)
    _reject_state_as_params(type(self).__name__, params, defaults)
    overrides: dict[str, Any] = dict(params) if params else {}
    if param_kwargs:
        duplicated = sorted(set(param_kwargs) & set(overrides))
        if duplicated:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: parameter(s) {duplicated} given twice — "
                f"once in params= and once as a keyword. Pass each parameter "
                f"exactly once (there is deliberately no precedence rule)."
            )
        overrides.update(param_kwargs)
    if overrides:
        unknown = set(overrides) - set(defaults)
        if unknown:
            # A misspelt parameter is nearly always one character off a
            # declared one, so name the closest match and spell the whole
            # constructor call rather than only listing what was valid.
            import difflib

            from tsdynamics.errors import InvalidParameterError

            cls_name = type(self).__name__
            guesses = {
                bad: (
                    [n for n in sorted(defaults) if n.lower() == bad.lower()]
                    or difflib.get_close_matches(bad, sorted(defaults), n=1, cutoff=0.5)
                )
                for bad in sorted(unknown)
            }
            fixed = {bad: g[0] for bad, g in guesses.items() if g}
            did_you_mean = (
                " Did you mean " + ", ".join(f"{b!r} -> {g!r}" for b, g in fixed.items()) + "?"
                if fixed
                else ""
            )
            call = ", ".join(
                f"{fixed.get(bad, bad)}={overrides[bad]!r}" for bad in sorted(unknown)
            )
            raise InvalidParameterError(
                f"{cls_name}: unknown parameter(s) {sorted(unknown)}. "
                f"Declared: {sorted(defaults)}.{did_you_mean}"
                + remedy(
                    f"{cls_name}({call})"
                    if fixed
                    else f"{cls_name}({sorted(defaults)[0]}={defaults[sorted(defaults)[0]]!r})",
                    lead="Spell it as:" if fixed else "The declared parameters are set like:",
                )
            )
        defaults.update(overrides)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    if resolved_dim is None:
        # Every downstream consumer does ``int(self.dim)``, so an undeclared
        # dimension used to surface as NumPy's ``int() argument must be ...
        # not 'NoneType'`` from deep inside resolve_ic -- which names neither
        # the class nor the attribute.  Writing a system class is the first
        # thing a new user does, so refuse it here, where the fix is one line.
        raise InvalidInputError(
            f"{type(self).__name__} does not declare its state-space dimension, "
            f"so the library cannot tell how many components its state has. "
            f"`dim` is the number of equations `_equations` returns."
            + remedy(
                f"class {type(self).__name__}({type(self).__bases__[0].__name__}):",
                "    dim = 3                       # the number of state components",
                '    variables = ("x", "y", "z")   # optional, names them',
                lead="Declare it on the class:",
            )
            + remedy(
                f"{type(self).__name__}(dim=3)",
                lead="...or pass it for this one instance:",
            )
        )
    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_explicit`` records whether the CURRENT
    # ``self.ic`` was chosen by the user (constructor / an explicit ``ic=``)
    # or merely auto-resolved (a random draw).  A user-chosen IC must never
    # be silently swapped for a random one — see ``DiscreteMap.run``.
    # ``np.array(..., copy=True)``, not ``asarray``: a float64 array passed in
    # would otherwise be *shared* with the caller, so mutating either side
    # silently moved the other's initial condition.  It is also what made
    # ``with_params`` (which forwards ``ic=self.ic``) hand back a system
    # aliasing the source's ic array.
    ic_arr = np.array(ic, dtype=float, copy=True) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)
    object.__setattr__(self, "_ic_explicit", ic is not None)

    # The private initial-condition RNG (see ``resolve_ic``).  Only the seed
    # is persisted; the Generator itself is a transient cache.
    object.__setattr__(self, "_ic_seed", None if seed is None else int(seed))
    object.__setattr__(self, "_ic_rng", None)

family property

family: str

"ode" | "dde" | "map" | "sde" — what kind of dynamics this is.

Replaces is_discrete, which could not tell a delay system from a stochastic one, so every caller that needed the distinction had to sniff the class.

info property

info: SystemInfo

Everything true about this system, in one printable record.

Absorbs reference, doi, known_lyapunov, field_labels and default_ic off the tab surface — they are facts, and a fact belongs in the record you print::

>>> print(Rossler().info)               # doctest: +SKIP
Roessler — 3-D continuous flow                     tsdynamics...Rossler
  equations   dx/dt = -y - z
  ...

copy

copy() -> SystemBase

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

The copy has its own independent params, so mutating the clone's parameters never affects the original. dim and field_shape are forwarded, so a system that took them from the constructor copies.

RETURNS DESCRIPTION
SystemBase

A fresh instance of the same subclass with copied params and ic.

See Also

copy : copy.copy(system) — the same independence, but it keeps the recorded meta (and does not re-run __init__).

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``, so mutating the clone's
    parameters never affects the original.  ``dim`` and ``field_shape`` are
    forwarded, so a system that took them from the constructor copies.

    Returns
    -------
    SystemBase
        A fresh instance of the same subclass with copied ``params`` and ``ic``.

    See Also
    --------
    __copy__ : ``copy.copy(system)`` — the same independence, but it *keeps*
        the recorded ``meta`` (and does not re-run ``__init__``).
    """
    return self._carry_ic_provenance(
        type(self)(
            params=cast(ParamSet, self.params).as_dict(),
            ic=self.ic.copy() if self.ic is not None else None,
            dim=cast(int, self.dim),
            field_shape=self.__dict__.get("_field_shape"),
        )
    )

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).run(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).run(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}
    # ``dim`` and ``field_shape`` are forwarded because a system that takes
    # them from the CONSTRUCTOR (the shape the docs teach you to write) has
    # nothing to re-read them from — measured, ``DimCtor(dim=2).with_params(k=2)``
    # raised "does not declare its state-space dimension", which broke
    # continuation and orbit diagrams for exactly those systems.
    return self._carry_ic_provenance(
        type(self)(
            params=new_p,
            ic=self.ic,
            dim=cast(int, self.dim),
            field_shape=self.__dict__.get("_field_shape"),
        )
    )

ContinuousSystem

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

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. Its signature is _equations(y, t, **params), and y is a state ACCESSOR, not an array: component i is y(i) — a call, which is also what lets a :class:~tsdynamics.families.delay.DelaySystem write a delayed access as y(i, t - tau). So neither y[0] nor x, y, z = y works here (a DiscreteMap's _step does take a plain state vector — that is the one place the two families differ)::

    class MyLorenz(ts.ContinuousSystem): params = {"sigma": 10.0, "rho": 28.0, "beta": 8 / 3} variables = ("x", "y", "z") dim = 3

       @staticmethod
       def _equations(y, t, sigma, rho, beta):
           x, u, z = y(0), y(1), y(2)          # NOT `x, u, z = y`
           return [sigma * (u - x), x * (rho - z) - u, x * u - beta * z]
    
  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, then JIT-compiled on first use (both memoised, so neither cost repeats); 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 The numerical kernel this system integrates with when the caller names none — the class-level counterpart of run(solver=...). Default "rk45"; declare "bdf" on a reliably stiff system. Spellings are normalised, so "RK45" and "dopri5" resolve to the same kernel, and system.info prints the resolved name.

Examples:

>>> lor = Lorenz()
>>> traj = lor.run(final_time=100, dt=0.01)
>>> t, y = traj.unpack()   # the two columns
>>> lor.sigma = 15.0     # change param — zero recompile cost
>>> traj2 = lor.run(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,
    seed: int | None = None,
    **param_kwargs: Any,
) -> 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.
        An ``ic`` given here is **explicit**: it is never silently replaced by
        a random draw (a diverging one raises instead).
    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.
    seed : int, optional
        Seed for the **initial-condition draw** used when neither ``ic`` nor
        ``system.info.default_ic`` supplies one.  The draw runs on a private
        :class:`numpy.random.Generator`, so it is reproducible *and* never
        touches the global ``numpy.random`` stream.  When omitted a fresh
        OS-entropy seed is drawn on first use and recorded on
        ``traj.meta["ic_seed"]``, so any run can be reproduced exactly by
        passing that value back as ``seed=``.
    **param_kwargs
        Per-parameter overrides given as plain keywords — ``Lorenz(sigma=12.0)``
        is exactly ``Lorenz(params={"sigma": 12.0})``.  Each name must be a
        declared parameter of this class; an unknown one raises and names the
        valid options.  A name given **both** here and inside ``params=`` is an
        error rather than a silent precedence rule.

    Raises
    ------
    InvalidParameterError
        If ``params`` or ``**param_kwargs`` contains a key that is not a
        declared parameter, or if the same parameter is given twice (once in
        ``params=`` and once as a keyword).

    Examples
    --------
    >>> Lorenz(sigma=12.0).sigma            # the natural spelling
    12.0
    >>> Lorenz(params={"sigma": 12.0}).sigma   # equivalent, still supported
    12.0
    """
    # Build ParamSet from class defaults + constructor overrides.  The two
    # override channels (``params=`` and free keywords) are merged here, and
    # both are validated against the *declared* parameter names.
    defaults = dict(type(self).params)
    _reject_state_as_params(type(self).__name__, params, defaults)
    overrides: dict[str, Any] = dict(params) if params else {}
    if param_kwargs:
        duplicated = sorted(set(param_kwargs) & set(overrides))
        if duplicated:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: parameter(s) {duplicated} given twice — "
                f"once in params= and once as a keyword. Pass each parameter "
                f"exactly once (there is deliberately no precedence rule)."
            )
        overrides.update(param_kwargs)
    if overrides:
        unknown = set(overrides) - set(defaults)
        if unknown:
            # A misspelt parameter is nearly always one character off a
            # declared one, so name the closest match and spell the whole
            # constructor call rather than only listing what was valid.
            import difflib

            from tsdynamics.errors import InvalidParameterError

            cls_name = type(self).__name__
            guesses = {
                bad: (
                    [n for n in sorted(defaults) if n.lower() == bad.lower()]
                    or difflib.get_close_matches(bad, sorted(defaults), n=1, cutoff=0.5)
                )
                for bad in sorted(unknown)
            }
            fixed = {bad: g[0] for bad, g in guesses.items() if g}
            did_you_mean = (
                " Did you mean " + ", ".join(f"{b!r} -> {g!r}" for b, g in fixed.items()) + "?"
                if fixed
                else ""
            )
            call = ", ".join(
                f"{fixed.get(bad, bad)}={overrides[bad]!r}" for bad in sorted(unknown)
            )
            raise InvalidParameterError(
                f"{cls_name}: unknown parameter(s) {sorted(unknown)}. "
                f"Declared: {sorted(defaults)}.{did_you_mean}"
                + remedy(
                    f"{cls_name}({call})"
                    if fixed
                    else f"{cls_name}({sorted(defaults)[0]}={defaults[sorted(defaults)[0]]!r})",
                    lead="Spell it as:" if fixed else "The declared parameters are set like:",
                )
            )
        defaults.update(overrides)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    if resolved_dim is None:
        # Every downstream consumer does ``int(self.dim)``, so an undeclared
        # dimension used to surface as NumPy's ``int() argument must be ...
        # not 'NoneType'`` from deep inside resolve_ic -- which names neither
        # the class nor the attribute.  Writing a system class is the first
        # thing a new user does, so refuse it here, where the fix is one line.
        raise InvalidInputError(
            f"{type(self).__name__} does not declare its state-space dimension, "
            f"so the library cannot tell how many components its state has. "
            f"`dim` is the number of equations `_equations` returns."
            + remedy(
                f"class {type(self).__name__}({type(self).__bases__[0].__name__}):",
                "    dim = 3                       # the number of state components",
                '    variables = ("x", "y", "z")   # optional, names them',
                lead="Declare it on the class:",
            )
            + remedy(
                f"{type(self).__name__}(dim=3)",
                lead="...or pass it for this one instance:",
            )
        )
    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_explicit`` records whether the CURRENT
    # ``self.ic`` was chosen by the user (constructor / an explicit ``ic=``)
    # or merely auto-resolved (a random draw).  A user-chosen IC must never
    # be silently swapped for a random one — see ``DiscreteMap.run``.
    # ``np.array(..., copy=True)``, not ``asarray``: a float64 array passed in
    # would otherwise be *shared* with the caller, so mutating either side
    # silently moved the other's initial condition.  It is also what made
    # ``with_params`` (which forwards ``ic=self.ic``) hand back a system
    # aliasing the source's ic array.
    ic_arr = np.array(ic, dtype=float, copy=True) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)
    object.__setattr__(self, "_ic_explicit", ic is not None)

    # The private initial-condition RNG (see ``resolve_ic``).  Only the seed
    # is persisted; the Generator itself is a transient cache.
    object.__setattr__(self, "_ic_seed", None if seed is None else int(seed))
    object.__setattr__(self, "_ic_rng", None)

reinit

reinit(
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    solver: str | None = None,
    rtol: float = DEFAULT_RTOL,
    atol: float = DEFAULT_ATOL,
    max_step: float | None = None,
    backend: str | None = None,
    **unknown: Any,
) -> None

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

solver= is the numerical kernel, spelled the same way run spells it (method= selects an estimator in v6) — and it is refused here with the same message run gives, because a keyword that teaches at one door and raises a bare TypeError at the next has taught nothing.

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

solver

Stepper configuration, as in :meth:run. max_step is stored and applied to every subsequent :meth:step.

TYPE: str | None DEFAULT: None

rtol

Stepper configuration, as in :meth:run. max_step is stored and applied to every subsequent :meth:step.

TYPE: str | None DEFAULT: None

atol

Stepper configuration, as in :meth:run. max_step is stored and applied to every subsequent :meth:step.

TYPE: str | None DEFAULT: None

max_step

Stepper configuration, as in :meth:run. max_step is stored and applied to every subsequent :meth:step.

TYPE: str | None DEFAULT: None

backend

Stepper configuration, as in :meth:run. max_step is stored and applied to every subsequent :meth:step.

TYPE: str | None DEFAULT: None

Notes

A reinit that raises leaves the system exactly as it was — self.ic is rolled back rather than left holding the initial condition that failed.

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,
    solver: str | None = None,
    rtol: float = DEFAULT_RTOL,
    atol: float = DEFAULT_ATOL,
    max_step: float | None = None,
    backend: str | None = None,
    **unknown: Any,
) -> None:
    """
    (Re)start the incremental stepper from state ``u`` at time ``t``.

    ``solver=`` is the numerical kernel, spelled the same way ``run`` spells
    it (``method=`` selects an *estimator* in v6) — and it is refused here
    with the same message ``run`` gives, because a keyword that teaches at
    one door and raises a bare ``TypeError`` at the next has taught nothing.

    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.
    solver, rtol, atol, max_step, backend
        Stepper configuration, as in :meth:`run`.  ``max_step`` is
        stored and applied to every subsequent :meth:`step`.

    Notes
    -----
    A ``reinit`` that raises leaves the system exactly as it was — ``self.ic``
    is rolled back rather than left holding the initial condition that failed.
    """
    if params:
        for k, v in params.items():
            self.params[k] = v
    # ``resolve_ic`` commits the resolved IC to ``self.ic`` before any of the
    # work below happens, and plenty below can raise (an unknown backend, an
    # unresolvable method, a tape that will not lower, an unavailable engine).
    # Without the guard a failed ``reinit`` latched the offending IC onto the
    # instance and every later, unrelated call silently started from it.
    reject_unknown_run_keywords(
        self,
        unknown,
        family="ode",
        accepted=("t", "params", "solver", "rtol", "atol", "max_step", "backend"),
        verb="reinit",
    )
    with self._ic_rollback():
        self._reinit_resolved(
            u,
            t=t,
            method=solver,
            rtol=rtol,
            atol=atol,
            max_step=max_step,
            backend=backend,
        )

step

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

Advance the system by n_or_dt and return the new state.

n_or_dt is a time increment, in the same unit as final_time; omitting it advances by _default_step_dt (0.01, the number every continuous family uses). A map counts iterations instead — that is the only thing the word means differently anywhere.

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:run 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 lyapunov_spectrum. 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 ``n_or_dt`` and return the new state.

    ``n_or_dt`` is a **time increment**, in the same unit as ``final_time``;
    omitting it advances by ``_default_step_dt`` (``0.01``, the number every
    continuous family uses).  A map counts **iterations** instead — that is
    the only thing the word means differently anywhere.

    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:`run` 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 ``lyapunov_spectrum``.  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,
            max_step=math.inf if self._step_max_step is None else float(self._step_max_step),
            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__,
        max_step=self._step_max_step,
    )
    # 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

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

This is the only route to the Jacobian as algebra, which is why it keeps a slot on a 19-name tab surface (CONTRACT.md §11, T2). :meth:jacobian evaluates ∂f/∂y to numbers at one state; this hands back the expressions, with the control parameters still symbols — so you can read the coupling structure off the matrix, take a further derivative, substitute a parameter to see where a term vanishes, or export the equations to another tool. Nothing else in the library exposes them: the engine consumes the lowered tape, and every analysis consumes numbers.

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

jacobian : the same matrix evaluated numerically at (u, t).

Examples:

>>> J = Lorenz().jacobian_sym()
>>> J[0]                                 # the row d(dx/dt)/dy
[-sigma, sigma, 0]
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).

    **This is the only route to the Jacobian as algebra**, which is why it
    keeps a slot on a 19-name tab surface (``CONTRACT.md`` §11, T2).
    :meth:`jacobian` evaluates ∂f/∂y to *numbers* at one state; this hands
    back the expressions, with the control parameters still symbols — so you
    can read the coupling structure off the matrix, take a further
    derivative, substitute a parameter to see where a term vanishes, or
    export the equations to another tool.  Nothing else in the library
    exposes them: the engine consumes the lowered tape, and every analysis
    consumes numbers.

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

    See Also
    --------
    jacobian : the same matrix evaluated numerically at ``(u, t)``.

    Examples
    --------
    >>> J = Lorenz().jacobian_sym()          # doctest: +SKIP
    >>> J[0]                                 # the row d(dx/dt)/dy
    [-sigma, sigma, 0]
    """
    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)

rhs

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

Evaluate the right-hand side :math:f(u, t) — the vector field itself.

The function :meth:jacobian is the derivative of. It had no public door until v6: system.jacobian(u, t) was public while the field it differentiates was reachable only as system._rhs_numeric(), which is a strange split — and a costly one, because three of the plotting layer's transforms (vector_field, flow_speed, streamlines) are this function, so a user had no sanctioned way to check the picture against the maths. Two independent readers and a blind tester reached for the private name; that is the definition of a missing door.

PARAMETER DESCRIPTION
u

State at which to evaluate the field.

TYPE: (array - like, shape(dim))

t

Time (it matters only for a non-autonomous system). Default 0.

TYPE: float DEFAULT: 0.0

RETURNS DESCRIPTION
(ndarray, shape(dim))

:math:du/dt at (u, t), with the system's current parameter values.

See Also

jacobian : The matrix of this function's partial derivatives.

Examples:

>>> import numpy as np, tsdynamics as ts
>>> lor = ts.systems.Lorenz()
>>> np.round(lor.rhs([1.0, 1.0, 1.0]), 6)
array([  0.      ,  26.      , -1.666667])

The Jacobian is this field's derivative, and a finite difference of one proves the other::

>>> u = np.array([1.0, 2.0, 3.0])
>>> h, col = 1e-6, 1
>>> step = np.zeros(3); step[col] = h
>>> numeric = (lor.rhs(u + step) - lor.rhs(u - step)) / (2 * h)
>>> bool(np.allclose(numeric, lor.jacobian(u)[:, col], atol=1e-6))
True
Source code in src/tsdynamics/families/continuous.py
def rhs(self, u: Any, t: float = 0.0) -> np.ndarray:
    r"""Evaluate the right-hand side :math:`f(u, t)` — the vector field itself.

    The function :meth:`jacobian` is the derivative *of*.  It had no public
    door until v6: ``system.jacobian(u, t)`` was public while the field it
    differentiates was reachable only as ``system._rhs_numeric()``, which is
    a strange split — and a costly one, because three of the plotting
    layer's transforms (``vector_field``, ``flow_speed``, ``streamlines``)
    **are** this function, so a user had no sanctioned way to check the
    picture against the maths.  Two independent readers and a blind tester
    reached for the private name; that is the definition of a missing door.

    Parameters
    ----------
    u : array-like, shape (dim,)
        State at which to evaluate the field.
    t : float, optional
        Time (it matters only for a non-autonomous system).  Default ``0``.

    Returns
    -------
    ndarray, shape (dim,)
        :math:`du/dt` at ``(u, t)``, with the system's **current** parameter
        values.

    See Also
    --------
    jacobian : The matrix of this function's partial derivatives.

    Examples
    --------
    >>> import numpy as np, tsdynamics as ts
    >>> lor = ts.systems.Lorenz()
    >>> np.round(lor.rhs([1.0, 1.0, 1.0]), 6)
    array([  0.      ,  26.      , -1.666667])

    The Jacobian is this field's derivative, and a finite difference of one
    proves the other::

    >>> u = np.array([1.0, 2.0, 3.0])
    >>> h, col = 1e-6, 1
    >>> step = np.zeros(3); step[col] = h
    >>> numeric = (lor.rhs(u + step) - lor.rhs(u - step)) / (2 * h)
    >>> bool(np.allclose(numeric, lor.jacobian(u)[:, col], atol=1e-6))
    True
    """
    dim = cast(int, self.dim)
    rhs_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(rhs_fn(arg), dtype=float).reshape(dim)

run

run(
    final_time: float = 100.0,
    dt: float | None = None,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    transient: float = 0.0,
    solver: str | None = None,
    rtol: float = DEFAULT_RTOL,
    atol: float = DEFAULT_ATOL,
    max_step: float | None = None,
    backend: str | None = None,
    seed: int | None = None,
    events: Any = None,
    **solver_options: Any,
) -> Trajectory

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

run is the trajectory verb — one word on flows, maps, delay and stochastic systems. integrate and trajectory were two more names for this method and are gone in v6.

PARAMETER DESCRIPTION
final_time

End of the integration window, in time units — the horizon word for a flow. (A map counts iterations instead and takes steps.) Default 100.0.

TYPE: float DEFAULT: 100.0

dt

Output sampling interval, in time units — the spacing of the returned grid. None (the default) means this family's _default_dt, which system.info prints under defaults. It does not set the integration accuracy: the internal stepper is adaptive and controlled by rtol/atol; interior samples are produced by the kernel's own continuous extension. Use rtol / atol to control accuracy and max_step to bound the internal step.

.. versionchanged:: 6.0 Before v6 the stepper was forced to land on every output sample, so a fine dt silently bought extra accuracy and a coarse one silently lost it. It no longer does. Kernels without a native continuous extension (everything but rk45 / tsit5 / dop853) still land on every sample.

TYPE: float DEFAULT: None

t0

Where the integration starts, in time units. Default 0.0. Allows warm restarts from a non-zero time (the IC is interpreted as the state at t0). Not to be confused with traj.after(t0), which cuts an already-recorded axis.

TYPE: float DEFAULT: 0.0

ic

Initial state at t0 — dim numbers, one per state component. Falls back to self.ic, then a U[0, 1)^dim draw. Passing it here does not change self.ic: one call, one run.

TYPE: array - like DEFAULT: None

solver

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 — the accuracy knob. Default :data:~tsdynamics._utils.tolerances.DEFAULT_RTOL / :data:~tsdynamics._utils.tolerances.DEFAULT_ATOL (1e-9 / 1e-12).

.. versionchanged:: 6.0 Tightened from 1e-6 / 1e-9. dt is now purely an output grid, so rtol is the only accuracy knob and the default had to carry the accuracy the old forced landing supplied for free. Measured median 1459x more accurate for 1.74x the cost; pass rtol=1e-6 for the pre-v6 trade.

TYPE: float DEFAULT: DEFAULT_RTOL

atol

Solver tolerances — the accuracy knob. Default :data:~tsdynamics._utils.tolerances.DEFAULT_RTOL / :data:~tsdynamics._utils.tolerances.DEFAULT_ATOL (1e-9 / 1e-12).

.. versionchanged:: 6.0 Tightened from 1e-6 / 1e-9. dt is now purely an output grid, so rtol is the only accuracy knob and the default had to carry the accuracy the old forced landing supplied for free. Measured median 1459x more accurate for 1.74x the cost; pass rtol=1e-6 for the pre-v6 trade.

TYPE: float DEFAULT: DEFAULT_RTOL

max_step

Upper bound on any single internal solver step, in time units. None (default) means no ceiling — the adaptive controller chooses freely from rtol/atol. Use it to stop an adaptive kernel stepping over a narrow feature (a thin resonance, a fast transient) in an otherwise smooth region, or to bound the detection resolution of an event march. Note this is a step size; the engine's max_steps is a step count and is not exposed here. A max_step far below the tolerance-driven natural step multiplies cost with no accuracy benefit and can trip :class:~tsdynamics.errors.StepBudgetError.

TYPE: float DEFAULT: None

backend

Where the ODE is integrated. Defaults to _default_backend ("jit"). The first two go through the shared engine seam (:func:tsdynamics._engine.run.integrate).

  • "jit" (default) — the Cranelift JIT: the lowered tape compiled to native code, with a process-wide compiled-evaluator cache, so the compile is paid once per distinct system rather than on every call.
  • "interp" — the SSA-tape interpreter. Same answers as the JIT, bit-for-bit; marginally faster on very small tapes, and the way to avoid the one-off compile.
  • "reference" — the dependency-light pure-Python/SciPy oracle. Not for production use: it is the independent cross-check the engine is validated against, and the wheel-free fallback.

.. versionchanged:: 6.0 The default moved from "interp" to "jit". Before v6 the JIT recompiled the whole tape on every FFI call, which made it slower than the interpreter for short runs; the v6 compiled-evaluator cache removed that per-call compile.

TYPE: ('jit', 'interp', 'reference') DEFAULT: "jit"

seed

Seed for the random initial-condition draw — the same meaning it has on :meth:DiscreteMap.run, on :meth:DelaySystem.run and on the constructor, so seed= reads identically on every family. (An SDE has a second source of randomness, so there seed= seeds the noise path as well as the draw.) It only bites when a draw happens: an explicit ic, an ic already resolved onto the system, and a class-level default_ic all take priority. The resolved seed is recorded on traj.meta["ic_seed"].

.. versionadded:: 6.0

TYPE: int DEFAULT: None

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

transient

Leading stretch of the run to discard, in time units (the same unit as final_time). The window is extended to transient + final_time and everything before t0 + transient is dropped, so the returned trajectory still spans final_time. One word on every family, in that family's own horizon unit: time here, iterations on a map (see :meth:DiscreteMap.run).

.. versionadded:: 6.0 Only the retired trajectory verb used to accept it, so a user who found run could not discard a transient at all.

TYPE: float DEFAULT: 0.0

RETURNS DESCRIPTION
Trajectory

A container of samples: len(traj) time points, iterating yields (t_i, y_i) pairs; traj.unpack() gives the two column arrays (t, y).

Source code in src/tsdynamics/families/continuous.py
def run(
    self,
    final_time: float = 100.0,
    dt: float | None = None,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    transient: float = 0.0,
    solver: str | None = None,
    rtol: float = DEFAULT_RTOL,
    atol: float = DEFAULT_ATOL,
    max_step: float | None = None,
    backend: str | None = None,
    seed: int | None = None,
    events: Any = None,
    **solver_options: Any,
) -> Trajectory:
    """
    Integrate the flow and return a :class:`~tsdynamics.families.Trajectory`.

    ``run`` is **the** trajectory verb — one word on flows, maps, delay and
    stochastic systems.  ``integrate`` and ``trajectory`` were two more names
    for this method and are gone in v6.

    Parameters
    ----------
    final_time : float
        End of the integration window, **in time units** — the horizon word
        for a flow.  (A map counts iterations instead and takes ``steps``.)
        Default 100.0.
    dt : float, optional
        Output sampling interval, **in time units** — the spacing of the
        returned grid.  ``None`` (the default) means this family's
        ``_default_dt``, which ``system.info`` prints under ``defaults``.  **It
        does not set the integration accuracy:** the internal stepper is
        adaptive and controlled by ``rtol``/``atol``; interior samples are
        produced by the kernel's own continuous extension.  Use ``rtol`` /
        ``atol`` to control accuracy and ``max_step`` to bound the internal
        step.

        .. versionchanged:: 6.0
            Before v6 the stepper was forced to land on every output sample,
            so a fine ``dt`` silently bought extra accuracy and a coarse one
            silently lost it.  It no longer does.  Kernels without a native
            continuous extension (everything but ``rk45`` / ``tsit5`` /
            ``dop853``) still land on every sample.
    t0 : float
        Where the integration **starts**, in time units. Default 0.0.
        Allows warm restarts from a non-zero time (the IC is interpreted as
        the state at ``t0``).  Not to be confused with ``traj.after(t0)``,
        which cuts an already-recorded axis.
    ic : array-like, optional
        Initial state at ``t0`` — ``dim`` numbers, one per state component.
        Falls back to ``self.ic``, then a ``U[0, 1)^dim`` draw.  **Passing
        it here does not change ``self.ic``**: one call, one run.
    solver : 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 — the accuracy knob.  Default
        :data:`~tsdynamics._utils.tolerances.DEFAULT_RTOL` /
        :data:`~tsdynamics._utils.tolerances.DEFAULT_ATOL` (``1e-9`` /
        ``1e-12``).

        .. versionchanged:: 6.0
            Tightened from ``1e-6`` / ``1e-9``.  ``dt`` is now purely an
            output grid, so ``rtol`` is the *only* accuracy knob and the
            default had to carry the accuracy the old forced landing supplied
            for free.  Measured median 1459x more accurate for 1.74x the
            cost; pass ``rtol=1e-6`` for the pre-v6 trade.
    max_step : float, optional
        Upper bound on any single internal solver step, in time units.
        ``None`` (default) means no ceiling — the adaptive controller chooses
        freely from ``rtol``/``atol``.  Use it to stop an adaptive kernel
        stepping over a narrow feature (a thin resonance, a fast transient)
        in an otherwise smooth region, or to bound the detection resolution
        of an event march.  Note this is a step *size*; the engine's
        ``max_steps`` is a step *count* and is not exposed here.  A
        ``max_step`` far below the tolerance-driven natural step multiplies
        cost with no accuracy benefit and can trip
        :class:`~tsdynamics.errors.StepBudgetError`.
    backend : {"jit", "interp", "reference"}, optional
        Where the ODE is integrated.  Defaults to ``_default_backend``
        (``"jit"``).  The first two go through the shared engine seam
        (:func:`tsdynamics._engine.run.integrate`).

        - ``"jit"`` (default) — the **Cranelift JIT**: the lowered tape
          compiled to native code, with a process-wide compiled-evaluator
          cache, so the compile is paid once per distinct system rather than
          on every call.
        - ``"interp"`` — the **SSA-tape interpreter**.  Same answers as the
          JIT, bit-for-bit; marginally faster on very small tapes, and the
          way to avoid the one-off compile.
        - ``"reference"`` — the dependency-light pure-Python/SciPy oracle.
          Not for production use: it is the independent cross-check the
          engine is validated against, and the wheel-free fallback.

        .. versionchanged:: 6.0
            The default moved from ``"interp"`` to ``"jit"``.  Before v6 the
            JIT recompiled the whole tape on every FFI call, which made it
            slower than the interpreter for short runs; the v6
            compiled-evaluator cache removed that per-call compile.
    seed : int, optional
        Seed for the **random initial-condition draw** — the same meaning it
        has on :meth:`DiscreteMap.run`, on :meth:`DelaySystem.run` and on the
        constructor, so ``seed=`` reads identically on every family.  (An SDE
        has a second source of randomness, so there ``seed=`` seeds the noise
        path as well as the draw.)  It only bites when a draw happens: an
        explicit ``ic``, an ``ic`` already resolved onto the system, and a
        class-level ``default_ic`` all take priority.  The resolved seed is
        recorded on ``traj.meta["ic_seed"]``.

        .. versionadded:: 6.0
    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"]``.
    transient : float, optional
        Leading stretch of the run to discard, in **time units** (the same
        unit as ``final_time``).  The window is extended to
        ``transient + final_time`` and everything before ``t0 + transient``
        is dropped, so the returned trajectory still spans ``final_time``.
        One word on every family, in that family's **own horizon unit**:
        time here, **iterations** on a map (see :meth:`DiscreteMap.run`).

        .. versionadded:: 6.0
            Only the retired ``trajectory`` verb used to accept it, so a
            user who found ``run`` could not discard a transient at all.

    Returns
    -------
    Trajectory
        A container of samples: ``len(traj)`` time points, iterating
        yields ``(t_i, y_i)`` pairs; ``traj.unpack()`` gives the two
        column arrays ``(t, y)``.
    """
    reject_unknown_run_keywords(self, solver_options, family="ode", accepted=_ODE_RUN_KEYWORDS)
    dt = self._default_dt if dt is None else dt
    transient = resolve_transient(transient, discrete=False)
    if transient > 0.0:
        traj = self.run(
            final_time=final_time + transient,
            dt=dt,
            t0=t0,
            ic=ic,
            solver=solver,
            rtol=rtol,
            atol=atol,
            max_step=max_step,
            backend=backend,
            seed=seed,
            events=events,
        )
        return traj.after(t0 + transient)
    if events is not None:
        return self._run_events(
            final_time=final_time,
            dt=dt,
            events=events,
            t0=t0,
            ic=ic,
            method=solver,
            rtol=rtol,
            atol=atol,
            max_step=max_step,
            backend=backend,
            seed=seed,
        )
    backend = backend if backend is not None else self._default_backend
    return self._dispatch(
        backend=backend,
        seed=seed,
        final_time=final_time,
        dt=dt,
        t0=t0,
        ic=ic,
        method=solver or self._default_method,
        rtol=rtol,
        atol=atol,
        max_step=max_step,
    )

DelaySystem

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

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, then JIT-compiled on first use (both memoised, so neither cost repeats). 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.

Tolerances

The DDE family keeps its own, looser default — :data:~tsdynamics._utils.tolerances.DDE_RTOL / :data:~tsdynamics._utils.tolerances.DDE_ATOL (1e-3 / 1e-3) — rather than the ODE :data:~tsdynamics._utils.tolerances.DEFAULT_RTOL / :data:~tsdynamics._utils.tolerances.DEFAULT_ATOL. The reason is not that the solver struggles when tightened (it does not: all six built-in DDEs complete at rtol=1e-12/atol=1e-15). It is that the tolerance is largely inert here: the method of steps lands on every output sample, so dt already bounds the internal step below the natural error. Measured over all six built-in DDEs at the default dt=0.02 to T=10, five return a bit-identical final state at rtol=1e-3 and at rtol=1e-9; the sixth (IkedaDelay, the only one whose step is tolerance-bound) costs 1.4x for a 6x accuracy gain. Tighten it explicitly if you need that.

.. versionchanged:: 6.0 The ODE default tightened to 1e-9/1e-12 to compensate for native dense output. DDEs never had dense output — the march always landed on every sample — so they lost nothing and this default is unchanged.

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:

>>> import tsdynamics as ts
>>> mg = MackeyGlass()
>>> hist = lambda s: [1.0 + 0.1 * np.sin(0.2 * s)]
>>> traj = mg.run(final_time=500, history=hist)
>>> exps = ts.analysis.lyapunov_spectrum(mg, k=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,
    seed: int | None = None,
    **param_kwargs: Any,
) -> 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.
        An ``ic`` given here is **explicit**: it is never silently replaced by
        a random draw (a diverging one raises instead).
    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.
    seed : int, optional
        Seed for the **initial-condition draw** used when neither ``ic`` nor
        ``system.info.default_ic`` supplies one.  The draw runs on a private
        :class:`numpy.random.Generator`, so it is reproducible *and* never
        touches the global ``numpy.random`` stream.  When omitted a fresh
        OS-entropy seed is drawn on first use and recorded on
        ``traj.meta["ic_seed"]``, so any run can be reproduced exactly by
        passing that value back as ``seed=``.
    **param_kwargs
        Per-parameter overrides given as plain keywords — ``Lorenz(sigma=12.0)``
        is exactly ``Lorenz(params={"sigma": 12.0})``.  Each name must be a
        declared parameter of this class; an unknown one raises and names the
        valid options.  A name given **both** here and inside ``params=`` is an
        error rather than a silent precedence rule.

    Raises
    ------
    InvalidParameterError
        If ``params`` or ``**param_kwargs`` contains a key that is not a
        declared parameter, or if the same parameter is given twice (once in
        ``params=`` and once as a keyword).

    Examples
    --------
    >>> Lorenz(sigma=12.0).sigma            # the natural spelling
    12.0
    >>> Lorenz(params={"sigma": 12.0}).sigma   # equivalent, still supported
    12.0
    """
    # Build ParamSet from class defaults + constructor overrides.  The two
    # override channels (``params=`` and free keywords) are merged here, and
    # both are validated against the *declared* parameter names.
    defaults = dict(type(self).params)
    _reject_state_as_params(type(self).__name__, params, defaults)
    overrides: dict[str, Any] = dict(params) if params else {}
    if param_kwargs:
        duplicated = sorted(set(param_kwargs) & set(overrides))
        if duplicated:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: parameter(s) {duplicated} given twice — "
                f"once in params= and once as a keyword. Pass each parameter "
                f"exactly once (there is deliberately no precedence rule)."
            )
        overrides.update(param_kwargs)
    if overrides:
        unknown = set(overrides) - set(defaults)
        if unknown:
            # A misspelt parameter is nearly always one character off a
            # declared one, so name the closest match and spell the whole
            # constructor call rather than only listing what was valid.
            import difflib

            from tsdynamics.errors import InvalidParameterError

            cls_name = type(self).__name__
            guesses = {
                bad: (
                    [n for n in sorted(defaults) if n.lower() == bad.lower()]
                    or difflib.get_close_matches(bad, sorted(defaults), n=1, cutoff=0.5)
                )
                for bad in sorted(unknown)
            }
            fixed = {bad: g[0] for bad, g in guesses.items() if g}
            did_you_mean = (
                " Did you mean " + ", ".join(f"{b!r} -> {g!r}" for b, g in fixed.items()) + "?"
                if fixed
                else ""
            )
            call = ", ".join(
                f"{fixed.get(bad, bad)}={overrides[bad]!r}" for bad in sorted(unknown)
            )
            raise InvalidParameterError(
                f"{cls_name}: unknown parameter(s) {sorted(unknown)}. "
                f"Declared: {sorted(defaults)}.{did_you_mean}"
                + remedy(
                    f"{cls_name}({call})"
                    if fixed
                    else f"{cls_name}({sorted(defaults)[0]}={defaults[sorted(defaults)[0]]!r})",
                    lead="Spell it as:" if fixed else "The declared parameters are set like:",
                )
            )
        defaults.update(overrides)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    if resolved_dim is None:
        # Every downstream consumer does ``int(self.dim)``, so an undeclared
        # dimension used to surface as NumPy's ``int() argument must be ...
        # not 'NoneType'`` from deep inside resolve_ic -- which names neither
        # the class nor the attribute.  Writing a system class is the first
        # thing a new user does, so refuse it here, where the fix is one line.
        raise InvalidInputError(
            f"{type(self).__name__} does not declare its state-space dimension, "
            f"so the library cannot tell how many components its state has. "
            f"`dim` is the number of equations `_equations` returns."
            + remedy(
                f"class {type(self).__name__}({type(self).__bases__[0].__name__}):",
                "    dim = 3                       # the number of state components",
                '    variables = ("x", "y", "z")   # optional, names them',
                lead="Declare it on the class:",
            )
            + remedy(
                f"{type(self).__name__}(dim=3)",
                lead="...or pass it for this one instance:",
            )
        )
    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_explicit`` records whether the CURRENT
    # ``self.ic`` was chosen by the user (constructor / an explicit ``ic=``)
    # or merely auto-resolved (a random draw).  A user-chosen IC must never
    # be silently swapped for a random one — see ``DiscreteMap.run``.
    # ``np.array(..., copy=True)``, not ``asarray``: a float64 array passed in
    # would otherwise be *shared* with the caller, so mutating either side
    # silently moved the other's initial condition.  It is also what made
    # ``with_params`` (which forwards ``ic=self.ic``) hand back a system
    # aliasing the source's ic array.
    ic_arr = np.array(ic, dtype=float, copy=True) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)
    object.__setattr__(self, "_ic_explicit", ic is not None)

    # The private initial-condition RNG (see ``resolve_ic``).  Only the seed
    # is persisted; the Generator itself is a transient cache.
    object.__setattr__(self, "_ic_seed", None if seed is None else int(seed))
    object.__setattr__(self, "_ic_rng", None)

reinit

reinit(
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    rtol: float | None = None,
    atol: float | None = None,
    backend: str | None = None,
    **unknown: 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, call :meth:run 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:run for a full trajectory.

PARAMETER DESCRIPTION
u

The constant past to restart from — dim numbers, held for every s ≤ 0.

TYPE: array - like DEFAULT: None

t

Must be 0.0: a delay system's clock is pinned to its history window [-tau_max, 0].

TYPE: float DEFAULT: None

params

Parameter overrides applied to this system in place.

TYPE: dict DEFAULT: None

rtol

Tolerances for the per-step march. Default :data:~tsdynamics._utils.tolerances.DDE_RTOL / DDE_ATOL.

TYPE: float DEFAULT: None

atol

Tolerances for the per-step march. Default :data:~tsdynamics._utils.tolerances.DDE_RTOL / DDE_ATOL.

TYPE: float DEFAULT: None

backend

Which evaluator drives the per-step march — the same word :meth:run takes, accepted here for the same reason :meth:ContinuousSystem.reinit accepts it: a stepping loop must be able to choose its engine. "reference" is refused (there is no pure-Python delay integrator).

TYPE: ('jit', 'interp') DEFAULT: "jit"

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,
    backend: str | None = None,
    **unknown: 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, call :meth:`run` 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:`run` for a full
    trajectory.

    Parameters
    ----------
    u : array-like, optional
        The **constant past** to restart from — ``dim`` numbers, held for
        every ``s ≤ 0``.
    t : float, optional
        Must be ``0.0``: a delay system's clock is pinned to its history
        window ``[-tau_max, 0]``.
    params : dict, optional
        Parameter overrides applied to this system **in place**.
    rtol, atol : float, optional
        Tolerances for the per-step march.  Default
        :data:`~tsdynamics._utils.tolerances.DDE_RTOL` / ``DDE_ATOL``.
    backend : {"jit", "interp"}, optional
        Which evaluator drives the per-step march — the same word
        :meth:`run` takes, accepted here for the same reason
        :meth:`ContinuousSystem.reinit` accepts it: a stepping loop must be
        able to choose its engine.  ``"reference"`` is refused (there is no
        pure-Python delay integrator).
    """
    reject_unknown_run_keywords(
        self,
        unknown,
        family="dde",
        accepted=("t", "params", "rtol", "atol", "backend"),
        verb="reinit",
    )
    if backend is not None:
        from tsdynamics._engine.run import resolve_backend

        if resolve_backend(backend) == "reference":
            raise InvalidParameterError(
                f"{type(self).__name__}.reinit(backend='reference'): there is no "
                "pure-Python delay integrator, so a DDE has no reference path."
                + remedy(f"{type(self).__name__}().reinit(backend='interp')")
            )
        self._step_backend = backend
    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)."
        )
    # ``resolve_ic`` commits the resolved IC to ``self.ic`` before anything
    # else, so a malformed ``u`` must leave the object exactly as it was —
    # the same contract ``integrate`` gets from ``_dispatch``.
    with self._ic_rollback():
        past_ic = self._resolve_ic(u)
    self._past_ic = past_ic
    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 n_or_dt time units and return the new state (forward-only).

n_or_dt is a time increment here, as on every continuous family; omitting it advances by _default_step_dt (0.01, the same number :class:~tsdynamics.families.ContinuousSystem and :class:~tsdynamics.families.StochasticSystem use — it was 0.1, a difference no mathematics asked for). A map counts iterations instead.

Source code in src/tsdynamics/families/delay.py
def step(self, n_or_dt: float | None = None) -> np.ndarray:
    """Advance by ``n_or_dt`` time units and return the new state (forward-only).

    ``n_or_dt`` is a **time increment** here, as on every continuous family;
    omitting it advances by ``_default_step_dt`` (``0.01``, the same number
    :class:`~tsdynamics.families.ContinuousSystem` and
    :class:`~tsdynamics.families.StochasticSystem` use — it was ``0.1``, a
    difference no mathematics asked for).  A map counts iterations instead.
    """
    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.run(
        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,
        backend=self._step_backend,
    )
    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()

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

run

run(
    final_time: float = 100.0,
    dt: float | None = None,
    *,
    ic: Any | None = None,
    history: History = None,
    transient: float = 0.0,
    solver: str | None = None,
    rtol: float | None = None,
    atol: float | None = None,
    backend: str | None = None,
    seed: int | None = None,
    **solver_options: Any,
) -> Trajectory

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

run is the trajectory verb — one word on flows, maps, delay and stochastic systems. integrate and trajectory were two more names for this method and are gone in v6.

The signature is closed. Before v6 this method accepted — and silently dropped — max_step, t0, events and outright typos: the run completed, the number was wrong, and nothing said so.

PARAMETER DESCRIPTION
final_time

End of the integration window, in time units — the horizon word for a flow. Default 100.0.

TYPE: float DEFAULT: 100.0

dt

Output sampling interval, in time units. None (the default) means this family's _default_dt, which system.info prints under defaults. The method of steps lands on every sample, so dt also bounds the internal step here (see rtol/atol).

TYPE: float DEFAULT: None

ic

The constant past, used when history is None — dim numbers, one per state component, held for every s ≤ 0. Falls back to self.ic, then a random draw.

TYPE: array - like DEFAULT: None

history

The past, as a function of the lag: h(s) is called with s ≤ 0 (s = 0 is the start of the run, s = -tau one delay before it) and must return dim numbers — the state at that moment of the pre-history. It is never called with s > 0. None (the default) uses a constant past equal to ic::

mg.run(final_time=500.0, history=lambda s: [1.0 + 0.1 * np.sin(0.2 * s)])

A constant past sitting exactly on a fixed point has no dynamics to grow, which is why the Lyapunov exponents of such a run are ~0.

TYPE: callable DEFAULT: None

rtol

Integration tolerances. Default :data:~tsdynamics._utils.tolerances.DDE_RTOL / :data:~tsdynamics._utils.tolerances.DDE_ATOL (both 1e-3) — the DDE family's own, deliberately looser than the ODE default because the method of steps lands on every output sample, so dt bounds the step and the tolerance is largely inert. See the class docstring for the measurement. Tightening is safe (no stall) if you need it.

TYPE: float DEFAULT: None

atol

Integration tolerances. Default :data:~tsdynamics._utils.tolerances.DDE_RTOL / :data:~tsdynamics._utils.tolerances.DDE_ATOL (both 1e-3) — the DDE family's own, deliberately looser than the ODE default because the method of steps lands on every output sample, so dt bounds the step and the tolerance is largely inert. See the class docstring for the measurement. Tightening is safe (no stall) if you need it.

TYPE: float DEFAULT: None

backend

Which evaluator drives the DDE engine. Defaults to _default_backend ("jit"). Both 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.

  • "jit" (default) — the Cranelift JIT, compiled once per distinct system and served from the process-wide compiled-evaluator cache.
  • "interp" — the SSA-tape interpreter; bit-for-bit identical, and the way to avoid the one-off compile.

Only constant delays lower; a state-dependent delay raises. backend="reference" is unsupported for DDEs (there is no pure-Python delay integrator).

.. versionchanged:: 6.0 The default moved from "interp" to "jit", once the v6 compiled-evaluator cache removed the JIT's per-call recompile.

TYPE: ('jit', 'interp') DEFAULT: "jit"

solver

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 name another explicit kernel with solver= if you need one.

TYPE: str DEFAULT: "rk45"

seed

Seed for the random initial-condition draw (the constant past when history is None) — the same meaning seed= has on every other family's trajectory producer and on the constructor. (An SDE has a second source of randomness, so there seed= seeds the noise path as well as the draw.) Inert unless a draw actually happens: an explicit ic, an already-resolved self.ic and a class-level default_ic all take priority. The resolved seed is recorded on traj.meta["ic_seed"].

.. versionadded:: 6.0

TYPE: int DEFAULT: None

transient

Leading stretch of the run to discard, in time units (the same unit as final_time). The window is extended to transient + final_time and everything before transient dropped. One word on every family, in that family's own horizon unit: time here, iterations on a map.

.. versionadded:: 6.0

TYPE: float DEFAULT: 0.0

RETURNS DESCRIPTION
Trajectory
Source code in src/tsdynamics/families/delay.py
def run(
    self,
    final_time: float = 100.0,
    dt: float | None = None,
    *,
    ic: Any | None = None,
    history: History = None,
    transient: float = 0.0,
    solver: str | None = None,
    rtol: float | None = None,
    atol: float | None = None,
    backend: str | None = None,
    seed: int | None = None,
    **solver_options: Any,
) -> Trajectory:
    """
    Integrate the delay system and return a :class:`~tsdynamics.families.Trajectory`.

    ``run`` is **the** trajectory verb — one word on flows, maps, delay and
    stochastic systems.  ``integrate`` and ``trajectory`` were two more names
    for this method and are gone in v6.

    The signature is **closed**.  Before v6 this method accepted — and
    silently dropped — ``max_step``, ``t0``, ``events`` and outright typos:
    the run completed, the number was wrong, and nothing said so.

    Parameters
    ----------
    final_time : float
        End of the integration window, **in time units** — the horizon word
        for a flow.  Default 100.0.
    dt : float, optional
        Output sampling interval, **in time units**.  ``None`` (the default)
        means this family's ``_default_dt``, which ``system.info`` prints
        under ``defaults``.  The method of steps lands on every sample, so
        ``dt`` also bounds the internal step here (see ``rtol``/``atol``).
    ic : array-like, optional
        The **constant** past, used when ``history`` is ``None`` — ``dim``
        numbers, one per state component, held for every ``s ≤ 0``.  Falls
        back to ``self.ic``, then a random draw.
    history : callable, optional
        The past, as a function of the **lag**: ``h(s)`` is called with
        ``s ≤ 0`` (``s = 0`` is the start of the run, ``s = -tau`` one
        delay before it) and must return ``dim`` numbers — the state at that
        moment of the pre-history.  It is never called with ``s > 0``.
        ``None`` (the default) uses a constant past equal to ``ic``::

            mg.run(final_time=500.0, history=lambda s: [1.0 + 0.1 * np.sin(0.2 * s)])

        A constant past sitting exactly on a fixed point has no dynamics to
        grow, which is why the Lyapunov exponents of such a run are ~0.
    rtol, atol : float
        Integration tolerances.  Default
        :data:`~tsdynamics._utils.tolerances.DDE_RTOL` /
        :data:`~tsdynamics._utils.tolerances.DDE_ATOL` (both ``1e-3``) — the
        DDE family's own, deliberately looser than the ODE default because
        the method of steps lands on every output sample, so ``dt`` bounds
        the step and the tolerance is largely inert.  See the class
        docstring for the measurement.  Tightening is *safe* (no stall) if
        you need it.
    backend : {"jit", "interp"}, optional
        Which evaluator drives the DDE engine.  Defaults to
        ``_default_backend`` (``"jit"``).  Both 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.

        - ``"jit"`` (default) — the Cranelift JIT, compiled once per distinct
          system and served from the process-wide compiled-evaluator cache.
        - ``"interp"`` — the SSA-tape interpreter; bit-for-bit identical, and
          the way to avoid the one-off compile.

        Only **constant** delays lower; a state-dependent delay raises.
        ``backend="reference"`` is unsupported for DDEs (there is no
        pure-Python delay integrator).

        .. versionchanged:: 6.0
            The default moved from ``"interp"`` to ``"jit"``, once the v6
            compiled-evaluator cache removed the JIT's per-call recompile.
    solver : 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 name another explicit
        kernel with ``solver=`` if you need one.
    seed : int, optional
        Seed for the **random initial-condition draw** (the constant past when
        ``history`` is ``None``) — the same meaning ``seed=`` has on every
        other family's trajectory producer and on the constructor.  (An SDE
        has a second source of randomness, so there ``seed=`` seeds the
        noise path as well as the draw.)  Inert unless a draw actually
        happens: an explicit ``ic``, an already-resolved ``self.ic`` and a
        class-level ``default_ic`` all take priority.  The resolved seed is
        recorded on ``traj.meta["ic_seed"]``.

        .. versionadded:: 6.0
    transient : float, optional
        Leading stretch of the run to discard, in **time units** (the same
        unit as ``final_time``).  The window is extended to
        ``transient + final_time`` and everything before ``transient``
        dropped.  One word on every family, in that family's **own horizon
        unit**: time here, **iterations** on a map.

        .. versionadded:: 6.0

    Returns
    -------
    Trajectory
    """
    reject_unknown_run_keywords(self, solver_options, family="dde", accepted=_DDE_RUN_KEYWORDS)
    dt = self._default_dt if dt is None else dt
    method = _resolve_dde_solver(self._default_method if solver is None else solver)
    transient = resolve_transient(transient, discrete=False)
    if transient > 0.0:
        traj = self.run(
            final_time + transient,
            dt,
            ic=ic,
            history=history,
            rtol=rtol,
            atol=atol,
            backend=backend,
            solver=method,
            seed=seed,
        )
        return traj.after(transient)
    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,
        seed=seed,
    )

DiscreteMap

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

Bases: SystemBase, ABC

Base class for discrete maps iterated on the engine.

Subclass contract
  1. Declare params = {...} and dim = N.
  2. Implement _step as a @staticmethod. Parameters arrive as positional arguments in the order they appear in the class-level params dict. It is :func:abc.abstractmethod and — since this class is an :class:abc.ABC, like every other family base — a subclass that omits it cannot be instantiated. (Before v6 DiscreteMap was the one family base that was not an ABC, so the abstract marker was inert and a missing kernel surfaced far downstream as a mystifying TapeCompileError from the lowering pass.)
  3. _jacobian is optional: by default the library differentiates _step symbolically, exactly as ContinuousSystem autogenerates its Jacobian from _equations. Write one only when the step cannot be traced (a Python if on the state) or when a one-sided slope on a discontinuity is the meaningful answer. Every catalogue map still ships a hand-written one, which the test suite uses as a cross-check oracle.
Iteration

iterate lowers _step to an in-process IR tape and runs the engine's native map loop, on a tape lowered and JIT-compiled once per system (both memoised). 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.run(steps=10_000)
>>> t_idx, X = traj.unpack()   # the two columns
>>> exps = h._lyapunov_spectrum(n=5_000)
>>> h_variant = h.with_params(a=1.2)
>>> traj2 = h_variant.run(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,
    seed: int | None = None,
    **param_kwargs: Any,
) -> 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.
        An ``ic`` given here is **explicit**: it is never silently replaced by
        a random draw (a diverging one raises instead).
    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.
    seed : int, optional
        Seed for the **initial-condition draw** used when neither ``ic`` nor
        ``system.info.default_ic`` supplies one.  The draw runs on a private
        :class:`numpy.random.Generator`, so it is reproducible *and* never
        touches the global ``numpy.random`` stream.  When omitted a fresh
        OS-entropy seed is drawn on first use and recorded on
        ``traj.meta["ic_seed"]``, so any run can be reproduced exactly by
        passing that value back as ``seed=``.
    **param_kwargs
        Per-parameter overrides given as plain keywords — ``Lorenz(sigma=12.0)``
        is exactly ``Lorenz(params={"sigma": 12.0})``.  Each name must be a
        declared parameter of this class; an unknown one raises and names the
        valid options.  A name given **both** here and inside ``params=`` is an
        error rather than a silent precedence rule.

    Raises
    ------
    InvalidParameterError
        If ``params`` or ``**param_kwargs`` contains a key that is not a
        declared parameter, or if the same parameter is given twice (once in
        ``params=`` and once as a keyword).

    Examples
    --------
    >>> Lorenz(sigma=12.0).sigma            # the natural spelling
    12.0
    >>> Lorenz(params={"sigma": 12.0}).sigma   # equivalent, still supported
    12.0
    """
    # Build ParamSet from class defaults + constructor overrides.  The two
    # override channels (``params=`` and free keywords) are merged here, and
    # both are validated against the *declared* parameter names.
    defaults = dict(type(self).params)
    _reject_state_as_params(type(self).__name__, params, defaults)
    overrides: dict[str, Any] = dict(params) if params else {}
    if param_kwargs:
        duplicated = sorted(set(param_kwargs) & set(overrides))
        if duplicated:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: parameter(s) {duplicated} given twice — "
                f"once in params= and once as a keyword. Pass each parameter "
                f"exactly once (there is deliberately no precedence rule)."
            )
        overrides.update(param_kwargs)
    if overrides:
        unknown = set(overrides) - set(defaults)
        if unknown:
            # A misspelt parameter is nearly always one character off a
            # declared one, so name the closest match and spell the whole
            # constructor call rather than only listing what was valid.
            import difflib

            from tsdynamics.errors import InvalidParameterError

            cls_name = type(self).__name__
            guesses = {
                bad: (
                    [n for n in sorted(defaults) if n.lower() == bad.lower()]
                    or difflib.get_close_matches(bad, sorted(defaults), n=1, cutoff=0.5)
                )
                for bad in sorted(unknown)
            }
            fixed = {bad: g[0] for bad, g in guesses.items() if g}
            did_you_mean = (
                " Did you mean " + ", ".join(f"{b!r} -> {g!r}" for b, g in fixed.items()) + "?"
                if fixed
                else ""
            )
            call = ", ".join(
                f"{fixed.get(bad, bad)}={overrides[bad]!r}" for bad in sorted(unknown)
            )
            raise InvalidParameterError(
                f"{cls_name}: unknown parameter(s) {sorted(unknown)}. "
                f"Declared: {sorted(defaults)}.{did_you_mean}"
                + remedy(
                    f"{cls_name}({call})"
                    if fixed
                    else f"{cls_name}({sorted(defaults)[0]}={defaults[sorted(defaults)[0]]!r})",
                    lead="Spell it as:" if fixed else "The declared parameters are set like:",
                )
            )
        defaults.update(overrides)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    if resolved_dim is None:
        # Every downstream consumer does ``int(self.dim)``, so an undeclared
        # dimension used to surface as NumPy's ``int() argument must be ...
        # not 'NoneType'`` from deep inside resolve_ic -- which names neither
        # the class nor the attribute.  Writing a system class is the first
        # thing a new user does, so refuse it here, where the fix is one line.
        raise InvalidInputError(
            f"{type(self).__name__} does not declare its state-space dimension, "
            f"so the library cannot tell how many components its state has. "
            f"`dim` is the number of equations `_equations` returns."
            + remedy(
                f"class {type(self).__name__}({type(self).__bases__[0].__name__}):",
                "    dim = 3                       # the number of state components",
                '    variables = ("x", "y", "z")   # optional, names them',
                lead="Declare it on the class:",
            )
            + remedy(
                f"{type(self).__name__}(dim=3)",
                lead="...or pass it for this one instance:",
            )
        )
    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_explicit`` records whether the CURRENT
    # ``self.ic`` was chosen by the user (constructor / an explicit ``ic=``)
    # or merely auto-resolved (a random draw).  A user-chosen IC must never
    # be silently swapped for a random one — see ``DiscreteMap.run``.
    # ``np.array(..., copy=True)``, not ``asarray``: a float64 array passed in
    # would otherwise be *shared* with the caller, so mutating either side
    # silently moved the other's initial condition.  It is also what made
    # ``with_params`` (which forwards ``ic=self.ic``) hand back a system
    # aliasing the source's ic array.
    ic_arr = np.array(ic, dtype=float, copy=True) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)
    object.__setattr__(self, "_ic_explicit", ic is not None)

    # The private initial-condition RNG (see ``resolve_ic``).  Only the seed
    # is persisted; the Generator itself is a transient cache.
    object.__setattr__(self, "_ic_seed", None if seed is None else int(seed))
    object.__setattr__(self, "_ic_rng", None)

jacobian

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

Return the tangent map J = df/dx at state u.

The same verb, the same signature and the same return type as a flow's :meth:~tsdynamics.families.ContinuousSystem.jacobian — a map simply has no t dependence, so the argument is accepted and ignored.

PARAMETER DESCRIPTION
u

The state to linearise at.

TYPE: (array - like, shape(dim))

t

Ignored; present so a caller holding any system can write sys.jacobian(u, t).

TYPE: float DEFAULT: 0.0

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

Examples:

>>> import tsdynamics as ts
>>> ts.systems.Henon().jacobian([0.5, 0.2]).shape
(2, 2)
Source code in src/tsdynamics/families/discrete.py
def jacobian(self, u: Any, t: float = 0.0) -> np.ndarray:
    """Return the tangent map ``J = df/dx`` at state ``u``.

    The same verb, the same signature and the same return type as a flow's
    :meth:`~tsdynamics.families.ContinuousSystem.jacobian` — a map simply
    has no ``t`` dependence, so the argument is accepted and ignored.

    Parameters
    ----------
    u : array-like, shape (dim,)
        The state to linearise at.
    t : float, optional
        Ignored; present so a caller holding *any* system can write
        ``sys.jacobian(u, t)``.

    Returns
    -------
    ndarray, shape (dim, dim)

    Examples
    --------
    >>> import tsdynamics as ts
    >>> ts.systems.Henon().jacobian([0.5, 0.2]).shape
    (2, 2)
    """
    state = np.asarray(u, dtype=float)
    jac = type(self)._jacobian(state, *cast(Any, self.params).as_tuple())
    n = int(cast(int, self.dim))
    return np.asarray(jac, dtype=float).reshape(n, n)

reinit

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

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

PARAMETER DESCRIPTION
u

The state to restart from — dim numbers. Resolved like ic= on :meth:run.

TYPE: array - like DEFAULT: None

t

The iteration count to restart the clock at, in iterations (a map has no time). Default 0.

TYPE: int DEFAULT: None

params

Parameter overrides applied to this system in place before the restart. (with_params is the non-mutating spelling.)

TYPE: dict DEFAULT: None

RAISES DESCRIPTION
InvalidParameterError

For any other keyword — refused by name. This used to end in a bare TypeError from the interpreter, so except InvalidParameterError caught a typo on a flow and missed the identical typo on a map.

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,
    **unknown: Any,
) -> None:
    """(Re)start stepping from state ``u`` at iteration count ``t``.

    Parameters
    ----------
    u : array-like, optional
        The state to restart from — ``dim`` numbers.  Resolved like ``ic=``
        on :meth:`run`.
    t : int, optional
        The iteration count to restart the clock at, **in iterations** (a
        map has no time).  Default 0.
    params : dict, optional
        Parameter overrides applied to this system **in place** before the
        restart.  (``with_params`` is the non-mutating spelling.)

    Raises
    ------
    InvalidParameterError
        For any other keyword — refused **by name**.  This used to end in a
        bare ``TypeError`` from the interpreter, so
        ``except InvalidParameterError`` caught a typo on a flow and missed
        the identical typo on a map.
    """
    reject_unknown_run_keywords(
        self,
        unknown,
        family="map",
        accepted=("t", "params"),
        verb="reinit",
    )
    if params:
        for k, v in params.items():
            self.params[k] = v
    # ``resolve_ic`` commits the resolved IC to ``self.ic`` up front, so a
    # malformed ``u`` (wrong length → a reshape ValueError) or a malformed
    # ``t`` must not leave a half-applied IC behind (the ``_dispatch`` /
    # ``iterate`` contract, applied to the stepping entry point).
    with self._ic_rollback():
        state = self._resolve_ic(u)
        n_now = int(t) if t is not None else 0
    self._state_now = state
    self._n_now = n_now

step

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

Advance n_or_dt iterations and return the new state.

n_or_dt is an iteration count here — the one family where it is not a time increment, because a map has no clock. 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_or_dt`` iterations and return the new state.

    ``n_or_dt`` is an **iteration count** here — the one family where it is
    not a time increment, because a map has no clock.  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() -> int

Return the current iteration count — an int, because it counts.

A flow's clock is a float time; a map's is a whole number of iterates, and reporting 7.0 where the docstring promised a count made the two readings of time() indistinguishable at a call site.

Source code in src/tsdynamics/families/discrete.py
def time(self) -> int:
    """Return the current iteration count — an ``int``, because it counts.

    A flow's clock is a float time; a map's is a whole number of iterates,
    and reporting ``7.0`` where the docstring promised a count made the two
    readings of ``time()`` indistinguishable at a call site.
    """
    return int(self._n_now)

run

run(
    steps: int = 1000,
    *,
    ic: Any | None = None,
    transient: int = 0,
    backend: str | None = None,
    seed: int | None = None,
    max_retries: int = 10,
    **solver_options: Any,
) -> Trajectory

Iterate the map and return a :class:~tsdynamics.families.Trajectory.

run is the trajectory verb — one word on flows, maps, delay and stochastic systems. iterate and trajectory were two more names for this method and are gone in v6.

A map's horizon is steps, a count of iterations, because a map has no continuous time; final_time is refused by name rather than reinterpreted. The positional form reads the same on both families: hen.run(1000), lor.run(100.0).

PARAMETER DESCRIPTION
steps

How many iterations to run — the horizon word for a map, counted in iterations. (A flow measures its horizon in time and takes final_time; a map has no clock, so final_time is refused by name.) Default 1000. The trajectory holds steps + 1 rows: the state it started from, then one per iteration — the same N + 1 a flow returns for final_time / dt == N.

TYPE: int DEFAULT: 1000

ic

Initial state — dim numbers, one per state component. Falls back to self.ic, then a random draw. Passing it here does not change self.ic: one call, one run.

TYPE: array - like DEFAULT: None

max_retries

Retry with a new random IC if divergence is detected — only when the initial condition was not chosen by the user. An explicit ic (here or on the constructor) that diverges raises instead of being silently swapped for a random one that traces a different orbit.

TYPE: int DEFAULT: 10

seed

Seed for the random-IC fallback and for the divergence retries, so an unseeded-IC run is reproducible. Equivalent to the constructor's seed=, and the same meaning it has on every other family; the resolved seed is recorded on traj.meta["ic_seed"].

TYPE: int DEFAULT: None

backend

Where the iteration runs. Defaults to _default_backend ("jit").

  • "jit" (default) — the Rust engine's native map loop driven by the Cranelift JIT: the lowered tape compiled to native code and served from a process-wide compiled-evaluator cache, so the compile is paid once per distinct system rather than on every call.
  • "interp" — the same native map loop driven by the SSA-tape interpreter. Bit-for-bit identical to the JIT; marginally faster on very small tapes, and the way to skip the one-off compile. Both require the compiled extension (:mod:tsdynamics._rust); until it is built they raise :class:~tsdynamics._engine.run.EngineNotAvailableError.
  • "reference" — the lowered next-state tape, iterated in pure Python. Not for production use: it is the dependency-light oracle the engine is validated against.

.. versionchanged:: 6.0 The default moved from "interp" to "jit", once the v6 compiled-evaluator cache removed the JIT's per-call recompile.

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: ('jit', 'interp', 'reference') DEFAULT: "jit"

transient

Leading stretch of the orbit to discard, in iterations (the same unit as steps). transient + steps iterations are run and the first transient dropped, so the returned trajectory still has steps samples. One word on every family, in that family's own horizon unit: iterations here, time on a flow (see :meth:ContinuousSystem.run).

.. versionadded:: 6.0 Only the retired trajectory verb used to accept it, so a user who found run could not discard a transient at all.

TYPE: int DEFAULT: 0

RETURNS DESCRIPTION
Trajectory

steps + 1 rows — the initial condition first, then one per iteration — with t = arange(steps + 1) (integer step indices, not float times), so y[k] is the state at t[k] = k.

.. versionchanged:: 6.0 The initial condition used to be dropped, so t[0] = 0 labelled :math:x_1 and traj["x"][n] was :math:x_{n+1} — measured, Logistic(r=2.8).run(steps=4, ic=[0.1]).y[0] was 0.252 (:math:f(0.1)) while Lorenz().run(ic=[1, 2, 3]).y[0] was [1, 2, 3]. Every cobweb therefore started one iterate late, and the two families disagreed about what a trajectory's first row means. They agree now.

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

InvalidParameterError

If steps < 1, or an unrecognised keyword is passed — including a flow keyword (final_time / dt / method / a tolerance), which a map has no meaning for. **kwargs exists only to catch those: nothing is silently dropped.

Source code in src/tsdynamics/families/discrete.py
def run(
    self,
    steps: int = 1000,
    *,
    ic: Any | None = None,
    transient: int = 0,
    backend: str | None = None,
    seed: int | None = None,
    max_retries: int = 10,
    **solver_options: Any,
) -> Trajectory:
    """
    Iterate the map and return a :class:`~tsdynamics.families.Trajectory`.

    ``run`` is **the** trajectory verb — one word on flows, maps, delay and
    stochastic systems.  ``iterate`` and ``trajectory`` were two more names
    for this method and are gone in v6.

    A map's horizon is ``steps``, a **count of iterations**, because a map
    has no continuous time; ``final_time`` is refused by name rather than
    reinterpreted.  The positional form reads the same on both families:
    ``hen.run(1000)``, ``lor.run(100.0)``.

    Parameters
    ----------
    steps : int
        How many iterations to **run** — the horizon word for a map, counted
        in **iterations**.  (A flow measures its horizon in time and takes
        ``final_time``; a map has no clock, so ``final_time`` is refused by
        name.)  Default 1000.  The trajectory holds ``steps + 1`` rows: the
        state it started from, then one per iteration — the same ``N + 1`` a
        flow returns for ``final_time / dt == N``.
    ic : array-like, optional
        Initial state — ``dim`` numbers, one per state component.  Falls
        back to ``self.ic``, then a random draw.  **Passing it here does not
        change ``self.ic``**: one call, one run.
    max_retries : int
        Retry with a new random IC if divergence is detected — only when the
        initial condition was **not chosen by the user**.  An explicit ``ic``
        (here *or* on the constructor) that diverges raises instead of being
        silently swapped for a random one that traces a different orbit.
    seed : int, optional
        Seed for the random-IC fallback and for the divergence retries, so an
        unseeded-IC run is reproducible.  Equivalent to the constructor's
        ``seed=``, and the same meaning it has on every other family; the
        resolved seed is recorded on ``traj.meta["ic_seed"]``.
    backend : {"jit", "interp", "reference"}, optional
        Where the iteration runs.  Defaults to ``_default_backend``
        (``"jit"``).

        - ``"jit"`` (default) — the Rust engine's native map loop driven by
          the **Cranelift JIT**: the lowered tape compiled to native code and
          served from a process-wide compiled-evaluator cache, so the compile
          is paid once per distinct system rather than on every call.
        - ``"interp"`` — the same native map loop driven by the **SSA-tape
          interpreter**.  Bit-for-bit identical to the JIT; marginally faster
          on very small tapes, and the way to skip the one-off compile.
          Both require the compiled extension (:mod:`tsdynamics._rust`);
          until it is built they raise
          :class:`~tsdynamics._engine.run.EngineNotAvailableError`.
        - ``"reference"`` — the lowered next-state tape, iterated in pure
          Python.  Not for production use: it is the dependency-light oracle
          the engine is validated against.

        .. versionchanged:: 6.0
            The default moved from ``"interp"`` to ``"jit"``, once the v6
            compiled-evaluator cache removed the JIT's per-call recompile.

        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`.
    transient : int, optional
        Leading stretch of the orbit to discard, in **iterations** (the same
        unit as ``steps``).  ``transient + steps`` iterations are run and the
        first ``transient`` dropped, so the returned trajectory still has
        ``steps`` samples.  One word on every family, in that family's
        **own horizon unit**: iterations here, **time** on a flow (see
        :meth:`ContinuousSystem.run`).

        .. versionadded:: 6.0
            Only the retired ``trajectory`` verb used to accept it, so a
            user who found ``run`` could not discard a transient at all.

    Returns
    -------
    Trajectory
        ``steps + 1`` rows — **the initial condition first**, then one per
        iteration — with ``t = arange(steps + 1)`` (integer step indices, not
        float times), so ``y[k]`` is the state at ``t[k] = k``.

        .. versionchanged:: 6.0
            The initial condition used to be dropped, so ``t[0] = 0`` labelled
            :math:`x_1` and ``traj["x"][n]`` was :math:`x_{n+1}` — measured,
            ``Logistic(r=2.8).run(steps=4, ic=[0.1]).y[0]`` was ``0.252``
            (:math:`f(0.1)`) while ``Lorenz().run(ic=[1, 2, 3]).y[0]`` was
            ``[1, 2, 3]``.  Every cobweb therefore started one iterate late,
            and the two families disagreed about what a trajectory's first
            row means.  They agree now.

    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).
    InvalidParameterError
        If ``steps < 1``, or an unrecognised keyword is passed — including a
        *flow* keyword (``final_time`` / ``dt`` / ``method`` / a tolerance),
        which a map has no meaning for.  ``**kwargs`` exists only to catch
        those: nothing is silently dropped.
    """
    reject_unknown_run_keywords(self, solver_options, family="map", accepted=_MAP_RUN_KEYWORDS)
    steps = int(steps)
    if steps < 1:
        raise InvalidParameterError(
            f"steps is a number of iterations, so it must be >= 1; got {steps}."
            + remedy(f"{type(self).__name__}().run(steps=1000)")
        )
    drop = int(resolve_transient(transient, discrete=True))
    backend = backend if backend is not None else self._default_backend

    with self._ic_rollback():
        traj = self._iterate_with_retries(
            steps=steps + drop, ic=ic, max_retries=max_retries, backend=backend, seed=seed
        )
    return traj[drop:] if drop > 0 else traj

StochasticSystem

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

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.run(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,
    seed: int | None = None,
    **param_kwargs: Any,
) -> 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.
        An ``ic`` given here is **explicit**: it is never silently replaced by
        a random draw (a diverging one raises instead).
    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.
    seed : int, optional
        Seed for the **initial-condition draw** used when neither ``ic`` nor
        ``system.info.default_ic`` supplies one.  The draw runs on a private
        :class:`numpy.random.Generator`, so it is reproducible *and* never
        touches the global ``numpy.random`` stream.  When omitted a fresh
        OS-entropy seed is drawn on first use and recorded on
        ``traj.meta["ic_seed"]``, so any run can be reproduced exactly by
        passing that value back as ``seed=``.
    **param_kwargs
        Per-parameter overrides given as plain keywords — ``Lorenz(sigma=12.0)``
        is exactly ``Lorenz(params={"sigma": 12.0})``.  Each name must be a
        declared parameter of this class; an unknown one raises and names the
        valid options.  A name given **both** here and inside ``params=`` is an
        error rather than a silent precedence rule.

    Raises
    ------
    InvalidParameterError
        If ``params`` or ``**param_kwargs`` contains a key that is not a
        declared parameter, or if the same parameter is given twice (once in
        ``params=`` and once as a keyword).

    Examples
    --------
    >>> Lorenz(sigma=12.0).sigma            # the natural spelling
    12.0
    >>> Lorenz(params={"sigma": 12.0}).sigma   # equivalent, still supported
    12.0
    """
    # Build ParamSet from class defaults + constructor overrides.  The two
    # override channels (``params=`` and free keywords) are merged here, and
    # both are validated against the *declared* parameter names.
    defaults = dict(type(self).params)
    _reject_state_as_params(type(self).__name__, params, defaults)
    overrides: dict[str, Any] = dict(params) if params else {}
    if param_kwargs:
        duplicated = sorted(set(param_kwargs) & set(overrides))
        if duplicated:
            from tsdynamics.errors import InvalidParameterError

            raise InvalidParameterError(
                f"{type(self).__name__}: parameter(s) {duplicated} given twice — "
                f"once in params= and once as a keyword. Pass each parameter "
                f"exactly once (there is deliberately no precedence rule)."
            )
        overrides.update(param_kwargs)
    if overrides:
        unknown = set(overrides) - set(defaults)
        if unknown:
            # A misspelt parameter is nearly always one character off a
            # declared one, so name the closest match and spell the whole
            # constructor call rather than only listing what was valid.
            import difflib

            from tsdynamics.errors import InvalidParameterError

            cls_name = type(self).__name__
            guesses = {
                bad: (
                    [n for n in sorted(defaults) if n.lower() == bad.lower()]
                    or difflib.get_close_matches(bad, sorted(defaults), n=1, cutoff=0.5)
                )
                for bad in sorted(unknown)
            }
            fixed = {bad: g[0] for bad, g in guesses.items() if g}
            did_you_mean = (
                " Did you mean " + ", ".join(f"{b!r} -> {g!r}" for b, g in fixed.items()) + "?"
                if fixed
                else ""
            )
            call = ", ".join(
                f"{fixed.get(bad, bad)}={overrides[bad]!r}" for bad in sorted(unknown)
            )
            raise InvalidParameterError(
                f"{cls_name}: unknown parameter(s) {sorted(unknown)}. "
                f"Declared: {sorted(defaults)}.{did_you_mean}"
                + remedy(
                    f"{cls_name}({call})"
                    if fixed
                    else f"{cls_name}({sorted(defaults)[0]}={defaults[sorted(defaults)[0]]!r})",
                    lead="Spell it as:" if fixed else "The declared parameters are set like:",
                )
            )
        defaults.update(overrides)
    object.__setattr__(self, "params", ParamSet(defaults))

    # dim: constructor arg > class attribute
    resolved_dim = dim if dim is not None else type(self).dim
    if resolved_dim is None:
        # Every downstream consumer does ``int(self.dim)``, so an undeclared
        # dimension used to surface as NumPy's ``int() argument must be ...
        # not 'NoneType'`` from deep inside resolve_ic -- which names neither
        # the class nor the attribute.  Writing a system class is the first
        # thing a new user does, so refuse it here, where the fix is one line.
        raise InvalidInputError(
            f"{type(self).__name__} does not declare its state-space dimension, "
            f"so the library cannot tell how many components its state has. "
            f"`dim` is the number of equations `_equations` returns."
            + remedy(
                f"class {type(self).__name__}({type(self).__bases__[0].__name__}):",
                "    dim = 3                       # the number of state components",
                '    variables = ("x", "y", "z")   # optional, names them',
                lead="Declare it on the class:",
            )
            + remedy(
                f"{type(self).__name__}(dim=3)",
                lead="...or pass it for this one instance:",
            )
        )
    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_explicit`` records whether the CURRENT
    # ``self.ic`` was chosen by the user (constructor / an explicit ``ic=``)
    # or merely auto-resolved (a random draw).  A user-chosen IC must never
    # be silently swapped for a random one — see ``DiscreteMap.run``.
    # ``np.array(..., copy=True)``, not ``asarray``: a float64 array passed in
    # would otherwise be *shared* with the caller, so mutating either side
    # silently moved the other's initial condition.  It is also what made
    # ``with_params`` (which forwards ``ic=self.ic``) hand back a system
    # aliasing the source's ic array.
    ic_arr = np.array(ic, dtype=float, copy=True) if ic is not None else None
    object.__setattr__(self, "ic", ic_arr)
    object.__setattr__(self, "_ic_explicit", ic is not None)

    # The private initial-condition RNG (see ``resolve_ic``).  Only the seed
    # is persisted; the Generator itself is a transient cache.
    object.__setattr__(self, "_ic_seed", None if seed is None else int(seed))
    object.__setattr__(self, "_ic_rng", None)

rhs

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

Return the drift f(u, t) — the deterministic skeleton's field.

The twin of :meth:jacobian, and named on the same principle: an SDE has two kernels, so "the right-hand side" is ambiguous by itself, and this one is the drift — the field whose Jacobian :meth:jacobian returns and whose equilibria the deterministic analyses find. The diffusion is not a vector field over (u, t) alone (it multiplies dW), so it has no place here.

PARAMETER DESCRIPTION
u

TYPE: (array - like, shape(dim))

t

TYPE: float DEFAULT: 0.0

RETURNS DESCRIPTION
(ndarray, shape(dim))

Examples:

>>> import tsdynamics as ts
>>> ou = ts.systems.OrnsteinUhlenbeck()
>>> ou.rhs([1.0]).shape
(1,)
Source code in src/tsdynamics/families/stochastic.py
def rhs(self, u: Any, t: float = 0.0) -> np.ndarray:
    """Return the **drift** ``f(u, t)`` — the deterministic skeleton's field.

    The twin of :meth:`jacobian`, and named on the same principle: an SDE has
    two kernels, so "the right-hand side" is ambiguous by itself, and this
    one is the drift — the field whose Jacobian :meth:`jacobian` returns and
    whose equilibria the deterministic analyses find.  The diffusion is not a
    vector field over ``(u, t)`` alone (it multiplies ``dW``), so it has no
    place here.

    Parameters
    ----------
    u : array-like, shape (dim,)
    t : float, optional

    Returns
    -------
    ndarray, shape (dim,)

    Examples
    --------
    >>> import tsdynamics as ts
    >>> ou = ts.systems.OrnsteinUhlenbeck()
    >>> ou.rhs([1.0]).shape
    (1,)
    """
    import symengine

    from tsdynamics._engine.symbols import state_time_symbols

    dim = int(cast(int, self.dim))
    y, t_sym = state_time_symbols()
    names = list(cast(Any, self.params))
    syms = {k: symengine.Symbol(k) for k in names}
    f = list(type(self)._drift(y, t_sym, **syms))
    state = [symengine.Symbol(f"_u{i}") for i in range(dim)]
    subs = {y(i): state[i] for i in range(dim)}
    args = state + [t_sym] + [syms[k] for k in names]
    fn = symengine.Lambdify(args, [symengine.sympify(e).subs(subs) for e in f])
    vals = [float(self.params[k]) for k in names]
    arg = np.concatenate([np.asarray(u, dtype=float).ravel(), [t], vals])
    return np.asarray(fn(arg), dtype=float).reshape(dim)

jacobian

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

Return the drift Jacobian d f/d u at state u.

An SDE has two kernels, so "the Jacobian" is ambiguous by itself; this one is the drift's, which is the object linear stability of the deterministic skeleton is read from. The diffusion's own derivative (Milstein's correction) is built by the lowering, not exposed here.

PARAMETER DESCRIPTION
u

TYPE: (array - like, shape(dim))

t

TYPE: float DEFAULT: 0.0

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

Examples:

>>> import tsdynamics as ts
>>> ts.systems.OrnsteinUhlenbeck().jacobian([1.0]).shape
(1, 1)
Source code in src/tsdynamics/families/stochastic.py
def jacobian(self, u: Any, t: float = 0.0) -> np.ndarray:
    """Return the **drift** Jacobian ``d f/d u`` at state ``u``.

    An SDE has two kernels, so "the Jacobian" is ambiguous by itself; this
    one is the drift's, which is the object linear stability of the
    deterministic skeleton is read from.  The diffusion's own derivative
    (Milstein's correction) is built by the lowering, not exposed here.

    Parameters
    ----------
    u : array-like, shape (dim,)
    t : float, optional

    Returns
    -------
    ndarray, shape (dim, dim)

    Examples
    --------
    >>> import tsdynamics as ts
    >>> ts.systems.OrnsteinUhlenbeck().jacobian([1.0]).shape
    (1, 1)
    """
    import symengine

    from tsdynamics._engine.symbols import state_time_symbols

    dim = int(cast(int, self.dim))
    y, t_sym = state_time_symbols()
    names = list(cast(Any, self.params))
    syms = {k: symengine.Symbol(k) for k in names}
    f = list(type(self)._drift(y, t_sym, **syms))
    rows = [[symengine.sympify(f[i]).diff(y(j)) for j in range(dim)] for i in range(dim)]
    # ``y(i)`` is a Function application, not a Symbol, and Lambdify needs
    # symbols — substitute after differentiating, exactly as the ODE family's
    # lambdified cache does.
    state = [symengine.Symbol(f"_u{i}") for i in range(dim)]
    subs = {y(i): state[i] for i in range(dim)}
    args = state + [t_sym] + [syms[k] for k in names]
    fn = symengine.Lambdify(
        args, [symengine.sympify(e).subs(subs) for row in rows for e in row]
    )
    vals = [float(self.params[k]) for k in names]
    arg = np.concatenate([np.asarray(u, dtype=float).ravel(), [t], vals])
    return np.asarray(fn(arg), dtype=float).reshape(dim, dim)

reinit

reinit(
    u: Any | None = None,
    *,
    t: float | None = None,
    params: dict[str, Any] | None = None,
    solver: str | None = None,
    seed: int | None = None,
    dt: float | None = None,
    **unknown: Any,
) -> 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

solver

"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,
    solver: str | None = None,
    seed: int | None = None,
    dt: float | None = None,
    **unknown: Any,
) -> 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.
    solver : 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``).
    """
    reject_unknown_run_keywords(
        self,
        unknown,
        family="sde",
        accepted=("t", "params", "solver", "seed", "dt"),
        verb="reinit",
    )
    if params:
        for k, v in params.items():
            self.params[k] = v
    t0 = float(t) if t is not None else 0.0
    # ``resolve_ic`` commits the IC before the method resolution and the tape
    # lowering below, either of which can raise — a failed ``reinit`` must
    # leave the object exactly as it was.
    with self._ic_rollback():
        ic_arr = self._resolve_ic(u)
        canon = self._resolve_method(solver)
        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 n_or_dt and return the new state.

n_or_dt is a time increment, in the same unit as final_time, and for an SDE it is also the noise scale (the increment is drawn ~ N(0, dt)); omitting it uses the stepper's dt — the dt= given to :meth:reinit, else _default_step_dt (0.01, the number every continuous family uses).

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 ``n_or_dt`` and return the new state.

    ``n_or_dt`` is a **time increment**, in the same unit as ``final_time``,
    and for an SDE it is also the noise scale (the increment is drawn
    ``~ N(0, dt)``); omitting it uses the stepper's ``dt`` — the ``dt=`` given
    to :meth:`reinit`, else ``_default_step_dt`` (``0.01``, the number every
    continuous family uses).

    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

run

run(
    final_time: float = 100.0,
    dt: float | None = None,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    transient: float = 0.0,
    solver: str | None = None,
    seed: int | None = None,
    backend: str | None = None,
    **solver_options: Any,
) -> Trajectory

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

run is the trajectory verb — one word on flows, maps, delay and stochastic systems. integrate and trajectory were two more names for this method and are gone in v6.

There is no rtol/atol here: an Ito scheme takes fixed dt steps and carries no embedded error estimate, so dt is both the discretisation and the output grid.

PARAMETER DESCRIPTION
final_time

End of the integration window, in time units — the horizon word for a flow.

TYPE: float DEFAULT: 100.0

dt

Fixed step size and output sampling interval, in time units. 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 — unlike a deterministic flow, coarsening it changes the path, not just the sampling.

TYPE: float DEFAULT: 0.02

t0

Where the integration starts, in time units (the IC is the state at t0).

TYPE: float DEFAULT: 0.0

ic

Initial state at t0 — dim numbers, one per state component. Falls back to self.ic, then a U[0, 1)^dim draw.

TYPE: array - like DEFAULT: None

solver

"euler_maruyama" (default, order 0.5) or "milstein" (order 1.0).

TYPE: str DEFAULT: None

seed

Makes the whole run reproducible: it seeds the noise realisation and the random initial-condition draw. (Every other family has only the second source of randomness, so there seed= means the draw alone — this is the one family where it also fixes the path.) Random if omitted; the resolved seed is recorded in traj.meta["seed"], and the IC seed in traj.meta["ic_seed"].

.. 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 ("jit", the Cranelift JIT served from the compiled-evaluator cache). "jit" / "interp" (the bit-for-bit identical SSA-tape interpreter) 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" pure-Python oracle 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.

.. versionchanged:: 6.0 The default moved from "interp" to "jit", once the v6 compiled-evaluator cache removed the JIT's per-call recompile.

TYPE: ('jit', 'interp', 'reference') DEFAULT: "jit"

transient

Leading stretch of the path to discard, in time units (the same unit as final_time). The window is extended to transient + final_time and everything before t0 + transient dropped. One word on every family, in that family's own horizon unit: time here, iterations on a map.

.. versionadded:: 6.0

TYPE: float DEFAULT: 0.0

RETURNS DESCRIPTION
Trajectory

A container of samples: len(traj) time points, iterating yields (t_i, y_i) pairs; traj.unpack() gives the two column arrays (t, y).

RAISES DESCRIPTION
InvalidParameterError

If an adaptive-solver keyword (rtol / atol / max_step) or events is passed. An SDE is integrated by a fixed-step scheme, so there is no error control to tune and no continuous extension to root-find a crossing on; the keyword is refused by name rather than silently ignored.

Source code in src/tsdynamics/families/stochastic.py
def run(
    self,
    final_time: float = 100.0,
    dt: float | None = None,
    *,
    t0: float = 0.0,
    ic: Any | None = None,
    transient: float = 0.0,
    solver: str | None = None,
    seed: int | None = None,
    backend: str | None = None,
    **solver_options: Any,
) -> Trajectory:
    """
    Integrate the SDE and return a :class:`~tsdynamics.families.Trajectory`.

    ``run`` is **the** trajectory verb — one word on flows, maps, delay and
    stochastic systems.  ``integrate`` and ``trajectory`` were two more names
    for this method and are gone in v6.

    There is no ``rtol``/``atol`` here: an Ito scheme takes fixed ``dt``
    steps and carries no embedded error estimate, so ``dt`` is both the
    discretisation and the output grid.

    Parameters
    ----------
    final_time : float, default 100.0
        End of the integration window, **in time units** — the horizon word
        for a flow.
    dt : float, default 0.02
        Fixed step size *and* output sampling interval, **in time units**.
        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`` — unlike a deterministic flow, coarsening it changes the path,
        not just the sampling.
    t0 : float, default 0.0
        Where the integration **starts**, in time units (the IC is the state
        at ``t0``).
    ic : array-like, optional
        Initial state at ``t0`` — ``dim`` numbers, one per state component.
        Falls back to ``self.ic``, then a ``U[0, 1)^dim`` draw.
    solver : str, optional
        ``"euler_maruyama"`` (default, order 0.5) or ``"milstein"``
        (order 1.0).
    seed : int, optional
        Makes the whole run reproducible: it seeds the **noise realisation**
        *and* the random initial-condition draw.  (Every other family has
        only the second source of randomness, so there ``seed=`` means the
        draw alone — this is the one family where it also fixes the path.)
        Random if omitted; the resolved seed is recorded in
        ``traj.meta["seed"]``, and the IC seed in ``traj.meta["ic_seed"]``.

        .. 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 : {"jit", "interp", "reference"}, optional
        Defaults to ``_default_backend`` (``"jit"``, the Cranelift JIT served
        from the compiled-evaluator cache).  ``"jit"`` / ``"interp"`` (the
        bit-for-bit identical SSA-tape interpreter) 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"`` pure-Python oracle 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.

        .. versionchanged:: 6.0
            The default moved from ``"interp"`` to ``"jit"``, once the v6
            compiled-evaluator cache removed the JIT's per-call recompile.
    transient : float, optional
        Leading stretch of the path to discard, in **time units** (the same
        unit as ``final_time``).  The window is extended to
        ``transient + final_time`` and everything before ``t0 + transient``
        dropped.  One word on every family, in that family's **own horizon
        unit**: time here, **iterations** on a map.

        .. versionadded:: 6.0

    Returns
    -------
    Trajectory
        A container of samples: ``len(traj)`` time points, iterating
        yields ``(t_i, y_i)`` pairs; ``traj.unpack()`` gives the two
        column arrays ``(t, y)``.

    Raises
    ------
    InvalidParameterError
        If an adaptive-solver keyword (``rtol`` / ``atol`` / ``max_step``) or
        ``events`` is passed.  An SDE is integrated by a **fixed-step**
        scheme, so there is no error control to tune and no continuous
        extension to root-find a crossing on; the keyword is refused by name
        rather than silently ignored.
    """
    reject_unknown_run_keywords(self, solver_options, family="sde", accepted=_SDE_RUN_KEYWORDS)
    dt = self._default_dt if dt is None else dt
    method = solver
    transient = resolve_transient(transient, discrete=False)
    if transient > 0.0:
        traj = self.run(
            final_time=final_time + transient,
            dt=dt,
            t0=t0,
            ic=ic,
            solver=method,
            seed=seed,
            backend=backend,
        )
        return traj.after(t0 + transient)
    # The SDE family does not route through ``_dispatch`` (the generic seam
    # cannot carry the noise seed), so it needs the IC rollback guard of its
    # own: a diverging / interrupted path must leave ``self.ic`` untouched.
    with self._ic_rollback():
        return self._integrate_resolved(
            final_time=final_time,
            dt=dt,
            t0=t0,
            ic=ic,
            method=method,
            seed=seed,
            backend=backend,
        )

System

Bases: Protocol

Structural type for steppable dynamical systems.

run

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

Produce a :class:~tsdynamics.families.Trajectory.

The one trajectory verb on every family and every wrapper. The horizon keyword follows the family (final_time for a flow, steps for a map); everything family-independent — ic, seed, backend, transient — is spelled identically everywhere.

Source code in src/tsdynamics/families/protocol.py
def run(self, *args: Any, **kwargs: Any) -> Trajectory:
    """Produce a :class:`~tsdynamics.families.Trajectory`.

    The one trajectory verb on every family and every wrapper.  The horizon
    keyword follows the family (``final_time`` for a flow, ``steps`` for a
    map); everything family-independent — ``ic``, ``seed``, ``backend``,
    ``transient`` — is spelled identically everywhere.
    """
    ...

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.
    """
    ...

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, **kwargs: Any) -> None

Restart the stepper from state u.

Source code in src/tsdynamics/families/protocol.py
def reinit(self, u: Any | None = None, **kwargs: Any) -> None:
    """Restart the stepper from state ``u``."""
    ...