Change ContinuousDS to use (t, u, du)

[Datseris]

Ooookay.

I succumb to the peer pressure. Time to change this...

[Datseris]

Also, this would make it so that the datasets from trajectory do not return the time vector.

on the other hand, the poincaresos from ChaosTools will also not return the time vector. but this is fiiiiine

[Datseris]

This will not allow for Jacobian to be computed automatically.

[Datseris]

Okay, we have to do this.

Changing the current code for the equations of motion to become (t, u, du) is trivial. I already do it anyways to create the ODESolver.

The problem is the Jacobian. Also, I want to see whether there is some way to take advantage of all the cool stuff of differentiation that is set up by Differential Equations.

What we need is a "variational integrator". An integrator that evolves forward in time both an orbit as well as deviation vectors. The equation for the deviation vectors is \dot{w} = J(u) w where of course J(u) must be computed always at the current state that is evolved in parallel with w. What we do now is the following:

function variational_integrator(ds::ContinuousDS, k::Int, t_final::Real,
    S::AbstractMatrix; diff_eq_kwargs = Dict())

    f! = ds.eom!
    jac! = ds.jacob!
    J = ds.J
    # the equations of motion `veom!` evolve the system and
    # k deviation vectors. Notice that the k deviation vectors
    # can also be considered a D×k matrix (which is the case
    # at `lyapunovs` function).
    # The e.o.m. for the system is f!(t, u , du).
    # The e.o.m. for the deviation vectors (tangent dynamics) are simply:
    # dY/dt = J(u) ⋅ Y
    # with J the Jacobian of the vector field at the current state
    # and Y being each of the k deviation vectors
    veom! = (t, u, du) -> begin
        us = view(u, :, 1)
        f!(view(du, :, 1), us)
        jac!(J, us)
        A_mul_B!(view(du, :, 2:k+1), J, view(u, :, 2:k+1))
    end

    varprob = ODEProblem(veom!, S, (zero(t_final), t_final))
    solver, newkw = get_solver(diff_eq_kwargs)
    vintegrator = init(varprob, solver; newkw..., save_everystep=false)
    return vintegrator
end

