:erge branch 'main' into stagger_step

CHANGELOG.md

@@ -1,3 +1,9 @@
+# v1.19
+
+- Global continuation can now be performed across any arbitrary curve
+  in parameter space. This is something completely novel
+  and we'll likely work on a paper on this!
+
 # v1.18
 
 This is a big release, with (hopefully) nothing breaking, but lots of deprecations!

Project.toml

@@ -2,8 +2,7 @@ name = "Attractors"
 uuid = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 authors = ["George Datseris <[email protected]>", "Kalel Rossi", "Alexandre Wagemakers"]
 repo = "https://github.com/JuliaDynamics/Attractors.jl.git"
-version = "1.18.5"
-
+version = "1.19.4"
 
 [deps]
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
@@ -33,7 +32,7 @@ Clustering = "0.15"
 ComplexityMeasures = "2.3, 3"
 Distances = "0.7, 0.8, 0.9, 0.10"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
-DynamicalSystemsBase = "3.0.3"
+DynamicalSystemsBase = "3.9.1"
 Makie = "≥ 0.19"
 Neighborhood = "0.2.2"
 Optim = "1.7"

docs/src/bfkit_comparison.jl

@@ -210,7 +210,7 @@ plot_attractors_curves(
 
 # ## Discussion and comparison
 
-# On the one hand, Attractors.jl found not only a single limit cycle,
+# Attractors.jl found not only a single limit cycle,
 # but also all system attractors, including chaotic ones.
 # It didn't require any initial guess regarding the limit cycle or its period,
 # but only a state space box that may contain attractors.
@@ -224,9 +224,17 @@ plot_attractors_curves(
 # And finally, Attractors.jl estimates a more general nonlocal measure of stability,
 # in the sense that if a set is nonlocally stable, it is guaranteed to be locally stable,
 # however the other way around isn't guaranteed.
-
-# On the other hand, traditional local continuation can track the unstable
-# branches, and automatically detect and label local bifurcations
+# Moreover, due to the orthogonality of finding and matching attractors,
+# as well as finding _all_ attractors, the global continuation of Attractors.jl
+# can continue along arbitrary user-defined curves in parameter space; not just
+# along a single parameter axis. This is possible because it is completely fine
+# for some attractors to stop existing during the global continuation,
+# while local continuation stops when attractors (and their unstable version)
+# stop existing.
+
+# Traditional local continuation can track the unstable
+# branches, and automatically detect and label local bifurcations,
+# both of which are not possible in Attractors.jl
 # (note we didn't bother to plot any of the detected bifurcations - if any are found).
 # In our experience having the local bifurcations is always useful.
 # Now, whether the unstable branches are of a limit cycle are useful or not,
@@ -248,6 +256,7 @@ plot_attractors_curves(
 # Using a slightly incorrect initial period guess of 20.0 instead of 19.0 also fails.
 # We imagine that this sensitivity would apply also to some other of the
 # several meta-parameters that enter a traditional continuation routine,
+# for example the thresholds and accuracy parameters related to Newton convergence,
 # but we didn't check further. This is exactly what we were alluding to
 # in the comparison we did in [Datseris2023](@cite), that traditional
 # local continuation "requires expertise and constant interventions".

docs/src/examples.md

@@ -2,7 +2,7 @@
 
 Note that the examples utilize some convenience plotting functions offered by Attractors.jl which come into scope when using `Makie` (or any of its backends such as `CairoMakie`), see the [visualization utilities](@ref) for more.
 
-## Newton's fractal (basins of 2D map)
+## Newton's fractal (basins of a 2D map)
 
 ```@example MAIN
 using Attractors

docs/src/index.md

@@ -10,12 +10,8 @@ using CairoMakie, Attractors
 
 ## Latest news
 
-- Brand new [tutorial](@ref tutorial) for the whole repo!
-- `continuation` has been renamed to `global_continuation`
-- Brand new continuation algorithm `AttractorSeedContinueMatch`
-  that generalizes existing infrastructure and makes Attractors.jl
-  more composable and more extendable.
-- New "matchers" and extendable interface for them.
+- Global continuation can now be performed across any arbitrary curve
+  in parameter space.
 - See the CHANGELOG.md (at the GitHub repo) for more!
 
 ## Getting started

docs/src/tutorial.jl

@@ -1,5 +1,13 @@
 # # [Attractors.jl Tutorial](@id tutorial)
 
+# ```@raw html
+# <video width="auto" controls loop>
+# <source src="../attracont.mp4" type="video/mp4">
+# </video>
+# ```
+
+#nb @doc Attractors
+
 # [`Attractors`](@ref) 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
@@ -13,6 +21,18 @@
 # but we won't cover anything else in this introductory tutorial.
 # See the [examples](@ref examples) page instead.
 
+
+# ### Package versions used
+
+import Pkg
+
+#nb # Activate an environment in the folder containing the notebook
+#nb Pkg.activate(dirname(@__DIR__))
+#nb Pkg.add(["DynamicalSystems", "CairoMakie", "GLMakie", "OrdinaryDiffEq", "BenchmarkTools"])
+
+Pkg.status(["Attractors", "CairoMakie", "OrdinaryDiffEq"])
+
+
 # ## Tutorial - copy-pasteable version
 
 # _Gotta go fast!_
@@ -294,7 +314,7 @@ pidx = 1 # index of the parameter
 # Then, we may call the [`global_continuation`](@ref) function.
 # We have to provide a continuation algorithm, which itself references an [`AttractorMapper`](@ref).
 # In this example we will re-use the `mapper` to create the "flagship product" of Attractors.jl
-# which is the generic [`AttractorSeedContinueMatch`](@ref).
+# which is the geenral [`AttractorSeedContinueMatch`](@ref).
 # 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).
@@ -486,6 +506,45 @@ 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))
+
+# 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
+)
+
+# and animate the result
+animate_attractors_continuation(
+    ds, attractors_cont, fractions_cont, pcurve;
+    savename = "curvecont.mp4"
+);
+
+# ```@raw html
+# <video width="auto" controls loop>
+# <source src="../curvecont.mp4" type="video/mp4">
+# </video>
+# ```
+
 # ## Conclusion and comparison with traditional local continuation
 
 # We've reached the end of the tutorial! Some aspects we haven't highlighted is

