Automatic Δt estimation for `basins_of_attraction`

[Datseris]

WIP.

[Datseris]

@awage here are my measurements for Δt with the new function I wrote in this PR:

• Lorenz84 with 60x60x60 = 0.01
• Magnetic pendulum with 200x200 = 0.01
• Thomas cyclical (3D, not poincare) with 250x250x250 = 0.05
• Lorenz96 with 50x50x50 = 0.01

So perhaps we need to re-consider what we think and what we say about Δt needed to be "time to cross 1 cell"? At least for Thomas and Magnetic Δt = 1.0 yields really good results even though it is so much larger than ...

[Datseris]

I also did a comparison regarding performance. What's the time when using Δt = 1, and when using the "recommended" from the algorithm I wrote:

• Thomas cyclical 3D: 0:02:33 versus 0:00:58

Problem: using the automated Δt of 0.13 yielded wrong results!!! It found more basins that normal:

Δt = 1:

Δt = 0.13:

I think the story isn't really totally clear to me at the moment. How can the Δt have such a massive impact on the output...?

[awage]

There is something to dig into here. What is the grid size of the latest simulations? 250x250x250?

[awage]

I can't reproduce these results. I have tried with different grids and time step and I always get 5 attractors. For example:

ds = Systems.thomas_cyclical(b = 0.1665)
xg = yg = zg = range(-6.0, 6.0; length = 80)
ChaosTools.automatic_Δt_basins(ds,(xg,yg,zg))
-> 0.25882088282464427
basins, attractors = basins_of_attraction((xg, yg, zg), ds; Δt = 0.2588)
Dict{Int16, Dataset{3, Float64}} with 5 entries:
  5 => 3-dimensional Dataset{Float64} with 1 points
  4 => 3-dimensional Dataset{Float64} with 129 points
  2 => 3-dimensional Dataset{Float64} with 344 points
  3 => 3-dimensional Dataset{Float64} with 343 points
  1 => 3-dimensional Dataset{Float64} with 1 points

EDIT: I got it with resolution 150 it stops working. I will start from there.

[awage]

Problem: using the automated Δt of 0.13 yielded wrong results!!! It found more basins that normal

@Datseris can you give me the parameters? I can't reproduce this bug. I always get 5 attractors.

[Datseris]

Sure, sorry for the late reply. I'm using:

b = 0.1665
system = :thomas_cyclical
xg = yg = zg = range(-6.0, 6.0; length = 100)
basin_kwargs = (mx_chk_hit_bas = 40, mx_chk_att = 5)

(notice its 3D basins)

the automated Δt finds Δt = 0.14032049283451245.

[awage]

There is something odd. I have tried with the parameters you gave me and I obtained a correct result for different time steps (0.14, 0.2 and 1.0). I suspect we are not using the same solver since I obtain Δt = 0.206 from automatic_Δt_basins.

Here is my code with the branch automatic_Δt_basins checked out:

using Revise
using DynamicalSystems
using Plots

