Stagger and step method for approximating chaotic saddles

[awage]

This is an implementation the stagger-and-step method [^Sweet2001] and [^Sala2016] to approximate the invariant
non-attracting set governing the chaotic transient dynamics of a system, namely the stable manifold of a chaotic saddle.

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.

The original method is implemented as well as an improvement by Sala et. al.

I added an example file for the docs.

[Datseris]

Sorry Alex, this is definitely on my todos to review. I am just super overwhelemd in the start of the new term and don't have much time. will review at some point!!!

[awage]

No pressure! This is a low priority PR. I did it just to learn a bit more on these techniques.

[Datseris]

Hi @awage thanks a lot for adding this! I don't have very critical comments at the moment, but I will also give a second review. Two main things I'd like is a plottable example (such as magnetic pendulum, newton or henon maps), and also, what happens in the main function if the max steps of the escape time are exceeded? The trajectory doesn't escape but you also don't do anything about it.

[Datseris]

You should put this file as a formal example in the Examples documentation page, and perhaps visualize the output. Is there a 2D or 3D example we can use so that the output can be visualized?

[Datseris]

The (projected) magnetic pendulum could be used for it as well perhaps?

[Datseris]

I don't understand why you do this extra step here.

[awage]

This is step part of the stagger_and_step. If the escape_time! is < Tm then we execute a step. If not, we look for a neighbor that fulfill the condition and after this we do one step.

[Datseris]

This function can be made significantly more performant by creating dedicated structs that store the dummy variables to reduce allocations, and are then used as function-like-objects. I can do this small change once the PR is ready to go.

[awage]

Thanks for the comments. I will start working this week on this PR.

[awage]

@Datseris I have adressed all your comments. I leave the plotable example in the comments below. It is the coupled Hénon map presented in the publications. If I have the time I will try a continuous system.

Note: In the test file I implemented a test to verify that the function escape_time! throws an error as intended. I am not sure that this is the way to test error handling.

using Attractors

# Coupled Hénon maps
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

# The region should not contain any attractors. 
R_min = [-4; -4.; -4.; -4.] 
R_max = [4.; 4.; 4.; 4.]


# Initial box and initial condition.
sampler, isinside = statespace_sampler(HRectangle(R_min,R_max))
x0 = sampler()
df = DeterministicIteratedMap(F!, x0) 
v = stagger_and_step!(df, x0, 10000, isinside; stagger_mode = :adaptive, δ = 1e-4, Tm = 10, max_steps = Int(1e5), δ₀ = 2.) 
v = StateSpaceSet(v)

using CairoMakie
scatter(v[:,1], v[:,3]; markersize = 3)

[Datseris]

Hi @awage I'm sorry I took so long to come back, life has been very intense for me lately. I still see in this PR that there are no changes in the docs folder. At minimum, the documentation string of the main function should be listed somewhere.

Thanks for attaching an example. My concern is that because this is 4D, how much sense does the visualization make? (can you post the output image please?)

[awage]

The projection on 2D of the saddle is like this:

This is the example published in the original paper.

[awage]

I have updated the docs. It seems to compile except for a problem with BifurcationToolkit.

[Datseris]

Seems like there are undefined references in the docs: https://github.com/JuliaDynamics/Attractors.jl/actions/runs/11835746833/job/32979169108?pr=148#step:5:174 (sorry for the late reply! but next time please check the automated integrations to resolve the errors)