Attractor & Basins Continuation

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 (Datseris et al., 2023).

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:

• RecurrencesFindAndMatch
• FeaturizeGroupAcrossParameter

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.

RecurrencesFindAndMatch <: AttractorsBasinsContinuation
RecurrencesFindAndMatch(mapper::AttractorsViaRecurrences; kwargs...)

A method for continuation as in (Datseris et al., 2023) 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.

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

Centroid(metric = Euclidean())

A distance that can be used in set_distance. The Centroid method returns the distance (according to metric) between the centroids (a.k.a. centers of mass) of the sets.

metric can be any function that takes in two static vectors are returns a positive definite number to use as a distance (and typically is a Metric from Distances.jl).

Hausdorff(metric = Euclidean())

A distance that can be used in set_distance. The Hausdorff distance is the greatest of all the distances from a point in one set to the closest point in the other set. The distance is calculated with the metric given to Hausdorff which defaults to Euclidean.

Hausdorff is 2x slower than StrictlyMinimumDistance, however it is a proper metric in the space of sets of state space sets.

StrictlyMinimumDistance([brute = false,] [metric = Euclidean(),])

A distance that can be used in set_distance. The StrictlyMinimumDistance returns the minimum distance of all the distances from a point in one set to the closest point in the other set. The distance is calculated with the given metric.

The brute::Bool argument switches the computation between a KDTree-based version, or brute force (i.e., calculation of all distances and picking the smallest one). Brute force performs better for datasets that are either large dimensional or have a small amount of points. Deciding a cutting point is not trivial, and is recommended to simply benchmark the set_distance function to make a decision.

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

set_distance(ssset1, ssset2 [, distance])

Calculate a distance between two StateSpaceSets, i.e., a distance defined between sets of points, as dictated by distance.

Possible distance types are:

• Centroid, which is the default, and 100s of times faster than the rest
• Hausdorff
• StrictlyMinimumDistance

setsofsets_distances(a₊, a₋ [, distance]) → distances

Calculate distances between sets of StateSpaceSets. Here a₊, a₋ are containers of StateSpaceSets, and the returned distances are dictionaries of distances. Specifically, distances[i][j] is the distance of the set in the i key of a₊ to the j key of a₋. Notice that distances from a₋ to a₊ are not computed at all (assumming symmetry in the distance function).

The distance can be as in set_distance, or it can be an arbitrary function that takes as input two state space sets and returns any positive-definite number as their "distance".

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), 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 are two unique keywords for match_continuation!:

• use_vanished::Bool
• retract_keys::Bool

If use_vanised = 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.

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

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

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 (Gelbrecht et al., 2020), besides the improvements of clustering accuracy and performance done by the developer team of Attractors.jl.