Instantaneous Lyapunov exponent for systems with parameter drift

[rusandris]

This is a non-autonomous version of Lyapunov exponent for systems with parameter drift, based on Dániel Jánosi, Tamás Tél, 2024 .

In non-autonomous systems with parameter drift, long time averages are less useful to assess chaoticity.

There are several different ways to address this, one is related to the ensemble approach.
Here, a new quantity called the Ensemble-averaged pairwise distance (EAPD) is used. An analog of the classical largest Lyapunov exponent can be extracted from the EAPD function ρ(t) : the local slope can be considered an instantaneous Lyapunov exponent.

This was previously discussed in #337 .

Example (low-precision version of Fig.11 a) from the article above):

using DynamicalSystems
using LinearAlgebra
using StatsBase
using Plots,LaTeXStrings


function duffing_drift(u0 = [0.1, 0.25]; ω = 1.0, β = 0.2, ε0 = 0.4, α=0.00045)
   return CoupledODEs(duffing_drift_rule, u0, [ω, β, ε0, α])
end

@inbounds function duffing_drift_rule(x, p, t)
           ω, β, ε0, α = p
           dx1 = x[2]
           dx2 = (ε0+α*t)*cos(ω*t) + x[1] - x[1]^3 - 2β * x[2]
           return SVector(dx1, dx2)
end

duffing = duffing_drift()
duffing_map = StroboscopicMap(duffing,2π)

init_states = randn(1000,2) #initalize ensemble
ρ,times = EAPD(duffing_map,init_states,100;Ttr=20,sliding_param_rate_index=4); #calc EAPD
lyapunov_instant(ρ,times;interval=1:20) #slope of first 20 time steps

pl = plot(times,ρ,ylabel=L"\rho(t)",lc=:gray10,lw=2,xlabel=L"t",legend=false,guidefontsize=20,dpi=300)

[Datseris]

Overall this looks fine, but cna you add the docstring and test before I do a proper review?

[Datseris]

can you also post the pictures as the artivcle is not open access?

[Datseris]

@rusandris is this PR still a draft or ready to review? Furthermore, are the test failures a result of this PR or unrelated?

[rusandris]

can you also post the pictures as the artivcle is not open access?

I'm not sure if I can just post the screenshots from the original article here, but anyways here are my version of Fig.11.

The slope values came out to be only approximately equal to the ones in the article (for an ensemble of N=10000), but trend-wise it seems to be okay.

[rusandris]

@rusandris is this PR still a draft or ready to review? Furthermore, are the test failures a result of this PR or unrelated?

I've added a test for the regular Henon map (autonomous) just to check if the instantaneous Lyapunov exponent coincides with the regular Lyap exp from the Benettin method (this is basically a sanity check that the time- and ensemble averages are the same).

Oops, I forgot to include my EAPD.jl tests in runtests.jl so the failures are definitely coming from somewhere else. I might need to look into those as well.

[Datseris]

The tests are good, but I am missing some clarity on what this function does.

I am also wondering why you do this "convergence to attractor" step first. Is it described in the paper?

[Datseris]

If the default value is normally distributed, why do we call it uniform? :D , maybe rename to perturbation_normal.

[rusandris]

That's funny, I guess it was uniform, then I changed it for some reason and forgot to rename the function

[Datseris]

Please make this a proper citation like in the docstrings of other functions of the library.

[Datseris]

Throughout this whole file, whenever there is a comma , it should be followed by a space. So all a,b,c needs to become a, b, c.

[Datseris]

There is a parameter index but there is no rate given as a keyword...? I am really confused. What does this keyword do?

[Datseris]

I don't understand what this function is supposed to do just by reading the docstring. We need to clarify some things: is the input ds expected to be autonomous or NON- ? And, is it mandatory that the system is in a parameter drift scenario, or there can be arbitrary non-autonomous behavior? If it can only be a drift, does the parameter drift happens inside ds or do you create it inside this new function? If it happens inside the ds, why do you need sliding_param_rate_index? And then you set this parameter to 0...?

Then inside the function I see this is done so that the ds reaches an attractor. But this is not said in the function.

I think this function is non-trivial and therefore warrants a ## Description section.

[rusandris]

You're right, I didn't explain the assumptions I'd made about the function itself.
I assumed that the parameter drift happens inside ds, meaning that the user already has a linearly drifting parameter in the definition. Let's take the system from above as an example here:

function duffing_drift(u0 = [0.1, 0.25]; ω = 1.0, β = 0.2, ε0 = 0.4, α=0.00045)
   return CoupledODEs(duffing_drift_rule, u0, [ω, β, ε0, α])
end

@inbounds function duffing_drift_rule(x, p, t)
           ω, β, ε0, α = p
           dx1 = x[2]
           dx2 = (ε0+α*t)*cos(ω*t) + x[1] - x[1]^3 - 2β * x[2]
           return SVector(dx1, dx2)
end

where ε parameter is replaced by the drifting term (ε0+α*t). I also assume that the drift rate α is in p. This is a simple way to accomplish drifting without extending the phase space.

I am also wondering why you do this "convergence to attractor" step first. Is it described in the paper?

Yes. You need some time to reach the snapshot attractor.
In practice, the convergence time is the time after which the averages of the phase space variables taken with respect to the ensembles are practically identical, although time dependent.
If init_states are randomly initialized (far from the attractor at the initial parameter), and there's no transient, the first few time steps cannot be used to calculate any reliable averages.

Because of this, there is a Ttr transient time option, which is spent on converging to the attractor at ε0, before the drift even starts happening. In other words, the drift needs to be switched off (drift rate α = 0) for this step. In order to switch it off, I need to know which parameter I need to set to 0, this is why sliding_param_rate_index
exists. This was the simplest way I could think of, but maybe there's a better way.

Now, what if you want to use this function in an autonomous setting (like I did in the tests with the Henon system)? Then you don't have a drift rate among the parameters, there is no need to switch off anything, in which case the default is sliding_param_rate_index=0 and the function still works.
Let me know if this makes sense.

As far as I know there is no standard way to create a drifting system. (Maybe here but I'm not sure.)

[Datseris]

It appears to me that CriticalTransitions.jl discusses a "drift" system in the opening image but nowhere in the stable docs is this mentioned or implemented. So we can safely ignore that.

Right, your approach is fine with me: we expect a parameter index pidx (mandatory input argument, and not optional keyword argument), and simply explain how this process works. that it is expected that p[pidx] multiplies time t and it is the "rate" parameter. For the non-drifiting henon map test just make a 3-parameter henon map where the 3rd parmaeter is not used anywhere. Then give pidx = 3.

Do you mind going ahead with adding the descriptions etc in the source code/docstring?

[Datseris]

cc @reykboerner or @ryandeeley I don't know if you are still working on CriticalTransitions.jl and what is the status for the drifting systems.

[Datseris]

Oh and @rusandris you shouldn't mess with the whole parameter container. Just do original_rate = current_parameter(ds, pidx) and then after setting it to zero do set_parameter!(ds, pidx, original_rate) to restore it.

[ryandeeley]

cc @reykboerner or @ryandeeley I don't know if you are still working on CriticalTransitions.jl and what is the status for the drifting systems.

Hi, the status regarding "drifting systems" is that @raphael-roemer and I have developed code for passing a user-defined ds::CoupledODEs (describing the system of interest in an autonomous form) and a rp::RateProtocol (describing the desired time-dependent forcing, or "drifting", one wishes to apply to the system) to a function RateSystem(), which returns a CoupledODEs describing the corresponding non-autonomous system. In particular the function RateSystem() creates a modified drift function f(u,p,t) from the information of ds.integ.f and rp. The code is available in the RateSystem branch of CriticalTransitions.jl (https://github.com/JuliaDynamics/CriticalTransitions.jl/tree/RateSystem), specifically in the file /test/ratesystem/RateSystemDraft2.jl, and the discussion surrounding this can be found in this pull request: JuliaDynamics/CriticalTransitions.jl#117. Documentation hasn't been written yet, I suppose in an ideal world that'll be done soon. I'm happy to help with anything if you'd like to implement this and cannot understand things just from the scripts.

[Datseris]

Thanks, that's good to know. @rusandris does this method you are implementing rely on a linear drift? I.e., it mandates a term p*t with t time? Or time dependence can be arbitrary?

For now we will proceed as is in this PR with pidx. In the future, when the drifting system is available as public API we can use it instead.

[rusandris]

Yes, it relies on linear drift, but it can be generalized further in the future if needed.