API

AttractorMapper(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:

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

For developers

AttractorMapper defines an extendable interface. A new type needs to subtype AttractorMapper and implement extract_attractors, id = mapper(u0) and the internal function Attractors.referenced_dynamical_system(mapper). 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) with ics a vector of initial conditions.

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

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 documentation, under the recurrences animation 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 distance of the current state to all attractors is computed. If any of these distances is < ε, then the label of the nearest attractor is returned. The distance is defined by the distance keyword.

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

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

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

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 AttractorsViaFeaturizing or aggregate_attractor_fractions.

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 = 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 <: GroupingConfig
GroupViaPairwiseComparison(; threshold::Real, metric...)

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

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

basins_fractions(
    mapper::AttractorMapper,
    ics::Union{AbstractVector, 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:

• an AbstractVector of initial conditions, in which case all are used. Typically this is a StateSpaceSet.
• a 0-argument function ics() that spits out random initial conditions. Then N random initial conditions are chosen. See statespace_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 case ics is a function.
• show_progress = true: Display a progress bar of the process.

basins_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 (Menck et al., 2013) 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.

basins_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 (Datseris and Wagemakers, 2022). 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.

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

statespace_sampler(region [, seed = 42]) → sampler, isinside

A function that facilitates sampling points randomly and uniformly in a state space region. It generates two functions:

• sampler is a 0-argument function that when called generates a random point inside a state space region. The point is always a Vector for type stability irrespectively of dimension. Generally, the generated point should be copied if it needs to be stored. (i.e., calling sampler() utilizes a shared vector) sampler is a thread-safe function.
• isinside is a 1-argument function that returns true if the given state space point is inside the region.

The region can be an instance of any of the following types (input arguments if not specified are vectors of length D, with D the state space dimension):

• HSphere(radius::Real, center): points inside the hypersphere (boundary excluded). Convenience method HSphere(radius::Real, D::Int) makes the center a D-long vector of zeros.
• HSphereSurface(radius, center): points on the hypersphere surface. Same convenience method as above is possible.
• HRectangle(mins, maxs): points in [min, max) for the bounds along each dimension.

The random number generator is always Xoshiro with the given seed.

statespace_sampler(grid::NTuple{N, AbstractRange} [, seed])

If given a grid that is a tuple of AbstractVectors, the minimum and maximum of the vectors are used to make an HRectangle region.

convergence_and_basins_fractions(mapper::AttractorMapper, ics)

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.

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

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

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

basin_entropy(basins::Array{Integer}, ε = size(basins, 1)÷10) -> Sb, Sbb

Return the basin entropy (Daza et al., 2016) 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. The size of the basins has to be divisible by ε.

Assuming that there are $N$ε-boxes that cover the basins, the basin entropy is estimated as (Daza et al., 2016)

\[S_b = \tfrac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{m_i}-p_{ij}\log(p_{ij})\]

where $m_i$ 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 ε.

basins_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. (Puy et al., 2021). 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 of Sbb

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

uncertainty_exponent(basins; kwargs...) -> ε, N_ε, α

Estimate the uncertainty exponent(Grebogi et al., 1983) 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.

test_wada_merge(basins, r) -> p

Test if the 2D array basins has the Wada property using the merging technique of (Daza et al., 2018).

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

edgetracking(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 StateSpaceSets, 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 bisection
• diverge_thresh = 1e-6: distance threshold for parallel integration
• u1: first initial state (defaults to first point in first entry of attractors)
• u2: second initial state (defaults to first point in second entry of attractors)
• maxiter = 100: maximum number of iterations before the algorithm stops
• abstol = 0.0: distance threshold for convergence of the updated edge state
• T_transient = 0.0: transient time before the algorithm starts saving the edge track
• tmax = Inf: maximum integration time of parallel trajectories until re-bisection
• Δt = 0.01: time step passed to step! when evolving the two trajectories
• ϵ_mapper = nothing: ϵ parameter in AttractorsViaProximity
• show_progress = true: if true, shows progress bar and information while running
• verbose = true: if false, silences print output and warnings while running
• kwargs...: additional keyword arguments to be passed to AttractorsViaProximity

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 (Battelino et al., 1988) and further described by (Skufca et al., 2006), the algorithm has been applied to, e.g., the laminar-turbulent boundary in plane Couette flow (Schneider et al., 2008), Wada basins (Wagemakers et al., 2020), as well as Melancholia states in conceptual (Mehling et al., 2023) and intermediate-complexity (Lucarini and Bódai, 2017) 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.

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

EdgeTrackingResults(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 between track1 and track2)
• track1::StateSpaceSet: the pseudo-trajectory tracking the edge within basin 1
• track2::StateSpaceSet: the pseudo-trajectory tracking the edge within basin 2
• time::Vector: time points of the above StateSpaceSets
• bisect_idx::Vector: indices of time at which a re-bisection occurred
• success::Bool: indicates whether the edge tracking has been successful or not

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

tipping_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 (Kaszás et al., 2019).

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 (Kaszás et al., 2019) 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.

minimal_fatal_shock(mapper::AttractorMapper, u0, search_area, algorithm; kw...)

Return the minimal fatal shock (also known as excitability threshold or stability 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 output mfs is a vector like u0.

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_fatal_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 the MFSBruteForce algorithm.
• target_id = nothing: when not nothing, 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 than the one it was originally in. This alternative basin is not returned, do mapper(u0 .+ mfs) if you need the ID.

The minimal fatal shock has many names. Many papers computed this quantity without explicitly naming it, or naming it something simple like "distance to the threshold". The first work that proposed the concept as a nonlocal stability quantifier was by (Klinshov et al., 2015) with the name "stability threshold". Here we use the name of (Halekotte and Feudel, 2020).

Our implementation is generic and works for any dynamical system, using either black box optimization or brute force searching approaches and the unique interface of Attractors.jl for mapping initial conditions to attractors. In contrast to (Klinshov et al., 2015) or (Halekotte and Feudel, 2020), our implementation does not place any assumptions on the nature of the dynamical system, or whether the basin boundaries are smooth.

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.

MFSBlackBoxOptim(; 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 not nothing, 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 of MFSBruteForce that can be used to provide an initial guess.
• bbkwargs = NamedTuple(): additional keyword arguments propagated to BlackBoxOptim.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.

MFSBruteForce(; 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.

Because this algorithm is based on hyperspheres, it assumes the Euclidean norm as the metric.

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

global_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 (Datseris et al., 2023).

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:

• AttractorSeedContinueMatch
• FeaturizeGroupAcrossParameter

Return

1. 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 the i-th parameter combination.
2. attractors_cont::Vector{Dict{Int, <:Any}}. The continued attractors. attractors_cont[i] is a dictionary mapping attractor ID to the attractor set at the i-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 from ics if ics is a function instead of set initial conditions.

AttractorSeedContinueMatch(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:

1. Seed initial conditions from previously found attractors
2. Propagate those forwards to "continue" previous attractors
3. Estimate basin fractions and potentially find new attractors
4. 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 (Datseris et al., 2023).

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 matcherMatchBySSSetDistance 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 (Datseris et al., 2023). This continuation method is still exported as RecurrencesFindAndMatch.

RecurrencesFindAndMatch <: GlobalContinuationAlgorithm
RecurrencesFindAndMatch(mapper::AttractorsViaRecurrences; kwargs...)

A method for global_continuation as in (Datseris et al., 2023) 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 to MatchBySSSetDistance.
• 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 (Datseris et al., 2023).

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 (Datseris et al., 2023), whereby attractors are found using the recurrences algorithm AttractorsViaRecurrences and they are then matched by their distance in state space MatchBySSSetDistance.

aggregate_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

1. fractions_cont: a vector of dictionaries mapping labels to basin fractions.
2. attractors_cont: a vector of dictionaries mapping labels to attractors. 1st and 2nd argument are exactly like the return values of global_continuation with RecurrencesFindAndMatch (or, they can be the return of basins_fractions).
3. featurizer: a 1-argument function to map an attractor into an appropriate feature to be grouped later. Features expected by GroupingConfig are SVector.
4. group_config: a subtype of GroupingConfig.
5. info_extraction: a function accepting a vector of features and returning a description of the features. I.e., exactly as in FeaturizeGroupAcrossParameter. The 5th argument is optional and defaults to the centroid of the features.

Return

1. aggregated_fractions: same as fractions_cont but now contains the fractions of the aggregated attractors.
2. aggregated_info: dictionary mapping the new labels of aggregated_fractions to the extracted information using info_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 StateSpaceSets and outputs a descriptor. E.g., info_extraction = vector -> mean(mean(x) for x in vector).

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

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 the attractors_cont contains in the return of global_continuation.

MatchBySSSetDistance(; 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 to setsofsets_distances.
• threshold = Inf: sets with distance larger than the threshold are guaranteed to not be mapped to each other.
• use_vanished = !isinf(threshold): value of the keyword use_vanished when used in match_sequentially!.

Description

In this matcher the values compared are StateSpaceSets 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.

MatchByBasinEnclosure(; 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 to AttractorsViaProximity. If nothing, it is estimated as a quarter of the minimum distance of centroids (in contrast to the default more accurate estimation in AttractorsViaProximity).
• Δt = 1, consecutive_lost_steps = 1000: also given to AttractorsViaProximity. 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 (Ritchie et al., 2023): attractors that quantify as "basin stable" are matched.

MatchByBasinOverlap(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 Arrays 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

IDMatcher

Supertype of all "matchers" that match can IDs labelling attractors. Currently available matchers:

• MatchBySSSetDistance
• MatchByBasinEnclosure
• MatchByBasinOverlap

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.

matching_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 StateSpaceSets representing attractors. Typically the +,- mean after and before some change of parameter of a dynamical system.

Keyword arguments

• ds: the dynamical system that generated a₊, a₋.
• p, pprev: the parameters corresponding to a₊, a₋. Both need to be iterables mapping parameter index to parameter value (such as Dict, Vector{Pair}, etc., so whatever can be given as input to DynamicalSystems.set_parameters!).
• next_id = next_free_id(a₊, a₋): the ID to give to values of a₊ that cannot be matched to a₋ 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.

matching_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 Arrays which represent basins of attraction and whose elements are the IDs.

See MatchByBasinOverlap for how matching works.

matching_map!(a₊, a₋, matcher; kw...) → rmap

Convenience function that first calls matching_map and then replaces the IDs in a₊ with this rmap.

match_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 (StateSpaceSets), 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 the p, pprev values given to matching_map. See global_continuation if you are unsure what this means.
• ds = nothing: propagated to matching_map.
• retract_keys::Bool = true: If true 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: If use_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 the matcher.

match_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 global continuation. continuation_quantity can for example be fractions_cont from global_continuation.

Centroid(metric = Euclidean())

A distance that can be used in set_distance. The Centroid method returns the distance (according to metric) between the centroids (a.k.a. centers of mass) of the sets.

metric can be any function that takes in two static vectors are returns a positive definite number to use as a distance (and typically is a Metric from Distances.jl).

Hausdorff(metric = Euclidean())

A distance that can be used in set_distance. The Hausdorff distance is the greatest of all the distances from a point in one set to the closest point in the other set. The distance is calculated with the metric given to Hausdorff which defaults to Euclidean.

Hausdorff is 2x slower than StrictlyMinimumDistance, however it is a proper metric in the space of sets of state space sets.

This metric only works for StateSpaceSets whose elements are SVectors.

StrictlyMinimumDistance([brute = false,] [metric = Euclidean(),])

A distance that can be used in set_distance. The StrictlyMinimumDistance returns the minimum distance of all the distances from a point in one set to the closest point in the other set. The distance is calculated with the given metric.

The brute::Bool argument switches the computation between a KDTree-based version, or brute force (i.e., calculation of all distances and picking the smallest one). Brute force performs better for datasets that are either large dimensional or have a small amount of points. Deciding a cutting point is not trivial, and is recommended to simply benchmark the set_distance function to make a decision.

If brute = false this metric only works for StateSpaceSets whose elements are SVectors.

set_distance(ssset1, ssset2 [, distance])

Calculate a distance between two StateSpaceSets, i.e., a distance defined between sets of points, as dictated by distance.

Possible distance types are:

• Centroid, which is the default, and 100s of times faster than the rest
• Hausdorff
• StrictlyMinimumDistance
• Any function f(A, B) that returns the distance between two state space sets A, B.

setsofsets_distances(a₊, a₋ [, distance]) → distances

Calculate distances between sets of StateSpaceSets. Here a₊, a₋ are containers of StateSpaceSets, and the returned distances are dictionaries of distances. Specifically, distances[i][j] is the distance of the set in the i key of a₊ to the j key of a₋. Distances from a₋ to a₊ are not computed at all, assumming symmetry in the distance function.

The distance can be anything valid for set_distance.

unique_keys(v::Iterator{<:AbstractDict})

Given a vector of dictionaries, return a sorted vector of the unique keys that are present across all dictionaries.

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

next_free_id(new::Dict, old::Dict)

Return the minimum key of the "new" dictionary that doesn't exist in the "old" dictionary.

plot_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 the Makie.scatter function that plots the attractors.

heatmap_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 and sckwargs as in plot_attractors.

shaded_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 plot
• maxit = maximum(iterations): clip the values of iterations 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.

plot_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 fraction
• separatorwidth = 1, separatorcolor = "white": adds a line separating the fractions if the style is :band
• axislegend_kwargs = (position = :lt,): propagated to axislegend if a legend is added
• series_kwargs = NamedTuple(): propagated to the band or scatterline plot
• Also all common plotting keywords.

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

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

animate_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 to trajectory 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 the Figure, the Axis, the "bar-like" axis containing the fractions, and the axislegend that adds the legend (if add_legend = true).
• add_legend = true: whether to display the axis legend.