(currently ds.eom! is of the form (du, u) but I want to change that to (t, u, du). The jacobian is of the form jacob!(J, u) which writes in-place. The function of the jacobian can be provided by the user, or be auto-differentiated, using:

J = zeros(T, D, D)
    jcf = ForwardDiff.JacobianConfig(eom!, du, state)
    ForwardDiff_jacob!(J, u) = ForwardDiff.jacobian!(
    J, eom!, du, u, jcf)

The problem is that this syntax works only when the equations of motion are 2-argument f!(du, u).

But in general, maybe there is some way to take advantage of all this huge code and the wrappers that diff eq has, and make such a variational integrator a reality with the syntax (t, u, du) and enabling auto-differentiation?

@ChrisRackauckas

[ChrisRackauckas]

You can just autodiff the closure (du,u)->f(t,u,du).

BTW, a variational integrator is different so don't use that term. It's like a symplectic method, except written directly on the Legrangian. I hope to have those as well:

https://en.wikipedia.org/wiki/Variational_integrator

[Datseris]

You can just autodiff the closure (du,u)->f(t,u,du)

But then what do you do when f(t, u, du) actually depends on time?

BTW, a variational integrator is different so don't use that term.

Damnit, here I go again on name hunts!

[ChrisRackauckas]

That's what I mean. You have to differentiate it at each time using the enclosed value of the current t.

[Datseris]

Ok I am not sure I understand. So you mean to do something like this:

# and Y being each of the k deviation vectors
veom! = (t, u, du) -> begin
    us = view(u, :, 1)
    f!(t, usm view(du, :, 1))
    # Here create a new closure
    c = (du,u)->f!(t,u,du)
    # Here perform a new autodifferentiation...?
    jcf = ForwardDiff.JacobianConfig(c, du, us)
    ForwardDiff.jacobian!(J, c, du, us, jcf)

    # Continue
    A_mul_B!(view(du, :, 2:k+1), J, view(u, :, 2:k+1))
end

Wouldn't that be impossibly slow?

Also, there isn't any way to do what I am trying to do, using the existing auto-differentiation handling of DiffEq, to re-use code?

[Datseris]

https://github.com/JuliaDiffEq/DiffEqDiffTools.jl/blob/master/src/function_wrappers.jl

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/derivative_wrappers.jl

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/rosenbrock_integrators.jl#L25-L59

I am posting the links for reference, so I don't lose them.

[ChrisRackauckas]

Wouldn't that be impossibly slow?

Well only because of how you did the JacConfig. jcf = ForwardDiff.JacobianConfig(c, du, us) should be defined outside of the eom and pre-cached since that's the cache allocation for ForwardDiff and you'll want to re-use that. Making a new closure each time like that should be fine since I believe Julia should elide the allocations for that, but if it doesn't then you just make a mutable callable type for that of course.

So the only extra allocation is to the JacobianConfig. You could in theory re-use the one from OrdinaryDiffEq.jl if you are using a stiff solver. However, most algorithms (like Tsit5 and Vern) don't actually need Jacobians, so this config doesn't exist. So in theory it could be more performant to do a check for specific methods and then grab that out of integrator.cache.jac_config, but in practice I don't think it's a big deal if your system isn't large and it would be a PITA to setup which algorithms have it and which don't (though we could set it up with a trait system, file an issue if you want it).

[Datseris]

I'll open a new branch today and test it. Really hope I can make it work... We'll see. Worst case scenario I just don't do the change from (du, u) to (t, u, du).

[Datseris]

I'll keep this open in case somebody else wants to solve it / do it in the future, but I can't.

I realized that I have no idea conceptually what a Jacobian even is in the case where one wants to treat time independently or "specially". In the case of e.g. the duffing oscillator:

du[1] = u[2]
du[2] = cos(ωt) - u[1] - u[1]^3 - u[2]

what I do with the knowledge I have is simply introduce a third dimension with du[3] = 1, and make t be u[3]. The Jacobian of this system has to be a 3x3 matrix.

I have no clue what needs to happen if one wants to treat the time "differently". The problem is that if you consider the Jacobian of only the 2D system where time is something special, then

J = [0   1;
      -1,  -1 -3 [u2]^2]

i.e. it does't have time dependence at all. What do you do with time then? Obviously the dynamics are not the same for different t values, yet this is not reflected in the 2D Jacobian...

[Datseris]

Alright, after some discussion with one of the best scientists on the field (U. Parlitz) I know how to tackle this. As long as the time dependence is "normal" (i.e. linearly coupled to the system), not caring about it for the jacobian is super fine. It is also easy to show this analytically, I don't know why I didn't really think about it...

@ChrisRackauckas I will adopt the (t, u, du) approach now. There is one thing I'd like to know.

Given that I have the equations of motion, DiffEq will anyways do a automatic differentiation for the Jacobian, no matter what?

Then

1. Given that I do have a function for the Jacobian, can I just pass it to the ODEProblem constructor?
2. If I do not have the above, is there a way I can get the function from an ODEProblem structure?

[Datseris]

If what I said was not clear, I am talking about the Jacobian of the equations of motion as if they were without time,

so for this:

du[1] = u[2]
du[2] = cos(ωt) - u[1] - u[1]^3 - u[2]

it would be this the same as the jacobian of this

du[1] = u[2]
du[2] = - u[1] - u[1]^3 - u[2]

[ChrisRackauckas]

Given that I have the equations of motion, DiffEq will anyways do a automatic differentiation for the Jacobian, no matter what?

No. Only solvers for stiff equations like Rosenbrock23 or TRBDF2 actually need Jacobians. All non-stiff solvers like Vern7 don't need Jacobians and don't have any of that setup.

Given that I do have a function for the Jacobian, can I just pass it to the ODEProblem constructor?

Yes.

http://docs.juliadiffeq.org/latest/features/performance_overloads.html#Declaring-Explicit-Jacobians-1

If I do not have the above, is there a way I can get the function from an ODEProblem structure?

No, err, not really. It's not in the ODEProblem and it's only setup with specific integration algorithms. There's a way to grab those components but I'm not sure that's helpful to you.