# Attractors, Basins, Tipping Points

## Finding Attractors

DynamicalSystems.jl defines a generic interface for finding attractors of dynamical systems. One first decides the instance of `GeneralizedDynamicalSystem`

they need. Then, an instance of `AttractorMapper`

is created from this dynamical system. This `mapper`

instance can be used to compute e.g., `basins_of_attraction`

, and the output can be further analyzed to get e.g., the `basin_entropy`

.

The example 2D basins of 4D system contains most of the functionality documented in this page.

`ChaosTools.AttractorMapper`

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

Subtypes of `AttractorMapper`

are structures that map initial conditions of `ds`

to attractors. Currently available mapping methods:

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:

`ChaosTools.AttractorsViaProximity`

— Type`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 `Dataset`

s 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).

`ChaosTools.AttractorsViaRecurrences`

— Type`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:

- 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. - The trajectory hits an already identified attractor
`mx_chk_att`

consecutive times: the initial condition is numbered with the attractor's number. - 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. - 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.

`ChaosTools.automatic_Δt_basins`

— Function`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.

`ChaosTools.AttractorsViaFeaturizing`

— Type`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 of attraction

Calculating basins of attraction, or their state space fractions, can be done with the functions `basins_fractions`

and `basins_of_attraction`

. See also the convenience `statespace_sampler`

and `match_attractors`

.

`ChaosTools.basins_fractions`

— Function`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 `Dataset`

s.

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.

`ChaosTools.basins_of_attraction`

— Function`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.

`ChaosTools.statespace_sampler`

— Function`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`

.

`ChaosTools.match_attractors!`

— Function`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 `Dataset`

s 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.

## Final state sensitivity / fractal boundaries

Several functions are provided related with analyzing the fractality of the boundaries of the basins of attraction:

`ChaosTools.basins_fractal_dimension`

— Function`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.

`ChaosTools.basin_entropy`

— Function`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 `ε`

.

`ChaosTools.basins_fractal_test`

— Function`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`

`ChaosTools.uncertainty_exponent`

— Function`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 points

`ChaosTools.tipping_probabilities`

— Function`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`

.

## Discrete system example

```
using DynamicalSystems
function newton_map(z, p, n)
z1 = z[1] + im*z[2]
dz1 = newton_f(z1, p[1])/newton_df(z1, p[1])
z1 = z1 - dz1
return SVector(real(z1), imag(z1))
end
newton_f(x, p) = x^p - 1
newton_df(x, p)= p*x^(p-1)
ds = DiscreteDynamicalSystem(newton_map, [0.1, 0.2], [3.0])
xg = yg = range(-1.5, 1.5; length = 400)
mapper = AttractorsViaRecurrences(ds, (xg, yg))
basins, attractors = basins_of_attraction(mapper; show_progress = false)
basins
```

```
400×400 Matrix{Int16}:
1 1 1 1 1 1 1 1 1 1 1 1 1 … 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 … 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2
⋮ ⋮ ⋮ ⋱ ⋮ ⋮
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 … 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
```

`attractors`

```
Dict{Int16, Dataset{2, Float64}} with 3 entries:
2 => 2-dimensional Dataset{Float64} with 1 points
3 => 2-dimensional Dataset{Float64} with 1 points
1 => 2-dimensional Dataset{Float64} with 1 points
```

Now let's plot this as a heatmap

```
using CairoMakie
# Set up some code for plotting attractors
function scatter_attractors!(ax, attractors)
for k ∈ keys(attractors)
x, y = columns(attractors[k])
scatter!(ax, attractors[k].data;
color = Cycled(k),
strokewidth = 3, strokecolor = :white
)
end
end
generate_cmap(n) = cgrad(Main.COLORS[1:n], n; categorical = true)
ids = sort!(unique(basins))
cmap = generate_cmap(length(ids))
fig, ax = heatmap(xg, yg, basins;
colormap = cmap, colorrange = (ids[1] - 0.5, ids[end]+0.5),
)
scatter_attractors!(ax, attractors)
fig
```

## Stroboscopic map example

This example targets periodically driven 2D continuous dynamical systems, like the Duffing oscillator:

```
using DynamicalSystems
ω=1.0; f = 0.2
ds = Systems.duffing([0.1, 0.25]; ω, f, d = 0.15, β = -1)
using OrdinaryDiffEq
diffeq = (alg = Vern9(), reltol = 1e-9)
smap = stroboscopicmap(ds, 2π/ω; diffeq)
```

```
Iterator of the stroboscopic map
rule f: duffing_rule
Period: 6.283185307179586
```

