re-open PR #83

[Datseris]

We continue here directly from #83 , I made a mistake and merged the previous one...

Closes #68

[Datseris]

@efosong I've added you to a team in JuliaDynamics so you can push at ChaosTools (and thus this PR directly). Sorry for the mistake but there is some work to be done here before you can get the issue bounty. Writing a docstring is an absolute must for example (as well as passing tests, which I am working on now).

[Datseris]

Okay, now tests pass (at least locally on my machine). They take quite some time, but this is expected given how many computations are necessary for the result.

I will spend some time implementing tests for the logistic map until we get the documentation pushed.

[efosong]

Thanks @Datseris 🙂

I'll write some documentation some time during the week.

I tried to write some discrete tests based on the Hénon map (as well as the "standard map") . It seemed to require a rather large number of iterations, which seems reasonable for discrete time systems, but I didn't really get any results I was happy with when I plotted them.

In some cases, I got convergence which I guess invalidates the 'limit cycle' based analysis in the paper so the results were a bit weird. I should've written this down what I observed when I was testing these, but didn't.

[Datseris]

@efosong can you please write down exactly what parameters you used for the Henon tests? For the standard map the method we are implementing here doesn't work. The standard map does not have attracting sets.

[Datseris]

I am testing with the henon map, whose behavior is pretty much established in the literature. The parameter b = 0.3 is fixed. The parameter a is varied which gives the following orbit diagram:

Now I run the code:

ds = Systems.henon()
a_reg = 0.8; a_cha = 1.4; a_reg2 = 1.0

for a in (a_reg, a_reg2, a_cha)
    set_parameter!(ds, 1, a)
    chaos_type, ν, C = predictability(ds; T_max = 40000)
    println("type = ", chaos_type, ". ν, C = ", round.((ν, C);digits=3))
end

(on the branch of this PR) and I get the results:

type = INDETERMINATE. ν, C = (NaN, 1.0)
type = PPC. ν, C = (0.013, 1.0)
type = SC. ν, C = (0.002, -0.029)

So the algorithm fails completely in the case of regular behavior. The NaN result shows that there has to be some problem with the implementation of the algorithm. But the result is also that a period-4 motion is "Partially Predictable Chaos" instead of "Laminar" behavior. This can't be correct, right?

@hwernecke when coming up with the algorithm, have you tested it with discrete systems? Is it our implementation that is at fault, or does it in general not perform well for discrete systems?

[hwernecke]

So the algorithm fails completely in the case of regular behavior. The NaN result shows that there has to be some problem with the implementation of the algorithm. But the result is also that a period-4 motion is "Partially Predictable Chaos" instead of "Laminar" behavior. This can't be correct, right?

@hwernecke when coming up with the algorithm, have you tested it with discrete systems? Is it our implementation that is at fault, or does it in general not perform well for discrete systems?

When developing the test we focused on continuous systems (ODE) with periodic and chaotic motion and only checked fixed point and discrete systems (maps) briefly. But now that you report the issue I remember what the problem is (sorry, could have happened earlier...).

Discrete Systems
In discrete systems there is basically two types of attractors:

• chaotic attractor: the test works as expected
• fixed points: here we encounter the problem
  Periodic motion in maps (equivalent to laminar flow in ODE) is a fixed point of higher order, i.e. the attracting state consists of several points (finite number!) that are approached and then visited one after the other.

Fixed points
That brings me back to fixed points in continuous systems. For the test fixed points in maps and ODE are equivalent. And for both of these the current implementation of the test will fail.
BUT: don't worry, once one understands, why the test fails it is rather easy to check for fixed points.

At a fixed point (and I only consider stable, i.e. attracting ones) the dynamics is characterized by a negative Lyapunov exponents, the largest of which is λmax<0. The distance d between two trajectories that have an initial distance δ can be approximated by the exponential evolution d(t) = δ exp(λmax t).
Thus, in the long-term limit t → ∞ the distance vanishes d → 0, which means that the fixed point is attracting/stable.

For the test that means:

• If one waits long enough before measuring the distance, one finds d = 0 (in theory) for arbitrary initial distance. In practice the accuracy is limit due to numerics and thus one finds a small finite distance close to the boundary of numerical precision, e.g. d ~ 10^{-19}.
  From that one would conclude that the scaling exponent vanishes ν = 0.
• If the convergence to the fixed point is slow or the evaluation of distance happens before having reached the boundary of numerical precision, one will most likely observe a seemingly linear scaling ν = 1.

Two ways out:

• One uses a separate fixed point detection, e.g. by checking whether the distance is close to numerical accuracy (works well in continuous systems).
• One detects a clearly negative largest Lyapunov exponent (works well in discrete systems).

