search

JuliaDynamics / ChaosTools.jl

Tools for the exploration of chaos and nonlinear dynamics
https://juliadynamics.github.io/DynamicalSystemsDocs.jl/chaostools/stable/
MIT License
187 stars 35 forks source link

Support for very high-dimensional systems in `lyapunovspectrum`/ usage of DDEs #255

Closed ikottlarz closed 2 years ago

ikottlarz commented 2 years ago

Issue description/ description of requested feature

I would like to use DynamicalSystems with a delay differential equation, specifically the Mackey-Glass equation. I'd like to calculate the Lyapunov spectrum of the system. DEEs are not supported in DynamicalSystems. So currently, I'm doing this by defining a DiscreteDynamicalSystem, in which I'm doing a "manual" integration step and keeping track of the past of the system in a vector. I can use this to calculate the Lyapunov spectrum for up to a system dimensionality of 90 (corresponding to a time delay of 9 with an integration step of 0.1 -> saving all steps up to 1 delay time). If I go larger, I get 

ERROR: LoadError: syntax: invalid syntax (memory-error out of gc handles)                                                                                                                                      
Stacktrace:                                                                                                                                                                                                    
  [1] top-level scope                                                                                                                                                                                          
    @ ~/Documents/master/complexity_entropy/test.jl:29                                                                                                                                                         
  [2] macro expansion                                                                                                                                                                                          
    @ ~/.julia/packages/StaticArrays/0bweZ/src/SMatrix.jl:35 [inlined]                                                                                                                                         
  [3] SArray                                                                                                                                                                                                   
    @ ~/.julia/packages/StaticArrays/0bweZ/src/SMatrix.jl:32 [inlined]                                                                                                                                         
  [4] macro expansion                                                                                                                                                                                          
    @ ~/.julia/packages/StaticArrays/0bweZ/src/arraymath.jl:51 [inlined]                                                                                                                                       
  [5] _rand                                                                                                                                                                                                    
    @ ~/.julia/packages/StaticArrays/0bweZ/src/arraymath.jl:43 [inlined]                                                                                                                                       
  [6] rand                                                                                                                                                                                                     
    @ ~/.julia/packages/StaticArrays/0bweZ/src/arraymath.jl:41 [inlined]                                                                                                                                       
  [7] rand(rng::Random.MersenneTwister, ::Type{SMatrix{100, 100, T, L} where {T, L}})                                                                                                                          
    @ Random /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Random/src/Random.jl:256                                                                                                    
  [8] rand(#unused#::Type{SMatrix{100, 100, T, L} where {T, L}})                                                                                                                                               
    @ Random /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Random/src/Random.jl:259                                                                                                    
  [9] orthonormal                                                                                                                                                                                              
    @ ~/.julia/packages/DelayEmbeddings/f4iJo/src/utils.jl:117 [inlined]                                                                                                                                       
 [10] lyapunovspectrum(ds::DiscreteDynamicalSystem{true, Vector{Float64}, 100, typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, DynamicalSystemsBase.var"#6#12"{ForwardDiff.JacobianConfig{Forwar\
dDiff.Tag{DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Float64}, Float64, 12, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.var"\
#5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int\
64, Int64, Float64, Float64}, Int64}, Float64}, Float64, 12}}}}, DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Vector{Float64}}, Matrix{Float64}, true\
}, N::Int64, k::Int64; kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:Ttr, :Δt), Tuple{Int64, Int64}}})                                                                       
    @ ChaosTools ~/.julia/packages/ChaosTools/Xe6aI/src/chaosdetection/lyapunovs.jl:70                                                                                                                         
 [11] top-level scope                                                                                                                                                                                          
    @ ~/Documents/master/complexity_entropy/test.jl:29                                                                                                                                                         
in expression starting at /home/inga/Documents/master/complexity_entropy/test.jl:29

From what I understand, this is due to the use of StaticArrays in the calculation of the spectrum, where

A very rough rule of thumb is that you should consider using a normal Array for arrays larger than 100 elements. ( Readme of StaticArrays.jl)

Maybe the usage of Base.Arrays could be allowed for large Systems? If not, is there another possible workaround I'm not seeing right now?

Also, are there any plans on maybe making a DelayDynamicalSystem or such available for DynamicalSystems?

MWE

This works fine for τ < 10 and throws the error stated above for larger values

using DynamicalSystems

function mackey_glass!(un, u, p, t)
    β, γ, n, Δt = p
    L = length(u)
    # last entry is x(t-τ)
    x_τ = u[L]
    # RK4
    k1 = β * x_τ / (1+x_τ^n) - γ*u[1]
    k2 = β * x_τ / (1+x_τ^n) - γ*(u[1]+0.5*k1)
    k3 = β * x_τ / (1+x_τ^n) - γ*(u[1]+0.5*k2)
    k4 = β * x_τ / (1+x_τ^n) - γ*(u[1]+k3)
    un[1] = u[1] + Δt * (k1/6 + k2/3 + k3/3 + k4/6)
    # keep history in memory
    un[2:L] = u[1:L-1] 
    return 
end

β = 2
γ = 1
n = 9.65
Δt = 0.1
τ = 10
N = 200000
u0 = zeros(Int(τ/Δt))
u0[1] = 1.
p = β, γ, n, Δt
ds = DiscreteDynamicalSystem(mackey_glass!, u0, p)
Lambdas = lyapunovspectrum(ds, N; Ttr=Int(1000/Δt), Δt=10)
Datseris commented 2 years ago

well, isn't that a face I haven't seen in a while :) Hi Inga!

Whether or not static arrays are used is user's choice. The algorithm uses normal Arrays if the system is "in-place". So something must be going wrong here, maybe a bug is at place. I'll try to identify it.

In the meantime, can you try again with providing a hard-coded Jacobian? Its easy here to write one analytically, and the automatic jacobian has strong drawbacks in very high dimensions.

Also, are there any plans on maybe making a DelayDynamicalSystem or such available for DynamicalSystems?

Well, sure, if one finds algorithms for them and decides to take the task of doing the PR that implements them ;) DDEs are "allowed" in DynamicalSystems.jl, they aren't just implemented. Reason is simple: when I started writing the library I didn't have any algorithms that really could use DDEs. (and as you have already realized above, you dont use the lyapunov algorithms for DDEs, you use it for discrete systems)

However, all of this can will change. After we implemented new types of dynamical systems, like stroboscopic maps, I realized DynamicalSystems.jl can benefit from yet another re-design. This is discussed in https://github.com/JuliaDynamics/DynamicalSystems.jl/issues/176 In this redesign, "dynamical systems" are just generic integrators. Many algorithms only require that, without caring if the integrator is a discrete or continuous, or perhaps even a delay system. Some algorithms do care. In any case, sure, DDEs are allowed, its just a matter of putting in the work. The redesign will make this work much less, but the redesign itself is some work :D

ikottlarz commented 2 years ago

Hi George :) When providing a jacobian, this works fine for τ >= 10*, so that's great, thanks a lot!

*EDIT: This works fine for τ >= 10 if I only calculate the k<90 largest lyapunov exponents, which I couldn't do before (see below). When trying to calculate all exponents, I get the same error as before. But this is fine for me since I now have a working solution for my specific problem.

I made another observation: When I tried to calculate only the k=20 largest lyapunov exponents using

Lambdas = lyapunovspectrum(ds, N, k; Ttr=Int(1000/Δt), Δt=10)

I got a DimensionMismatch for k != dimension(ds) when I did not provide a jacobian (here for a 30-dimensional system):

ERROR: LoadError: DimensionMismatch("tried to assign 30×30 array to 30×20 destination")                                                                                                                        
Stacktrace:                                                                                                                                                                                                    
  [1] throw_setindex_mismatch(X::LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}, I::Tuple{Int64, Int64})                                                                                                 
    @ Base ./indices.jl:193                                                                                                                                                                                    
  [2] setindex_shape_check                                                                                                                                                                                     
    @ ./indices.jl:253 [inlined]                                                                                                                                                                               
  [3] _unsafe_setindex!(::IndexLinear, ::Matrix{Float64}, ::LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}}, ::Base.Slice{Base.OneTo{Int64}}, ::UnitRange{Int64})                                         
    @ Base ./multidimensional.jl:896                                                                                                                                                                           
  [4] _setindex!                                                                                                                                                                                               
    @ ./multidimensional.jl:887 [inlined]                                                                                                                                                                      
  [5] setindex!                                                                                                                                                                                                
    @ ./abstractarray.jl:1267 [inlined]                                                                                                                                                                        
  [6] set_deviations!                                                                                                                                                                                          
    @ ~/.julia/packages/DynamicalSystemsBase/jey1Q/src/core/discrete.jl:53 [inlined]                                                                                                                           
  [7] lyapunovspectrum(integ::DynamicalSystemsBase.MinimalDiscreteIntegrator{true, Matrix{Float64}, 630, DynamicalSystemsBase.var"#20#23"{20, typeof(mackey_glass!), Matrix{Float64}, ForwardDiff.JacobianConf\
ig{ForwardDiff.Tag{DynamicalSystemsBase.var"#19#22"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}}, Float64}, Float64, 10, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.v\
ar"#19#22"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}}, Float64}, Float64, 10}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.var"#19#22"{typeof(mackey_glass!), Tuple{Int64\
, Int64, Float64, Float64}}, Float64}, Float64, 10}}}}}, Tuple{Int64, Int64, Float64, Float64}}, N::Int64, Δt::Int64, Ttr::Int64, show_progress::Bool)                                                         
    @ ChaosTools ~/.julia/packages/ChaosTools/Xe6aI/src/chaosdetection/lyapunovs.jl:103                                                                                                                        
  [8] lyapunovspectrum(ds::DiscreteDynamicalSystem{true, Vector{Float64}, 30, typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, DynamicalSystemsBase.var"#6#12"{ForwardDiff.JacobianConfig{Forward\
Diff.Tag{DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Float64}, Float64, 10, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.var"#\
5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Float64}, Float64, 10}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int6\
4, Int64, Float64, Float64}, Int64}, Float64}, Float64, 10}}}}, DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Vector{Float64}}, Matrix{Float64}, true}\
, N::Int64, Q0::SMatrix{30, 20, Float64, 600}; Ttr::Int64, Δt::Int64, u0::Vector{Float64}, show_progress::Bool, diffeq::NamedTuple{(), Tuple{}}, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, Named\
Tuple{(), Tuple{}}})                                                                                                                                                                                           
    @ ChaosTools ~/.julia/packages/ChaosTools/Xe6aI/src/chaosdetection/lyapunovs.jl:89                                                                                                                         
  [9] lyapunovspectrum(ds::DiscreteDynamicalSystem{true, Vector{Float64}, 30, typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, DynamicalSystemsBase.var"#6#12"{ForwardDiff.JacobianConfig{Forward\
Diff.Tag{DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Float64}, Float64, 10, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.var"#\
5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Float64}, Float64, 10}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int6\
4, Int64, Float64, Float64}, Int64}, Float64}, Float64, 10}}}}, DynamicalSystemsBase.var"#5#11"{typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, Int64}, Vector{Float64}}, Matrix{Float64}, true}\
, N::Int64, k::Int64; kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:Ttr, :Δt), Tuple{Int64, Int64}}})                                                                        
    @ ChaosTools ~/.julia/packages/ChaosTools/Xe6aI/src/chaosdetection/lyapunovs.jl:70                                                                                                                         
 [10] macro expansion                                                                                                                                                                                          
    @ ./show.jl:955 [inlined]                                                                                                                                                                                  
 [11] top-level scope                                                                                                                                                                                          
    @ ~/Documents/master/complexity_entropy/src/complexity_entropy/pipelines/data_generation/mackey_glass.jl:82                                                                                                
in expression starting at /home/inga/Documents/master/complexity_entropy/src/complexity_entropy/pipelines/data_generation/mackey_glass.jl:66

~Maybe both problems are related to ForwardDiff?~ Since the first problem still occurs with a "manual" jacobian, this is probably not related to ForwardDiff

Well, sure, if one finds algorithms for them and decides to take the task of doing the PR that implements them ;)

I was (naively?) thinking that maybe one could write something that generalizes what I've done, just using a DiscreteDynamicalSystem? Although that of course wouldn't use any of the nice and sophisticated integrators in DifferentialEquations.jl, and probably has some other drawbacks that I'm not seeing here since I don't know that much about DDEs.

Datseris commented 2 years ago

this sounds definitely like a bug, i'll try to pinpoint it when I have more time and better vision. In the meantime if you find the source of it of course please do post it here.

awage commented 2 years ago

Hi @ikottlarz , I've seen that your using a discretization of your delay equation. I wonder if you really need a RK4 scheme, if you look closely at your system you can factor a lot of terms, it reduces to: 

function mackey_glass!(un, u, p, t)
    β, γ, n, Δt = p
    L = length(u)
    # last entry is x(t-τ)
    x_τ = u[L]
    # RK4
    K = (β * x_τ / (1+x_τ^n)) - γ*u[1]
    un[1] = u[1] + Δt * K * SOME CONSTANT DEPENDING ON γ
    # keep history in memory
    un[2:L] = u[1:L-1] 
    return 
end

That will speed up considerably your simulations.

ikottlarz commented 2 years ago

Hi @awage, thank you for your comment!

One thing: Actually, the expression I posted above is not entirely correct. The correct expression would be

[...]
    k1 = Δt * (β * x_τ / (1+x_τ^n) - γ*u[1])
    k2 = Δt * (β * x_τ / (1+x_τ^n) - γ*(u[1]+0.5*k1))
    k3 = Δt * (β * x_τ / (1+x_τ^n) - γ*(u[1]+0.5*k2))
    k4 = Δt * (β * x_τ / (1+x_τ^n) - γ*(u[1]+k3))
    un[1] = u[1] + k1/6 + k2/3 + k3/3 + k4/6

or, equivalently

    k1 = β * x_τ / (1+x_τ^n) - γ*u[1]
    k2 = β * x_τ / (1+x_τ^n) - γ*(u[1]+0.5*Δt*k1)
    k3 = β * x_τ / (1+x_τ^n) - γ*(u[1]+0.5*Δt*k2)
    k4 = β * x_τ / (1+x_τ^n) - γ*(u[1]+Δt*k3)
    un[1] = u[1] + Δt * (k1/6 + k2/3 + k3/3 + k4/6)

(I had originally used the first version, then, while posting here, thought "oh right, why don't I factor out the Δt?", not realizing that I had also put it into the k_i expressions, then found while this also integrates the system without error, it is not a MG time series)

But in any case, you're right, I can factor out terms. I did this and am now left with

function mackey_glass!(un, u, p, t)
    β, γ, n, Δt = p
    L = length(u)
    # last entry is x(t-τ)
    x_τ = u[L]
    # RK4 with factorizing of terms
    K = β * x_τ / (1+x_τ^n) - γ*u[1]
    C = 0.5*Δt*γ
    un[1] = u[1] + Δt/6 * K * (6 - 2C * (3 - C * (2 - C)))
    # keep history in memory
    un[2:L] = u[1:L-1] 
    return 
end

Although it seems to me (without having actually quantified it yet) that this version is not that much faster than what I had before.

awage commented 2 years ago

The performance bottleneck is probably in the Jacobian evaluation then.

Off topic: since all variables are delayed version of the first one, what do the lyapunov spectrum look like? You need all the exponents?

ikottlarz commented 2 years ago

The performance bottleneck is probably in the Jacobian evaluation then.

Yeah I think so too. It did get significantly faster once I provided a Jacobian myself.

Off topic: since all variables are delayed version of the first one, what do the lyapunov spectrum look like? You need all the exponents?

I'm mostly interested in the KY-dimension of the system and for the parameter's I'm using, D_KY~τ (at least for the delays I've considered so far). Right now, I'm calculating only the k=τ+5 largest exponents, which seems to suffice for the delays I've considered. Here's a plot of the spectrum and the KY-dim: image

ikottlarz commented 2 years ago

Update: I was originally going for delays up to 25. If I go even larger, the code breaks again at τ = 27 with the original error, but interestingly, with a different stack trace: 

...
τ = 27
k = τ+5
Lambdas = lyapunovspectrum(ds, T, k; Ttr=Int(1000/Δt), Δt=10)

yields

ERROR: LoadError: syntax: invalid syntax (memory-error out of gc handles)
Stacktrace:
  [1] top-level scope
    @ ~/Documents/master/complexity_entropy/src/complexity_entropy/pipelines/data_generation/mackey_glass.jl:64
  [2] macro expansion
    @ ~/.julia/packages/StaticArrays/NQjQM/src/SMatrix.jl:35 [inlined]
  [3] SArray
    @ ~/.julia/packages/StaticArrays/NQjQM/src/SMatrix.jl:32 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/StaticArrays/NQjQM/src/arraymath.jl:51 [inlined]
  [5] _rand
    @ ~/.julia/packages/StaticArrays/NQjQM/src/arraymath.jl:43 [inlined]
  [6] rand
    @ ~/.julia/packages/StaticArrays/NQjQM/src/arraymath.jl:41 [inlined]
  [7] rand(rng::Random.MersenneTwister, ::Type{SMatrix{270, 32, T, L} where {T, L}})
    @ Random /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Random/src/Random.jl:256
  [8] rand(#unused#::Type{SMatrix{270, 32, T, L} where {T, L}})
    @ Random /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Random/src/Random.jl:259
  [9] orthonormal
    @ ~/.julia/packages/DelayEmbeddings/f4iJo/src/utils.jl:117 [inlined]
 [10] lyapunovspectrum(ds::DiscreteDynamicalSystem{true, Vector{Float64}, 270, typeof(mackey_glass!), Tuple{Int64, Int64, Float64, Float64}, typeof(mackey_glass_jac!), Matrix{Float64}, false}, N::Int64, k::Int64; kwargs::Base.Iterators.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:Ttr, :Δt), Tuple{Int64, Int64}}})
    @ ChaosTools ~/.julia/packages/ChaosTools/Xe6aI/src/chaosdetection/lyapunovs.jl:70
 [11] top-level scope
    @ ~/Documents/master/complexity_entropy/src/complexity_entropy/pipelines/data_generation/mackey_glass.jl:88
in expression starting at /home/inga/Documents/master/complexity_entropy/src/complexity_entropy/pipelines/data_generation/mackey_glass.jl:64

So, I think this one is due to the orthonormal function, in which StaticArrays are always used. It says in the docstring 

"""
    orthonormal(D, k) -> ws
Return a matrix `ws` with `k` columns, each being
an `D`-dimensional orthonormal vector.

Always return `SMatrix` for stability reasons.
"""

I don't understand enough about these stability reasons to know if it would be a good idea to allow for other types, or I should find some other solution.

Datseris commented 2 years ago

oh shiiiiiiiiiiiiiiiiiiiiiiiiiiiit

yeah, yeah, you are right, this is the bug. Thanks so much for finding it out!!!!! Yes, we need to modify this function into the version

orthonormal([T,] D, k)

with T deciding the return type Matrix or SMatrix. Then it returns normal matrices. This version should be used in the inplace format always actually. Daaaaaaaaaamn this is a blast from the past, I wrote this function literally 5 years ago.

ikottlarz commented 2 years ago

Alright, I made a pull request in DelayEmbeddings.jl which I think would solve this