Attractors.jl Tutorial

Attractors is a component of the DynamicalSystems.jl library. This tutorial will walk you through its main functionality. That is, given a DynamicalSystem instance, find all its attractors and their basins of attraction. Then, continue these attractors, and their stability properties, across a parameter value. It also offers various functions that compute nonlocal stability properties for an attractor, any of which can be used in the continuation to quantify stability.

Besides this main functionality, there are plenty of other stuff, like for example edgestate or basins_fractal_dimension, but we won't cover anything else in this introductory tutorial. See the examples page instead.

Package versions used

import Pkg

Pkg.status(["Attractors", "CairoMakie", "OrdinaryDiffEq"])

Status `~/work/Attractors.jl/Attractors.jl/docs/Project.toml`
  [f3fd9213] Attractors v1.22.1 `~/work/Attractors.jl/Attractors.jl`
  [13f3f980] CairoMakie v0.12.15
  [1dea7af3] OrdinaryDiffEq v6.89.0

Tutorial - copy-pasteable version

Gotta go fast!

using Attractors, CairoMakie, OrdinaryDiffEq
## Define key input: a `DynamicalSystem`
function modified_lorenz_rule(u, p, t)
    x, y, z = u; a, b = p
    dx = y - x
    dy = - x*z + b*abs(z)
    dz = x*y - a
    return SVector(dx, dy, dz)
end
p0 = [5.0, 0.1] # parameters
u0 = [-4.0, 5, 0] # state
diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9, dt = 0.01) # solver options
ds = CoupledODEs(modified_lorenz_rule, u0, p0; diffeq)

## Define key input: an `AttractorMaper` that finds
## attractors of a `DynamicalSystem`
grid = (
    range(-15.0, 15.0; length = 150), # x
    range(-20.0, 20.0; length = 150), # y
    range(-20.0, 20.0; length = 150), # z
)
mapper = AttractorsViaRecurrences(ds, grid;
    consecutive_recurrences = 1000,
    consecutive_lost_steps = 100,
)

## Find attractors and their basins of attraction state space fraction
## by randomly sampling initial conditions in state sapce
sampler, = statespace_sampler(grid)
algo = AttractorSeedContinueMatch(mapper)
fs = basins_fractions(mapper, sampler)
attractors = extract_attractors(mapper)

## found two attractors: one is a limit cycle, the other is chaotic
## visualize them
plot_attractors(attractors)

## continue all attractors and their basin fractions across any arbigrary
## curve in parameter space using a global continuation algorithm
algo = AttractorSeedContinueMatch(mapper)
params(θ) = [1 => 5 + 0.5cos(θ), 2 => 0.1 + 0.01sin(θ)]
pcurve = params.(range(0, 2π; length = 101))
fractions_cont, attractors_cont = global_continuation(
	algo, pcurve, sampler; samples_per_parameter = 1_000
)

## and visualize the results
fig = plot_basins_attractors_curves(
	fractions_cont, attractors_cont, A -> minimum(A[:, 1]), pcurve; add_legend = false
)

Input: a `DynamicalSystem`

The key input for most functionality of Attractors.jl is an instance of a DynamicalSystem. If you don't know how to make a DynamicalSystem, you need to consult the main tutorial of the DynamicalSystems.jl library. For this tutorial we will use a modified Lorenz-like system with equations

\[\begin{align*} \dot{x} & = y - x \\ \dot{y} &= -x*z + b*|z| \\ \dot{z} &= x*y - a \\ \end{align*}\]

which we define in code as

using Attractors # part of `DynamicalSystems`, so it re-exports functionality for making them!
using OrdinaryDiffEq # for accessing advanced ODE Solvers

function modified_lorenz_rule(u, p, t)
    x, y, z = u; a, b = p
    dx = y - x
    dy = - x*z + b*abs(z)
    dz = x*y - a
    return SVector(dx, dy, dz)
end

p0 = [5.0, 0.1] # parameters
u0 = [-4.0, 5, 0] # state
diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9, dt = 0.01) # solver options
ds = CoupledODEs(modified_lorenz_rule, u0, p0; diffeq)

3-dimensional CoupledODEs
 deterministic: true
 discrete time: false
 in-place:      false
 dynamic rule:  modified_lorenz_rule
 ODE solver:    Vern9
 ODE kwargs:    (abstol = 1.0e-9, reltol = 1.0e-9, dt = 0.01)
 parameters:    [5.0, 0.1]
 time:          0.0
 state:         [-4.0, 5.0, 0.0]

