Overhaul ContinuousDS

[Datseris]

EDIT: THIS IS THE NEW AND CURRENT PLAN, IN PROGRESS BY ME IN A BRANCH.

Example at the end!

• Rename ContinuousDS to ContinuousDynamicalSystem with alias CDS.

• CDS will have 2 fields: prob::ODEProblem and jacobian (which is a function).

  • At a later point we might consider having tangent dynamics function as part of it as well.

• Define two functions: makeintegrator(CDS) and maketangentintegrator(CDS) (names to be decided), that return integrators about the system or about the system+tangent space.

• Internally, all functions of DynamicalSystemsBase.jl and ChaosTools.jl (and any other packages of the ecosystem) will use step!(integrator, Δt). An argument like dt will be given to functions that actually care about dt.

• Define a function step!(integrator, Δt) that steps until a span of Δt has passed.

• All high-level functions and internals for both DynamicalSystemsBase.jl and ChaosTools.jl will not care about the tspan property of ODEProblem.

• Be sure that current time and parameters is passed properly into the Jacobian function from ForwardDiff.

• Define type aliases for CDS and DDS, e.g. CDS{ODE, A,B,C} = DS{ODE, A, B, C} where {ODE<:ODEProblem, A, B, C}. This will reduce integrator constructors by half!

• Provide high level constructors CDS(eom, state, p [, jacobian]) and CDS(prob [, jacobian]). The default tspan will be made (0.0, Inf).

• Parameterize CDS with IIP (already part of the ODEProblem).

• All of the above are meant to work for 2 cases: both in-place (iip) and out-of-place versions. This is a lot of work but the performance benefits of SVector for small systems are worth it.

• Create an oop Lorenz63 model and implement tests, similar to the tests here: https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/blob/master/test/discrete_types.jl

• Create a ParallelIntegrator (this has already started in the source code) that evolves pairs of vectors of states.

• All high level functions that want to use a CDS, will first call makeintegrator(CDS) and immediatelly pass that to a lower-level function that uses the integrator and step!. This is very important for parallelization stuff.

Here is the struct:

struct ContinuousDynamicalSystem{IIP, IAD, ODE<:ODEProblem, JAC}
    prob::ODEProblem
    jacobian::JAC
end

where IIP is "isinplace" and IAD is "is auto differentiated".

If the equations for the system are out of place, the Jacobian must be as well. The same for in-place. Supporting all 4 possible variants is a lot of work that I cannot do.

EDIT MORE: Try to combine DISCRETE + Continuous into DynamicalSystem, parameterized by the ODEProblem type. Meaning:

struct DynamicalSystem{
        DEP<:DEProblem, # problem type also gives dispatch purposes
        IIP, # is in place , for dispatch purposes
        JAC, # jacobian function (either user-provided or FD)
        IAD} # is auto differentiated? Only for constructing TangentEvolver
    prob::DEP
    jacobian::JAC
end

[Datseris]

@ChrisRackauckas if you have any suggestions that will make the above easier / etc. please go ahead!

[ChrisRackauckas]

Nah, that looks like how I'd do it. You do have to code a lot of paths in two ways to account for this. Overall, the out-of-place form is usually better for what you're doing I'd think (<100 ODEs).

[Datseris]

Alright, I wish somebody will want to work on this :) I think by viewing how I did it in #23 it can make things very clear!

[Datseris]

@ChrisRackauckas My plan here is to make a continuous dynamical system have a ODEIntegrator as a field instead of a ODEProblem. The reason is that it is much much faster to just reinit! an integrator than to create a new one each time!!!

What do you think?

[ChrisRackauckas]

Well, that's only true for in-place systems. If you're using out-of-place, the cache is essentially nothing. YMMV but that's something to think about.

[tkf]

