Basin boundaries

basins_fractal_dimension(BoA::ArrayBasinsOfAttraction; kwargs...) -> V_ε, N_ε, d
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.

The second function signature exists for backwards compatibility.

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(BoA::ArrayBasinsOfAttraction{Integer}, ε = size(BoA.basins, 1)÷10) -> Sb, Sbb
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.

The second function signature exists for backwards compatibility.

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

basin_entropy(BoA::SampledBasinsOfAttraction{<:Integer}, n = 100; kw...) -> Sb, Sbb
basin_entropy(points::AbstractVector, labels::AbstractVector, n = 100; kw...) -> Sb, Sbb

Return the basin entropy Sb and basin boundary entropy Sbb of sampled basins using a nearest-neighbor estimate. For each sampled point, a local neighborhood is formed by its n nearest neighbors, and the probabilities $p_{ij}$ are estimated from the relative frequencies of basin IDs in this neighborhood.

Keyword arguments

• metric = Euclidean(): the metric to use when constructing the KDTree.

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(BoA::ArrayBasinsOfAttraction; kwargs...) -> ε, N_ε, α
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.

The second function signature exists for backwards compatibility.

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(BoA::ArrayBasinsOfAttraction,r) -> p
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).

The second function signature exists for backwards compatibility.

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 BasinMapProximity 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 = 1e-9: 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.
• show_progress = true: if true, shows progress bar and information while running.
• verbose = true: if false, silences print output and warnings while running.
• kw...: additional keyword arguments to be passed to BasinMapProximity. We strongly recommend to either pass in a high Ttr, or a very small ε, (smaller than the default estimated by BasinMapProximity) to avoid transient parts of trajectories wrongly being classified as having converged.

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., 2024) 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 BasinMap 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.

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, bmap::BasinMap; 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 bmap 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 bmap must be an BasinMap of subtype BasinMapProximity or BasinMapRecurrences.

stagger_and_step(ds::DynamicalSystem, x0, N::Int, 
    isinside::Function; kwargs...) -> trajectory

Implement the stagger-and-step method (Sweet et al., 2001) and (Sala et al., 2016) to approximate the invariant non-attracting set governing the chaotic transient dynamics of a system, namely the stable manifold of a chaotic saddle.

Before starting the algorithm, we must search a region of the phase space defined by a membership function isinside that fulfils three conditions:

1. The function isinside(x) must return true if the point x is inside the chosen bounded region and false otherwise. See statespace_sampler as a helper to construct this
2. The region contains no attractors.
3. The stable manifold of the saddle must be inside the defined region.

Given the dynamical system ds and a initial guess x0, the algorithm returns a StateSpaceState named trajectory of N points from the phase space close to the stable manifold. If we set one of this point as an initial condition of ds, the trajectory escape from the region after at least Tm steps of ds. The search is stochastic and depends on the parameter δ defining a (small) neighborhood of search.

Keyword arguments

• δ = 1e-10: A small number constraining the random search around a particular point. The interpretation of this number will depend on the probability distribution chosen for the sampling (see stagger_mode).

• Tm = 30: The minimum number of steps of ds before the trajectory escapes from the bounding box defined by isinside.

• max_steps = 10^5: The search for a new candidate point may fail at some point. If the search fails after max_steps, a new initial point is set and the method starts from a new point.

• max_escape_time = 10000: If the trajectory stays in the defined region after max_escape_time steps, there is probably an attractor in the region and the algorithm will throw an error.

• stagger_mode = :exp: There are several ways to produce candidate points x to fulfill the condition T(x) > Tm. The available methods are:

  • :exp: A candidate is sampled from a truncated exponential distribution in a random direction u around the current x such that x_c = x + u*r. r = 10^-s with s taken from a uniform distribution in [-15, δ].

  • :unif: The next candidate is x_c = x + u*r with r taken from a uniform distribution [0,δ] and u a random direction around x.

  • :adaptive: The next candidate is x_c = x + u*r with r drawn from a gaussian distribution with variance δ and mean zero. The variance is adapted according to a free parameter γ such that: δ = δ/γ if no candidate is found and δ = δ*γ when it succeeds.

• γ = 1.1: It is free parameter for the adaptive stagger method :adaptive.

• δ₀ = 1.0: This is the radius for the first stagger trajectory search. The algorithm looks for a point sufficiently close to the saddle before switching to the stagger-and-step routine. The search radius must be large enough to find a suitable initial candidate. To type δ₀ use \delta<TAB>\_0<TAB>.

• rng::AbstractRNG = Xoshiro(): Random number generator. Use this for reproducibility.

Description

The method relies on the stagger-and-step algorithm that searches initial conditions close to the saddle with escapes time T(x_n) > Tm. The function T represents the time at which the trajectory with initial condition x_n steps out from a region defined by the user (see the argument isinside). This time is a discrete or continuous variable depending on the dynamical system ds. In the case of a continuous dynamical system, a fixed timestep Δt is recommended.

Given the dynamical mapping F, if the step x_{n+1} = F(x_n) fulfills the condition T(x_{n+1}) > Tm we accept this next point, this is the step part of the method. If not, the method search randomly the next point in a neighborhood following a given probability distribution, this is the stagger part. The stagger process sometimes fails to find a new candidate and a new starting point of the trajectory is chosen within the defined region. See the keyword argument stagger_mode for the different available methods.

The method produces a pseudo-trajectory of N points δ-close to the stable manifold of the chaotic saddle.