Extension of basin general to n-dim

[awage]

I found a way to extend to n-dimension. It seems to work with the magnetic pendulum quite fine. I am amazed how good Julia handle n-dim arrays and CartesianIndex is a smashing. It was not much to change after all.

I tried in 4D and it worked but computational load explodes quickly.

ds = Systems.magnetic_pendulum(γ=1, d=0.3, α=0.2, ω=0.5, N=3)
integ = integrator(ds, u0=[0,0,0,0], reltol=1e-9)
xg=range(-2,2,length=100)
yg=range(-2,2,length=100)
zg=range(-1,1,length=10)
bsn , att = basins_general([xg, yg, zg], ds; dt=1., idxs=1:3, complete_state=[0])

I have still some tests to do and also I have to check the complete_state functionality. It may be broken.

[Datseris]

We should test this PR with some system that has attractors spanning all three dimensions. For example: https://www.sciencedirect.com/science/article/pii/S0375960114002850?via%3Dihub

Furthermore we should try to address #183 as much as possible in this PR. I will be doing some performance improvements over the weekend.

[Datseris]

@awage keep those performance annotations coming. Tomorrow I will work on the performance of this PR and will incorporate fixes for all comments.

[Datseris]

@awage I'm working on this right now. One thing that we should definitely improve before merging is the code readability. One thing we should definitely improve before merging is the procedure! function:

1. It is very large, 100 lines of code. It should be split up into more, smaller, atomic functions, each being up to 20 lines of code.
2. It's name does not convey it's purpose, so I will rename it to something that includes the "basins"

[Datseris]

While working on the source code I have found extremely confusing the following: what the states u represent. Typically u as a symbol is used to represent the full system state. In the source, it is not at all for me clear when u means the full state, or the projection of the state on the sub-space denoted by grid. @awage do you think you can go through and clarify this? I don't think I can do it.

Let's do the following:

• u is used exclusively to denote the full system state, as is done in the rest of the library.
• y can be used to denote the projected state, which is of smaller dimensionality, and same dimensionality as the grid.

[awage]

Ok, I will work on these issues next week. Since the algorithm is a finite state machine maybe there is way to define states like :search_attractor :atractor_found and so on. I will try to find a way.

[Datseris]

I've addressed most of the concerns, did a lot of performance improvements, and increased code clarity a lot. What remains is:

• Use symbol y for projected state and u for full state
• Replace all literal constants like 60 or 20*60 with named keyword arguments.
• Address the TODO markings I've left in the source code.

Have a look and let me know what can be improved. Even though I've tried to make the source be as identical as originally regarding output, this was really hard to control, as I've done some really large changes. As a result, now the e.g. Henon map test fails. It says that count(basin .== 1) is 4261 instead of 4700 and that count(basin .== -1) is 5739 instead of 5730. I couldn't pinpoint exactly which change made this disparity.

[awage]

We should test this PR with some system that has attractors spanning all three dimensions. For example: https://www.sciencedirect.com/science/article/pii/S0375960114002850?via%3Dihub

This system is a pain in the neck with its non-attracting torus, the algorithm only detect pieces of the torus and the attractor. There is more gentle system with attractors in 3D:

http://inaesp.org/PublicJG/freire_et_al_CHAOS2008.pdf (https://doi.org/10.1063/1.2953589)

basins_general works just fine:

[Datseris]

Awesome, you can put this system in DynamicalSystemsBase.jl if you'd like!

[awage]

I have complety changed the algorithm. I think it is much less cryptic right now. It will be easier to maintain and to extend. In the meantime I solved some problems in the original algorithm since they were cutting corners.

The idea now is that you are in particular state and depending on the cell you hit you have a transition. Then depending on an internal counter of the state you may have actions or transition (I left a schematic drawing). There are still some minor issues to solve.

It took some serious brain time to put it that way!

EDIT: there are some transitions missing on the picture.

[Datseris]