ext/plotting.jl

@@ -395,7 +395,13 @@ end
 # Videos
 ##########################################################################################
 function Attractors.animate_attractors_continuation(
-        ds::DynamicalSystem, attractors_cont, fractions_cont, prange, pidx;
+    ds::DynamicalSystem, attractors_cont, fractions_cont, prange, pidx; kw...)
+    pcurve = [[pidx => p] for p in prange]
+    return animate_attractors_continuation(ds, attractors_cont, fractions_cont, pcurve; kw...)
+end
+
+function Attractors.animate_attractors_continuation(
+        ds::DynamicalSystem, attractors_cont, fractions_cont, pcurve;
         savename = "attracont.mp4", access = SVector(1, 2),
         limits = auto_attractor_lims(attractors_cont, access),
         framerate = 4, markersize = 10,
@@ -406,10 +412,11 @@ function Attractors.animate_attractors_continuation(
         T = 100,
         figure = NamedTuple(), axis = NamedTuple(), fracaxis = NamedTuple(),
         legend = NamedTuple(),
+        add_legend = length(ukeys) ≤ 6
     )
     length(access) ≠ 2 && error("Need two indices to select two dimensions of `ds`.")
     K = length(ukeys)
-    fig = Figure(figure...)
+    fig = Figure(; figure...)
     ax = Axis(fig[1,1]; limits, axis...)
     fracax = Axis(fig[1,2]; width = 50, limits = (0,1,0,1), ylabel = "fractions",
         yaxisposition = :right, fracaxis...
@@ -423,19 +430,21 @@ function Attractors.animate_attractors_continuation(
     for k in ukeys
         plotf!(ax, att_obs[k]; color = (colors[k], 0.75), label = "$k", markersize, marker = markers[k])
     end
-    axislegend(ax; legend...)
+    if add_legend
+        axislegend(ax; legend...)
+    end
 
     # setup fractions axis
     heights = Observable(fill(0.1, K))
     barcolors = [colors[k] for k in ukeys]
     barplot!(fracax, fill(0.5, K), heights; width = 1, gap = 0, stack=1:K, color = barcolors)
 
-    record(fig, savename, eachindex(prange); framerate) do i
-        p = prange[i]
-        ax.title = "p = $p"
+    record(fig, savename, eachindex(pcurve); framerate) do i
+        p = pcurve[i]
+        ax.title = "p: $p" # TODO: Add compat printing here.
         attractors = attractors_cont[i]
         fractions = fractions_cont[i]
-        set_parameter!(ds, pidx, p)
+        set_parameters!(ds, p)
         heights[] = [get(fractions, k, 0) for k in ukeys]
 
         for (k, att) in attractors

src/basins/basins_utilities.jl

@@ -1,3 +1,46 @@
+# It works for all mappers that define a `basins_fractions` method.
+"""
+    basins_of_attraction(mapper::AttractorMapper, grid::Tuple) → basins, attractors
+
+Compute the full basins of attraction as identified by the given `mapper`,
+which includes a reference to a [`DynamicalSystem`](@ref) and return them
+along with (perhaps approximated) found attractors.
+
+`grid` is a tuple of ranges defining the grid of initial conditions that partition
+the state space into boxes with size the step size of each range.
+For example, `grid = (xg, yg)` where `xg = yg = range(-5, 5; length = 100)`.
+The grid has to be the same dimensionality as the state space expected by the
+integrator/system used in `mapper`. E.g., a [`ProjectedDynamicalSystem`](@ref)
+could be used for lower dimensional projections, etc. A special case here is
+a [`PoincareMap`](@ref) with `plane` being `Tuple{Int, <: Real}`. In this special
+scenario the grid can be one dimension smaller than the state space, in which case
+the partitioning happens directly on the hyperplane the Poincaré map operates on.
+
+`basins_of_attraction` function is a convenience 5-lines-of-code wrapper which uses the
+`labels` returned by [`basins_fractions`](@ref) and simply assigns them to a full array
+corresponding to the state space partitioning indicated by `grid`.
+
+See also [`convergence_and_basins_of_attraction`](@ref).
+"""
+function basins_of_attraction(mapper::AttractorMapper, grid::Tuple; kwargs...)
+    basins = zeros(Int32, map(length, grid))
+    I = CartesianIndices(basins)
+    A = StateSpaceSet([generate_ic_on_grid(grid, i) for i in vec(I)])
+    fs, labels = basins_fractions(mapper, A; kwargs...)
+    attractors = extract_attractors(mapper)
+    vec(basins) .= vec(labels)
+    return basins, attractors
+end
+
+# Type-stable generation of an initial condition given a grid array index
+@generated function generate_ic_on_grid(grid::NTuple{B, T}, ind) where {B, T}
+    gens = [:(grid[$k][ind[$k]]) for k=1:B]
+    quote
+        Base.@_inline_meta
+        @inbounds return SVector{$B, Float64}($(gens...))
+    end
+end
+
 """
     basins_fractions(basins::AbstractArray [,ids]) → fs::Dict