ds = Systems.thomas_cyclical(b = 0.1665)
xg = yg = zg = range(-6.0, 6.0; length = 100)
dt = ChaosTools.automatic_Δt_basins(ds, (xg,yg,zg))
basins, attractors = basins_of_attraction((xg, yg, zg), ds; Δt = dt, show_progress = true, mx_chk_hit_bas = 40, mx_chk_att = 5, mx_chk_fnd_att = 100)
heatmap(basins[:,:,50]')

[Datseris]

Oh yeah, I am using the solver we used in the paper, i.e. Vern9!

The Δt you obtain from the ChaosTools function does not depend on the solver. Have a look at the source code. It actually evaluates the vector field at random points and uses the derivative to derive the timescale. It picks points at random hence the result fluctuates. We also need to keep in mind that the step! we use inside the basins_of_attraction function is approximate, not exact. It does not step for exactly Δt. And it also cannot step for less that the actual internal integrator step, which is decided adaptively.

[awage]

I still cannot reproduce the figure. I get the right result for different time steps (which is kind of comforting). There is something I am not doing right but I don't see what. Here is my code:

using Revise
using DynamicalSystems
using Plots
using OrdinaryDiffEq
ds = Systems.thomas_cyclical(b = 0.1665)
xg = yg = zg = range(-6.0, 6.0; length = 100)
dt = ChaosTools.automatic_Δt_basins(ds, (xg,yg,zg))
default_diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9, maxiters = Int(1e12))
dt = 0.14032049283451245
basins, attractors = basins_of_attraction((xg, yg, zg), ds; Δt = dt, show_progress = true, diffeq = default_diffeq, mx_chk_hit_bas = 40, mx_chk_att = 5, mx_chk_fnd_att = 100)
heatmap(basins[:,:,50]')

[Datseris]

Alright, sorry for taking so long. Here is the code that gives 7 attractors for the system:

using DynamicalSystems

b = 0.1665
ds = Systems.thomas_cyclical(; b)
basin_kwargs = (mx_chk_hit_bas = 40, mx_chk_att = 5)
using OrdinaryDiffEq: Vern9
default_diffeq = (alg = Vern9(), reltol = 1e-9, abstol = 1e-9, maxiters = Int(1e12))

xg = yg = zg = range(-6.0, 6.0; length = 100)
grid = (xg, yg, zg)

basins, attractors = basins_of_attraction(
    grid, ds; diffeq = default_diffeq, Δt = 0.14, basin_kwargs..., 
)

length(attractors)

Notice that Δt = 0.14 is the automatically found one.

[awage]

I have pasted and ran the code and I have the good results:

Can you check what is the default value of mx_chk_fnd_att in your implementation?

[Datseris]

My implementation is the one that is now (correctly pushed) here on this branch. I see:

as the default parameters.

Oh man, this is very weird that you're not getting what I am getting! I am not lying I swear I get 7 attractors :D :D :D

[awage]

Wait I have just pulled the latest code. I am checking now.

[Datseris]

Well pulling the latest code shouldn't matter, because I didn't change any of the internals. I only changed the default Δt value in this PR. But in the code we both run, we anyway prescribe a Δt... I'm also looking at this in parallel...

[awage]

OK, I got the seven attractors with the latest code. But on an earlier version I got 5. I will track the bug tomorrow.

[awage]

Got it!!

The parameter mx_chk_loc_att should be higher. I bumped it to 100 and now there are 5 attractors. At any rate these parameters (mx_chk_fnd_att and mx_chk_loc_att ) should be quite high, between 100 and 200.

[Datseris]

Awesome, thanks. I'll bump both of them to 100 in this PR then.

[Datseris]

This finally resolved the difficulties with the Lorenz84 system. Also, I noticed that mx_chk_fnd_att=200, mx_chk_loc_att=50 is good enough. In fact, it matters not how large mx_chk_loc_att is, correct? It has no bearing on the algorithm because it does not color any new cells if mx_chk_fnd_att is already sufficiently large. So mx_chk_loc_att only decides how many attractor points the algorithm finds and puts into the output datasets attractors.

@awage is this correct? Anyways, we now have the beautiful:

scatter points: points found by our algorithm. lines: integrated trajectories.

[awage]

Maybe I should clarify the docstring and/or change the names of the constants:

• mx_chk_fnd_att : this is the number of recurrent visits that must be achieved in order to claim that we have an attractor. However there is no coloring at this stage.
• mx_chk_loc_att : the coloring starts here. Each time that the program finds a unvisited cell or a visited cell it turns it into an attractor cell. This is repeated until we visit mx_chk_loc_att consecutive attractor cells.

It doesn't hurt to increase mx_chk_loc_att. It does color cells as attractors each time it visits a new one.

EDIT: Also it has a marginal impact on the performances.

[Datseris]

Okay, well this was NOT clear to me, because I thought that all cells found during the mx_chk_fnd_att phase are also counted as attractor cells. Can you please make this clear in the paper as well?

[awage]

The reason why you cannot do this is that some cells marked as visited are part of the transient. So you must discard all the visited cells and start over coloring only the attractor cells. I will rephrase this part.

[Datseris]

Yeah, once again, there is another default value to consider... Here is the code:

    ds = Systems.magnetic_pendulum(γ=1, d=0.2, α=0.2, ω=0.8, N=3)
    xg = range(-2,2,length=100)
    yg = range(-2,2,length=100)
    complete_state(y) = SVector(0.0, 0.0)
    basin, attractors = basins_of_attraction((xg,yg), ds;
    idxs=1:2, complete_state, show_progress = false, Δt = 0.02)

For Δt = 1 it is the version before this PR. Putting Δt = nothing is the version of this PR, which does the automated estimation, which yields internally Δt = 0.02. Unfortunately, the results of basins is wrong with this Δt:

Now, if I use the keyword mx_chk_hit_bas = 100, I get the correct result:

So it seems (and this makes perfect logical sense) that most of our algorithm keywords depend on both the grid spacing as well as the choice of Δt... Should we also increase the default mx_chk_hit_bas ?

[Datseris]

actually the above image is still incorrect. Do you see that inner yellow point on the top-right corner of the Teal basin of attraction on the left side? It should be completely yellow but it is contaminated with Teal. I've increased the Δt automatic algorithm to yield 4 times what it used to yield and the problem goes away. Perhaps there are problems created if the trajectory does not leave the cell per step (which sometimes doesn't happen with current Δt)

anyway here is the finally correct version:

[Datseris]

aaaand there are more problems... The "zoom in" functionality no longer works. Continuing from the above code and doing

    xg = yg = range(-2, -1.9, length = 50)
    bas, att = basins_of_attraction((xg,yg), ds;
    idxs = 1:2, complete_state, show_progress = false, mx_chk_hit_bas = 20,
    attractors = attractors, mx_chk_lost = 1000, ε = 1e-3, Δt = 1.0)

(i.e., just running the tests in the basins_tests.jl file)

It doesn't work. All basins get labelled -1. What the hell, how is this even possible! I am using the same Δt as "before this PR", yet the process fails completely.

[awage]

Sorry, I was doing the same. Fixing tests.

[awage]

actually the above image is still incorrect. Do you see that inner yellow point on the top-right corner of the Teal basin of attraction on the left side? It should be completely yellow but it is contaminated with Teal. I've increased the Δt automatic algorithm to yield 4 times what it used to yield and the problem goes away. Perhaps there are problems created if the trajectory does not leave the cell per step (which sometimes doesn't happen with current Δt)

It is better to increase Δt rather than picking a small value. If the trajectory skips a few cells it doesn't matter much. It is sometimes important (Lorenz84 for instance). 4 to 10 times the minimum value sounds reasonable.

[Datseris]

Okay, I make it 10 times the automated and will put an according note to the docs.