dysts API Reference¶

The API reference for the dysts repository: https://github.com/williamgilpin/dysts

Base Classes¶

Dynamical systems in Python

class dysts.base.BaseDyn(metadata_path: str | None = None, metadata: dict[str, Any] | None = None, required_fields: tuple[str | tuple[str, ...], ...] = ('parameters', ('dimension', 'embedding_dimension')), **extra_metadata)¶

A base class for dynamical systems

has_jacobian() → bool¶: Check if the subclass has implemented the _jac method.

static load_system_metadata(system_name: str, data_path: str) → dict[str, Any]¶

Load data from a JSON file

Returns None if the system name is not found

static load_trajectory(data_path: str, system_name: str, return_times=False, standardize=False)¶

Load a precomputed trajectory for the dynamical system

Parameters:

data_path (str) – Path to the data file, of format {name}.json.gz
standardize (bool) – Standardize the output time series.
return_times (bool) – Whether to return the timepoints at which the solution was computed

Returns:

A T x D trajectory tpts, sol (ndarray): T x 1 timepoint array, and T x D trajectory

Return type:

sol (ndarray)

make_trajectory(*args, **kwargs)¶: Make a trajectory for the dynamical system

transform_ic(transform_fn: Callable[[ndarray, Any], ndarray | None]) → bool¶: Updates the initial condition via a transform function

transform_params(transform_fn: Callable[[str, ndarray, Any], ndarray | None]) → bool¶: Updates the current parameter list via a transform function

A discrete map dynamical system class

make_trajectory(n: int, init_cond: ndarray | None = None, inverse: bool = False, return_times: bool = False, standardize: bool = False, **kwargs) → ndarray | tuple[ndarray, ndarray]¶

Generate a fixed-length trajectory with default parameters and initial condition(s).

Parameters:

n (int) – The length of each trajectory.
init_cond (Optional[np.ndarray]) – Initial conditions. If None, use default.
inverse (bool) – Whether to reverse the trajectory.
return_times (bool) – Whether to return the timepoints of the solution.
standardize (bool) – Whether to standardize the output time series.

Returns:

If return_times is False, returns the trajectory. If return_times is True, returns a tuple of (timepoints, trajectory).

Return type:

Union[np.ndarray, Tuple[np.ndarray, np.ndarray]]

rhs(X)¶: The right hand side of a dynamical map

rhs_inv(Xp)¶: The inverse of the right hand side of a dynamical map

A continuous dynamical system base class

jac(X, t)¶: The Jacobian of the dynamical system

make_trajectory(n: int, dt: float = 0.001, init_cond: ndarray | None = None, resample: bool = True, pts_per_period: int = 100, return_times: bool = False, standardize: bool = False, postprocess: bool = True, noise: float = 0.0, timescale: str = 'Fourier', method: str = 'Radau', random_seed: int = 0, rtol: float = 1e-12, atol: float = 1e-12, verbose: bool = False, **kwargs) → ndarray | tuple[ndarray, ndarray] | None | tuple[None, None]¶

Generate a fixed-length trajectory for the dynamical system.

Parameters:

n – Total number of trajectory points.
dt – Time step for integration. Defaults to 1e-3 or system’s default if set.
init_cond – Initial conditions. If None, uses system’s default.
resample – Whether to resample trajectory to match dominant Fourier components.
pts_per_period – Number of points per period if resampling.
return_times – If True, return time points along with trajectory.
standardize – If True, standardize the output time series.
postprocess – If True, apply coordinate conversions and rescalings.
noise – Stochasticity level in integrated dynamics (corresponds to Brownian motion).
timescale – Timescale for resampling. “Fourier” (default) or “Lyapunov”.
method – Integration method.
random_seed – Seed for random number generation.
rtol – Relative tolerance for integration.
atol – Absolute tolerance for integration.
verbose – Whether to trigger warnings.
**kwargs – Additional arguments for integration routine.

Returns:

np.ndarray: T x D trajectory array. If return_times is True:

Tuple[np.ndarray, np.ndarray]: T x 1 time point array and T x D trajectory array.

None: If no complete trajectories are found.

Return type:

If return_times is False

Raises:

ValueError – If an invalid timescale is provided.

rhs(X, t)¶: The right hand side of a dynamical equation

A delayed differential equation object. Uses a exposed fork of ddeint The delay timescale is assumed to be the “tau” field. The embedding dimensions are set to a default value, but delay equations are infinite dimensional.

delayed_rhs(X: Callable[[float], Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]], t: float, tau: float) → Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | bool | int | float | complex | str | bytes | _NestedSequence[bool | int | float | complex | str | bytes]¶: The right hand side of a dynamical equation

Generate a fixed-length trajectory with default timestep, parameters, and initial conditions.

Parameters:

n (int) – the total number of trajectory points
method (str) – Not used. Currently Euler is the only option here
noise (float) – The amplitude of brownian forcing
resample (bool) – whether to resample trajectories to have matching dominant Fourier components
pts_per_period (int) – if resampling, the number of points per period
standardize (bool) – Standardize the output time series.
timescale (str) – The timescale to use for resampling. “Fourier” (default) uses the dominant significant Fourier timescale, estimated using the periodogram of the system and surrogates. “Lyapunov” uses the Lypunov timescale of the system.
return_times (bool) – Whether to return the timepoints at which the solution was computed
postprocess (bool) – Whether to apply coordinate conversions and other domain-specific rescalings to the integration coordinates
embedding_dim (int) – Optionally augment solution with delay embedded trajectories. If equation is multi-dimensional with dimension d, then this will return a flattened delay embedding of shape (n, d*embedding_dim) where each consecutive (n, d) block will be a trajectory for a different delay parameter.
history_function (callable) – Function for specifying past conditions i.e. points for t < 0