For stroboscopic maps, we strongly recommend using a higher precision integrator from OrdinaryDiffEq.jl. Now we can compute the basins of `smap`

just like before.

```
xg = yg = range(-2.2, 2.2; length = 200)
mapper = AttractorsViaRecurrences(smap, (xg, yg); diffeq)
basins, attractors = basins_of_attraction(mapper; show_progress = false)
basins
```

```
200×200 Matrix{Int16}:
1 2 1 1 1 1 1 1 1 2 2 1 2 … 1 1 2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 2 1 2 2 1 1 1 2 1 2 2 2 1 1 2 2 2 2 1 2
1 1 2 1 1 1 1 2 2 1 2 1 2 1 1 1 1 1 1 2 2 2 1 1 2
1 1 1 2 2 2 1 1 1 2 1 2 2 2 2 2 1 1 1 1 2 2 2 2 2
2 2 1 2 2 2 1 2 2 2 2 2 2 1 1 1 1 2 2 1 1 1 1 2 2
2 2 2 2 2 2 2 2 1 1 1 1 2 … 1 1 1 1 1 1 1 2 2 1 1 1
2 2 2 1 1 2 1 1 2 2 2 2 1 1 1 1 1 2 2 1 1 1 1 1 1
2 2 2 2 1 2 2 2 1 1 2 1 2 2 1 1 1 1 1 2 1 1 1 2 1
2 2 2 2 1 2 2 2 1 2 2 2 1 1 2 2 1 1 1 1 1 2 2 1 2
1 1 2 2 1 1 2 1 1 1 1 1 1 1 1 1 2 2 1 2 1 1 1 2 1
⋮ ⋮ ⋮ ⋱ ⋮ ⋮
2 1 2 1 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2
1 2 2 1 2 2 2 2 2 2 1 1 1 1 1 1 1 1 2 2 1 1 2 2 1
1 2 1 1 2 2 2 1 1 2 2 2 2 1 1 1 2 2 2 1 2 2 1 2 2
1 1 2 2 2 1 1 2 2 2 2 1 1 1 1 1 2 1 2 1 1 2 1 1 1
2 2 2 2 2 1 2 1 1 1 2 1 2 … 1 1 1 2 2 2 2 2 2 1 2 2
2 1 1 2 2 1 1 1 2 2 2 1 1 1 1 2 1 2 2 2 2 2 1 2 1
2 2 2 2 1 2 2 1 2 1 2 2 2 2 2 2 2 2 2 2 2 2 1 2 1
1 2 2 2 2 2 1 1 2 1 2 2 1 2 2 2 1 1 1 2 2 1 1 1 2
2 1 2 2 2 2 2 2 2 2 2 1 2 1 1 1 2 1 1 1 2 1 2 1 2
```

And visualize the result as a heatmap, scattering the found attractors via scatter.

```
using CairoMakie
ids = sort!(unique(basins))
cmap = generate_cmap(length(ids))
fig, ax = heatmap(xg, yg, basins;
colormap = cmap, colorrange = (ids[1] - 0.5, ids[end]+0.5),
)
scatter_attractors!(ax, attractors)
fig
```

## 2D basins of 4D system

In this section we will calculate the basins of attraction of the four-dimensional magnetic pendulum. We know that the attractors of this system are all individual fixed points on the (x, y) plane so we will only compute the basins there. We can also use this opportunity to highlight a different method, the `AttractorsViaProximity`

which works when we already know where the attractors are.

### Computing the basins

```
ds = Systems.magnetic_pendulum(d=0.2, α=0.2, ω=0.8, N=3)
psys = projected_integrator(ds, [1, 2], [0.0, 0.0])
```

```
Integrator of a projected system
rule f: MagneticPendulum
projection: [1, 2]
dimension: 2
complete state: [0.0, 0.0]
```

For this systems we know the attractors are close to the magnet positions, so we can just do

```
attractors = Dict(i => Dataset([ds.f.magnets[i]]) for i in 1:3)
mapper = AttractorsViaProximity(psys, attractors)
```

```
AttractorsViaProximity
rule f: MagneticPendulum
type: ProjectedIntegrator
ε: 0.8660254037844386
Δt: 1
Ttr: 100
attractors: Dict{Int64, Dataset{2, Float64}} with 3 entries:
2 => 2-dimensional Dataset{Float64} with 1 points
3 => 2-dimensional Dataset{Float64} with 1 points
1 => 2-dimensional Dataset{Float64} with 1 points
```

and as before, get the basins of attraction