Finding attractors

There are two major methods for finding attractors in dynamical systems. Explanation of how they work is in their respective docs.

You can consult (Datseris et al., 2023) for a comparison between the two.

As far as the user is concerned, both algorithms are part of the same interface, and can be used in the same way. The interface is extendable as well, and works as follows.

First, we create an instance of such an "attractor finding algorithm", which we call AttractorMapper. For example, AttractorsViaRecurrences requires a tesselated grid of the state space to search for attractors in. It also allows the user to tune some meta parameters, but in our example they are already tuned for the dynamical system at hand. So we initialize

grid = (
    range(-10.0, 10.0; length = 150), # x
    range(-15.0, 15.0; length = 150), # y
    range(-15.0, 15.0; length = 150), # z
)

mapper = AttractorsViaRecurrences(ds, grid;
    consecutive_recurrences = 1000, attractor_locate_steps = 1000,
    consecutive_lost_steps = 100,
)

AttractorsViaRecurrences
 system:      CoupledODEs
 grid:        (-10.0:0.1342281879194631:10.0, -15.0:0.20134228187919462:15.0, -15.0:0.20134228187919462:15.0)
 attractors:  Dict{Int64, StateSpaceSet{3, Float64, SVector{3, Float64}}}()

This mapper can map any initial condition to the corresponding attractor ID, for example

mapper([-4.0, 5, 0])

while

mapper([4.0, 2, 0])

the fact that these two different initial conditions got assigned different IDs means that they converged to a different attractor. The attractors are stored in the mapper internally, to obtain them we use the function

attractors = extract_attractors(mapper)

Dict{Int64, StateSpaceSet{3, Float64, SVector{3, Float64}}} with 2 entries:
  2 => 3-dimensional StateSpaceSet{Float64} with 320 points
  1 => 3-dimensional StateSpaceSet{Float64} with 934 points

In Attractors.jl, all information regarding attractors is always a standard Julia Dict, which maps attractor IDs (positive integers) to the corresponding quantity. Here the quantity are the attractors themselves, represented as StateSpaceSet.

We can visualize them with the convenience plotting function

using CairoMakie
plot_attractors(attractors)

(this convenience function is a simple loop over scattering the values of the attractors dictionary)

In our example system we see that for the chosen parameters there are two coexisting attractors: a limit cycle and a chaotic attractor. There may be more attractors though! We've only checked two initial conditions, so we could have found at most two attractors! However, it can get tedious to manually iterate over initial conditions, which is why this mapper is typically given to higher level functions for finding attractors and their basins of attraction. The simplest one is basins_fractions. Using the mapper, it finds "all" attractors of the dynamical system and reports the state space fraction each attractors attracts. The search is probabilistic, so "all" attractors means those that at least one initial condition converged to.

We can provide explicitly initial conditions to basins_fraction, however it is typically simpler to provide it with with a state space sampler instead: a function that generates random initial conditions in the region of the state space that we are interested in. Here this region coincides with grid, so we can simply do:

sampler, = statespace_sampler(grid)

sampler() # random i.c.

3-element Vector{Float64}:
 -2.645638264032062
  1.047314055517088
  7.278008671940473

sampler() # another random i.c.

3-element Vector{Float64}:
  8.679171398887526
 11.06137754367522
 -2.9299278174265613

and finally call

fs = basins_fractions(mapper, sampler)

Dict{Int64, Float64} with 2 entries:
  2 => 0.357
  1 => 0.643

The returned fs is a dictionary mapping each attractor ID to the fraction of the state space the corresponding basin occupies. With this we can confirm that there are (likely) only two attractors and that both attractors are robust as both have sufficiently large basin fractions.

To obtain the full basins, which is computationally much more expensive, use basins_of_attraction.

You can use alternative algorithms in basins_fractions, see the documentation of AttractorMapper for possible subtypes. AttractorMapper defines an extendable interface and can be enriched with other methods in the future!

Different Attractor Mapper

Attractors.jl utilizes composable interfaces throughout its functionality. In the above example we used one particular method to find attractors, via recurrences in the state space. An alternative is AttractorsViaFeaturizing.

For this method, we need to provide a "featurizing" function that given an trajectory (which is likely an attractor), it returns some features that will hopefully distinguish different attractors in a subsequent grouping step. Finding good features is typically a trial-and-error process, but for our system we already have some good features:

using Statistics: mean

function featurizer(A, t) # t is the time vector associated with trajectory A
    xmin = minimum(A[:, 1])
    ycen = mean(A[:, 2])
    return SVector(xmin, ycen)
end

featurizer (generic function with 1 method)

from which we initialize

mapper2 = AttractorsViaFeaturizing(ds, featurizer; Δt = 0.1)

AttractorsViaFeaturizing
 system:      CoupledODEs
 Ttr:         100.0
 Δt:          0.1
 T:           100.0
 group via:   GroupViaClustering
 featurizer:  featurizer

AttractorsViaFeaturizing allows for a third input, which is a "grouping configuration", that dictates how features will be grouped into attractors, as features are extracted from (randomly) sampled state space trajectories. In this tutorial we leave it at its default value, which is clustering using the DBSCAN algorithm. The keyword arguments are meta parameters which control how long to integrate each initial condition for, and what sampling time, to produce a trajectory A given to the featurizer function. Because one of the two attractors is chaotic, we need denser sampling time than the default.

We can use mapper2 exactly as mapper:

fs2 = basins_fractions(mapper2, sampler)

attractors2 = extract_attractors(mapper2)

plot_attractors(attractors2)

This mapper also found the attractors, but we should warn you: this mapper is less robust than AttractorsViaRecurrences. One of the reasons for this is that AttractorsViaFeaturizing is not auto-terminating. For example, if we do not have enough transient integration time, the two attractors will get confused into one:

mapper3 = AttractorsViaFeaturizing(ds, featurizer; Ttr = 10, Δt = 0.1)
fs3 = basins_fractions(mapper3, sampler)
attractors3 = extract_attractors(mapper3)
plot_attractors(attractors3)

On the other hand, the downside of AttractorsViaRecurrences is that it can take quite a while to converge for chaotic high dimensional systems.

Global continuation

If you have heard before the word "continuation", then you are likely aware of the traditional continuation-based bifurcation analysis (CBA) offered by many software, such as AUTO, MatCont, and in Julia BifurcationKit.jl. Here we offer a completely different kind of continuation called global continuation.

The traditional continuation analysis continues the curves of individual fixed points (and under some conditions limit cycles) across the joint state-parameter space and tracks their local (linear) stability. This approach needs to manually be "re-run" for every individual branch of fixed points or limit cycles. The global continuation in Attractors.jl finds all attractors, including chaotic or quasiperiodic ones, in the whole of the state space (that it searches in), without manual intervention. It then continues all of these attractors concurrently along a parameter axis. Additionally, the global continuation tracks a nonlocal stability property which by default is the basin fraction.

This is a fundamental difference. Because all attractors are simultaneously tracked across the parameter axis, the user may arbitrarily estimate any property of the attractors and how it varies as the parameter varies. A more detailed comparison between these two approaches can be found in (Datseris et al., 2023). See also the comparison page in our docs that attempts to do the same analysis of our Tutorial with traditional continuation software.

To perform the continuation is extremely simple. First, we decide what parameter, and what range, to continue over:

prange = 4.5:0.01:6
pidx = 1 # index of the parameter

Then, we may call the global_continuation function. We have to provide a continuation algorithm, which itself references an AttractorMapper. In this example we will re-use the mapper to create the "flagship product" of Attractors.jl which is the geenral AttractorSeedContinueMatch. This algorithm uses the mapper to find all attractors at each parameter value and from the found attractors it continues them along a parameter axis using a seeding process (see its documentation string). Then, it performs a "matching" step, ensuring a "continuity" of the attractor label across the parameter axis. For now we ignore the matching step, leaving it to the default value. We'll use the mapper we created above and define

ascm = AttractorSeedContinueMatch(mapper)

