Edgetracking algorithm

.gitignore

@@ -3,4 +3,5 @@ Manifest.toml
 *.scss
 *.css
 .vscode
-*style.jl
+*style.jl
+.DS_Store

Project.toml

@@ -2,7 +2,7 @@ name = "Attractors"
 uuid = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 authors = ["George Datseris <datseris.george@gmail.com>", "Kalel Rossi", "Alexandre Wagemakers"]
 repo = "https://github.com/JuliaDynamics/Attractors.jl.git"
-version = "1.13.6"
+version = "1.14.0"
 
 [deps]
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"

docs/make.jl

@@ -15,7 +15,6 @@ pages = [
     "references.md",
 ]
 
-
 import Downloads
 Downloads.download(
     "https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/build_docs_with_style.jl",
@@ -31,5 +30,5 @@ bib = CitationBibliography(
 )
 
 build_docs_with_style(pages, Attractors, DynamicalSystemsBase, StateSpaceSets;
-    expandfirst = ["index.md"], bib,
+    expandfirst = ["index.md"], bib
 )

docs/refs.bib

@@ -158,4 +158,65 @@ @article{Halekotte2020
   author = {Lukas Halekotte and Ulrike Feudel},
   title = {Minimal fatal shocks in multistable complex networks},
   journal = {Scientific Reports}
+}
+
+@article{Battelino1988,
+  title={Multiple coexisting attractors, basin boundaries and basic sets},
+  author={Battelino, Peter M and Grebogi, Celso and Ott, Edward and Yorke, James A and Yorke, Ellen D},
+  journal={Physica D: Nonlinear Phenomena},
+  volume={32},
+  number={2},
+  pages={296--305},
+  year={1988},
+  publisher={Elsevier}
+}
+
+@article{Skufca2006,
+  title={Edge of chaos in a parallel shear flow},
+  author={Skufca, Joseph D and Yorke, James A and Eckhardt, Bruno},
+  journal={Physical review letters},
+  volume={96},
+  number={17},
+  pages={174101},
+  year={2006},
+  publisher={APS}
+}
+
+@article{Schneider2008,
+  title={Laminar-turbulent boundary in plane Couette flow},
+  author={Schneider, Tobias M and Gibson, John F and Lagha, Maher and De Lillo, Filippo and Eckhardt, Bruno},
+  journal={Physical Review E},
+  volume={78},
+  number={3},
+  pages={037301},
+  year={2008},
+  publisher={APS}
+}
+
+@article{Wagemakers2020,
+  title={The saddle-straddle method to test for Wada basins},
+  author={Wagemakers, Alexandre and Daza, Alvar and Sanju{\'a}n, Miguel AF},
+  journal={Communications in Nonlinear Science and Numerical Simulation},
+  volume={84},
+  pages={105167},
+  year={2020},
+  publisher={Elsevier}
+}
+
+@article{Lucarini2017,
+  title={Edge states in the climate system: exploring global instabilities and critical transitions},
+  author={Lucarini, Valerio and B{\'o}dai, Tam{\'a}s},
+  journal={Nonlinearity},
+  volume={30},
+  number={7},
+  pages={R32},
+  year={2017},
+  publisher={IOP Publishing}
+}
+
+@article{Mehling2023,
+  title={Limits to predictability of the asymptotic state of the Atlantic Meridional Overturning Circulation in a conceptual climate model},
+  author={Mehling, Oliver and B{\"o}rner, Reyk and Lucarini, Valerio},
+  journal={arXiv preprint arXiv:2308.16251},
+  year={2023}
 }

docs/src/basins.md

@@ -29,6 +29,15 @@ basins_fractal_test
 uncertainty_exponent
 ```
 
+## Edge tracking and edge states
+The edge tracking algorithm allows to locate and construct so-called edge states (also referred to as *Melancholia states*) embedded in the basin boundary separating different basins of attraction. These could be saddle points, unstable periodic orbits or chaotic saddles. The general idea is that these sets can be found because they act as attractors when restricting to the basin boundary.
+
+```@docs
+edgetracking
+EdgeTrackingResults
+bisect_to_edge
+```
+
 ## Tipping points
 This page discusses functionality related with tipping points in dynamical systems with known rule. If instead you are interested in identifying tipping points in measured timeseries, have a look at [TransitionIndicators.jl](https://github.com/JuliaDynamics/TransitionIndicators.jl).
 

docs/src/continuation.md

@@ -24,6 +24,9 @@ RecurrencesFindAndMatch
 
 ```@docs
 match_statespacesets!
+Centroid
+Hausdorff
+StrictlyMinimumDistance
 replacement_map
 set_distance
 setsofsets_distances

docs/src/examples.md

@@ -321,7 +321,7 @@ fig
 
 To achieve even better results for this kind of problematic systems than with previuosly introduced `Irregular Grids`  we provide a functionality to construct `Subdivision Based Grids` in which
 one can obtain more coarse or dense structure not only along some axis but for a specific regions where the state space flow has
-significantly different speed. [`subdivided_based_grid`](@ref) enables automatic evaluation of velocity vectors for regions of originally user specified
+significantly different speed. [`subdivision_based_grid`](@ref) enables automatic evaluation of velocity vectors for regions of originally user specified
 grid to further treat those areas as having more dense or coarse structure than others.
 
 ```@example MAIN
@@ -536,7 +536,7 @@ plot_basins_curves(aggregated_fractions, prange;
 
 ## Trivial featurizing and grouping for basins fractions
 
-This is a rather trivial example showcasing the usage of [`AttractorsViaFeaturizing`](@ref). Let us use once again the magnetic pendulum example. For it, we have a really good idea of what features will uniquely describe each attractor: the last points of a trajectory (which should be very close to the magnetic the trajectory converged to). To provide this information to the [`AttractorsviaFeaturizing`](@ref) we just create a julia function that returns this last point
+This is a rather trivial example showcasing the usage of [`AttractorsViaFeaturizing`](@ref). Let us use once again the magnetic pendulum example. For it, we have a really good idea of what features will uniquely describe each attractor: the last points of a trajectory (which should be very close to the magnetic the trajectory converged to). To provide this information to the [`AttractorsViaFeaturizing`](@ref) we just create a julia function that returns this last point
 
 ```@example MAIN
 using Attractors
@@ -871,3 +871,92 @@ animate_attractors_via_recurrences(mapper, u0s)
 <source src="https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/videos/attractors/recurrence_algorithm.mp4?raw=true" type="video/mp4">
 </video>
 ```
+
+## Edge tracking
+
+To showcase how to run the [`edgetracking`](@ref) algorithm, let us use it to find the
+saddle point of the bistable FitzHugh-Nagumo (FHN) model, a two-dimensional ODE system
+originally conceived to represent a spiking neuron.
+We define the system in the following form:
+
+```@example MAIN
+function fitzhugh_nagumo(u,p,t)
+    x, y = u
+    eps, beta = p
+    dx = (x - x^3 - y)/eps
+    dy = -beta*y + x
+    return SVector{2}([dx, dy])
+end
+
+params = [0.1, 3.0]
+ds = CoupledODEs(fitzhugh_nagumo, ones(2), params, diffeq=(;alg = Vern9(), reltol=1e-11))
+```
+
+Now, we can use Attractors.jl to compute the fixed points and basins of attraction of the
+FHN model.
+
+```@example MAIN
+xg = yg = range(-1.5, 1.5; length = 201)
+grid = (xg, yg)
+mapper = AttractorsViaRecurrences(ds, grid; sparse=false)
+basins, attractors = basins_of_attraction(mapper)
+
+for i in 1:length(attractors)
+    println(attractors[i][1])
+end
+```
+
+The `basins_of_attraction` function found three fixed points: the two stable nodes of the
+system (labeled A and B) and the saddle point at the origin. The saddle is an unstable
+equilibrium and can therefore not be found by simulation, but we can find it using the
+[`edgetracking`](@ref) algorithm. For illustration, let us initialize the algorithm from
+two initial conditions `init1` and `init2` (which must belong to different basins
+of attraction, see figure below).
+
+```julia
+attractors_AB = Dict(1 => attractors[1], 2 => attractors[2])
+init1, init2 = [-1.0, -1.0], [-1.0, 0.2]
+```
+
+Now, we run the edge tracking algorithm:
+
+```@example MAIN
+bisect_thresh, diverge_thresh, Δt, abstol = 1e-3, 2e-3, 1e-5, 1e-3
+et = edgetracking(ds, attractors_AB; u1=init1, u2=init2,
+    bisect_thresh, diverge_thresh, Δt, abstol)
+
+et.edge[end]
+```
+
+The algorithm has converged to the origin (up to the specified accuracy) where the saddle
+is located. The figure below shows how the algorithm has iteratively tracked along the basin
+boundary from the two initial conditions (red points) to the saddle (green square). Points
+of the edge track (orange) at which a re-bisection occured are marked with a white border.
+The figure also depicts two trajectories (blue) intialized on either side of the basin
+boundary at the first bisection point. We see that these trajectories follow the basin
+boundary for a while but then relax to either attractor before reaching the saddle. By
+counteracting the instability of the saddle, the edge tracking algorithm instead allows to
+track the basin boundary all the way to the saddle, or edge state.
+
+```@example MAIN
+traj1 = trajectory(ds, 2, et.track1[et.bisect_idx[1]], Δt=1e-5)
+traj2 = trajectory(ds, 2, et.track2[et.bisect_idx[1]], Δt=1e-5)
+
+fig = Figure()
+ax = Axis(fig[1,1], xlabel="x", ylabel="y")
+heatmap_basins_attractors!(ax, grid, basins, attractors, add_legend=false, labels=Dict(1=>"Attractor A", 2=>"Attractor B", 3=>"Saddle"))
+lines!(ax, traj1[1][:,1], traj1[1][:,2], color=:dodgerblue, linewidth=2, label="Trajectories")
+lines!(ax, traj2[1][:,1], traj2[1][:,2], color=:dodgerblue, linewidth=2)
+lines!(ax, et.edge[:,1], et.edge[:,2], color=:orange, linestyle=:dash)
+scatter!(ax, et.edge[et.bisect_idx,1], et.edge[et.bisect_idx,2], color=:white, markersize=15, marker=:circle, zorder=10)
+scatter!(ax, et.edge[:,1], et.edge[:,2], color=:orange, markersize=11, marker=:circle, zorder=10, label="Edge track")
+scatter!(ax, [-1.0,-1.0], [-1.0, 0.2], color=:red, markersize=15, label="Initial conditions")
+xlims!(ax, -1.2, 1.1); ylims!(ax, -1.3, 0.8)
+axislegend(ax, position=:rb)
+fig        
+```
+
+In this simple two-dimensional model, we could of course have found the saddle directly by
+computing the zeroes of the ODE system. However, the edge tracking algorithm allows finding
+edge states also in high-dimensional and chaotic systems where a simple computation of
+unstable equilibria becomes infeasible.

docs/src/index.md

@@ -10,9 +10,7 @@ using CairoMakie, Attractors
 
 ## Latest news
 
-- Our paper on the global stability analysis framework offered by Attractors.jl ([`continuation`](@ref)) and the novel continuation offered by [`RecurrencesFindAndMatch`](@ref) is published as a _Featured Article_ in Chaos (https://pubs.aip.org/aip/cha/article/33/7/073151/2904709/Framework-for-global-stability-analysis-of) and has been featured in the AIP publishing showcase (https://www.growkudos.com/publications/10.1063%25252F5.0159675/reader)
-- New function [`minimal_fatal_shock`](@ref)
-- New function [`match_continuation!`](@ref) which improves the matching during a continuation process where attractors disappear and reappear.
+- New functions [`edgetracking`](@ref) and [`bisect_to_edge`](@ref) added that implement an **edge tracking algorithm** to find saddles or *edge states* in dynamical systems, also when they are unstable chaotic sets.
 
 ## Outline of Attractors.jl
 

docs/src/visualization.md

@@ -11,7 +11,7 @@ heatmap_basins_attractors(grid, basins, attractors; kwargs...)
 heatmap_basins_attractors!(ax, grid, basins, attractors; kwargs...)
 ```
 
-## Common plotting keywords
+## [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.

src/Attractors.jl

@@ -18,6 +18,7 @@ include("dict_utils.jl")
 include("mapping/attractor_mapping.jl")
 include("basins/basins.jl")
 include("continuation/basins_fractions_continuation_api.jl")
+include("boundaries/edgetracking.jl")
 include("deprecated.jl")
 include("tipping/tipping.jl")