Toolbox module

Module holding model’s tools (analysis, Lyapunov vectos, …).

Lyapunov module

Module with the classes of multi-thread the computation of the various Lyapunov vectors and exponents. Integrate using the Runge-Kutta method defined in the integrate module. See [LYAP-KP12] for more details on the Lyapunov vectors theoretical framework.

Module classes

References

LYAP-BGGS80

Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn. Lyapunov characteristic exponents for smooth dynamical systems and for hamiltonian systems; a method for computing all of them. part 1: theory. Meccanica, 15(1):9–20, 1980. URL: https://link.springer.com/article/10.1007/BF02128236.

LYAP-DPV21(1,2,3)

Jonathan Demaeyer, Stephen G Penny, and Stéphane Vannitsem. Identifying efficient ensemble perturbations for initializing subseasonal-to-seasonal prediction. arXiv preprint arXiv:2109.07979, 2021. URL: https://arxiv.org/abs/2109.07979.

LYAP-ER85(1,2,3)

J-P Eckmann and David Ruelle. Ergodic theory of chaos and strange attractors. The theory of chaotic attractors, pages 273–312, 1985. URL: https://link.springer.com/chapter/10.1007/978-0-387-21830-4_17.

LYAP-GPT+07(1,2,3)

Francesco Ginelli, Pietro Poggi, Alessio Turchi, Hugues Chaté, Roberto Livi, and Antonio Politi. Characterizing dynamics with covariant lyapunov vectors. Physical review letters, 99(13):130601, 2007. URL: https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.130601.

LYAP-KP12(1,2,3,4,5,6,7)

Pavel V Kuptsov and Ulrich Parlitz. Theory and computation of covariant lyapunov vectors. Journal of nonlinear science, 22(5):727–762, 2012. URL: https://link.springer.com/article/10.1007/s00332-012-9126-5.

class qgs.toolbox.lyapunov.ClvProcess(processID, func, func_jac, b, c, a, ics_queue, clv_queue, noise_pert)[source]

Bases: Process

CovariantLyapunovsEstimator’s workers class. Allows to multi-thread Lyapunov vectors estimation.

Parameters

  • processID (int) – Number identifying the worker.

  • func (callable) – Numba-jitted function to integrate assigned to the worker.

  • func_jac (callable) – Numba-jitted Jacobian matrix function to integrate assigned to the worker.

  • b (ndarray, optional) – Vector of coefficients \(b_i\) of the Runge-Kutta method .

  • c (ndarray, optional) – Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .

  • a (ndarray, optional) – Vector of coefficients \(a_i\) of the Runge-Kutta method .

  • ics_queue (multiprocessing.JoinableQueue) – Queue to which the worker ask for initial conditions and parameters input.

  • clv_queue (multiprocessing.Queue) – Queue to which the worker returns the estimation results.

processID

Number identifying the worker.

Type

int

func

Numba-jitted function to integrate assigned to the worker.

Type

callable

func_jac

Numba-jitted Jacobian matrix function to integrate assigned to the worker.

Type

callable

b

Vector of coefficients \(b_i\) of the Runge-Kutta method .

Type

ndarray

c

Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .

Type

ndarray

a

Vector of coefficients \(a_i\) of the Runge-Kutta method .

Type

ndarray

run()[source]

Main worker computing routine. Perform the estimation with the fetched initial conditions and parameters.

class qgs.toolbox.lyapunov.CovariantLyapunovsEstimator(num_threads=None, b=None, c=None, a=None, number_of_dimensions=None, noise_pert=0.0, method=0)[source]

Bases: object

Class to compute the Covariant Lyapunov vectors (CLVs) and exponents along a trajectory of a dynamical system

\[\dot{\boldsymbol{x}} = \boldsymbol{f}(t, \boldsymbol{x})\]

with a set of LyapProcess and a specified Runge-Kutta method. The tangent linear model must also be provided. I.e. one must provide the linearized ODEs

\[\dot{\boldsymbol{\delta x}} = \boldsymbol{\mathrm{J}}(t, \boldsymbol{x}) \cdot \boldsymbol{\delta x}\]

