Reference

Metrics¶

Forecast-quality and dynamical metrics for reservoir models, with first-class support for reporting horizons in Lyapunov times — the natural clock of a chaotic attractor — rather than raw integration steps. A valid-prediction time in Lyapunov times (\(\text{VPT}_{\Lambda t} = \text{steps} \cdot dt \cdot \lambda_\mathrm{max}\)) is hardware- and integration-step-independent, so it is comparable across systems and papers.

The conversions, the canonical \(\lambda_\mathrm{max}\) tables, valid_prediction_time, and reservoir_lyapunov are pure PyTorch/NumPy and always available. The series- and System-based paths of largest_lyapunov need the optional resdag[dynamics] extra (tsdynamics, Python ≥ 3.12) — see Work · Lyapunov time for the prose walkthrough.

Forecast accuracy¶

metrics ¶

Reservoir-computing metrics¶

Forecast-quality and dynamical metrics for reservoir models, with first-class support for reporting horizons in Lyapunov times — the natural clock of a chaotic attractor — instead of raw integration steps.

Lyapunov bookkeeping

lyapunov_time, steps_per_lyapunov_time, steps_to_lyapunov, lyapunov_to_steps Pure, dependency-free conversions between step counts and Lyapunov times. CANONICAL_LAMBDA_MAX, MACKEY_GLASS_LAMBDA Literature largest-Lyapunov-exponent tables for resdag's built-in generators. largest_lyapunov Estimate the largest Lyapunov exponent of a generator name (canonical), a measured series, or a tsdynamics System. reservoir_lyapunov Native-torch Benettin/QR estimate of a trained reservoir's own maximal Lyapunov exponent (an edge-of-chaos / Echo-State-Property diagnostic).

Forecast accuracy

valid_prediction_time Contiguous valid forecast horizon (NRMSE below a threshold), reportable in steps or Lyapunov times.

Optional dependency

The series- and System-based Lyapunov estimators require the optional tsdynamics package (pip install resdag[dynamics], Python >= 3.12); every tsdynamics import is lazy and guarded. Everything else — the conversions, the canonical constants, the reservoir Jacobian estimate, and the forecast metrics — is pure PyTorch/NumPy and always available.

Examples:

>>> from resdag.metrics import steps_to_lyapunov, valid_prediction_time
>>> round(steps_to_lyapunov(357, dt=0.02, lambda_max=0.9056), 2)
6.47

valid_prediction_time ¶

valid_prediction_time(pred: Tensor, truth: Tensor, *, threshold: float = 0.3, metric: str = 'nrmse', units: str = 'steps', dt: float | None = None, lambda_max: float | None = None, dataset: str | None = None, **kw: float) -> float

Valid Prediction Time: contiguous steps a forecast stays accurate.

Per-step error is the NRMSE — the RMSE across feature channels, each channel normalised by the per-dimension standard deviation of truth over the (batch, time) axes — reduced across the batch with the median. The VPT is the number of steps, counting from the start, whose error stays strictly below threshold (0 if the very first step is already over).

This matches :func:resdag.hpo.losses.forecast_horizon exactly (that copy is kept out of the optional-HPO import path); see the module docstring.

PARAMETER	DESCRIPTION
`pred`	Forecast and ground truth of shape `(batch, time, dim)`. TYPE: `Tensor`
`truth`	Forecast and ground truth of shape `(batch, time, dim)`. TYPE: `Tensor`
`threshold`	NRMSE threshold below which a step counts as "valid". TYPE: `float` DEFAULT: `0.3`
`metric`	Per-step error metric. Only NRMSE is implemented (the reservoir-computing standard); the argument exists for forward compatibility. TYPE: `'nrmse'` DEFAULT: `"nrmse"`
`units`	Return the horizon in raw steps, or converted to Lyapunov times (`steps * dt * lambda_max`). `"lyapunov"` needs `dt` and `lambda_max`, or a known `dataset` name to infer them. TYPE: `('steps', 'lyapunov')` DEFAULT: `"steps"`
`dt`	Integration step, required for `units="lyapunov"` unless inferable from `dataset`. TYPE: `float` DEFAULT: `None`
`lambda_max`	Largest Lyapunov exponent, required for `units="lyapunov"` unless inferable from a chaotic `dataset`. TYPE: `float` DEFAULT: `None`
`dataset`	A built-in generator name (`"lorenz"`, `"rossler"`, `"mackey_glass"`, ...) used to look up `dt` / `lambda_max` when they are not given explicitly. TYPE: `str` DEFAULT: `None`
`**kw`	Extra keys forwarded to inference; `tau=` selects the Mackey-Glass exponent. TYPE: `float` DEFAULT: `{}`

RETURNS	DESCRIPTION
`float`	The VPT: an integer step count as a `float` when `units="steps"`, or the horizon in Lyapunov times when `units="lyapunov"`.

RAISES	DESCRIPTION
`ValueError`	If `metric` is not `"nrmse"`, `units` is unknown, or `units="lyapunov"` cannot resolve `dt` / `lambda_max`.

Examples:

>>> import torch
>>> pred = torch.zeros(1, 100, 3)
>>> truth = torch.zeros(1, 100, 3)
>>> valid_prediction_time(pred, truth)  # identical -> full horizon
100.0

Source code in src/resdag/metrics/forecast.py

def valid_prediction_time(
    pred: torch.Tensor,
    truth: torch.Tensor,
    *,
    threshold: float = 0.3,
    metric: str = "nrmse",
    units: str = "steps",
    dt: float | None = None,
    lambda_max: float | None = None,
    dataset: str | None = None,
    **kw: float,
) -> float:
    """Valid Prediction Time: contiguous steps a forecast stays accurate.

    Per-step error is the **NRMSE** — the RMSE across feature channels, each
    channel normalised by the per-dimension standard deviation of ``truth`` over
    the ``(batch, time)`` axes — reduced across the batch with the **median**.
    The VPT is the number of steps, counting from the start, whose error stays
    strictly below ``threshold`` (``0`` if the very first step is already over).

    This matches :func:`resdag.hpo.losses.forecast_horizon` exactly (that copy
    is kept out of the optional-HPO import path); see the module docstring.

    Parameters
    ----------
    pred, truth : torch.Tensor
        Forecast and ground truth of shape ``(batch, time, dim)``.
    threshold : float, default=0.3
        NRMSE threshold below which a step counts as "valid".
    metric : {"nrmse"}, default="nrmse"
        Per-step error metric.  Only NRMSE is implemented (the reservoir-computing
        standard); the argument exists for forward compatibility.
    units : {"steps", "lyapunov"}, default="steps"
        Return the horizon in raw steps, or converted to Lyapunov times
        (``steps * dt * lambda_max``).  ``"lyapunov"`` needs ``dt`` and
        ``lambda_max``, or a known ``dataset`` name to infer them.
    dt : float, optional
        Integration step, required for ``units="lyapunov"`` unless inferable
        from ``dataset``.
    lambda_max : float, optional
        Largest Lyapunov exponent, required for ``units="lyapunov"`` unless
        inferable from a chaotic ``dataset``.
    dataset : str, optional
        A built-in generator name (``"lorenz"``, ``"rossler"``,
        ``"mackey_glass"``, ...) used to look up ``dt`` / ``lambda_max`` when
        they are not given explicitly.
    **kw
        Extra keys forwarded to inference; ``tau=`` selects the Mackey-Glass
        exponent.

    Returns
    -------
    float
        The VPT: an integer step count as a ``float`` when ``units="steps"``, or
        the horizon in Lyapunov times when ``units="lyapunov"``.

    Raises
    ------
    ValueError
        If ``metric`` is not ``"nrmse"``, ``units`` is unknown, or
        ``units="lyapunov"`` cannot resolve ``dt`` / ``lambda_max``.

    Examples
    --------
    >>> import torch
    >>> pred = torch.zeros(1, 100, 3)
    >>> truth = torch.zeros(1, 100, 3)
    >>> valid_prediction_time(pred, truth)  # identical -> full horizon
    100.0
    """
    if metric != "nrmse":
        raise ValueError(f"Unknown metric {metric!r}; only 'nrmse' is implemented.")
    if units not in ("steps", "lyapunov"):
        raise ValueError(f"Unknown units {units!r}; expected 'steps' or 'lyapunov'.")

    p = pred.detach().cpu().double().numpy()
    t = truth.detach().cpu().double().numpy()
    # Compare over the overlapping prefix if the two differ in length.
    n = min(p.shape[1], t.shape[1])
    p, t = p[:, :n], t[:, :n]

    scale = t.std(axis=(0, 1), keepdims=True)
    scale = np.where(scale == 0.0, 1.0, scale)
    err = np.sqrt((((p - t) / scale) ** 2).mean(axis=2))  # (batch, time)
    e_t = np.median(err, axis=0)  # (time,)

    below = e_t < threshold
    if below.size == 0 or not below[0]:
        steps = 0
    elif (~below).any():
        steps = int(np.argmax(~below))
    else:
        steps = int(below.size)

    if units == "steps":
        return float(steps)

    dt_r, lam_r = _infer_dt_lambda(dataset, dt, lambda_max, **kw)
    return steps_to_lyapunov(steps, dt_r, lam_r)

Lyapunov-time conversions¶

metrics ¶

Reservoir-computing metrics¶

Forecast-quality and dynamical metrics for reservoir models, with first-class support for reporting horizons in Lyapunov times — the natural clock of a chaotic attractor — instead of raw integration steps.

Lyapunov bookkeeping

lyapunov_time, steps_per_lyapunov_time, steps_to_lyapunov, lyapunov_to_steps Pure, dependency-free conversions between step counts and Lyapunov times. CANONICAL_LAMBDA_MAX, MACKEY_GLASS_LAMBDA Literature largest-Lyapunov-exponent tables for resdag's built-in generators. largest_lyapunov Estimate the largest Lyapunov exponent of a generator name (canonical), a measured series, or a tsdynamics System. reservoir_lyapunov Native-torch Benettin/QR estimate of a trained reservoir's own maximal Lyapunov exponent (an edge-of-chaos / Echo-State-Property diagnostic).

