`Attractors.Attractors`

— Module**Attractors.jl**

A Julia module for

- finding attractors of arbitrary dynamical systems
- finding their basins of attraction or the state space fractions of the basins
- analyzing global stability of attractors (also called non-local stability or resilience)
- "continuing" the attractors and their basins over a parameter range
- finding the basin boundaries 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. 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.

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

.`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 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.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 Examples 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 <: AttractorsBasinsContinuation
FeaturizeGroupAcrossParameter(mapper::AttractorsViaFeaturizing; kwargs...)
```

A method for `continuation`

. It uses the featurizing approach discussed in `AttractorsViaFeaturizing`

and hence requires an instance of that mapper as an input. When used in `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.

**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.`par_weight = 0`

: See below the section on MCBB.

**MCBB special version**

If the chosen grouping method is `GroupViaClustering`

, the additional keyword `par_weight::Real`

can be used. If it is ≠ 0, the distance matrix between features obtains an extra weight that is proportional to the distance `par_weight*|p[i] - p[j]|`

between the parameters used when extracting features. The range of parameters is normalized to 0-1 such that the largest distance in the parameter space is 1. The normalization is done because the feature space is also (by default) normalized to 0-1.

This version of the algorithm is the original "MCBB" continuation method described in Gelbrecht2020, besides the improvements of clustering accuracy and performance done by the developer team of Attractors.jl.

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

`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 initial conditions to IDs, then instead extend the
**internal function**`feature_to_group(feature, config)`

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

penalty value for the objective function, allows to adjust optimization algorithm to find the minimal fatal shock,`default = 1000.0`

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

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

— TypeAlias for `RecurrencesFindAndMatch`

.

`Attractors.RecurrencesFindAndMatch`

— Type```
RecurrencesFindAndMatch <: AttractorsBasinsContinuation
RecurrencesFindAndMatch(mapper::AttractorsViaRecurrences; kwargs...)
```

A method for `continuation`

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

) and the "matching attractors" functionality offered by `match_continuation!`

.

You can use `RAFM`

as an alias.

**Keyword arguments**

`distance = Centroid(), threshold = Inf, use_vanished = !isinf(threshold)`

: propagated to`match_continuation!`

.`info_extraction = identity`

: A function that takes as an input an attractor (`StateSpaceSet`

) and outputs whatever information should be stored. It is used to return the`attractors_info`

in`continuation`

. Note that the same attractors that are stored in`attractors_info`

are also used to perform the matching in`match_continuation!`

, hence this keyword should be altered with care.`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**

At the first parameter slice of the continuation process, attractors and their fractions are found as described in the `AttractorsViaRecurrences`

mapper using recurrences in state space. At each subsequent parameter slice, new attractors are found by seeding initial conditions from the previously found attractors and then running these initial conditions through the recurrences algorithm of the `mapper`

. Seeding initial conditions close to previous attractors accelerates the main bottleneck of `AttractorsViaRecurrences`

, which is finding the attractors.

After the special initial conditions are mapped to attractors, attractor basin fractions are computed by sampling random initial conditions using the provided `sampler`

in `continuation`

) and mapping them to attractors using the `AttractorsViaRecurrences`

mapper. 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". Once the basins fractions are computed, the parameter is incremented again and we perform the steps as before.

This process continues for all parameter values. After all parameters are exhausted, the found attractors (and their fractions) are "matched" to the previous ones. I.e., their *IDs are changed*, so that attractors that are "similar" to those at a previous parameter get assigned the same ID. Matching is done by the `match_continuation!`

function and is an *orthogonal* step. This means, that if you don't like the initial outcome of the matching process, you may call `match_continuation!`

again on the outcome with different matching-related keywords. You do not need to compute attractors and basins again!

Matching is a very sophisticated process that can be understood in detail by reading the docstrings of `match_statespacesets!`

first and then `match_continuation!`

. Here is a short summary: attractors from previous and current parameter are matched based on their "distance". By default this is distance in state space, but any measure of "distance" may be used, such as the distance between Lyapunov spectra. Matching prioritizes new->old pairs with smallest distance: once these are matched the next available new->old pair with smallest distance is matched, until all new/old attractors have been matched. The `threshold`

keyword establishes that attractors with distance > `threshold`

do not get matched. Additionally, use `use_vanished = true`

if you want to include as matching candidates attractors that have vanished during the continuation process.

`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._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_curves, attractors_info, 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 `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_curves`

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

: a vector of dictionaries mapping labels to attractors. 1st and 2nd argument are exactly like the return values of`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_curves`

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_info, fractions_curves, 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 `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_info, fractions_curves`

are the output of `continuation`

while `ds, prange, pidx`

are the input to `continuation`

.

**Keyword arguments**

`savename = "test.mp4"`

: name of video output file`framerate = 4`

: framerate of video output`markersize = 10`

`Δt, T`

: propagated to`trajectory`

for evolving an initial condition sampled from an attractor- Also all common plotting keywords.