where \(\boldsymbol{\mathrm{J}} = \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{x}}\) is the Jacobian matrix of \(\boldsymbol{f}\).

Parameters

  • num_threads (None or int, optional) – Number of LyapProcess workers (threads) to use. If None, use the number of machine’s cores available. Default to None.

  • b (None or ndarray, optional) – Vector of coefficients \(b_i\) of the Runge-Kutta method . If None, use the classic RK4 method coefficients. Default to None.

  • c (None or ndarray, optional) – Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method . If None, use the classic RK4 method coefficients. Default to None.

  • a (None or ndarray, optional) – Vector of coefficients \(a_i\) of the Runge-Kutta method . If None, use the classic RK4 method coefficients. Default to None.

  • number_of_dimensions (None or int, optional) – Allow to hardcode the dynamical system dimension. If None, evaluate the dimension from the callable func. Default to None.

  • method (int, optional) –

    Allow to select the method used to compute the CLVs. Presently can be 0 or 1:

    • 0: Uses the method of Ginelli et al. [LYAP-GPT+07]. Suitable for a trajectory not too long (depends on the memory available).

    • 1: Uses the method of the intersection of the subspace spanned by the BLVs and FLVs described in [LYAP-ER85] and [LYAP-KP12] (see also [LYAP-DPV21], Appendix A). Suitable for longer trajectories (uses less memory).

    Default to 0, i.e. Ginelli et al. algorithm.

  • noise_pert (float, optional) – Noise perturbation amplitude parameter of the diagonal of the R matrix in the QR decomposition during the Ginelli step. Mainly done to avoid ill-conditioned matrices near tangencies (see [LYAP-KP12]). Default to 0 (no perturbation). Only apply if using the Ginelli et al. algorithm, i.e. if method=0.

num_threads

Number of LyapProcess workers (threads) to use.

Type

int

b

Vector of coefficients \(b_i\) of the Runge-Kutta method .

Type

ndarray

c

Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .

Type

ndarray

a

Vector of coefficients \(a_i\) of the Runge-Kutta method .

Type

ndarray

n_dim

Dynamical system dimension.

Type

int

n_vec

The number of Lyapunov vectors to compute.

Type

int

n_traj

The number of trajectories (initial conditions) computed at the last estimation performed by the estimator.

Type

int

n_records

The number of saved states of the last estimation performed by the estimator.

Type

int

ic

Store the estimator initial conditions.

Type

ndarray

func

Last function \(\boldsymbol{f}\) used by the estimator.

Type

callable

func_jac

Last Jacobian matrix function \(\boldsymbol{J}\) used by the estimator.

Type

callable

method

Select the method used to compute the CLVs:

  • 0: Uses the method of Ginelli et al. [LYAP-GPT+07]. Suitable for a trajectory not too long (depends on the memory available).

  • 1: Uses the method of the intersection of the subspaces spanned by the BLVs and FLVs described in [LYAP-ER85] and [LYAP-KP12] (see also [LYAP-DPV21], Appendix A). Suitable for longer trajectories (uses less memory).

Type

int

noise_pert

Noise perturbation parameter of the diagonal of the matrix resulting from the backpropagation during the Ginelli step. Mainly done to avoid ill-conditioned matrices near tangencies (see [LYAP-KP12]). Only apply if using the Ginelli et al. algorithm, i.e. if method=0.

Type

float

compute_clvs(t0, ta, tb, tc, dt, mdt, ic=None, write_steps=1, n_vec=None, method=None, backward_vectors=False, forward_vectors=False)[source]

Estimate the Covariant Lyapunov Vectors (CLVs) along a given trajectory, always integrating the said trajectory forward in time from ic at t0 to time tc. Return the CLVs between ta and tb. The result of the estimation can be obtained afterward by calling get_clvs().