AttractorSeedContinueMatch{AttractorsViaRecurrences{CoupledODEs{false, 3, OrdinaryDiffEqCore.ODEIntegrator{OrdinaryDiffEqVerner.Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, false, SVector{3, Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{SVector{3, Float64}}, SciMLBase.ODESolution{Float64, 2, Vector{SVector{3, Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, Nothing, SciMLBase.ODEProblem{SVector{3, Float64}, Tuple{Float64, Float64}, false, Vector{Float64}, SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.modified_lorenz_rule), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqVerner.Vern9{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.modified_lorenz_rule), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{SVector{3, Float64}}, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, Nothing, OrdinaryDiffEqVerner.Vern9ConstantCache, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.modified_lorenz_rule), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqVerner.Vern9ConstantCache, OrdinaryDiffEqCore.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEqCore.PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Bool, SciMLBase.CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Tuple{}, Tuple{}, Tuple{}}, SVector{3, Float64}, Float64, Nothing, OrdinaryDiffEqCore.DefaultInit, Nothing}, Vector{Float64}}, Attractors.BasinsInfo{3, Attractors.RegularGrid{3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, Float64, SVector{3, Float64}, Attractors.SparseArray{Int64, 3}}, Attractors.RegularGrid{3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol, Symbol}, @NamedTuple{consecutive_recurrences::Int64, attractor_locate_steps::Int64, consecutive_lost_steps::Int64}}}, MatchBySSSetDistance{Centroid{Euclidean}, Float64}, typeof(Attractors._default_seeding)}(AttractorsViaRecurrences
 system:      CoupledODEs
 grid:        (-10.0:0.1342281879194631:10.0, -15.0:0.20134228187919462:15.0, -15.0:0.20134228187919462:15.0)
 attractors:  Dict{Int64, StateSpaceSet{3, Float64, SVector{3, Float64}}}(2 => 3-dimensional StateSpaceSet{Float64} with 320 points, 1 => 3-dimensional StateSpaceSet{Float64} with 934 points)
, MatchBySSSetDistance{Centroid{Euclidean}, Float64}(Centroid{Euclidean}(Euclidean(0.0)), Inf, false), Attractors._default_seeding)

and call

fractions_cont, attractors_cont = global_continuation(
	ascm, prange, pidx, sampler; samples_per_parameter = 1_000
)

([Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0), Dict(1 => 1.0)  …  Dict(2 => 0.9201596806387226, 3 => 0.07984031936127745), Dict(2 => 1.0), Dict(2 => 1.0), Dict(2 => 1.0), Dict(2 => 1.0), Dict(2 => 1.0), Dict(2 => 1.0), Dict(2 => 1.0), Dict(2 => 1.0), Dict(2 => 1.0)], Dict{Int64, StateSpaceSet{3, Float64, SVector{3, Float64}}}[Dict(1 => 3-dimensional StateSpaceSet{Float64} with 832 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 749 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 748 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 797 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 831 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 666 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 601 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 785 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 698 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 716 points)  …  Dict(2 => 3-dimensional StateSpaceSet{Float64} with 489 points, 3 => 3-dimensional StateSpaceSet{Float64} with 410 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 503 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 457 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 449 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 429 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 361 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 344 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 337 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 336 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 336 points)])

the output is given as two vectors. Each vector is a dictionary mapping attractor IDs to their basin fractions, or their state space sets, respectively. Both vectors have the same size as the parameter range. For example, the attractors at the 34-th parameter value are:

attractors_cont[34]

Dict{Int64, StateSpaceSet{3, Float64, SVector{3, Float64}}} with 1 entry:
  1 => 3-dimensional StateSpaceSet{Float64} with 492 points

There is a fantastic convenience function for animating the attractors evolution, that utilizes things we have already defined:

animate_attractors_continuation(
    ds, attractors_cont, fractions_cont, prange, pidx;
);

Hah, how cool is that! The attractors pop in and out of existence like out of nowhere! It would be incredibly difficult to find these attractors in traditional continuation software where a rough estimate of the period is required! (It would also be too hard due to the presence of chaos for most of the parameter values, but that's another issue!)

Now typically a continuation is visualized in a 2D plot where the x axis is the parameter axis. We can do this with the convenience function:

fig = plot_basins_attractors_curves(
	fractions_cont, attractors_cont, A -> minimum(A[:, 1]), prange,
)

In the top panel are the basin fractions, by default plotted as stacked bars. Bottom panel is a visualization of the tracked attractors. The argument A -> minimum(A[:, 1]) is simply a function that maps an attractor into a real number for plotting.

Different matching procedures

By default attractors are matched by their distance in state space. The default matcher is MatchBySSSetDistance, and is given implicitly as a default 2nd argument when creating AttractorSeedContinueMatch. But like anything else in Attractors.jl, "matchers" also follow a well-defined and extendable interface, see IDMatchers for that.

Let's say that the default matching that we chose above isn't desirable. For example, one may argue that the attractor that pops up at the end of the continuation should have been assigned the same ID as attractor 1, because they are both to the left (see the video above). In reality one wouldn't really request that, because looking the video of attractors above shows that the attractors labelled "1", "2", and "3" are all completely different. But we argue here for example that "3" should have been the same as "1".

Thankfully during a global continuation the "matching" step is completely separated from the "finding and continuing" step. If we don't like the initial matching, we can call match_sequentially! with a new instance of a matcher, and match again, without having to recompute the attractors and their basin fractions. For example, using this matcher:

matcher = MatchBySSSetDistance(use_vanished = true)

MatchBySSSetDistance{Centroid{Euclidean}, Float64}(Centroid{Euclidean}(Euclidean(0.0)), Inf, true)

will compare a new attractor with the latest instance of attractors with a given ID that have ever existed, irrespectively if they exist in the current parameter or not. This means, that the attractor "3" would in fact be compared with both attractor "2" and "1", even if "1" doesn't exist in the parameter "3" started existing at. And because "3" is closer to "1" than to "2", it will get matched to attractor "1" and get the same ID.

Let's see this in action:

attractors_cont2 = deepcopy(attractors_cont)

match_sequentially!(attractors_cont2, matcher)

fig = plot_attractors_curves(
	attractors_cont2, A -> minimum(A[:, 1]), prange,
)

and as we can see, the new attractor at the end of the parameter range got assigned the same ID as the original attractor "1". For more ways of matching attractors see IDMatcher.

Enhancing the continuation

The biggest strength of Attractors.jl is that it is not an isolated software. It is part of DynamicalSystems.jl. Here, we will use the full power of DynamicalSystems.jl and enrich the above continuation with various other measures of nonlocal stability, in particular Lyapunov exponents and the minimal fatal shock. First, let's plot again the continuation and label some things or clarity

fig = plot_basins_attractors_curves(
	fractions_cont, attractors_cont, A -> minimum(A[:, 1]), prange; add_legend = false
)

ax1 = content(fig[2,1])

ax1.ylabel = "min(A₁)"

fig

First, let's estimate the maximum Lyapunov exponent (MLE) for all attractors, using the lyapunov function that comes from the ChaosTools.jl submodule.

using ChaosTools: lyapunov

lis = map(enumerate(prange)) do (i, p) # loop over parameters
    set_parameter!(ds, pidx, p) # important! We use the dynamical system!
    attractors = attractors_cont[i]
    # Return a dictionary mapping attractor IDs to their MLE
    Dict(k => lyapunov(ds, 10000.0; u0 = A[1]) for (k, A) in attractors)
end

151-element Vector{Dict{Int64, Float64}}:
 Dict(1 => 0.11063040059085028)
 Dict(1 => 0.08701055206081301)
 Dict(1 => 0.09438114119318956)
 Dict(1 => 0.09431940240352951)
 Dict(1 => 0.08537075228426165)
 Dict(1 => 0.04030058499652367)
 Dict(1 => 8.51054506401255e-5)
 Dict(1 => 0.08943243194967629)
 Dict(1 => 0.08542254292087191)
 Dict(1 => 0.08295323081952044)
 ⋮
 Dict(2 => 0.0002349420201453823)
 Dict(2 => 0.00016701888055929077)
 Dict(2 => 0.00033949274914120946)
 Dict(2 => 0.00015453123165271555)
 Dict(2 => -4.7922207539117585e-5)
 Dict(2 => 6.78146409068883e-5)
 Dict(2 => 9.898874500937138e-5)
 Dict(2 => -0.00012005471291116607)
 Dict(2 => -4.197578470353276e-6)

The above map loop may be intimidating if you are a beginner, but it is really just a shorter way to write a for loop for our example. We iterate over all parameters, and for each we first update the dynamical system with the correct parameter, and then extract the MLE for each attractor. map just means that we don't have to pre-allocate a new vector before the loop; it creates it for us.

We can visualize the LE with the other convenience function plot_continuation_curves!,

ax2 = Axis(fig[3, 1]; ylabel = "MLE")
plot_continuation_curves!(ax2, lis, prange; add_legend = false)

fig

This reveals crucial information for tha attractors, whether they are chaotic or not, that we would otherwise obtain only by visualizing the system dynamics at every single parameter. The story we can see now is that the dynamics start with a limit cycle (0 Lyapunov exponent), go into bi-stability of chaos and limit cycle, then there is only one limit cycle again, and then a chaotic attractor appears again, for a second bistable regime.

The last piece of information to add is yet another measure of nonlocal stability: the minimal fatal shock (MFS), which is provided by minimal_fatal_shock. The code to estimate this is similar with the map block for the MLE. Here however we re-use the created mapper, but now we must not forget to reset it inbetween parameter increments:

using LinearAlgebra: norm
search_area = collect(extrema.(grid ./ 2)) # smaller search = faster results
search_algorithm = MFSBlackBoxOptim(max_steps = 1000, guess = ones(3))

mfss = map(enumerate(prange)) do (i, p)
    set_parameter!(ds, pidx, p)
    reset_mapper!(mapper) # reset so that we don't have to re-initialize
    # We need a special clause here: if there is only 1 attractor,
    # then there is no MFS. It is undefined. We set it to `NaN`,
    # which conveniently, will result to nothing being plotted by Makie.
    attractors = attractors_cont[i]
    if length(attractors) == 1
        return Dict(k => NaN for (k, A) in attractors)
    end
    # otherwise, compute the actual MFS from the first point of each attractor
    Dict(k =>
        norm(minimal_fatal_shock(mapper, A[1], search_area, search_algorithm))
        for (k, A) in attractors
    )
end

151-element Vector{Dict{Int64, Float64}}:
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 Dict(1 => NaN)
 ⋮
 Dict(2 => NaN)
 Dict(2 => NaN)
 Dict(2 => NaN)
 Dict(2 => NaN)
 Dict(2 => NaN)
 Dict(2 => NaN)
 Dict(2 => NaN)
 Dict(2 => NaN)
 Dict(2 => NaN)

In a real application we wouldn't use the first point of each attractor, as the first point is completely random on the attractor (at least, for the [AttractorsViaRecurrences] mapper we use here). We would do this by examining the whole A object in the above block instead of just using A[1]. But this is a tutorial so we don't care!

Right, so now we can visualize the MFS with the rest of the other quantities:

ax3 = Axis(fig[4, 1]; ylabel = "MFS", xlabel = "parameter")
plot_continuation_curves!(ax3, mfss, prange; add_legend = false)

# make the figure prettier
for ax in (ax1, ax2,); hidexdecorations!(ax; grid = false); end
resize!(fig, 500, 500)
fig

Continuation along arbitrary parameter curves

One of the many advantages of the global continuation is that we can choose what parameters to continue over. We can provide any arbitrary curve in parameter space. This is possible because (1) finding and matching attractors are two completely orthogonal steps, and (2) it is completely fine for attractors to dissapear (and perhaps re-appear) during a global continuation.

#For example, we can probe an elipsoid defined as

params(θ) = [1 => 5 + 0.5cos(θ), 2 => 0.1 + 0.01sin(θ)]
pcurve = params.(range(0, 2π; length = 101))

101-element Vector{Vector{Pair{Int64, Float64}}}:
 [1 => 5.5, 2 => 0.1]
 [1 => 5.499013364214136, 2 => 0.10062790519529313]
 [1 => 5.496057350657239, 2 => 0.10125333233564304]
 [1 => 5.491143625364344, 2 => 0.10187381314585725]
 [1 => 5.4842915805643155, 2 => 0.10248689887164855]
 [1 => 5.475528258147577, 2 => 0.10309016994374948]
 [1 => 5.4648882429441255, 2 => 0.10368124552684678]
 [1 => 5.45241352623301, 2 => 0.10425779291565074]
 [1 => 5.438153340021932, 2 => 0.10481753674101715]
 [1 => 5.422163962751007, 2 => 0.10535826794978997]
 ⋮
 [1 => 5.438153340021932, 2 => 0.09518246325898286]
 [1 => 5.45241352623301, 2 => 0.09574220708434927]
 [1 => 5.4648882429441255, 2 => 0.09631875447315323]
 [1 => 5.475528258147577, 2 => 0.09690983005625053]
 [1 => 5.4842915805643155, 2 => 0.09751310112835145]
 [1 => 5.491143625364344, 2 => 0.09812618685414276]
 [1 => 5.496057350657239, 2 => 0.09874666766435695]
 [1 => 5.499013364214136, 2 => 0.09937209480470688]
 [1 => 5.5, 2 => 0.1]

here each component maps the parameter index to its value. We can just give this pcurve to the global continuation, using the same mapper and continuation algorithm, but adjusting the matching process so that vanished attractors are kept in "memory"

matcher = MatchBySSSetDistance(use_vanished = true)

ascm = AttractorSeedContinueMatch(mapper, matcher)

fractions_cont, attractors_cont = global_continuation(
	ascm, pcurve, sampler; samples_per_parameter = 1_000
)

([Dict(2 => 0.73, 1 => 0.27), Dict(2 => 0.7594810379241517, 1 => 0.2405189620758483), Dict(2 => 0.7295409181636726, 1 => 0.27045908183632733), Dict(2 => 0.7415169660678643, 1 => 0.25848303393213573), Dict(2 => 0.7385229540918163, 1 => 0.26147704590818366), Dict(2 => 0.7415169660678643, 1 => 0.25848303393213573), Dict(2 => 0.7544910179640718, 1 => 0.24550898203592814), Dict(2 => 0.7315369261477046, 1 => 0.2684630738522954), Dict(2 => 0.7564870259481038, 1 => 0.2435129740518962), Dict(2 => 0.7355289421157685, 1 => 0.26447105788423153)  …  Dict(2 => 0.7025948103792415, 1 => 0.29740518962075846), Dict(2 => 0.6976047904191617, 1 => 0.3023952095808383), Dict(2 => 0.718562874251497, 1 => 0.281437125748503), Dict(2 => 0.7524950099800399, 1 => 0.24750499001996007), Dict(2 => 0.7215568862275449, 1 => 0.27844311377245506), Dict(2 => 0.7395209580838323, 1 => 0.26047904191616766), Dict(2 => 0.7405189620758483, 1 => 0.25948103792415167), Dict(2 => 0.7265469061876247, 1 => 0.27345309381237526), Dict(2 => 0.7325349301397206, 1 => 0.26746506986027946), Dict(2 => 0.7445109780439122, 1 => 0.2554890219560878)], Dict{Int64, StateSpaceSet{3, Float64, SVector{3, Float64}}}[Dict(2 => 3-dimensional StateSpaceSet{Float64} with 518 points, 1 => 3-dimensional StateSpaceSet{Float64} with 340 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 480 points, 1 => 3-dimensional StateSpaceSet{Float64} with 313 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 563 points, 1 => 3-dimensional StateSpaceSet{Float64} with 326 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 618 points, 1 => 3-dimensional StateSpaceSet{Float64} with 332 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 632 points, 1 => 3-dimensional StateSpaceSet{Float64} with 323 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 646 points, 1 => 3-dimensional StateSpaceSet{Float64} with 344 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 658 points, 1 => 3-dimensional StateSpaceSet{Float64} with 320 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 652 points, 1 => 3-dimensional StateSpaceSet{Float64} with 345 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 695 points, 1 => 3-dimensional StateSpaceSet{Float64} with 304 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 669 points, 1 => 3-dimensional StateSpaceSet{Float64} with 316 points)  …  Dict(2 => 3-dimensional StateSpaceSet{Float64} with 647 points, 1 => 3-dimensional StateSpaceSet{Float64} with 315 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 643 points, 1 => 3-dimensional StateSpaceSet{Float64} with 312 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 591 points, 1 => 3-dimensional StateSpaceSet{Float64} with 354 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 605 points, 1 => 3-dimensional StateSpaceSet{Float64} with 329 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 573 points, 1 => 3-dimensional StateSpaceSet{Float64} with 325 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 537 points, 1 => 3-dimensional StateSpaceSet{Float64} with 328 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 483 points, 1 => 3-dimensional StateSpaceSet{Float64} with 341 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 509 points, 1 => 3-dimensional StateSpaceSet{Float64} with 335 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 510 points, 1 => 3-dimensional StateSpaceSet{Float64} with 324 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 487 points, 1 => 3-dimensional StateSpaceSet{Float64} with 339 points)])

and animate the result

animate_attractors_continuation(
    ds, attractors_cont, fractions_cont, pcurve;
    savename = "curvecont.mp4"
);

Conclusion and comparison with traditional local continuation

We've reached the end of the tutorial! Some aspects we haven't highlighted is how most of the infrastructure of Attractors.jl is fully extendable. You will see this when reading the documentation strings of key structures like AttractorMapper. All documentation strings are in the API page. See the examples page for more varied applications. And lastly, see the comparison page in our docs that attempts to do the same analysis of our Tutorial with traditional local continuation and bifurcation analysis software showing that (at least for this example) using Attractors.jl is clearly beneficial over the alternatives.