Nonlocal Stability

StabilityQuantifiersAccumulator(bmap::BasinMap [, extras]; kwargs...)

A special data structure that allows mapping initial conditions to their basins while at the same time calculating many stability quantifiers in the most efficient way possible. bmap is any instance of a BasinMap, although for BasinMapFeaturizeGroup the convergence times won't make sense.

The accummulator records several quantifiers of stability (or resilience) defined in (Morr et al., 2026), and a few more added after the publication, including the intermingledness of basins of attraction (Datseris et al., 2026), see list below. However, it also allows computing any additional user-defined quantifier that is a function of the attractors and/or their basins of attraction via the extras argument, see the Extra quantifiers section below.

StabilityQuantifiersAccumulator can be used as any BasinMap with library functions such as basins_fractions. After mapping all initial conditions to their basins, the finalize_accumulator function should be called which will return a dictionary of all stability quantifiers estimated by the accumulator, where each entry maps the stability quantifier description (String) to a dictionary mapping attractor IDs to the stability quantifier value. Calling reset_mapper!(accumulator) cleans up all accumulated quantifiers. This functionality was developed as part of (Morr et al., 2026) and has now been extended to work for any BasinMap, current or future.

Using with global_continuation: Since StabilityQuantifiersAccumulator is formally an BasinMap, it can be used with global_continuation. Simply give it as a bmap input to AttractorSeedContinueMatch and then call global_continuation. Then, the field .quantifiers of GlobalContinuationOutput will be filled with all the quantifiers estimated by the accumulator at each parameter. These will also be matched just like the attractors and their fractions are matched.

Use stability_quantifiers_along_continuation for continuation of stability quantifiers computed on the basis of an BasinMapProximity bmap from already found attractors. This is useful to do for quantifiers related to the convergence time, which is defined more rirogously and is estimated more accurately for a proximity bmap.

Keyword arguments

• finite_time = 1.0: Finite time horizon considered for the computation of the finite_time_basin_stability. Initial coditions with a convergence time larger than finite_time are not considered to be in the respective finite time basin. Convergence time is determined by the bmap.
• weighting_distribution::Distribution: Distribution of uncertain initial conditions used for example in the computation of basin_stability. By default it is a uniform distribution everywhere in the state space.* distance = Centroid(): Set-distance for how to compute the distance between an initial condition u0 and an attractor A. Estimated via set_distance([u0], A, distance).
• idistances = [Euclidean()]: A vector of point distances given to intermingedness for calculating the intermingledness of basins of attraction. If this is empty, intermingledness is not calculated.
• n_basin_entropy = 100: number of nearest neighbors used for sampled-basin entropy.
• basin_entropy = true: if true, compute sampled-basin entropy quantifiers; if false, skip them because they can be expensive.

Description

StabilityQuantifiersAccumulator efficiently uses a single id = bmap(u0) call to accumulate information for many different stability quantifiers corresponding to each attractor of the dynamical system. It accumulates all these different quantifiers when different initial conditions are mapped through it. After enough u0s have been given to the accumulator, they can be finalized (comput maxima or averages) using finalize!(accumulator).

You can extent this functionality by adding new stability quantifiers as long as their estimation can be done on the basis of the three quantities accumulated: the basin dictionary mapping initial conditions to attractors, the distance of each initial condition to each attractor, and the convergence time of each initial condition.

The following stability quantifiers are estimated for each attractor (and the returned dictionary maps strings with the names of the quantifiers to the dictionaries containing the quantifier values for each attractor):

Local (fixed point) stability quantifiers

These quantifiers apply only to fixed point attractors. Their value is NaN if an attractor is not a fixed point (length(A) > 1). If an unstable fixed point attractor is recorded (due to an initial condition starting there for example), a value Inf is assigned to all quantifiers.

• characteristic_return_time: The reciprocal of the largest real part of the eigenvalues of the Jacobian matrix at the fixed point (continuous time), or 1/|log(ρ)| where ρ is the spectral radius (discrete time).
• reactivity: The largest growth rate of the linearized system at the fixed point (continuous time), or ρ - 1 where ρ is the spectral radius (discrete time). See also (Krakovská et al., 2024).
• maximal_amplification: The maximal (with respect to disturbances) amplification of the linearized system at the attractor over all time (continuous time). Always 1 for discrete time systems, as the linearized map does not amplify disturbances beyond their initial norm when the fixed point is stable.
• maximal_amplification_time: The time at which the maximal amplification occurs (continuous time). Always 0 for discrete time systems.

Nonlocal stability quantifiers

The information for these nonlocal stability quantifiers is accumulated while initial conditions are mapped to attractors. Afterwards it is aggregated according to the probability density weighting_distribution when calling finalize_accumulator!.