Parameters

  • t0 (float) – Initial time of the time integration. Corresponds to the initial condition’s ic time.

  • ta (float) – Define the time span between t0 and ta of the first part of the algorithm, which obtain the convergence to the Backward Lyapunov vectors (initialization of the Benettin algorithm).

  • tb (float) – Define the time span between ta and tb where the Covariant Lyapunov Vectors are computed.

  • tc (float) – Final time of the time integration algorithm. Define the time span between tb and tc where, depending on the value of method, the convergence to the Forward Lyapunov Vectors or to the Covariant Lyapunov Vectors (thanks to the Ginelli steps) is obtained.

  • dt (float) – Timestep of the integration.

  • mdt (float) – Micro-timestep to integrate the tangent linear equation between the nonlinear system dt timesteps. Should be smaller or equal to dt.

  • ic (None or ndarray(float), optional) –

    Initial conditions of the system. Can be a 1D or a 2D array:

    • 1D: Provide a single initial condition. Should be of shape (n_dim,) where n_dim = \(\mathrm{dim}(\boldsymbol{x})\).

    • 2D: Provide an ensemble of initial condition. Should be of shape (n_traj, n_dim) where n_dim = \(\mathrm{dim}(\boldsymbol{x})\), and where n_traj is the number of initial conditions.

    If None, use the initial conditions stored in ic. If then ic is None, use a zero initial condition. Default to None.

  • write_steps (int, optional) – Save the state of the integration in memory every write_steps steps. The other intermediary steps are lost. It determines the size of the returned objects. Default is 1. Set to 0 to return only the final state.

  • n_vec (int, optional) – The number of Lyapunov vectors to compute. Should be smaller or equal to n_dim.

  • method (int, optional) –

    Allow to select the method used to compute the CLVs. Presently can be 0 or 1:

    • 0: Uses the method of Ginelli et al. [LYAP-GPT+07]. Suitable for a trajectory not too long (depends on the memory available).

    • 1: Uses the method of the intersection of the subspace spanned by the BLVs and FLVs described in [LYAP-ER85] and [LYAP-KP12] (see also [LYAP-DPV21], Appendix A). Suitable for longer trajectories (uses less memory).

    Use the Ginelli et al. algorithm if not provided.

  • backward_vectors (bool, optional) – Store also the computed Backward Lyapunov vectors between ta and tb. Only applies if method=1. Does not store the BLVs if not provided.

  • forward_vectors (bool, optional) – Store also the computed Forward Lyapunov vectors between ta and tb. Only applies if method=1. Does not store the FLVs if not provided.

get_blvs()[source]

Returns the BLVs obtained during the previous CLVs estimation.

Returns

time, traj, exponents, vectors – The result of the estimation:

  • time: Time at which the state of the system was saved. Array of shape (n_records,).

  • traj: Saved dynamical system states. 3D array of shape (n_traj, n_dim, n_records). If n_traj = 1, a 2D array of shape (n_dim, n_records) is returned instead.

  • exponents: Saved estimates of the local Lyapunov exponents along the trajectory. 3D array of shape (n_traj, n_vec, n_records). If n_traj = 1, a 2D array of shape (n_vec, n_records) is returned instead.

  • vectors: Saved estimates of the local Lyapunov vectors along the trajectory. Depending on the input initial conditions, it is maximum a 4D array of shape (n_traj, n_dim, n_vec, n_records). If one of the dimension is 1, it is squeezed.

Return type

ndarray

Warning

The BLVs are only available if method is set to 1.

get_clvs()[source]

Returns the result of the previous CLVs estimation.

Returns

time, traj, exponents, vectors – The result of the estimation:

  • time: Time at which the state of the system was saved. Array of shape (n_records,).

  • traj: Saved dynamical system states. 3D array of shape (n_traj, n_dim, n_records). If n_traj = 1, a 2D array of shape (n_dim, n_records) is returned instead.

  • exponents: Saved estimates of the local Lyapunov exponents along the trajectory. 3D array of shape (n_traj, n_vec, n_records). If n_traj = 1, a 2D array of shape (n_vec, n_records) is returned instead.

  • vectors: Saved estimates of the local Lyapunov vectors along the trajectory. Depending on the input initial conditions, it is maximum a 4D array of shape (n_traj, n_dim, n_vec, n_records). If one of the dimension is 1, it is squeezed.