Beautiful. Just beautiful. I'll do another round of review when I have the time. What I'm thinking is that now, with the approach we have, how much of it is really based on the work of Nussa/Yorke? Probably better to make a nice figure of the state machine and use this as docs instead of saying "we implement the method of Nussa/Yorke"?

[awage]

The basic idea is the same:

• you have a grid and you keep track of the events in an array
• you make sure to keep the transient times to the attractors as short as possible.

The internals of the algorithm are now very different. I will document the algorithm as clearly as possible later this week.

There is something I need to solve. The :bas_hit state is only valid for 2D maps. It assumes that when you hit the same basin several times in a row you will end to the attractor of that basin. This is not always reliable in dimension above 2. But for 2D it is really quick so it is worth keeping it.

[awage]

@Datseris : I have a question about the functionality of basins_general. You left some comments that suggest that the grid should have the dimension of the phase space of the dynamical system. It is not necessary, you can always use a grid in a lower dimension projected space. For example:

ds = Systems.magnetic_pendulum(γ=1, d=0.3, α=0.2, ω=0.5, N=3)
integ = integrator(ds, u0=[0,0,0,0], reltol=1e-9)
xg=range(-2,2,length=100)
yg=range(-2,2,length=100)
zg=range(-1,1,length=10)
bsn,att = basins_general((xg, yg, zg), ds; dt=1., idxs=1:3, complete_state=[0])

This should work. Maybe you want to cleave the functionalities: on the one hand use basins_2D for projection and basins_general for n-dim cases. This is a design choice, since we are not limited by the algorithm.

[Datseris]

You left some comments that suggest that the grid should have the dimension of the phase space of the dynamical system

Oh sorry, this was never my intention. Maybe I misspoke in the comment. Sure, the grid can have smaller dimensionality.

My idea is that basins_general becomes the basins function that has all important documentation, while basins_2D will become just a convenience that also allows for poincare maps and stroboscopic maps.

I also suggest that we rename basins_general to basins_of_attraction or attractor_basins. Just basins would be nice, but I'd prefer to leave this name unbinded so that users can use it for the actual basins matrix.

[awage]

I think we can get everything working with basins_general. It seems that maps are working and Poincaré maps also with dispatch on basins_general. The interface would be a unique function that in most case doesn't need much keywords arguments.

The only systems left are stroboscopic maps. But it can be solved if we allow to step exactly at T.

What is your catch?

[Datseris]

You are right. Everything can be done with a single function. For stroboscopic maps, a user would provide diffeq = (alg = Tsit5(), dt = T, adaptive = false) to ensure the stepping is one period precicely.

[Datseris]

todo for me: we can make complete_state a functor to avoid re-creating states all the time.

[awage]

From my end it is OK and closes #183 and #182 . But I detected a problem with the tag propagation of diffeq. I tried to add a callback on an equation but the program doesn't respond. I don't know if it should be addressed in this PR:

using DifferentialEquations
 using DynamicalSystems
function forced_pendulum(u, p, t)
    @inbounds begin
    d = p[1]; F = p[2]; omega = p[3]
    du1 = u[2]
    du2 = -d*u[2] - sin(u[1])+ F*cos(omega*t)
    return SVector{2}(du1, du2)
    end
end

# We have to define a callback to wrap the phase in [-π,π]
function affect!(integrator)
    if integrator.u[1] < 0
        integrator.u[1] += 2*π
    else
        integrator.u[1] -= 2*π
    end
end
condition(u,t,integrator) = (integrator.u[1] < -π  || integrator.u[1] > π)
cb = DiscreteCallback(condition,affect!)

F = 1.66 ; ω = 1.; d=0.2;
ds=ContinuousDynamicalSystem(forced_pendulum, rand(2),[d, F, ω])
xg = range(-pi,pi,length=100); yg = range(-2.,4.,length=100)
bsn,att = basins_of_attraction((xg, yg), ds; T=2*pi/ω, diffeq = ( save_everystep=false, reltol=1e-9, callback =cb))

[Datseris]