dysts.base.staticjit(func: Callable) → Callable¶: Decorator to apply numba’s njit decorator to a static method

Utilities¶

Helper utilities for working with time series arrays. This module is intended to have no dependencies on the rest of the package.

dysts.utils.utils.cartesian_to_polar(x, y)¶: Convert 2D cartesian coordinates to polar coordinates

dysts.utils.utils.find_characteristic_timescale(y, k=1, window=True)¶: Find the k leading characteristic timescales in a time series using the power spectrum..

dysts.utils.utils.find_psd(y, window=True)¶: Find the power spectrum of a signal

dysts.utils.utils.find_significant_frequencies(sig, window=True, fs=1, n_samples=100, significance_threshold=0.95, surrogate_method='rp', show=False, return_amplitudes=False)¶

Find power spectral frequencies that are significant in a signal, by comparing the appearance of a peak with its appearance in randomly-shuffled surrogates.

If no significant freqencies are detected, the significance floor is lowered

Parameters:

window (bool) – Whether to window the signal before taking the FFT
thresh (float) – The number of standard deviations above mean to be significant
fs (int) – the sampling frequency
n_samples (int) – the number of surrogates to create
show (bool) – whether to show the psd of the signal and the surrogate

Returns:

The frequencies overrated in the dataset amps (ndarray): the amplitudes of the PSD at the identified frequencies

Return type:

freqs (ndarray)

dysts.utils.utils.find_slope(x, y)¶

Given two vectors or arrays, compute the best fit slope using an analytic formula. For arrays is computed along the last axis.

Parameters:

x (ndarray) – (N,) or (M, N)
y (ndarray) – (N,) or (M, N)

Returns:

the values of the slope for each of the last dimensions

Return type:

b (ndarray)

dysts.utils.utils.freq_from_autocorr(sig, fs=1)¶

Estimate frequency using autocorrelation

Parameters:

sig (ndarray) – A univariate signal
fs (int) – The sampling frequency

Returns:

The dominant frequency

Return type:

out (float)

References

Modified from the following https://gist.github.com/endolith/255291

dysts.utils.utils.freq_from_fft(sig, fs=1)¶

Estimate frequency of a signal from the peak of the power spectrum

Parameters:

sig (ndarray) – A univariate signal
fs (int) – The sampling frequency

Returns:

The dominant frequency

Return type:

out (float)

References

Modified from the following https://gist.github.com/endolith/255291

dysts.utils.utils.jac_fd(func0, y0, eps=0.001, m=1, method='central', verbose=False)¶

Calculate numerical jacobian of a function with respect to a reference value

Parameters:

func (callable) – a vector-valued function
y0 (ndarray) – a point around which to take the gradient
eps (float) – the step size for the finite difference calculation

Returns:

a numerical estimate of the Jacobian about that point

Return type:

jac (ndarray)

dysts.utils.utils.logarithmic_n(min_n, max_n, factor)¶

Return a list of values by successively multiplying a minimum value min_n by a factor > 1 until a maximum value max_n is reached. Non-integer results are rounded down. Based on a similar function in the nolds Python library.

Parameters:

min_n (float) – minimum value (must be < max_n)
max_n (float) – maximum value (must be > min_n)
factor (float) – factor used to increase min_n (must be > 1)

Returns:

min_n, min_n * factor, min_n * factor^2, … min_n * factor^i < max_n without duplicates

Return type:

list of integers

dysts.utils.utils.make_epsilon_ball(pt, n, eps=1e-05, random_state=None)¶

Uniformly sample a fixed-radius ball of points around a given point via using Muller’s method

Parameters:

pt (ndarray) – The center of the sampling
n (int) – The number of points to sample
eps (float) – The radius of the ball
random_state (int) – Initialize the random number generator

Returns:

The set of randomly-sampled points

Return type:

out (ndarray)

dysts.utils.utils.make_surrogate(data, method='rp')¶

Parameters:

data (ndarray) – A one-dimensional time series
method (str) – “rs” or rp”

Returns:

A single random surrogate time series

Return type:

surr_data (ndarray)

dysts.utils.utils.min_data_points_rosenstein(emb_dim, lag, trajectory_len, min_tsep)¶

Adapted from the nolds Python library: https://github.com/CSchoel/nolds/blob/master/nolds/measures.py

Helper function that calculates the minimum number of data points required to use lyap_r. Note that none of the required parameters may be set to None. :param emb_dim: Embedding dimension :type emb_dim: int :param lag: Lag between time series values :type lag: int :param trajectory_len: Length of trajectory to follow :type trajectory_len: int :param min_tsep: Minimum temporal separation :type min_tsep: int

Returns:: minimum number of data points required to call lyap_r with the given parameters
Return type:: int

dysts.utils.utils.nan_fill(a, axis=0, mode='forward')¶

Fill NaN values in an array using a specified method

Parameters:

a (np.ndarray) – the array to fill
axis (int) – the axis along which to fill NaN values
mode (str) – the method to use for filling NaN values

Returns:

the array with NaN values filled

Return type:

a_filled (np.ndarray)

dysts.utils.utils.nanmean_trimmed(arr, percentile=0.5, axis=None)¶

Compute the trimmed mean of an array along a specified axis, ignoring NaNs.

Parameters:

arr (np.array) – Input array
percentile (float) – The fraction of data to be trimmed. Should be between 0 and 1. Default is 0.5.
axis (int) – Axis along which to compute the trimmed mean. If None (the default), compute over the whole array.

Returns:

Trimmed mean of the input array along the specified axis.

dysts.utils.utils.pad_axis(arr, d, axis=-1, padding=0)¶: Pad axis of arr with a constant padding to a desired shape

dysts.utils.utils.pad_to_shape(arr, target_shape)¶

Given an array, and a target shape, pad the dimensions in order to reach the desired shape Currently, if the rank of the array is lower than the target shape, singleton dimensions are appended to the rank

Parameters:

arr (ndarray) – The array to pad.
target_shape (iterable) – The desired shape.

Returns:

The padded array,

Return type:

arr (ndarray)

dysts.utils.utils.polar_to_cartesian(r, th)¶: Convert polar coordinates to 2D Cartesian coordinates

dysts.utils.utils.rowwise_euclidean(x, y)¶: Computes the euclidean distance across rows

dysts.utils.utils.standardize_ts(a, scale=1.0)¶

Standardize an array along dimension -2 For dimensions with zero variance, divide by one instead of zero

Parameters:

a (ndarray) – a matrix containing a time series or batch of time series with shape (T, D) or (B, T, D)
scale (float) – the number of standard deviations by which to scale

Returns:

A standardized time series with the same shape as: the input

Return type:

ts_scaled (ndarray)

Native python utilities

class dysts.utils.native_utils.ComputationHolder(func, *args, timeout: float = 10, **kwargs)¶

A wrapper class to force a computation to stop after a timeout.

Parameters

func (callable): the function to run args (tuple): the arguments to pass to the function kwargs (dict): the keyword arguments to pass to the function timeout (int): the timeout in seconds. If None is passed, the computation

will run indefinitely until it finishes.

Example

>>> def my_func():
...     while True:
...         print("hello")
...         time.sleep(8)
>>> ch = ComputationHolder(my_func, timeout=3)
>>> ch.run()
hello
hello
hello
None

dysts.utils.native_utils.convert_json_to_gzip(fpath, encoding='utf-8', delete_original=False)¶

Convert a json file to a gzip file in a format that can be easily read by the dysts package. By default, the gzip file will be saved with the same name and in the same directory as the json file, but with a “.gz” extension.

Parameters:

fpath (str) – Path to the json file to be converted
encoding (str) – Encoding to use when writing the gzip file
delete_original (bool) – Whether to delete the original json file after conversion. Default is False.

Returns:

None

dysts.utils.native_utils.has_module(module_name: str) → bool¶: Check if a module is installed

dysts.utils.native_utils.num_unspecified_params(func) → int¶

Count required parameters that haven’t been specified through partial or defaults.

Parameters:: func – Function or partial function to inspect
Returns:: Number of parameters that must be specified when calling the function

Utilities for integration

class dysts.utils.integration_utils.dde(f, jac=None)¶

This class overwrites a few functions of scipy.integrate.ode to allow for updates of the pseudo-variable Y between each integration step.

integrate(t, step=False, relax=False)¶

Find y=y(t), set y as an initial condition, and return y.

Parameters:

t (float) – The endpoint of the integration step.
step (bool) – If True, and if the integrator supports the step method, then perform a single integration step and return. This parameter is provided in order to expose internals of the implementation, and should not be changed from its default value in most cases.
relax (bool) – If True and if the integrator supports the run_relax method, then integrate until t_1 >= t and return. relax is not referenced if step=True. This parameter is provided in order to expose internals of the implementation, and should not be changed from its default value in most cases.

Returns:

y – The integrated value at t

Return type:

float

set_initial_value(Y)¶: Set initial conditions y(t) = y.

class dysts.utils.integration_utils.ddeVar(g, tc=0)¶

The instances of this class are special function-like variables which store their past values in an interpolator and can be called for any past time: Y(t), Y(t-d). Very convenient for the integration of DDEs.

update(t, Y)¶: Add one new (ti,yi) to the interpolator

dysts.utils.integration_utils.ddeint(func, g, tt, fargs=None)¶

Solves Delay Differential Equations

Similar to scipy.integrate.odeint. Solves a Delay differential Equation system (DDE) defined by

Y(t) = g(t) for t<0 Y’(t) = func(Y,t) for t>= 0

Where func can involve past values of Y, like Y(t-d).

Parameters:

func – a function Y,t,args -> Y’(t), where args is optional. The variable Y is an instance of class ddeVar, which means that it is called like a function: Y(t), Y(t-d), etc. Y(t) returns either a number or a numpy array (for multivariate systems).
g – The ‘history function’. A function g(t)=Y(t) for t<0, g(t) returns either a number or a numpy array (for multivariate systems).
tt – The vector of times [t0, t1, …] at which the system must be solved.
fargs – Additional arguments to be passed to parameter func, if any.