It's related to what I commented in the discrete DS issue #21 (comment) but I prefer to make DS types small and immutable structs; almost like a thin wrapper around DiffEqBase.DiscreteProblem and DiffEqBase.ODEProblem for packing up phase and tangent space dynamics (and other thing? like the definition of the domain?).

For example, including integrator into the DS type is bad if you want to send the serialized DS struct over TCP (say) for doing parameter sweep of slightly different parameters and/or initial conditions using multiple machines. You'd just want to send DS definitions, not temporal variables.

[Datseris]

You can't change the system state if it is immutable, which means you have to create a new one every time you want to change state and you can't evolve in place either .

…

On Feb 8, 2018 5:21 AM, "Takafumi Arakaki" ***@***.***> wrote: It's related to what I commented in the discrete DS issue #21 (comment) <#21 (comment)> but I prefer to make DS types small and immutable structs; almost like a thin wrapper around DiffEqBase.DiscreteProblem and DiffEqBase.ODEProblem for packing up phase and tangent space dynamics (and other thing? like the definition of the domain?). For example, including integrator into the DS type is bad if you want to send the serialized DS struct over TCP (say) for doing parameter sweep of slightly different parameters and/or initial conditions using multiple machines. You'd just want to send DS definitions, not temporal variables. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#20 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ASwgYbYknsEZULl9HTcNRn242_Wj7X_Eks5tSnYtgaJpZM4R4nVL> .

[tkf]

@Datseris Please read #21 (comment) first

[Datseris]

Yeah Makes 0 sense to have an "integrator" and "solved" or discrete systems. Scientifically speaking , they are solved from the moment they are defined.

…

On Feb 8, 2018 5:57 AM, "Takafumi Arakaki" ***@***.***> wrote: @Datseris <https://github.com/datseris> Please read #21 (comment) <#21 (comment)> first — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#20 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ASwgYSuxyy6di4JfGu2jPZ5M1Q3d7NTrks5tSn6xgaJpZM4R4nVL> .

[Datseris]

Edit: "solver"

…

On Feb 8, 2018 6:07 AM, "George Datseris" ***@***.***> wrote: Yeah Makes 0 sense to have an "integrator" and "solved" or discrete systems. Scientifically speaking , they are solved from the moment they are defined. On Feb 8, 2018 5:57 AM, "Takafumi Arakaki" ***@***.***> wrote: > @Datseris <https://github.com/datseris> Please read #21 (comment) > <#21 (comment)> > first > > — > You are receiving this because you were mentioned. > Reply to this email directly, view it on GitHub > <#20 (comment)>, > or mute the thread > <https://github.com/notifications/unsubscribe-auth/ASwgYSuxyy6di4JfGu2jPZ5M1Q3d7NTrks5tSn6xgaJpZM4R4nVL> > . >

[tkf]

Another issue not mentioned in the plan:

Currently, I think t = 0 is always used in some places like evaluating Jacobian using ForwardDiff. This restriction can be removed.

Also, what ODEProblem.tspan means for the DS type can be made more meaningful by adding a "contract" in the API that it is used for inter-orthonormalization-interval in Lyapunov analysis (or something similar). In other words, users can specify the tspan that is comparable/shorter than putative "Lyapunov times" (or just some kind of timescale of the system) to avoid too many QR decompositions, etc. This also is nice for downstream libraries implementing something like attractor visualizations, orbit diagram, etc., since you can know some kind of timescale of the system from the DS type. You can make the default to such functions nicer by using "100x tspan" or something. Having default tspan = (0.0, 100.0) is rather meaningless since you don't know what the time scale of the system is. But "100x the characteristic timescale" sounds a bit more reasonable.

There is an additional nice side-effect because it works nice with periodic drive. User can use an integer multiple of periods for tspan[2] - tspan[1] then if you do

reinit!(tangent_system_integrator, u0)
solve!(tangent_system_integrator)

to solve the system for a inter-orthonormalization-interval, then the ODEProblem.f gets the t values without loss of precision with very long iterations. I'm not sure if it really matters, but some chaotic system has super long transient/convergence so probably some researchers care.

