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.continuation
— Functioncontinuation(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
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 thei
-th parameter.attractors_info <: Vector{Dict{Int, <:Any}}
. Information about the attractors.attractors_info[i]
is a dictionary mapping attractor ID to information about the attractor at thei
-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 fromics
.
Recurrences continuation (best)
Attractors.RecurrencesSeededContinuation
— TypeRecurrencesSeededContinuation <: 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 tomatch_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 theattractors_info
inrsc
.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 of10
seeds is done per attractor.
Matching attractors
Attractors.match_attractor_ids!
— Functionmatch_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₊
.
Attractors.replacement_map
— Functionreplacement_map(a₊, a₋; distance = Centroid(), threshold = Inf) → rmap
Return a dictionary mapping keys in a₊
to new keys in a₋
, as explained in match_attractor_ids!
.
StateSpaceSets.set_distance
— Functionset_distance(ssset1, ssset2 [, distance])
Calculate a distance between two StateSpaceSet
s, 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 restHausdorff
StrictlyMinimumDistance
StateSpaceSets.setsofsets_distances
— Functionsetsofsets_distances(a₊, a₋ [, distance]) → distances
Calculate distances between sets of StateSpaceSet
s. Here a₊, a₋
are containers of StateSpaceSet
s, 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!
— Functionmatch_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).
Aggregating attractors and fractions
Attractors.aggregate_attractor_fractions
— Functionaggregate_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
fractions_curves
: a vector of dictionaries mapping labels to basin fractions.attractors_info
: a vector of dictionaries mapping labels to attractors. 1st and 2nd argument are exactly like the return values ofcontinuation
withRecurrencesSeededContinuation
(or, they can be the return ofbasins_fractions
).featurizer
: a 1-argument function to map an attractor into an appropriate feature to be grouped later. Features expected byGroupingConfig
areSVector
.group_config
: a subtype ofGroupingConfig
.info_extraction
: a function accepting a vector of features and returning a description of the features. I.e., exactly as inGroupAcrossParameterContinuation
. The 5th argument is optional and defaults to the centroid of the features.
Return
aggregated_fractions
: same asfractions_curves
but now contains the fractions of the aggregated attractors.aggregated_info
: dictionary mapping the new labels ofaggregated_fractions
to the extracted information usinginfo_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 StateSpaceSet
s and outputs a descriptor. E.g., info_extraction = vector -> mean(mean(x) for x in vector)
.
Grouping continuation
Attractors.GroupAcrossParameterContinuation
— TypeGroupAcrossParameterContinuation(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.
- Gelbrecht2021Maximilian Gelbrecht et al 2021, Monte Carlo basin bifurcation analysis, New J. Phys.22 03303