Attractors, Basins, Tipping Points

AttractorMapper(ds::GeneralizedDynamicalSystem, args...; kwargs...) → mapper

Subtypes of AttractorMapper are structures that map initial conditions of ds to attractors. Currently available mapping methods:

• AttractorsViaProximity
• AttractorsViaRecurrences
• AttractorsViaFeaturizing

All AttractorMapper subtypes can be used with basins_fractions or basins_of_attraction.

In addition, some mappers can be called as a function of an initial condition:

label = mapper(u0)

and this will on the fly compute and return the label of the attractor u0 converges at. The mappers that can do this are:

• AttractorsViaProximity
• AttractorsViaRecurrences

AttractorsViaProximity(ds::DynamicalSystem, attractors::Dict [, ε]; kwargs...)

Map initial conditions to attractors based on whether the trajectory reaches ε-distance close to any of the user-provided attractors. They have to be in a form of a dictionary mapping attractor labels to Datasets containing the attractors.

The system gets stepped, and at each step the minimum distance to all attractors is computed. If any of these distances is < ε, then the label of the nearest attractor is returned. If an ε::Real is not provided by the user, a value is computed automatically as half of the minimum distance between all attractors. This operation can be expensive for large attractor datasets.

Because in this method the attractors are already known to the user, the method can also be called supervised.

Keywords

• Ttr = 100: Transient time to first evolve the system for before checking for proximity.
• Δt = 1: Integration step time (only valid for continuous systems).
• horizon_limit = 1e3: If the maximum distance of the trajectory from any of the given attractors exceeds this limit, it is assumed that the trajectory diverged (gets labelled as -1).
• mx_chk_lost = 1000: If the integrator has been stepped this many times without coming ε-near to any attractor, it is assumed that the trajectory diverged (gets labelled as -1).
• diffeq = NamedTuple(): Keywords propagated to DifferentialEquations.jl (only valid for continuous systems).

AttractorsViaRecurrences(ds::GeneralizedDynamicalSystem, grid::Tuple; kwargs...)

Map initial conditions to attractors by identifying attractors on the fly based on recurrences in the state space, as outlined by Datseris & Wagemakers^{[Datseris2022]}. Works for any case encapsulated by GeneralizedDynamicalSystem.

grid is a tuple of ranges partitioning the state space so that a finite state machine can operate on top of it. For example grid = (xg, yg) where xg = yg = range(-5, 5; length = 100) for a two-dimensional system. The grid has to be the same dimensionality as the state space, use a projected_integrator if you want to search for attractors in a lower dimensional subspace.

Keyword Arguments

• Δt: Approximate time step of the integrator, which is 1 for discrete systems. For continuous systems, an automatic value is calculated using automatic_Δt_basins.
• diffeq = NamedTuple(): Keyword arguments propagated to integrator. Only valid for ContinuousDynamicalSystem. It is recommended to choose high accuracy solvers for this application, e.g. diffeq = (alg=Vern9(), reltol=1e-9, abstol=1e-9).
• mx_chk_att = 2: A parameter that sets the maximum checks of consecutives hits of an attractor before deciding the basin of the initial condition.
• mx_chk_hit_bas = 10: Maximum check of consecutive visits of the same basin of attraction. This number can be increased for higher accuracy.
• mx_chk_fnd_att = 100: Maximum check of unnumbered cell before considering we have an attractor. This number can be increased for higher accuracy.
• mx_chk_loc_att = 100: Maximum check of consecutive cells marked as an attractor before considering that we have all the available pieces of the attractor.
• mx_chk_lost = 20: Maximum check of iterations outside the defined grid before we consider the orbit lost outside. This number can be increased for higher accuracy.
• horizon_limit = 1e6: If the norm of the integrator state reaches this limit we consider that the orbit diverges.

Description

An initial condition given to an instance of AttractorsViaRecurrences is iterated based on the integrator corresponding to ds. A recurrence in the state space means that the trajectory has converged to an attractor. This is the basis for finding attractors.