I hope that helps to overcome the misery of fixed points. And again: I am sorry that I completely forgot about this issue...

[Datseris]

@hwernecke perfect, thank you very much for your reply. I am happy because this is also what I have been thinking as the answer to this "problem".

From my perspective the best way forward is to simply deduce the result for discrete systems based on $λmax$. Let's say if $λmax < -1e-2$ we return immediately LAM as a result. In general I am fine with this since I believe your algorithm is more useful for distinguishing partially predictable chaos from strong chaos. Identifying laminar phases is much easier and there are many ways for a researcher to confirm that their system is laminar.

@efosong I will take care of the internals of the function to handle discrete systems or very small exponents. In the meantime, can you please work on the docstring? There is also one more thing remaining as a test: to bring in one more continuous system to check. We could use either of the systems Hendrik suggested in his comments in #83 , or we could also use the second system that Hendrik uses in the paper (which is a bit simpler than the neural model).

Then this PR is done!

[efosong]

@Datseris I've added a docstring. I tried to copy the style of some of the other docstrings. Please can you take a look and tell me if you think it's clear/any changes need to be made?

I'm also trying to implement the 3 neuron network in Wernecke et al 2018. I think I have it working, but I think my choice of initial conditions is quite different to the paper and perhaps causing it to converge to a different attractor. @hwernecke please can you let me know how you initialised x and b in this paper? In particular I'm trying to replicate figure 4 and get qualitatively similar behaviour but slightly different results --- my initialisation is random, but if you chose the initial conditions deliberately it would be good to know so I can be more confident in my comparison.

I could put a pull request in to add this model as a predefined system (with an arbitrary number of neurons) if you'd like @Datseris.

[Datseris]

@Datseris I've added a docstring. I tried to copy the style of some of the other docstrings. Please can you take a look and tell me if you think it's clear/any changes need to be made?

Awesome. Sure, I'm on it!

I could put a pull request in to add this model as a predefined system (with an arbitrary number of neurons) if you'd like @Datseris.

Definitely, seems like a plan!

[Datseris]

@efosong you have done a wonderful job here. The code is clear, understandable and transparent. The documentation string is excellent (I only did some minor changes, mostly to account for the above discussion; we return LAM for negative exponents), and the method seems to work fine for the Lorenz tests.

Your contribution is now certainly enough for you to get the bounty. Of course, if you want to contribute tests for the neural system, this will be welcomed. But I just want to point out that you don't have to if you don't want to! Let me know if you still want to do this, otherwise I'll merge.

(I am now writing tests for the Henon map!)

[Datseris]

Just a word of caution: I've changed "laminar" to "regular" throughout. In every other place in DynamicalSystems.jl documentation I use the word regular instead of laminar, so it is nice to have it fitting more.

[efosong]

Thanks @Datseris and @hwernecke for helping with this (my first open source) PR.

@Datseris: I think I might add tests for this neuron model at some point, but perhaps that can be added in a separate PR some time in the future so I can get this PR checked off my Todo list. If you're happy to merge the pull request now, then please do it 🙂

I ran into a little trouble with using my model with a parallel_integrator (for some reason, when I pass a vector of states type Array{Array{Float64, 1}, 1}, it interprets it as a matrix initial condition instead) and unless it's an easy fix I don't think I have time right now to dig into it.

At some point in the future I also hope to add some parallelisation as @SebastianM-C suggested, but this will come in another PR.

I might also take a look at some of the other issues at some point (e.g. #77, #84) and hopefully I'll be able to contribute a bit more independently for those ones!

[Datseris]

@efosong I am very happy to see that you are considering to contribute more!

I ran into a little trouble with using my model with a parallel_integrator (for some reason, when I pass a vector of states type Array{Array{Float64, 1}, 1}, it interprets it as a matrix initial condition instead) and unless it's an easy fix I don't think I have time right now to dig into it.

If you open an issue at DynamicalSystemsBase with a MWE (or if you PR the system you defined you can mention the error there) I can help you with it.

[Datseris]

@efosong don't forget to request the bounty!

[hwernecke]

@hwernecke please can you let me know how you initialised x and b in this paper? In particular I'm trying to replicate figure 4 and get qualitatively similar behaviour but slightly different results --- my initialisation is random, but if you chose the initial conditions deliberately it would be good to know so I can be more confident in my comparison.

@efosong Your plot looks great!
Indeed for the choice of parameters you used there are two coexisting attractors. In Fig. 4 of the paper we presented a trajectory from the other one. But you can compare your results with Fig. 6(a) panel (iii), where trajectories from both attractors are shown.
If you wish to find the other attractor, you have to try different initial conditions. We never checked the basin of attraction carefully.