Return type

ndarray

get_flvs()[source]

Returns the FLVs obtained during the previous CLVs estimation.

Returns

time, traj, exponents, vectors – The result of the estimation:

  • time: Time at which the state of the system was saved. Array of shape (n_records,).

  • traj: Saved dynamical system states. 3D array of shape (n_traj, n_dim, n_records). If n_traj = 1, a 2D array of shape (n_dim, n_records) is returned instead.

  • exponents: Saved estimates of the local Lyapunov exponents along the trajectory. 3D array of shape (n_traj, n_vec, n_records). If n_traj = 1, a 2D array of shape (n_vec, n_records) is returned instead.

  • vectors: Saved estimates of the local Lyapunov vectors along the trajectory. Depending on the input initial conditions, it is maximum a 4D array of shape (n_traj, n_dim, n_vec, n_records). If one of the dimension is 1, it is squeezed.

Return type

ndarray

Warning

The FLVs are only available if method is set to 1.

set_bca(b=None, c=None, a=None, ic_init=True)[source]

Set the coefficients of the Runge-Kutta method and restart the estimator.

Parameters

  • b (None or ndarray, optional) – Vector of coefficients \(b_i\) of the Runge-Kutta method . If None, does not reinitialize these coefficients.

  • c (None or ndarray, optional) – Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method . If None, does not reinitialize these coefficients.

  • a (None or ndarray, optional) – Vector of coefficients \(a_i\) of the Runge-Kutta method . If None, does not reinitialize these coefficients.

  • ic_init (bool, optional) – Re-initialize or not the initial conditions of the estimator. Default to True.

set_func(f, fjac)[source]

Set the Numba-jitted function \(\boldsymbol{f}\) and Jacobian matrix function \(\boldsymbol{\mathrm{J}}\) to integrate.

Parameters

  • f (callable) – The Numba-jitted function \(\boldsymbol{f}\). Should have the signature f(t, x) where x is the state value and t is the time.

  • fjac (callable) – The Numba-jitted Jacobian matrix function \(\boldsymbol{J}\). Should have the signature J(t, x) where x is the state value and t is the time.

Warning

This function restarts the estimator!

set_noise_pert(noise_pert)[source]

Set the noise perturbation noise_pert parameter.

Parameters

noise_pert (float, optional) – Noise perturbation amplitude parameter of the diagonal of the R matrix in the QR decomposition during the Ginelli step. Mainly done to avoid ill-conditioned matrices near tangencies (see [LYAP-KP12]). Only apply if using the Ginelli et al. algorithm, i.e. if method is 0.

start()[source]

Start or restart the workers (threads) of the estimator.

Warning

If the estimator was not previously terminated, it will be terminated first in the case of a restart.

terminate()[source]

Stop the workers (threads) and release the resources of the estimator.

class qgs.toolbox.lyapunov.LyapProcess(processID, func, func_jac, b, c, a, ics_queue, lyap_queue)[source]

Bases: Process

LyapunovsEstimator’s workers class. Allows to multi-thread Lyapunov vectors estimation.

Parameters

  • processID (int) – Number identifying the worker.

  • func (callable) – Numba-jitted function to integrate assigned to the worker.

  • func_jac (callable) – Numba-jitted Jacobian matrix function to integrate assigned to the worker.

  • b (ndarray, optional) – Vector of coefficients \(b_i\) of the Runge-Kutta method .

  • c (ndarray, optional) – Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .

  • a (ndarray, optional) – Vector of coefficients \(a_i\) of the Runge-Kutta method .

  • ics_queue (multiprocessing.JoinableQueue) – Queue to which the worker ask for initial conditions and parameters input.

  • lyap_queue (multiprocessing.Queue) – Queue to which the worker returns the estimation results.

processID

Number identifying the worker.

Type

int

func

Numba-jitted function to integrate assigned to the worker.

Type

callable

func_jac

Numba-jitted Jacobian matrix function to integrate assigned to the worker.

Type

callable

b

Vector of coefficients \(b_i\) of the Runge-Kutta method .