Forecast accuracy

valid_prediction_time Contiguous valid forecast horizon (NRMSE below a threshold), reportable in steps or Lyapunov times.

Optional dependency

The series- and System-based Lyapunov estimators require the optional tsdynamics package (pip install resdag[dynamics], Python >= 3.12); every tsdynamics import is lazy and guarded. Everything else — the conversions, the canonical constants, the reservoir Jacobian estimate, and the forecast metrics — is pure PyTorch/NumPy and always available.

Examples:

>>> from resdag.metrics import steps_to_lyapunov, valid_prediction_time
>>> round(steps_to_lyapunov(357, dt=0.02, lambda_max=0.9056), 2)
6.47

lyapunov_time ¶

lyapunov_time(lambda_max: float) -> float

Return the Lyapunov time :math:t_\lambda = 1 / \lambda_\mathrm{max}.

One Lyapunov time is the characteristic timescale over which infinitesimally close trajectories of a chaotic system separate by a factor of :math:e.

PARAMETER	DESCRIPTION
`lambda_max`	Largest Lyapunov exponent (must be strictly positive for a chaotic system). TYPE: `float`

RETURNS	DESCRIPTION
`float`	The Lyapunov time :math:`1 / \lambda_\mathrm{max}`, in the same time units as `lambda_max` is inverse to.

RAISES	DESCRIPTION
`ValueError`	If `lambda_max` is not strictly positive.

Source code in src/resdag/metrics/lyapunov.py

def lyapunov_time(lambda_max: float) -> float:
    """Return the Lyapunov time :math:`t_\\lambda = 1 / \\lambda_\\mathrm{max}`.

    One Lyapunov time is the characteristic timescale over which infinitesimally
    close trajectories of a chaotic system separate by a factor of :math:`e`.

    Parameters
    ----------
    lambda_max : float
        Largest Lyapunov exponent (must be strictly positive for a chaotic
        system).

    Returns
    -------
    float
        The Lyapunov time :math:`1 / \\lambda_\\mathrm{max}`, in the same time
        units as ``lambda_max`` is inverse to.

    Raises
    ------
    ValueError
        If ``lambda_max`` is not strictly positive.
    """
    if lambda_max <= 0.0:
        raise ValueError(f"lambda_max must be positive to define a Lyapunov time, got {lambda_max}")
    return 1.0 / lambda_max

steps_per_lyapunov_time ¶

steps_per_lyapunov_time(dt: float, lambda_max: float) -> float

Return the number of integration steps in one Lyapunov time.

A convenience for turning a horizon-in-steps into a horizon-in-Lyapunov-times (divide the step count by this) and vice versa.

PARAMETER	DESCRIPTION
`dt`	Integration step (time between consecutive samples). TYPE: `float`
`lambda_max`	Largest Lyapunov exponent (strictly positive). TYPE: `float`

RETURNS	DESCRIPTION
`float`	`1 / (dt * lambda_max)` — the number of samples spanning one Lyapunov time :math:`t_\lambda`.

RAISES	DESCRIPTION
`ValueError`	If `dt` or `lambda_max` is not strictly positive.

Source code in src/resdag/metrics/lyapunov.py

def steps_per_lyapunov_time(dt: float, lambda_max: float) -> float:
    """Return the number of integration steps in one Lyapunov time.

    A convenience for turning a horizon-in-steps into a horizon-in-Lyapunov-times
    (divide the step count by this) and vice versa.

    Parameters
    ----------
    dt : float
        Integration step (time between consecutive samples).
    lambda_max : float
        Largest Lyapunov exponent (strictly positive).

    Returns
    -------
    float
        ``1 / (dt * lambda_max)`` — the number of samples spanning one Lyapunov
        time :math:`t_\\lambda`.

    Raises
    ------
    ValueError
        If ``dt`` or ``lambda_max`` is not strictly positive.
    """
    if dt <= 0.0:
        raise ValueError(f"dt must be positive, got {dt}")
    if lambda_max <= 0.0:
        raise ValueError(f"lambda_max must be positive, got {lambda_max}")
    return 1.0 / (dt * lambda_max)

steps_to_lyapunov ¶

steps_to_lyapunov(steps: float, dt: float, lambda_max: float) -> float

Convert a step count to Lyapunov times :math:\Lambda t.

This is the headline conversion: VPT_in_Lyapunov = steps * dt * lambda_max.

PARAMETER	DESCRIPTION
`steps`	Number of integration steps (e.g. a valid-prediction time). TYPE: `float`
`dt`	Integration step. TYPE: `float`
`lambda_max`	Largest Lyapunov exponent (strictly positive). TYPE: `float`

RETURNS	DESCRIPTION
`float`	The horizon expressed in Lyapunov times, `steps * dt * lambda_max`.

RAISES	DESCRIPTION
`ValueError`	If `dt` or `lambda_max` is not strictly positive.

Examples:

A Lorenz-63 valid-prediction time of 357 steps at dt = 0.02:

>>> round(steps_to_lyapunov(357, dt=0.02, lambda_max=0.9056), 2)
6.47

Source code in src/resdag/metrics/lyapunov.py

def steps_to_lyapunov(steps: float, dt: float, lambda_max: float) -> float:
    """Convert a step count to Lyapunov times :math:`\\Lambda t`.

    This is the headline conversion: ``VPT_in_Lyapunov = steps * dt *
    lambda_max``.

    Parameters
    ----------
    steps : float
        Number of integration steps (e.g. a valid-prediction time).
    dt : float
        Integration step.
    lambda_max : float
        Largest Lyapunov exponent (strictly positive).

    Returns
    -------
    float
        The horizon expressed in Lyapunov times, ``steps * dt * lambda_max``.

    Raises
    ------
    ValueError
        If ``dt`` or ``lambda_max`` is not strictly positive.

    Examples
    --------
    A Lorenz-63 valid-prediction time of 357 steps at ``dt = 0.02``:

    >>> round(steps_to_lyapunov(357, dt=0.02, lambda_max=0.9056), 2)
    6.47
    """
    if dt <= 0.0:
        raise ValueError(f"dt must be positive, got {dt}")
    if lambda_max <= 0.0:
        raise ValueError(f"lambda_max must be positive, got {lambda_max}")
    return steps * dt * lambda_max

lyapunov_to_steps ¶

lyapunov_to_steps(lyap: float, dt: float, lambda_max: float) -> int

Convert a horizon in Lyapunov times back to a whole number of steps.

The inverse of :func:steps_to_lyapunov, rounded to the nearest integer step (horizons are discrete step counts).

PARAMETER	DESCRIPTION
`lyap`	Horizon in Lyapunov times :math:`\Lambda t`. TYPE: `float`
`dt`	Integration step. TYPE: `float`
`lambda_max`	Largest Lyapunov exponent (strictly positive). TYPE: `float`

RETURNS	DESCRIPTION
`int`	`round(lyap / (dt * lambda_max))` — the equivalent step count.

RAISES	DESCRIPTION
`ValueError`	If `dt` or `lambda_max` is not strictly positive.

Source code in src/resdag/metrics/lyapunov.py

def lyapunov_to_steps(lyap: float, dt: float, lambda_max: float) -> int:
    """Convert a horizon in Lyapunov times back to a whole number of steps.

    The inverse of :func:`steps_to_lyapunov`, rounded to the nearest integer
    step (horizons are discrete step counts).

    Parameters
    ----------
    lyap : float
        Horizon in Lyapunov times :math:`\\Lambda t`.
    dt : float
        Integration step.
    lambda_max : float
        Largest Lyapunov exponent (strictly positive).

    Returns
    -------
    int
        ``round(lyap / (dt * lambda_max))`` — the equivalent step count.

    Raises
    ------
    ValueError
        If ``dt`` or ``lambda_max`` is not strictly positive.
    """
    if dt <= 0.0:
        raise ValueError(f"dt must be positive, got {dt}")
    if lambda_max <= 0.0:
        raise ValueError(f"lambda_max must be positive, got {lambda_max}")
    return int(round(lyap / (dt * lambda_max)))

Canonical exponents¶

metrics ¶

Reservoir-computing metrics¶

Forecast-quality and dynamical metrics for reservoir models, with first-class support for reporting horizons in Lyapunov times — the natural clock of a chaotic attractor — instead of raw integration steps.

Lyapunov bookkeeping

lyapunov_time, steps_per_lyapunov_time, steps_to_lyapunov, lyapunov_to_steps Pure, dependency-free conversions between step counts and Lyapunov times. CANONICAL_LAMBDA_MAX, MACKEY_GLASS_LAMBDA Literature largest-Lyapunov-exponent tables for resdag's built-in generators. largest_lyapunov Estimate the largest Lyapunov exponent of a generator name (canonical), a measured series, or a tsdynamics System. reservoir_lyapunov Native-torch Benettin/QR estimate of a trained reservoir's own maximal Lyapunov exponent (an edge-of-chaos / Echo-State-Property diagnostic).

Forecast accuracy

valid_prediction_time Contiguous valid forecast horizon (NRMSE below a threshold), reportable in steps or Lyapunov times.

Optional dependency

The series- and System-based Lyapunov estimators require the optional tsdynamics package (pip install resdag[dynamics], Python >= 3.12); every tsdynamics import is lazy and guarded. Everything else — the conversions, the canonical constants, the reservoir Jacobian estimate, and the forecast metrics — is pure PyTorch/NumPy and always available.

Examples:

>>> from resdag.metrics import steps_to_lyapunov, valid_prediction_time
>>> round(steps_to_lyapunov(357, dt=0.02, lambda_max=0.9056), 2)
6.47

CANONICAL_LAMBDA_MAX `module-attribute` ¶

CANONICAL_LAMBDA_MAX: dict[str, float] = {'lorenz': 0.9056, 'rossler': 0.0714}

MACKEY_GLASS_LAMBDA `module-attribute` ¶

