`Attractors.Attractors`

— Module**Attractors.jl**

Attractors.jl is a Julia package for

- Finding all attractors, and all types of attractors, of arbitrary dynamical systems. An extendable interface allows for new algorithms for finding attractors.
- Finding their basins of attraction or the state space fractions of the basins. This includes finding exit basins (divergence to infinity).
- Analyzing nonlocal stability of attractors (also called global stability or resilience).
- Performing
**global continuation**of attractors and their basins (or other measures of stability), over a parameter range. Global continuation is a new, cutting-edge type of continuation that offers several advantages over traditional local continuation (AUTO, MatCont, BifurcationKit.jl, etc.), see the comparison in our docs. - Finding the basin boundaries and edges states and analyzing their fractal properties.
- Tipping points related functionality for systems with known dynamic rule.
- And more!

It can be used as a standalone package, or as part of DynamicalSystems.jl.

To install it, run `import Pkg; Pkg.add("Attractors")`

.

All further information is provided in the documentation, which you can either find online or build locally by running the `docs/make.jl`

file.

*Previously, Attractors.jl was part of ChaosTools.jl*

`Attractors.AttractorMapper`

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

Subtypes of `AttractorMapper`

are structures that map initial conditions of `ds`

to attractors. The found attractors are stored inside the `mapper`

, and can be obtained by calling `attractors = extract_attractors(mapper)`

.

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:

`AttractorsViaProximity`

`AttractorsViaRecurrences`

`AttractorsViaFeaturizing`

with the`GroupViaHistogram`

configuration.

**For developers**

`AttractorMapper`

defines an extendable interface. A new type needs to implement `extract_attractors`

and `id = mapper(u0)`

. From these, everything else in the entire rest of the library just works! If it is not possible to implement `id = mapper(u0)`

, then instead extend `basins_fractions(mapper, ics)`

.

`Attractors.AttractorSeedContinueMatch`

— Type`AttractorSeedContinueMatch(mapper, matcher = MatchBySSSetDistance(); seeding)`

A global continuation method for `global_continuation`

. `mapper`

is any subtype of `AttractorMapper`

which implements `extract_attractors`

, i.e., it finds the actual attractors. `matcher`

is a configuration of how to match attractor IDs, and at the moment can only be an instance of `MatchBySSSetDistance`

.

**Description**

This is a general/composable global continuation method based on a 4-step process:

- Seed initial conditions from previously found attractors
- Propagate those forwards to "continue" previous attractors
- Estimate basin fractions and potentially find new attractors
- Match attractors

**Step 0 - Finding initial attractors**

At the first parameter slice of the global continuation process, attractors and their fractions are found using the given `mapper`

and `basins_fractions`

. See the `mapper`

documentation and `AttractorMapper`

for details on how this works. Then, from the second parameter onwards the continuation occurs.

**Step 1 - Seeding initial conditions**

Initial conditions can be seeded from previously found attractors. This is controlled by the `seeding`

keyword, which must be a function that given a `StateSpaceSet`

(an attractor), it returns an iterator of initial conditions. By default the first point of an attractor is provided as the only seed.

Seeding can be turned off by providing the dummy function `seeding = A -> []`

, i.e., it always returns an empty iterator and hence no seeds and we skip to step 2.

**Step 2 - Continuing the seeds**

The dynamical system referenced by the `mapper`

is now set to the new parameter value. The seeds are run through the `mapper`

to converge to attractors at the new parameter value. Seeding initial conditions close to previous attractors increases the probability that if an attractor continues to exist in the new parameter, it is found. Additionally, for some `mappers`

this seeding process improves the accuracy as well as performance of finding attractors, see e.g. discussion in Datseris2023.

This seeding works for any `mapper`

, regardless of if they can map individual initial conditions with the `mapper(u0)`

syntax! If this syntax isn't supported, steps 2 and 3 are done together.

**Step 3 - Estimate basins fractions**

After the special seeded initial conditions are mapped to attractors, attractor basin fractions are computed by sampling additional initial conditions using the provided `ics`

in `global_continuation`

. I.e., exactly as in `basins_fractions`

. Naturally, during this step new attractors may be found, besides those found using the "seeding from previous attractors".

**Step 4 - Matching**

Normally the ID an attractor gets assigned is somewhat a random integer. Therefore, to ensure a logical output of the global continuation process, attractors need to be "matched". This means: attractor and fractions must have their *IDs changed*, so that attractors that are "similar" to those at a previous parameter get assigned the same ID.

What is "similar enough" is controlled by the `matcher`

input. The default `matcher`

`MatchBySSSetDistance`

matches sets which have small distance in state space. The matching algorithm itself can be quite involved, so read the documentation of the `matcher`

for how matching works.

A note on matching: the `MatchBySSSetDistance`

can also be used after the continuation is completed, as it only requires as input the state space sets (attractors), without caring at which parameter each attractor exists at. If you don't like the final matching output, you may use a different instance of `MatchBySSSetDistance`

and call `match_sequentially!`

again on the output, without having to recompute the whole global continuation!

**Step 5 - Finish**

After matching the parameter is incremented. Steps 1-4 repeat until all parameter values are exhausted.

**Further note**

This global continuation method is a generalization of the "RAFM" continuation described in Datseris2023. This continuation method is still exported as `RecurrencesFindAndMatch`

.

`Attractors.AttractorsViaFeaturizing`

— Type```
AttractorsViaFeaturizing(
ds::DynamicalSystem, featurizer::Function,
grouping_config = GroupViaClustering(); kwargs...
)
```

Initialize a `mapper`

that maps initial conditions to attractors using a featurizing and grouping approach. This is a supercase of the featurizing and clustering approach that is utilized by bSTAB Stender2021 and MCBB Gelbrecht2020. See `AttractorMapper`

for how to use the `mapper`

. This `mapper`

also allows the syntax `mapper(u0)`

but only if the `grouping_config`

is *not* `GroupViaClustering`

.

`featurizer`

is a function `f(A, t)`

that takes as an input an integrated trajectory `A::StateSpaceSet`

and the corresponding time vector `t`

and returns a vector `v`

of features describing the trajectory. For better performance, it is strongly recommended that `v isa SVector{<:Real}`

.

`grouping_config`

is an instance of any subtype of `GroupingConfig`

and decides how features will be grouped into attractors, see below.

See also the intermediate functions `extract_features`

and `group_features`

, which can be utilized when wanting to work directly with features.

**Keyword arguments**

`T=100, Ttr=100, Δt=1`

: Propagated to`DynamicalSystems.trajectory`

for integrating an initial condition to yield`A, t`

.`threaded = true`

: Whether to run the generation of features over threads by integrating trajectories in parallel.

**Description**

The trajectory `X`

of an initial condition is transformed into features. Each feature is a number useful in *characterizing the attractor* the initial condition ends up at, and **distinguishing it from other attractors**. Example features are the mean or standard deviation of some the dimensions of the trajectory, the entropy of some of the dimensions, the fractal dimension of `X`

, or anything else you may fancy.

All feature vectors (each initial condition = 1 feature vector) are then grouped using one of the sevaral available grouping configurations. Each group is assumed to be a unique attractor, and hence each initial condition is labelled according to the group it is part of. 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, and the right way to group them, in contrast to `AttractorsViaRecurrences`

.

Attractors are stored and can be accessed with `extract_attractors`

, however it should be clear that this mapper never actually finds attractors. They way we store attractors is by picking the first initial condition that belongs to the corresponding "attractor group", and then recording its trajectory with the same arguments `T, Ttr, Δt`

. This is stored as the attractor, but of course there is no guarantee that this is actually an attractor.

