Tools for the exploration of chaos and nonlinear dynamics

PlaneCrossing(plane, dir) → z

Create a struct that can be called as a function z(u) that returns the signed distance of state u from the hyperplane plane (positive means in front of the hyperplane). See poincaresos for what plane can be (tuple or vector).

_ac_period(v, t; ϵ = 0.2, L = length(v)÷10)

Use the autocorrelation function (AC). The value where the AC first comes back close to 1 is the period of the signal. The keyword L = length(v)÷10 denotes the length of the AC (thus, given the default setting, this method will fail if there less than 10 periods in the signal). The keyword ϵ = 0.2 means that 1-ϵ counts as "1" for the AC.

_ls_period(v, t;
            minimum_period::Real = 2 * (t[end] - t[1]) / length(t),
            maximum_period::Real = (t[end] - t[1]) / 1.5,
            kwargs... # see the documentation of LombScargle.jl for these
    )

Uses the Lomb-Scargle algorithm to compute a periodogram. The advantage of the Lomb-Scargle method is that it does not require an equally sampled dataset, and it performs well on undersampled datasets. Constraints have been set on the period, since Lomb-Scargle tends to have false peaks at very low frequencies. That being said, it's a very flexible method. It is extremely customizable, and the keyword arguments that can be passed to it are given in the documentation.

_mt_period(v, t;
                nw = 4, ntapers = DSP.ceil(2nw)-1,
                window = DSP.dpss(length(s), nw, ntapers), kwargs...)

The multitaper method reduces estimation bias by obtaining multiple independent estimates from the same sample. Data tapers are then windowed and the power spectra are obtained. nw is the time-bandwidth product, and ntapers is the number of tapers. If window is not specified, the signal is tapered with ntapers discrete prolate spheroidal sequences with time-bandwidth product nw. Each sequence is equally weighted; adaptive multitaper is not (yet) supported. If window is specified, each column is applied as a taper. The sum of periodograms is normalized by the total sum of squares of window.

_periodogram_period(v, t; kwargs...)

Use the fast Fourier transform to compute a periodogram (power-spectrum) of the given data. Data must be evenly sampled.

_zc_period(v, t; line = mean(v))

Find the zero crossings of the data, and use the average difference between zero crossings as the period. This is a naïve implementation, with no interpolation; however, it's useful as a sanity check.

The keyword line controls where the "zero point" is.

The implementation of the function was inspired by this gist, and has been modified for performance and to support arbitrary time grids.

basin_fractions(basins::Array) → fs::Dict

Calculate the fraction of the basins of attraction encoded in basins. The elements of basins are integers, enumerating the attractor that the entry of basins converges to. Return a dictionary that maps attractor IDs to their relative fractions.

In^[Menck2013] the authors use these fractions to quantify the stability of a basin of attraction, and specifically how it changes when a parameter is changed.

basins_2D(xg, yg, integ; kwargs...) → basins, attractors

Compute an estimate of the basins of attraction of a "two dimensional system" of the plane onto itself according to the method of Nusse & Yorke^[Yorke1997]. The dynamical system can be:

• An actual 2D DiscreteDynamicalSystem or ContinuousDynamicalSystem.
• 2D poincaré map of a 3D ContinuousDynamicalSystem.
• A 2D stroboscopic map, i.e. a periodically forced 2D ContinuousDynamicalSystem.

For a higher-dimensional dynamical systems, use basins_general.

integ is an istance of an integrator, not a DynamicalSystem. This includes the output of poincaremap. See documentation online for examples for all cases! xg, yg are 1-dimensional ranges that define the grid of the initial conditions to test. The output basins is a matrix on the grid (xg, yg), see below for details. The output attractors is a dictionary whose keys correspond to the attractor number and the values contains the points of the attractors found on the map. Notice that for some attractors this list may be incomplete.

See also match_attractors!, basin_fractions, tipping_probabilities.

Keyword Arguments

• T : Period of the stroboscopic map, in case integ is an integrator of a 2D continuous dynamical system with periodic time forcing.
• mc_att = 10: A parameter that sets the maximum checks of consecutives hits of an attractor before deciding the basin of the initial condition.
• mc_bas = 10 : Maximum check of consecutive visits of the same basin of attraction. This number can be increased for higher accuracy.
• mc_unmb = 60 : Maximum check of unnumbered cell before considering we have an attractor. This number can be increased for higher accuracy.

Description

This method identifies the attractors and their basins of attraction on the grid without prior knowledge about the system. At the end of a successfull computation the function returns a matrix coding the basins of attraction and a dictionary with all attractors found.

basins has the following organization:

• The basins are coded in sequential order from 1 to the number of attractors Na .
• If the trajectory diverges or converges to an attractor outside the defined grid it is numbered -1

attractors has the following organization:

• The keys of the dictionary correspond to the number of the attractor.
• The value associated to a key is a Dataset with the guessed location of the attractor on the state space.

The method starts by picking the first available initial condition on the plane not yet numbered. The dynamical system is then iterated until one of the following conditions happens:

1. The trajectory hits a known attractor already numbered mc_att consecutive times: the initial condition is numbered with the corresponding number.
2. The trajectory diverges or hits an attractor outside the defined grid: the initial condition is set to -1
3. The trajectory hits a known basin mc_bas times in a row: the initial condition belongs to that basin and is numbered accordingly.
4. The trajectory hits mc_unmb times in a row an unnumbered cell: it is considered an attractor and is labelled with a new number.

Regarding performace, this method is at worst as fast as tracking the attractors. In most cases there is a signicative improvement in speed.

basins_general(xg, yg, ds::DynamicalSystem; kwargs...) -> basin, attractors

Compute an estimate of the basins of attraction of a higher-dimensional dynamical system ds on a projection of the system dynamics on a two-dimensional plane.

Like basins_2D, xg, yg are ranges defining the grid of initial conditions on the plane. Refer to basins_2D for more details regarding the algorithm. Notice that to use the efficient algorithm of basins_2D we have to project the dynamics on a 2D plane. There are edge cases where the system may have two attractors that are close on the plane but are far apart in another dimension. They could be collapsed or confused into the same attractor. This is a drawback of this method.

This function can be used to make attractor basins of higher dimension via the complete_state keyword. E.g. to make 3D basins you can make many 2D basins slices and concatenate them. For example:

zg = 0:0.01:1 # the range defining the z part of the grid
bs, as = [], []
for z ∈ zg
    b, a = basins_general(xg, yg, ds; complete_state = [z, 0.0])
    push!(bs, b); push!(as, a)
end
# use `match_attractors!` to match potential basins, then do:
basins_3D = cat(bs...; dims = 3)

Keyword Arguments

• dt = 1: Approximate time step of the integrator. It is recommended to use values ≥ 1.
• idxs = 1:2: This vector selects the two variables of the system that will define the "plane" the dynamics will be projected into.
• complete_state = zeros(D-2): This argument allows setting the remaining variables of the dynamical system state on each planar initial condition x, y. It can be either a vector of length D-2, or a function f(x, y) that returns a vector of length D-2.
• mc_att, mc_bas, mc_unmb: As in basins_2D.
• diffeq...: Keyword arguments propagated to integrator.

boxed_correlationsum(data, εs, r0 = maximum(εs); q = 2 , M = size(data, 2), w = 0)
boxed_correlationsum(data; q = 2 , M = size(data, 2), w = 0)

Estimate the box assisted q-order correlation sum^[Kantz2003] out of a dataset data for radii εs by splitting the data into boxes of size r0 beforehand.

Description

C_q(ε) is calculated for every ε ∈ εs and each of the boxes to then be summed up afterwards. If M is unequal to the dimension of the data, only the first M dimensions are considered for the box distribution (this is called the prism-assisted version). The method of splitting the data into boxes was implemented according to Theiler^{[Theiler1987]}. If only a dataset is given the radii εs and boxsize r0 are chosen by calculating estimate_r0_buenoorovio.