The word "distance" here refers to the distance established by the distance keyword.

• mean_convergence_time: The convergence time is determined by the bmap using convergence_time. The mean is computed with respect to the weighting_distribution.
• maximal_convergence_time: The maximal convergence time of initial conditions to the attractor. Only initial conditions with non-zero probability under weighting_distribution are considered.
• median_convergence_time: The median convergence time of initial conditions under the weighting_distribution.
• mean_convergence_pace: The mean convergence pace of initial conditions to the attractor. Similar to the mean convergence time, except that each convergence time is divided by the distance of the respective initial condition to the attractor.
• maximal_convergence_pace: The maximal convergence pace of initial conditions to the attractor. Only initial conditions with non-zero probability under weighting_distribution are considered.
• median_convergence_pace: The median convergence pace of initial conditions under the weighting_distribution.
• minimal_critical_shock_magnitude: The minimal distance of the attractor to the closest non-zero probability point (under weighting_distribution) in a basin of attraction of a different attractor. If only a single attractor exists, the value Inf is assigned. Distance is computed via the distance keyword (default Centroid()), so for an aggregated attractor this is the distance from the non-basin point to the centroid of the merged attractor set, which may not be a meaningful state-space point. To instead obtain the minimal distance over all attractor points to any other-basin point, use e.g. distance = StrictlyMinimumDistance().
• maximal_noncritical_shock_magnitude: The distance of the attractor to the furthest non-zero probability point (under weighting_distribution) of its own basin of attraction. If only a single attractor exists, the value Inf is assigned.
• mean_noncritical_shock_magnitude: same as above but computing the mean under weighting_distribution instead of maximum distance.
• basin_fraction: The fraction of initial conditions that converge to the attractor.
• basin_stability: The fraction of initial conditions that converge to the attractor, weighted by weighting_distribution. For the default value of weighting_distribution this is identical to basin_fraction.
• finite_time_basin_stability: The fraction of initial conditions that converge to the attractor within the time horizon finite_time, weighted by weighting_distribution.
• intermingledness$(i): intermingledness of the basins of attraction corresponding to the i-th distance function given to the idistances keyword. Multiple entries may be be produced, each ending with the number i. Because intermingledness is expensive to compute for a large number of initial conditions, you can disable this by providing an empty vector to idistances. See intermingedness for more information.
• basin_entropy and boundary_basin_entropy: estimated from the sampled basins with the nearest-neighbor method (if basin_entropy = true).

Extra quantifiers

The accumulator produces SampledBasinsOfAttraction instances internally. Thus, any function of such objects can be estimated at the end of the accumulation (at the call of finalize_accumulator). To enhance the output of the accumulator with arbitrary quantifiers you can provide the optional argument extras. It is a dictionary mapping additional quantifier names (as Strings) to functions f(sboa, ds). Instead of a dictionary, you can provide an iterable of functions f, in which case the function name is used as the dict's string key.

Each function f takes as an input the sampled basins and the dynamical system and returns a dictionary mapping attractor IDs to an output (anything works, but real numbers make most sense). Each provided function produces an additional entry in the output dictionary of the accumulator, with name provided by its key in the dictionary.

Here's an example extras taken from one of the online examples:

function extra_function(sboa::SampledBasinsOfAttraction, ds::DynamicalSystem)
    ids = extract_basins(sboa)
    u0s = extract_domain(sboa)
    att = extract_attractors(sboa) # we don't need the attractors for this example
    out = Dict{Int, Float64}() # our function must return a Dict
    for i in eachindex(ids)
        id = ids[i]
        out[id] = max(out[id], u0s[i][2])
    end
    return out
end

extras = Dict("maxv" => extra_function)

finalize_accumulator(accumulator::StabilityQuantifiersAccumulator)

Return a dictionary mapping stability quantifiers (strings) to dictionaries mapping attractor IDs to corresponding quantifier values.

The attractors themselves are those of the accumulator's bmap and can be obtained with extract_attractors(accumulator).

See StabilityQuantifiersAccumulator for more.

To compute stability quantifiers for aggregated groups of attractors, merge them first with aggregate_continuation and pass the merged attractors to stability_quantifiers_along_continuation.

stability_quantifiers_along_continuation(
    ds::DynamicalSystem, attractors_cont, pcurve, sampler;
    kw...
)

Perform a global continuation of all stability quantifiers estimated by StabilityQuantifiersAccumulator using the found attractors of a previous call to global_continuation using the ds. The inputs pcurve, sampler are the same as in global_continuation (i.e., a vector and an InitialConditionsSampler).

This method is special because it always creates a BasinMapProximity for the attractors at a given point along the global continuation, and then estimates the stability quantifiers using StabilityQuantifiersAccumulator and the proximity map.