A finite state machine (FSM) follows the trajectory in the state space, and constantly maps it to the given grid. The FSM decides when an initial condition has successfully converged into an attractor. An array, internally called "basins", stores the state of the FSM on the grid, according to the indexing system described in ^{[Datseris2022]}. As the system is integrated more and more, the information of the "basins" becomes richer and richer with more identified attractors or with grid cells that belong to basins of already found attractors. Notice that only in the special method basins_of_attraction(mapper::AttractorsViaRecurrences) the information of the attraction or exit basins is utilized. In other functions like basins_fractions only the attractor locations are utilized.

The iteration of a given initial condition continues until one of the following happens:

1. The trajectory hits mx_chk_fnd_att times in a row grid cells previously visited: it is considered that an attractor is found and is labelled with a new number.
2. The trajectory hits an already identified attractor mx_chk_att consecutive times: the initial condition is numbered with the attractor's number.
3. The trajectory hits a known basin mx_chk_hit_bas times in a row: the initial condition belongs to that basin and is numbered accordingly.
4. The trajectory spends mx_chk_lost steps outside the defined grid or the norm of the integrator state becomes > than horizon_limit: the initial condition is set to -1.

automatic_Δt_basins(integ, grid; N = 5000) → Δt

Calculate an optimal Δt value for basins_of_attraction. This is done by evaluating the dynamic rule f (vector field) at N randomly chosen points of the grid. The average f is then compared with the diagonal length of a grid cell and their ratio provides Δt.

Notice that Δt should not be too small which happens typically if the grid resolution is high. It is okay for basins_of_attraction if the trajectory skips a few cells. But if Δt is too small the default values for all other keywords such as mx_chk_hit_bas need to be increased drastically.

Also, Δt that is smaller than the internal step size of the integrator will cause a performance drop.

AttractorsViaFeaturizing(ds::DynamicalSystem, featurizer::Function; kwargs...) → mapper

Initialize a mapper that maps initial conditions to attractors using the featurizing and clustering method of ^{[Stender2021]}. See AttractorMapper for how to use the mapper.

featurizer is a function that takes as an input an integrated trajectory A::Dataset and the corresponding time vector t and returns a Vector{<:Real} of features describing the trajectory.

Keyword arguments

Integration

• T=100, Ttr=100, Δt=1, diffeq=NamedTuple(): Propagated to trajectory.

Feature extraction and classification

• attractors_ic = nothing Enables supervised version, see below. If given, must be a Dataset of initial conditions each leading to a different attractor.
• min_neighbors = 10: (unsupervised method only) minimum number of neighbors (i.e. of similar features) each feature needs to have in order to be considered in a cluster (fewer than this, it is labeled as an outlier, -1).
• clust_method_norm=Euclidean(): metric to be used in the clustering.
• clustering_threshold = 0.0: Maximum allowed distance between a feature and the cluster center for it to be considered inside the cluster. Only used when clust_method = "kNN_thresholded".
• clust_method = clustering_threshold > 0 ? "kNN_thresholded" : "kNN": (supervised method only) which clusterization method to apply. If "kNN", the first-neighbor clustering is used. If "kNN_thresholded", a subsequent step is taken, which considers as unclassified (label -1) the features whose distance to the nearest template is above the clustering_threshold.
• rescale_features = true: (unsupervised method): if true, rescale each dimension of the

extracted features separately into the range [0,1].

• optimal_radius_method = silhouettes (unsupervised method): the method used to determine

the optimal radius for clustering features in the unsupervised method. The silhouettes method chooses the radius that maximizes the average silhouette values of clusters, and is an iterative optimization procedure that may take some time to execute. The elbow method chooses the the radius according to the elbow (knee, highest-derivative method) (see optimal_radius_dbscan_elbow for details), and is quicker though possibly leads to worse clustering.

Description

The trajectory X of each initial condition is transformed into a vector of features. Each feature is a number useful in characterizing the attractor the initial condition ends up at, and distinguishing it from other attrators. Example features are the mean or standard deviation of one of the of the timeseries of the trajectory, the entropy of the first two dimensions, the fractal dimension of X, or anything else you may fancy. The vectors of features are then used to identify to which attractor each trajectory belongs (i.e. in which basin of attractor each initial condition is in). The method thus relies on the user having at least some basic idea about what attractors to expect in order to pick the right features, in contrast to AttractorsViaRecurrences.

