Attractors.AttractorsModule

Attractors.jl

docsdev docsstable Paper CI codecov Package Downloads

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.AttractorMapperType
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:

Attractors.AttractorsViaFeaturizingType
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 trajectoryA::StateSpaceSetand the corresponding time vectortand returns a vectorvof features describing the trajectory. For better performance, it is strongly recommended thatv 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.AttractorsViaProximityType
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 StateSpaceSets 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 StateSpaceSets. 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.AttractorsViaRecurrencesType
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:

  1. A tuple of sorted AbstractRanges for a regular grid.

Example is grid = (xg, yg) where xg = yg = range(-5, 5; length = 100) for a two-dimensional system.

  1. A tuple of sorted AbstractVectors 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).

  1. 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.EdgeTrackingResultsType
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 StateSpaceSets
  • 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.FeaturizeGroupAcrossParameterMethod
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.GroupViaClusteringType
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.GroupViaHistogramType
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.GroupViaNearestFeatureType
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.GroupViaPairwiseComparisonType
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.GroupingConfigType

GroupingConfig

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:

  1. Make a new type and subtype GroupingConfig.
  2. 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.
  3. 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.
  4. Include the new grouping file in the grouping/all_grouping_configs.jl and list it in this documentation string.
Attractors.MFSBlackBoxOptimType
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.MFSBruteForceType
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.RecurrencesFindAndMatchType
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.SubdivisionBasedGridType
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.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_fractionsFunction
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

  1. fractions_curves: a vector of dictionaries mapping labels to basin fractions.
  2. 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).
  3. featurizer: a 1-argument function to map an attractor into an appropriate feature to be grouped later. Features expected by GroupingConfig are SVector.
  4. group_config: a subtype of GroupingConfig.
  5. 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

  1. aggregated_fractions: same as fractions_curves but now contains the fractions of the aggregated attractors.
  2. 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_continuationFunction
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_basinsMethod
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_entropyFunction
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_dimensionMethod
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_testMethod
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_fractionsFunction
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_fractionsMethod
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_attractionMethod
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_attractionMethod
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_edgeMethod
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.

Info

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_assignmentMethod

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

Attractors.continuationFunction
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

  1. 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.
  2. 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.crude_initial_radiusMethod

This function generates a random perturbation of the initial point u0 within specified "searcharea" 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 totaliterations times and returns the best perturbation found.

Attractors.edgetrackingMethod
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 StateSpaceSets, 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.edgetrackingMethod
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_attractorsMethod
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_featuresMethod
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_groupMethod
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!Method

Main 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_featuresMethod
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.iterations_to_convergeFunction
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_forceMethod

This 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_shockMethod
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_silhouetteMethod

Find 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_originalMethod

Find 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_curvesFunction
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_curvesFunction
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.retract_keys_to_consecutiveMethod
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_heatmapFunction
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.

Attractors.silhouettes_newMethod

Calculate 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_gridMethod
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 AbstractRanges). 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_probabilitiesMethod
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_exponentMethod
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_keysMethod
unique_keys(v::Iterator{<:AbstractDict})

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