MACKEY_GLASS_LAMBDA: dict[int, float] = {17: 0.0086, 23: 0.0097, 30: 0.0071}

Estimators¶

metrics ¶

Reservoir-computing metrics¶

Forecast-quality and dynamical metrics for reservoir models, with first-class support for reporting horizons in Lyapunov times — the natural clock of a chaotic attractor — instead of raw integration steps.

Lyapunov bookkeeping

lyapunov_time, steps_per_lyapunov_time, steps_to_lyapunov, lyapunov_to_steps Pure, dependency-free conversions between step counts and Lyapunov times. CANONICAL_LAMBDA_MAX, MACKEY_GLASS_LAMBDA Literature largest-Lyapunov-exponent tables for resdag's built-in generators. largest_lyapunov Estimate the largest Lyapunov exponent of a generator name (canonical), a measured series, or a tsdynamics System. reservoir_lyapunov Native-torch Benettin/QR estimate of a trained reservoir's own maximal Lyapunov exponent (an edge-of-chaos / Echo-State-Property diagnostic).

Forecast accuracy

valid_prediction_time Contiguous valid forecast horizon (NRMSE below a threshold), reportable in steps or Lyapunov times.

Optional dependency

The series- and System-based Lyapunov estimators require the optional tsdynamics package (pip install resdag[dynamics], Python >= 3.12); every tsdynamics import is lazy and guarded. Everything else — the conversions, the canonical constants, the reservoir Jacobian estimate, and the forecast metrics — is pure PyTorch/NumPy and always available.

Examples:

>>> from resdag.metrics import steps_to_lyapunov, valid_prediction_time
>>> round(steps_to_lyapunov(357, dt=0.02, lambda_max=0.9056), 2)
6.47

largest_lyapunov ¶

largest_lyapunov(source: Any, *, dt: float | None = None, from_: str | None = None, **kw: Any) -> float

Estimate the largest Lyapunov exponent of source.

A single entry point that dispatches on the type of source:

generator name (:class:str, e.g. "lorenz", "mackey_glass") — returns the canonical literature value from :data:CANONICAL_LAMBDA_MAX / :data:MACKEY_GLASS_LAMBDA (keyed by tau= for Mackey-Glass). No dependency.
time series (:class:numpy.ndarray / :class:torch.Tensor) — a diagnostic estimate via tsdynamics.lyapunov_from_data. Data-driven exponents are noisy; you must opt in with from_="data" so a series is never silently mistaken for accurate. Needs resdag[dynamics].
tsdynamics System — a high-accuracy Benettin estimate via tsdynamics.lyapunov_spectrum. Needs resdag[dynamics].

PARAMETER	DESCRIPTION
`source`	What to estimate the exponent of; see above. TYPE: `str, numpy.ndarray, torch.Tensor, or tsdynamics System`
`dt`	Integration step. Required for the series path; forwarded to the System path if given. Ignored for a generator name. TYPE: `float` DEFAULT: `None`
`from_`	Set to `"data"` to force the noisy data-driven estimate for an array source. Without it, passing a raw array raises to prevent silently trusting a diagnostic number. TYPE: `('data', None)` DEFAULT: `"data"`
`**kw`	Extra arguments: `tau=` selects the Mackey-Glass canonical value; otherwise forwarded to the underlying `tsdynamics` estimator (e.g. `method=`, `fit=` for data; `final_time=`, `transient=` for a System). TYPE: `Any` DEFAULT: `{}`

RETURNS	DESCRIPTION
`float`	The estimated largest Lyapunov exponent (inverse natural-time units).

RAISES	DESCRIPTION
`ImportError`	If a `tsdynamics` path is taken but the package is not installed.
`TypeError`	If an array is passed without `from_="data"`.
`KeyError`	If a generator name has no canonical value.

Examples:

>>> largest_lyapunov("lorenz")
0.9056
>>> largest_lyapunov("mackey_glass", tau=17)
0.0086

Source code in src/resdag/metrics/lyapunov.py

