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]

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.

[rusandris]

Tests seem to be ok, the error that keeps popping up comes from elsewhere:

Lorenz stable: Test Failed at /home/runner/work/ChaosTools.jl/ChaosTools.jl/test/chaosdetection/lyapunovs.jl:67
  Expression: ≈(lyapunov(ds, 1000, Ttr = 100), 0, atol = 0.001)
   Evaluated: -0.0015200001782929603 ≈ 0 (atol=0.001)

I suggest this needs to be dealt with in a separate PR.

[Datseris]

Thanks @rusandris ; Let's finish this PR then! You need to:

1. Add the new method to the documentation.
2. Add a changelog entry.
3. Increase the minor version of the Project.toml file!

[codecov-commenter]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 62.21%. Comparing base (928bdb0) to head (9dca714).
⚠️ Report is 15 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #345      +/-   ##
==========================================
+ Coverage   54.86%   62.21%   +7.35%     
==========================================
  Files          22       26       +4     
  Lines         904     1064     +160     
==========================================
+ Hits          496      662     +166     
+ Misses        408      402       -6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Datseris]

dimensionless? why? it has the dimensions of the state space.

[rusandris]

It was meant in a sense that your variables in your state space vectors need to be dimensionless in order the distance to make sense. But I don't think it needs that theoretical detail there, so I took it out for now.

[rusandris]

Thanks @rusandris ; Let's finish this PR then! You need to:

1. Add the new method to the documentation.

I've added to lyapunovs.md, and it looks fine for a first version (I managed to build it locally), but I couldn't figure out how to render paper references properly yet

[Datseris]

Sorry for being late, I had some personal life problems. I am thinking whether it is possible to combine this with the RateSystem of JuliaDynamics/CriticalTransitions.jl#117 which creates a RateSystem that could be used directly here. But to do this, I must first ask @rusandris : does this functionality here only work / make sense if the forcing is linear? Or the methodology would work with any kind of forcing?

[rusandris]

does this functionality here only work / make sense if the forcing is linear? Or the methodology would work with any kind of forcing?

In the original review paper the parameter drift was chosen to be linear for simplicity. The method itself - calculating the ensemble averaged distance between trajectory pairs - doesn't depend on this linearity. You just have to be careful when taking the slope of ρ(t), as it only can be considered as Lyapunov exponent if the distance between pairs is typically small. And that depends on your step size, not the functional form of the drift. I might be wrong but that's how I understand it right now.

Looking at the code, the function ensemble_averaged_pairwise_distance does use the assumption that the drift is linear. Particularly, pidx parameter is for the index of the drifting rate. But that can be generalized if needed, especially if we want to make this compatible with the RateSystem.