w is the Theiler window.

The function is explicitly optimized for q = 2 and becomes quite slow for q ≠ 2.

See correlationsum for the definition of C_q and also data_boxing.

boxed_correlationsum_2(boxes, contents, data, εs; w = 0)

For a vector of boxes and the indices of their contents inside of data, calculate the classic correlationsum of a radius or multiple radii εs. w is the Theiler window, for explanation see boxed_correlationsum.

boxed_correlationsum_q(boxes, contents, data, εs, q; w = 0)

For a vector of boxes and the indices of their contents inside of data, calculate the q-order correlationsum of a radius or radii εs. w is the Theiler window, for explanation see boxed_correlationsum.

boxregion(as, bs) -> sampler, restraining

Define a box in $\mathbb{R}^d$ with edges the as and bs and then return two functions: sampler, which generates a random initial condition in that box and restraining that returns true if a given state is in the box.

broomhead_king(s::AbstractVector, d::Int) -> U, S, Vtr

Return the Broomhead-King coordinates of a timeseries s by performing svd on high-dimensional embedding if s with dimension d with minimum delay.

Description

Broomhead and King coordinates is an approach proposed in ^{[Broomhead1987]} that applies the Karhunen–Loève theorem to delay coordinates embedding with smallest possible delay.

The function performs singular value decomposition on the d-dimensional matrix $X$ of $s$,

\[X = \frac{1}{\sqrt{N}}\left( \begin{array}{cccc} x_1 & x_2 & \ldots & x_d \\ x_2 & x_3 & \ldots & x_{d+1}\\ \vdots & \vdots & \vdots & \vdots \\ x_{N-d+1} & x_{N-d+2} &\ldots & x_N \end{array} \right) = U\cdot S \cdot V^{tr}.\]

where $x := s - \bar{s}$. The columns of $U$ can then be used as a new coordinate system, and by considering the values of the singular values $S$ you can decide how many columns of $U$ are "important". See the documentation page for example application.

capacitydim(args...) = generalizeddim(args...; α = 0)

correlationsum(X, ε::Real; w = 0, norm = Euclidean(), q = 2) → C_q(ε)

Calculate the q-order correlation sum of X (Dataset or timeseries) for a given radius ε and norm.

The function boxed_correlationsum is faster and should be preferred over this one.

Description

The correlation sum is done using the formula:

\[C_2(\epsilon) = \frac{2}{(N-w)(N-w-1)}\sum_{i=1}^{N}\sum_{j=1+w+i}^{N} B(||X_i - X_j|| < \epsilon)\]

for q=2 and

\[C_q(\epsilon) = \left[\frac{1}{\alpha} \sum_{i=w+1}^{N-w}\left[\sum_{j:|i-j| > w} B(||X_i - X_j|| < \epsilon)\right]^{q-1}\right]^{1/(q-1)}\]

where

\[\alpha = (N-2w)(N-2w-1)^{(q-1)}\]

for q≠2, where $N$ is its length and $B$ gives 1 if the argument is true. w is the Theiler window. If ε is a vector its values have to be ordered. See the article of Grassberger for the general definition ^{[Grassberger2007]} and the book "Nonlinear Time Series Analysis" ^[Kantz2003], Ch. 6, for a discussion around w and choosing best values and Ch. 11.3 for the explicit definition of the q-order correlationsum.

correlationsum(X, εs::AbstractVector; w, norm, q) → C_q(ε)

If εs is a vector, C_q is calculated for each ε ∈ εs. If also q=2, some strong optimizations are done, but this requires the allocation a matrix of size N×N. If this is larger than your available memory please use instead:

[correlationsum(..., ε) for ε in εs]

See grassberger for more. See also takens_best_estimate.

correlationsum_fixedmass(data, max_j; metric = Euclidean(), M = length(data)) → rs, ys

A fixed mass algorithm for the calculation of the fractal dimension $\Delta$ according to [^Grassberger 1988] with max_j the maximum number of neighbours that should be considered for the calculation, M defines the number of points considered for the calculation, default is the whole data set.

Implements

\[\Delta \overline{\log r^{(j)}} \sim Ψ(j) - \log N\]

where $\Psi(j) = \frac{\text{d} \log Γ(j)}{\text{d} j}$, rs = $\overline{\log r^{(j)}}$ and ys = $\Psi(j) - \log N$.

$\Delta$ can be computed by using linear_region(rs, ys).

[^Grassberger 1988]: Peter Grassberger (1988) Finite sample Corrections to Entropy and Dimension Estimates, Physics Letters A 128(6-7)

data_boxing(data, r0, M = size(data, 2))

Distribute the data points into boxes of size r0. Return box positions and the contents of each box as two separate vectors. Implemented according to the paper by Theiler^{[Theiler1987]} improving the algorithm by Grassberger and Procaccia^{[Grassberger1983]}. If M is smaller than the dimension of the data, only the first M dimensions are considered for the distribution into boxes.

See also: boxed_correlationsum.

Return the true distance of u from u0 according to metric defined by ε.

draw_basin!(xg, yg, integ, iter_f!::Function, reinit_f!::Function)

Compute an estimate of the basin of attraction on a two-dimensional plane. This is a low level function, for higher level functions see: basins_2D, basins_general

Arguments

• xg, yg : 1-dim range vector that defines the grid of the initial conditions to test.
• integ : integrator handle of the dynamical system.
• iter_f! : function that iterates the map or the system, see step! from DifferentialEquations.jl and

examples for a Poincaré map of a continuous system.

• reinit_f! : function that sets the initial condition to test on a two dimensional projection of the phase space.

dyca(data, eig_thresold) -> eigenvalues, proj_mat, projected_data

Compute the Dynamical Component analysis (DyCA) of the given data^[Uhl2018] used for dimensionality reduction.

Return the eigenvalues, projection matrix, and reduced-dimension data (which are just data*proj_mat).

Description

Dynamical Component Analysis (DyCA) is a method to detect projection vectors to reduce the dimensionality of multi-variate, high-dimensional deterministic datasets. Unlike methods like PCA or ICA that make a stochasticity assumption, DyCA relies on a determinacy assumption on the time-series and is based on the solution of a generalized eigenvalue problem. After choosing an appropriate eigenvalue threshold and solving the eigenvalue problem, the obtained eigenvectors are used to project the high-dimensional dataset onto a lower dimension. The obtained eigenvalues measure the quality of the assumption of linear determinism for the investigated data. Furthermore, the number of the generalized eigenvalues with a value of approximately 1.0 are a measure of the number of linear equations contained in the dataset. This property is useful in detecting regions with highly deterministic parts in the time-series and also as a preprocessing step for reservoir computing of high-dimensional spatio-temporal data.

The generalised eigenvalue problem we solve is:

\[C_1 C_0^{-1} C_1^{\top} \bar{u} = \lambda C_2 \bar{u} \]

where $C_0$ is the correlation matrix of the data with itself, $C_1$ the correlation matrix of the data with its derivative, and $C_2$ the correlation matrix of the derivative of the data with itself. The eigenvectors $\bar{u}$ with eigenvalues approximately 1 and their $C_1^{-1} C_2 u$ counterpart, form the space where the data is projected onto.

estimate_boxsizes(A::Dataset; kwargs...) → εs

Return k exponentially spaced values: εs = base .^ range(lower + w, upper + z; length = k), that are a good estimate for sizes ε that are used in calculating a Fractal Dimension. It is strongly recommended to regularize input dataset A before using this function.

Let d₋ be the minimum pair-wise distance in A and d₊ the maximum distance along each of the dimensions of A. Then lower = log(base, d₋) and upper = log(base, d₊). Because by default w=1, z=-1, the returned sizes are an order of mangitude larger than the minimum distance, and an order of magnitude smaller than the maximum distance.

Keywords