cf. This is how I do van der Pol with periodic drive:
https://github.com/tkf/LyapunovExponents.jl/blob/master/src/examples/van_der_pol.jl

Actually, the code above using reinit! is not ideal for ODE since (I suppose) the ODE integrator cannot reuse some cached information as it cannot identify tspan[1] with tspan[2]. But I wonder if it can be treated in DifferentialEquations.jl side. (I'll open an issue there.)

[Datseris]

Also, what ODEProblem.tspan means for the DS type can be made more meaningful by adding a "contract" in the API that it is used for inter-orthonormalization-interval in Lyapunov analysis (or something similar). In other words, users can specify the tspan that is comparable/shorter than putative "Lyapunov times" (or just some kind of timescale of the system) to avoid too many QR decompositions, etc

Yeah, but we can also do this as an argument to the top level function though. And it will be less restricting this way imo. The plan is to always use step! internally, and progress your time accordingly up to a user-given dt.

(EDIT: with this approach you only reinit! at the end of the function)

[tkf]

Just a quick thought. Maybe we can merge ODEs and maps:

struct DynamicalSystem{IIP, IAD, DEP<:DEProblem, JAC}
    prob::DEP
    jacobian::JAC
end

This way, supporting random dynamical system (JuliaDynamics/ChaosTools.jl#32) would be transparent by accepting prob to be SDEProblem or RODEProblem.

[tkf]

If the equations for the system are out of place, the Jacobian must be as well. The same for in-place.

If you make this guarantee (which makes sense), you don't need IIP for the DS type since it's already can be retrieved from the prob. Making the type trait for the DS type makes sense even in that case though:

import DiffEqBase: isinplace
isinplace(ds::ContinuousDynamicalSystem) = isinplace(ds.prob)

[tkf]

Define two functions: makeintegrator(CDS) and maketangentintegrator(CDS) (names to be decided), that return integrators about the system or about the system+tangent space.

Re: name, maybe just use DiffEqBase.init?

DiffEqBase.init(ds) = init(ds, Val{:phase})    # default to phase
DiffEqBase.init(ds, ::Type{Val{:phase}})       # phase integrator
DiffEqBase.init(ds, ::Type{Val{:tangent}})     # phase-tangent integrator

...or :variational or :augmented instead of :tangent.

@RainerEngelken Do you know what the standard term (or any other terms, for brainstorming) for the "derived" dynamical system? By that, I mean the ODE/map coupling the original system and the linearized system for tangent vector evolution.

I use "tangent" too but I think it's better to reserve that for the ODE/map for evolving only for tangent vectors. For example, you might want to solve the phase space dynamics first and solve tangent dynamics several times given the solution of that trajectory (by using sol(t) to get the phase space state). One example of the case you'd want to solve the phase-space dynamics separately is the CLV algorithm by Kuptsov & Parlitz (2012). I'm not arguing that the "pure-tangent integrator" has to be added in the code right now or anything. It's just that it may get into our way in the future if we don't avoid it right now.

[Datseris]

If you make this guarantee (which makes sense), you don't need IIP for the DS type since it's already can be retrieved from the prob.

True, but it is still useful for dispatch purposes. From internal tests I found it useful, and it doesn't affect any high or low level API so I will still use it.

Just a quick thought. Maybe we can merge ODEs and maps:

struct DynamicalSystem{IIP, IAD, DEP<:DEProblem, JAC}
    prob::DEP
    jacobian::JAC
end

Wonderful idea. I only hope it will work out in the end. I will try it.

Re: name, maybe just use DiffEqBase.init?

DiffEqBase.init(ds) = init(ds, Val{:phase}) # default to phase
DiffEqBase.init(ds, ::Type{Val{:phase}}) # phase integrator
DiffEqBase.init(ds, ::Type{Val{:tangent}}) # phase-tangent integrator

I am sorry, but I cannot agree with this. I heavily dislike using Val for dispatch in APIs.
We can reconsider the names, but I don't see me using Val.

[tkf]

FYI, you can do

DiffEqBase.init(ds, kind::Symbol) = init(ds, Val{kind})

to increase readability; i.e., users can do init(ds, :tangent) when they don't need performance for initialization (which is like almost always). The users who has to make the integrater in a tight loop can still use the Val version. Though I'm not super opinionated about this.

[Datseris]

At the moment I am using:

integrator(ds) -> returns appropriate integrator
tangent_integrator(ds) -> returns appropriate tangent integrator

I don't want add methods to e.g. ODEIntegrator as this would limit extendability to e.g. JuliaDynamics/ChaosTools.jl#32

What do you think about this syntax? It works and is readable for both kind of continuous or discrete systems.

(By the way I have succeeded in merging both discrete and continuous into one DynamicalSystem!)

[tkf]

I agree that it should not be ODEIntegrator, since it can be discrete DS or stochastic. There has to be a unified API to turn a DS type into a integrator.

The direction of integrator and tangent_integrator sounds good but: (1) "integrator" for discrete DS sounds wired, though I can live with it (you just have to think about appropriate measure space in your mind!). BTW, this is why I mentioned DiffEqBase.init; its semantics is already established and the name is agnostic about continuous or discrete. (2) I think the use of "tangent" has to be re-considered: #20 (comment)

[Datseris]

I am using integrator even for discrete, because in DiffEq, discrete problems are evolved using ODEIntegrator... :D
init may be reconsidered.

Alright, namings can be decided at the end, I don't really care about them atm.

btw I will need your help when I am done with everything, to go over the PR. It is a massive upgrade and one pair of eyes alone (mine) can't be trusted.

[tkf]

I am using integrator even for discrete

I know, but the libraries downstream to DynamicalSystemsBase does not have to see the name directly. Also, SDE/RODE/DDE/PDE can be used in the future. As you already noticed, avoid using ODEIntegrator helps.

About this:

Define a function step!(integrator, Δt) that steps until a span of Δt has passed.

Checkout SciML/OrdinaryDiffEq.jl#253
It may (partially) go into DifferentialEquations.jl

[tkf]

If you make this guarantee (which makes sense), you don't need IIP for the DS type since it's already can be retrieved from the prob.

True, but it is still useful for dispatch purposes.

How about using "type alias"?:

using DiffEqBase

struct DynamicalSystem{DEP<:DEProblem, JAC}
    prob::DEP
    jacobian::JAC

    DynamicalSystem(prob::DEP, jacobian::JAC = nothing) where {DEP, JAC} =
        new{DEP, JAC}(prob, jacobian)
end

const IPContinuousDS = DynamicalSystem{DEP} where {
    uType, tType, DEP <: ODEProblem{uType, tType, true}}

const NIPContinuousDS = DynamicalSystem{DEP} where {
    uType, tType, DEP <: ODEProblem{uType, tType, false}}

ds_ip = DynamicalSystem(ODEProblem{true}((a...) -> nothing, [0.0], (0.0, 1.0)))
ds_nip = DynamicalSystem(ODEProblem{false}((a...) -> nothing, [0.0], (0.0, 1.0)))

ds_ip :: IPContinuousDS
ds_nip :: NIPContinuousDS

f(::IPContinuousDS) = true
f(::NIPContinuousDS) = false
@show f(ds_ip)
@show f(ds_nip)

Since DiscreteProblem <: AbstractODEProblem (lol), swapping ODEProblem with AbstractODEProblem defines IP/NIP version of generic DS. But probably we shouldn't be using this inheritance since it won't work anyway with SDE/RODE/DDE/PDE. It would be nice if DEProblem{isinplace} was defined, but probably there are good reasons.

[Datseris]

Let's continue this on #21