```
xg = yg = range(-4, 4; length=150)
basins, = basins_of_attraction(mapper, (xg, yg); show_progress = false)
ids = sort!(unique(basins))
cmap = generate_cmap(length(ids))
fig, ax = heatmap(xg, yg, basins;
colormap = cmap, colorrange = (ids[1] - 0.5, ids[end]+0.5),
)
scatter_attractors!(ax, attractors)
fig
```

### Computing the uncertainty exponent

Let's now calculate the `uncertainty_exponent`

for this system as well. The calculation is straightforward:

```
ε, f_ε, α = uncertainty_exponent(basins)
fig, ax = lines(log.(ε), log.(f_ε))
ax.title = "α = $(round(α; digits=3))"
fig
```

The actual uncertainty exponent is the slope of the curve (α) and indeed we get an exponent near 0 as we know a-priory the basins have fractal boundaries for the magnetic pendulum.

### Computing the tipping probabilities

We will compute the tipping probabilities using the magnetic pendulum's example as the "before" state. For the "after" state we will change the `γ`

parameter of the third magnet to be so small, its basin of attraction will virtually disappear. As we don't know *when* the basin of the third magnet will disappear, we switch the attractor finding algorithm back to `AttractorsViaRecurrences`

.

```
ds = Systems.magnetic_pendulum(d=0.2, α=0.2, ω=0.8, N=3, γs = [1.0, 1.0, 0.1])
psys = projected_integrator(ds, [1, 2], [0.0, 0.0]; diffeq = (reltol = 1e-9,))
mapper = AttractorsViaRecurrences(psys, (xg, yg); Δt = 1)
basins_after, attractors_after = basins_of_attraction(
mapper; show_progress = false
)
# matching attractors is important!
match_attractors!(basins, attractors, basins_after, attractors_after)
# now plot
ids = sort!(unique(basins_after))
cmap = generate_cmap(length(ids))
fig, ax = heatmap(xg, yg, basins_after;
colormap = cmap, colorrange = (ids[1] - 0.5, ids[end]+0.5),
)
scatter_attractors!(ax, attractors_after)
fig
```

`P = tipping_probabilities(basins, basins_after)`

```
3×2 Matrix{Float64}:
0.5 0.5
0.450115 0.549885
0.550155 0.449845
```

As you can see `P`