def largest_lyapunov(
    source: Any,
    *,
    dt: float | None = None,
    from_: str | None = None,
    **kw: Any,
) -> float:
    """Estimate the largest Lyapunov exponent of ``source``.

    A single entry point that dispatches on the type of ``source``:

    * **generator name** (:class:`str`, e.g. ``"lorenz"``, ``"mackey_glass"``) —
      returns the canonical literature value from :data:`CANONICAL_LAMBDA_MAX`
      / :data:`MACKEY_GLASS_LAMBDA` (keyed by ``tau=`` for Mackey-Glass).
      **No dependency.**
    * **time series** (:class:`numpy.ndarray` / :class:`torch.Tensor`) — a
      *diagnostic* estimate via ``tsdynamics.lyapunov_from_data``.  Data-driven
      exponents are noisy; you must opt in with ``from_="data"`` so a series is
      never silently mistaken for accurate.  **Needs** ``resdag[dynamics]``.
    * **tsdynamics System** — a high-accuracy Benettin estimate via
      ``tsdynamics.lyapunov_spectrum``.  **Needs** ``resdag[dynamics]``.

    Parameters
    ----------
    source : str, numpy.ndarray, torch.Tensor, or tsdynamics System
        What to estimate the exponent of; see above.
    dt : float, optional
        Integration step.  Required for the series path; forwarded to the
        System path if given.  Ignored for a generator name.
    from_ : {"data", None}, optional
        Set to ``"data"`` to force the noisy data-driven estimate for an array
        source.  Without it, passing a raw array raises to prevent silently
        trusting a diagnostic number.
    **kw
        Extra arguments: ``tau=`` selects the Mackey-Glass canonical value;
        otherwise forwarded to the underlying `tsdynamics` estimator (e.g.
        ``method=``, ``fit=`` for data; ``final_time=``, ``transient=`` for a
        System).

    Returns
    -------
    float
        The estimated largest Lyapunov exponent (inverse natural-time units).

    Raises
    ------
    ImportError
        If a `tsdynamics` path is taken but the package is not installed.
    TypeError
        If an array is passed without ``from_="data"``.
    KeyError
        If a generator name has no canonical value.

    Examples
    --------
    >>> largest_lyapunov("lorenz")
    0.9056
    >>> largest_lyapunov("mackey_glass", tau=17)
    0.0086
    """
    if isinstance(source, str):
        return _lambda_from_name(source, **kw)

    if isinstance(source, torch.Tensor):
        source = source.detach().cpu().numpy()

    if isinstance(source, np.ndarray):
        if from_ != "data":
            raise TypeError(
                "Estimating a Lyapunov exponent from a raw data series is a noisy "
                "diagnostic; pass from_='data' (and dt=...) to opt in explicitly, "
                "or use a generator name / tsdynamics System for an accurate value."
            )
        if dt is None:
            raise ValueError("dt is required to estimate a Lyapunov exponent from data.")
        return _lambda_from_series(source, dt=dt, **kw)

    # Anything else is assumed to be a tsdynamics System (dispatched lazily so we
    # never import tsdynamics just to check the type).
    return _lambda_from_system(source, dt=dt, **kw)

reservoir_lyapunov ¶

reservoir_lyapunov(source: Any, warmup_inputs: Tensor | tuple[Tensor, ...], *, horizon: int = 1000, k: int = 1, transient: int = 100, seed: int | None = 0) -> float | ndarray

Largest (or top-k) Lyapunov exponent of a driven ESN reservoir.

Estimates the maximal Lyapunov exponent of the reservoir's own dynamics via Benettin's algorithm with QR reorthonormalisation of the analytic tangent map, iterated along an input-driven reservoir trajectory. This is a standard reservoir-computing quality / edge-of-chaos diagnostic: a healthy ESN sits just below the edge, with a maximal exponent that is negative but close to zero (the Echo State Property requires the driven dynamics to contract, so trajectories from different initial states converge).

The reservoir map is the leaky-ESN update

.. math::

h_t = (1-\alpha) h_{t-1} + \alpha\, f(W_\mathrm{fb} x_t
      + W_\mathrm{in} u_t + W_\mathrm{hh} h_{t-1} + b),

whose Jacobian with respect to the state is

.. math::