Type

ndarray

c

Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .

Type

ndarray

a

Vector of coefficients \(a_i\) of the Runge-Kutta method .

Type

ndarray

run()[source]

Main worker computing routine. Perform the estimation with the fetched initial conditions and parameters.

class qgs.toolbox.lyapunov.LyapunovsEstimator(num_threads=None, b=None, c=None, a=None, number_of_dimensions=None)[source]

Bases: object

Class to compute the Forward and Backward Lyapunov vectors and exponents along a trajectory of a dynamical system

\[\dot{\boldsymbol{x}} = \boldsymbol{f}(t, \boldsymbol{x})\]

with a set of LyapProcess and a specified Runge-Kutta method. The tangent linear model must also be provided. I.e. one must provide the linearized ODEs

\[\dot{\boldsymbol{\delta x}} = \boldsymbol{\mathrm{J}}(t, \boldsymbol{x}) \cdot \boldsymbol{\delta x}\]

where \(\boldsymbol{\mathrm{J}} = \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{x}}\) is the Jacobian matrix of \(\boldsymbol{f}\). The method used to compute the Lyapunov vectors is the one introduced by Benettin et al. [LYAP-BGGS80].

Parameters

  • num_threads (None or int, optional) – Number of LyapProcess workers (threads) to use. If None, use the number of machine’s cores available. Default to None.

  • b (None or ndarray, optional) – Vector of coefficients \(b_i\) of the Runge-Kutta method . If None, use the classic RK4 method coefficients. Default to None.

  • c (None or ndarray, optional) – Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method . If None, use the classic RK4 method coefficients. Default to None.

  • a (None or ndarray, optional) – Vector of coefficients \(a_i\) of the Runge-Kutta method . If None, use the classic RK4 method coefficients. Default to None.

  • number_of_dimensions (None or int, optional) – Allow to hardcode the dynamical system dimension. If None, evaluate the dimension from the callable func. Default to None.

num_threads

Number of LyapProcess workers (threads) to use.

Type

int

b

Vector of coefficients \(b_i\) of the Runge-Kutta method .

Type

ndarray

c

Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .

Type

ndarray

a

Vector of coefficients \(a_i\) of the Runge-Kutta method .

Type

ndarray

n_dim

Dynamical system dimension.

Type

int

n_vec

The number of Lyapunov vectors to compute.

Type

int

n_traj

The number of trajectories (initial conditions) computed at the last estimation performed by the estimator.

Type

int

n_records

The number of saved states of the last estimation performed by the estimator.

Type

int

ic

Store the estimator initial conditions.

Type

ndarray

func

Last function \(\boldsymbol{f}\) used by the estimator.

Type

callable

func_jac

Last Jacobian matrix function \(\boldsymbol{J}\) used by the estimator.

Type

callable

compute_lyapunovs(t0, tw, t, dt, mdt, ic=None, write_steps=1, n_vec=None, forward=False, adjoint=False, inverse=False)[source]

Estimate the Lyapunov vectors using the Benettin algorithm along a given trajectory, always integrating the said trajectory forward in time from ic at t0 to time t. The result of the estimation can be obtained afterward by calling get_lyapunovs(). If forward is True, it yields the Forward Lyapunov Vectors (FLVs) between t0 and tw, otherwise, returns the Backward Lyapunov Vectors (BLVs) between tw and t.