Examples

We will solve the delayed Lotka-Volterra system defined as

For t < 0: x(t) = 1+t y(t) = 2-t

For t >= 0: dx/dt = 0.5* ( 1- y(t-d) ) dy/dt = -0.5* ( 1- x(t-d) )

The delay d is a tunable parameter of the model.

import numpy as np
from ddeint import ddeint

def model(XY,t,d):
    x, y = XY(t)
    xd, yd = XY(t-d)
    return np.array([0.5*x*(1-yd), -0.5*y*(1-xd)])

g = lambda t : np.array([1+t,2-t]) # 'history' at t<0
tt = np.linspace(0,30,20000) # times for integration
d = 0.5 # set parameter d
yy = ddeint(model,g,tt,fargs=(d,)) # solve the DDE !

dysts.utils.integration_utils.generate_ic_ensemble(model, tpts0, n_samples, frac_perturb_param=0.1, frac_transient=0.1, ic_range=None, random_state=0)¶

Generate an ensemble of trajectories with random initial conditions, labelled by different initial conditions

Parameters:

model (callable_) – function defining the numerical derivative
tpts (ndarray) – the timesteps over which to run the simulation
n_samples (int) – the number of different initial conditons
frac_perturb_param (float) – the amount to perturb the ic by
frac_transient (float) – the fraction of time for the time series to settle onto the attractor
ic_range (list) – a starting value for the initial conditions
random_state (int) – the seed for the random number generator

Returns:

A set of initial conditions

Return type:

all_samples (array)

dysts.utils.integration_utils.integrate_dyn(f: Callable, ic: ndarray, tvals: ndarray, noise: float | ndarray = 0.0, dtval=None, **kwargs)¶

Given the RHS of a dynamical system, integrate the system noise > 0 requires the Python library sdeint (assumes Brownian noise)

Parameters:

f (callable) – The right hand side of a system of ODEs.
ic (ndarray) – the initial conditions
tvals (ndarray) – times points at which to evaluate the solution
noise (float or iterable) – The amplitude of the Langevin forcing term. If a vector is passed, this will be different for each dynamical variable
dtval (float) – The starting integration timestep. This will be the exact timestep for fixed-step integrators, or stochastic integration.
kwargs (dict) – Arguments passed to scipy.integrate.solve_ivp.

Returns:

The integrated trajectory

Return type:

sol (ndarray)

dysts.utils.integration_utils.integrate_weiner(f, noise_amp, ic, tvals)¶

Given the RHS of a dynamical system, integrate the system assuming Brownian noise Requires the Python library sdeint

Parameters:

f (callable) – the right hand side of a system of ODE
noise_amp (float) – the amplitude of the Langevin forcing term

dysts.utils.integration_utils.resample_timepoints(model, ic, tpts, cutoff=10000, pts_per_period=100)¶

Given a differential equation, initial condition, and a set of integration points, determine a new set of timepoints that scales to the periodicity of the model

Parameters:

model (callable) – the right hand side of a set of ODEs
ic (list) – the initial conditions
tpts (array) – the timepoints over which to integrate
pts_per_period (int) – the number of timepoints to sample in each period

Returns:

The resampled timepoints

Return type:

new_timepoints (ndarray)

Analysis¶

Functions that act on DynSys or DynMap objects

dysts.analysis.calculate_lyapunov_exponent(traj1, traj2, dt=1.0)¶

Calculate the lyapunov exponent of two multidimensional trajectories using simple linear regression based on the log-transformed separation of the trajectories.

Parameters:

traj1 (np.ndarray) – trajectory 1 with shape (n_timesteps, n_dimensions)
traj2 (np.ndarray) – trajectory 2 with shape (n_timesteps, n_dimensions)
dt (float) – time step between timesteps

Returns:

lyapunov exponent

Return type:

float

dysts.analysis.compute_timestep(system: DynSys, total_length: int = 40000, transient_fraction: float = 0.2, num_iters: int = 20, pts_per_period: int = 1000, timescale: str = 'Fourier') → tuple[ndarray, ndarray]¶

Given a dynamical system object, find the integration timestep based on the largest signficant frequency

Parameters:

model (DynSys) – A dynamical systems object.
total_length (int) – The total trajectory length to use to determine timescales.
transient_fraction (float) – The fraction of a trajectory to discard as transient
num_iters (int) – The number of refinements to the timestep
pts_per_period (int) – The target integration timestep relative to the signal.
visualize (bool) – Whether to plot timestep versus time, in order to identify problems with the procedure

Returns: dt (float): The best integration timestep period (float, optional): The dominant timescale in the signal

dysts.analysis.corr_gpdim(traj1, traj2, register=False, standardize=False, **kwargs)¶

Given two multivariate time series, estimate their similarity using the cross Grassberger-Procaccia dimension

This quantity is defined as the cross-correlation between the two time series normalized by the product of the Grassberger-Procaccia dimension of each time series. np.sqrt(<ox oy> / ox oy)

Parameters:

traj1 (np.array) – T x D, where T is the number of timepoints, and D is the number of dimensions
traj2 (np.array) – T x D, where T is the number of timepoints, and D is the number of dimensions
register (bool) – Whether to register the two time series before computing the cross-correlation
standardize (bool) – Whether to standardize the time series before computing the cross-correlation

Returns:

The cross-correlation between the two time series

Return type:

float

dysts.analysis.estimate_powerlaw(data0)¶

Given a 1D array of continuous-valued data, estimate the power law exponent using the maximum likelihood estimator proposed by Clauset, Shalizi, Newman (2009).

Parameters:: data0 (np.ndarray) – An array of continuous-valued data
Returns:: The estimated power law exponent
Return type:: float

dysts.analysis.find_lyapunov_exponents(model, traj_length, pts_per_period=500, tol=1e-08, min_tpts=10, **kwargs)¶

Given a dynamical system, compute its spectrum of Lyapunov exponents. :param model: the right hand side of a differential equation, in format

func(X, t)

Parameters:

traj_length (int) – the length of each trajectory used to calulate Lyapunov exponents
pts_per_period (int) – the sampling density of the trajectory
kwargs – additional keyword arguments to pass to the model’s make_trajectory method

Returns:

A list of computed Lyapunov exponents

Return type:

final_lyap (ndarray)

References

Christiansen & Rugh (1997). Computing Lyapunov spectra with continuous: Gram-Schmidt orthonormalization

Example

>>> import dysts
>>> model = dysts.Lorenz()
>>> lyap = dysts.find_lyapunov_exponents(model, 1000, pts_per_period=1000)
>>> print(lyap)

dysts.analysis.gp_dim(data, y_data=None, rvals=None, nmax=100)¶

Estimate the Grassberger-Procaccia dimension for a numpy array using the empirical correlation integral.

Parameters:

data (np.array) – T x D, where T is the number of datapoints/timepoints, and D is the number of features/dimensions
y_data (np.array, Optional) – A second dataset of shape T2 x D, for computing cross-correlation.
rvals (np.array) – A list of radii
nmax (int) – The number of points at which to evaluate the correlation integral

Returns:

The discrete bins at which the correlation integral is: estimated

corr_sum (np.array): The estimates of the correlation integral at each bin

Return type:

rvals (np.array)

dysts.analysis.gpdistance(traj1, traj2, standardize=True, register=False, **kwargs)¶

Given two multivariate time series, estimate their similarity using the cross Grassberger-Procaccia distance

Parameters:

traj1 (np.array) – T x D, where T is the number of timepoints, and D is the number of dimensions
traj2 (np.array) – T x D, where T is the number of timepoints, and D is the number of dimensions
register (bool) – Whether to register the two time series before computing the cross-correlation
standardize (bool) – Whether to standardize the time series before computing the cross-correlation

dysts.analysis.kaplan_yorke_dimension(spectrum0)¶: Calculate the Kaplan-Yorke dimension, given a list of Lyapunov exponents

dysts.analysis.logarithmic_n(min_n, max_n, factor)¶

Return a list of values by successively multiplying a minimum value min_n by a factor > 1 until a maximum value max_n is reached. Non-integer results are rounded down. Based on a similar function in the nolds Python library.

Parameters:

min_n (float) – minimum value (must be < max_n)
max_n (float) – maximum value (must be > min_n)
factor (float) – factor used to increase min_n (must be > 1)

Returns:

min_n, min_n * factor, min_n * factor^2, … min_n * factor^i < max_n without duplicates

Return type:

list of integers

dysts.analysis.max_lyapunov_exponent(eq: DynSys, max_walltime: float, rtol: float = 0.001, atol: float = 1e-10, n_samples: int = 1000, traj_length: int = 5000, **kwargs)¶

Calculate the lyapunov spectrum of the system using a naive method based on the log-transformed separation of the trajectories over time.

Parameters:

eq (dysts.DynSys) – equation to calculate the lyapunov spectrum of
rtol (float) – relative tolerance for the separation of the trajectories
atol (float) – absolute tolerance for the separation of the trajectories
n_samples (int) – number of initial conditions to sample
traj_length (int) – length of the trajectories to sample. This should be long enough to ensure that most trajectories explore the attractor.
max_walltime (float) – maximum walltime in seconds to spend on the calculation of a given trajectory. If the calculation takes longer than this, the trajectory is discarded and a new one is sampled.
**kwargs – keyword arguments to pass to the sample / make_trajectory method of the dynamical equation

Returns:

largest lyapunov exponent

Return type:

float

Example

>>> import dysts
>>> eq = dysts.Lorenz()
>>> max_lyap = dysts.lyapunov_exponent_naive(eq)

dysts.analysis.rowwise_euclidean(x, y)¶: Computes the euclidean distance across rows

dysts.analysis.sample_initial_conditions(model, points_to_sample, traj_length=1000, pts_per_period=30)¶

Generate a random sample of initial conditions from a dynamical system

Parameters:

model (callable) – the right hand side of a differential equation, in format func(X, t)
points_to_sample (int) – the number of random initial conditions to sample
traj_length (int) – the total length of the reference trajectory from which points are drawn
pts_per_period (int) – the sampling density of the trajectory

Returns:

The points with shape (points_to_sample, d)

Return type:

sample_points (ndarray)

Systems¶

Utilities for the implemented systems

dysts.systems.compute_trajectory_statistics(n: int, subset: Sequence[str], datapath: PathLike | None = None, **kwargs) → dict[str, dict[str, ndarray[tuple[int, ...], dtype[float64]]]]¶