J_t = (1-\alpha) I + \alpha\, \mathrm{diag}\!\big(f'(\mathrm{pre}_t)\big)\, W_\mathrm{hh},

with :math:f' = 1 - \tanh^2 for the (default) tanh activation. The driving inputs come from warmup_inputs; the trajectory is looped over them if horizon exceeds their length, so a short synchronisation window still yields a long tangent-space average.

PARAMETER	DESCRIPTION
`source`	The reservoir to analyse, or a model containing exactly one ESN reservoir. TYPE: `BaseReservoirLayer or ESNModel`
`warmup_inputs`	Driving inputs of shape `(batch, time, features)`. A bare tensor is the feedback stream; a tuple is `(feedback, driver, ...)`. Only the first sample of the batch is used (the exponent is a property of the dynamics, not the batch). TYPE: `torch.Tensor or tuple of torch.Tensor`
`horizon`	Number of tangent-map iterations (after the transient). Longer gives a tighter estimate. TYPE: `int` DEFAULT: `1000`
`k`	Number of leading exponents to return. `k=1` returns the scalar maximal exponent; `k>1` returns the top-`k` as a `(k,)` array. TYPE: `int` DEFAULT: `1`
`transient`	Leading steps whose tangent growth is discarded so the estimate reflects the attractor rather than the initial state. Clamped so at least a few accumulation steps remain. TYPE: `int` DEFAULT: `100`
`seed`	Seed for the random initial tangent frame. Defaults to `0` so the estimate is fully deterministic (the converged leading exponent is robust to the frame regardless). Pass `None` to draw from the process-global RNG instead. TYPE: `int or None` DEFAULT: `0`

RETURNS	DESCRIPTION
`float or ndarray`	The maximal Lyapunov exponent (`k=1`) or the top-`k` spectrum as a `(k,)` array, in inverse-step units (per reservoir update).

RAISES	DESCRIPTION
`ValueError`	If `source` has no (or more than one) ESN reservoir, or `k` exceeds the reservoir size.
`NotImplementedError`	If the reservoir activation is not `tanh` (the only activation whose analytic derivative is wired here).

Notes

Exponents are per reservoir step. Multiply by 1/dt to compare with a natural-time attractor exponent from :data:CANONICAL_LAMBDA_MAX.

Examples:

>>> import torch
>>> from resdag.layers import ESNLayer
>>> res = ESNLayer(100, feedback_size=3, spectral_radius=0.9, seed=0)
>>> feedback = torch.randn(1, 300, 3)
>>> lam = reservoir_lyapunov(res, feedback, horizon=500)
>>> lam < 0.0  # a stable, ESP-satisfying reservoir contracts
True

Source code in src/resdag/metrics/lyapunov.py

def reservoir_lyapunov(
    source: Any,
    warmup_inputs: torch.Tensor | tuple[torch.Tensor, ...],
    *,
    horizon: int = 1000,
    k: int = 1,
    transient: int = 100,
    seed: int | None = 0,
) -> float | np.ndarray:
    """Largest (or top-``k``) Lyapunov exponent of a driven ESN reservoir.

    Estimates the maximal Lyapunov exponent of the reservoir's own dynamics via
    **Benettin's algorithm with QR reorthonormalisation** of the analytic
    tangent map, iterated along an input-driven reservoir trajectory.  This is a
    standard reservoir-computing quality / *edge-of-chaos* diagnostic: a healthy
    ESN sits just below the edge, with a maximal exponent that is negative but
    close to zero (the Echo State Property requires the driven dynamics to
    contract, so trajectories from different initial states converge).

    The reservoir map is the leaky-ESN update

    .. math::

        h_t = (1-\\alpha) h_{t-1} + \\alpha\\, f(W_\\mathrm{fb} x_t
              + W_\\mathrm{in} u_t + W_\\mathrm{hh} h_{t-1} + b),

    whose Jacobian with respect to the state is

    .. math::

        J_t = (1-\\alpha) I + \\alpha\\, \\mathrm{diag}\\!\\big(f'(\\mathrm{pre}_t)\\big)\\, W_\\mathrm{hh},

    with :math:`f' = 1 - \\tanh^2` for the (default) ``tanh`` activation.  The
    driving inputs come from ``warmup_inputs``; the trajectory is looped over
    them if ``horizon`` exceeds their length, so a short synchronisation window
    still yields a long tangent-space average.

    Parameters
    ----------
    source : BaseReservoirLayer or ESNModel
        The reservoir to analyse, or a model containing exactly one ESN
        reservoir.
    warmup_inputs : torch.Tensor or tuple of torch.Tensor
        Driving inputs of shape ``(batch, time, features)``.  A bare tensor is
        the feedback stream; a tuple is ``(feedback, driver, ...)``.  Only the
        first sample of the batch is used (the exponent is a property of the
        dynamics, not the batch).
    horizon : int, default=1000
        Number of tangent-map iterations (after the transient).  Longer gives a
        tighter estimate.
    k : int, default=1
        Number of leading exponents to return.  ``k=1`` returns the scalar
        maximal exponent; ``k>1`` returns the top-``k`` as a ``(k,)`` array.
    transient : int, default=100
        Leading steps whose tangent growth is discarded so the estimate reflects
        the attractor rather than the initial state.  Clamped so at least a few
        accumulation steps remain.
    seed : int or None, default=0
        Seed for the random initial tangent frame.  Defaults to ``0`` so the
        estimate is fully deterministic (the converged leading exponent is robust
        to the frame regardless).  Pass ``None`` to draw from the process-global
        RNG instead.

    Returns
    -------
    float or numpy.ndarray
        The maximal Lyapunov exponent (``k=1``) or the top-``k`` spectrum as a
        ``(k,)`` array, in **inverse-step units** (per reservoir update).

    Raises
    ------
    ValueError
        If ``source`` has no (or more than one) ESN reservoir, or ``k`` exceeds
        the reservoir size.
    NotImplementedError
        If the reservoir activation is not ``tanh`` (the only activation whose
        analytic derivative is wired here).

    Notes
    -----
    Exponents are per **reservoir step**.  Multiply by ``1/dt`` to compare with
    a natural-time attractor exponent from :data:`CANONICAL_LAMBDA_MAX`.

    Examples
    --------
    >>> import torch
    >>> from resdag.layers import ESNLayer
    >>> res = ESNLayer(100, feedback_size=3, spectral_radius=0.9, seed=0)
    >>> feedback = torch.randn(1, 300, 3)
    >>> lam = reservoir_lyapunov(res, feedback, horizon=500)
    >>> lam < 0.0  # a stable, ESP-satisfying reservoir contracts
    True
    """
    from resdag.layers.cells import ESNCell

    reservoir = _resolve_reservoir(source)
    # ``_resolve_reservoir`` guarantees a cell that owns ``weight_hh`` (an ESN
    # cell); narrow the type so the analytic-Jacobian access below is checked.
    cell = cast("ESNCell", reservoir.cell)

    activation = getattr(cell, "activation", "tanh")
    if activation != "tanh":
        raise NotImplementedError(
            f"reservoir_lyapunov only supports the 'tanh' activation "
            f"(analytic derivative 1 - tanh^2); got {activation!r}."
        )

    inputs: tuple[torch.Tensor, ...]
    if isinstance(warmup_inputs, torch.Tensor):
        inputs = (warmup_inputs,)
    else:
        inputs = tuple(warmup_inputs)
    if not inputs:
        raise ValueError("warmup_inputs must contain at least the feedback stream.")

    weight_hh: torch.Tensor = cell.weight_hh
    native_dtype: torch.dtype = weight_hh.dtype
    device = weight_hh.device
    # Work in float64 for a numerically stable QR accumulation regardless of the
    # reservoir's own dtype.
    W_hh = weight_hh.detach().to(torch.float64)
    N = int(W_hh.shape[0])

    if not 1 <= k <= N:
        raise ValueError(f"k must be in [1, reservoir_size={N}], got {k}.")

    alpha = float(getattr(cell, "leak_rate", 1.0))

    # Use the first sample of the batch; keep a leading singleton batch axis so
    # the cell's single-step API (project_inputs / forward) accepts it.
    per_step_inputs = [x[:1].to(device=device, dtype=native_dtype) for x in inputs]
    seq_len = per_step_inputs[0].shape[1]
    total_steps = transient + horizon

    # Synchronise the reservoir state on the driving inputs, then read it back so
    # the tangent iteration starts on the attractor.  Route through the public
    # forward so state layout / validation is handled by the layer.
    reservoir.reset_state(batch_size=1)
    with torch.no_grad():
        reservoir(*per_step_inputs)
    state = reservoir.get_state()
    if state is None:  # pragma: no cover - forward just materialised the state
        raise RuntimeError("Reservoir state was not materialised by the warmup pass.")
    h = state.to(torch.float64).reshape(N)  # (N,)

    # Orthonormal tangent frame (N, k); columns tracked through the flow.  The
    # converged leading exponent is robust to the initial frame, but seeding it
    # (default ``seed=0``) makes the whole estimate deterministic instead of
    # dependent on the process-global RNG.  ``seed=None`` draws from the global
    # RNG (non-deterministic).
    generator: torch.Generator | None = None
    if seed is not None:
        generator = torch.Generator(device=device)
        generator.manual_seed(int(seed))
    Q = torch.linalg.qr(torch.randn(N, k, dtype=torch.float64, device=device, generator=generator))[
        0
    ]
    log_sums = torch.zeros(k, dtype=torch.float64, device=device)
    n_accum = 0

    identity = torch.eye(N, dtype=torch.float64, device=device)

    with torch.no_grad():
        for step in range(total_steps):
            # Cycle through the driving window if the horizon outruns it.
            t = step % seq_len
            slices = [x[:, t, :] for x in per_step_inputs]  # each (1, features)
            # Pre-activation in the cell's native dtype, then promote.
            projected = cell.project_inputs(slices).to(torch.float64).reshape(N)  # (N,)
            pre = projected + W_hh @ h  # (N,)
            new_h = torch.tanh(pre)
            if alpha != 1.0:
                new_h = (1.0 - alpha) * h + alpha * new_h

            # Analytic Jacobian J = (1-a) I + a diag(f'(pre)) W_hh.
            fprime = 1.0 - torch.tanh(pre) ** 2  # (N,)
            J = alpha * (fprime[:, None] * W_hh)
            if alpha != 1.0:
                J = J + (1.0 - alpha) * identity

            # Evolve the tangent frame and reorthonormalise (Benettin/QR).
            Z = J @ Q  # (N, k)
            Q, R = torch.linalg.qr(Z)
            diag_r = torch.abs(torch.diagonal(R))
            diag_r = torch.clamp(diag_r, min=1e-300)  # guard log(0)
            if step >= transient:
                log_sums += torch.log(diag_r)
                n_accum += 1

            h = new_h

    if n_accum == 0:  # pragma: no cover - horizon >= 1 guarantees accumulation
        raise ValueError("No accumulation steps; increase horizon relative to transient.")

    exponents = (log_sums / n_accum).cpu().numpy()  # per-step exponents
    if k == 1:
        return float(exponents[0])
    return exponents

Metrics¶

Forecast accuracy¶

metrics ¶

Reservoir-computing metrics¶

valid_prediction_time ¶

Lyapunov-time conversions¶

metrics ¶

Reservoir-computing metrics¶

lyapunov_time ¶

steps_per_lyapunov_time ¶

steps_to_lyapunov ¶

lyapunov_to_steps ¶

Canonical exponents¶

metrics ¶

Reservoir-computing metrics¶

CANONICAL_LAMBDA_MAX module-attribute ¶

MACKEY_GLASS_LAMBDA module-attribute ¶

Estimators¶

metrics ¶

Reservoir-computing metrics¶

largest_lyapunov ¶

reservoir_lyapunov ¶

CANONICAL_LAMBDA_MAX `module-attribute` ¶

MACKEY_GLASS_LAMBDA `module-attribute` ¶