Basin maps

BasinMap(ds::DynamicalSystem, args...; kwargs...) → bmap

Subtypes of BasinMap are structures that map initial conditions of ds to unique IDs. These IDs (typically) correspond to the basins of attraction and corresponding attractors of ds. The found attractors are stored inside the bmap, and can be obtained by calling attractors = extract_attractors(bmap) when applicable.

Currently available mapping methods:

• BasinMapProximity
• BasinMapRecurrences
• BasinMapFeaturizeGroup

All BasinMap subtypes can be used with basins_fractions or basins_of_attraction.

In addition, most basin maps can be called as a function of an initial condition:

label = bmap(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:

• BasinMapProximity
• BasinMapRecurrences
• BasinMapFeaturizeGroup with the GroupViaHistogram or GroupViaNearestFeature configurations.

See also StabilityQuantifiersAccumulator that extends this interface to accelerate estimation of stability quantifiers.

For developers

BasinMap defines an extendable interface. A new type needs to subtype BasinMap and implement the following:

• extract_attractors
• the internal function Attractors.referenced_dynamical_system(bmap).
• the internal function Attractors.can_map_individual_ic(bmap)::Bool.
• and if possible, id = bmap(u0)

From these, everything else in the rest of the library just works!

If it is not possible to implement id = bmap(u0), then instead extend the function basins_fractions_grouped(bmap, ics, progress, labels), where ics is always a Vector of initial conditions, progress is a preinitialized progress bar, and labels is a preinitialized container of labels. If !isempty(labels), then its full length must be filled with the ids corresponding to the first N entries of ics.

extract_attractors(BoA::BasinsOfAttraction) → attractors

Return the attractors component of a BasinsOfAttraction object, a dictionary mapping attractor labels to attractors represented as StateSpaceSet's.

extract_attractors(bmap::AttractorsMapper) → attractors

Return a dictionary mapping label IDs to attractors found by the bmap. This function should be called after mapping initial conditions with bmap (e.g., calling basins_fractions) so that the attractors have actually been found first.

For developing a new bmap: extend the internal function _extract_attractors.

convergence_time(bmap::BasinMap) → t

Return the approximate time the bmap took to converge to an attractor. This function should be called just right after bmap(u0) was called with u0 the initial condition of interest. Hence it is only valid with BasinMap 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.

BasinMapRecurrences(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.
2. 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).
3. 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_recurrences. 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 stop_at_Δt = true.
• stop_at_Δt = false: Only used if the input dynamical system is CoupledODEs. It is given as a third input to step!. Value true 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. stop_at_Δt = true should also be used if you want an accurate estimate of convergence_time.

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 BasinMapRecurrences 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 documentation, under the recurrences animation page.

automatic_Δt_recurrences(ds::DynamicalSystem, grid; kw...) → Δt

Calculate an optimal Δt value for BasinMapRecurrences. This is done by evaluating the dynamic rule f (vector field) at N randomly chosen points within the bounding box of the grid. The aggregated f is then compared with the average diagonal length of a grid cell and their ratio provides Δt, multiplied by factor.

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.

Keyword arguments

• N = 5000: random samples in the grid to use.
• agg = mean: function used to average f. median is another option, and another one is v -> quantile(v, 0.1) for some low quantile.
• factor = 10: factor multiplying the ratio of diagonal/speed.
• norm = LinearAlgebra.norm: function used to obtain speed from state space velocity.

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 BasinMapRecurrences. 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 BasinMapRecurrences being stuck in some state space regions with a small motion speed and false identification of attractors.

BasinMapProximity(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, which have to be in a form of a dictionary mapping attractor labels to StateSpaceSets containing the attractors.

Keywords

• ε = nothing: Distance below which a trajectory has converged to an attractor, see below. Type \varepsilon<TAB> to input ε.
• Ttr = 0: Transient time to first evolve the system for before checking for proximity.
• Δt = 1: Step time given to step!.
• stop_at_Δt = false: Third argument 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 = 10000: If the ds has been stepped this many times without coming ε-near to any attractor, it is assumed that the trajectory diverged (gets labelled as -1).
• distance = StrictlyMinimumDistance(): Distance function for evaluating the distance between the trajectory end-point and the given attractors. Can be anything given to set_distance.

Description

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

attractors do not have to be "true" attractors. Any arbitrary sets in the state space can be provided.

If an ε::Real is not provided by the user, a value is computed automatically as 1/10th 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.

The convergence_time is Inf if an initial condition has not converged. As such, the convergence time is always a float type even for discrete time systems.

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

Initialize a bmap 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 BasinMap for how to use the bmap. This bmap also allows the syntax bmap(u0) if the grouping_config is is notGroupViaClustering or GroupViaPairwiseComparison.

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 DynamicalSystems.trajectory for integrating an initial condition to yield A, t.
• container: Also propagated to DynamicalSystems.trajectory, but here its default value is Vector if dimension(ds) ≥ 32.
• 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 BasinMapRecurrences.

Attractors are stored and can be accessed with extract_attractors, however it should be clear that this bmap 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.

Note that the convergence time of this bmap is always the constant T + Ttr.

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

GroupingConfig

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

Currently available grouping configurations are:

• GroupViaClustering
• GroupViaNearestFeature
• GroupViaHistogram
• GroupViaPairwiseComparison

For developers

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 features to group index, then instead extend the internal functionfeature_to_group(feature, config). This will also allow doing id = bmap(u0) with BasinMapFeaturizeGroup.
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 BasinMapFeaturizeGroup. 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 BasinMapFeaturizeGroup. 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...) <: GroupingConfig

Initialize a struct that contains instructions on how to group features in group_features (or BasinMapFeaturizeGroup). GroupViaNearestFeature accepts a template, which is a dictionary mapping unique labels to unique feature vectors. Then, feature vectors are labelled according to the feature vector in templates that is closest (the label is the key of the closest template).

templates can also be just a vector of feature vectors.

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.
• use_svector = true: Convert template vectors to SVector.
• use_kdtree = true: If true, use a KDTree for nearest-neighbor lookup. If false, use brute-force search over all template vectors.

GroupViaPairwiseComparison <: GroupingConfig
GroupViaPairwiseComparison(; threshold::Real, metric...)

Initialize a struct that contains instructions on how to group features in BasinMapFeaturizeGroup. 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 = 0.1: A real number defining the maximum distance two features can have to be considered in the same cluster - above the threshold, features are different. This value simply needs to be large enough to differentiate clusters. A good value for threshold depends on the feature variability within a cluster the chosen metric, and whether features are rescaled. See description below for more.
• 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 for the default metric. threshold, however, it should be avoid for when there is only one attractor for the system, because it leads to it being wrongly classified as many different attractors at the same location).

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 scales as O(N) in memory and performance with N the number of features. This is a huge difference versus the O(N^2) of GroupViaClustering.

extract_features(bmap::BasinMapFeaturizeGroup, 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 bmap. Keyword N is ignored if ics isa StateSpaceSet.