Compute mean and std for given trajectory list :param n: The number of timepoints to integrate :type n: int :param subset: A list of system names. Defaults to all continuous systems.

Can also pass in sys_class as a kwarg to specify other system classes.

Parameters:

datapath (str) – Path to save the computed statistics to a JSON file
kwargs – Additional keyword arguments passed to the integration routine

dysts.systems.get_attractor_list(sys_class: str = 'continuous', exclude: list[str] = []) → list[str]¶

Get names of implemented dynamical systems

Parameters:

sys_class – class of systems to get the name of - must be one of [‘continuous’, ‘continuous_no_delay’, ‘delay’, ‘discrete’]
exclude – list of systems to exclude, optional

Returns:

Sorted list of systems belonging to sys_class

dysts.systems.get_system_data(sys_class: str = 'continuous', exclude: list[str] = []) → dict[str, Any]¶

Get system data from the dedicated system class json files

Parameters:

sys_class – class of systems to get the name of - must be one of [‘continuous’, ‘continuous_no_delay’, ‘delay’, ‘discrete’]
exclude – list of systems to exclude, optional

Returns:

Data from json file filtered by sys_class as a dict

dysts.systems.make_trajectory_ensemble(n: int, use_tqdm: bool = True, use_multiprocessing: bool = False, subset: Sequence[str] | Sequence[BaseDyn] | None = None, event_fns: Sequence[Callable] | None = None, silent_errors: bool = False, **kwargs) → dict[str, ndarray[tuple[int, ...], dtype[float64]] | None]¶

Integrate multiple dynamical systems with identical settings

Parameters:

n (int) – The number of timepoints to integrate
use_tqdm (bool) – Whether to use a progress bar
use_multiprocessing (bool) – Not yet implemented.
subset (list) – A list of system names or BaseDyn (e.g. custom dynamical systems). Defaults to all continuous systems. Can also pass in sys_class as a kwarg to specify other system classes.
event_fns (list) –
A list of functions that can take the signatures:
1. (system: BaseDyn) -> (event_fn: Callable[[float, Array], float])
2. (event_fn: Callable[[float, Array], float])
3. () -> (event_fn: Callable[[float, Array], float])
If multiprocessing is used, the event functions must be of type 1 or 3 to avoid state persistence across processes.
silent_errors (bool) – Whether to fail silently if an error occurs
kwargs (dict) – Integration options passed to each system’s make_trajectory() method

Returns:

A dictionary containing trajectories for each system

Return type:

all_sols (dict)

Sampling¶

Example sampling functions for dysts

class dysts.sampling.BaseSampler(random_seed: int = 0)¶

Base class for any sampling-based transformations

set_rng(rng: Generator) → None¶: Required for multiprocessing

class dysts.sampling.GaussianInitialConditionSampler(random_seed: int = 0, scale: float = 0.0001, verbose: bool = False)¶: Sample gaussian perturbations for each initial condition in a given system list

class dysts.sampling.GaussianParamSampler(random_seed: int = 0, scale: float = 0.01, verbose: bool = False)¶

Sample gaussian perturbations for system parameters

Note

This is a MWE of a parameter transform
Other parameter transforms should follow this dataclass template

Parameters:

random_seed – for random sampling
scale – std (isotropic) of gaussian used for sampling

class dysts.sampling.OnAttractorInitCondSampler(random_seed: int = 0, reference_traj_length: int = 4096, reference_traj_transient: int = 500, trajectory_cache: dict[str, ~numpy.ndarray[tuple[int, ...], ~numpy.dtype[~numpy.float64]]] = <factory>, verbose: bool = False, events: list[~typing.Callable] | None = None)¶

Sample points from the attractor of a system

Subtleties:

This is slow, it requires integrating each system with its default parameters before sampling from the attractor.
The sampled initial conditions from this sampler are necessarily tied to the attractor defined by the default parameters.

Parameters:

reference_traj_length – Length of the reference trajectory to use for sampling ic on attractor.
reference_traj_transient – Transient length to ignore for the reference trajectory
trajectory_cache – Cache of reference trajectories for each system.
events – events to pass to solve_ivp

Metrics¶

Metrics for comparing two time series. These metrics are included to faciliate benchmarking of the algorithms in this package while reducing dependencies.

For more exhaustive sets of metrics, use the external tslearn, darts, or sktime libraries.

class dysts.metrics.GaussianMixture(means, covariances, weights=None)¶

A Gaussian Mixture Model class.

Parameters:

means (list) – A list of means for each component of the GMM.
covariances (list) – A list of covariance matrices for each component of the GMM.
weights (list) – A list of weights for each component of the GMM.

means¶

An array of means for each component of the GMM.

Type:: np.ndarray

covariances¶

An array of covariance matrices for each component of the GMM.

Type:: np.ndarray

weights¶

An array of weights for each component of the GMM.

Type:: np.ndarray

n_components¶

The number of components in the GMM.

Type:: int

gaussians¶

A list of multivariate_normal objects, one for each component

Type:: list

sample(n_samples=1)¶

Draw samples from the Gaussian Mixture Model.

Parameters:: n_samples (int) – The number of samples to draw.
Returns:: An array of shape (n_samples, ndim) containing the drawn samples.
Return type:: samples (np.ndarray)

dysts.metrics.are_broadcastable(shape1: tuple[int, ...], shape2: tuple[int, ...]) → bool¶