The algorithm of^{[Stender2021]} that we use has two versions to do this. If the attractors are not known a-priori the unsupervised versions should be used. Here, the vectors of features of each initial condition are mapped to an attractor by analysing how the features are clustered in the feature space. Using the DBSCAN algorithm, we identify these clusters of features, and consider each cluster to represent an attractor. Features whose attractor is not identified are labeled as -1. If each feature spans different scales of magnitude, rescaling them into the same [0,1] interval can bring significant improvements in the clustering in case the Euclidean distance metric is used.

In the supervised version, the attractors are known to the user, who provides one initial condition for each attractor using the attractors_ic keyword. The algorithm then evolves these initial conditions, extracts their features, and uses them as templates representing the attrators. Each trajectory is considered to belong to the nearest template (based on the distance in feature space). Notice that the functionality of this version is similar to AttractorsViaProximity. Generally speaking, the AttractorsViaProximity is superior. However, if the dynamical system has extremely high-dimensionality, there may be reasons to use the supervised method of this featurizing algorithm instead.

Parallelization note

The trajectories in this method are integrated in parallel using Threads. To enable this, simply start Julia with the number of threads you want to use.

basins_fractions(mapper::AttractorMapper, ics::Union{Dataset, Function}; kwargs...)

Approximate the state space fractions fs of the basins of attraction of a dynamical stystem by mapping initial conditions to attractors using mapper (which contains a reference to a GeneralizedDynamicalSystem). The fractions are simply the ratios of how many initial conditions ended up at each attractor.

Initial conditions to use are defined by ics. It can be:

• a Dataset of initial conditions, in which case all are used.
• a 0-argument function ics() that spits out random initial conditions. Then N random initial conditions are chosen. See statespace_sampler to generate such functions.

The returned arguments are fs. If ics is a Dataset then the labels of each initial and roughly approximated attractors are also returned: fs, labels, attractors.

The output fs is a dictionary whose keys are the labels given to each attractor (always integers enumerating the different attractors), and the values are their respective fractions. The label -1 is given to any initial condition where mapper could not match to an attractor (this depends on the mapper type). attractors has the same structure, mapping labels to Datasets.

See AttractorMapper for all possible mapper types.

Keyword arguments

• N = 1000: Number of random initial conditions to generate in case ics is a function.
• show_progress = true: Display a progress bar of the process.

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

Calculate the state space 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 (i.e., like the output of basins_of_attraction). 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_of_attraction(mapper::AttractorMapper, grid::Tuple) → basins, attractors

Compute the full basins of attraction as identified by the given mapper, which includes a reference to a GeneralizedDynamicalSystem and return them along with (perhaps approximated) found attractors.

grid is a tuple of ranges defining the grid of initial conditions that partition the state space into boxes with size the step size of each range. For example, grid = (xg, yg) where xg = yg = range(-5, 5; length = 100). The grid has to be the same dimensionality as the state space expected by the integrator/system used in mapper. E.g., a projected_integrator could be used for lower dimensional projections, etc. A special case here is a poincaremap with plane being Tuple{Int, <: Real}. In this special scenario the grid can be one dimension smaller than the state space, in which case the partitioning happens directly on the hyperplane the Poincaré map operates on.

basins_of_attraction function is a convenience 5-lines-of-code wrapper which uses the labels returned by basins_fractions and simply assings them to a full array corresponding to the state space partitioning indicated by grid.

basins_of_attraction(mapper::AttractorsViaRecurrences; show_progress = true)

This is a special method of basins_of_attraction that using recurrences does exactly what is described in the paper by Datseris & Wagemakers^{[Datseris2022]}. By enforcing that the internal grid of mapper is the same as the grid of initial conditions to map to attractors, the method can further utilize found exit and attraction basins, making the computation faster as the grid is processed more and more.

statespace_sampler(rng = Random.GLOBAL_RNG; kwargs...) → sampler, isinside

Convenience function that creates two functions. sampler is a 0-argument function that generates random points inside a state space region defined by the keywords. isinside is a 1-argument function that decides returns true if the given state space point is inside that region.

The regions can be:

• Rectangular box, with edges min_bounds and max_bounds. The sampling of the points inside the box is decided by the keyword method which can be either "uniform" or "multgauss".
• Sphere, of spheredims dimensions, radius radius and centered on center.

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 (i.e, output of basins_of_attraction). Typically the +,- mean after and before some change of parameter for a system.

