Finding Attractors

AttractorMapper(ds::DynamicalSystem, args...; kwargs...) → mapper

Subtypes of AttractorMapper are structures that map initial conditions of ds to attractors. Currently available mapping methods:

• AttractorsViaProximity
• AttractorsViaRecurrences
• AttractorsViaFeaturizing

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 the GroupViaHistogram configuration.

AttractorsViaRecurrences(ds::DynamicalSystem, grid::Tuple; kwargs...)

Map initial conditions of ds to attractors by identifying attractors on the fly based on recurrences in the state space, as outlined by Datseris & Wagemakers^{[Datseris2022]}.

grid is a tuple of ranges partitioning the state space so that a finite state machine can operate on top of it. For example grid = (xg, yg) where xg = yg = range(-5, 5; length = 100) for a two-dimensional system. 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 interval representation of the state space grid. If true, uses a sparse array, whose memory usage is in general more efficient than a regular array obtained with sparse=false. In practice, the sparse representation should always be preferred when searching for basins_fractions. Only for very low dimensional systems and for computing the full basins_of_attraction the non-sparse version should be used.

Time evolution configuration

• Ttr = 0: Skip a transient before the recurrence routine begins.
• Δt: Approximate integration time step (second argument of the step! function). The keyword Dt can also be used instead if Δ (\Delta) is not accessible. It is 1 for discrete time systems. For continuous systems, an automatic value is calculated using automatic_Δt_basins. For very fine grids, this can become very small, much smaller than the typical integrator internal step size in case of adaptive integrators. In such cases, use force_non_adaptive = true.
• force_non_adaptive = false: Only used if the input dynamical system is CoupledODEs. If true the additional keywords adaptive = false, dt = Δt are given as diffeq to the CoupledODEs. This means that adaptive integration is turned off and Δt is used as the ODE integrator timestep. This is useful in (1) very fine grids, and (2) if some of the attractors are limit cycles. We have noticed that in this case the integrator timestep becomes commensurate with the limit cycle period, leading to incorrectly counting the limit cycle as more than one attractor.

Finite state machine configuration

• mx_chk_att = 2: Μaximum checks of consecutives hits of an existing attractor cell before declaring convergence to that existing attractor.
• mx_chk_hit_bas = 10: Maximum check of consecutive visits of the same basin of attraction before declaring convergence to an existing attractor.
• mx_chk_fnd_att = 100: Maximum check of consecutive visits to a previously visited unlabeled cell before declaring we have found a new attractor.
• mx_chk_loc_att = 100: Maximum check of consecutive visits to cells marked as a new attractor, during the attractor identification phase, before declaring we that we have identified the new attractor with sufficient cells.
• 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, after an attractor is found. If true, each visited cell will only store a point once, which is desirable for fixed points and limit cycles. If false, at least mx_chk_loc_att points are stored per attractor, leading to more densely stored attractors, which may be desirable for instance in chaotic attractors.
• mx_chk_lost = 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.
• mx_chk_safety = 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 haulting for innappropriate grids, where a found attractor may intersect in the same cell with a new attractor the orbit traces (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. A recurrence in the state space 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 the given grid. The FSM decides when an initial condition has successfully converged into an attractor. An array, internally called "basins", stores the state of the FSM on the grid, according to the indexing system described in ^{[Datseris2022]}. As the system is integrated more and more, the information of the "basins" becomes richer and richer with more identified attractors or with grid cells that belong to basins of already found attractors. Notice that only in the special method basins_of_attraction(mapper::AttractorsViaRecurrences) the information of the attraction or exit basins is utilized. In other functions like basins_fractions only the attractor locations are utilized, as the basins themselves are not stored.

The iteration of a given initial condition continues until one of the following happens:

• The trajectory hits mx_chk_fnd_att times in a row grid cells previously visited: it is considered that an attractor is found and is labelled with a new ID. Then, iteration continues a bit more until we have identified the attractor with sufficient accuracy, i.e., until mx_chk_loc_att cells with the new ID have been visited.
• The trajectory hits an already identified attractor mx_chk_att consecutive times: the initial condition is numbered with the attractor's ID.
• The trajectory hits a known basin mx_chk_hit_bas times in a row: the initial condition belongs to that basin and is numbered accordingly. Notice that basins are stored and used only when sparse = false.
• The trajectory spends mx_chk_lost steps outside the defined grid or the norm of the integrator state becomes > than horizon_limit: the initial condition's label is set to -1.
• If none of the above happens, the initial condition is labelled -1 after and mx_chk_safety integrator steps.

automatic_Δt_basins(ds::DynamicalSystem, grid; N = 5000) → Δt

Calculate an optimal Δt value for basins_of_attraction. This is done by evaluating the dynamic rule f (vector field) at N randomly chosen points of the grid. The average f is then compared with the 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 for basins_of_attraction if the trajectory skips a few cells. But if Δt is too small the default values for all other keywords such as mx_chk_hit_bas need to be increased drastically.

Also, Δt that is smaller than the internal step size of the integrator will cause a performance drop.

AttractorsViaProximity(ds::DynamicalSystem, attractors::Dict [, ε]; kwargs...)

Map initial conditions to attractors based on whether the trajectory reaches ε-distance close to any of the user-provided attractors. They have to be in a form of a dictionary mapping attractor labels to StateSpaceSets containing the attractors.

The system gets stepped, and at each step the minimum distance to all attractors is computed. If any of these distances is < ε, then the label of the nearest attractor is returned.

If an ε::Real is not provided by the user, a value is computed automatically as half of the minimum distance between all attractors. This operation can be expensive for large attractor StateSpaceSets. If length(attractors) == 1, then ε becomes 1/10 of the diagonal of the box containing the attractor. If length(attractors) == 1 and the attractor is a single point, an error is thrown.

Because in this method the attractors are already known to the user, the method can also be called supervised.

Keywords

• Ttr = 100: Transient time to first evolve the system for before checking for proximity.
• Δt = 1: Step time given to step!.
• horizon_limit = 1e3: If the maximum distance of the trajectory from any of the given attractors exceeds this limit, it is assumed that the trajectory diverged (gets labelled as -1).
• mx_chk_lost = 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).