has size 3×2, as after the change only 2 attractors have been identified in the system (3 still exist but our state space discretization isn't accurate enough to find the 3rd because it has such a small basin). Also, the first row of `P`

is 50% probability to each other magnet, as it should be due to the system's symmetry.

### Computing the basin fractions via featurizing

We can also compute the fraction of initial conditions that go to each of the three attractors for the magnetic pendulum via featurizing using `AttractorsViaFeaturizing`

because here it is trivial to find the appropriate features that distinguish the attractors: simply the final value of the x and y coordinates of the trajectory.

Let's also just use `basins_fractions`

with random sampling instead of computing a full array of basins of attractions, just for some variety.

So let's define a random sampler and the featurizing function

```
using Random
rng = Random.MersenneTwister(1234)
sampler, _ = statespace_sampler(rng; min_bounds=[-4, -4], max_bounds=[4, 4])
featurizer(A, t) = A[end]
```

`featurizer (generic function with 1 method)`

which gives

```
ds = Systems.magnetic_pendulum(d=0.2, α=0.2, ω=0.8, N=3)
psys = projected_integrator(ds, [1, 2], [0.0, 0.0]; diffeq = (reltol = 1e-9,))
mapper = AttractorsViaFeaturizing(psys, featurizer)
fs = basins_fractions(mapper, sampler; N = 1000, show_progress = false)
```

```
Dict{Int64, Float64} with 3 entries:
2 => 0.341
3 => 0.306
1 => 0.353
```

As it should, we get about 33% fraction for each attractor.

## 3D basins

To showcase the true power of `AttractorsViaRecurrences`

we need to use a system whose attractors span higher-dimensional space. An example is

`ds = Systems.thomas_cyclical(b = 0.1665)`

```
3-dimensional continuous dynamical system
state: [1.0, 0.0, 0.0]
rule f: thomas_rule
in-place? false
jacobian: thomas_jacob
parameters: [0.1665]
```

which, for this parameter, contains 5 coexisting attractors. 3 of them are entangled periodic orbits that span across all three dimensions, and the remaining 2 are fixed points.

To compute the basins we define a three-dimensional grid and call on it `basins_of_attraction`

.

```
# This computation takes about an hour
xg = yg = zg = range(-6.0, 6.0; length = 251)
mapper = AttractorsViaRecurrences(ds, (xg, yg, zg))
basins, attractors = basins_of_attraction(mapper)
attractors
```

```
Dict{Int16, Dataset{3, Float64}} with 5 entries:
5 => 3-dimensional Dataset{Float64} with 1 points
4 => 3-dimensional Dataset{Float64} with 379 points
6 => 3-dimensional Dataset{Float64} with 1 points
2 => 3-dimensional Dataset{Float64} with 538 points
3 => 3-dimensional Dataset{Float64} with 537 points
1 => 3-dimensional Dataset{Float64} with 1 points
```

The basins of attraction are very complicated. We can try to visualize them by animating the 2D slices at each z value, to obtain:

Then, we visualize the attractors to obtain:

In the animation above, the scattered points are the attractor values the function `AttractorsViaRecurrences`

found by itself. Of course, for the periodic orbits these points are incomplete. Once the function's logic understood we are on an attractor, it stops computing. However, we also simulated lines, by evolving initial conditions colored appropriately with the basins output.

The animation was produced with the code:

```
using GLMakie
fig = Figure()
display(fig)
ax = fig[1,1] = Axis3(fig; title = "found attractors")
cmap = cgrad(:dense, 6; categorical = true)
for i in keys(attractors)
tr = attractors[i]
markersize = length(attractors[i]) > 10 ? 2000 : 6000
marker = length(attractors[i]) > 10 ? :circle : :rect
scatter!(ax, columns(tr)...; markersize, marker, transparency = true, color = cmap[i])
j = findfirst(isequal(i), basins)
x = xg[j[1]]
y = yg[j[2]]
z = zg[j[3]]
tr = trajectory(ds, 100, SVector(x,y,z); Ttr = 100)
lines!(ax, columns(tr)...; linewidth = 1.0, color = cmap[i])
end
a = range(0, 2π; length = 200) .+ π/4
record(fig, "cyclical_attractors.mp4", 1:length(a)) do i
ax.azimuth = a[i]
end
```

## Poincaré map example

In the previous example we saw that this system has periodic attractors. In the Poincaré map these periodic attractors become points. We can use any instance of `AttractorMapper`

and `poincaremap`

to find basins of attraction on the Poincaré surface of section.

```
ds = Systems.thomas_cyclical(b = 0.1665)
xg = yg = range(-6.0, 6.0; length = 100)
pmap = poincaremap(ds, (3, 0.0);
rootkw = (xrtol = 1e-8, atol = 1e-8), diffeq = (reltol=1e-9,)
)
mapper = AttractorsViaRecurrences(pmap, (xg, yg))
basins, attractors = basins_of_attraction(mapper; show_progress = false)
attractors
```

```
Dict{Int16, Dataset{3, Float64}} with 3 entries:
2 => 3-dimensional Dataset{Float64} with 1 points
3 => 3-dimensional Dataset{Float64} with 3 points
1 => 3-dimensional Dataset{Float64} with 1 points
```

`attractors[1]`

```
3-dimensional Dataset{Float64} with 1 points
1.83899 -4.15575 -2.08263e-11
```

Looks good so far, but let's plot it as well:

```
ids = sort!(unique(basins))
cmap = generate_cmap(length(ids))
fig, ax = heatmap(xg, yg, basins;
colormap = cmap, colorrange = (ids[1] - 0.5, ids[end]+0.5),
)
scatter_attractors!(ax, attractors)
fig
```

This aligns perfectly with the video we produced above.

- Datseris2022G. Datseris and A. Wagemakers,
*Effortless estimation of basins of attraction*, Chaos 32, 023104 (2022) - Stender2021Stender & Hoffmann, bSTAB: an open-source software for computing the basin stability of multi-stable dynamical systems
- Menck2013Menck, Heitzig, Marwan & Kurths. How basin stability complements the linear stability paradigm. Nature Physics, 9(2), 89–92
- Datseris2022G. Datseris and A. Wagemakers,
*Effortless estimation of basins of attraction*, Chaos 32, 023104 (2022) - Daza2016A. Daza, A. Wagemakers, B. Georgeot, D. Guéry-Odelin and M. A. F. Sanjuán, Basin entropy: a new tool to analyze uncertainty in dynamical systems, Sci. Rep., 6, 31416, 2016.
- Grebogi1983C. Grebogi, S. W. McDonald, E. Ott and J. A. Yorke, Final state sensitivity:
- Kaszás2019Kaszás, Feudel & Tél. Tipping phenomena in typical dynamical systems subjected to parameter drift. Scientific Reports, 9(1)