`Attractors.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 `StateSpaceSet`

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 `StateSpaceSet`

s. If `length(attractors) == 1`

, then `ε`

becomes 1/10 of the diagonal of the box containing the attractor. If `length(attractors) == 1`

and the attractor is a single point, an error is thrown.

**Keywords**

`Ttr = 100`

: Transient time to first evolve the system for before checking for proximity.`Δt = 1`

: Step time given to`step!`

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

).`consecutive_lost_steps = 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`

).

`Attractors.AttractorsViaRecurrences`

— Type`AttractorsViaRecurrences(ds::DynamicalSystem, grid; kwargs...)`

Map initial conditions of `ds`

to attractors by identifying attractors on the fly based on recurrences in the state space, as outlined in Datseris2022. However, the Description section below for has a more accurate (and simpler) exposition to the algorithm than the paper.

`grid`

is instructions for partitioning the state space into finite-sized cells so that a finite state machine can operate on top of it. Possibilities are:

- A tuple of sorted
`AbstractRange`

s for a regular grid.

Example is `grid = (xg, yg)`

where `xg = yg = range(-5, 5; length = 100)`

for a two-dimensional system.

- A tuple of sorted
`AbstractVector`

s for an irregular grid, for example

`grid = (xg, yg)`

with `xg = range(0, 10.0^(1/2); length = 200).^2, yg = range(-5, 5; length = 100)`

.

- An instance of the special grid type

`SubdivisionBasedGrid`

, which can be created either manually or by using `subdivision_based_grid`

. This automatically analyzes and adapts grid discretization levels in accordance with state space flow speed in different regions.

The grid has to be the same dimensionality as the state space, use a `ProjectedDynamicalSystem`

if you want to search for attractors in a lower dimensional subspace.

**Keyword arguments**

`sparse = true`

: control the storage type of the state space grid. If true, uses a sparse array, whose memory usage is in general more efficient than a regular array obtained with`sparse=false`

. In practice, the sparse representation should always be preferred when searching for`basins_fractions`

. Only for very low dimensional systems and for computing the full`basins_of_attraction`

the non-sparse version should be used.

**Time evolution configuration**

`Ttr = 0`

: Skip a transient before the recurrence routine begins.`Δt`

: Approximate integration time step (second argument of the`step!`

function). The keyword`Dt`

can also be used instead if`Δ`

(`\Delta`

) is not accessible. It is`1`

for discrete time systems. For continuous systems, an automatic value is calculated using`automatic_Δt_basins`

. For very fine grids, this can become very small, much smaller than the typical integrator internal step size in case of adaptive integrators. In such cases, use`force_non_adaptive = true`

.`force_non_adaptive = false`

: Only used if the input dynamical system is`CoupledODEs`

. If`true`

the additional keywords`adaptive = false, dt = Δt`

are given as`diffeq`

to the`CoupledODEs`

. This means that adaptive integration is turned off and`Δt`

is used as the ODE integrator timestep. This is useful in (1) very fine grids, and (2) if some of the attractors are limit cycles. We have noticed that in this case the integrator timestep becomes commensurate with the limit cycle period, leading to incorrectly counting the limit cycle as more than one attractor.

**Finite state machine configuration**

`consecutive_recurrences = 100`

: Number of consecutive visits to previously visited unlabeled cells (i.e., recurrences) required before declaring we have converged to a new attractor. This number tunes the accuracy of converging to attractors and should generally be high (and even higher for chaotic systems).`attractor_locate_steps = 1000`

: Number of subsequent steps taken to locate accurately the new attractor after the convergence phase is over. Once`attractor_locate_steps`

steps have been taken, the new attractor has been identified with sufficient accuracy and iteration stops. This number can be very high without much impact to overall performance.`store_once_per_cell = true`

: Control if multiple points in state space that belong to the same cell are stored or not in the attractor, when a new attractor is found. If`true`

, each visited cell will only store a point once, which is desirable for fixed points and limit cycles. If`false`

then`attractor_locate_steps`

points are stored per attractor, leading to more densely stored attractors, which may be desirable for instance in chaotic attractors.`consecutive_attractor_steps = 2`

: Μaximum checks of consecutives hits of an existing attractor cell before declaring convergence to that existing attractor.`consecutive_basin_steps = 10`

: Number of consecutive visits of the same basin of attraction required before declaring convergence to an existing attractor. This is ignored if`sparse = true`

, as basins are not stored internally in that case.`consecutive_lost_steps = 20`

: Maximum check of iterations outside the defined grid before we declare the orbit lost outside and hence assign it label`-1`

.`horizon_limit = 1e6`

: If the norm of the integrator state reaches this limit we declare that the orbit diverged to infinity.`maximum_iterations = Int(1e6)`

: A safety counter that is always increasing for each initial condition. Once exceeded, the algorithm assigns`-1`

and throws a warning. This clause exists to stop the algorithm never halting for inappropriate grids. It may happen when a newly found attractor orbit intersects in the same cell of a previously found attractor (which leads to infinite resetting of all counters).

**Description**

An initial condition given to an instance of `AttractorsViaRecurrences`

is iterated based on the integrator corresponding to `ds`

. Enough recurrences in the state space (i.e., a trajectory visited a region it has visited before) 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 a cell in the given `grid`

. The grid cells store information: they are empty, visited, basins, or attractor cells. The state of the FSM is decided based on the cell type and the previous state of the FSM. Whenever the FSM recurs its state, its internal counter is increased, otherwise it is reset to 0. Once the internal counter reaches a threshold, the FSM terminates or changes its state. The possibilities for termination are the following:

- The trajectory hits
`consecutive_recurrences`

times in a row previously visited cells: it is considered that an attractor is found and is labelled with a new ID. Then, iteration continues for`attractor_locate_steps`

steps. Each cell visited in this period stores the "attractor" information. Then iteration terminates and the initial condition is numbered with the attractor's ID. - The trajectory hits an already identified attractor
`consecutive_attractor_steps`

consecutive times: the initial condition is numbered with the attractor's basin ID. - The trajectory hits a known basin
`consecutive_basin_steps`

times in a row: the initial condition belongs to that basin and is numbered accordingly. Notice that basins are stored and used only when`sparse = false`

otherwise this clause is ignored. - The trajectory spends
`consecutive_lost_steps`

steps outside the defined grid or the norm of the dynamical system state becomes > than`horizon_limit`

: the initial condition is labelled`-1`

. - If none of the above happens, the initial condition is labelled
`-1`

after`maximum_iterations`

steps.

There are some special internal optimizations and details that we do not describe here but can be found in comments in the source code. (E.g., a special timer exists for the "lost" state which does not interrupt the main timer of the FSM.)

A video illustrating how the algorithm works can be found in the online documentation, under the recurrences animation page.

`Attractors.EdgeTrackingResults`

— Type`EdgeTrackingResults(edge, track1, track2, time, bisect_idx)`

Data type that stores output of the `edgetracking`

algorithm.

**Fields**

`edge::StateSpaceSet`

: the pseudo-trajectory representing the tracked edge segment (given by the average in state space between`track1`

and`track2`

)`track1::StateSpaceSet`

: the pseudo-trajectory tracking the edge within basin 1`track2::StateSpaceSet`

: the pseudo-trajectory tracking the edge within basin 2`time::Vector`

: time points of the above`StateSpaceSet`

s`bisect_idx::Vector`

: indices of`time`

at which a re-bisection occurred`success::Bool`

: indicates whether the edge tracking has been successful or not

`Attractors.FeaturizeGroupAcrossParameter`

— Method```
FeaturizeGroupAcrossParameter <: GlobalContinuationAlgorithm
FeaturizeGroupAcrossParameter(mapper::AttractorsViaFeaturizing; kwargs...)
```

A method for `global_continuation`

. It uses the featurizing approach discussed in `AttractorsViaFeaturizing`

and hence requires an instance of that mapper as an input. When used in `global_continuation`

, features are extracted and then grouped across a parameter range. Said differently, all features of all initial conditions across all parameter values are put into the same "pool" and then grouped as dictated by the `group_config`

of the mapper. After the grouping is finished the feature label fractions are distributed to each parameter value they came from.

This continuation method is based on, but strongly generalizes, the approaches in the papers Gelbrecht2020 and Stender2021.

**Keyword arguments**

`info_extraction::Function`

a function that takes as an input a vector of feature-vectors (corresponding to a cluster) and returns a description of the cluster. By default, the centroid of the cluster is used. This is what the`attractors_cont`

contains in the return of`global_continuation`

.

`Attractors.GlobalContinuationAlgorithm`

— Type`GlobalContinuationAlgorithm`

Supertype of all algorithms used in `global_continuation`

. Each algorithm typically references an `AttractorMapper`

, as well as contains more information for how to continue/track/match attractors across a parameter range.

See `global_continuation`

for more.

`Attractors.GroupViaClustering`

— Type`GroupViaClustering(; kwargs...)`

Initialize a struct that contains instructions on how to group features in `AttractorsViaFeaturizing`

. `GroupViaClustering`

clusters features into groups using DBSCAN, similar to the original work by bSTAB Stender2021 and MCBB Gelbrecht2020. Several options on clustering are available, see keywords below.

The defaults are a significant improvement over existing literature, see Description.

**Keyword arguments**

`clust_distance_metric = Euclidean()`

: A metric to be used in the clustering. It can be any function`f(a, b)`

that returns the distance between real-valued vectors`a, b`

. All metrics from Distances.jl can be used here.`rescale_features = true`

: if true, rescale each dimension of the extracted features separately into the range`[0,1]`

. This typically leads to more accurate clustering.`min_neighbors = 10`

: minimum number of neighbors (i.e. of similar features) each feature needs to have, including counting its own self, in order to be considered in a cluster (fewer than this, it is labeled as an outlier,`-1`

).`use_mmap = false`

: whether to use an on-disk map for creating the distance matrix of the features. Useful when the features are so many where a matrix with side their length would not fit to memory.

**Keywords for optimal radius estimation**

`optimal_radius_method::Union{Real, String} = "silhouettes_optim"`

: if a real number, it is the radius used to cluster features. Otherwise, it determines the method used to automatically determine that radius. Possible values are:`"silhouettes"`

: Performs a linear (sequential) search for the radius that maximizes a statistic of the silhouette values of clusters (typically the mean). This can be chosen with`silhouette_statistic`

. The linear search may take some time to finish. To increase speed, the number of radii iterated through can be reduced by decreasing`num_attempts_radius`

(see its entry below).`"silhouettes_optim"`

: Same as`"silhouettes"`

but performs an optimized search via Optim.jl. It's faster than`"silhouettes"`

, with typically the same accuracy (the search here is not guaranteed to always find the global maximum, though it typically gets close).`"knee"`

: chooses the the radius according to the knee (a.k.a. elbow, highest-derivative method) and is quicker, though generally leading to much worse clustering. It requires that`min_neighbors`

> 1.

`num_attempts_radius = 100`

: number of radii that the`optimal_radius_method`

will try out in its iterative procedure. Higher values increase the accuracy of clustering, though not necessarily much, while always reducing speed.`silhouette_statistic::Function = mean`

: statistic (e.g. mean or minimum) of the silhouettes that is maximized in the "optimal" clustering. The original implementation in Stender2021 used the`minimum`

of the silhouettes, and typically performs less accurately than the`mean`

.`max_used_features = 0`

: if not`0`

, it should be an`Int`

denoting the max amount of features to be used when finding the optimal radius. Useful when clustering a very large number of features (e.g., high accuracy estimation of fractions of basins of attraction).

**Description**

The DBSCAN clustering algorithm is used to automatically identify clusters of similar features. Each feature vector is a point in a feature space. Each cluster then basically groups points that are closely packed together. Closely packed means that the points have at least `min_neighbors`

inside a ball of radius `optimal_radius`

centered on them. This method typically works well if the radius is chosen well, which is not necessarily an easy task. Currently, three methods are implemented to automatically estimate an "optimal" radius.

**Estimating the optimal radius**

The default method is the **silhouettes method**, which includes keywords `silhouette`

and `silhouette_optim`

. Both of them search for the radius that optimizes the clustering, meaning the one that maximizes a statistic `silhouette_statistic`

(e.g. mean value) of a quantifier for the quality of each cluster. This quantifier is the silhouette value of each identified cluster. A silhouette value measures how similar a point is to the cluster it currently belongs to, compared to the other clusters, and ranges from -1 (worst matching) to +1 (ideal matching). If only one cluster is found, the assigned silhouette is zero. So for each attempted radius in the search the clusters are computed, their silhouettes calculated, and the statistic of these silhouettes computed. The algorithm then finds the radius that leads to the maximum such statistic. For `optimal_radius_method = "silhouettes"`

, the search is done linearly, from a minimum to a maximum candidate radius for `optimal_radius_method = "silhouettes"`

; `optimal_radius_method = silhouettes_optim`

, it is done via an optimized search performed by Optim.jl which is typically faster and with similar accuracy. A third alternative is the`"elbow"`

method, which works by calculating the distance of each point to its k-nearest-neighbors (with `k=min_neighbors`

) and finding the distance corresponding to the highest derivative in the curve of the distances, sorted in ascending order. This distance is chosen as the optimal radius. It is described in Ester1996 and Schubert2017. It typically performs considerably worse than the `"silhouette"`

methods.

`Attractors.GroupViaHistogram`

— Type`GroupViaHistogram(binning::FixedRectangularBinning)`

Initialize a struct that contains instructions on how to group features in `AttractorsViaFeaturizing`

. `GroupViaHistogram`

performs a histogram in feature space. Then, all features that are in the same histogram bin get the same label. The `binning`

is an instance of `FixedRectangularBinning`

from ComplexityMeasures.jl. (the reason to not allow `RectangularBinning`

is because during continuation we need to ensure that bins remain identical).

`Attractors.GroupViaNearestFeature`

— Type`GroupViaNearestFeature(templates; kwargs...)`

Initialize a struct that contains instructions on how to group features in `AttractorsViaFeaturizing`

. `GroupViaNearestFeature`

accepts a `template`

, which is a vector of features. Then, generated features from initial conditions in `AttractorsViaFeaturizing`

are labelled according to the feature in `templates`

that is closest (the label is the index of the closest template).

`templates`

can be a vector or dictionary mapping keys to templates. Internally all templates are converted to `SVector`

for performance. Hence, it is strongly recommended that both `templates`

and the output of the `featurizer`

function in `AttractorsViaFeaturizing`

return `SVector`

types.

**Keyword arguments**

`metric = Euclidean()`

: metric to be used to quantify distances in the feature space.`max_distance = Inf`

: Maximum allowed distance between a feature and its nearest template for it to be assigned to that template. By default,`Inf`

guarantees that a feature is assigned to its nearest template regardless of the distance. Features that exceed`max_distance`

to their nearest template get labelled`-1`

.

`Attractors.GroupViaPairwiseComparison`

— Type`GroupViaPairwiseComparison(; threshold::Real, kwargs...)`

Initialize a struct that contains instructions on how to group features in `AttractorsViaFeaturizing`

. `GroupViaPairwiseComparison`

groups features and identifies clusters by considering the pairwise distance between features. It can be used as an alternative to the clustering method in `GroupViaClustering`

, having the advantage that it is simpler, typically faster and uses less memory.

**Keyword arguments**

`threshold`

(mandatory): A real number defining the maximum distance two features can be to be considered in the same cluster - above the threshold, features are different. This value simply needs to be large enough to differentiate clusters.`metric = Euclidean()`

: A function`metric(a, b)`

that returns the distance between two features`a`

and`b`

, outputs of`featurizer`

. Any`Metric`

from Distances.jl can be used here.`rescale_features = true`

: if true, rescale each dimension of the extracted features separately into the range`[0,1]`

. This typically leads to more accurate grouping.

**Description**

This algorithm assumes that the features are well-separated into distinct clouds, with the maximum radius of the cloud controlled by `threshold`

. Since the systems are deterministic, this is achievable with a good-enough `featurizer`

function, by removing transients, and running the trajectories for sufficiently long. It then considers that features belong to the same attractor when their pairwise distance, computed using `metric`

, is smaller than or equal to `threshold`

, and that they belong to different attractors when the distance is bigger. Attractors correspond to each grouping of similar features. In this way, the key parameter `threshold`

is basically the amount of variation permissible in the features belonging to the same attractor. If they are well-chosen, the value can be relatively small and does not need to be fine tuned.

The `threshold`

should achieve a balance: one one hand, it should be large enough to account for variations in the features from the same attractor - if it's not large enough, the algorithm will find duplicate attractors. On the other hand, it should be small enough to not group together features from distinct attractors. This requires some knowledge of how spread the features are. If it's too big, the algorithm will miss some attractors, as it groups 2+ distinct attractors together. Therefore, as a rule of thumb, one can repeat the procedure a few times, starting with a relatively large value and reducing it until no more attractors are found and no duplicates appear.

The method uses relatively little memory, as it only stores vectors whose size is on order of the number of attractors of the system.

`Attractors.GroupingConfig`

— TypeGroupingConfig

Supertype for configuration structs on how to group features together. Used in several occasions such as `AttractorsViaFeaturizing`

or `aggregate_attractor_fractions`

.

Currently available grouping configurations are:

**For developers**

`GroupingConfig`

defines an extendable interface. The only thing necessary for a new grouping configuration is to:

- Make a new type and subtype
`GroupingConfig`

. - If the grouping allows for mapping individual features to group index, then instead extend the
**internal function**`feature_to_group(feature, config)`

. This will also allow doing`id = mapper(u0)`

with`AttractorsViaFeaturizing`

. - Else, extend the function
`group_features(features, config)`

. You could still extend`group_features`

even if (2.) is satisfied, if there are any performance benefits. - Include the new grouping file in the
`grouping/all_grouping_configs.jl`

and list it in this documentation string.

`Attractors.IDMatcher`

— Type`IDMatcher`

Supertype of all "matchers" that match can IDs labelling attractors. Currently available matchers:

Matchers implement an extendable interface based on the function `matching_map`

. This function is used by the higher level function `match_sequentially!`

, which can be called after any call to a global continuation to match attractors differently, if the matching used originally during the continuation was not the best.

`Attractors.MFSBlackBoxOptim`

— Type`MFSBlackBoxOptim(; kwargs...)`

The black box derivative-free optimization algorithm used in `minimal_fatal_shock`

.

**Keyword arguments**

`guess = nothing`

: a initial guess for the minimal fatal shock given to the optimization algorithm. If not`nothing`

,`random_algo`

below is ignored.`max_steps = 10000`

: maximum number of steps for the optimization algorithm.`penalty = 1000.0`

: penalty value for the objective function for perturbations that do not lead to a different basin of attraction. This value is added to the norm of the perturbation and its value should be much larger than the typical sizes of the basins of attraction.`print_info`

: boolean value, if true, the optimization algorithm will print information on the evaluation steps of objective function,`default = false`

.`random_algo = MFSBruteForce(100, 100, 0.99)`

: an instance of`MFSBruteForce`

that can be used to provide an initial guess.`bbkwargs = NamedTuple()`

: additional keyword arguments propagated to`BlackBoxOptim.bboptimize`

for selecting solver, accuracy, and more.

**Description**

The algorithm uses BlackBoxOptim.jl and a penalized objective function to minimize. y function used as a constraint function. So, if we hit another basin during the search we encourage the algorithm otherwise we punish it with some penalty. The function to minimize is (besides some details):

```
function mfs_objective(perturbation, u0, mapper, penalty)
dist = norm(perturbation)
if mapper(u0 + perturbation) == mapper(u0)
# penalize if we stay in the same basin:
return dist + penalty
else
return dist
end
end
```

Using an initial guess can be beneficial to both performance and accuracy, which is why the output of a crude `MFSBruteForce`

is used to provide a guess. This can be disabled by either passing a `guess`

vector explicitly or by giving `nothing`

as `random_algo`

.

`Attractors.MFSBruteForce`

— Type`MFSBruteForce(; kwargs...)`

The brute force randomized search algorithm used in `minimal_fatal_shock`

.

It consists of two steps: random initialization and sphere radius reduction. On the first step, the algorithm generates random perturbations within the search area and records the perturbation that leads to a different basin but with the smallest magnitude. With this obtained perturbation it proceeds to the second step. On the second step, the algorithm generates random perturbations on the surface of the hypersphere with radius equal to the norm of the perturbation found in the first step. It reduces the radius of the hypersphere and continues searching for the better result with a smaller radius. Each time a better result is found, the radius is reduced further.

The algorithm records the perturbation with smallest radius that leads to a different basin.

**Keyword arguments**

`initial_iterations = 10000`

: number of random perturbations to try in the first step of the algorithm.`sphere_iterations = 10000`

: number of steps while initializing random points on hypersphere and decreasing its radius.`sphere_decrease_factor = 0.999`

factor by which the radius of the hypersphere is decreased (at each step the radius is multiplied by this number). Number closer to 1 means more refined accuracy

`Attractors.MatchByBasinEnclosure`

— Type`MatchByBasinEnclosure(; kw...) <: IDMatcher`

A matcher that matches attractors by whether they are enclosed in the basin of a new attractor or not.

**Keyword arguments**

`ε = nothing`

: distance threshold given to`AttractorsViaProximity`

. If`nothing`

, it is estimated as a quarter of the minimum distance of centroids (in contrast to the default more accurate estimation in`AttractorsViaProximity`

).`Δt = 1, consecutive_lost_steps = 1000`

: also given to`AttractorsViaProximity`

. We have not yet decided what should happen to attractors that did not converge to one of the current attractors within this number of steps. At the moment they get assigned the next available free ID but this may change in future releases.`distance = Centroid()`

: metric to estimate distances between state space sets in case there are co-flowing attractors, see below.`seeding = A -> A[end]`

: how to select a point from the attractor to see if it is enclosed in the basin of a new attractor.

**Description**

An attractor `A₋`

is a set in a state space that occupies a particular region (or, a single point, if it is a fixed point). This region is always within the basin of attraction of said attractor. When the parameter of the dynamical system is incremented, the attractors `A₊`

in the new parameter have basins that may have changed in shape and size.

The new attractor `A₊`

is "matched" (i.e., has its ID changed) to the old attractor `A₋`

attractor if `A₋`

is located inside the basin of attraction of `A₊`

. To see if `A₋`

is in the basin of `A₊`

, we first pick a point from `A₋`

using the `seeding`

keyword argument. By default this is the last point on the attractor, but it could be anything else, including the centroid of the attractor (`mean(A)`

). This point is given as an initial condition to an `AttractorsViaProximity`

mapper that maps initial conditions to the `₊`

attractors when the trajectories from the initial conditions are `ε`

-close to the `₊`

attractors.

There can be the situation where multiple `₋`

attractors converge to the same `₊`

attractor, which we call "coflowing attractors". In this scenario matching is prioritized for the `₋`

attractor that is closest to the `₊`

in terms of state space set distance, which is estimated with the `distance`

keyword, which can be anything `MatchBySSSetDistance`

accepts. The closest `₊`

attractor gets the ID of the `₋`

closest attractor that converge to it.

Basin enclosure is a concept similar to "basin instability" in Ritchie2023.

`Attractors.MatchByBasinOverlap`

— Type`MatchByBasinOverlap(threshold = Inf)`

A matcher that matches IDs given full basins of attraction.

**Description**

This matcher cannot be used in with the generic global continuation method of `AttractorSeedContinueMatch`

. This matcher matches IDs of attractors whose basins of attraction before and after `b₋, b₊`

have the most overlap (in pixels). This overlap is normalized in 0-1 (with 1 meaning 100% of a basin in `b₋`

is overlaping with some other basin in `b₊`

). Therefore, the values this matcher compares are *full basins of attraction*, not attractors themselves (hence why it can't be given to `AttractorSeedContinueMatch`

). Rather, you may use this matcher with `matching_map`

.

The `threshold`

can dissallow matching between basins that do not have enough overlap. Basins whose overlap is less than `1/threshold`

are guaranteed to get assined different IDs. For example: for `threshold = 2`

basins that have ≤ 50% overlap get different IDs guaranteed. By default, there is no threshold.

The information of the basins of attraction is typically an `Array`

, i.e., the direct output of `basins_of_attraction`

. For convenience, as well as backwards compatibility, when using `matching_map`

with this mapper you may provide two `Array`

s `b₊, b₋`

representing basins of attraction after and before, and the conversion to dictionaries will happen internally as it is supposed to. To replace the `IDs`

in `b₊`

given the replacement map just call `replace!(b₊, rmap...)`

, or use the in-place version `matching_map!`

directly.

A lower-level input for this matcher in `matching_map`

can be dictionaries mapping IDs to vectors of cartesian indices, where the indices mean which parts of the state space belong to which ID

`Attractors.MatchBySSSetDistance`

— Type`MatchBySSSetDistance(; distance = Centroid(), threshold = Inf, use_vanished = false)`

A matcher type that matches IDs by the distance of their corresponding state space sets.

**Keyword arguments**

`distance = Centroid()`

: distance to match by, given to`setsofsets_distances`

.`threshold = Inf`

: sets with distance larger than the`threshold`

are guaranteed to not be mapped to each other.`use_vanished = !isinf(threshold)`

: value of the keyword`use_vanished`

when used in`match_sequentially!`

.

**Description**

In this matcher the values compared are `StateSpaceSet`

s which in most cases represent attractors in the state space, but may also represent any other set such as a group of features.

Here is how this matcher works: (recall in this conversation that sets/attractors are stored in dictionaries, mapping keys/IDs to the sets, and we want to match keys in the "new" dictionary (`a₊`

) to those in the "old" dictionary (`a₋`

)).

The distance between all possible pairs of sets between the "old" sets and "new" sets is computed as a formal distance between sets. This is controlled by the `distance`

option, itself given to the lower-level `setsofsets_distances`

function, so `distance`

can be whatever that function accepts. That is, one of `Centroid`

, `Hausdorff`

, `StrictlyMinimumDistance`

, or any arbitrary user-provided function `f`

that given two sets `f(A, B)`

it returns a positive number (their distance).

Sets (in particular, their corresponding IDs) are then matched according to this distance. First, all possible ID pairs (old, new) are sorted according to the distance of their corresponding sets. The pair with smallest distance is matched. IDs in matched pairs are removed from the matching pool to ensure a unique mapping. Then, the next pair with least remaining distance is matched, and the process repeats until all pairs are exhausted.

Additionally, you can provide a `threshold`

value. If the distance between two sets is larger than this `threshold`

, then it is guaranteed that the two sets will get assigned different ID in the replacement map, and hence, the set in `a₊`

gets the next available integer as its ID.

`Attractors.SubdivisionBasedGrid`

— Type`SubdivisionBasedGrid(grid::NTuple{D, <:AbstractRange}, lvl_array::Array{Int, D})`

Given a coarse `grid`

tesselating the state space, construct a `SubdivisionBasedGrid`

based on the given level array `lvl_array`

that should have the same dimension as `grid`

. The level array has non-negative integer values, with 0 meaning that the corresponding cell of the coarse `grid`

should not be subdivided any further. Value `n > 0`

means that the corresponding cell will be subdivided in total `2^n`

times (along each dimension), resulting in finer cells within the original coarse cell.

`Attractors.RAFM`

— FunctionAlias for `RecurrencesFindAndMatch`

.

`Attractors.RecurrencesFindAndMatch`

— Method```
RecurrencesFindAndMatch <: GlobalContinuationAlgorithm
RecurrencesFindAndMatch(mapper::AttractorsViaRecurrences; kwargs...)
```

A method for `global_continuation`

as in Datseris2023 that is based on the recurrences algorithm for finding attractors (`AttractorsViaRecurrences`

) and then matching them according to their state space distance.

**Keyword arguments**

`distance = Centroid(), threshold = Inf`

: passed to`MatchBySSSetDistance`

.`seeds_from_attractor`

: A function that takes as an input an attractor and returns an iterator of initial conditions to be seeded from the attractor for the next parameter slice. By default, we sample only the first stored point on the attractor.

**Description**

`RecurrencesFindAndMatch`

is a wrapper type. It is has been generalized by `AttractorSeedContinueMatch`

. It is still exported for backwards compatibility and to have a clear reference to the original algorithm developed in Datseris2023.

The source code of `RecurrencesFindAndMatch`

is trival: it takes the given mapper, it initializes a `MatchBySSSetDistance`

, and along with `seeds_from_attractor`

it makes the `AttractorSeedContinueMatch`

instance. This is the process described in Datseris2023, whereby attractors are found using the recurrences algorithm `AttractorsViaRecurrences`

and they are then matched by their distance in state space `MatchBySSSetDistance`

.

`Attractors._rescale_to_01`

— MethodDo "min-max" rescaling of vector of feature vectors so that its values span `[0,1]`

.

`Attractors.additive_dict_merge!`

— Method`additive_dict_merge!(d1::Dict, d2::Dict)`

Merge keys and values of `d2`

into `d1`

additively: the values of the same keys are added together in `d1`

and new keys are given to `d1`

as-is.

`Attractors.aggregate_attractor_fractions`

— Function```
aggregate_attractor_fractions(
fractions_cont, attractors_cont, featurizer, group_config [, info_extraction]
)
```

Aggregate the already-estimated curves of fractions of basins of attraction of similar attractors using the same pipeline used by `GroupingConfig`

. The most typical application of this function is to transform the output of a `global_continuation`

with `RecurrencesFindAndMatch`

so that similar attractors, even across parameter space, are grouped into one "attractor". Thus, the fractions of their basins are aggregated.

You could also use this function to aggregate attractors and their fractions even in a single parameter configuration, i.e., using the output of `basins_fractions`

.

This function is useful in cases where you want the accuracy and performance of `AttractorsViaRecurrences`

, but you also want the convenience of "grouping" similar attractrors like in `AttractorsViaFeaturizing`

for presentation or analysis purposes. For example, a high dimensional model of competition dynamics across multispecies may have extreme multistability. After finding this multistability however, one may care about aggregating all attractors into two groups: where a given species is extinct or not. This is the example highlighted in our documentation, in Extinction of a species in a multistable competition model.

**Input**

`fractions_cont`

: a vector of dictionaries mapping labels to basin fractions.`attractors_cont`

: a vector of dictionaries mapping labels to attractors. 1st and 2nd argument are exactly like the return values of`global_continuation`

with`RecurrencesFindAndMatch`

(or, they can be the return of`basins_fractions`

).`featurizer`

: a 1-argument function to map an attractor into an appropriate feature to be grouped later. Features expected by`GroupingConfig`

are`SVector`

.`group_config`

: a subtype of`GroupingConfig`

.`info_extraction`

: a function accepting a vector of features and returning a description of the features. I.e., exactly as in`FeaturizeGroupAcrossParameter`

. The 5th argument is optional and defaults to the centroid of the features.

**Return**

`aggregated_fractions`

: same as`fractions_cont`

but now contains the fractions of the aggregated attractors.`aggregated_info`

: dictionary mapping the new labels of`aggregated_fractions`

to the extracted information using`info_extraction`

.

**Clustering attractors directly**

*(this is rather advanced)*

You may also use the DBSCAN clustering approach here to group attractors based on their state space distance (the `set_distance`

) by making a distance matrix as expected by the DBSCAN implementation. For this, use `identity`

as `featurizer`

, and choose `GroupViaClustering`

as the `group_config`

with `clust_distance_metric = set_distance`

and provide a numerical value for `optimal_radius_method`

when initializing the `GroupViaClustering`

, and also, for the `info_extraction`

argument, you now need to provide a function that expects a *vector of StateSpaceSets* and outputs a descriptor. E.g.,

`info_extraction = vector -> mean(mean(x) for x in vector)`

.`Attractors.animate_attractors_continuation`

— Function```
animate_attractors_continuation(
ds::DynamicalSystem, attractors_cont, fractions_cont, prange, pidx;
kwargs...
)
```

Animate how the found system attractors and their corresponding basin fractions change as the system parameter is increased. This function combines the input and output of the `global_continuation`

function into a video output.

The input dynamical system `ds`

is used to evolve initial conditions sampled from the found attractors, so that the attractors are better visualized. `attractors_cont, fractions_cont`

are the output of `global_continuation`

while `ds, prange, pidx`

are the input to `global_continuation`

.

**Keyword arguments**

`savename = "attracont.mp4"`

: name of video output file.`framerate = 4`

: framerate of video output.`Δt, T`

: propagated to`trajectory`

for evolving an initial condition sampled from an attractor.- Also all common plotting keywords.
`figure, axis, fracaxis`

: named tuples propagated as keyword arguments to the creation of the`Figure`

, the`Axis`

and the "bar-like" axis containing the fractions.

`Attractors.automatic_Δt_basins`

— Method`automatic_Δt_basins(ds::DynamicalSystem, 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 within the bounding box of the grid. The average `f`

is then compared with the average 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 if the trajectory skips a few cells. Also, `Δt`

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

`Attractors.basin_entropy`

— Function`basin_entropy(basins::Array, ε = 20) -> Sb, Sbb`

Compute the basin entropy Daza2016 `Sb`

and basin boundary entropy `Sbb`

of the given `basins`

of attraction by considering `ε`

boxes along each dimension.

**Description**

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

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 arbitrarily 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 `ε`

.

`Attractors.basins_fractal_dimension`

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

Estimate the fractal dimension `d`

of the boundary between basins of attraction using a box-counting algorithm for the boxes that contain at least two different basin IDs.

**Keyword arguments**

`range_ε = 2:maximum(size(basins))÷20`

is the range of sizes of the box to test (in pixels).

**Description**

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 output `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.

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.

`Attractors.basins_fractal_test`

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

Perform an automated test to decide if the boundary of the basins has fractal structures based on the method of Puy et al. Puy2021. Return `test_res`

(`:fractal`

or `:smooth`

) and the mean basin boundary entropy.

**Keyword arguments**

`ε = 20`

: size of the box to compute the basin boundary entropy.`Ntotal = 1000`

: number of balls to test in the boundary for the computation of`Sbb`

**Description**

The test "looks" at the basins with a magnifier of size `ε`

at random. If what we see in the magnifier looks like a smooth boundary (onn 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 boxes 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.

`Attractors.basins_fractions`

— Function`basins_fractions(basins::AbstractArray [,ids]) → 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. Optionally you may give a vector of `ids`

to calculate the fractions of only the chosen ids (by default `ids = unique(basins)`

).

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. For this, see `global_continuation`

.

`Attractors.basins_fractions`

— Method```
basins_fractions(
mapper::AttractorMapper,
ics::Union{StateSpaceSet, Function};
kwargs...
)
```

Approximate the state space fractions `fs`

of the basins of attraction of a dynamical system by mapping initial conditions to attractors using `mapper`

(which contains a reference to a `DynamicalSystem`

). 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
`StateSpaceSet`

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.

**Return**

The function will always return `fractions`

, which is a dictionary whose keys are the labels given to each attractor (always integers enumerating the different attractors), and whose values are the respective basins fractions. The label `-1`

is given to any initial condition where `mapper`

could not match to an attractor (this depends on the `mapper`

type).

If `ics`

is a `StateSpaceSet`

the function will also return `labels`

, which is a *vector*, of equal length to `ics`

, that contains the label each initial condition was mapped to.

See `AttractorMapper`

for all possible `mapper`

types, and use `extract_attractors`

(after calling `basins_fractions`

) to extract the stored attractors from the `mapper`

. See also `convergence_and_basins_fractions`

.

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

`Attractors.basins_of_attraction`

— Method`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 `DynamicalSystem`

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 `ProjectedDynamicalSystem`

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 assigns them to a full array corresponding to the state space partitioning indicated by `grid`

.

See also `convergence_and_basins_of_attraction`

.

`Attractors.basins_of_attraction`

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

`Attractors.bisect_to_edge`

— Method`bisect_to_edge(pds::ParallelDynamicalSystem, mapper::AttractorMapper; kwargs...) -> u1, u2`

Finds the basin boundary between two states `u1, u2 = current_states(pds)`

by bisecting along a straight line in phase space. The states `u1`

and `u2`

must belong to different basins.

Returns a triple `u1, u2, success`

, where `u1, u2`

are two new states located on either side of the basin boundary that lie less than `bisect_thresh`

(Euclidean distance in state space) apart from each other, and `success`

is a Bool indicating whether the bisection was successful (it may fail if the `mapper`

maps both states to the same basin of attraction, in which case a warning is raised).

**Keyword arguments**

`bisect_thresh = 1e-7`

: The maximum (Euclidean) distance between the two returned states.

**Description**

`pds`

is a `ParallelDynamicalSystem`

with two states. The `mapper`

must be an `AttractorMapper`

of subtype `AttractorsViaProximity`

or `AttractorsViaRecurrences`

.

If the straight line between `u1`

and `u2`

intersects the basin boundary multiple times, the method will find one of these intersection points. If more than two attractors exist, one of the two returned states may belong to a different basin than the initial conditions `u1`

and `u2`

. A warning is raised if the bisection involves a third basin.

`Attractors.cluster_assignment`

— MethodUtil function for `classify_features`

. Returns the assignment vector, in which the i-th component is the cluster index of the i-th feature.

`Attractors.convergence_and_basins_fractions`

— Method`convergence_and_basins_fractions(mapper::AttractorMapper, ics::StateSpaceSet)`

An extension of `basins_fractions`

. Return `fs, labels, convergence`

. The first two are as in `basins_fractions`

, and `convergence`

is a vector containing the time each initial condition took to converge to its attractor. Only usable with mappers that support `id = mapper(u0)`

.

See also `convergence_time`

.

**Keyword arguments**

`show_progress = true`

: show progress bar.

`Attractors.convergence_and_basins_of_attraction`

— Method`convergence_and_basins_of_attraction(mapper::AttractorMapper, grid)`

An extension of `basins_of_attraction`

. Return `basins, attractors, convergence`

, with `basins, attractors`

as in `basins_of_attraction`

, and `convergence`

being an array with same shape as `basins`

. It contains the time each initial condition took to converge to its attractor. It is useful to give to `shaded_basins_heatmap`

.

See also `convergence_time`

.

**Keyword arguments**

`show_progress = true`

: show progress bar.

`Attractors.convergence_time`

— Function`convergence_time(mapper::AttractorMapper) → t`

Return the approximate time the `mapper`

took to converge to an attractor. This function should be called just right after `mapper(u0)`

was called with `u0`

the initial condition of interest. Hence it is only valid with `AttractorMapper`

subtypes that support this syntax.

Obtaining the convergence time is computationally free, so that `convergence_and_basins_fractions`

can always be used instead of `basins_fractions`

.

`Attractors.crude_initial_radius`

— MethodThis function generates a random perturbation of the initial point `u0`

within specified "search*area" and checks if it is in the same basin of attraction. It does so by generating a random vector of length dim and then adding it to u0. If the perturbation is not in the same basin of attraction, it calculates the norm of the perturbation and compares it to the best perturbation found so far. If the norm is smaller, it updates the best perturbation found so far. It repeats this process total*iterations times and returns the best perturbation found.

`Attractors.edgetracking`

— Method`edgetracking(ds::DynamicalSystem, attractors::Dict; kwargs...)`

Track along a basin boundary in a dynamical system `ds`

with two or more attractors in order to find an *edge state*. Results are returned in the form of `EdgeTrackingResults`

, which contains the pseudo-trajectory `edge`

representing the track on the basin boundary, along with additional output (see below).

The system's `attractors`

are specified as a `Dict`

of `StateSpaceSet`

s, as in `AttractorsViaProximity`

or the output of `extract_attractors`

. By default, the algorithm is initialized from the first and second attractor in `attractors`

. Alternatively, the initial states can be set via keyword arguments `u1`

, `u2`

(see below). Note that the two initial states must belong to different basins of attraction.

**Keyword arguments**

`bisect_thresh = 1e-7`

: distance threshold for bisection`diverge_thresh = 1e-6`

: distance threshold for parallel integration`u1`

: first initial state (defaults to first point in first entry of`attractors`

)`u2`

: second initial state (defaults to first point in second entry of`attractors`

)`maxiter = 100`

: maximum number of iterations before the algorithm stops`abstol = 0.0`

: distance threshold for convergence of the updated edge state`T_transient = 0.0`

: transient time before the algorithm starts saving the edge track`tmax = Inf`

: maximum integration time of parallel trajectories until re-bisection`Δt = 0.01`

: time step passed to`step!`

when evolving the two trajectories`ϵ_mapper = nothing`

:`ϵ`

parameter in`AttractorsViaProximity`

`show_progress = true`

: if true, shows progress bar and information while running`verbose = true`

: if false, silences print output and warnings while running`kwargs...`

: additional keyword arguments to be passed to`AttractorsViaProximity`

**Description**

The edge tracking algorithm is a numerical method to find an *edge state* or (possibly chaotic) saddle on the boundary between two basins of attraction. Introduced by Battelino1988 and further described by Skufca2006, the algorithm has been applied to, e.g., the laminar-turbulent boundary in plane Couette flow Schneider2008, Wada basins Wagemakers2020, as well as Melancholia states in conceptual Mehling2023 and intermediate-complexity Lucarini2017 climate models. Relying only on forward integration of the system, it works even in high-dimensional systems with complicated fractal basin boundary structures.

The algorithm consists of two main steps: bisection and tracking. First, it iteratively bisects along a straight line in state space between the intial states `u1`

and `u2`

to find the separating basin boundary. The bisection stops when the two updated states are less than `bisect_thresh`

(Euclidean distance in state space) apart from each other. Next, a `ParallelDynamicalSystem`

is initialized from these two updated states and integrated forward until the two trajectories diverge from each other by more than `diverge_thresh`

(Euclidean distance). The two final states of the parallel integration are then used as new states `u1`

and `u2`

for a new bisection, and so on, until a stopping criterion is fulfilled.

Two stopping criteria are implemented via the keyword arguments `maxiter`

and `abstol`

. Either the algorithm stops when the number of iterations reaches `maxiter`

, or when the state space position of the updated edge point changes by less than `abstol`

(in Euclidean distance) compared to the previous iteration. Convergence below `abstol`

happens after sufficient iterations if the edge state is a saddle point. However, the edge state may also be an unstable limit cycle or a chaotic saddle. In these cases, the algorithm will never actually converge to a point but (after a transient period) continue populating the set constituting the edge state by tracking along it.

A central idea behind this algorithm is that basin boundaries are typically the stable manifolds of unstable sets, namely edge states or saddles. The flow along the basin boundary will thus lead to these sets, and the iterative bisection neutralizes the unstable direction of the flow away from the basin boundary. If the system possesses multiple edge states, the algorithm will find one of them depending on where the initial bisection locates the boundary.

**Output**

Returns a data type `EdgeTrackingResults`

containing the results.

Sometimes, the AttractorMapper used in the algorithm may erroneously identify both states `u1`

and `u2`

with the same basin of attraction due to being very close to the basin boundary. If this happens, a warning is raised and `EdgeTrackingResults.success = false`

.

`Attractors.edgetracking`

— Method`edgetracking(pds::ParallelDynamicalSystem, mapper::AttractorMapper; kwargs...)`

Low-level function for running the edge tracking algorithm, see `edgetracking`

for a description, keyword arguments and output type.

`pds`

is a `ParallelDynamicalSystem`

with two states. The `mapper`

must be an `AttractorMapper`

of subtype `AttractorsViaProximity`

or `AttractorsViaRecurrences`

.

`Attractors.extract_attractors`

— Method`extract_attractors(mapper::AttractorsMapper) → attractors`

Return a dictionary mapping label IDs to attractors found by the `mapper`

. This function should be called after calling `basins_fractions`

with the given `mapper`

so that the attractors have actually been found first.

For `AttractorsViaFeaturizing`

, the attractors are only stored if the mapper was called with pre-defined initial conditions rather than a sampler (function returning initial conditions).

`Attractors.extract_features`

— Method`extract_features(mapper, ics; N = 1000, show_progress = true)`

Return a vector of the features of each initial condition in `ics`

(as in `basins_fractions`

), using the configuration of `mapper::AttractorsViaFeaturizing`

. Keyword `N`

is ignored if `ics isa StateSpaceSet`

.

`Attractors.feature_to_group`

— Method`feature_to_group(feature::AbstractVector, group_config::GroupingConfig) → group_label`

Map the given feature vector to its group label (integer). This is an internal function. It is strongly recommended that `feature isa SVector`

.

`Attractors.finite_state_machine!`

— MethodMain procedure which performs one step of the finite state machine. Directly implements the algorithm of Datseris & Wagemakers 2021, see the flowchart (Figure 2). I.e., it obtains current state, updates state with respect to the cell code the trajectory it is currently on, performs the check of counters, and sets the f.s.m. to the new state, or increments counter if same state.

The basins and attractors are coded in the array with odd numbers for the basins and even numbers for the attractors. The attractor `2n`

has the corresponding basin `2n+1`

. This codification is changed when the basins and attractors are returned to the user. Diverging trajectories and the trajectories staying outside the grid are coded with -1.

The label `1`

(initial value) outlined in the paper is `0`

here instead. The function returns `0`

unless the FSM has terminated its operation.

`Attractors.global_continuation`

— Function`global_continuation(gca::GlobalContinuationAlgorithm, prange, pidx, ics; kwargs...)`

Find and continue attractors (or representations of attractors) and the fractions of their basins of attraction across a parameter range. `global_continuation`

is the central function of the framework for global stability analysis illustrated in Datseris2023.

The global continuation algorithm typically references an `AttractorMapper`

which is used to find the attractors and basins of a dynamical system. Additional arguments that control how to continue/track/match attractors across a parameter range are given when creating `gca`

.

The basin fractions and the attractors (or some representation of them) are continued across the parameter range `prange`

, for the parameter of the system with index `pidx`

(any index valid in `DynamicalSystems.set_parameter!`

can be used).

`ics`

are the initial conditions to use when globally sampling the state space. Like in `basins_fractions`

it can be either a set vector of initial conditions, or a 0-argument function that generates random initial conditions.

Possible subtypes of `GlobalContinuationAlgorithm`

are:

**Return**

`fractions_cont::Vector{Dict{Int, Float64}}`

. The fractions of basins of attraction.`fractions_cont[i]`

is a dictionary mapping attractor IDs to their basin fraction at the`i`

-th parameter.`attractors_cont::Vector{Dict{Int, <:Any}}`

. Information about the attractors.`attractors_cont[i]`

is a dictionary mapping attractor ID to information about the attractor at the`i`

-th parameter. The type of information stored depends on the chosen global continuation type, but typically it is the attractors themselves as`StateSpaceSet`

s.

**Keyword arguments**

`show_progress = true`

: display a progress bar of the computation.`samples_per_parameter = 100`

: amount of initial conditions sampled at each parameter from`ics`

if`ics`

is a function instead of set initial conditions.

`Attractors.group_features`

— Method`group_features(features, group_config::GroupingConfig) → labels`

Group the given iterable of "features" (anything that can be grouped, typically vectors of real numbers) according to the configuration and return the labels (vector of equal length as `features`

). See `GroupingConfig`

for possible grouping configuration configurations.

`Attractors.haussdorff_distance`

— Method`haussdorff_distance(M1, M2) -> hd`

Compute the Hausdorff distance between two binary matrices of the same size. First a distance matrix is computed using `distance_matrix`

for each matrix M1 and M2. The entries of this matrix are the distance in L1 metric to the closest 0 pixel in the initial matrix. The distance being 0 if it belongs to the basin with id 0.

`Attractors.heatmap_basins_attractors`

— Function`heatmap_basins_attractors(grid, basins, attractors; kwargs...)`

Plot a heatmap of found (2-dimensional) `basins`

of attraction and corresponding `attractors`

, i.e., the output of `basins_of_attraction`

.

**Keyword arguments**

- All the common plotting keywords.

`Attractors.match_sequentially!`

— Method`match_sequentially!(dicts::Vector{Dict{Int, Any}}, matcher::IDMatcher; kw...)`

Match the `dicts`

, a vector of dictionaries mapping IDs (integers) to values, according to the given `matcher`

by sequentially applying the `matching_map`

function to all elements of `dicts`

besides the first one.

In the context of Attractors.jl `dicts`

are typically dictionaries mapping IDs to attractors (`StateSpaceSet`

s), however the function is generic and would work for any values that `matcher`

works with.

Return `rmaps`

, which is a vector of dictionaries. `rmaps[i]`

contains the `matching_map`

for `attractors[i+1]`

, i.e., the pairs of `old => new`

IDs.

**Keyword arguments**

`prange = eachindex(attractors)`

: the range of parameters from which to extract the`p, pprev`

values given to`matching_map`

.`pidx, ds`

: both propagated to`matching_map`

and are`nothing`

by default.`retract_keys::Bool = true`

: If`true`

at the end the function will "retract" keys (i.e., make the integers smaller integers) so that all unique IDs are the 1-incremented positive integers. E.g., if the IDs where 1, 6, 8, they will become 1, 2, 3. The special ID -1 is unaffected by this.`use_vanished = false`

: If`use_vanised = true`

, then IDs (and their corresponding sets) that existed before but have vanished are kept in "memory" when it comes to matching: the current dictionary values (the attractor sets) are compared to the latest instance of all values that have ever existed, each with a unique ID, and get matched to their closest ones. The value of this keyword is obtained from the`matcher`

.

`Attractors.match_sequentially!`

— Method`match_sequentially!(continuation_quantity::Vector{Dict}, rmaps::Vector{Dict})`

Do the same as in `match_sequentially!`

above, now given the vector of matching maps, and for any arbitrary quantity that has been tracked in the global*continuation. `continuation*quantity`can for example be`

fractions*cont from [global*continuation`](@ref).

`Attractors.matching_map!`

— Method`matching_map!(a₊, a₋, matcher; kw...) → rmap`

Convenience function that first calls `matching_map`

and then replaces the IDs in `a₊`

with this `rmap`

.

`Attractors.matching_map`

— Method`matching_map(b₊::AbstractArray, b₋::AbstractArray, matcher::MatchByBasinOverlap)`

Special case of `matching_map`

where instead of having as input dictionaries mapping IDs to values, we have `Array`

s which represent basins of attraction and whose elements are the IDs.

See `MatchByBasinOverlap`

for how matching works.

`Attractors.matching_map`

— Method```
matching_map(
a₊::Dict, a₋::Dict, matcher;
ds::DynamicalSystem, p, pprev, pidx, next_id
) → rmap
```

Given dictionaries `a₊, a₋`

mapping IDs to values, return a *replacement map*: a dictionary mapping the IDs (keys) in dictionary `a₊`

to IDs (keys) in dictionary `a₋`

, so that so that values in `a₊`

that are the "closest" to values in `a₋`

get assigned the same key as in `a₋`

. In this way keys of `a₊`

are "matched" to keys of `a₋`

. Use `swap_dict_keys`

to apply `rmap`

to `a₊`

or to other dictionaries with same keys as `a₊`

.

How matching happens, i.e., how "closeness" is defined, depends on the algorithm `matcher`

.

The values contained in `a₊, a₋`

can be anything supported by `matcher`

. Within Attractors.jl they are typically `StateSpaceSet`

s representing attractors. Typically the +,- mean after and before some change of parameter of a dynamical system.

**Keyword arguments**

`ds`

: the dynamical system that generated`a₊, a₋`

.`p, pprev, pidx`

: the parameter values corresponding to`a₊, a₋`

and the index the parameter has in`ds`

.`next_id = next_free_id(a₊, a₋)`

: the ID to give to values of`a₊`

that cannot be matched to`a₋`

and hence must obtain a new unique ID.

Some matchers like `MatchBySSSetDistance`

do not utilize `ds, p, pprev, pidx`

in any way while other matchers like `MatchByBasinEnclosure`

do, and those require expliticly giving values to `ds, p, pprev, pidx`

as their default values is just `nothing`

.

`Attractors.mfs_brute_force`

— MethodThis function works on the results obtained by `crude_initial_radius`

. It starts from the best shock found so far and tries to find a better one by continuously reducing the radius of the sphere on the surface of which it generates random perturbations. If perturbation with the same basin of attraction is found, it updates the best shock found so far and reduces the radius of the sphere. It repeats this process total_iterations times and returns the best perturbation found.

`Attractors.minimal_fatal_shock`

— Method`minimal_fatal_shock(mapper::AttractorMapper, u0, search_area, algorithm; kw...)`

Return the *minimal fatal shock* `mfs`

(also known as *excitability threshold*) for the initial point `u0`

according to the specified `algorithm`

given a `mapper`

that satisfies the `id = mapper(u0)`

interface (see `AttractorMapper`

if you are not sure which mappers do that). The `mapper`

contains a reference to a `DynamicalSystem`

. The options for `algorithm`

are: `MFSBruteForce`

or `MFSBlackBoxOptim`

. For high dimensional systems `MFSBlackBoxOptim`

is likely more accurate.

The `search_area`

dictates the state space range for the search of the `mfs`

. It can be a 2-tuple of (min, max) values, in which case the same values are used for each dimension of the system in `mapper`

. Otherwise, it can be a vector of 2-tuples, each for each dimension of the system. The search area is defined w.r.t. to `u0`

(i.e., it is the search area for perturbations of `u0`

).

An alias to `minimal_fata_shock`

is `excitability_threshold`

.

**Keyword arguments**

`metric = LinearAlgebra.norm`

: a metric function that gives the norm of a perturbation vector. This keyword is ignored for the`MFSBruteForce`

algorithm.`target_id = nothing`

: when not`nothing`

, it should be an integer or a vector of integers corresponding to target attractor label(s). Then, the MFS is estimated based only on perturbations that lead to the target attractor(s).

**Description**

The minimal fatal shock is defined as the smallest-norm perturbation of the initial point `u0`

that will lead it a different basin of attraction. It is inspired by the paper "Minimal fatal shocks in multistable complex networks" Halekotte2020, however the implementation here is generic: it works for *any* dynamical system.

The *excitability threshold* is a concept nearly identical, however, instead of looking for a perturbation that simply brings us out of the basin, we look for the smallest perturbation that brings us into specified basin(s). This is enabled via the keyword `target_id`

.

`Attractors.next_free_id`

— Method`next_free_id(new::Dict, old::Dict)`

Return the minimum key of the "new" dictionary that doesn't exist in the "old" dictionary.

`Attractors.optimal_radius_dbscan_knee`

— MethodFind the optimal radius ϵ of a point neighborhood for use in DBSCAN through the elbow method (knee method, highest derivative method).

`Attractors.optimal_radius_dbscan_silhouette`

— MethodFind the optimal radius ε of a point neighborhood to use in DBSCAN, the unsupervised clustering method for `AttractorsViaFeaturizing`

. Iteratively search for the radius that leads to the best clustering, as characterized by quantifiers known as silhouettes. Does a linear (sequential) search.

`Attractors.optimal_radius_dbscan_silhouette_optim`

— MethodSame as `optimal_radius_dbscan_silhouette`

, but find minimum via optimization with Optim.jl.

`Attractors.optimal_radius_dbscan_silhouette_original`

— MethodFind the optimal radius ε of a point neighborhood to use in DBSCAN, the unsupervised clustering method for `AttractorsViaFeaturizing`

. The basic idea is to iteratively search for the radius that leads to the best clustering, as characterized by quantifiers known as silhouettes.

`Attractors.overwrite_dict!`

— Method`overwrite_dict!(old::Dict, new::Dict)`

In-place overwrite the `old`

dictionary for the key-value pairs of the `new`

.

`Attractors.plot_attractors`

— Function`plot_attractors(attractors::Dict{Int, StateSpaceSet}; kwargs...)`

Plot the attractors as a scatter plot.

**Keyword arguments**

- All the common plotting keywords. Particularly important is the
`access`

keyword. `sckwargs = (strokewidth = 0.5, strokecolor = :black,)`

: additional keywords propagated to the`Makie.scatter`

function that plots the attractors.

`Attractors.plot_attractors_curves`

— Function`plot_attractors_curves(attractors_cont, attractor_to_real, prange = 1:length(); kwargs...)`

Same as in `plot_basins_curves`

but visualizes the attractor dependence on the parameter instead of their fraction. The function `attractor_to_real`

takes as input a `StateSpaceSet`

(attractor) and returns a real number so that it can be plotted versus the parameter axis. See also `plot_basins_attractors_curves`

.

Same keywords as `plot_basins_curves`

. See also `plot_continuation_curves`

.

`Attractors.plot_basins_attractors_curves`

— Function```
plot_basins_attractors_curves(
fractions_cont, attractors_cont, a2rs [, prange]
kwargs...
)
```

Convenience combination of `plot_basins_curves`

and `plot_attractors_curves`

in a multi-panel plot that shares legend, colors, markers, etc. This function allows `a2rs`

to be a `Vector`

of functions, each mapping attractors into real numbers. Below the basins fractions plot, one additional panel is created for each entry in `a2rs`

. `a2rs`

can also be a single function, in which case only one panel is made.

`Attractors.plot_basins_curves`

— Function`plot_basins_curves(fractions_cont, prange = 1:length(); kwargs...)`

Plot the fractions of basins of attraction versus a parameter range, i.e., visualize the output of `global_continuation`

. See also `plot_basins_attractors_curves`

and `plot_continuation_curves`

.

**Keyword arguments**

`style = :band`

: how to visualize the basin fractions. Choices are`:band`

for a band plot with cumulative sum = 1 or`:lines`

for a lines plot of each basin fraction`separatorwidth = 1, separatorcolor = "white"`

: adds a line separating the fractions if the style is`:band`

`axislegend_kwargs = (position = :lt,)`

: propagated to`axislegend`

if a legend is added`series_kwargs = NamedTuple()`

: propagated to the band or scatterline plot- Also all common plotting keywords.

`Attractors.plot_continuation_curves`

— Function`plot_continuation_curves(continuation_info, prange = 1:length(); kwargs...)`

Same as in `plot_basins_curves`

but visualizes any arbitrary quantity characterizing the continuation. Hence, the `continuation_info`

is of exactly the same format as `fractions_cont`

: a vector of dictionaries, each dictionary mapping attractor IDs to real numbers. `continuation_info`

is meant to accompany `attractor_info`

in `plot_attractors_curves`

. To produce `continuation_info`

from `attractor_info`

you can do something like:

```
continuation_info = map(attractors_cont) do attractors
Dict(k => f(A) for (k, A) in attractors)
end
```

with `f`

your function of interest that returns a real number.

`Attractors.recurrences_map_to_label!`

— Method`recurrences_map_to_label!(bsn_nfo::BasinsInfo, ds, u0; kwargs...) -> ic_label`

Return the label of the attractor that the initial condition `u0`

converges to, or `-1`

if it does not convergence anywhere (e.g., divergence to infinity or exceeding `maximum_iterations`

).

Notice the numbering system `cell_label`

is as in `finite_state_machine!`

. Even numbers are attractors, odd numbers are basins.

`Attractors.reset_mapper!`

— Method`reset_mapper!(mapper::AttractorsMapper)`

Reset all accumulated information of `mapper`

so that it is as if it has just been initialized. Useful in `for`

loops that loop over a parameter of the dynamical system stored in `mapper`

.

`Attractors.retract_keys_to_consecutive`

— Method`retract_keys_to_consecutive(v::Vector{<:Dict}) → rmap`

Given a vector of dictionaries with various positive integer keys, retract all keys so that consecutive integers are used. So if the dictionaries have overall keys 2, 3, 42, then they will transformed to 1, 2, 3.

Return the replacement map used to replace keys in all dictionaries with `swap_dict_keys!`

.

As this function is used in attractor matching in `global_continuation`

it skips the special key `-1`

.

`Attractors.shaded_basins_heatmap`

— Function`shaded_basins_heatmap(grid, basins, attractors, iterations; kwargs...)`

Plot a heatmap of found (2-dimensional) `basins`

of attraction and corresponding `attractors`

. A matrix `iterations`

with the same size of `basins`

must be provided to shade the color according to the value of this matrix. A small value corresponds to a light color and a large value to a darker tone. This is useful to represent the number of iterations taken for each initial condition to converge. See also `convergence_time`

to store this iteration number.

**Keyword arguments**

`show_attractors = true`

: shows the attractor on plot`maxit = maximum(iterations)`

: clip the values of`iterations`

to

the value `maxit`

. Useful when there are some very long iterations and keep the range constrained to a given interval.

- All the common plotting keywords.

`Attractors.silhouettes_new`

— MethodCalculate silhouettes. A bit slower than the implementation in `Clustering.jl`

but seems to be more robust. The latter seems to be incorrect in some cases.

`Attractors.subdivision_based_grid`

— Method`subdivision_based_grid(ds::DynamicalSystem, grid; maxlevel = 4, q = 0.99)`

Construct a grid structure `SubdivisionBasedGrid`

that can be directly passed as a grid to `AttractorsViaRecurrences`

. The input `grid`

is an originally coarse grid (a tuple of `AbstractRange`

s). The state space speed is evaluate in all cells of the `grid`

. Cells with small speed (when compared to the "max" speed) resultin in this cell being subdivided more. To avoid problems with spikes in the speed, the `q`

-th quantile of the velocities is used as the "max" speed (use `q = 1`

for true maximum). The subdivisions in the resulting grid are clamped to at most value `maxlevel`

.

This approach is designed for *continuous time* systems in which different areas of the state space flow may have significantly different velocity. In case of originally coarse grids, this may lead `AttractorsViaRecurrences`

being stuck in some state space regions with a small motion speed and false identification of attractors.

`Attractors.swap_dict_keys!`

— Method`swap_dict_keys!(d::Dict, matching_map::Dict)`

Swap the keys of a dictionary `d`

given a `matching_map`

which maps old keys to new keys. Also ensure that a swap can happen at most once, e.g., if input `d`

has a key `4`

, and `rmap = Dict(4 => 3, 3 => 2)`

, then the key `4`

will be transformed to `3`

and not further to `2`

.

`Attractors.test_wada_merge`

— Method`test_wada_merge(basins, r) -> p`

Test if the 2D array `basins`

has the Wada property using the merging technique of Daza2018.

**Description**

The technique consists in computing the generalized basins of each attractor. These new basins are formed with on of the basins and the union of the other basins. A new boundary is defined by these two objects. The algorithm then computes the distance between each boundaries of these basins pairwise. If all the boundaries are within some distance `r`

, there is a unique boundary separating the basins and we have the wada property. The algorithm returns the maximum proportion of pixels of a boundary with distance strictly greater than `r`

from another boundary.

If `p == 0`

, we have the Wada property for this value of `r`

. If `p > 0`

, the criteria to decide if the basins are Wada is left to the user. Numerical inaccuracies may be responsible for a small percentage of points with distance larger than `r`

`Attractors.tipping_probabilities`

— Method`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 Kaszas2019.

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 Kaszas2019 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`

.

`Attractors.uncertainty_exponent`

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

Estimate the uncertainty exponentGrebogi1983 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 output `α`

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.

`Attractors.unique_keys`

— Method`unique_keys(v::Iterator{<:AbstractDict})`

Given a vector of dictionaries, return a sorted vector of the unique keys that are present across all dictionaries.