Stagger and step method for approximating chaotic saddles

docs/refs.bib

@@ -243,6 +243,29 @@ @Article{Ritchie2023
   DOI = {10.5194/esd-14-669-2023}
 }
 
+
+
+@article{sweet2001,
+  title={Stagger-and-step method: Detecting and computing chaotic saddles in higher dimensions},
+  author={Sweet, David and Nusse, Helena E and Yorke, James A},
+  journal={Physical Review Letters},
+  volume={86},
+  number={11},
+  pages={2261},
+  year={2001},
+  publisher={APS}
+}
+
+@article{sala2016,
+  title={Searching chaotic saddles in high dimensions},
+  author={Sala, Matteo and Leit{\~a}o, Jorge C and Altmann, Eduardo G},
+  journal={Chaos: An Interdisciplinary Journal of Nonlinear Science},
+  volume={26},
+  number={12},
+  year={2016},
+  publisher={AIP Publishing}
+}
+
 @ARTICLE{Klinshov2015,
   title     = "Stability threshold approach for complex dynamical systems",
   author    = "Klinshov, Vladimir V and Nekorkin, Vladimir I and Kurths,

docs/src/api.md

@@ -98,13 +98,14 @@ uncertainty_exponent
 test_wada_merge
 ```
 
-## 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.
+## Edge tracking, edge states and chaotic saddles
+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. Another technique to get a pseudo trajectory close to a saddle is the stagger-and-step method that requires little information on the dynamical system. 
 
 ```@docs
 edgetracking
 EdgeTrackingResults
 bisect_to_edge
+stagger_and_step
 ```
 
 ## Tipping points
@@ -234,4 +235,4 @@ plot_basins_attractors_curves
 
 ```@docs
 animate_attractors_continuation
-```
+```

docs/src/examples.md

@@ -730,7 +730,6 @@ 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"))
@@ -749,3 +748,64 @@ In this simple two-dimensional model, we could of course have found the saddle d
 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.
+
+
+## Invariant saddle of a dynamical system 
+
+The stagger-and-step method approximates the invariant 
+non-attracting set governing the chaotic transient dynamics 
+of a system, namely the stable manifold of a chaotic saddle. 
+
+Given the dynamical system `ds` and a initial guess `x0` in a 
+region *with no attractors*, the algorithm provides `N` points 
+close to the  stable manifold that escape from the region 
+after at least `Tm` steps of `ds`. 
+
+We first set the dynamical system, in our case we set up two coupled Hénon map
+that are known to have a chaotic saddle that generates chaotic transients
+before the trajectories escape: 
+
+```@example MAIN
+function F!(du, u ,p, n)
+    x,y,u,v = u
+    A = 3; B = 0.3; C = 5.; D = 0.3; k = 0.4; 
+    du[1] = A - x^2 + B*y + k*(x-u)
+    du[2] = x
+    du[3] = C - u^2 + D*v + k*(u-x)
+    du[4] = u
+    return 
+end
+ds = DeterministicIteratedMap(F!, zeros(4)) 
+```
+
+Next we define a region in the phase space that should not contain attractors. 
+Using this region we also define a `sampler` and a membership function `isinside`: 
+
+```@example MAIN
+R_min = [-4; -4.; -4.; -4.] 
+R_max = [4.; 4.; 4.; 4.]
+sampler, isinside = statespace_sampler(HRectangle(R_min,R_max))
+```
+
+And we are ready! We can now call the function `stagger_and_step` with an initial 
+condition `x0`: 
+
+```@example MAIN
+x0 = sampler()
+v = stagger_and_step!(ds, x0, 10000, isinside; stagger_mode = :adaptive, δ = 1e-4, Tm = 10, max_steps = Int(1e5), δ₀ = 2.) 
+```
+The `stagger_mode` keyword select the type of search in the 
+phase space to stick close to the saddle at each step. 
+The mode `:adaptive` adapts the radius of the stochastic
+search as a function of the success of the search process. 
+
+Finally we can represent a projection of the chaotic saddle found in this 
+example: 
+
+```@example MAIN
+fig = Figure()
+ax = Axis(fig[1,1], xlabel="x", ylabel="y")
+v = StateSpaceSet(v)
+scatter!(ax, v[:,1], v[:,3]; markersize = 3)
+fig
+```

src/Attractors.jl

@@ -19,6 +19,7 @@ include("basins/basins.jl")
 include("continuation/basins_fractions_continuation_api.jl")
 include("matching/matching_interface.jl")
 include("boundaries/edgetracking.jl")
+include("boundaries/stagger_and_step.jl")
 include("deprecated.jl")
 include("tipping/tipping.jl")
 

src/boundaries/stagger_and_step.jl

@@ -0,0 +1,241 @@
+export stagger_and_step!, stagger_trajectory!
+using LinearAlgebra: norm
+using Random
+using ProgressMeter
+
+
+"""
+    stagger_trajectory!(ds, x0, Tm, isinside; kwargs...) -> xi
+
+Return a point `xi` which _guarantees_ `T(xi) > 
+Tm` with a random walk search around the initial coordinates 
+`x0`. `T(xi)` is the escape time of the initial condition `xi`.  
+In the case where the algorithm cannot find this point the algorithm 
+returns nothing and the parameters must be adjusted to find the suitable point.
+This is an auxiliary function for [`stagger_and_step`](@ref).
+Keyword arguments and definitions are identical for both functions. 
+
+The initial search radius `δ₀` is big, `δ₀ = 1.0` by default.
+"""
+function stagger_trajectory!(ds, x0, Tm, isinside; δ₀ = 1., stagger_mode = :exp,
+        max_steps = Int(1e5), γ = 1.1, max_escape_time = 10000)
+    T = escape_time!(ds, x0, isinside; max_escape_time)
+    xi = copy(x0) 
+    while !(T > Tm)  # we must have T > Tm at each step 
+        xi, T = stagger!(ds, xi, δ₀, T, isinside; γ, stagger_mode, max_steps)
+        if T < 0
+            return nothing
+        end 
+
+    end
+    return xi
+end
+
+
+"""
+    stagger_and_step!(ds::DynamicalSystem, x0, N::Int, 
+        isinside::Function; kwargs...) -> trajectory
+
+Implement the stagger-and-step method 
+[Sweet2001](@cite) and [Sala2016](@cite) to approximate the invariant 
+non-attracting set governing the chaotic transient dynamics 
+of a system, namely the stable manifold of a chaotic saddle. 
+
+Given the dynamical system `ds` and a initial guess `x0` in a 
+region *with no attractors* defined by the membership function 
+`isinside`, the algorithm provides `N` points close to the 
+stable manifold that escape from the region after at least `Tm` 
+steps of `ds`. The search is stochastic and depends on the 
+parameter `δ` defining a (small) neighborhood of search. 
+
+The function `isinside(x)` returns `true` if the point `x` is 
+inside the chosen bounded region and `false` otherwise. See 
+[`statespace_sampler`](@ref) to construct this function.
+
+## Keyword arguments
+* `δ = 1e-10`: A small number constraining the random
+  search around a particular point. The interpretation of this 
+  number will depend on the distribution chosen for the 
+  generation (see `stagger_mode`).
+
+* `Tm = 30`: The minimum number of iterations of `ds` before 
+  the trajectory escapes from the bounding box defined by
+  `isinside`. 
+
+* `max_steps = 10^5`: The search for a new candidate point may 
+  fail at some point. If the search fails after `max_steps`, 
+  a new initial point is set and the method starts from a new
+  point.
+
+* `max_escape_time = 10000`: If the trajectory stays in the 
+  defined region after `max_escape_time` steps, there is probably
+  an attractor in the region and the algorithm will fail. 
+
+* `stagger_mode = :exp`: There are several ways to produce 
+  candidate points `x` that have to fulfill the condition 
+  `T(x) > Tm`. The available methods are: 
+
+    * `:exp`: An candidate following an truncated exponential 
+      distribution in a random direction `u` around the current
+      `x` such that `x_c = x + u*r`. `r = 10^-s` with `s` taken
+      from a uniform distribution in [-15, δ]. This mode fails
+      often but still manage to provide long enough stretch of
+      trajectories. 
+
+    * `:unif`: The next candidate is `x_c = x + u*r` with `r` 
+      taken from a uniform distribution [0,δ]. 
+    * `:adaptive`: The next candidate is `x_c = x + u*r` with 
+      `r` drawn from a gaussian distribution with variance  δ.
+      The variance changes according to a free parameter `γ`  
+      such that `δ = δ/γ` if no candidate is found and `δ = δ*γ`
+      when it succeeds. `γ = 1.1`: It is the free parameter for 
+      the adaptive stagger method. 
+* `δ₀ = 1.0`: This is the radius for the first stagger 
+  trajectory search. The algorithm looks for a point 
+  sufficiently close to the saddle before switching to the 
+  stagger-and-step routine. The search radius must be large 
+  enough to find a suitable initial. 
+
+## Description 
+
+The method relies on the stagger-and-step algorithm that 
+computes points close to the saddle that escapes in a time 
+`T(x_n) > Tm`. The function `T` represents the escape time 
+from a region defined by the user (see the argument 
+`isinside`).
+
+Given the dynamical mapping `F`, if the iteration `x_{n+1} = 
+F(x_n)` respects the condition `T(x_{n+1}) > Tm` we accept 
+this next point, this is the _step_ part of the method. If 
+not, the method search randomly the next point in a 
+neighborhood following a given probability distribution, this 
+is the _stagger_ part. This part sometimes fails to find a new 
+candidate and a new starting point of the trajectory is chosen 
+within the defined region. See the keyword argument 
+`stagger_mode` for the different available methods.   
+
+The method produces a pseudo-trajectory of `N` points δ-close
+to the stable manifold of the chaotic saddle. 
+
+"""
+function stagger_and_step!(ds::DynamicalSystem, x0, N::Int, isinside::Function; δ = 1e-10, Tm  = 30, 
+    γ = 1.1, max_steps = Int(1e5), max_escape_time = 10000, stagger_mode = :exp, δ₀ = 1., 
+    show_progress = true)
+
+    progress = ProgressMeter.Progress(
+        N; desc = "Saddle estimation: ", dt = 1.0
+    )
+    xi = stagger_trajectory!(ds, x0, Tm, isinside; δ₀, stagger_mode = :unif, max_steps) 
+    if isnothing(xi)
+        error("Cannot find a stagger trajectory. Choose a different starting point or 
+                search radius δ₀.")
+    end
+
+    v = Vector{typeof(current_state(ds))}(undef, N)
+    v[1] = xi
+    for n in 1:N
+        show_progress && ProgressMeter.update!(progress, n)
+        if escape_time!(ds, xi, isinside; max_escape_time) > Tm
+            reinit!(ds, xi; t0 = 0)
+        else
+            xp, Tp = stagger!(ds, xi, δ, Tm, isinside; stagger_mode, max_steps, γ)
+            # The stagger step may fail. We reinitiate the algorithm from a new initial condition.
+            if Tp < 0
+                xp = stagger_trajectory!(ds, x0, Tm, isinside; δ₀, stagger_mode = :exp, max_steps, γ) 
+                if isnothing(xp)
+                    error("Cannot find a stagger trajectory. Choose a different starting 
+                          point or search radius δ₀.")
+                end
+                δ = 0.1
+            end
+            reinit!(ds, xp; t0 = 0)
+        end 
+        step!(ds)
+        xi = copy(current_state(ds))
+        v[n] = xi
+    end
+    return v
+end
+
+
+
+function escape_time!(ds, x0, isinside; max_escape_time = 10000) 
+    x = copy(x0) 
+    set_state!(ds,x)
+    reinit!(ds, x; t0 = 0)
+    k = 1; 
+    while isinside(x) 
+        if k > max_escape_time 
+            error("The trajectory did not escape for ", k, " steps, you probably have \
+an attractor in the defined region.")
+        end
+        step!(ds)
+        x = current_state(ds)
+        k += 1
+    end
+    return current_time(ds)
+end
+
+function rand_u(δ, n; stagger_mode = :exp)
+    if stagger_mode == :exp 
+        a = -log10(δ)
+        s = (15-a)*rand() + a
+        u = (rand(n).- 0.5)
+        u = u/norm(u)
+        return u*10^-s
+    elseif stagger_mode == :unif
+        s = δ*rand()
+        u = (rand(n).- 0.5)
+        u = u/norm(u)
+        return u*s
+    elseif stagger_mode == :adaptive
+        s = δ*randn()
+        u = (rand(n).- 0.5)
+        u = u/norm(u)
+        return u*s
+    end
+end
+
+"""
+    stagger!(ds, x0, δ, Tm, isinside; kwargs...) -> stagger_point
+
+This function searches a new candidate in a neighborhood of x0 with a random search 
+depending on some distribution. If the search fails it returns a negative time.
+"""
+function stagger!(ds, x0, δ, Tm, isinside; max_steps = Int(1e6), γ = 1.1, stagger_mode = :exp, 
+        verbose = false, max_escape_time = 10000)
+    Tp = 0; xp = zeros(length(x0)); k = 1; 
+    T0 = escape_time!(ds, x0, isinside; max_escape_time)
+    if !isinside(x0)
+        error("x0 must be in grid")
+    end
+    while Tp ≤ Tm 
+        xp = x0 .+ rand_u(δ,length(x0); stagger_mode)
+
+        if k > max_steps 
+           if verbose 
+               @warn "Stagger search fails, δ is too small or T is too large. 
+                We reinitiate the algorithm
+               "
+           end
+           return 0,-1
+        end
+        Tp = escape_time!(ds, xp, isinside; max_escape_time)
+        if stagger_mode == :adaptive
+            # We adapt the variance of the search
+            # if the alg. can't find a candidate
+            if Tp < T0
+                δ = δ/γ
+            elseif Tp == T0
+                δ = δ*γ
+            end
+            if Tp == Tm
+                # The adaptive alg. accepts T == Tp
+                return xp, Tp
+            end
+        end
+        k = k + 1
+    end
+    return xp, Tp
+end
+