Thanks a lot for this great work. I'll have a look over the weekend, fix as much as I can, and merge and tag. Remaining issues (e.g. if I can't fix the callback) will be opened as new issues.

[Datseris]

@awage I am doing a performance improvement now. One thing I've noticed is that during the re-initialization process we were constantly creating new vectors. So I created a helper type that stores two dummy vectors and does in-place creation of a new state. What I'm having difficulty is extending this to the Poincare Map. Why is it that we're doing completely different handling for it?

[Datseris]

ok I got what's going on, I'll push a fix soon

[awage]

@awage I am doing a performance improvement now. One thing I've noticed is that during the re-initialization process we were constantly creating new vectors. So I created a helper type that stores two dummy vectors and does in-place creation of a new state. What I'm having difficulty is extending this to the Poincare Map. Why is it that we're doing completely different handling for it?

Pmap has its own integrator, it is a bit special.

[Datseris]

Alright, I'm finally done. I was a bit discontent with how much more complexity the Poincare map was adding to the code. I decided to eliminate the complexity, in favor of treating the "state" of the Poincare map as the full system state. I will add a bit more info about this at the docs. I have to open a PR at DynamicalSystems.jl for the new documentation.

[Datseris]

Hm, something may be off. I'm finding 4 attractors for the duffing oscillator:

four colors from

using DynamicalSystems, PyPlot
ω=1.0; f = 0.2
ds = Systems.duffing([0.1, 0.25]; ω, f, d = 0.15, β = -1)
basins, attractors = basins_of_attraction((xg, yg), ds; T=2π/ω)
pcolormesh(xg, yg, basins')

[Datseris]

I've also noticed that the uncertainty_exponent is still limited to a 2D grid. Is there any fundamental limitation that keeps it there, or we can also improve this algorithm to be arbitrary dimensional?

[Datseris]

I have the documentation PR ready. https://github.com/JuliaDynamics/DynamicalSystems.jl/pull/143

So after we resolve the the multiple basins for duffing I merge and tag!

[awage]

Hm, something may be off. I'm finding 4 attractors for the duffing oscillator:

Something awkward happens with SimpleATsit5() again. It does not seems related to DynamicalSystems, MWE:

using  DynamicalSystems, Plots, SimpleDiffEq, OrdinaryDiffEq
ds = Systems.duffing([0.1, 0.25]; ω = 1., f = 0.2, d = 0.15, β = -1)
T = 2π/1.
integ = integrator(ds, alg=SimpleATsit5())
integ2 = integrator(ds, alg=Tsit5())
integ.u = integ2.u = [-0.2, 0.1]
v=Vector{SVector{2,Float64}}()
v2=Vector{SVector{2,Float64}}()
for k = 1:700
    step!(integ, T, true)
    step!(integ2, T, true)
    push!(v,integ.u)
    push!(v2,integ2.u)
end
v = Dataset(v); v2 = Dataset(v2)
plot(v[:,2],seriestype=:scatter)
plot!(v2[:,2],seriestype=:scatter)

If you look at the picture the problem is not that they go to different attractor (it could happen since it is chaotic and we have different solvers). The problem is the very long transient before getting to the steady state. It is maybe the stepping problem that pops up again.

[Datseris]

I'll simply modify the documentation PR to instruct to users to use a higher accuracy solver than the default one.

[Datseris]

This is good to go and merge from my end, if you don't have any final remarks.

[awage]

I've also noticed that the uncertainty_exponent is still limited to a 2D grid. Is there any fundamental limitation that keeps it there, or we can also improve this algorithm to be arbitrary dimensional?

Well, this is something to dig. In the original algorithm they take "balls" in one dimension. The reason is that you need "balls" of radius ε to measure correctly the uncertainty exponent α. When the resolution is different along axes you may have arrays of different size along the dimensions. It causes trouble.

On the other hand, I did not find this solution satisfactory since you may have a bias due to the direction of the axis relative to the boundary. I will try to find where they say that it is independent of the chosen ball.

[awage]

Thank you very much @Datseris! for me everything fine.