Basins of Attraction

basins_fractions(
    mapper::AttractorMapper,
    ics::Union{StateSpaceSet, Function};
    kwargs...
)

Approximate the state space fractions fs of the basins of attraction of a dynamical system by mapping initial conditions to attractors using mapper (which contains a reference to a DynamicalSystem). The fractions are simply the ratios of how many initial conditions ended up at each attractor.

Initial conditions to use are defined by ics. It can be:

• a StateSpaceSet of initial conditions, in which case all are used.
• a 0-argument function ics() that spits out random initial conditions. 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 continuation.

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

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.

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.

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::StateSpaceSet)

An extension of basins_fractions. Return fs, labels, convergence. The first two are as in basins_fractions, and convergence is a vector containing the time each initial condition took to converge to its attractor. Only usable with mappers that support id = mapper(u0).

See also convergence_time.

Keyword arguments

• show_progress = true: show progress bar.

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, ε = 20) -> Sb, Sbb

Compute the basin entropy (Daza et al., 2016) Sb and basin boundary entropy Sbb of the given basins of attraction by considering ε boxes along each dimension.

Description

First, the input basins is divided regularly into n-dimensional boxes of side ε (along all dimensions). Then Sb is simply the average of the Gibbs entropy computed over these boxes. The function returns the basin entropy Sb as well as the boundary basin entropy Sbb. The later is the average of the entropy only for 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 attractors. 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 Sbb basin_entropy for Ntotal random boxes of radius ε. If the computed value is equal to theoretical value of a smooth boundary (taking into account statistical errors and biases) then we decide that we have a smooth boundary. Notice that the response test_res may depend on the chosen ball radius ε. For larger size, we may observe structures for smooth boundary and we obtain a different answer.

The output test_res is a symbol describing the nature of the basin and the output Sbb is the estimated value of the boundary basin entropy with the sampling method.

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 mfs (also known as excitability threshold) for the initial point u0 according to the specified algorithm given a mapper that satisfies the id = mapper(u0) interface (see AttractorMapper if you are not sure which mappers do that). The mapper contains a reference to a DynamicalSystem. The options for algorithm are: MFSBruteForce or MFSBlackBoxOptim. For high dimensional systems MFSBlackBoxOptim is likely more accurate.

The search_area dictates the state space range for the search of the mfs. It can be a 2-tuple of (min, max) values, in which case the same values are used for each dimension of the system in mapper. Otherwise, it can be a vector of 2-tuples, each for each dimension of the system. The search area is defined w.r.t. to u0 (i.e., it is the search area for perturbations of u0).

An alias to minimal_fata_shock is excitability_threshold.

Keyword arguments

• metric = LinearAlgebra.norm: a metric function that gives the norm of a perturbation vector. This keyword is ignored for 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. It is inspired by the paper "Minimal fatal shocks in multistable complex networks" (Halekotte and Feudel, 2020), however the implementation here is generic: it works for any dynamical system.

The excitability threshold is a concept nearly identical, however, instead of looking for a perturbation that simply brings us out of the basin, we look for the smallest perturbation that brings us into specified basin(s). This is enabled via the keyword target_id.

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.

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