• w = 1, z = -1, k = 20 : as explained above.
• base = MathConstants.e : the base used in the log function.

estimate_period(v::Vector, method, t=0:length(v)-1; kwargs...)

Estimate the period of the signal v, with accompanying time vector t, using the given method.

If t is an AbstractArray, then it is iterated through to ensure that it's evenly sampled (if necessary for the algorithm). To avoid this, you can pass any AbstractRange, like a UnitRange or a LinRange, which are defined to be evenly sampled.

Methods requiring evenly sampled data

These methods are faster, but some are error-prone.

• :periodogram or :pg: Use the fast Fourier transform to compute a periodogram (power-spectrum) of the given data. Data must be evenly sampled.

• :multitaper or mt: The multitaper method reduces estimation bias by using multiple independent estimates from the same sample. Data tapers are then windowed and the power spectra are obtained. Available keywords follow: nw is the time-bandwidth product, and ntapers is the number of tapers. If window is not specified, the signal is tapered with ntapers discrete prolate spheroidal sequences with time-bandwidth product nw. Each sequence is equally weighted; adaptive multitaper is not (yet) supported. If window is specified, each column is applied as a taper. The sum of periodograms is normalized by the total sum of squares of window.

• :autocorrelation or :ac: Use the autocorrelation function (AC). The value where the AC first comes back close to 1 is the period of the signal. The keyword L = length(v)÷10 denotes the length of the AC (thus, given the default setting, this method will fail if there less than 10 periods in the signal). The keyword ϵ = 0.2 (\epsilon) means that 1-ϵ counts as "1" for the AC.

Methods not requiring evenly sampled data

These methods tend to be slow, but versatile and low-error.

• :lombscargle or :ls: Use the Lomb-Scargle algorithm to compute a periodogram. The advantage of the Lomb-Scargle method is that it does not require an equally sampled dataset and performs well on undersampled datasets. Constraints have been set on the period, since Lomb-Scargle tends to have false peaks at very low frequencies. That being said, it's a very flexible method. It is extremely customizable, and the keyword arguments that can be passed to it are given in the documentation.

• :zerocrossing or :zc: Find the zero crossings of the data, and use the average difference between zero crossings as the period. This is a naïve implementation, with only linear interpolation; however, it's useful as a sanity check. The keyword line controls where the "crossing point" is. It deffaults to mean(v).

For more information on the periodogram methods, see the documentation of DSP.jl and LombScargle.jl.

estimate_r0_buenoorovio(X, M = size(X, 2))

Estimates a reasonable size for boxing the time series X proposed by Bueno-Orovio and Pérez-García^[Bueno2007] before calculating the correlation dimension as presented by Theiler^{[Theiler1983]}. If instead of boxes, prisms are chosen everything stays the same but M is the dimension of the prism. To do so the dimension ν is estimated by running the algorithm by Grassberger and Procaccia^{[Grassberger1983]} with √N points where N is the number of total data points. An effective size ℓ of the attractor is calculated by boxing a small subset of size N/10 into boxes of sidelength r_ℓ and counting the number of filled boxes η_ℓ.

\[\ell = r_\ell \eta_\ell ^{1/\nu}\]

The optimal number of filled boxes η_opt is calculated by minimising the number of calculations.

\[\eta_\textrm{opt} = N^{2/3}\cdot \frac{3^\nu - 1}{3^M - 1}^{1/2}.\]

M is the dimension of the data or the number of edges on the prism that don't span the whole dataset.

Then the optimal boxsize $r_0$ computes as

\[r_0 = \ell / \eta_\textrm{opt}^{1/\nu}.\]

estimate_r0_theiler(data)

Estimates a reasonable size for boxing the data before calculating the correlation dimension proposed by Theiler^{[Theiler1987]}. To do so the dimension is estimated by running the algorithm by Grassberger and Procaccia^{[Grassberger1983]} with √N points where N is the number of total data points. Then the optimal boxsize $r_0$ computes as

\[r_0 = R (2/N)^{1/\nu}\]

where $R$ is the size of the chaotic attractor and $\nu$ is the estimated dimension.

exit_entry_times(ds, u₀, εs, T; diffeq...) → exits, entries

Collect exit and entry times for a ball/box centered at u₀ with radii εs (see below), in the state space of the given discrete dynamical system (function not yet available for continuous systems). Return the exit and (re-)entry return times to the set(s), where each of these is a vector containing all collected times for the respective ε-radius set, for ε ∈ εs.

Use transit_return(exits, entries) to transform the output into transit and return times, and see also mean_return_times for both continuous and discrete systems.

Description

Transit time statistics are important for the transport properties of dynamical systems^[Meiss1997] and can even be connected with the fractal dimension of chaotic sets^[Boev2014].

The current algorithm collects exit and re-entry times to given sets in the state space, which are centered at u₀ (algorithm always starts at u₀ and the initial state of ds is irrelevant). εs is always a Vector.

The sets around u₀ are nested hyper-spheres of radius ε ∈ εs, if each entry of εs is a real number. The sets can also be hyper-rectangles (boxes), if each entry of εs is a vector itself. Then, the i-th box is defined by the space covered by u0 .± εs[i] (thus the actual box size is 2εs[i]!).

The reason to input multiple εs at once is purely for performance.

For discrete systems, exit time is recorded immediatelly after exitting of the set, and re-entry is recorded immediatelly on re-entry. This means that if an orbit needs 1 step to leave the set and then it re-enters immediatelly on the next step, the return time is 1. For continuous systems high-order interpolation is done to accurately record the time of exactly crossing the ε-ball/box.

expansionentropy(ds::DynamicalSystem, sampler, restraining; kwargs...)

Calculate the expansion entropy^[Hunt2015] of ds, in the restraining region $S$ defined by restraining, by estimating the slope of the biggest linear region of the curve $\log E_{t0+T, t0}(f, S)$ versus $T$ (using linear_region). This is an approximation of the expansion entropy $H_0$, according to^[Hunt2015].

sampler is a 0-argument function that generates a random initial condition (a sample) of ds. restraining is a 1-argument function restraining(u) that given the state u it returns true if the state is inside the restraining region $S$.

Use boxregion for an easy way to define sampler and restraining on a multidimension box.

Keyword Arguments

• N = 1000 : Number of samples taken at each batch (same as $N$ of [1]).
• steps = 40 : The maximal steps for which the system will be run.
• Ttr = 0 : Transient time to evolve each initial condition before starting to comute $E$. This is t0 of [1] and of the following notation.
• batches = 100 : Number of batches to run the calculation, see below.
• diffeq... : Other keywords are propagated to the solvers of DifferentialEquations.jl.

Description

N samples are initialized and propagated forwards in time (along with their tangent space). At every time $t$ in [t0+dt, t0+2dt, ... t0+steps*dt] we calculate $H$:

\[H[t] = \log E_{t0+T, t0}(f, S),\]

with

