Attractors.Attractors
— ModuleAttractors.jl
Attractors.jl is a Julia package for
- Finding all attractors, and all types of attractors, of arbitrary dynamical systems. An extendable interface allows for new algorithms for finding attractors.
- Finding their basins of attraction or the state space fractions of the basins. This includes finding exit basins (divergence to infinity).
- Analyzing nonlocal stability of attractors (also called global stability or resilience).
- Performing global continuation of attractors and their basins (or other measures of stability), over a parameter range. Global continuation is a new, cutting-edge type of continuation that offers several advantages over traditional local continuation (AUTO, MatCont, BifurcationKit.jl, etc.), see the comparison in our docs.
- Finding the basin boundaries and edges states and analyzing their fractal properties.
- Tipping points related functionality for systems with known dynamic rule.
- And more!
It can be used as a standalone package, or as part of DynamicalSystems.jl.
To install it, run import Pkg; Pkg.add("Attractors")
.
All further information is provided in the documentation, which you can either find online or build locally by running the docs/make.jl
file.
Previously, Attractors.jl was part of ChaosTools.jl
Attractors.AttractorMapper
— TypeAttractorMapper(ds::DynamicalSystem, args...; kwargs...) → mapper
Subtypes of AttractorMapper
are structures that map initial conditions of ds
to attractors. The found attractors are stored inside the mapper
, and can be obtained by calling attractors = extract_attractors(mapper)
.
Currently available mapping methods:
All AttractorMapper
subtypes can be used with basins_fractions
or basins_of_attraction
.
In addition, some mappers can be called as a function of an initial condition:
label = mapper(u0)
and this will on the fly compute and return the label of the attractor u0
converges at. The mappers that can do this are:
AttractorsViaProximity
AttractorsViaRecurrences
AttractorsViaFeaturizing
with theGroupViaHistogram
configuration.
For developers
AttractorMapper
defines an extendable interface. A new type needs to implement extract_attractors
and id = mapper(u0)
. From these, everything else in the entire rest of the library just works! If it is not possible to implement id = mapper(u0)
, then instead extend basins_fractions(mapper, ics)
.
Attractors.AttractorSeedContinueMatch
— TypeAttractorSeedContinueMatch(mapper, matcher = MatchBySSSetDistance(); seeding)
A global continuation method for global_continuation
. mapper
is any subtype of AttractorMapper
which implements extract_attractors
, i.e., it finds the actual attractors. matcher
is a configuration of how to match attractor IDs, see IDMatcher
for more options.
Description
This is a general/composable global continuation method based on a 4-step process:
- Seed initial conditions from previously found attractors
- Propagate those forwards to "continue" previous attractors
- Estimate basin fractions and potentially find new attractors
- Match attractors
Step 0 - Finding initial attractors
At the first parameter slice of the global continuation process, attractors and their fractions are found using the given mapper
and basins_fractions
. See the mapper
documentation and AttractorMapper
for details on how this works. Then, from the second parameter onwards the continuation occurs.
Step 1 - Seeding initial conditions
Initial conditions can be seeded from previously found attractors. This is controlled by the seeding
keyword, which must be a function that given a StateSpaceSet
(an attractor), it returns an iterator of initial conditions. By default the first point of an attractor is provided as the only seed.
Seeding can be turned off by providing the dummy function seeding = A -> []
, i.e., it always returns an empty iterator and hence no seeds and we skip to step 2.
Step 2 - Continuing the seeds
The dynamical system referenced by the mapper
is now set to the new parameter value. The seeds are run through the mapper
to converge to attractors at the new parameter value. Seeding initial conditions close to previous attractors increases the probability that if an attractor continues to exist in the new parameter, it is found. Additionally, for some mappers
this seeding process improves the accuracy as well as performance of finding attractors, see e.g. discussion in Datseris2023.
This seeding works for any mapper
, regardless of if they can map individual initial conditions with the mapper(u0)
syntax! If this syntax isn't supported, steps 2 and 3 are done together.
Step 3 - Estimate basins fractions
After the special seeded initial conditions are mapped to attractors, attractor basin fractions are computed by sampling additional initial conditions using the provided ics
in global_continuation
. 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".
Step 4 - Matching
Normally the ID an attractor gets assigned is somewhat a random integer. Therefore, to ensure a logical output of the global continuation process, attractors need to be "matched". This means: attractor and fractions must have their IDs changed, so that attractors that are "similar" to those at a previous parameter get assigned the same ID.
What is "similar enough" is controlled by the matcher
input. The default matcher
MatchBySSSetDistance
matches sets which have small distance in state space. The matching algorithm itself can be quite involved, so read the documentation of the matcher
for how matching works.
A note on matching: the MatchBySSSetDistance
can also be used after the continuation is completed, as it only requires as input the state space sets (attractors), without caring at which parameter each attractor exists at. If you don't like the final matching output, you may use a different instance of MatchBySSSetDistance
and call match_sequentially!
again on the output, without having to recompute the whole global continuation!
Step 5 - Finish
After matching the parameter is incremented. Steps 1-4 repeat until all parameter values are exhausted.
Further note
This global continuation method is a generalization of the "RAFM" continuation described in Datseris2023. This continuation method is still exported as RecurrencesFindAndMatch
.
Attractors.AttractorsViaFeaturizing
— TypeAttractorsViaFeaturizing(
ds::DynamicalSystem, featurizer::Function,
grouping_config = GroupViaClustering(); kwargs...
)
Initialize a mapper
that maps initial conditions to attractors using a featurizing and grouping approach. This is a supercase of the featurizing and clustering approach that is utilized by bSTAB Stender2021 and MCBB Gelbrecht2020. See AttractorMapper
for how to use the mapper
. This mapper
also allows the syntax mapper(u0)
but only if the grouping_config
is not GroupViaClustering
.
featurizer
is a function f(A, t)
that takes as an input an integrated trajectory A::StateSpaceSet
and the corresponding time vector t
and returns a vector v
of features describing the trajectory. For better performance, it is strongly recommended that v isa SVector{<:Real}
.
grouping_config
is an instance of any subtype of GroupingConfig
and decides how features will be grouped into attractors, see below.
See also the intermediate functions extract_features
and group_features
, which can be utilized when wanting to work directly with features.
Keyword arguments
T=100, Ttr=100, Δt=1
: Propagated toDynamicalSystems.trajectory
for integrating an initial condition to yieldA, t
.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 feature 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 are stored and can be accessed with extract_attractors
, however it should be clear that this mapper never actually finds attractors. They way we store attractors is by picking the first initial condition that belongs to the corresponding "attractor group", and then recording its trajectory with the same arguments T, Ttr, Δt
. This is stored as the attractor, but of course there is no guarantee that this is actually an attractor.
Attractors.AttractorsViaProximity
— TypeAttractorsViaProximity(ds::DynamicalSystem, attractors::Dict [, ε]; kwargs...)
Map initial conditions to attractors based on whether the trajectory reaches ε
-distance close to any of the user-provided attractors
. They have to be in a form of a dictionary mapping attractor labels to StateSpaceSet
s containing the attractors.
The system gets stepped, and at each step the minimum distance to all attractors is computed. If any of these distances is < ε
, then the label of the nearest attractor is returned.
If an ε::Real
is not provided by the user, a value is computed automatically as half of the minimum distance between all attractors. This operation can be expensive for large StateSpaceSet
s. If length(attractors) == 1
, then ε
becomes 1/10 of the diagonal of the box containing the attractor. If length(attractors) == 1
and the attractor is a single point, an error is thrown.
Keywords
Ttr = 100
: Transient time to first evolve the system for before checking for proximity.Δt = 1
: Step time given tostep!
.horizon_limit = 1e3
: If the maximum distance of the trajectory from any of the given attractors exceeds this limit, it is assumed that the trajectory diverged (gets labelled as-1
).consecutive_lost_steps = 1000
: If the integrator has been stepped this many times without comingε
-near to any attractor, it is assumed that the trajectory diverged (gets labelled as-1
).
Attractors.AttractorsViaRecurrences
— TypeAttractorsViaRecurrences(ds::DynamicalSystem, grid; kwargs...)
Map initial conditions of ds
to attractors by identifying attractors on the fly based on recurrences in the state space, as outlined in Datseris2022. However, the Description section below for has a more accurate (and simpler) exposition to the algorithm than the paper.
grid
is instructions for partitioning the state space into finite-sized cells so that a finite state machine can operate on top of it. Possibilities are:
- A tuple of sorted
AbstractRange
s for a regular grid.
Example is grid = (xg, yg)
where xg = yg = range(-5, 5; length = 100)
for a two-dimensional system.
- A tuple of sorted
AbstractVector
s for an irregular grid, for example
grid = (xg, yg)
with xg = range(0, 10.0^(1/2); length = 200).^2, yg = range(-5, 5; length = 100)
.
- An instance of the special grid type
SubdivisionBasedGrid
, which can be created either manually or by using subdivision_based_grid
. This automatically analyzes and adapts grid discretization levels in accordance with state space flow speed in different regions.
The grid has to be the same dimensionality as the state space, use a ProjectedDynamicalSystem
if you want to search for attractors in a lower dimensional subspace.
Keyword arguments
sparse = true
: control the storage type of the state space grid. If true, uses a sparse array, whose memory usage is in general more efficient than a regular array obtained withsparse=false
. In practice, the sparse representation should always be preferred when searching forbasins_fractions
. Only for very low dimensional systems and for computing the fullbasins_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 thestep!
function). The keywordDt
can also be used instead ifΔ
(\Delta
) is not accessible. It is1
for discrete time systems. For continuous systems, an automatic value is calculated usingautomatic_Δ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, useforce_non_adaptive = true
.force_non_adaptive = false
: Only used if the input dynamical system isCoupledODEs
. Iftrue
the additional keywordsadaptive = false, dt = Δt
are given asdiffeq
to theCoupledODEs
. 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. Onceattractor_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. Iftrue
, each visited cell will only store a point once, which is desirable for fixed points and limit cycles. Iffalse
thenattractor_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 ifsparse = 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 forattractor_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 whensparse = 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 > thanhorizon_limit
: the initial condition is labelled-1
. - If none of the above happens, the initial condition is labelled
-1
aftermaximum_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 documentation, under the recurrences animation page.
Attractors.EdgeTrackingResults
— TypeEdgeTrackingResults(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 betweentrack1
andtrack2
)track1::StateSpaceSet
: the pseudo-trajectory tracking the edge within basin 1track2::StateSpaceSet
: the pseudo-trajectory tracking the edge within basin 2time::Vector
: time points of the aboveStateSpaceSet
sbisect_idx::Vector
: indices oftime
at which a re-bisection occurredsuccess::Bool
: indicates whether the edge tracking has been successful or not
Attractors.FeaturizeGroupAcrossParameter
— MethodFeaturizeGroupAcrossParameter <: GlobalContinuationAlgorithm
FeaturizeGroupAcrossParameter(mapper::AttractorsViaFeaturizing; kwargs...)
A method for global_continuation
. It uses the featurizing approach discussed in AttractorsViaFeaturizing
and hence requires an instance of that mapper as an input. When used in global_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.
This continuation method is based on, but strongly generalizes, the approaches in the papers Gelbrecht2020 and Stender2021.
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. This is what theattractors_cont
contains in the return ofglobal_continuation
.
Attractors.GlobalContinuationAlgorithm
— TypeGlobalContinuationAlgorithm
Supertype of all algorithms used in global_continuation
. Each algorithm typically references an AttractorMapper
, as well as contains more information for how to continue/track/match attractors across a parameter range.
See global_continuation
for more.
Attractors.GroupViaClustering
— TypeGroupViaClustering(; 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 functionf(a, b)
that returns the distance between real-valued vectorsa, 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 withsilhouette_statistic
. The linear search may take some time to finish. To increase speed, the number of radii iterated through can be reduced by decreasingnum_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 thatmin_neighbors
> 1.
num_attempts_radius = 100
: number of radii that theoptimal_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 theminimum
of the silhouettes, and typically performs less accurately than themean
.max_used_features = 0
: if not0
, it should be anInt
denoting the max amount of features to be used when finding the optimal radius. Useful when clustering a very large number of features (e.g., high accuracy estimation of fractions of basins of attraction).
Description
The DBSCAN clustering algorithm is used to automatically identify clusters of similar features. Each feature vector is a point in a feature space. Each cluster then basically groups points that are closely packed together. Closely packed means that the points have at least min_neighbors
inside a ball of radius optimal_radius
centered on them. This method typically works well if the radius is chosen well, which is not necessarily an easy task. Currently, three methods are implemented to automatically estimate an "optimal" radius.
Estimating the optimal radius
The default method is the silhouettes method, which includes keywords silhouette
and silhouette_optim
. Both of them search for the radius that optimizes the clustering, meaning the one that maximizes a statistic silhouette_statistic
(e.g. mean value) of a quantifier for the quality of each cluster. This quantifier is the silhouette value of each identified cluster. A silhouette value measures how similar a point is to the cluster it currently belongs to, compared to the other clusters, and ranges from -1 (worst matching) to +1 (ideal matching). If only one cluster is found, the assigned silhouette is zero. So for each attempted radius in the search the clusters are computed, their silhouettes calculated, and the statistic of these silhouettes computed. The algorithm then finds the radius that leads to the maximum such statistic. For optimal_radius_method = "silhouettes"
, the search is done linearly, from a minimum to a maximum candidate radius for optimal_radius_method = "silhouettes"
; optimal_radius_method = silhouettes_optim
, it is done via an optimized search performed by Optim.jl which is typically faster and with similar accuracy. A third alternative is the"elbow"
method, which works by calculating the distance of each point to its k-nearest-neighbors (with k=min_neighbors
) and finding the distance corresponding to the highest derivative in the curve of the distances, sorted in ascending order. This distance is chosen as the optimal radius. It is described in Ester1996 and Schubert2017. It typically performs considerably worse than the "silhouette"
methods.
Attractors.GroupViaHistogram
— TypeGroupViaHistogram(binning::FixedRectangularBinning)
Initialize a struct that contains instructions on how to group features in AttractorsViaFeaturizing
. GroupViaHistogram
performs a histogram in feature space. Then, all features that are in the same histogram bin get the same label. The binning
is an instance of FixedRectangularBinning
from ComplexityMeasures.jl. (the reason to not allow RectangularBinning
is because during continuation we need to ensure that bins remain identical).
Attractors.GroupViaNearestFeature
— TypeGroupViaNearestFeature(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 exceedmax_distance
to their nearest template get labelled-1
.
Attractors.GroupViaPairwiseComparison
— TypeGroupViaPairwiseComparison(; 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 functionmetric(a, b)
that returns the distance between two featuresa
andb
, outputs offeaturizer
. AnyMetric
from Distances.jl can be used here.rescale_features = true
: if true, rescale each dimension of the extracted features separately into the range[0,1]
. This typically leads to more accurate grouping.
Description
This algorithm assumes that the features are well-separated into distinct clouds, with the maximum radius of the cloud controlled by threshold
. Since the systems are deterministic, this is achievable with a good-enough featurizer
function, by removing transients, and running the trajectories for sufficiently long. It then considers that features belong to the same attractor when their pairwise distance, computed using metric
, is smaller than or equal to threshold
, and that they belong to different attractors when the distance is bigger. Attractors correspond to each grouping of similar features. In this way, the key parameter threshold
is basically the amount of variation permissible in the features belonging to the same attractor. If they are well-chosen, the value can be relatively small and does not need to be fine tuned.
The threshold
should achieve a balance: one one hand, it should be large enough to account for variations in the features from the same attractor - if it's not large enough, the algorithm will find duplicate attractors. On the other hand, it should be small enough to not group together features from distinct attractors. This requires some knowledge of how spread the features are. If it's too big, the algorithm will miss some attractors, as it groups 2+ distinct attractors together. Therefore, as a rule of thumb, one can repeat the procedure a few times, starting with a relatively large value and reducing it until no more attractors are found and no duplicates appear.
The method uses relatively little memory, as it only stores vectors whose size is on order of the number of attractors of the system.
Attractors.GroupingConfig
— TypeGroupingConfig
Supertype for configuration structs on how to group features together. Used in several occasions such as AttractorsViaFeaturizing
or aggregate_attractor_fractions
.
Currently available grouping configurations are:
For developers
GroupingConfig
defines an extendable interface. The only thing necessary for a new grouping configuration is to:
- Make a new type and subtype
GroupingConfig
. - If the grouping allows for mapping individual features to group index, then instead extend the internal function
feature_to_group(feature, config)
. This will also allow doingid = mapper(u0)
withAttractorsViaFeaturizing
. - Else, extend the function
group_features(features, config)
. You could still extendgroup_features
even if (2.) is satisfied, if there are any performance benefits. - Include the new grouping file in the
grouping/all_grouping_configs.jl
and list it in this documentation string.
Attractors.IDMatcher
— TypeIDMatcher
Supertype of all "matchers" that match can IDs labelling attractors. Currently available matchers:
Matchers implement an extendable interface based on the function matching_map
. This function is used by the higher level function match_sequentially!
, which can be called after any call to a global continuation to match attractors differently, if the matching used originally during the continuation was not the best.
Attractors.MFSBlackBoxOptim
— TypeMFSBlackBoxOptim(; 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 notnothing
,random_algo
below is ignored.max_steps = 10000
: maximum number of steps for the optimization algorithm.penalty = 1000.0
: penalty value for the objective function for perturbations that do not lead to a different basin of attraction. This value is added to the norm of the perturbation and its value should be much larger than the typical sizes of the basins of attraction.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 ofMFSBruteForce
that can be used to provide an initial guess.bbkwargs = NamedTuple()
: additional keyword arguments propagated toBlackBoxOptim.bboptimize
for selecting solver, accuracy, and more.
Description
The algorithm uses BlackBoxOptim.jl and a penalized objective function to minimize. y function used as a constraint function. So, if we hit another basin during the search we encourage the algorithm otherwise we punish it with some penalty. The function to minimize is (besides some details):
function mfs_objective(perturbation, u0, mapper, penalty)
dist = norm(perturbation)
if mapper(u0 + perturbation) == mapper(u0)
# penalize if we stay in the same basin:
return dist + penalty
else
return dist
end
end
Using an initial guess can be beneficial to both performance and accuracy, which is why the output of a crude MFSBruteForce
is used to provide a guess. This can be disabled by either passing a guess
vector explicitly or by giving nothing
as random_algo
.
Attractors.MFSBruteForce
— TypeMFSBruteForce(; 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.MatchByBasinEnclosure
— TypeMatchByBasinEnclosure(; kw...) <: IDMatcher
A matcher that matches attractors by whether they are enclosed in the basin of a new attractor or not.
Keyword arguments
ε = nothing
: distance threshold given toAttractorsViaProximity
. Ifnothing
, it is estimated as a quarter of the minimum distance of centroids (in contrast to the default more accurate estimation inAttractorsViaProximity
).Δt = 1, consecutive_lost_steps = 1000
: also given toAttractorsViaProximity
. We have not yet decided what should happen to attractors that did not converge to one of the current attractors within this number of steps. At the moment they get assigned the next available free ID but this may change in future releases.distance = Centroid()
: metric to estimate distances between state space sets in case there are co-flowing attractors, see below.seeding = A -> A[end]
: how to select a point from the attractor to see if it is enclosed in the basin of a new attractor.
Description
An attractor A₋
is a set in a state space that occupies a particular region (or, a single point, if it is a fixed point). This region is always within the basin of attraction of said attractor. When the parameter of the dynamical system is incremented, the attractors A₊
in the new parameter have basins that may have changed in shape and size.
The new attractor A₊
is "matched" (i.e., has its ID changed) to the old attractor A₋
attractor if A₋
is located inside the basin of attraction of A₊
. To see if A₋
is in the basin of A₊
, we first pick a point from A₋
using the seeding
keyword argument. By default this is the last point on the attractor, but it could be anything else, including the centroid of the attractor (mean(A)
). This point is given as an initial condition to an AttractorsViaProximity
mapper that maps initial conditions to the ₊
attractors when the trajectories from the initial conditions are ε
-close to the ₊
attractors.
There can be the situation where multiple ₋
attractors converge to the same ₊
attractor, which we call "coflowing attractors". In this scenario matching is prioritized for the ₋
attractor that is closest to the ₊
in terms of state space set distance, which is estimated with the distance
keyword, which can be anything MatchBySSSetDistance
accepts. The closest ₊
attractor gets the ID of the ₋
closest attractor that converge to it.
Basin enclosure is a concept similar to "basin (in)stability" in Ritchie2023: attractors that quantify as "basin stable" are matched.
Attractors.MatchByBasinOverlap
— TypeMatchByBasinOverlap(threshold = Inf)
A matcher that matches IDs given full basins of attraction.
Description
This matcher cannot be used in with the generic global continuation method of AttractorSeedContinueMatch
. This matcher 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% of a basin in b₋
is overlaping with some other basin in b₊
). Therefore, the values this matcher compares are full basins of attraction, not attractors themselves (hence why it can't be given to AttractorSeedContinueMatch
). Rather, you may use this matcher with matching_map
.
The threshold
can dissallow matching between basins that do not have enough overlap. Basins whose overlap is less than 1/threshold
are guaranteed to get assined different IDs. For example: for threshold = 2
basins that have ≤ 50% overlap get different IDs guaranteed. By default, there is no threshold.
The information of the basins of attraction is typically an Array
, i.e., the direct output of basins_of_attraction
. For convenience, as well as backwards compatibility, when using matching_map
with this mapper you may provide two Array
s b₊, b₋
representing basins of attraction after and before, and the conversion to dictionaries will happen internally as it is supposed to. To replace the IDs
in b₊
given the replacement map just call replace!(b₊, rmap...)
, or use the in-place version matching_map!
directly.
A lower-level input for this matcher in matching_map
can be dictionaries mapping IDs to vectors of cartesian indices, where the indices mean which parts of the state space belong to which ID
Attractors.MatchBySSSetDistance
— TypeMatchBySSSetDistance(; distance = Centroid(), threshold = Inf, use_vanished = false)
A matcher type that matches IDs by the distance of their corresponding state space sets.
Keyword arguments
distance = Centroid()
: distance to match by, given tosetsofsets_distances
.threshold = Inf
: sets with distance larger than thethreshold
are guaranteed to not be mapped to each other.use_vanished = !isinf(threshold)
: value of the keyworduse_vanished
when used inmatch_sequentially!
.
Description
In this matcher the values compared are StateSpaceSet
s which in most cases represent attractors in the state space, but may also represent any other set such as a group of features.
Here is how this matcher works: (recall in this conversation that sets/attractors are stored in dictionaries, mapping keys/IDs to the sets, and we want to match keys in the "new" dictionary (a₊
) to those in the "old" dictionary (a₋
)).
The distance between all possible pairs of sets between the "old" sets and "new" sets is computed as a formal distance between sets. This is controlled by the distance
option, itself given to the lower-level setsofsets_distances
function, so distance
can be whatever that function accepts. That is, one of Centroid
, Hausdorff
, StrictlyMinimumDistance
, or any arbitrary user-provided function f
that given two sets f(A, B)
it returns a positive number (their distance).
Sets (in particular, their corresponding IDs) are then matched according to this distance. First, all possible ID pairs (old, new) are sorted according to the distance of their corresponding sets. The pair with smallest distance is matched. IDs 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 sets is larger than this threshold
, then it is guaranteed that the two sets will get assigned different ID in the replacement map, and hence, the set in a₊
gets the next available integer as its ID.
Attractors.SubdivisionBasedGrid
— TypeSubdivisionBasedGrid(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.RAFM
— FunctionAlias for RecurrencesFindAndMatch
.
Attractors.RecurrencesFindAndMatch
— MethodRecurrencesFindAndMatch <: GlobalContinuationAlgorithm
RecurrencesFindAndMatch(mapper::AttractorsViaRecurrences; kwargs...)
A method for global_continuation
as in Datseris2023 that is based on the recurrences algorithm for finding attractors (AttractorsViaRecurrences
) and then matching them according to their state space distance.
Keyword arguments
distance = Centroid(), threshold = Inf
: passed toMatchBySSSetDistance
.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
RecurrencesFindAndMatch
is a wrapper type. It is has been generalized by AttractorSeedContinueMatch
. It is still exported for backwards compatibility and to have a clear reference to the original algorithm developed in Datseris2023.
The source code of RecurrencesFindAndMatch
is trival: it takes the given mapper, it initializes a MatchBySSSetDistance
, and along with seeds_from_attractor
it makes the AttractorSeedContinueMatch
instance. This is the process described in Datseris2023, whereby attractors are found using the recurrences algorithm AttractorsViaRecurrences
and they are then matched by their distance in state space MatchBySSSetDistance
.
Attractors._rescale_to_01
— MethodDo "min-max" rescaling of vector of feature vectors so that its values span [0,1]
.
Attractors.additive_dict_merge!
— Methodadditive_dict_merge!(d1::Dict, d2::Dict)
Merge keys and values of d2
into d1
additively: the values of the same keys are added together in d1
and new keys are given to d1
as-is.
Attractors.aggregate_attractor_fractions
— Functionaggregate_attractor_fractions(
fractions_cont, attractors_cont, 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 a global_continuation
with RecurrencesFindAndMatch
so that similar attractors, even across parameter space, are grouped into one "attractor". Thus, the fractions of their basins are aggregated.
You could also use this function to aggregate attractors and their fractions even in a single parameter configuration, i.e., using the output of basins_fractions
.
This function is useful in cases where you want the accuracy and performance of AttractorsViaRecurrences
, but you also want the convenience of "grouping" similar attractrors like in AttractorsViaFeaturizing
for presentation or analysis purposes. For example, a high dimensional model of competition dynamics across multispecies may have extreme multistability. After finding this multistability however, one may care about aggregating all attractors into two groups: where a given species is extinct or not. This is the example highlighted in our documentation, in Extinction of a species in a multistable competition model.
Input
fractions_cont
: a vector of dictionaries mapping labels to basin fractions.attractors_cont
: a vector of dictionaries mapping labels to attractors. 1st and 2nd argument are exactly like the return values ofglobal_continuation
withRecurrencesFindAndMatch
(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 inFeaturizeGroupAcrossParameter
. The 5th argument is optional and defaults to the centroid of the features.
Return
aggregated_fractions
: same asfractions_cont
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)
.
Attractors.animate_attractors_continuation
— Functionanimate_attractors_continuation(
ds::DynamicalSystem, attractors_cont, fractions_cont, pcurve;
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 global_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_cont, fractions_cont
are the output of global_continuation
while ds, pcurve
are the input to global_continuation
.
Keyword arguments
savename = "attracont.mp4"
: name of video output file.framerate = 4
: framerate of video output.Δt, T
: propagated totrajectory
for evolving an initial condition sampled from an attractor.- Also all common plotting keywords.
figure, axis, fracaxis, legend
: named tuples propagated as keyword arguments to the creation of theFigure
, theAxis
, the "bar-like" axis containing the fractions, and theaxislegend
that adds the legend (ifadd_legend = true
).add_legend = true
: whether to display the axis legend.
Attractors.automatic_Δt_basins
— Methodautomatic_Δt_basins(ds::DynamicalSystem, grid; N = 5000) → Δt
Calculate an optimal Δt
value for basins_of_attraction
. This is done by evaluating the dynamic rule f
(vector field) at N
randomly chosen points within the bounding box of the grid. The average f
is then compared with the average diagonal length of a grid cell and their ratio provides Δt
.
Notice that Δt
should not be too small which happens typically if the grid resolution is high. It is okay if the trajectory skips a few cells. Also, Δt
that is smaller than the internal step size of the integrator will cause a performance drop.
Attractors.basin_entropy
— Methodbasin_entropy(basins::Array, ε = 20) -> Sb, Sbb
Return the basin entropy Daza2016 Sb
and basin boundary entropy Sbb
of the given basins
of attraction by considering ε
-sized boxes along each dimension.
Description
First, the n-dimensional input basins
is divided regularly into n-dimensional boxes of side ε
. If ε
is an integer, the same size is used for all dimensions, otherwise ε
can be a tuple with the same size as the dimensions of basins
. Assuming that there are $N$ ε
-boxes that cover the basins
, the basin entropy is estimated as Daza2016
\[S_b = \tfrac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{m_i}-p_{ij}\log(p_{ij})\]
where $m$ is the number of unique IDs (integers of basins
) in box $i$ and $p_{ij}$ is the relative frequency (probability) to obtain ID $j$ in the $i$ box (simply the count of IDs $j$ divided by the total in the box).
Sbb
is the boundary basin entropy. This follows the same definition as $S_b$, but now averaged over only only 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 unique IDs in basins
. In this case the boundary is intermingled: for a given initial condition we can find another initial condition that lead to another basin arbitrarily close. It provides also a simple criterion for fractality: if the boundary basin entropy Sbb
is above log(2)
then we have a fractal boundary. It doesn't mean that basins with values below cannot have a fractal boundary, for a more precise test see basins_fractal_test
. An important feature of the basin entropy is that it allows comparisons between different basins using the same box size ε
.
Attractors.basins_fractal_dimension
— Methodbasins_fractal_dimension(basins; kwargs...) -> V_ε, N_ε, d
Estimate the fractal dimension d
of the boundary between basins of attraction using a box-counting algorithm for the boxes that contain at least two different basin IDs.
Keyword arguments
range_ε = 2:maximum(size(basins))÷20
is the range of sizes of the box to test (in pixels).
Description
The output N_ε
is a vector with the number of the balls of radius ε
(in pixels) that contain at least two initial conditions that lead to different attractors. V_ε
is a vector with the corresponding size of the balls. The output d
is the estimation of the box-counting dimension of the boundary by fitting a line in the log.(N_ε)
vs log.(1/V_ε)
curve. However it is recommended to analyze the curve directly for more accuracy.
It is the implementation of the popular algorithm of the estimation of the box-counting dimension. The algorithm search for a covering the boundary with N_ε
boxes of size ε
in pixels.
Attractors.basins_fractal_test
— Methodbasins_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 ofSbb
Description
The test "looks" at the basins with a magnifier of size ε
at random. If what we see in the magnifier looks like a smooth boundary (onn average) we decide that the boundary is smooth. If it is not smooth we can say that at the scale ε
we have structures, i.e., it is fractal.
In practice the algorithm computes the boundary basin entropy Sbb
basin_entropy
for Ntotal
random boxes of radius ε
. If the computed value is equal to theoretical value of a smooth boundary (taking into account statistical errors and biases) then we decide that we have a smooth boundary. Notice that the response test_res
may depend on the chosen ball radius ε
. For larger size, we may observe structures for smooth boundary and we obtain a different answer.
The output test_res
is a symbol describing the nature of the basin and the output Sbb
is the estimated value of the boundary basin entropy with the sampling method.
Attractors.basins_fractions
— Functionbasins_fractions(basins::AbstractArray [,ids]) → 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. Optionally you may give a vector of ids
to calculate the fractions of only the chosen ids (by default ids = unique(basins)
).
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 global_continuation
.
Attractors.basins_fractions
— Methodbasins_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. ThenN
random initial conditions are chosen. Seestatespace_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 a 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
. See also convergence_and_basins_fractions
.
Keyword arguments
N = 1000
: Number of random initial conditions to generate in caseics
is a function.show_progress = true
: Display a progress bar of the process.
Attractors.basins_of_attraction
— Methodbasins_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
.
See also convergence_and_basins_of_attraction
.
Attractors.basins_of_attraction
— Methodbasins_of_attraction(mapper::AttractorsViaRecurrences; show_progress = true)
This is a special method of basins_of_attraction
that using recurrences does exactly what is described in the paper by Datseris & Wagemakers Datseris2022. By enforcing that the internal grid of mapper
is the same as the grid of initial conditions to map to attractors, the method can further utilize found exit and attraction basins, making the computation faster as the grid is processed more and more.
Attractors.bisect_to_edge
— Methodbisect_to_edge(pds::ParallelDynamicalSystem, mapper::AttractorMapper; kwargs...) -> u1, u2
Finds the basin boundary between two states u1, u2 = current_states(pds)
by bisecting along a straight line in phase space. The states u1
and u2
must belong to different basins.
Returns a triple u1, u2, success
, where u1, u2
are two new states located on either side of the basin boundary that lie less than bisect_thresh
(Euclidean distance in state space) apart from each other, and success
is a Bool indicating whether the bisection was successful (it may fail if the mapper
maps both states to the same basin of attraction, in which case a warning is raised).
Keyword arguments
bisect_thresh = 1e-7
: The maximum (Euclidean) distance between the two returned states.
Description
pds
is a ParallelDynamicalSystem
with two states. The mapper
must be an AttractorMapper
of subtype AttractorsViaProximity
or AttractorsViaRecurrences
.
If the straight line between u1
and u2
intersects the basin boundary multiple times, the method will find one of these intersection points. If more than two attractors exist, one of the two returned states may belong to a different basin than the initial conditions u1
and u2
. A warning is raised if the bisection involves a third basin.
Attractors.cluster_assignment
— MethodUtil function for classify_features
. Returns the assignment vector, in which the i-th component is the cluster index of the i-th feature.
Attractors.convergence_and_basins_fractions
— Methodconvergence_and_basins_fractions(mapper::AttractorMapper, ics::StateSpaceSet)
An extension of basins_fractions
. Return fs, labels, convergence
. The first two are as in basins_fractions
, and convergence
is a vector containing the time each initial condition took to converge to its attractor. Only usable with mappers that support id = mapper(u0)
.
See also convergence_time
.
Keyword arguments
show_progress = true
: show progress bar.
Attractors.convergence_and_basins_of_attraction
— Methodconvergence_and_basins_of_attraction(mapper::AttractorMapper, grid)
An extension of basins_of_attraction
. Return basins, attractors, convergence
, with basins, attractors
as in basins_of_attraction
, and convergence
being an array with same shape as basins
. It contains the time each initial condition took to converge to its attractor. It is useful to give to shaded_basins_heatmap
.
See also convergence_time
.
Keyword arguments
show_progress = true
: show progress bar.
Attractors.convergence_time
— Functionconvergence_time(mapper::AttractorMapper) → t
Return the approximate time 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.
Obtaining the convergence time is computationally free, so that convergence_and_basins_fractions
can always be used instead of basins_fractions
.
Attractors.crude_initial_radius
— MethodThis 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.edgetracking
— Methodedgetracking(ds::DynamicalSystem, attractors::Dict; kwargs...)
Track along a basin boundary in a dynamical system ds
with two or more attractors in order to find an edge state. Results are returned in the form of EdgeTrackingResults
, which contains the pseudo-trajectory edge
representing the track on the basin boundary, along with additional output (see below).
The system's attractors
are specified as a Dict
of StateSpaceSet
s, as in AttractorsViaProximity
or the output of extract_attractors
. By default, the algorithm is initialized from the first and second attractor in attractors
. Alternatively, the initial states can be set via keyword arguments u1
, u2
(see below). Note that the two initial states must belong to different basins of attraction.
Keyword arguments
bisect_thresh = 1e-7
: distance threshold for bisectiondiverge_thresh = 1e-6
: distance threshold for parallel integrationu1
: first initial state (defaults to first point in first entry ofattractors
)u2
: second initial state (defaults to first point in second entry ofattractors
)maxiter = 100
: maximum number of iterations before the algorithm stopsabstol = 0.0
: distance threshold for convergence of the updated edge stateT_transient = 0.0
: transient time before the algorithm starts saving the edge tracktmax = Inf
: maximum integration time of parallel trajectories until re-bisectionΔt = 0.01
: time step passed tostep!
when evolving the two trajectoriesϵ_mapper = nothing
:ϵ
parameter inAttractorsViaProximity
show_progress = true
: if true, shows progress bar and information while runningverbose = true
: if false, silences print output and warnings while runningkwargs...
: additional keyword arguments to be passed toAttractorsViaProximity
Description
The edge tracking algorithm is a numerical method to find an edge state or (possibly chaotic) saddle on the boundary between two basins of attraction. Introduced by Battelino1988 and further described by Skufca2006, the algorithm has been applied to, e.g., the laminar-turbulent boundary in plane Couette flow Schneider2008, Wada basins Wagemakers2020, as well as Melancholia states in conceptual Mehling2023 and intermediate-complexity Lucarini2017 climate models. Relying only on forward integration of the system, it works even in high-dimensional systems with complicated fractal basin boundary structures.
The algorithm consists of two main steps: bisection and tracking. First, it iteratively bisects along a straight line in state space between the intial states u1
and u2
to find the separating basin boundary. The bisection stops when the two updated states are less than bisect_thresh
(Euclidean distance in state space) apart from each other. Next, a ParallelDynamicalSystem
is initialized from these two updated states and integrated forward until the two trajectories diverge from each other by more than diverge_thresh
(Euclidean distance). The two final states of the parallel integration are then used as new states u1
and u2
for a new bisection, and so on, until a stopping criterion is fulfilled.
Two stopping criteria are implemented via the keyword arguments maxiter
and abstol
. Either the algorithm stops when the number of iterations reaches maxiter
, or when the state space position of the updated edge point changes by less than abstol
(in Euclidean distance) compared to the previous iteration. Convergence below abstol
happens after sufficient iterations if the edge state is a saddle point. However, the edge state may also be an unstable limit cycle or a chaotic saddle. In these cases, the algorithm will never actually converge to a point but (after a transient period) continue populating the set constituting the edge state by tracking along it.
A central idea behind this algorithm is that basin boundaries are typically the stable manifolds of unstable sets, namely edge states or saddles. The flow along the basin boundary will thus lead to these sets, and the iterative bisection neutralizes the unstable direction of the flow away from the basin boundary. If the system possesses multiple edge states, the algorithm will find one of them depending on where the initial bisection locates the boundary.
Output
Returns a data type EdgeTrackingResults
containing the results.
Sometimes, the AttractorMapper used in the algorithm may erroneously identify both states u1
and u2
with the same basin of attraction due to being very close to the basin boundary. If this happens, a warning is raised and EdgeTrackingResults.success = false
.
Attractors.edgetracking
— Methodedgetracking(pds::ParallelDynamicalSystem, mapper::AttractorMapper; kwargs...)
Low-level function for running the edge tracking algorithm, see edgetracking
for a description, keyword arguments and output type.
pds
is a ParallelDynamicalSystem
with two states. The mapper
must be an AttractorMapper
of subtype AttractorsViaProximity
or AttractorsViaRecurrences
.
Attractors.extract_attractors
— Methodextract_attractors(mapper::AttractorsMapper) → attractors
Return a dictionary mapping label IDs to attractors found by the mapper
. This function should be called after calling basins_fractions
with the given mapper
so that the attractors have actually been found first.
For AttractorsViaFeaturizing
, the attractors are only stored if the mapper was called with pre-defined initial conditions rather than a sampler (function returning initial conditions).
Attractors.extract_features
— Methodextract_features(mapper, ics; N = 1000, show_progress = true)
Return a vector of the features of each initial condition in ics
(as in basins_fractions
), using the configuration of mapper::AttractorsViaFeaturizing
. Keyword N
is ignored if ics isa StateSpaceSet
.
Attractors.feature_to_group
— Methodfeature_to_group(feature::AbstractVector, group_config::GroupingConfig) → group_label
Map the given feature vector to its group label (integer). This is an internal function. It is strongly recommended that feature isa SVector
.
Attractors.finite_state_machine!
— MethodMain procedure which performs one step of the finite state machine. Directly implements the algorithm of Datseris & Wagemakers 2021, see the flowchart (Figure 2). I.e., it obtains current state, updates state with respect to the cell code the trajectory it is currently on, performs the check of counters, and sets the f.s.m. to the new state, or increments counter if same state.
The basins and attractors are coded in the array with odd numbers for the basins and even numbers for the attractors. The attractor 2n
has the corresponding basin 2n+1
. This codification is changed when the basins and attractors are returned to the user. Diverging trajectories and the trajectories staying outside the grid are coded with -1.
The label 1
(initial value) outlined in the paper is 0
here instead. The function returns 0
unless the FSM has terminated its operation.
Attractors.global_continuation
— Methodglobal_continuation(gca::GlobalContinuationAlgorithm, prange, pidx, ics; kwargs...)
global_continuation(gca::GlobalContinuationAlgorithm, pcurve, ics; kwargs...)
Find and continue attractors (or representations of attractors) and the fractions of their basins of attraction across a parameter range. global_continuation
is the central function of the framework for global stability analysis illustrated in Datseris2023.
The global continuation algorithm typically references an AttractorMapper
which is used to find the attractors and basins of a dynamical system. Additional arguments that control how to continue/track/match attractors across a parameter range are given when creating gca
.
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 DynamicalSystems.set_parameter!
can be used). In contrast to traditional continuation (see online Tutorial for a comparison), global continuation can be performed over arbitrary user-defined curves in parameter space. The second call signature with pcurve
allows for this possibility. In this case pcurve
is a vector of iterables, where each itereable maps parameter indices to parameter values. These iterables can be dictionaries, named tuples, Vector{Pair}
, etc., and the sequence of the iterables defines a curve in parameter space. In fact, the version with prange, pidx
simply defines pcurve = [[pidx => p] for p in prange]
and calls the second method.
ics
are the initial conditions to use when globally sampling the state space. Like in basins_fractions
it can be either a set vector of initial conditions, or a 0-argument function that generates random initial conditions.
Possible subtypes of GlobalContinuationAlgorithm
are:
Return
fractions_cont::Vector{Dict{Int, Float64}}
. The fractions of basins of attraction.fractions_cont[i]
is a dictionary mapping attractor IDs to their basin fraction at thei
-th parameter combination.attractors_cont::Vector{Dict{Int, <:Any}}
. The continued attractors.attractors_cont[i]
is a dictionary mapping attractor ID to the attractor set at thei
-th parameter combination.
Keyword arguments
show_progress = true
: display a progress bar of the computation.samples_per_parameter = 100
: amount of initial conditions sampled at each parameter combination fromics
ifics
is a function instead of set initial conditions.
Attractors.group_features
— Methodgroup_features(features, group_config::GroupingConfig) → labels
Group the given iterable of "features" (anything that can be grouped, typically vectors of real numbers) according to the configuration and return the labels (vector of equal length as features
). See GroupingConfig
for possible grouping configuration configurations.
Attractors.haussdorff_distance
— Methodhaussdorff_distance(M1, M2) -> hd
Compute the Hausdorff distance between two binary matrices of the same size. First a distance matrix is computed using distance_matrix
for each matrix M1 and M2. The entries of this matrix are the distance in L1 metric to the closest 0 pixel in the initial matrix. The distance being 0 if it belongs to the basin with id 0.
Attractors.heatmap_basins_attractors
— Functionheatmap_basins_attractors(grid, basins, attractors; kwargs...)
Plot a heatmap of found (2-dimensional) basins
of attraction and corresponding attractors
, i.e., the output of basins_of_attraction
.
Keyword arguments
- All the common plotting keywords.
Attractors.match_sequentially!
— Methodmatch_sequentially!(dicts::Vector{Dict{Int, Any}}, matcher::IDMatcher; kw...)
Match the dicts
, a vector of dictionaries mapping IDs (integers) to values, according to the given matcher
by sequentially applying the matching_map
function to all elements of dicts
besides the first one.
In the context of Attractors.jl dicts
are typically dictionaries mapping IDs to attractors (StateSpaceSet
s), however the function is generic and would work for any values that matcher
works with.
Return rmaps
, which is a vector of dictionaries. rmaps[i]
contains the matching_map
for attractors[i+1]
, i.e., the pairs of old => new
IDs.
Keyword arguments
pcurve = nothing
: the curve of parameters along which the continuation occured, from which to extract thep, pprev
values given tomatching_map
. Seeglobal_continuation
if you are unsure what this means.ds = nothing
: propagated tomatching_map
.retract_keys::Bool = true
: Iftrue
at the end the function 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.use_vanished = false
: Ifuse_vanised = true
, then IDs (and their corresponding sets) that existed before but have vanished are kept in "memory" when it comes to matching: the current dictionary values (the attractor sets) are compared to the latest instance of all values that have ever existed, each with a unique ID, and get matched to their closest ones. The value of this keyword is obtained from thematcher
.
Attractors.match_sequentially!
— Methodmatch_sequentially!(continuation_quantity::Vector{Dict}, rmaps::Vector{Dict})
Do the same as in match_sequentially!
above, now given the vector of matching maps, and for any arbitrary quantity that has been tracked in the globalcontinuation. `continuationquantitycan for example be
fractionscontfrom [
globalcontinuation`](@ref).
Attractors.matching_map!
— Methodmatching_map!(a₊, a₋, matcher; kw...) → rmap
Convenience function that first calls matching_map
and then replaces the IDs in a₊
with this rmap
.
Attractors.matching_map
— Methodmatching_map(b₊::AbstractArray, b₋::AbstractArray, matcher::MatchByBasinOverlap)
Special case of matching_map
where instead of having as input dictionaries mapping IDs to values, we have Array
s which represent basins of attraction and whose elements are the IDs.
See MatchByBasinOverlap
for how matching works.
Attractors.matching_map
— Methodmatching_map(
a₊::Dict, a₋::Dict, matcher;
ds::DynamicalSystem, p, pprev, next_id
) → rmap
Given dictionaries a₊, a₋
mapping IDs to values, return a replacement map: a dictionary mapping the IDs (keys) in dictionary a₊
to IDs (keys) in dictionary a₋
, so that so that values in a₊
that are the "closest" to values in a₋
get assigned the same key as in a₋
. In this way keys of a₊
are "matched" to keys of a₋
. Use swap_dict_keys
to apply rmap
to a₊
or to other dictionaries with same keys as a₊
.
How matching happens, i.e., how "closeness" is defined, depends on the algorithm matcher
.
The values contained in a₊, a₋
can be anything supported by matcher
. Within Attractors.jl they are typically StateSpaceSet
s representing attractors. Typically the +,- mean after and before some change of parameter of a dynamical system.
Keyword arguments
ds
: the dynamical system that generateda₊, a₋
.p, pprev
: the parameters corresponding toa₊, a₋
. Both need to be iterables mapping parameter index to parameter value (such asDict, Vector{Pair}
, etc., so whatever can be given as input toDynamicalSystems.set_parameters!
).next_id = next_free_id(a₊, a₋)
: the ID to give to values ofa₊
that cannot be matched toa₋
and hence must obtain a new unique ID.
Some matchers like MatchBySSSetDistance
do not utilize ds, p, pprev
in any way while other matchers like MatchByBasinEnclosure
do, and those require expliticly giving values to ds, p, pprev
as their default values is just nothing
.
Attractors.mfs_brute_force
— MethodThis function works on the results obtained by crude_initial_radius
. It starts from the best shock found so far and tries to find a better one by continuously reducing the radius of the sphere on the surface of which it generates random perturbations. If perturbation with the same basin of attraction is found, it updates the best shock found so far and reduces the radius of the sphere. It repeats this process total_iterations times and returns the best perturbation found.
Attractors.minimal_fatal_shock
— Methodminimal_fatal_shock(mapper::AttractorMapper, u0, search_area, algorithm; kw...)
Return the minimal fatal shock mfs
(also known as excitability threshold) 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
. For 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
).
An alias to minimal_fata_shock
is excitability_threshold
.
Keyword arguments
metric = LinearAlgebra.norm
: a metric function that gives the norm of a perturbation vector. This keyword is ignored for theMFSBruteForce
algorithm.target_id = nothing
: when notnothing
, it should be an integer or a vector of integers corresponding to target attractor label(s). Then, the MFS is estimated based only on perturbations that lead to the target attractor(s).
Description
The minimal fatal shock is defined as the 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.
The excitability threshold is a concept nearly identical, however, instead of looking for a perturbation that simply brings us out of the basin, we look for the smallest perturbation that brings us into specified basin(s). This is enabled via the keyword target_id
.
Attractors.next_free_id
— Methodnext_free_id(new::Dict, old::Dict)
Return the minimum key of the "new" dictionary that doesn't exist in the "old" dictionary.
Attractors.optimal_radius_dbscan_knee
— MethodFind the optimal radius ϵ of a point neighborhood for use in DBSCAN through the elbow method (knee method, highest derivative method).
Attractors.optimal_radius_dbscan_silhouette
— MethodFind the optimal radius ε of a point neighborhood to use in DBSCAN, the unsupervised clustering method for AttractorsViaFeaturizing
. Iteratively search for the radius that leads to the best clustering, as characterized by quantifiers known as silhouettes. Does a linear (sequential) search.
Attractors.optimal_radius_dbscan_silhouette_optim
— MethodSame as optimal_radius_dbscan_silhouette
, but find minimum via optimization with Optim.jl.
Attractors.optimal_radius_dbscan_silhouette_original
— MethodFind the optimal radius ε of a point neighborhood to use in DBSCAN, the unsupervised clustering method for AttractorsViaFeaturizing
. The basic idea is to iteratively search for the radius that leads to the best clustering, as characterized by quantifiers known as silhouettes.
Attractors.overwrite_dict!
— Methodoverwrite_dict!(old::Dict, new::Dict)
In-place overwrite the old
dictionary for the key-value pairs of the new
.
Attractors.plot_attractors
— Functionplot_attractors(attractors::Dict{Int, StateSpaceSet}; kwargs...)
Plot the attractors as a scatter plot.
Keyword arguments
- All the common plotting keywords. Particularly important is the
access
keyword. sckwargs = (strokewidth = 0.5, strokecolor = :black,)
: additional keywords propagated to theMakie.scatter
function that plots the attractors.
Attractors.plot_attractors_curves
— Functionplot_attractors_curves(attractors_cont, attractor_to_real [, prange]; kw...)
Same as in plot_basins_curves
but visualize the attractor dependence on the parameter(s) instead of their basin 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. See also plot_basins_attractors_curves
.
Same keywords as plot_basins_curves
. See also plot_continuation_curves
.
Attractors.plot_basins_attractors_curves
— Functionplot_basins_attractors_curves(
fractions_cont, attractors_cont, a2rs [, prange]
kwargs...
)
Convenience combination of plot_basins_curves
and plot_attractors_curves
in a multi-panel plot that shares legend, colors, markers, etc. This function allows a2rs
to be a Vector
of functions, each mapping attractors into real numbers. Below the basins fractions plot, one additional panel is created for each entry in a2rs
. a2rs
can also be a single function, in which case only one panel is made.
Attractors.plot_basins_curves
— Functionplot_basins_curves(fractions_cont [, prange]; kw...)
Plot the fractions of basins of attraction versus a parameter range/curve, i.e., visualize the output of global_continuation
. See also plot_basins_attractors_curves
and plot_continuation_curves
.
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 fractionseparatorwidth = 1, separatorcolor = "white"
: adds a line separating the fractions if the style is:band
axislegend_kwargs = (position = :lt,)
: propagated toaxislegend
if a legend is addedseries_kwargs = NamedTuple()
: propagated to the band or scatterline plot- Also all common plotting keywords.
Attractors.plot_continuation_curves
— Functionplot_continuation_curves(continuation_info [, prange]; kwargs...)
Same as in plot_basins_curves
but visualize any arbitrary quantity characterizing the continuation. Hence, the continuation_info
is of exactly the same format as fractions_cont
: a vector of dictionaries, each dictionary mapping attractor IDs to real numbers. continuation_info
is meant to accompany attractor_info
in plot_attractors_curves
. To produce continuation_info
from attractor_info
you can do something like:
continuation_info = map(attractors_cont) do attractors
Dict(k => f(A) for (k, A) in attractors)
end
with f
your function of interest that returns a real number.
Attractors.recurrences_map_to_label!
— Methodrecurrences_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.reset_mapper!
— Methodreset_mapper!(mapper::AttractorsMapper)
Reset all accumulated information of mapper
so that it is as if it has just been initialized. Useful in for
loops that loop over a parameter of the dynamical system stored in mapper
.
Attractors.retract_keys_to_consecutive
— Methodretract_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 global_continuation
it skips the special key -1
.
Attractors.shaded_basins_heatmap
— Functionshaded_basins_heatmap(grid, basins, attractors, iterations; kwargs...)
Plot a heatmap of found (2-dimensional) basins
of attraction and corresponding attractors
. A matrix iterations
with the same size of basins
must be provided to shade the color according to the value of this matrix. A 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 condition to converge. See also convergence_time
to store this iteration number.
Keyword arguments
show_attractors = true
: shows the attractor on plotmaxit = maximum(iterations)
: clip the values ofiterations
to
the value maxit
. Useful when there are some very long iterations and keep the range constrained to a given interval.
- All the common plotting keywords.
Attractors.silhouettes_new
— MethodCalculate silhouettes. A bit slower than the implementation in Clustering.jl
but seems to be more robust. The latter seems to be incorrect in some cases.
Attractors.subdivision_based_grid
— Methodsubdivision_based_grid(ds::DynamicalSystem, grid; maxlevel = 4, q = 0.99)
Construct a grid structure SubdivisionBasedGrid
that can be directly passed as a grid to AttractorsViaRecurrences
. The input grid
is an originally coarse grid (a tuple of AbstractRange
s). The state space speed is evaluate in all cells of the grid
. Cells with small speed (when compared to the "max" speed) resultin in this cell being subdivided more. To avoid problems with spikes in the speed, the q
-th quantile of the velocities is used as the "max" speed (use q = 1
for true maximum). The subdivisions in the resulting grid are clamped to at most value maxlevel
.
This approach is designed for continuous time systems in which different areas of the state space flow may have significantly different velocity. In case of originally coarse grids, this may lead AttractorsViaRecurrences
being stuck in some state space regions with a small motion speed and false identification of attractors.
Attractors.swap_dict_keys!
— Methodswap_dict_keys!(d::Dict, matching_map::Dict)
Swap the keys of a dictionary d
given a matching_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.test_wada_merge
— Methodtest_wada_merge(basins, r) -> p
Test if the 2D array basins
has the Wada property using the merging technique of Daza2018.
Description
The technique consists in computing the generalized basins of each attractor. These new basins are formed with on of the basins and the union of the other basins. A new boundary is defined by these two objects. The algorithm then computes the distance between each boundaries of these basins pairwise. If all the boundaries are within some distance r
, there is a unique boundary separating the basins and we have the wada property. The algorithm returns the maximum proportion of pixels of a boundary with distance strictly greater than r
from another boundary.
If p == 0
, we have the Wada property for this value of r
. If p > 0
, the criteria to decide if the basins are Wada is left to the user. Numerical inaccuracies may be responsible for a small percentage of points with distance larger than r
Attractors.tipping_probabilities
— Methodtipping_probabilities(basins_before, basins_after) → P
Return the tipping probabilities of the computed basins before and after a change in the system parameters (or time forcing), according to the definition of Kaszas2019.
The input basins
are integer-valued arrays, where the integers enumerate the attractor, e.g. the output of basins_of_attraction
.
Description
Let $\mathcal{B}_i(p)$ denote the basin of attraction of attractor $A_i$ at parameter(s) $p$. Kaszás et al Kaszas2019 define the tipping probability from $A_i$ to $A_j$, given a parameter change in the system of $p_- \to p_+$, as
\[P(A_i \to A_j | p_- \to p_+) = \frac{|\mathcal{B}_j(p_+) \cap \mathcal{B}_i(p_-)|}{|\mathcal{B}_i(p_-)|}\]
where $|\cdot|$ is simply the volume of the enclosed set. The value of $P(A_i \to A_j | p_- \to p_+)$ is P[i, j]
. The equation describes something quite simple: what is the overlap of the basin of attraction of $A_i$ at $p_-$ with that of the attractor $A_j$ at $p_+$. If basins_before, basins_after
contain values of -1
, corresponding to trajectories that diverge, this is considered as the last attractor of the system in P
.
Attractors.uncertainty_exponent
— Methoduncertainty_exponent(basins; kwargs...) -> ε, N_ε, α
Estimate the uncertainty exponentGrebogi1983 of the basins of attraction. This exponent is related to the final state sensitivity of the trajectories in the phase space. An exponent close to 1
means basins with smooth boundaries whereas an exponent close to 0
represent completely fractalized basins, also called riddled basins.
The output N_ε
is a vector with the number of the balls of radius ε
(in pixels) that contain at least two initial conditions that lead to different attractors. The output α
is the estimation of the uncertainty exponent using the box-counting dimension of the boundary by fitting a line in the log.(N_ε)
vs log.(1/ε)
curve. However it is recommended to analyze the curve directly for more accuracy.
Keyword arguments
range_ε = 2:maximum(size(basins))÷20
is the range of sizes of the ball to test (in pixels).
Description
A phase space with a fractal boundary may cause a uncertainty on the final state of the dynamical system for a given initial condition. A measure of this final state sensitivity is the uncertainty exponent. The algorithm probes the basin of attraction with balls of size ε
at random. If there are a least two initial conditions that lead to different attractors, a ball is tagged "uncertain". f_ε
is the fraction of "uncertain balls" to the total number of tries in the basin. In analogy to the fractal dimension, there is a scaling law between, f_ε ~ ε^α
. The number that characterizes this scaling is called the uncertainty exponent α
.
Notice that the uncertainty exponent and the box counting dimension of the boundary are related. We have Δ₀ = D - α
where Δ₀
is the box counting dimension computed with basins_fractal_dimension
and D
is the dimension of the phase space. The algorithm first estimates the box counting dimension of the boundary and returns the uncertainty exponent.
Attractors.unique_keys
— Methodunique_keys(v::Iterator{<:AbstractDict})
Given a vector of dictionaries, return a sorted vector of the unique keys that are present across all dictionaries.