In basins_of_attraction 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.

basins_fractal_dimension(basins; kwargs...) -> V_ε, N_ε ,d

Estimate the Fractal Dimension d of the boundary between basins of attraction using the box-counting algorithm.

The output N_ε is a vector with the number of the balls of radius ε (in pixels) that contain at least two initial conditions that lead to different attractors. V_ε is a vector with the corresponding size of the balls. The ouput d is the estimation of the box-counting dimension of the boundary by fitting a line in the log.(N_ε) vs log.(1/V_ε) curve. However it is recommended to analyze the curve directly for more accuracy.

Keyword arguments

• range_ε = 2:maximum(size(basins))÷20 is the range of sizes of the ball to test (in pixels).

Description

It is the implementation of the popular algorithm of the estimation of the box-counting dimension. The algorithm search for a covering the boundary with N_ε boxes of size ε in pixels.

basin_entropy(basins, ε = 20) -> Sb, Sbb

This algorithm computes the basin entropy Sb of the basins of attraction. First, the input basins is divided regularly into n-dimensional boxes of side ε (along all dimensions). Then Sb is simply the average of the Gibbs entropy computed over these boxes. The function returns the basin entropy Sb as well as the boundary basin entropy Sbb. The later is the average of the entropy only for boxes that contains at least two different basins, that is, for the boxes on the boundary.

The basin entropy is a measure of the uncertainty on the initial conditions of the basins. It is maximum at the value log(n_att) being n_att the number of attractors. In this case the boundary is intermingled: for a given initial condition we can find another initial condition that lead to another basin arbitriraly close. It provides also a simple criterion for fractality: if the boundary basin entropy Sbb is above log(2) then we have a fractal boundary. It doesn't mean that basins with values below cannot have a fractal boundary, for a more precise test see basins_fractal_test. An important feature of the basin entropy is that it allows comparisons between different basins using the same box size ε.

basins_fractal_test(basins; ε = 20, Ntotal = 1000) -> test_res, Sbb

This is an automated test to decide if the boundary of the basins has fractal structures. The bottom line is to look at the basins with a magnifier of size ε at random in basins. If what we see in the magnifier looks like a smooth boundary (in average) we decide that the boundary is smooth. If it is not smooth we can say that at the scale ε we have structures, i.e., it is fractal.

In practice the algorithm computes the boundary basin entropy Sbb basin_entropy for Ntotal random balls of radius ε. If the computed value is equal to theoretical value of a smooth boundary (taking into account statistical errors and biases) then we decide that we have a smooth boundary. Notice that the response test_res may depend on the chosen ball radius ε. For larger size, we may observe structures for smooth boundary and we obtain a different answer.

The output test_res is a symbol describing the nature of the basin and the output Sbb is the estimated value of the boundary basin entropy with the sampling method.

^[Puy2021] Andreu Puy, Alvar Daza, Alexandre Wagemakers, Miguel A. F. Sanjuán. A test for fractal boundaries based on the basin entropy. Commun Nonlinear Sci Numer Simulat, 95, 105588, 2021.

Keyword arguments

• ε = 20: size of the ball for the test of basin. The result of the test may change with the size.
• Ntotal = 1000: number of balls to test in the boundary for the computation of Sbb

uncertainty_exponent(basins; kwargs...) -> ε, N_ε ,α

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 completely fractalized basins, also called riddled basins.

The output N_ε is a vector with the number 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 using the box-counting dimension of the boundary by fitting a line in the log.(N_ε) vs log.(1/ε) curve. However it is recommended to analyze the curve directly for more accuracy.

Keyword arguments

• range_ε = 2:maximum(size(basins))÷20 is the range of sizes of the ball to test (in pixels).

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 Δ₀ = D - α where Δ₀ is the box counting dimension computed with basins_fractal_dimension and D is the dimension of the phase space. The algorithm first estimates the box counting dimension of the boundary and returns the uncertainty exponent.

An obstruction to predictability, Physics Letters A, 99, 9, 1983

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^{[Kaszás2019]}.

The input basins are integer-valued arrays, where the integers enumerate the attractor, e.g. the output of basins_of_attraction.

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.