\[E_{t0+T, t0}(f, S) = \frac 1 N \sum_{i'} G(Df_{t0+t, t0}(x_i))\]

(using same notation as ^[Hunt2015]). In principle $E$ is the average largest possible growth ratio within the restraining region (sampled by the initial conditions). The summation is only over $x_i$ that stay inside the region $S$ defined by the boolean function restraining. This process is done by the expansionentropy_sample function.

Then, this is repeated for batches amount of times, as recommended in^[Hunt2015]. From all these batches, the mean and std of $H$ is computed at every time point. This is done by the expansionentropy_batch function. When plotted versus $t$, these create the curves and error bars of e.g. Figs 2, 3 of [1].

This function expansionentropy simply returns the slope of the biggest linear region of the curve $H$ versus $t$, which approximates the expansion entropy $H_0$. It is therefore recommended to use expansionentropy_batch directly and evaluate the result yourself, as this step is known to be inaccurate for non-chaotic systems (where $H$ fluctuates strongly around 0).

expansionentropy_batch(ds, sampler, restraining; kwargs...)

Run expansionentropy_samplebatch times, and return times, mean(H), std(H) for all resulting H, see expansionentropy.

Accepts the same arguments as expansionentropy.

expansionentropy_sample(ds, sampler, restraining; kwargs...)

Return times, H for one sample of ds (see expansionentropy). Accepts the same argumets as expansionentropy, besides batches.

find_neighborboxes_2(index, boxes, contents) → indices

For an index into boxes all neighbouring boxes beginning from the current one are searched. If the found box is indeed a neighbour, the contents of that box are added to indices.

find_neighborboxes_q(index, boxes, contents, q) → indices

For an index into boxes all neighbouring boxes are searched. If the found box is indeed a neighbour, the contents of that box are added to indices.

gali(ds::DynamicalSystem, tmax, k::Int | Q0; kwargs...) -> GALI_k, t

Compute $\text{GALI}_k$^[Skokos2007] for a given k up to time tmax. Return $\text{GALI}_k(t)$ and time vector $t$.

The third argument, which sets the order of gali, can be an integer k, or a matrix with its columns being the deviation vectors (then k = size(Q0)[2]). In the first case random orthonormal vectors are chosen.

Keyword Arguments

• threshold = 1e-12 : If GALI_k falls below the threshold iteration is terminated.
• dt = 1 : Time-step between deviation vector normalizations. For continuous systems this is approximate.
• u0 : Initial state for the system. Defaults to get_state(ds).
• diffeq... : Keyword arguments propagated into init of DifferentialEquations.jl. See trajectory for examples. Only valid for continuous systems.

Description

The Generalized Alignment Index, $\text{GALI}_k$, is an efficient (and very fast) indicator of chaotic or regular behavior type in $D$-dimensional Hamiltonian systems ($D$ is number of variables). The asymptotic behavior of $\text{GALI}_k(t)$ depends critically on the type of orbit resulting from the initial condition. If it is a chaotic orbit, then

\[\text{GALI}_k(t) \sim \exp\left[\sum_{j=1}^k (\lambda_1 - \lambda_j)t \right]\]

with $\lambda_j$ being the j-th Lyapunov exponent (see lyapunov, lyapunovspectrum). If on the other hand the orbit is regular, corresponding to movement in $d$-dimensional torus with $1 \le d \le D/2$ then it holds

\[\text{GALI}_k(t) \sim \begin{cases} \text{const.}, & \text{if} \;\; 2 \le k \le d \; \; \text{and} \; \;d > 1 \\ t^{-(k - d)}, & \text{if} \;\; d < k \le D - d \\ t^{-(2k - D)}, & \text{if} \;\; D - d < k \le D \end{cases}\]

Traditionally, if $\text{GALI}_k(t)$ does not become less than the threshold until tmax the given orbit is said to be chaotic, otherwise it is regular.

Our implementation is not based on the original paper, but rather in the method described in^{[Skokos2016b]}, which uses the product of the singular values of $A$, a matrix that has as columns the deviation vectors.

Performance Notes

This function uses a tangent_integrator. For loops over initial conditions and/or parameter values one should use the low level method that accepts an integrator, and reinit! it to new initial conditions. See the "advanced documentation" for info on the integrator object. The low level method is

ChaosTools.gali(tinteg, tmax, dt, threshold)

(section 5.3.1 and ref. [85] therein), Lecture Notes in Physics 915, Springer (2016)

generalized_dim(dataset [, sizes]; q = 1, base = MathConstants.e) -> D_α

Return the α order generalized dimension of the dataset, by calculating the genentropy for each ε ∈ sizes.

The case of α=0 is often called "capacity" or "box-counting" dimension.

Description

The returned dimension is approximated by the (inverse) power law exponent of the scaling of the genentropy versus the box size ε, where ε ∈ sizes.

Calling this function performs a lot of automated steps:

1. A vector of box sizes is decided by calling sizes = estimate_boxsizes(dataset), if sizes is not given.
2. For each element of sizes the appropriate entropy is calculated, through h = genentropy.(Ref(dataset), sizes; α, base). Let x = -log.(sizes).
3. The curve h(x) is decomposed into linear regions, using linear_regions(x, h).
4. The biggest linear region is chosen, and a fit for the slope of that region is performed using the function linear_region, which does a simple linear regression fit using linreg. This slope is the return value of generalized_dim.

By doing these steps one by one yourself, you can adjust the keyword arguments given to each of these function calls, refining the accuracy of the result.

The following aliases are provided:

• α = 0 : boxcounting_dim, capacity_dim
• α = 1 : information_dim

grassberger(data, εs = estimate_boxsizes(data); kwargs...) → D_C

Use the method of Grassberger and Proccacia^{[Grassberger1983]}, and the correction by Theiler^{[Theiler1986]}, to estimate the correlation dimension D_C of the given data.

This function does something extrely simple:

cm = correlationsum(data, εs; kwargs...)
return linear_region(log.(sizes), log(cm))[2]

i.e. it calculates correlationsum for various radii and then tries to find a linear region in the plot of the log of the correlation sum versus log(ε). See generalized_dim for a more thorough explanation.

See also takens_best_estimate.

informationdim(args...) = generalizeddim(args...; α = 1)

inner_correlationsum_2(indices_X, indices_Y, data, εs; norm = Euclidean(), w = 0)

Calculates the classic correlation sum for values X inside a box, considering Y consisting of all values in that box and the ones in neighbouring boxes for all distances ε ∈ εs calculated by norm. To obtain the position of the values in the original time series data, they are passed as indices_X and indices_Y.

w is the Theiler window. Each index to the original array is checked for the distance of the compared index. If this absolute value is not higher than w its element is not used in the calculation of the correlationsum.

See also: correlationsum

inner_correlationsum_q(indices_X, indices_Y, data, εs, q::Real; norm = Euclidean(), w = 0)

Calculates the q-order correlation sum for values X inside a box, considering Y consisting of all values in that box and the ones in neighbouring boxes for all distances ε ∈ εs calculated by norm. To obtain the position of the values in the original time series data, they are passed as indices_X and indices_Y.

w is the Theiler window. The first and last w points of this data set are not used by themselves to calculate the correlationsum.

See also: correlationsum

Return true if state is outside ε-ball

kaplanyorke_dim(λs::AbstractVector)

Calculate the Kaplan-Yorke dimension, a.k.a. Lyapunov dimension^[Kaplan1970].

Description

The Kaplan-Yorke dimension is simply the point where cumsum(λs) becomes zero (interpolated):

\[ D_{KY} = k + \frac{\sum_{i=1}^k \lambda_i}{|\lambda_{k+1}|},\quad k = \max_j \left[ \sum_{i=1}^j \lambda_i > 0 \right].\]

If the sum of the exponents never becomes negative the function will return the length of the input vector.

Useful in combination with lyapunovspectrum.

kernelprob(X, ε; norm = Euclidean()) → p::Probabilities

Associate each point in X (Dataset or timesries) with a probability p using the "kernel estimation" (also called "nearest neighbor kernel estimation" and other names):

\[p_j = \frac{1}{N}\sum_{i=1}^N B(||X_i - X_j|| < \epsilon)\]

where $N$ is its length and $B$ gives 1 if the argument is true.

See also genentropy and correlationsum. kernelprob is equivalent with probabilities(X, NaiveKernel(ε, TreeDistance(norm))).

lambdamatrix(λ, inds::Vector{Int}, sings) -> Λk

Return the matrix $\mathbf{\Lambda}_k$ used to create a new dynamical system with some unstable fixed points turned to stable in the function periodicorbits.

Arguments

1. λ<:Real : the multiplier of the $C_k$ matrix, with 0<λ<1.
2. inds::Vector{Int} : The ith entry of this vector gives the row of the nonzero element of the ith column of $C_k$.
3. sings::Vector{<:Real} : The element of the ith column of $C_k$ is +1 if signs[i] > 0 and -1 otherwise (sings can also be Bool vector).

Calling lambdamatrix(λ, D::Int) creates a random $\mathbf{\Lambda}_k$ by randomly generating an inds and a signs from all possible combinations. The collections of all these combinations can be obtained from the function lambdaperms.

Description

Each element of indsmust be unique such that the resulting matrix is orthogonal and represents the group of special reflections and permutations.

Deciding the appropriate values for λ, inds, sings is not trivial. However, in ref.^[Pingel2000] there is a lot of information that can help with that decision. Also, by appropriately choosing various values for λ, one can sort periodic orbits from e.g. least unstable to most unstable, see^{[Diakonos1998]} for details.

lambdaperms(D) -> indperms, singperms

Return two collections that each contain all possible combinations of indices (total of $D!$) and signs (total of $2^D$) for dimension D (see lambdamatrix).

linear_region(x, y; kwargs...) -> ((ind1, ind2), slope)

Call linear_regions and identify and return the largest linear region and its slope. The region starts and stops at x[ind1:ind2].

The keywords dxi, tol are propagated as-is to linear_regions. The keyword ignore_saturation = true ignores saturation that (typically) happens at the final points of the curve y(x), where the curve flattens out.

linear_regions(x, y; dxi::Int = 1, tol = 0.25) -> (lrs, tangents)

Identify regions where the curve y(x) is linear, by scanning the x-axis every dxi indices sequentially (e.g. at x[1] to x[5], x[5] to x[10], x[10] to x[15] and so on if dxi=5).

If the slope (calculated via linear regression) of a region of width dxi is approximatelly equal to that of the previous region, within tolerance tol, then these two regions belong to the same linear region.

Return the indices of x that correspond to linear regions, lrs, and the correcttangents at each region (obtained via a second linear regression at each accumulated region).

linreg(x, y) -> a, b

Perform a linear regression to find the best coefficients so that the curve: z = a + b*x has the least squared error with y.

local_growth_rates(ds, points::Dataset; S=100, Δt=5, kwargs...) → λlocal

Compute the exponential local growth rate(s) of perturbations of the dynamical system ds for initial conditions given in points. For each initial condition u ∈ points, S total perturbations are created and evolved for time Δt. The exponential local growth rate is defined simply by log(g/g0)/Δt with g0 the initial pertrubation size and g the size after Δt. Thus, λlocal is a matrix of size (length(points), S).

This function is a modification of lyapunov. It uses the full nonlinear dynamics to evolve the perturbations, but does not do any re-scaling, thus allowing probing state and time dependence of perturbation growth. The actual growth is given by exp(λlocal * Δt).

The output of this function is sometimes referred as "Nonlinear Local Lyapunov Exponent".

Keyword Arguments

• perturbation: If given, it should be a function perturbation(ds, u, j) that outputs a pertrubation vector (preferrably SVector) given the system, current initial condition u and the counter j ∈ 1:S. If not given, a random perturbation is generated with norm given by the keyword e = 1e-6.
• diffeq...: Keywords propagated to the solvers of DifferentialEquations.jl.

lyapunov(ds::DynamicalSystem, Τ; kwargs...) -> λ

Calculate the maximum Lyapunov exponent λ using a method due to Benettin ^{[Benettin1976]}, which simply evolves two neighboring trajectories (one called "given" and one called "test") while constantly rescaling the test one. T denotes the total time of evolution (should be Int for discrete systems).

See also lyapunovspectrum, local_growth_rates.

Keyword Arguments

• u0 = get_state(ds) : Initial condition.
• Ttr = 0 : Extra "transient" time to evolve the trajectories before starting to measure the expontent. Should be Int for discrete systems.
• d0 = 1e-9 : Initial & rescaling distance between the two neighboring trajectories.
• upper_threshold = 1e-6 : Upper distance threshold for rescaling.
• lower_threshold = 1e-12 : Lower distance threshold for rescaling (in order to be able to detect negative exponents).
• dt = 1 : Time of evolution between each check of distance exceeding the thresholds. For continuous systems this is approximate.
• inittest = (u1, d0) -> u1 .+ d0/sqrt(D) : A function that given (u1, d0) initializes the test state with distance d0 from the given state u1 (D is the dimension of the system). This function can be used when you want to avoid the test state appearing in a region of the phase-space where it would have e.g. different energy or escape to infinity.
• diffeq... : Keyword arguments propagated into init of DifferentialEquations.jl. See trajectory for examples. Only valid for continuous systems.

Description

Two neighboring trajectories with initial distance d0 are evolved in time. At time $t_i$ their distance $d(t_i)$ either exceeds the upper_threshold, or is lower than lower_threshold, which initializes a rescaling of the test trajectory back to having distance d0 from the given one, while the rescaling keeps the difference vector along the maximal expansion/contraction direction: $u_2 \to u_1+(u_2−u_1)/(d(t_i)/d_0)$.

The maximum Lyapunov exponent is the average of the time-local Lyapunov exponents

\[\lambda = \frac{1}{t_{n} - t_0}\sum_{i=1}^{n} \ln\left( a_i \right),\quad a_i = \frac{d(t_{i})}{d_0}.\]

Performance Notes

This function uses a parallel_integrator. For loops over initial conditions and/or parameter values one should use the low level method that accepts an integrator, and reinit! it to new initial conditions. See the "advanced documentation" for info on the integrator object. The low level method is

lyapunov(pinteg, T, Ttr, dt, d0, ut, lt)

lyapunovspectrum(ds::DynamicalSystem, N [, k::Int | Q0]; kwargs...) -> λs

Calculate the spectrum of Lyapunov exponents ^{[Lyapunov1992]} of ds by applying a QR-decomposition on the parallelepiped matrix N times. Return the spectrum sorted from maximum to minimum.

The third argument k is optional, and dictates how many lyapunov exponents to calculate (defaults to dimension(ds)). Instead of passing an integer k you can pass a pre-initialized matrix Q0 whose columns are initial deviation vectors (then k = size(Q0)[2]).

See also lyapunov, local_growth_rates.

Keyword Arguments

• u0 = get_state(ds) : State to start from.
• Ttr = 0 : Extra "transient" time to evolve the system before application of the algorithm. Should be Int for discrete systems. Both the system and the deviation vectors are evolved for this time.
• dt = 1 : Time of individual evolutions between successive orthonormalization steps. For continuous systems this is approximate.
• diffeq... : Keyword arguments propagated into init of DifferentialEquations.jl. See trajectory for examples. Only valid for continuous systems.

Description

The method we employ is "H2" of ^[Geist1990], originally stated in ^{[Benettin1980]}. The deviation vectors defining a D-dimensional parallepiped in tangent space are evolved using the tangent dynamics of the system. A QR-decomposition at each step yields the local growth rate for each dimension of the parallepiped. The growth rates are then averaged over N successive steps, yielding the lyapunov exponent spectrum (at each step the parallepiped is re-normalized).

Performance Notes

This function uses a tangent_integrator. For loops over initial conditions and/or parameter values one should use the low level method that accepts an integrator, and reinit! it to new initial conditions. See the "advanced documentation" for info on the integrator object. The low level method is

lyapunovspectrum(tinteg, N, dt::Real, Ttr::Real)

If you want to obtain the convergence timeseries of the Lyapunov spectrum, use the method

ChaosTools.lyapunovspectrum_convergence(tinteg, N, dt, Ttr)

(not exported).

match_attractors!(b₋, a₋, b₊, a₊, [, method = :distance])

Attempt to match the attractors in basins/attractors b₊, a₊ with those at b₋, a₋. b is an array whose values encode the attractor ID, while a is a dictionary mapping IDs to Datasets containing the attractors (e.g. output of basins_general). Typically the +,- mean after and before some change of parameter for a system.

In basins_general different attractors get assigned different IDs, however which attractor gets which ID is somewhat arbitrary, and computing the basins of the same system for slightly different parameters could label the "same" attractors (at the different parameters) with different IDs. match_attractors! tries to "match" them by modifying the attractor IDs.

The modification of IDs is always done on the b, a that have less attractors.

method decides the matching process:

• method = :overlap matches attractors whose basins before and after have the most overlap (in pixels).
• method = :distance matches attractors whose state space distance the smallest.

matrix_gradient(matrix)

Compute the gradient of a matrix along axis=1.

Description

Compute the gradient using second order accurate central differences in the interior points and first order accurate one-sides differences at the boundaries. We find the standard second order approximation by using:

\[\hat{f}_i^{(1)} = \frac{f(x_{i+1}-f(x_{i-1})}{2h} + O(h^2)\]

The returned gradient matrix hence has the same shape as the input array. Here we compute the gradient along axis=1 (row-wise), so to compute gradient along axis=2 (column-wise), the tranpose of the input matrix must be given.

maximalexpansion(M)

Calculates the maximal expansion rate of M, i.e. the product of all singular values of M that are greater than 1. In the notation of [1], it is the function $G$.

mean_return_times(ds::DynamicalSystem, u₀, εs, T; kwargs...) → τ, c

Return the mean return times τ, as well as the amount of returns c, for subsets of the state space of ds defined by u₀, εs. The ds is evolved for a maximum of T time. This function behaves similarly to exit_entry_times and thus see that one for the meaning of u₀ and εs.

This function supports both discrete and continuous systems, however the optimizations done in discrete systems (where all nested ε-sets are checked at the same time), are not done here yet, which leads to disproportionally lower performance since each ε-related set is checked individually from start.

Continuous systems allow for the following keywords:

• i=10 How many points to interpolate the trajectory in-between steps to find candidate crossing regions.
• dmin If the trajectory is at least dmin distance away from u0, the algorithm that checks for crossings of the ε-set is not initiated. By default obtains the a value 4 times as large as the radius of the maximum ε-set.
• diffeq... All extra keywords are propagated to solvers of DifferentialEquations.jl, see trajectory. It is strongly recommended to use high accuracy keywords here, e.g. alg = Vern9(), reltol = 1e-12, abstol = 1e-12, maxiters = typemax(Int).

For continuous systems T, i, dmin can be vectors with same size as εs, to help increase accuracy of small ε.

minimum_pairwise_distance(A::Dataset, metric = Euclidean())

Return min_d, min_pair: the minimum pairwise distance of all points in the dataset, and the corresponding point pair.

modified mean square displacement

molteno_boxing(data::Dataset; k0 = 10) → (probs, εs)

Distribute the data into boxes whose size is halved in each step. Stop if the average number of points per filled box falls below the threshold k0.

Return probs, a vector of Propabilities for different box sizes and the corresponding box sizes εs.

Description

Project the data onto the whole interval of numbers that is covered by UInt64. This projected data is distributed into boxes whose size decreases by factor 2 in each step. For each box that contains more than one point 2^D new boxes are created where D is the dimension of the data.

The process of dividing the data into new boxes stops when the number of points over the number of filled boxes falls below k0. The box sizes εs are calculated and returned together with the probs.

See^{[Molteno1993]} for more.

molteno_dim(data::Dataset; k0 = 10, q = 1.0, base = ℯ)

Calculate the generalized dimension using the algorithm for box division defined by Molteno^{[Molteno1993]}.

Description

Divide the data into boxes with each new box having half the side length of the former box using molteno_boxing. Break if the number of points over the number of filled boxes falls below k0. Then the generalized dimension can be calculated by using genentropy to calculate the sum over the logarithm also considering possible approximations and fitting this to the logarithm of one over the boxsize using linear_region.

This algorithm is faster than the traditional approach of using probabilities, but it is only suited for low dimensional data since it divides each box into 2^D new boxes if D is the dimension. For large D this leads to low numbers of box divisions before the threshold is passed and the divison stops and as a result to a low number of data points to fit the dimension to and thereby a poor estimate.

molteno_subboxes(box, data, iteration)

Divide a box containing indices into data to 2^D smaller boxes and sort the points contained in the former box into the new boxes. Implemented according to Molteno^[Molteno]. Sorts the elements of the former box into the smaller boxes using cheap bit shifting and & operations on the value of data at each box element. iteration determines which bit of the array should be shifted to the last position.

normalize_deviations!(tinteg)

Normalize (in-place) the deviation vectors of the tangent integrator.

numericallyapunov(R::Dataset, ks;  refstates, w, distance, ntype)

Return E = [E(k) for k ∈ ks], where E(k) is the average logarithmic distance between states of a neighborhood that are evolved in time for k steps (k must be integer). Typically R is the result of delay coordinates of a single timeseries.

Keyword Arguments

• refstates = 1:(length(R) - ks[end]) : Vector of indices that notes which states of the reconstruction should be used as "reference states", which means that the algorithm is applied for all state indices contained in refstates.
• w::Int = 1 : The Theiler window.
• ntype = NeighborNumber(1) : The neighborhood type. Either NeighborNumber or WithinRange. See Neighborhoods for more info.
• distance::Metric = Cityblock() : The distance function used in the logarithmic distance of nearby states. The allowed distances are Cityblock() and Euclidean(). See below for more info. The metric for finding neighbors is always the Euclidean one.

Description

If the dataset exhibits exponential divergence of nearby states, then it should hold

\[E(k) \approx \lambda\cdot k \cdot \Delta t + E(0)\]

for a well defined region in the k axis, where $\lambda$ is the approximated maximum Lyapunov exponent. $\Delta t$ is the time between samples in the original timeseries. You can use linear_region with arguments (ks .* Δt, E) to identify the slope (= $\lambda$) immediatelly, assuming you have choosen sufficiently good ks such that the linear scaling region is bigger than the saturated region.

The algorithm used in this function is due to Parlitz^[Skokos2016], which itself expands upon Kantz ^[Kantz1994]. In sort, for each reference state a neighborhood is evaluated. Then, for each point in this neighborhood, the logarithmic distance between reference state and neighborhood state(s) is calculated as the "time" index k increases. The average of the above over all neighborhood states over all reference states is the returned result.

If the Metric is Euclidean() then use the Euclidean distance of the full D-dimensional points (distance $d_E$ in ref.^[Skokos2016]). If however the Metric is Cityblock(), calculate the absolute distance of only the first elements of the m+k and n+k points of R (distance $d_F$ in ref.^[Skokos2016], useful when R comes from delay embedding).

orbitdiagram(ds::DiscreteDynamicalSystem, i, p_index, pvalues; kwargs...)

Compute the orbit diagram (also called bifurcation diagram) of the given system, saving the i variable(s) for parameter values pvalues. The p_index specifies which parameter of the equations of motion is to be changed.

i can be Int or AbstractVector{Int}. If i is Int, returns a vector of vectors. Else it returns vectors of vectors of vectors. Each entry are the points at each parameter value.

Keyword Arguments

• Ttr::Int = 1000 : Transient steps; each orbit is evolved for Ttr first before saving output.
• n::Int = 100 : Amount of points to save for each initial condition.
• dt = 1 : Stepping time. Changing this will give you the orbit diagram of the dt order map.
• u0 = nothing : Specify an initial state. If nothing, the previous state after each parameter is used to seed the new initial condition at the new parameter (with the very first state being the system's state). This makes convergence to the attractor faster, necessitating smaller Ttr. Otherwise u0 can be a standard state, or a vector of states, so that a specific state is used for each parameter.
• ulims = (-Inf, Inf) : only record system states within ulims (only valid if i isa Int).

See also poincaresos and produce_orbitdiagram.

periodicorbits(ds::DiscreteDynamicalSystem,
               o, ics [, λs, indss, singss]; kwargs...) -> FP

Find fixed points FP of order o for the map ds using the algorithm due to Schmelcher & Diakonos^{[Schmelcher1997]}. ics is a collection of initial conditions (container of vectors) to be evolved.

Optional Arguments

The optional arguments λs, indss, singssmust be containers of appropriate values, besides λs which can also be a number. The elements of those containers are passed to: lambdamatrix(λ, inds, sings), which creates the appropriate $\mathbf{\Lambda}_k$ matrix. If these arguments are not given, a random permutation will be chosen for them, with λ=0.001.

Keyword Arguments

• maxiters::Int = 100000 : Maximum amount of iterations an i.c. will be iterated before claiming it has not converged.
• disttol = 1e-10 : Distance tolerance. If the 2-norm of a previous state with the next one is ≤ disttol then it has converged to a fixed point.
• inftol = 10.0 : If a state reaches norm(state) ≥ inftol it is assumed that it has escaped to infinity (and is thus abandoned).
• roundtol::Int = 4 : The found fixed points are rounded to roundtol digits before pushed into the list of returned fixed points FP, if they are not already contained in FP. This is done so that FP doesn't contain duplicate fixed points (notice that this has nothing to do with disttol). Turn this to typemax(Int) to get the full precision of the algorithm.

Description

The algorithm used can detect periodic orbits by turning fixed points of the original map ds to stable ones, through the transformation

\[\mathbf{x}_{n+1} = \mathbf{x}_n + \mathbf{\Lambda}_k\left(f^{(o)}(\mathbf{x}_n) - \mathbf{x}_n\right)\]

with $f$ = eom. The index $k$ counts the various possible $\mathbf{\Lambda}_k$.

Performance Notes

All initial conditions are evolved for all$\mathbf{\Lambda}_k$ which can very quickly lead to long computation times.

permentropy(x, m = 3; τ = 1, base = Base.MathConstants.e)

Compute the permutation entropy^[Brandt2002] of given order m from the x timeseries.

This method is equivalent with

genentropy(x, SymbolicPermutation(; m, τ); base)

permentropy_old(x::AbstractVector, order [, interval=1]; base = Base.MathConstants.e)

Compute the permutation entropy^[Brandt2002] of given order from the x timeseries.

Optionally, interval can be specified to use x[t0:interval:t1] when calculating permutation of the sliding windows between t0 and t1 = t0 + interval * (order - 1).

Optionally use base for the logarithms.

poincaremap!(integ, plane_distance, planecrossing, Tmax, idxs, rootkw)

Low level function that actual performs the algorithm of finding the next crossing of the Poincaré surface of section. Return the state at the section or nothing if evolved for more than Tmax without any crossing.

poincaremap(ds::ContinuousDynamicalSystem, plane, Tmax=1e6; kwargs...) → pmap

Return a map (integrator) that produces iterations over the Poincaré map of the ds. This map is defined as the sequence of points on the Poincaré surface of section. See poincaresos for details on plane and all other kwargs.

You can progress the map one step on the section by calling step!(pmap), which also returns the next state crossing the hyperplane. You can also set the integrator to start from a new state u, which doesn't have to be on the hyperplane, by using reinit!(pmap, u) and then calling step! as normally.

Notice: The argument Tmax exists so that the integrator can terminate instead of being evolved for infinite time, to avoid cases where iteration would continue forever for ill-defined hyperplanes or for convergence to fixed points. Once the system has been evolved for more than Tmax, step!(pmap) will always return nothing.

Example:

ds = Systems.rikitake([0.,0.,0.], μ = 0.47, α = 1.0)
pmap = poincaremap(ds, (3,0.), Tmax=200.)
next_state_on_psos = step!(pmap)
# Change initial condition
reinit!(pmap, [1.0, 0, 0])
next_state_on_psos = step!(pmap)

poincaresos(A::Dataset, plane; kwargs...)

Calculate the Poincaré surface of section of the given dataset with the given plane by performing linear interpolation betweeen points that sandwich the hyperplane.

Argument plane and keywords direction, warning, idxs are the same as above.

poincaresos(ds::ContinuousDynamicalSystem, plane, tfinal = 1000.0; kwargs...)

Calculate the Poincaré surface of section^[Tabor1989] of the given system with the given plane. The system is evolved for total time of tfinal. Return a Dataset of the points that are on the surface of section. See also poincaremap for the map version.

If the state of the system is $\mathbf{u} = (u_1, \ldots, u_D)$ then the equation defining a hyperplane is

\[a_1u_1 + \dots + a_Du_D = \mathbf{a}\cdot\mathbf{u}=b\]

where $\mathbf{a}, b$ are the parameters of the hyperplane.

In code, plane can be either:

• A Tuple{Int, <: Number}, like (j, r) : the plane is defined as when the j variable of the system equals the value r.
• A vector of length D+1. The first D elements of the vector correspond to $\mathbf{a}$ while the last element is $b$.

This function uses ds and higher order interpolation from DifferentialEquations.jl to create a high accuracy estimate of the section. See also produce_orbitdiagram.

Keyword Arguments

• direction = -1 : Only crossings with sign(direction) are considered to belong to the surface of section. Positive direction means going from less than $b$ to greater than $b$.
• idxs = 1:dimension(ds) : Optionally you can choose which variables to save. Defaults to the entire state.
• Ttr = 0.0 : Transient time to evolve the system before starting to compute the PSOS.
• u0 = get_state(ds) : Specify an initial state.
• warning = true : Throw a warning if the Poincaré section was empty.
• rootkw = (xrtol = 1e-6, atol = 1e-6) : A NamedTuple of keyword arguments passed to find_zero from Roots.jl.
• diffeq... : All other extra keyword arguments are propagated into init of DifferentialEquations.jl. See trajectory for examples.

predictability(ds::DynamicalSystem; kwargs...) -> chaos_type, ν, C

Determine whether ds displays strongly chaotic, partially-predictable chaotic or regular behaviour, using the method by Wernecke et al. described in^{[Wernecke2017]}.

Return the type of the behavior, the cross-distance scaling coefficient ν and the correlation coefficient C. Typical values for ν, C and chaos_type are given in Table 2 of^{[Wernecke2017]}:

Keyword Arguments

• Ttr = 200 : Extra "transient" time to evolve the system before sampling from the trajectory. Should be Int for discrete systems.
• T_sample = 1e4 : Time to evolve the system for taking samples. Should be Int for discrete systems.
• n_samples = 500 : Number of samples to take for use in calculating statistics.
• λ_max = lyapunov(ds, 5000) : Value to use for largest Lyapunov exponent for finding the Lyapunov prediction time. If it is less than zero a regular result is returned immediatelly.
• d_tol = 1e-3 : tolerance distance to use for calculating Lyapunov prediction time.
• T_multiplier = 10 : Multiplier from the Lyapunov prediction time to the evaluation time.
• T_max = Inf : Maximum time at which to evaluate trajectory distance. If the internally computed evaluation time is larger than T_max, stop at T_max instead.
• δ_range = 10.0 .^ (-9:-6) : Range of initial condition perturbation distances to use to determine scaling ν.
• diffeq... : Keyword arguments propagated into init of DifferentialEquations.jl. See trajectory for examples. Only valid for continuous systems.

Description

Samples points from a trajectory of the system to be used as initial conditions. Each of these initial conditions is randomly perturbed by a distance δ, and the trajectories for both the original and perturbed initial conditions are computed to the 'evaluation time' T.

The average (over the samples) distance and cross-correlation coefficient of the state at time T is computed. This is repeated for a range of δ (defined by δ_range), and linear regression is used to determine how the distance and cross-correlation scale with δ, allowing for identification of chaos type.

The evaluation time T is calculated as T = T_multiplier*Tλ, where the Lyapunov prediction time Tλ = log(d_tol/δ)/λ_max. This may be very large if the λ_max is small, e.g. when the system is regular, so this internally computed time T can be overridden by a smaller T_max set by the user.

Performance Notes

For continuous systems, it is likely that the maxiters used by the integrators needs to be increased, e.g. to 1e9. This is part of the diffeq kwargs. In addition, be aware that this function does a lot of internal computations. It is operating in a different speed than e.g. lyapunov.

produce_orbitdiagram(
	ds::ContinuousDynamicalSystem, plane, i::Int, p_index, pvalues; kwargs...
)

Produce an orbit diagram (sometimes wrongly called bifurcation diagram) for the i variable(s) of the given continuous system by computing Poincaré surfaces of section using plane for the given parameter values (see poincaresos).

i can be Int or AbstractVector{Int}. If i is Int, returns a vector of vectors. Else it returns a vector of vectors of vectors. Each entry are the points at each parameter value.

Keyword Arguments

• printparams::Bool = false : Whether to print the parameter used during computation in order to keep track of running time.
• direction, warning, Ttr, rootkw, diffeq... : Propagated into poincaresos.
• u0 = nothing : Specify an initial state. If nothing, the previous state after each parameter is used to seed the new initial condition at the new parameter (with the very first state being the system's state). This makes convergence to the attractor faster, necessitating smaller Ttr. Otherwise u0 can be a standard state, or a vector of states, so that a specific state is used for each parameter.

Description

For each parameter, a PSOS reduces the system from a flow to a map. This then allows the formal computation of an "orbit diagram" for the i variable of the system, just like it is done in orbitdiagram.

The parameter change is done as p[p_index] = value taking values from pvalues and thus you must use a parameter container that supports this (either Array, LMArray, dictionary or other).

See also poincaresos, orbitdiagram.

real_to_uint64(data::Dataset{D,T}) where {D, T}

Calculate maximum and minimum value of data to then project the values onto $[0, M - ε(M)]$ where $\epsilon$ is the precision of the used Type and $M$ is the maximum value of the UInt64 type.

reinit_integ_idxs!(integ, y, idxs, u, remidxs)

reinit! given integrator by setting its idxs entries of the state as y, and the remidxs ones as u.

takens_best_estimate(X, εmax, metric = Chebyshev(),εmin = 0) → D_C, D_C_95u, D_C_95l

Use the so-called "Takens' best estimate" ^{[Takens1985][Theiler1988]} method for estimating the correlation dimension D_C and the upper (D_C_95u) and lower (D_C_95l) confidence limit for the given dataset X.

The original formula is

\[D_C \approx \frac{C(\epsilon_\text{max})}{\int_0^{\epsilon_\text{max}}(C(\epsilon) / \epsilon) \, d\epsilon}\]

where $C$ is the correlationsum and $\epsilon_\text{max}$ is an upper cutoff. Here we use the later expression

\[D_C \approx - \frac{1}{\eta},\quad \eta = \frac{1}{(N-1)^*}\sum_{[i, j]^*}\log(||X_i - X_j|| / \epsilon_\text{max})\]

where the sum happens for all $i, j$ so that $i < j$ and $||X_i - X_j|| < \epsilon_\text{max}$. In the above expression, the bias in the original paper has already been corrected, as suggested in ^{[Borovkova1999]}.

The confidence limits are estimated from the log-likelihood function by finding the values of D_C where the function has fallen by 2 from its maximum, see e.g. ^[Barlow] chapter 5.3 Because the CLT does not apply (no independent measurements), the limits are not neccesarily symmetric.

According to ^{[Borovkova1999]}, introducing a lower cutoff εmin can make the algorithm more stable (no divergence), this option is given but defaults to zero.

If X comes from a delay coordinates embedding of a timseries x, a recommended value for $\epsilon_\text{max}$ is std(x)/4.

testchaos01(φ::Vector [, cs, N0]) -> chaotic?

Perform the so called "0-1" test for chaos introduced by Gottwald and Melbourne^{[Gottwald2016]} on the timeseries φ. Return true if φ is chaotic, false otherwise.

Description

This method tests if the given timeseries is chaotic or not by transforming it into a two-dimensional diffusive process. If the timeseries is chaotic, the mean square displacement of the process grows as sqrt(length(φ)), while it stays constant if the timeseries is regular. The implementation here computes K, the correlation coefficient (median of Kc for c ∈ cs), and simply checks if K > 0.5.

If you want to access the various Kc you should call the method testchaos01(φ, c::Real, N0) which returns Kc.

cs defaults to 3π/5*rand(10) + π/4 and N0, the length of the two-dimensional process, is N0 = length(φ)/10.

Notice that for data sampled from continous dynamical systems, some care must be taken regarding the values of cs, see^{[Gottwald2016]}.

tipping_probabilities(basins_before, basins_after) → P

Return the tipping probabilities of the computed basins before and after a change in the system parameters (or time forcing), according to the definition of ^..

The input basins are integer-valued arrays, where the integers enumerate the attractor. They can be of any dimensionality provided that size(basins_before) == size(basins_after). Typically they are 2D, as the output of [basins2D](@ref) or [basinsgeneral`](@ref)

Description

Let $\mathcal{B}_i(p)$ denote the basin of attraction of attractor $A_i$ at parameter(s) $p$. Kaszás et al^{[Kaszás2019]} define the tipping probability from $A_i$ to $A_j$, given a parameter change in the system of $p_- \to p_+$, as

\[P(A_i \to A_j | p_- \to p_+) = \frac{|\mathcal{B}_j(p_+) \cap \mathcal{B}_i(p_-)|}{|\mathcal{B}_i(p_-)|}\]

where $|\cdot|$ is simply the volume of the enclosed set. The value of $P(A_i \to A_j | p_- \to p_+)$ is P[i, j]. The equation describes something quite simple: what is the overlap of the basin of attraction of $A_i$ at $p_-$ with that of the attractor $A_j$ at $p_+$. If basins_before, basins_after contain values of -1, corresponding to trajectories that diverge, this is considered as the last attractor of the system in P.

transit_return(exits, entries) → transit, return

Convert the output of exit_entry_times to the transit and return times.

uncertainty_exponent(xg, yg, basins::Matrix; kwargs...) -> ε, f_ε ,α

Estimate the uncertainty exponent^{[Grebogi1983]} of the basins of attraction. This exponent is related to the final state sensitivity of the trajectories in the phase space. An exponent close to 1 means basins with smooth boundaries whereas an exponent close to 0 represent complety fractalized basins, also called riddled basins.

xg, yg are 1-dim ranges that defines the grid of the initial conditions. basins is the matrix containing the information of the basin, i.e. the output of basins_2D or basins_general.

The output f_ε is a vector with the fraction of the balls of radius ε (in pixels) that contain at least two initial conditions that lead to different attractors. The ouput α is the estimation of the uncertainty exponent of the basins of attraction by fitting a line in the log.(f_ε) vs log.(ε) curve, however it is recommended to analyze the curve directly for more accuracy.

Keyword arguments

• precision = 1e-5 is the variance of the estimator of the uncertainty function. Values between 1e-7 and 1e-5 brings reasonable results.
• max_ε = floor(Int, length(xg)/20) is the maximum size in pixels of the ball to test.

Description

A phase space with a fractal boundary may cause a uncertainty on the final state of the dynamical system for a given initial condition. A measure of this final state sensitivity is the uncertainty exponent. The algorithm probes the basin of attraction with balls of size ε at random. If there are a least two initial conditions that lead to different attractors, a ball is tagged "uncertain". f_ε is the fraction of "uncertain balls" to the total number of tries in the basin. In analogy to the fractal dimension, there is a scaling law between, f_ε ~ ε^α. The number that characterizes this scaling is called the uncertainty exponent α.

Notice that the uncertainty exponent and the box counting dimension of the boundary are related. We have Δ₀ = 2 - α where Δ₀ is the box counting dimension, see Fractal Dimension

Return the signed distance of state to ε-ball (negative means inside ball)