There are two reasons to use this method:

1. You are interested in quantifiers related to the convergence time, which is defined more rirogously and is estimated more accurately for a proximity map.
2. You want more control over the values of ε, finite_time, weighting_distribution, all of which are allowed to be Vectors with the same length as pcurve. (they can always be functions)

Keyword arguments

• ε = nothing: given to BasinMapProximity.
• proximity_mapper_options = NamedTuple(): extra keywords for BasinMapProximity.
• distance, finite_time, weighting_distribution: given to StabilityQuantifiersAccumulator.
• basin_entropy = true: given to StabilityQuantifiersAccumulator to enable/disable sampled-basin entropy quantifiers.
• n_basin_entropy = 100: given to StabilityQuantifiersAccumulator as the number of nearest neighbors used by sampled-basin entropy quantifiers.

Aggregating attractors

This function computes stability quantifiers for whatever attractors it is given. To obtain quantifiers for aggregated groups of attractors, first merge them with aggregate_continuation and pass the resulting agg_attractors_cont here: each merged group is then treated as a single attractor, so all quantifiers (including those that need the raw basin data, like medians and critical shock magnitudes) are computed correctly for the group, with IDs that stay consistent along the parameter axis.

minimal_critical_shock(bmap::BasinMap, u0, search_area, algorithm; kw...)

Return the minimal critical shock for the initial point u0 according to the specified algorithm given a bmap that satisfies the id = bmap(u0) interface (see BasinMap if you are not sure which mappers do that). The output mfs is a vector like u0.

The bmap contains a reference to a DynamicalSystem. The options for algorithm are: MCSBruteForce or MCSBlackBoxOptim. For high dimensional systems MCSBlackBoxOptim 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 bmap. 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_critical_shock is excitability_threshold. Other names for the concept are or stability threshold or minimal fatal shock.

Keyword arguments

• metric = LinearAlgebra.norm: a metric function that gives the norm of a perturbation vector. This keyword is ignored for the MCSBruteForce 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 critical 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 bmap(u0 .+ mfs) if you need the ID.

The minimal critical 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.

MCSBlackBoxOptim(; kwargs...)

The black box derivative-free optimization algorithm used in minimal_critical_shock.

Keyword arguments

• guess = nothing: a initial guess for the minimal critical 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 = MCSBruteForce(100, 100, 0.99): an instance of MCSBruteForce 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, bmap, penalty)
    dist = norm(perturbation)
    if bmap(u0 + perturbation) == bmap(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 MCSBruteForce is used to provide a guess. This can be disabled by either passing a guess vector explicitly or by giving nothing as random_algo.

MCSBruteForce(; kwargs...)

The brute force randomized search algorithm used in minimal_critical_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.
• seed = rand(1:10000): seed for the random number generator used when sampling random perturbations.

intermingledness(points::StatesSpaceSet, labels [, distance]; kw...)

Calculate the intermingledness (Datseris et al., 2026) of the points which have been divided into groups (typically attractors) indicated by labels. Return a dictionary mapping unique labels to their intermingledness. See also boundary_intermingledness.

The optional distance = Euclidean() argument is a function dictating how to estimate distances between points. A vector of distance functions can also be given as distance, in which case a vector of intermingledness is returned corresponding to each distance function. In (Datseris et al., 2026) intermingledness was estimated individually per dimension of points, which you can achieve by e.g.,

using Distances: WeightedEuclidean
D = dimension(points)
weights = [(1:D .== i) for i in 1:D]
distances = WeightedEuclidean.(weights)

The summarizer = mean keyword argument dictates how to summarize the intermingedness statistic across other groups.

Description

Intermingledness is a way to quantify similarity or dissimilarity between the different groups that points was grouped in. For example, points can be initial conditions fed into basin_fractions, and labels the output. Or, points can be feature vectors and labels the output of group_features.

Intermingledness is effectively the ratio of the pairwise-averaged intra-group distance divided by the pairwise-averaged inter-group distance. If it is 1, points are as close to points in their own group as they are to points in other groups. See (Datseris et al., 2026) for examples using intermingledness and the detailed definition or honestly, just look at the source code, it is only 10 lines!

boundary_intermingledness(points::StatesSpaceSet, labels [, distance]; kw...)

Calculate the boundary intermingledness (Datseris et al., 2026) of the points which have been divided into groups (typically attractors) indicated by labels. Return a dictionary mapping unique labels to their b-intermingledness. See also boundary_intermingledness.

The summarizer = mean keyword argument dictates how to summarize the b-intermingedness statistic across other groups.

Description

Boundary interminglendess is an alternative to intermingledness focused on the boundary set, see appendix C of (Datseris et al., 2026) for more.