AttractorsViaFeaturizing(
    ds::DynamicalSystem, featurizer::Function,
    grouping_config = GroupViaClustering(); kwargs...
)

Initialize a mapper that maps initial conditions to attractors using a featurizing and grouping approach. This is a supercase of the featurizing and clustering approach that is utilized by bSTAB^{[Stender2021]} and MCBB^{[Gelbrecht2021]}. See AttractorMapper for how to use the mapper. This mapper also allows the syntax mapper(u0) but only if the grouping_config is not GroupViaClustering.

featurizer is a function f(A, t) that takes as an input an integrated trajectoryA::StateSpaceSetand the corresponding time vectortand returns a vectorvof features describing the trajectory. For better performance, it is strongly recommended thatv isa SVector{<:Real}`.

grouping_config is an instance of any subtype of GroupingConfig and decides how features will be grouped into attractors, see below.

See also the intermediate functions extract_features and group_features, which can be utilized when wanting to work directly with features.

Keyword arguments

• T=100, Ttr=100, Δt=1: Propagated to trajectory.
• threaded = true: Whether to run the generation of features over threads by integrating trajectories in parallel.

Description

The trajectory X of an initial condition is transformed into features. Each feature is a number useful in characterizing the attractor the initial condition ends up at, and distinguishing it from other attractors. Example features are the mean or standard deviation of some the dimensions of the trajectory, the entropy of some of the dimensions, the fractal dimension of X, or anything else you may fancy.

All feature vectors (each initial condition = 1 vector) are then grouped using one of the sevaral available grouping configurations. Each group is assumed to be a unique attractor, and hence each initial condition is labelled according to the group it is part of. The method thus relies on the user having at least some basic idea about what attractors to expect in order to pick the right features, and the right way to group them, in contrast to AttractorsViaRecurrences.

GroupingConfig

Supertype for configuration structs on how to group features together. Used in several occasions such as AttractorsViaFeaturizing or aggregate_attractor_fractions.

Currently available grouping configurations are:

• GroupViaClustering
• GroupViaNearestFeature
• GroupViaHistogram

GroupingConfig defines an extendable interface. The only thing necessary for a new grouping configuration is to:

1. Make a new type and subtype GroupingConfig.
2. If the grouping allows for mapping individual initial conditions to IDs, then instead extend the internal function feature_to_group(feature, config). This will allow doing id = mapper(u0) with AttractorsViaFeaturizing.
3. Else, extend the function group_features(features, config). You could still extend group_features even if (2.) is satisfied, if there are any performance benefits.
4. Include the new grouping file in the grouping/all_grouping_configs.jl and list it in this documentation string.

GroupViaClustering(; kwargs...)

Initialize a struct that contains instructions on how to group features in AttractorsViaFeaturizing. GroupViaClustering clusters features into groups using DBSCAN, similar to the original work by bSTAB^{[Stender2021]} and MCBB^{[Gelbrecht2021]}. Several options on clustering are available, see keywords below.

The defaults are a significant improvement over existing literature, see Description.

Keyword arguments

• clust_distance_metric = Euclidean(): A metric to be used in the clustering. It can be any function f(a, b) that returns the distance between vectors a, b. All metrics from Distances.jl can be used here.
• rescale_features = true: if true, rescale each dimension of the extracted features separately into the range [0,1]. This typically leads to more accurate clustering.
• min_neighbors = 10: minimum number of neighbors (i.e. of similar features) each feature needs to have, including counting its own self, in order to be considered in a cluster (fewer than this, it is labeled as an outlier, -1).
• use_mmap = false: whether to use an on-disk map for creating the distance matrix of the features. Useful when the features are so many where a matrix with side their length would not fit to memory.

Keywords for optimal radius estimation

• optimal_radius_method::Union{Real, String} = "silhouettes_optim": if a real number, it is the radius used to cluster features. Otherwise, it determines the method used to automatically determine that radius. Possible values are:
  • "silhouettes": Performs a linear (sequential) search for the radius that maximizes a statistic of the silhouette values of clusters (typically the mean). This can be chosen with silhouette_statistic. The linear search may take some time to finish. To increase speed, the number of radii iterated through can be reduced by decreasing num_attempts_radius (see its entry below).
  • "silhouettes_optim": Same as "silhouettes" but performs an optimized search via Optim.jl. It's faster than "silhouettes", with typically the same accuracy (the search here is not guaranteed to always find the global maximum, though it typically gets close).
  • "knee": chooses the the radius according to the knee (a.k.a. elbow, highest-derivative method) and is quicker, though generally leading to much worse clustering. It requires that min_neighbors > 1.
• num_attempts_radius = 100: number of radii that the optimal_radius_method will try out in its iterative procedure. Higher values increase the accuracy of clustering, though not necessarily much, while always reducing speed.
• silhouette_statistic::Function = mean: statistic (e.g. mean or minimum) of the silhouettes that is maximized in the "optimal" clustering. The original implementation in ^{[Stender2021]} used the minimum of the silhouettes, and typically performs less accurately than the mean.
• max_used_features = 0: if not 0, it should be an Int denoting the max amount of features to be used when finding the optimal radius. Useful when clustering a very large number of features (e.g., high accuracy estimation of fractions of basins of attraction).

Description

The DBSCAN clustering algorithm is used to automatically identify clusters of similar features. Each feature vector is a point in a feature space. Each cluster then basically groups points that are closely packed together. Closely packed means that the points have at least min_neighbors inside a ball of radius optimal_radius centered on them. This method typically works well if the radius is chosen well, which is not necessarily an easy task. Currently, three methods are implemented to automatically estimate an "optimal" radius.

Estimating the optimal radius

The default method is the silhouettes method, which includes keywords silhouette and silhouette_optim. Both of them search for the radius that optimizes the clustering, meaning the one that maximizes a statistic silhouette_statistic (e.g. mean value) of a quantifier for the quality of each cluster. This quantifier is the silhouette value of each identified cluster. A silhouette value measures how similar a point is to the cluster it currently belongs to, compared to the other clusters, and ranges from -1 (worst matching) to +1 (ideal matching). If only one cluster is found, the assigned silhouette is zero. So for each attempted radius in the search the clusters are computed, their silhouettes calculated, and the statistic of these silhouettes computed. The algorithm then finds the radius that leads to the maximum such statistic. For optimal_radius_method = "silhouettes", the search is done linearly, from a minimum to a maximum candidate radius for optimal_radius_method = "silhouettes"; optimal_radius_method = silhouettes_optim, it is done via an optimized search performed by Optim.jl which is typically faster and with similar accuracy. A third alternative is the"elbow" method, which works by calculating the distance of each point to its k-nearest-neighbors (with k=min_neighbors) and finding the distance corresponding to the highest derivative in the curve of the distances, sorted in ascending order. This distance is chosen as the optimal radius. It is described in ^{[Kriegel1996]} and ^{[Schubert2017]}. It typically performs considerably worse than the "silhouette" methods.

GroupViaHistogram(binning::FixedRectangularBinning)

Initialize a struct that contains instructions on how to group features in AttractorsViaFeaturizing. GroupViaHistogram performs a histogram in feature space. Then, all features that are in the same histogram bin get the same label. The binning is an instance of FixedRectangularBinning from ComplexityMeasures.jl. (the reason to not allow RectangularBinning is because during continuation we need to ensure that bins remain identical).

GroupViaNearestFeature(templates; kwargs...)

Initialize a struct that contains instructions on how to group features in AttractorsViaFeaturizing. GroupViaNearestFeature accepts a template, which is a vector of features. Then, generated features from initial conditions in AttractorsViaFeaturizing are labelled according to the feature in templates that is closest (the label is the index of the closest template).

templates can be a vector or dictionary mapping keys to templates. Internally all templates are converted to SVector for performance. Hence, it is strongly recommended that both templates and the output of the featurizer function in AttractorsViaFeaturizing return SVector types.

Keyword arguments

• metric = Euclidean(): metric to be used to quantify distances in the feature space.
• max_distance = Inf: Maximum allowed distance between a feature and its nearest template for it to be assigned to that template. By default, Inf guarantees that a feature is assigned to its nearest template regardless of the distance. Features that exceed max_distance to their nearest template get labelled -1.

group_features(features, group_config::GroupingConfig) → labels

Group the given vector of feature vectors according to the configuration and return the labels (vector of equal length as features). See AttractorsViaFeaturizing for possible configurations.

extract_features(mapper, ics; N = 1000, show_progress = true)

Return a vector of the features of each initial condition in ics (as in basins_fractions), using the configuration of mapper::AttractorsViaFeaturizing. Keyword N is ignored if ics isa StateSpaceSet.