pmap get_state returns interpolated state

awage commented 3 years ago

This is an attempt to solve the problem of the Poincaré maps. In this PR the state returned by get_state is the state interpolated on the Poincaré map.

There is a small impact on the performance that I still don't know how to solve in step! (there is a if else). Tests are Ok.

Datseris commented 3 years ago

What's the performance regression when calling poincaresos as in the test suite?

awage commented 3 years ago

No impact on poincaresos but a small dent in performance for step! :

Master: 
julia> @btime psos = poincaresos(ds, (2, 0.0), 1000.0, direction = +1)
1.921 ms (2828 allocations: 253.23 KiB)

julia> @btime  for i in 1:100000; step!(pmap); end
  7.273 ms (500000 allocations: 15.26 MiB)

PR: 
  julia> @btime psos = poincaresos(ds, (2, 0.0), 1000.0, direction = +1)
  1.913 ms (2828 allocations: 253.23 KiB)

  julia> @btime  for i in 1:100000; step!(pmap); end
  7.614 ms (500000 allocations: 15.26 MiB)

Datseris commented 3 years ago

Hm this change in allocations is huge. You don't see its impact on performance because I assume you are using a low dimensional out of place system to do the calculations. For higher dimensional systems this would make a big difference.

In the first place, does this fix the problems we had with the basins of attraction?

awage commented 3 years ago

The change is big from poincaresos to step!. This PR doesn't add any allocation when compared to master version.

This fix the problem with basin_of_attraction. What I though was a fantastic example of basins with the Rikitake ODE is in fact lame and boring:

ds = Systems.rikitake(μ = 0.47, α = 1.0)
plane = (3, 0.0)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6;
    idxs = 1:2, rootkw = (xrtol = 1e-12, atol = 1e-12), reltol=1e-9
)
xg = yg = range(-6.,6.,length=200)
basin, attractors = basins_of_attraction((xg, yg), pmap)

attractors[1]
3-dimensional Dataset{Float64} with 2 points
 -2.20231   1.05136   1.59021e-9
  2.36865  -1.33907  -4.63543e-10

The problem with this parameter is that there is long transient before we get to the attractor. It fools the algorithm that thinks he has found different attractors when you compute basisns without pmap. It can be avoided with mx_chk_bas or mx_chk_att but pmap requires less tweaking. This Rikitake is a tough guy.

imagen

Datseris commented 3 years ago

we should have a simpler example I think for the poincare map. My suggestion would be the Thomas_Cyclical system, for same parameter as the animations we produced and for z = 0 as the hyerplane. The poincare section will have three fixed points (the intersection of the three attracting periodic orbits). And the basis of these three points will be fractal. Can you please apply the poincare-basins to the thomas system? Let's be sure we have something working first, and then look at what should be fixed.

awage commented 3 years ago

I spotted a problem. step! returns nothing when the map is not well defined while get_state always returns a vector. We should be able to detect this in the algorithm. I have to rethink this PR. pmap is a headache but I think it is a valuable tool.

awage commented 3 years ago

While I am not very satisfied with the solution it does not affect the performance and fix the problem of tracking the trajectory on the plane.

Master: 
@btime psos = poincaresos(ds, (2, 0.0), 1000.0, direction = +1)
  1.721 ms (2828 allocations: 253.23 KiB)

@btime  for i in 1:100000; step!(pmap); end
 6.384 ms (500000 allocations: 15.26 MiB)

PR: 
@btime psos = poincaresos(ds, (2, 0.0), 1000.0, direction = +1)
  1.718 ms (2828 allocations: 253.23 KiB)

@btime  for i in 1:100000; step!(pmap); end
  6.324 ms (500000 allocations: 15.26 MiB)

Thomas_Cyclical system is very sharp with pmap imagen

awage commented 3 years ago

Also I found a bug in check_next_state! that should be pushed in master. If not in this PR in another.

Datseris commented 3 years ago

Feel free to push the bugfix here.

@awage can you please provide a summary of the changes you did, and more importantly, why you did them? This will make it easier for me to review. By looking at the code, seems like the behavior of PoincareMap was changed, so it would be good to give a rationale.

awage commented 3 years ago

Hi! Here is a quick summary of the reason of the changes:

step! returns nothing when the trajectory does not return to the map after Tmax. But if get_state returns nothing we have a problem with view in get_projected_state that is expecting a vector.
The trick is to return a vector of NaN in get_state if the orbit does not return. This way we can check in _identify_basin_of_cell! if the orbit is off the plane (get_state is used to get the full state of the system).
This has no impact on the behavior of pmap since step! and poincaresos will use the low level function and return the interpolated projected state.

I think the fundamental problem is that step! does not behave in the same way for PoincaŕeMap and for ContinuousDynamicalSystems. In the second case you never get nothing and there is no type conflict.

Another solution is to change the behavior of step! in order to be consistent with get_state so that they return both NaN or something characteristic that makes clear that the trajectory does not return to the map.

Datseris commented 3 years ago

The tests we have now are unfeasible. Just the test of the poincare map on the Thomas system takes 6 minutes, which is unsustainable; we need to come up with smaller tests, which is what I am doing now. Once I have resolved this, we can merge. Unfortuantely, I found a really weird behavior. If I make the grid of Thomas coarser, the process fails:

    ds = Systems.thomas_cyclical(b = 0.1665)
    xg = yg = range(-6.0, 6.0; length = 50)
    pmap = poincaremap(ds, (3, 0.0); idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9)
    basin, attractors = basins_of_attraction((xg,yg), pmap; show_progress = true)

basin has -1 entries. In fact, all entries are -1.

EDIT: Forget what I said. I did something wrong which I will fix now. The current branch, as-is, works fine.

Datseris commented 3 years ago

BTW I think the current solution is just fine, I don't think it is more useful to make the return type a vector of NaN or something like that.

Datseris commented 3 years ago

Okay, I'm happy with the final version. One thing I found clunky and unintuitive is why the PoincareMap was having projection indices, given that it operates on the full system state. It made the remaining code confusing, and the interaction with the basins of attraction even more confusing. For example, the old code would find as attractors 3-dimensional points, even with the PoincareMap having idxs = 1:2. I thought that was weird. So I completely removed idxs from PoincareMap.

The other thing I did was change the logic of Tmax in PoincareMap. Previously, it would count time from the creation of the map. Now it instead counts time from the start of step!, which I think is more sensible and will better allow longer computation times without worrying whther the iteration exceeded Tmax even for a well-defined hyperplane.

awage commented 3 years ago

Awesome!!

I found an example to compare the output of the algorithm with a published result (Sprott again). I am preparing a PR with this system:

using DynamicalSystems
@inline @inbounds function dl_lorenz(u, p, t)
    R = p[1];
    x = u[1]; y = u[2]; z = u[3];
    dx = y - x
    dy = -x*z
    dz = x*y - R
    return SVector{3}(dx, dy, dz)
end
R=4.7;
ds = ContinuousDynamicalSystem(dl_lorenz, rand(3), [R])
xg = yg = range(-10.,10.,length=500)
pmap = poincaremap(ds, (3, 0.), Tmax=1e2)
basin,attractors = basins_of_attraction((xg,yg), pmap)
plot(xg,yg,basin',seriestype=:heatmap)

Original in J. C. Sprott, Simplest Chaotic Flows with Involutional Symmetries, Int. Jour. Bifurcation and Chaos 24, 1450009 (2014):

imagen

Computed with DynamicalSystems.jl (30 seconds for 500x500!!) imagen

JuliaDynamics / ChaosTools.jl

pmap get_state returns interpolated state #190