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

Map initial conditions of ds to attractors by identifying attractors on the fly based on recurrences in the state space, as outlined in (Datseris and Wagemakers, 2022). 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:

1. A tuple of sorted AbstractRanges for a regular grid.

Example is grid = (xg, yg) where xg = yg = range(-5, 5; length = 100) for a two-dimensional system.

1. A tuple of sorted AbstractVectors 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).

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

• 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. Once attractor_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. If true, each visited cell will only store a point once, which is desirable for fixed points and limit cycles. If false then attractor_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 if sparse = 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 for attractor_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 when sparse = 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 > than horizon_limit: the initial condition is labelled -1.
• If none of the above happens, the initial condition is labelled -1 after maximum_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 Examples page.

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

SubdivisionBasedGrid(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.

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

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

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

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 (Stender and Hoffmann, 2021) and MCBB (Gelbrecht et al., 2020). 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 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
• GroupViaPairwiseComparison

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 (Stender and Hoffmann, 2021) and MCBB (Gelbrecht et al., 2020). 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 real-valued 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 (Stender and Hoffmann, 2021) 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 (Ester et al., 1996) and (Schubert et al., 2017). 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.

GroupViaPairwiseComparison(; 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 function metric(a, b) that returns the distance between two features a and b, outputs of featurizer. Any Metric 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.

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.