New generic continuation type and new matcher interface (#132)

* add new plot_continuation_curves visualization

* move recurrences animation to hidden page

* add enriching continuaton example

* reset comment

* clarify plotting utility

* remove invalid zorder attribute

* add ChaosTool to project docs

* WIP on tutorial

* write a tutorial!

* make a dedicated tutorial page

* simpler getting started statement

* put literate in project

* fix build errors

* make autolimits a bit larger

* add new generalized continuation type

* wip im tired and stop here for now

* wip more wip

* strcture matcher api around replacement map

* finish docstring of replacement map to be general

* port over the continuation matching function

* reference the retract keys

* port over docstring of MatchSSS

* delete old file

* Update src/continuation/continuation_afam.jl

Co-authored-by: Alexandre Wagemakers <[email protected]>

* Update src/matching/sssdistance.jl

Co-authored-by: Alexandre Wagemakers <[email protected]>

* improve docstring

* finish up things (maybe?)

* rename to match sequentially

* deprecate info extraction

* fix tests and name exports

* implement retracted

* correct code so that tutorial works

* update doc links

* finish basin overlap (TODO: deprecation and porting)

* deprecate old basins overlap function

* finish some deprecations

* finish basin overlap code and tests

* fix aggregation tests and remove obsolete test

* global rename `fractions_curves` to `fractions_cont`

* glboal rename attractors_info to attractors_cont

* push basin_enclosure

* rename to IDMatcher

* do not promote IDMatcher, but keep concrete types for clarity

* rename supertype to `GlobalContinuationAlgorithm`

* small improvement

* deprecate continuation to global_continuation

* deprecate part 2

* compilable module

* rename to AttractorSeedContinueMatch

* more updates

* add changelog

* export reset mapper

* make Featurizing always store the "attractors"

* add changelog

* add test for extract attractors for hte mapper

* add tests for generalized seed continue with featurizing grouping

* improve docstring of new continuation

* finish incredible docstring of new continuation

* reference find and match

* finish extendable matching interface

* more reasonable statement in tutorial

* correct docstring

* address minor comment

* more minor address

* add matching by basin enclosure

* addres minor

* correct ehader for basin enclosure

* add error for currently missing coflowing logic

* fix dumb seeding mistake in generic continuation

* export IDmatcher

* add comment on the logic new to old

* correct changelog

---------

Co-authored-by: Alexandre Wagemakers <[email protected]>

CHANGELOG.md

@@ -1,8 +1,22 @@
 # v1.18
 
+This is a big release, with (hopefully) nothing breaking, but lots of deprecations!
+
+## New stuff
+
+- New central Tutorial for Attractors.jl. It also highlights how to enrich a continuation output with other measures of nonlocal stability.
+- New global continuation algorithm that generalizes RAFM: `AttractorSeedContinueMatch`.
+- There is now an extendable API for "matchers", ways to match state spaces sets across a continuation. See `IDMatcher` for the API.
 - New plotting function `plot_continuation_curves!` to add additional information to the `plot_basins_attractors_curves` type plots.
-- New example in the docs for enriching a continuation output with other examples of nonlocal stability.
-- New exported function `reset_mapper!` to clear all stored information in an `AttractorsViaRecurrences` instance. Used in the aforementioned new example.
+- New exported function `reset_mapper!` to clear all stored information in an `AttractorMapper`.
+- `AttractorsViaFeaturizing` now always stores the attractors and implements `extract_attractors`.
+
+
+## Deprecations and renaming
+
+- Function `continuation` has been deprecated for `global_continuation` in preparation for a future where both local/linear/tradiational continuation as well as our "attractors and basins continuation" are both provided by DynamicalSystems.jl.
+- `match_continuation!` has been deprecated for `match_sequentially!`.
+- Option `par_weight` is deprecated in `FeaturizeGroupAcrossParameter`. Part of the developer team (`@Datseris`, `@KalelR`) discussed this an concluded that `par_weight` doesn't make much scientific sense to include. Since it obfuscates the code and the documentation, it is no longer documented but still exported. It will be unavailable in the next breaking release.
 
 # v1.17
 

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.17.1"
+version = "1.18.0"
 
 
 [deps]

docs/src/api.md

@@ -7,9 +7,9 @@ Attractors.jl defines a generic interface for finding attractors of dynamical sy
 
 ```@docs
 AttractorMapper
+extract_attractors
 ```
 
-
 ### Recurrences
 
 ```@docs
@@ -62,13 +62,13 @@ extract_features
 
 
 ## Basins of attraction
+
 Calculating basins of attraction, or their state space fractions, can be done with the functions:
 - [`basins_fractions`](@ref)
 - [`basins_of_attraction`](@ref).
 
 ```@docs
 basins_fractions
-extract_attractors
 basins_of_attraction
 statespace_sampler
 ```
@@ -122,10 +122,16 @@ MFSBlackBoxOptim
 MFSBruteForce
 ```
 
-## Continuation API
+## Global continuation
 
 ```@docs
-continuation
+global_continuation
+```
+
+### General seeding-based continuation
+
+```@docs
+AttractorSeedContinueMatch
 ```
 
 ### Recurrences continuation
@@ -146,29 +152,34 @@ aggregate_attractor_fractions
 FeaturizeGroupAcrossParameter
 ```
 
-## Matching continuation output
+## Matching attractors
+
+Matching attractors follow an extendable interface based on [`IDMatcher`](@ref).
+The available matchers are:
+
+```@docs
+MatchBySSSetDistance
+MatchByBasinOverlap
+```
+
+### Matching interface
 
-One of the most novel, and most powerful features of Attractors.jl is its
-ability to "match" the continuation result: i.e., functionality that decides
-which attractor gets which label at each parameter value, so that there is a
-"continuity" across the parameter axis.
+```@docs
+IDMatcher
+replacement_map
+replacement_map!
+match_sequentially!
+Attractors._match_attractors
+```
 
-The matching process is entirely orthogonal to the continuation, something
-completely novel in continuation software. This means, that if you don't like
-the way the matching worked in the first time you estimated the attractors and their
-basins, you don't have to re-compute them! You can simply relabel them using the
-[`match_continuation!`](@ref) function!
+### Low-level distance functions
 
 ```@docs
-match_statespacesets!
 Centroid
 Hausdorff
 StrictlyMinimumDistance
-replacement_map
 set_distance
 setsofsets_distances
-match_continuation!
-match_basins_ids!
 ```
 
 ## Visualization utilities
@@ -185,6 +196,7 @@ heatmap_basins_attractors!(ax, grid, basins, attractors; kwargs...)
 ```
 
 ### [Common plotting keywords](@id common_plot_kwargs)
+
 Common keywords for plotting functions in Attractors.jl are:
 
 - `ukeys`: the basin ids (unique keys, vector of integers) to use. By default all existing keys are used.

docs/src/examples.md

@@ -409,7 +409,7 @@ fig
 
 ## Basin fractions continuation in the magnetic pendulum
 
-Perhaps the simplest application of [`continuation`](@ref) is to produce a plot of how the fractions of attractors change as we continuously change the parameter we changed above to calculate tipping probabilities.
+Perhaps the simplest application of [`global_continuation`](@ref) is to produce a plot of how the fractions of attractors change as we continuously change the parameter we changed above to calculate tipping probabilities.
 
 ### Computing the fractions
 
@@ -434,12 +434,12 @@ sampler, = statespace_sampler(region, 1234)
 # continue attractors and basins:
 # `Inf` threshold fits here, as attractors move smoothly in parameter space
 rsc = RecurrencesFindAndMatch(mapper; threshold = Inf)
-fractions_curves, attractors_info = continuation(
+fractions_cont, attractors_cont = global_continuation(
     rsc, prange, pidx, sampler;
     show_progress = false, samples_per_parameter = 100
 )
 # Show some characteristic fractions:
-fractions_curves[[1, 50, 101]]
+fractions_cont[[1, 50, 101]]
 ```
 
 
@@ -448,7 +448,7 @@ We visualize them using a predefined function that you can find in `docs/basins_
 
 ```@example MAIN
 # careful; `prange` isn't a vector of reals!
-plot_basins_curves(fractions_curves, γγ)
+plot_basins_curves(fractions_cont, γγ)
 ```
 
 
@@ -467,7 +467,7 @@ function real_number_repr(attractor)
 end
 
 for (i, γ) in enumerate(γγ)
-    for (k, attractor) in attractors_info[i]
+    for (k, attractor) in attractors_cont[i]
         scatter!(ax, γ, real_number_repr(attractor); color = Cycled(k))
     end
 end
@@ -504,12 +504,12 @@ sampler, = statespace_sampler(grid, 1234)
 # initialize mapper
 mapper = AttractorsViaRecurrences(ds, grid; recurrences_kwargs...)
 # perform continuation of attractors and their basins
-continuation = RecurrencesFindAndMatch(mapper; threshold = Inf)
-fractions_curves, attractors_info = continuation(
-    continuation, prange, pidx, sampler;
+alg = RecurrencesFindAndMatch(mapper; threshold = Inf)
+fractions_cont, attractors_cont = global_continuation(
+    alg, prange, pidx, sampler;
     show_progress = true, samples_per_parameter
-);
-plot_basins_curves(fractions_curves, prange; separatorwidth = 1)
+)
+plot_basins_curves(fractions_cont, prange; separatorwidth = 1)
 ```
 
 ![](https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/videos/attractors/multispecies_competition_fractions.png)
@@ -537,7 +537,7 @@ isextinct(A, idx = unitidxs) = all(a -> a <= 1e-2, A[:, idx])
 groupingconfig = GroupViaClustering(; min_neighbors=1, optimal_radius_method=0.5)
 
 aggregated_fractions, aggregated_info = aggregate_attractor_fractions(
-    fractions_curves, attractors_info, featurizer, groupingconfig
+    fractions_cont, attractors_cont, featurizer, groupingconfig
 )
 
 plot_basins_curves(aggregated_fractions, prange;
@@ -620,10 +620,10 @@ featurizer(a, t) = a[end]
 clusterspecs = GroupViaClustering(optimal_radius_method = "silhouettes", max_used_features = 200)
 mapper = AttractorsViaFeaturizing(ds, featurizer, clusterspecs; T = 20, threaded = true)
 gap = FeaturizeGroupAcrossParameter(mapper; par_weight = 1.0)
-fractions_curves, clusters_info = continuation(
+fractions_cont, clusters_info = global_continuation(
     gap, rrange, ridx, sampler; show_progress = false
 )
-fractions_curves
+fractions_cont
 ```
 
 Looking at the information of the "attractors" (here the clusters of the grouping procedure) does not make it clear which label corresponds to which kind of attractor, but we can look at the:
@@ -647,9 +647,9 @@ using ComplexityMeasures: ValueHistogram, FixedRectangularBinning, probabilities
 # you decide the binning for the histogram, but for a valid estimation of
 # distances, all histograms must have exactly the same bins, and hence be
 # computed with fixed ranges, i.e., using the `FixedRectangularBinning`
-const binning = FixedRectangularBinning(range(-5, 5; length = 11))
 
 function histogram_featurizer(A, t)
+    binning = FixedRectangularBinning(range(-5, 5; length = 11))
     ms = mean.(columns(A)) # vector of mean of each variable
     p = probabilities(ValueHistogram(binning), ms) # this is the histogram
     return vec(p) # because Distances.jl doesn't know `Probabilities`

docs/src/tutorial.jl

@@ -30,7 +30,7 @@
 # ```
 
 # which we define in code as
-using Attractors
+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)
@@ -78,7 +78,7 @@ mapper = AttractorsViaRecurrences(ds, grid;
     consecutive_lost_steps = 100,
 )
 
-# This `mapper` can map any initial condition `u` to the corresponding
+# This `mapper` can map any initial condition to the corresponding
 # attractor ID, for example
 
 mapper([-4.0, 5, 0])
@@ -87,7 +87,7 @@ mapper([-4.0, 5, 0])
 
 mapper([4.0, 2, 0])
 
-# the fact that these two different conditions got assigned different IDs means
+# 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
@@ -151,18 +151,23 @@ fs = basins_fractions(mapper, sampler)
 # [`AttractorMapper`](@ref) defines an extendable interface and can be enriched
 # with other methods in the future!
 
-# ## Continuation
+# ## 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](https://github.com/bifurcationkit/BifurcationKit.jl). Here we offer a completely different kind of continuation called **attractors & basins 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](https://github.com/bifurcationkit/BifurcationKit.jl).
+# Here we offer a completely different kind of continuation called **global continuation**.
 
 # A direct comparison of the two approaches is not truly possible, because they do different things.
-# The traditional linearized continuation analysis continues the curves of individual fixed
-# points across the joint state-parameter space. The attractor and basins continuation first
-# finds all attractors at all parameter values and then _matches_ appropriately similar
-# attractors across different parameters, giving the illusion of continuing them individually.
-
-# This is a fundamental difference. With our approach, one finds all attractors
-# (or almost all, for insufficiently dense sampling). And because all attractors are simultaneously
+# The traditional continuation analysis continues the curves of individual fixed
+# points across the joint state-parameter space and tracks their _local (linear) stability_.
+# The global continuation in Attractors.jl finds all attractors, including chaotic ones,
+# in the whole of the state space (that it searches in), and continues all of these attractors
+# concurrently along a parameter axis.
+# Additionally, this 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 [Datseris2023](@cite).
@@ -173,39 +178,37 @@ fs = basins_fractions(mapper, sampler)
 prange = 4.7:0.02:6
 pidx = 1 # index of the parameter
 
-# Then, we may call the [`continuation`](@ref) function.
+# 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 a [`RecurrencesFindAndMatch`](@ref) continuation algorithm.
 # This algorithm uses the `mapper` to find all attractors at each parameter value.
 # Then, it performs a "matching" step, ensuring a "continuity" of the attractor
 # label across the parameter axis. You can read the docstring for more details,
 # as this algorithm is quite sophisticated!
 
-# For now we can use all of its default options which are reliable 99% of the time
+# For now we can use all of its default options which are reliable most of the time
 
 rafm = RecurrencesFindAndMatch(mapper)
 
 # and call
 
-fractions_curves, attractors_info = continuation(
+fractions_cont, attractors_cont = global_continuation(
 	rafm, prange, pidx, sampler; samples_per_parameter = 1_000
 )
 
-attractors_info
-
 # the output is given as two vectors. Each vector is a dictionary
-# mapping attractor IDs to their fractions, or their state space sets, respectively.
+# 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_info[34]
+attractors_cont[34]
 
 # There is a fantastic convenience function for animating
 # the attractors evolution, that utilizes things we have
 # already defined:
 
 animate_attractors_continuation(
-    ds, attractors_info, fractions_curves, prange, pidx;
+    ds, attractors_cont, fractions_cont, prange, pidx;
 );
 
 # ```@raw html
@@ -223,7 +226,7 @@ animate_attractors_continuation(
 # parameter axis. We can do this with the convenience function:
 
 fig = plot_basins_attractors_curves(
-	fractions_curves, attractors_info, A -> minimum(A[:, 1]), prange,
+	fractions_cont, attractors_cont, A -> minimum(A[:, 1]), prange,
 )
 
 # In the top panel are the basin fractions, by default plotted as stacked bars.
@@ -232,15 +235,13 @@ fig = plot_basins_attractors_curves(
 # an attractor into a real number for plotting.
 # We can provide more functions to visualize other aspects of the attractors:
 
-using Statistics: std
-
 a2rs = [
     A -> minimum(A[:, 1]),
     A -> log(length(A)), # proxy for "complexity"
 ]
 
 fig = plot_basins_attractors_curves(
-	fractions_curves, attractors_info, a2rs, prange; add_legend = false
+	fractions_cont, attractors_cont, a2rs, prange; add_legend = false
 )
 
 ax1, ax2 = content.((fig[2,1], fig[3,1]))
@@ -265,7 +266,7 @@ 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_info[i]
+    attractors = attractors_cont[i]
     Dict(k => lyapunov(ds, 2000.0; u0 = A[1]) for (k, A) in attractors)
 end
 
@@ -300,7 +301,7 @@ mfss = map(enumerate(prange)) do (i, p)
     ## 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_info[i]
+    attractors = attractors_cont[i]
     if length(attractors) == 1
         return Dict(k => NaN for (k, A) in attractors)
     end