Attractor & Basins Continuation

A new kind of continuation

If you have heard before the word "continuation", then you are likely aware of the standard (linearized) continuation offered by many software, in Julia particularly the BifurcationKit.jl.

Here we offer a completely different kind of continuation called attractors & basins continuation.

A direct comparison of the two approaches is not truly possible, because they do different things. The traditional linearized continuation analysis continues the curves of individual fixed points across the joint state-parameter space. The attractor and basins continuation first finds all attractors at all parameter values and then matches appropriately similar attractors across different parameters, giving the illusion of continuing them individually. Additionally, the curves of stable fixed points in the joint parameter space is only a small by-product of the attractor basins continuation, and the main information is the basin fractions and how these change in the parameter space.

A more detailed comparison of the methods will be provided soon in a paper we are writing (will put link here, stay tuned!).

Continuation API

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.

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

Recurrences continuation (best)

Attractors.RecurrencesSeededContinuationType
RecurrencesSeededContinuation <: AttractorsBasinsContinuation
RecurrencesSeededContinuation(mapper::AttractorsViaRecurrences; kwargs...)

A method for rsc. TODO: Cite our preprint here.

Description

At the first parameter slice 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 piping 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 attractors are found, their fractions are computed by running new initial conditions through the AttractorsViaRecurrences mapper. This process continues until all parameter values are exhausted and for each parameter value the attractors and their fractions are found.

Then, the different attractors across parameters are matched so that they have the same ID. The matching process is based on distances between attractors. The function that computes these distances is setsofsets_distances and the matching function is match_attractor_ids! (please read those docstrings as well).

At each parameter slice beyond the first, the new attractors are matched to the previous attractors found in the previous parameter value by a direct call to the match_attractor_ids! function. Hence, the matching of attractors here works "slice by slice" on the parameter axis and the attractors that are closest to each other (in state space, but for two different parameter values) get assigned the same label.

Keyword arguments

  • distance, threshold: propagated to match_attractor_ids!.
  • 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 rsc.
  • 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 some points from existing attractors according to how many points the attractors themselves contain. A maximum of 10 seeds is done per attractor.
source

Matching attractors

Attractors.match_attractor_ids!Function
match_attractor_ids!(a₊::AbstractDict, a₋; distance = Centroid(), threshold = Inf)

Given dictionaries a₊, a₋ mapping IDs to attractors (StateSpaceSet instances), match attractor IDs in dictionary a₊ so that its attractors 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.

Description

When finding attractors and their fractions in DynamicalSystems.jl, different attractors get assigned different IDs. However which attractor gets which ID is somewhat arbitrary. Finding the attractors of the same system for slightly different parameters could label "similar" attractors (at the different parameters) with different IDs. match_attractors_ids! tries to "match" them by modifying the attractor IDs, i.e., the keys of the given dictionaries.

The matching happens according to the output of the setsofsets_distances function with the keyword distance. distance` can be whatever that function accepts. Attractors are then match according to distance, with unique mapping. The closest attractors (before and after) are mapped to each other, and are removed from the matching pool, and then the next pair with least remaining distance is matched, and so on.

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

source
StateSpaceSets.setsofsets_distancesFunction
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 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. However, distance can also be any arbitrary user function that takes as input two state space sets and returns any positive-definite number as their "distance".

Attractors.match_basins_ids!Function
match_basins_ids!(b₊::AbstractArray, b₋; threshold = Inf)

Similar to match_attractor_ids! 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).

source

Aggregating attractors and fractions

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 RecurrencesSeededContinuation 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 dimesional 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 RecurrencesSeededContinuation (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 GroupAcrossParameterContinuation. 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).

source

Grouping continuation

Attractors.GroupAcrossParameterContinuationType
GroupAcrossParameterContinuation(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.

source