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
LyapunovsEstimator
to estimate the Backward and Forward Lyapunov Vectors (BLVs and FLVs) along a trajectoryCovariantLyapunovsEstimator
to estimate the Covariant Lyapunov Vectors (CLVs) along a trajectory
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:
multiprocessing.context.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.
- b
Vector of coefficients \(b_i\) of the Runge-Kutta method .
- Type
- c
Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .
- Type
- a
Vector of coefficients \(a_i\) of the Runge-Kutta method .
- Type
- 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
- b
Vector of coefficients \(b_i\) of the Runge-Kutta method .
- Type
- c
Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .
- Type
- a
Vector of coefficients \(a_i\) of the Runge-Kutta method .
- Type
- n_traj
The number of trajectories (initial conditions) computed at the last estimation performed by the estimator.
- Type
- 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
- 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
- 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 thenic
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
). Ifn_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
). Ifn_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
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
). Ifn_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
). Ifn_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
- 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
). Ifn_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
). Ifn_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
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)
wherex
is the state value andt
is the time.fjac (callable) – The Numba-jitted Jacobian matrix function \(\boldsymbol{J}\). Should have the signature
J(t, x)
wherex
is the state value andt
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.
- class qgs.toolbox.lyapunov.LyapProcess(processID, func, func_jac, b, c, a, ics_queue, lyap_queue)[source]
Bases:
multiprocessing.context.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.
- b
Vector of coefficients \(b_i\) of the Runge-Kutta method .
- Type
- c
Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .
- Type
- a
Vector of coefficients \(a_i\) of the Runge-Kutta method .
- Type
- 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
- b
Vector of coefficients \(b_i\) of the Runge-Kutta method .
- Type
- c
Matrix of coefficients \(c_{i,j}\) of the Runge-Kutta method .
- Type
- a
Vector of coefficients \(a_i\) of the Runge-Kutta method .
- Type
- n_traj
The number of trajectories (initial conditions) computed at the last estimation performed by the estimator.
- Type
- 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 thenic
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
). Ifn_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
). Ifn_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
- 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)
wherex
is the state value andt
is the time.fjac (callable) – The Numba-jitted Jacobian matrix function \(\boldsymbol{J}\). Should have the signature
J(t, x)
wherex
is the state value andt
is the time.
Warning
This function restarts the estimator!