Parameters

  • t0 (float) – Initial time of the time integration. Corresponds to the initial condition’s ic time.

  • tw (float) – Time at which the algorithm start to store the Lyapunov vectors. Define thus also the transient before the which the Lyapunov vectors are considered as having not yet converged. Must be between t0 and t.

  • t (float) – Final time of the time integration. Corresponds to the final condition.

  • dt (float) – Timestep of the integration.

  • mdt (float) – Micro-timestep to integrate the tangent linear equation between the nonlinear system dt timesteps. Should be smaller or equal to dt.

  • ic (None or ndarray(float), optional) –

    Initial conditions of the system. Can be a 1D or a 2D array:

    • 1D: Provide a single initial condition. Should be of shape (n_dim,) where n_dim = \(\mathrm{dim}(\boldsymbol{x})\).

    • 2D: Provide an ensemble of initial condition. Should be of shape (n_traj, n_dim) where n_dim = \(\mathrm{dim}(\boldsymbol{x})\), and where n_traj is the number of initial conditions.

    If None, use the initial conditions stored in ic. If then ic is None, use a zero initial condition. Default to None.

  • forward (bool, optional) – If True, yield the Forward Lyapunov Vectors (FLVs) between t0 and tw. If False, yield the Backward Lyapunov Vectors (BLVs) between tw and t. Default to False, i.e. Backward Lyapunov Vectors estimation.

  • adjoint (bool, optional) – If true, integrate the tangent \(\dot{\boldsymbol{\delta x}} = \boldsymbol{\mathrm{J}}(t, \boldsymbol{x}) \cdot \boldsymbol{\delta x}\) , else, integrate the adjoint linear model \(\dot{\boldsymbol{\delta x}} = \boldsymbol{\mathrm{J}}^T(t, \boldsymbol{x}) \cdot \boldsymbol{\delta x}\). Integrate the tangent model by default.

  • inverse (bool, optional) – Whether or not to invert the Jacobian matrix \(\boldsymbol{\mathrm{J}}(t, \boldsymbol{x}) \rightarrow \boldsymbol{\mathrm{J}}^{-1}(t, \boldsymbol{x})\). False by default.

  • write_steps (int, optional) – Save the state of the integration in memory every write_steps steps. The other intermediary steps are lost. It determines the size of the returned objects. Default is 1. Set to 0 to return only the final state.

  • n_vec (int, optional) – The number of Lyapunov vectors to compute. Should be smaller or equal to n_dim.

get_lyapunovs()[source]

Returns the result of the previous Lyapunov vectors estimation.

Returns

time, traj, exponents, vectors – The result of the estimation:

  • time: Time at which the state of the system was saved. Array of shape (n_records,).

  • traj: Saved dynamical system states. 3D array of shape (n_traj, n_dim, n_records). If n_traj = 1, a 2D array of shape (n_dim, n_records) is returned instead.

  • exponents: Saved estimates of the local Lyapunov exponents along the trajectory. 3D array of shape (n_traj, n_vec, n_records). If n_traj = 1, a 2D array of shape (n_vec, n_records) is returned instead.

  • vectors: Saved estimates of the local Lyapunov vectors along the trajectory. Depending on the input initial conditions, it is maximum a 4D array of shape (n_traj, n_dim, n_vec, n_records). If one of the dimension is 1, it is squeezed.

Return type

ndarray

set_bca(b=None, c=None, a=None, ic_init=True)[source]

Set the coefficients of the Runge-Kutta method and restart the estimator.

Parameters

  • b (None or ndarray, optional) – Vector of coefficients \(b_i\) of the Runge-Kutta method . If None, does not reinitialize these coefficients.

  • c (None or ndarray, optional) – Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method . If None, does not reinitialize these coefficients.

  • a (None or ndarray, optional) – Vector of coefficients \(a_i\) of the Runge-Kutta method . If None, does not reinitialize these coefficients.

  • ic_init (bool, optional) – Re-initialize or not the initial conditions of the estimator. Default to True.

set_func(f, fjac)[source]

Set the Numba-jitted function \(\boldsymbol{f}\) and Jacobian matrix function \(\boldsymbol{\mathrm{J}}\) to integrate.

Parameters

  • f (callable) – The Numba-jitted function \(\boldsymbol{f}\). Should have the signature f(t, x) where x is the state value and t is the time.

  • fjac (callable) – The Numba-jitted Jacobian matrix function \(\boldsymbol{J}\). Should have the signature J(t, x) where x is the state value and t is the time.

Warning

This function restarts the estimator!

start()[source]

Start or restart the workers (threads) of the estimator.

Warning

If the estimator was not previously terminated, it will be terminated first in the case of a restart.

terminate()[source]

Stop the workers (threads) and release the resources of the estimator.