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.19.4 `~/work/Attractors.jl/Attractors.jl`
[13f3f980] CairoMakie v0.12.7
[1dea7af3] OrdinaryDiffEq v6.87.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 a parameter axis
# using a global continuation algorithm
algo = AttractorSeedContinueMatch(mapper)
prange = 4.7:0.02:6
pidx = 1
fractions_cont, attractors_cont = global_continuation(
algo, prange, pidx, sampler; samples_per_parameter = 1_000
)
# and visualize the results
fig = plot_basins_attractors_curves(
fractions_cont, attractors_cont, A -> minimum(A[:, 1]), prange; 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}}()
This mapper
can map any initial condition to the corresponding attractor ID, for example
mapper([-4.0, 5, 0])
1
while
mapper([4.0, 2, 0])
2
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}} with 2 entries:
2 => 3-dimensional StateSpaceSet{Float64} with 314 points
1 => 3-dimensional StateSpaceSet{Float64} with 893 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)
![Example block output](0bad23b6.png)
(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}:
0.6596581918663791
-2.8083611073063413
-13.617364590734098
sampler() # another random i.c.
3-element Vector{Float64}:
-8.9860464345904
11.107959562212084
-7.715979904448993
and finally call
fs = basins_fractions(mapper, sampler)
Dict{Int64, Float64} with 2 entries:
2 => 0.349
1 => 0.651
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)
![Example block output](e4ca82f2.png)
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)
![Example block output](66b9aec2.png)
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
1
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, OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.OrdinaryDiffEqVerner.Vern9{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.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}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEq.OrdinaryDiffEqVerner.Vern9{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.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}, Vector{SVector{3, Float64}}, Vector{Float64}, Vector{Vector{SVector{3, Float64}}}, Nothing, OrdinaryDiffEq.OrdinaryDiffEqVerner.Vern9ConstantCache, Nothing}, SciMLBase.DEStats, 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}, OrdinaryDiffEq.OrdinaryDiffEqVerner.Vern9ConstantCache, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.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, OrdinaryDiffEq.DefaultInit, Nothing}, Vector{Float64}}, Attractors.BasinsInfo{3, Attractors.RegularGrid{3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, 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}}(2 => 3-dimensional StateSpaceSet{Float64} with 314 points, 1 => 3-dimensional StateSpaceSet{Float64} with 893 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.937125748502994, 3 => 0.06287425149700598), 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}}[Dict(1 => 3-dimensional StateSpaceSet{Float64} with 843 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 746 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 826 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 784 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 736 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 654 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 540 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 768 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 754 points), Dict(1 => 3-dimensional StateSpaceSet{Float64} with 758 points) … Dict(2 => 3-dimensional StateSpaceSet{Float64} with 501 points, 3 => 3-dimensional StateSpaceSet{Float64} with 447 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 499 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 466 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 463 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 427 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 364 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 340 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 326 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 336 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 348 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}} with 1 entry:
1 => 3-dimensional StateSpaceSet{Float64} with 489 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,
)
![Example block output](c0f51df8.png)
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,
)
![Example block output](b6218b98.png)
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
![Example block output](be7bb423.png)
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.11457746082750221)
Dict(1 => 0.06960907516347392)
Dict(1 => 0.0899495823729251)
Dict(1 => 0.0961200667822121)
Dict(1 => 0.07961090186853706)
Dict(1 => 0.04210473114090952)
Dict(1 => 0.00019841781930877444)
Dict(1 => 0.092435145085838)
Dict(1 => 0.08511765729253534)
Dict(1 => 0.07863522548623969)
⋮
Dict(2 => -4.034668545983373e-5)
Dict(2 => 0.00016751728947577275)
Dict(2 => 0.00014981730511878652)
Dict(2 => 0.00029320515854832095)
Dict(2 => 0.0002371189379029226)
Dict(2 => 4.4418025507002096e-5)
Dict(2 => 9.74838967161415e-5)
Dict(2 => 8.189699640859993e-5)
Dict(2 => 4.436708377727187e-5)
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
![Example block output](bb58d465.png)
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
![Example block output](9887b9b2.png)
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.732, 1 => 0.268), Dict(2 => 0.7115768463073853, 1 => 0.2884231536926148), Dict(2 => 0.7345309381237525, 1 => 0.2654690618762475), Dict(2 => 0.7534930139720559, 1 => 0.2465069860279441), Dict(2 => 0.7275449101796407, 1 => 0.27245508982035926), Dict(2 => 0.7534930139720559, 1 => 0.2465069860279441), Dict(2 => 0.7395209580838323, 1 => 0.26047904191616766), Dict(2 => 0.7554890219560878, 1 => 0.24451097804391217), Dict(2 => 0.7345309381237525, 1 => 0.2654690618762475), Dict(2 => 0.7564870259481038, 1 => 0.2435129740518962) … Dict(2 => 0.7155688622754491, 1 => 0.2844311377245509), Dict(2 => 0.717564870259481, 1 => 0.282435129740519), Dict(2 => 0.7514970059880239, 1 => 0.24850299401197604), Dict(2 => 0.7345309381237525, 1 => 0.2654690618762475), Dict(2 => 0.7145708582834331, 1 => 0.28542914171656686), Dict(2 => 0.7335329341317365, 1 => 0.26646706586826346), Dict(2 => 0.7285429141716567, 1 => 0.2714570858283433), Dict(2 => 0.7265469061876247, 1 => 0.27345309381237526), Dict(2 => 0.7375249500998003, 1 => 0.2624750499001996), Dict(2 => 0.7654690618762475, 1 => 0.2345309381237525)], Dict{Int64, StateSpaceSet{3, Float64}}[Dict(2 => 3-dimensional StateSpaceSet{Float64} with 492 points, 1 => 3-dimensional StateSpaceSet{Float64} with 328 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 502 points, 1 => 3-dimensional StateSpaceSet{Float64} with 323 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 556 points, 1 => 3-dimensional StateSpaceSet{Float64} with 324 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 635 points, 1 => 3-dimensional StateSpaceSet{Float64} with 336 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 635 points, 1 => 3-dimensional StateSpaceSet{Float64} with 342 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 658 points, 1 => 3-dimensional StateSpaceSet{Float64} with 323 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 677 points, 1 => 3-dimensional StateSpaceSet{Float64} with 316 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 657 points, 1 => 3-dimensional StateSpaceSet{Float64} with 320 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 627 points, 1 => 3-dimensional StateSpaceSet{Float64} with 335 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 710 points, 1 => 3-dimensional StateSpaceSet{Float64} with 327 points) … Dict(2 => 3-dimensional StateSpaceSet{Float64} with 684 points, 1 => 3-dimensional StateSpaceSet{Float64} with 313 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 632 points, 1 => 3-dimensional StateSpaceSet{Float64} with 322 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 651 points, 1 => 3-dimensional StateSpaceSet{Float64} with 341 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 597 points, 1 => 3-dimensional StateSpaceSet{Float64} with 321 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 546 points, 1 => 3-dimensional StateSpaceSet{Float64} with 327 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 530 points, 1 => 3-dimensional StateSpaceSet{Float64} with 313 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 488 points, 1 => 3-dimensional StateSpaceSet{Float64} with 317 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 494 points, 1 => 3-dimensional StateSpaceSet{Float64} with 301 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 524 points, 1 => 3-dimensional StateSpaceSet{Float64} with 327 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 485 points, 1 => 3-dimensional StateSpaceSet{Float64} with 330 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.