`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) → 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. For this, see `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 *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`

.

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

.

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

— Function`continuation(abc::AttractorsBasinsContinuation, prange, pidx, ics; kwargs...)`

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

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

The continuation type `abc`

is a subtype of `AttractorsBasinsContinuation`

and contains an `AttractorMapper`

. The mapper contains information on how to find the attractors and basins of a dynamical system. Additional arguments and keyword arguments given when creating `abc`

further tune the continuation and how attractors are matched across different parameter values.

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 `set_parameter!`

can be used).

`ics`

is a 0-argument function generating initial conditions for the dynamical system (as in `basins_fractions`

).

Possible subtypes of `AttractorsBasinsContinuation`

are:

**Return**

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

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

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

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

. Information about the attractors.`attractors_info[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 continuation type.

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

.

`Attractors.convergence_and_basins_of_attraction`

— Method`convergence_and_basins_of_attraction(mapper::AttractorsViaRecurrences, grid) -> basins, attractors, iterations`

An extension of `basins_of_attraction`

. This version also returns `iterations`

, which is the number of iterations each initial condition took to converge to the attractor. Useful to give to `shaded_basins_heatmap`

.

**Keyword arguments**

`show_progress = true`

: show progress bar

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

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

Group the given vector of feature vectors according to the configuration and return the labels (vector of equal length as `features`

). See `AttractorsViaFeaturizing`

for possible configurations.

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

— Function`iterations_to_converge(mapper::AttractorMapper) → n::Int`

Return the number of iterations 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.

`Attractors.match_basins_ids!`

— Method`match_basins_ids!(b₊::AbstractArray, b₋; threshold = Inf)`

Similar to `match_statespacesets!`

but operate on basin arrays instead (the arrays typically returned by `basins_of_attraction`

).

This method 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% overlap of pixels). The `threshold`

in this case is compared to the inverse of the overlap (so, for `threshold = 2`

attractors that have less than 50% overlap get different IDs guaranteed).

`Attractors.match_continuation!`

— Method`match_continuation!(fractions_curves::Vector{<:Dict}, attractors_info::Vector{<:Dict}; kwargs...)`

Loop over all entries in the given arguments (which are typically the direct outputs of `continuation`

with `RecurrencesFindAndMatch`

), and match the attractor IDs in both the attractors container and the basins fractions container. This means that we loop over each entry of the vectors (skipping the first), and in each entry we attempt to match the current dictionary keys to the keys of the previous dictionary using `match_statespacesets!`

.

The keywords `distance, threshold`

are propagated to `match_statespacesets!`

. However, there is a unique keyword for `match_continuation!`

: `use_vanished::Bool`

. If `true`

, then attractors that existed before but have vanished are kept in "memory" when it comes to matching: the new attractors are compared to the latest instance of all attractors that have ever existed, and get matched to their closest ones as per `match_statespacesets!`

. If `false`

, vanished attractors are ignored. Note that in this case new attractors that cannot be matched to any previous attractors will get an appropriately incremented ID. E.g., if we started with three attractors, and attractor 3 vanished, and at some later parameter value we again have three attractors, the new third attractor will *not* have ID 3, but 4 (i.e., the next available ID).

By default `use_vanished = !isinf(threshold)`

and since the default value for `threshold`

is `Inf`

, `use_vanished`

is `false`

.

The last keyword is `retract_keys = true`

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

`match_continuation!(attractors_info::Vector{<:Dict}; kwargs...)`

This is a convenience method that only uses and modifies the state space set dictionary container without the need for a basins fractions container.

`Attractors.match_statespacesets!`

— Method`match_statespacesets!(a₊::AbstractDict, a₋; distance = Centroid(), threshold = Inf)`

Given dictionaries `a₊, a₋`

mapping IDs to `StateSpaceSet`

instances, match the IDs in dictionary `a₊`

so that its sets that are the closest to those in dictionary `a₋`

get assigned the same key as in `a₋`

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

Return the replacement map, a dictionary mapping old keys of `a₊`

to the new ones that they were mapped to. You can obtain this map, without modifying the dictionaries, by calling the `replacement_map`

function directly.

**Keyword arguments**

`distance = Centroid()`

: given to`setsofsets_distances`

.`threshold = Inf`

: attractors with distance larger than the`threshold`

are guaranteed to not be mapped to each other.

**Description**

The distance between all possible pairs of sets between the "old" and "new" containers is computed using `setsofsets_distances`

with the keyword `distance`

. `distance`

can be whatever that function accepts, i.e., one of `Centroid, Hausdorff, StrictlyMinimumDistance`

or any arbitrary user- provided function that given two sets it returns a positive number (their distance). State space sets are then matched according to this distance. First, all possible pairs (old, new, distance) are sorted according to their distance. The pair with smallest distance is matched. Sets 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 attractors is larger than this `threshold`

, then it is guaranteed that the attractors will get assigned different key in the dictionary `a₊`

(which is the next available integer).

`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) → mfs`

Return the minimal fatal shock `mfs`

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`

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

).

**Description**

The minimal fatal shock is defined as the smallest (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.

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

— Function`plot_attractors_curves(attractors_info, 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.

Same keywords as `plot_basins_curves`

.

`Attractors.plot_basins_attractors_curves`

— Function```
plot_basins_attractors_curves(
fractions_curves, attractors_info, attractor_to_real [, prange]
kwargs...
)
```

Convenience combination of `plot_basins_curves`

and `plot_attractors_curves`

in a two-panel plot that shares legend, colors, markers, etc.

`Attractors.plot_basins_curves`

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

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

.

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

— Method`replacement_map(a₊, a₋; distance = Centroid(), threshold = Inf) → rmap`

Return a dictionary mapping keys in `a₊`

to new keys in `a₋`

, as explained in `match_statespacesets!`

.

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

. An matrix `iterations`

with the same size of `basins`

must be provided to shade the color according to the value of this matrix. An 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 conditions. See also [`iterations_to_converge`

]@ref 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, replacement_map::Dict)`

Swap the keys of a dictionary `d`

given the `replacement_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.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.