Convenience function continuation_series (#163)

* fix ignoring of accesses in attractors plot

* allow 3D access to make 3D axis

* make the continuation series function general usfeful

* improve docstring of globa lcontinuation

* rename file to global contiinuation

* make public the continuation series

* mention the new function in the docs

* change recurrences keyword to stop at ΔTr

* remove bifurcation kit limit in compilation

* make convergence time of proximity mapper more safe

* better handling of converge time

* better docstring

* base ref value

* limit both BFKit and KrylovKit

* fix missing attractors

* fix missing ds reference

* even more limiting of versions

* Update CHANGELOG.md

CHANGELOG.md

@@ -1,3 +1,11 @@
+# v1.23
+
+- New convenience function `continuation_series` that was already used internally
+  when plotting, but now also made public.
+- The keyword `force_non_adaptive` of `AttractorsViaRecurrences` is deprecated in favor
+  of the keyword `stop_at_Δt` which has the same practical impact with simpler implementation.
+- Keyword `stop_at_Δt` added to `AttractorsViaProximity`.
+
 # v1.22
 
 - The function `plot_continuation_curves` has been changed and now plots using

Project.toml

@@ -2,7 +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.22.1"
+version = "1.23.0"
 
 
 [deps]

docs/Project.toml

@@ -13,6 +13,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 PredefinedDynamicalSystems = "31e2f376-db9e-427a-b76e-a14f56347a14"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 
 [compat]
-BifurcationKit = "0.3"
+BifurcationKit = "0.3.6"
+KrylovKit = "0.8.1"

docs/src/api.md

@@ -126,6 +126,8 @@ MFSBruteForce
 
 ```@docs
 global_continuation
+GlobalContinuationAlgorithm
+continuation_series
 ```
 
 ### General seeding-based continuation

docs/src/tutorial.jl

@@ -79,14 +79,15 @@ Pkg.status(["Attractors", "CairoMakie", "OrdinaryDiffEq"])
 # ## 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))
+# angles = range(0, 2π; length = 101)
+# pcurve = params.(angles)
 # 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
+# 	fractions_cont, attractors_cont, A -> minimum(A[:, 1]), angles; add_legend = false
 # )
 # ```
 
@@ -284,9 +285,10 @@ plot_attractors(attractors3)
 # 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](https://github.com/bifurcationkit/BifurcationKit.jl).
+# These software perform **local continuation**.
 # Here we offer a completely different kind of continuation called **global continuation**.
 
-# The traditional continuation analysis continues the curves of individual _fixed
+# The local continuation 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
@@ -298,23 +300,23 @@ plot_attractors(attractors3)
 # 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
+# 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 [Datseris2023](@cite).
 # See also the [comparison page](@ref bfkit_comparison) 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:
+# To perform a global continuation is surprisingly simple. First, we decide what parameter,
+# and what range of that parameter, to continue over:
 
 prange = 4.5:0.01:6
 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 geenral [`AttractorSeedContinueMatch`](@ref).
+# which is the general [`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).
@@ -337,8 +339,10 @@ fractions_cont, attractors_cont = global_continuation(
 
 attractors_cont[34]
 
-# There is a fantastic convenience function for animating
-# the attractors evolution, that utilizes things we have
+# If you want to transform the output to the alternative format of
+# a dictionary of vectors, use [`continuation_series`](@ref).
+# You typically don't have to though, because there is a fantastic convenience function
+# for animating the attractors evolution, that utilizes things we have
 # already defined:
 
 animate_attractors_continuation(
@@ -352,7 +356,7 @@ animate_attractors_continuation(
 # ```
 
 # 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
+# It can be 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!)
 

ext/plotting.jl

@@ -33,10 +33,11 @@ end
 ##########################################################################################
 # Attractors
 ##########################################################################################
-function Attractors.plot_attractors(a; kw...)
+function Attractors.plot_attractors(a; access = SVector(1, 2), kw...)
     fig = Figure()
-    ax = Axis(fig[1,1])
-    plot_attractors!(ax, a; kw...)
+    AX = length(access) == 2 ? Axis : Axis3
+    ax = AX(fig[1,1])
+    plot_attractors!(ax, a; access, kw...)
     return fig
 end
 
@@ -52,8 +53,7 @@ function Attractors.plot_attractors!(ax, attractors;
     for k in ukeys
         k ∉ keys(attractors) && continue
         A = attractors[k]
-        x, y = columns(A[:, access])
-        scatter!(ax, x, y;
+        scatter!(ax, A[:, access];
             color = (colors[k], 0.9), markersize = 20,
             marker = markers[k],
             label = "$(labels[k])",
@@ -237,14 +237,14 @@ function Attractors.plot_basins_curves!(ax, fractions_cont, prange = 1:length(fr
     if style == :band
         # transform to cumulative sum
         for j in 2:length(bands)
-            bands[j] .+= bands[j-1]
+            bands[k[j]] .+= bands[k[j-1]]
         end
         for (j, k) in enumerate(ukeys)
             if j == 1
-                l, u = 0, bands[j]
+                l, u = 0, bands[k]
                 l = fill(0f0, length(u))
             else
-                l, u = bands[j-1], bands[j]
+                l, u = bands[k-1], bands[k]
             end
             band!(ax, prange, l, u;
                 color = colors[k], label = "$(labels[k])", series_kwargs...
@@ -255,8 +255,8 @@ function Attractors.plot_basins_curves!(ax, fractions_cont, prange = 1:length(fr
         end
         ylims!(ax, 0, 1)
     elseif style == :lines
-        for (j, k) in enumerate(ukeys)
-            scatterlines!(ax, prange, bands[j];
+        for k in ukeys
+            scatterlines!(ax, prange, bands[k];
                 color = colors[k], label = "$(labels[k])", marker = markers[k],
                 markersize = 10, linewidth = 3, series_kwargs...
             )
@@ -270,16 +270,6 @@ function Attractors.plot_basins_curves!(ax, fractions_cont, prange = 1:length(fr
     return
 end
 
-function continuation_series(continuation_info, defval, ukeys = unique_keys(continuation_info))
-    bands = [zeros(length(continuation_info)) for _ in ukeys]
-    for i in eachindex(continuation_info)
-        for (j, k) in enumerate(ukeys)
-            bands[j][i] = get(continuation_info[i], k, defval)
-        end
-    end
-    return bands
-end
-
 function Attractors.plot_attractors_curves(attractors_cont, attractor_to_real, prange = 1:length(attractors_cont); kwargs...)
     fig = Figure()
     ax = Axis(fig[1,1])

src/Attractors.jl

@@ -16,7 +16,7 @@ using Reexport
 include("dict_utils.jl")
 include("mapping/attractor_mapping.jl")
 include("basins/basins.jl")
-include("continuation/basins_fractions_continuation_api.jl")
+include("continuation/global_continuation_api.jl")
 include("matching/matching_interface.jl")
 include("boundaries/edgetracking.jl")
 include("deprecated.jl")

src/continuation/basins_fractions_continuation_api.jl → src/continuation/global_continuation_api.jl

@@ -1,4 +1,4 @@
-export global_continuation, GlobalContinuationAlgorithm
+export global_continuation, GlobalContinuationAlgorithm, continuation_series
 
 """
     GlobalContinuationAlgorithm
@@ -17,7 +17,29 @@ abstract type GlobalContinuationAlgorithm end
     global_continuation(gca::GlobalContinuationAlgorithm, pcurve, ics; kwargs...)
 
 Find and continue attractors (or representations of attractors)
-and the fractions of their basins of attraction across a parameter range.
+and the fractions of their basins of attraction across a parameter range `pcurve`
+by sampling given initial conditions `ics` according to algorithm `gca`.
+
+Return:
+
+1. `fractions_cont::Vector{Dict{Int, Float64}}`. The fractions of basins of attraction.
+   `fractions_cont[i]` is a dictionary mapping attractor IDs to their basin fraction
+   at the `i`-th parameter combination.
+2. `attractors_cont::Vector{Dict{Int, <:Any}}`. The continued attractors.
+   `attractors_cont[i]` is a dictionary mapping attractor ID to the
+   attractor set at the `i`-th parameter combination.
+
+See the function [`continuation_series`](@ref) if you wish to transform the output
+to an alternative format.
+
+## Keyword arguments
+
+- `show_progress = true`: display a progress bar of the computation.
+- `samples_per_parameter = 100`: amount of initial conditions sampled at each parameter
+  combination from `ics` if `ics` is a function instead of set initial conditions.
+
+## Description
+
 `global_continuation` is the central function of the framework for global stability analysis
 illustrated in [Datseris2023](@cite).
 
@@ -46,28 +68,33 @@ Possible subtypes of `GlobalContinuationAlgorithm` are:
 
 - [`AttractorSeedContinueMatch`](@ref)
 - [`FeaturizeGroupAcrossParameter`](@ref)
-
-## Return
-
-1. `fractions_cont::Vector{Dict{Int, Float64}}`. The fractions of basins of attraction.
-   `fractions_cont[i]` is a dictionary mapping attractor IDs to their basin fraction
-   at the `i`-th parameter combination.
-2. `attractors_cont::Vector{Dict{Int, <:Any}}`. The continued attractors.
-   `attractors_cont[i]` is a dictionary mapping attractor ID to the
-   attractor set at the `i`-th parameter combination.
-
-## Keyword arguments
-
-- `show_progress = true`: display a progress bar of the computation.
-- `samples_per_parameter = 100`: amount of initial conditions sampled at each parameter
-  combination from `ics` if `ics` is a function instead of set initial conditions.
 """
 function global_continuation(alg::GlobalContinuationAlgorithm, prange::AbstractVector, pidx, sampler; kw...)
     # everything is propagated to the curve setting
     pcurve = [[pidx => p] for p in prange]
     return global_continuation(alg, pcurve, sampler; kw...)
 end
 
+
+"""
+    continuation_series(continuation_info, defval = NaN)
+
+Transform a continuation quantity (a vector of dictionaries, each dictionary
+mapping attractor IDs to real numbers) to a dictionary of vectors where
+the `k` entry is the series of the continuation quantity corresponding to
+attractor with ID `k`. `defval` denotes the value to assign in the series
+if an attractor does not exist at this particular series index.
+"""
+function continuation_series(continuation_info::AbstractVector{<:AbstractDict}, defval = NaN, ukeys = unique_keys(continuation_info))
+    bands = Dict(k => zeros(length(continuation_info)) for k in ukeys)
+    for i in eachindex(continuation_info)
+        for k in ukeys
+            bands[k][i] = get(continuation_info[i], k, defval)
+        end
+    end
+    return bands
+end
+
 include("continuation_ascm_generic.jl")
 include("continuation_recurrences.jl")
 include("continuation_grouping.jl")

src/dict_utils.jl

@@ -114,4 +114,4 @@ function next_free_id(keys₊, keys₋)
     s = setdiff(keys₊, keys₋)
     nextid = isempty(s) ? maximum(keys₋) + 1 : minimum(s)
     return nextid
-end
+end