Check if two numpy arrays are broadcastable.

Parameters:

shape1 (tuple[int, ...]) – The shape of the first array
shape2 (tuple[int, ...]) – The shape of the second array

Returns:

True if the arrays are broadcastable, False otherwise

Return type:

bool

Examples

>>> are_broadcastable((1, 2, 3), (4, 5, 6))
False
>>> are_broadcastable((1, 2, 3), (3,))
True

dysts.metrics.average_hellinger_distance(ts_true: ndarray, ts_gen: ndarray, num_freq_bins: int = 100) → float¶

Compute the average Hellinger distance between power spectra of two multivariate time series.

Parameters:

ts_true (np.ndarray) – True time series, shape (n_samples, n_dimensions).
ts_gen (np.ndarray) – Generated time series, shape (n_samples, n_dimensions).
num_freq_bins (int) – Number of frequency bins to use in FFT for power spectrum.

Returns:

Average Hellinger distance across all dimensions.

Return type:

avg_dh

References

Mikhaeil et al. Advances in Neural Information Processing Systems, 35:: 11297–11312, December 2022.

dysts.metrics.calculate_season_error(y_past, m, time_dim=-1)¶

Calculate the mean absolute error between the forward and backward slices of the past data.

Parameters:

y_past (np.ndarray) – The past data
m (int) – The season length
time_dim (int) – The dimension of the time series

Returns:

The mean absolute error

Return type:

float

Examples

>>> calculate_season_error(y_past, m)

dysts.metrics.coefficient_of_variation(y_true, y_pred, eps=1e-10)¶

Coefficient of Variation of the root mean squared error relative to the mean of the true values

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The Coefficient of Variation

Return type:

float

dysts.metrics.compute_metrics(y_true: ndarray, y_pred: ndarray, verbose: bool = False, include: list[str] | None = None, batch_axis: int | None = None) → dict[str, float]¶

Compute multiple time series metrics

Parameters:

y_true (np.ndarray) – The true values of shape (…, T, …)
y_pred (np.ndarray) – The predicted values of shape (…, T, …)
verbose (bool) – Whether to print the computed metrics. Default is False.
include (optional, list) – The metrics to include. Default is None, which computes all metrics. Otherwise, specify a list of metrics to compute.
batch_axis (optional, int) – The axis to treat as the batch dimension. Default is None, which adds a singleton batchdimension at axis 0

Returns:

A dictionary containing the computed metrics

Return type:

dict

Raises:

ValueError – If the batch dimension is not the same for y_true and y_pred
ValueError – If the shapes of y_true and y_pred are not broadcastable

Examples

>>> compute_metrics(y_true, y_pred)
>>> compute_metrics(y_true, y_pred, batch_axis=1)
>>> compute_metrics(y_true, y_pred, include=["mse", "mae"])
>>> compute_metrics(y_true, y_pred, include=["mse", "mae"], batch_axis=1)

dysts.metrics.dtw(y_true, y_pred)¶

Compute the Dynamic Time Warping (DTW) distance between two time series.

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The DTW distance

Return type:

float

dysts.metrics.estimate_kl_divergence(true_orbit, generated_orbit, n_samples=1000, sigma_scale: float | None = 1.0) → float¶

Estimate KL divergence between observed and generated orbits using Gaussian Mixture Models (GMMs).

Parameters:

observed_orbit (np.ndarray) – Observed orbit points, with shape (T, N) where T is the number of time steps and N is the dimensionality.
generated_orbit (np.ndarray) – Generated orbit points, with shape (T, N) where T is the number of time steps and N is the dimensionality.
n_samples (int) – Number of Monte Carlo samples.
sigma_scale (float) – Variance parameter for the GMMs. If None, the variance is estimated locally for each point.

Returns:

Estimated KL divergence

Return type:

float

References

Hess, Florian, et al. “Generalized teacher forcing for learning chaotic dynamics.” Proceedings of the 40th International Conference on Machine Learning. 2023.

Development:: Rank-order (copula) transform each orbit coordinate, in order to reduce sensitivity to spacing among time series.

dysts.metrics.find_sigma(dists, tol=1e-06)¶

Given a list of distances to k nearest neighbors, find the sigma for each point

Parameters:

dists (np.ndarray) – A matrix of shape (k,)
tol (float) – The tolerance for the sigma

Returns:

The sigma np.ndarray: The transformed distances

Return type:

float

Examples

>>> sigma, dists_transformed = find_sigma(dists)
>>> sigma, dists_transformed = find_sigma(dists, tol=1e-4)

Raises:: ValueError – If the distance matrix is not a 2D array

dysts.metrics.hellinger_distance(p, q, axis=0)¶

Compute the Hellinger distance between two distributions.

Parameters:

p (np.ndarray) – The first distribution
q (np.ndarray) – The second distribution
axis (int) – The axis to sum over

Returns:

The Hellinger distance

Return type:

float

dysts.metrics.horizoned_metric(y_true, y_pred, metric, *args, horizon=None, **kwargs)¶

Compute a metric over a range of horizons

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values
metric (callable) – The metric function
*args – Additional arguments to pass to the metric function
horizon (int) – The maximum horizon to compute the metric over. If None, the horizon is set to the length of the time series
**kwargs – Additional keyword arguments to pass to the metric function

Returns:

The metric values at each horizon

Return type:

np.ndarray

dysts.metrics.kendall(y_true, y_pred)¶

Kendall-Tau Correlation. Returns dimensionwise mean for multivariate time series of shape (T, D)

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The Kendall-Tau Correlation

Return type:

float

dysts.metrics.mae(y_true, y_pred)¶

Mean Absolute Error

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The MAE

Return type:

float

dysts.metrics.mape(y_true, y_pred, eps=1e-10)¶

The Mean Absolute Percentage Error

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The MAPE

Return type:

float

dysts.metrics.marre(y_true, y_pred, eps=1e-10)¶

Mean Absolute Ranged Relative Error

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The MARRE

Return type:

float

dysts.metrics.mase(y, yhat, y_train=None, m=1, time_dim=-1, eps=1e-10)¶

The mean absolute scaled error.

Parameters:

y (ndarray) – The true values.
yhat (ndarray) – The predicted values.
y_train (ndarray) – The training values.
m (int) – The season length, which is the number of time steps that are skipped when computing the denominator. Default is 1.
time_dim (int) – The dimension of the time series. Default is -1.

Returns:

The MASE error

Return type:

mase_val (float)

dysts.metrics.mse(y_true, y_pred)¶

Mean Squared Error

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The MSE

Return type:

float

dysts.metrics.msis(y, yhat_lower, yhat_upper, y_obs, m, time_dim=-1, a=0.05, eps=1e-10)¶

The mean scaled interval score.

Parameters:

y (np.ndarray) – An array containing the true values.
yhat_lower – An array containing the a% quantile of the predicted distribution.
yhat_upper – An array containing the (1-a)% quantile of the predicted distribution.
y_obs – An array containing the training values.
m – The season length.
a – A scalar in [0, 1] specifying the quantile window to evaluate.

Returns:

The scalar MSIS.

dysts.metrics.mutual_information(y_true, y_pred)¶

Mutual Information. Returns dimensionwise mean for multivariate time series of shape (T, D). Computes the mutual information separately for each dimension and returns the mean.

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The Mutual Information

Return type:

float

dysts.metrics.nrmse(y_true, y_pred, eps=1e-10, scale=None)¶

Normalized Root Mean Squared Error

Parameters:

y_true (np.ndarray) – True values of shape (T, D)
y_pred (np.ndarray) – Predicted values of shape (T, D)
eps (float) – Small value to avoid division by zero
scale (np.ndarray) – Standard deviation of the true values of shape (D,). If None, the standard deviation is computed from the true values.

Returns:

NRMSE

Return type:

float

dysts.metrics.ope(y_true, y_pred)¶

Optimality Percentage Error

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The OPE

Return type:

float

dysts.metrics.pearson(y_true, y_pred)¶

Pearson Correlation. Returns dimensionwise mean for multivariate time series of shape (T, D)

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The Pearson Correlation

Return type:

float

dysts.metrics.r2_score(y_true, y_pred)¶

The R2 Score

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The R2 Score

Return type:

float

dysts.metrics.relu(x)¶

Rectified Linear Unit

Parameters:: x (np.ndarray) – The input array
Returns:: The output array
Return type:: np.ndarray

dysts.metrics.rmse(x, y)¶

Root Mean Squared Error

Parameters:

x (np.ndarray) – The true values
y (np.ndarray) – The predicted values

Returns:

The RMSE

Return type:

float

dysts.metrics.rmsle(y_true, y_pred, eps=1e-10)¶

Root Mean Squared Log Error. In case of negative values, the series is shifted to the positive domain.

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The RMSLE

Return type:

float

dysts.metrics.simplex_neighbors(X, metric='euclidean', k=20, tol=1e-06)¶

Compute the distance between points in a dataset using the simplex distance metric.

Parameters:

X (np.ndarray) – dataset of shape (n, d)
Y (np.ndarray) – dataset of shape (m, d)
metric (str) – distance metric to use
k (int) – number of nearest neighbors to use in the distance calculation
tol (float) – tolerance for the distance calculation

Returns:

weights matrix of shape (n,) idx (np.ndarray): index matrix of shape (n, k) sigmas (np.ndarray): sigmas matrix of shape (n,)

Return type:

wgts (np.ndarray)

Examples

>>> wgts, idx, sigmas = simplex_neighbors(X)
>>> wgts, idx, sigmas = simplex_neighbors(X, metric="cosine")
>>> wgts, idx, sigmas = simplex_neighbors(X, k=10)
>>> wgts, idx, sigmas = simplex_neighbors(X, tol=1e-4)

dysts.metrics.smape(x, y, eps=1e-10, scaled=True)¶: Symmetric mean absolute percentage error

dysts.metrics.spearman(y_true, y_pred)¶

Spearman Correlation. Returns dimensionwise mean for multivariate time series of shape (T, D)

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The Spearman Correlation

Return type:

float

dysts.metrics.wape(y_true, y_pred, eps=1e-10)¶

Weighted Absolute Percentage Error

Parameters:

y_true (np.ndarray) – The true values
y_pred (np.ndarray) – The predicted values

Returns:

The WAPE

Return type:

float

dysts API Reference¶

Base Classes¶

Utilities¶

Analysis¶

Systems¶

Sampling¶

Metrics¶

Indices and tables¶

dysts

Navigation

Related Topics