Add the JuliaDiffEq Common Interface

[ChrisRackauckas]

Hey,
I was wondering if you would like to include the code to plug into the JuliaDiffEq "common interface". This common interface is a type-based interface which makes all of the differential equations solvers have a common solver signature. The full discussion is here:

SciML/Roadmap#5

Since the discussion is long, let me give a quick summary. There is a problem type which defines the Cauchy problem to solve:

type ODEProblem{uType,tType,isinplace,F} <: AbstractODEProblem{uType,tType,isinplace,F}
  f::F
  u0::uType
  tspan::Vector{tType}
end

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/common_interface/src/problems.jl#L26

Each algorithm is a subtype of AbstractODEAlgorithm and thus the command solve(prob,alg;kwargs....) dispatches to the appropriate solvers from the appropriate packages. A similar setup exists for DAE problems (as exemplified by Sundials.jl and now DASKR.jl).

There are many advantages to adding this interface. For one, people already know and use this API, meaning that it will be easy for users to pick up your integrators as use them as a backend simply by changing the algorithm choice. Secondly, packages can write generic algorithms on the common interface, which then allows any backend solver to be used. For example, DiffEqParamEstim.jl performs parameter estimation using an ODE solver, and any ODE solver which implements the common interface can be used due to dispatch magic (you just pass in the algorithm type you wish to use). Another example is DiffEqSensitivity.jl which uses whichever backend to perform sensitivity analysis. Other tools include DiffEqUncertainty.jl which will be able to perform uncertainty quantification on numerical methods (once common callback interfaces are sorted out), and DiffEqDevTools which provides methods to test / benchmark any solver which implements the common interface (which is used to make these benchmarks).

You can see the full set of algorithms we already have collected here (this is missing DASKR.jl and LSODE.jl which have recently implemented this same API as well):

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/solvers/ode_solve.html

It would be really nice to have geometric integrators! Let me know if you need starter code.

[michakraus]

This sounds very interesting. I would be very happy to provide that interface. However, there are a few things I need to understand first.

I made a few design decisions which result in a somewhat different structure of GeomDAE compared to DifferentialEquations. The main reason for this is performance. GeomDAE is essentially allocation-free, which can lead to about an order of magnitude better performance for fixed time step integrators compared to ODE and DifferentialEquations [*].

For example, the time span is not part of the ODE problem but of the solution. The reason for this is, that in general I want to run simulations for several hundred millions of time steps, but I do not want to keep all of those in memory. I want to run for, say, a million time steps, then save, then run another million time steps, etc... now this particular issue is easily fixed, but there are others which might be more complicated.

This is not to say that these things are cast in stone, I just have to understand what needs to be changed, if it makes sense to change it, and what would be the best way to do so. It also seems that these differences are not too severe, so this might be rather easily doable.

One thing I am not yet sure about is how your solution data structures look like. I think they are a bit different from mine. Of course, one can always copy data, but that might not be desirable. Although, I am also not yet entirely satisfied with my solution data structures. So I will have a look at yours and see if they could be unified and what that would take.

[*] Admittedly, I benchmarked this only for a few relatively small problems and only for RK4. For bigger problems and different integrators the outcome is probably different. So the above statement is not meant to be a general one. A typical outcome (here for the pendulum) looks like this:

[ChrisRackauckas]

Few comments:

GeomDAE is essentially allocation-free, which can lead to about an order of magnitude better performance for fixed time step integrators compared to ODE and DifferentialEquations [*].

DifferentialEquations.jl is allocation-free if you aren't saving the extra dense output stuff. So I am skeptical that you could get an order of magnitude faster than it for RK4, though from my benchmarks you will see that I did get that much faster than ODE.jl. Was this with save_timeseries and dense toggled? Could I see a copy of those benchmarks?

For example, the time span is not part of the ODE problem but of the solution. The reason for this is, that in general I want to run simulations for several hundred millions of time steps, but I do not want to keep all of those in memory. I want to run for, say, a million time steps, then save, then run another million time steps, etc... now this particular issue is easily fixed, but there are others which might be more complicated.

The save_timeseries option defaults to true, but if set to false the methods will not save every timestep. To save specific times, saveat will save by interpolation (i.e. the adaptive timestepper will not be constrained and it will use interpolations to save at the specified values) and tstops will constrain the solver to hit certain points. This is explained in the options common to all solvers:

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/basics/common_solver_opts.html

For example, the time span is not part of the ODE problem but of the solution. The reason for this is, that in general I want to run simulations for several hundred millions of time steps, but I do not want to keep all of those in memory. I want to run for, say, a million time steps, then save, then run another million time steps, etc... now this particular issue is easily fixed, but there are others which might be more complicated.

This can be done with infinite timespans and callbacks. Admittedly we're still working on callbacks (just OrdinaryDiffEq.jl implements them), but this simple to do.

One thing I am not yet sure about is how your solution data structures look like. I think they are a bit different from mine. Of course, one can always copy data, but that might not be desirable. Although, I am also not yet entirely satisfied with my solution data structures. So I will have a look at yours and see if they could be unified and what that would take.

The formal specification is here:

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/basics/solution.html

You can subtype AbstractODESolution and provide that interface, or simply save to the solution type defined in DiffEqBase:

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/solutions/ode_solutions.jl#L3

That solution type uses Vector{Vector{T}}. I chose this because it ends up being faster in many cases and allows for ODEs to change sizes (something that comes up in my research). That said, if your subtype of AbstractODESolution implements the same solution handling interface, then the things like plot recipes, parameter estimation, uncertainty quantification, Monte Carlo methods, etc. will all work seamlessly, so it's not necessary to use the same solution type.

Of course, one can always copy data, but that might not be desirable.

Even if you save in some other form (matrix?), that shouldn't be necessary. You could just create a Vector{Vector{T}} via views.

[ChrisRackauckas]

[*] Admittedly, I benchmarked this only for a few relatively small problems and only for RK4. For bigger problems and different integrators the outcome is probably different. So the above statement is not meant to be a general one. A typical outcome (here for the pendulum) looks like this:

Can I see the code you used? What's saved in each case?

[michakraus]

Can I see the code you used? What's saved in each case?

Sure. I am cleaning up the notebook right now. But what I am doing is essentially this:

function f_de(t, x)
    [x[2],  sin(x[1])]
end

function run_pendulum_de(Δt=Δt, ntime=ntime)
    solve(ODEProblem(f_de, x0, (0., Δt*ntime)), RK4(), dt=Δt, dense=false)
end

In both cases, GeomDAE and DifferentialEquations, I save every time step. More on this in the reply to your previous (longer) post. Stay tuned.

[michakraus]

Let me comment only on these performance issues. The other issues need some more time. I am a bit busy this and next week and do not really have time to get into this interface business right now. But I would like to work on this.

GeomDAE is essentially allocation-free, which can lead to about an order of magnitude better performance for fixed time step integrators compared to ODE and DifferentialEquations [*].

DifferentialEquations.jl is allocation-free if you aren't saving the extra dense output stuff. So I am skeptical that you could get an order of magnitude faster than it for RK4, though from my benchmarks you will see that I did get that much faster than ODE.jl.

Ok, very good. I was surprised that DifferentialEquations is so much slower and supposed that I might not run it the right way, but so far didn't really have an incentive to track this down. On the other hand, my premise was to run it the (more or less) simplest possible way, i.e., prescribing an equation, initial conditions, time step and an integration method, without getting too deep into technicalities, e.g., I assume that when I choose RK4, time steps are automatically fixed.

I added the full Jupyter notebook to the repository (examples/Performance_Comparison_Pendulum_1D.ipynb). With respect to the previous version, I added adaptive=false, which interestingly seems to yield some gain (a factor ~2 or so).

Was this with save_timeseries and dense toggled? Could I see a copy of those benchmarks?

I had dense=false but did not set save_timeseries to false as I did want to save every time step.

[mauro3]

It would be good to get feedback on what doesn't work for GeomDAE.jl in DiffEqBase.jl, so we can make adjustments.

[ChrisRackauckas]

Thanks for sharing the code. I tested it out right now and see what happened.

First of all, you didn't use the pre-allocated function with DifferentialEquations.jl. In your example, you could've just passed DifferentialEquations.jl

function f_de(t, x, fx)
    fx[1] = x[2]
    fx[2] = sin(x[1])
end

as well. Instead you passed:

function f_de(t, x)
    [x[2],  sin(x[1])]
end

That's a large part of the performance problem, and I'm not sure why you didn't use a pre-allocated function there.

Also, I put a patch in for OrdinaryDiffEq to ignore callbacks functions when they aren't doing anything because for some reason they weren't inlining properly:

JuliaLang/METADATA.jl#7181

For some reason the tags haven't gone through (that's unusually long!). So these timings are on master for OrdinaryDiffEq.jl until that goes through. It's probably about 10x slower on the current release due to this bug. Then there's also a patch on RecursiveArrayTools.jl dealing with copy vs deepcopy instances (it's complicated, basically, for simple cases like Array solutions should copy, but there was a spot where it should deepcopy).

So for a fair timing, you should checkout master on those two packages (or wait for tags) and use the pre-allocated function.

Fair Benchmark

I put together a fair benchmark where the two just solve the same function instead of different functions. This is taken from your examples folder.

What I timed was the following: (master on OrdinaryDiffEq.jl and RecursiveArrayTools.jl due to the changes mentioned above)

#Setup
using GeomDAE
using PyPlot
using DifferentialEquations

const Δt    = 0.1
const ntime = 100000
const neps  = 1E-14
const nmax  = 20

function f(t, x, fx)
    fx[1] = x[2]
    fx[2] = sin(x[1])
    nothing
end

x0 = [acos(0.4), 0.0]

# GeomDAE
ode = ODE(f, x0)
tableau = getTableauERK4()
int = Integrator(ode, tableau, Δt)
sol = Solution(ode, Δt, ntime)
integrate!(int, sol)
#set_initial_conditions!(sol, ode)
@time integrate!(int, sol)

# DifferentialEquations
prob = ODEProblem(f,x0,(0.0,ntime*Δt))
@time sol2 = solve(prob,RK4(),dt=Δt,save_timeseries=true,dense=false)

and got the timings:

# GeomDAE.jl
0.051890 seconds (4 allocations: 160 bytes)
0.052013 seconds (4 allocations: 160 bytes)
0.051918 seconds (4 allocations: 160 bytes)

# OrdinaryDiffEq.jl (DifferentialEquations.jl's main explicit methods)
0.041300 seconds (100.20 k allocations: 13.166 MB)
0.040307 seconds (100.20 k allocations: 13.166 MB)
0.040174 seconds (100.20 k allocations: 13.166 MB)

A few things to note:

First of all, the common interface has in the works "initialized versions" (still some details to work out, which is why it's not implemented yet). For example:

integrator = init(prob,RK4(),dt=Δt,save_timeseries=true,dense=false)
sol = solve(integrator)

Initialization actually doesn't cost all that much though (it's what makes the allocations look like so much more, but the one time cost actually doesn't take much time), so you don't really gain much. push! works just fine because it's ammortized.

(I know @mauro3 will be curious about this. Here the initialization cost was ~0.001 seconds, so about 1%)

This time difference is pretty much what I expected though, because this is all due to the costs of using a tableau as a vector instead of specifically written functions on constants (takes more time to write, but the results are faster since it can optimize a lot more and it doesn't have to access arrays which is a costly step). It's not as big as the usual cost because you get some gains somewhere, likely by treating dt as a constant whereas OrdinaryDiffEq.jl lets it change each step.

Interesting though.

[ChrisRackauckas]

But benchmarks aside (I knew there had to be some issue with the benchmark, that's why I got curious), after playing with it, it seems you can easily implement the common interface like:

function solve(prob,alg;kwargs...)
  ode = ODE(prob.f, prob.u0)
  int = Integrator(ode, alg.tableau, dt)
  sol = Solution(ode, dt, tspan[2]÷dt)
  integrate!(int, sol)
  ... # Solution handling is a little trickier since the solution types are different
end

which is then called like:

alg = GeomDAE(tableau = getTableauERK4())
sol = solve(prob,alg;dt=...)

The only somewhat issue is the ntimes argument handling. GeomDAE.jl doesn't do adaptive timestepping it seems? Why not just use an end time? You can always calculate an end time, and since the problem type is mutable you can always iteratively solve by doing:

sol = solve(prob,alg;kwargs...) # solve once
prob.tspan = (prob.tspan(2),prob.tspan(2)+dt*ntimes)
sol = solve(prob,alg;kwargs...) # solve on the next interval...

so using times isn't limiting. However, using number of iterations is limiting when considering adaptive timestepping.

[michakraus]

First of all, you didn't use the pre-allocated function with DifferentialEquations.jl. [...] That's a large part of the performance problem, and I'm not sure why you didn't use a pre-allocated function there.

I am not sure either. I setup this comparison probably like half a year ago and I think I just took the example from the DifferentialEquations docs/tutorial. Fixing this helps a bit, but not yet too much due to the bug you mentioned.

This time difference is pretty much what I expected though, because this is all due to the costs of using a tableau as a vector instead of specifically written functions on constants (takes more time to write, but the results are faster since it can optimize a lot more and it doesn't have to access arrays which is a costly step). It's not as big as the usual cost because you get some gains somewhere, likely by treating dt as a constant whereas OrdinaryDiffEq.jl lets it change each step.

This is very interesting. I will look into this. Maybe I can borrow your solution.

[ChrisRackauckas]

I am not sure either. I setup this comparison probably like half a year ago and I think I just took the example from the DifferentialEquations docs/tutorial. Fixing this helps a bit, but not yet too much due to the bug you mentioned.

Ahh yes, half a year ago, things were very different (I think that would've been about the time I released, right? I don't think I had ODE solvers at that point).

[michakraus]

Right now, GeomDAE.jl doesn't support adaptive time stepping. However, I have been thinking about how to implement this. In the adaptive case, I would not use ntimes to determine the number of time steps but the number of points to save the solution w.r.t. the initial time step. For the non-adaptive case I also have nsave to determine how often to save. That probably needs to be unified in a better way when adding the adaptive case, so I can also change it altogether. The biggest question here is how to prescribe when to save the solution for output. I will have another look at how you are doing this.

An interesting point in that respect is that you mentioned that you use interpolation for getting the output values at the prescribed times. For me, this is not an option as that would inherently not preserve structure (symplecticity, momenta, etc.). Instead, I was thinking about always using the actual integrator to hit the required times of the output array, i.e., if t < timeseries[i] and t+dt > timeseries[i], then set dt = timeseries[i] - t, where t is the current time and timeseries is the array of output times.

Another thought about this is the following: While in principle it wouldn't be a problem to add the final time to the ODE problem, I was thinking of this in terms of an initial value problem, which is to say that it doesn't have a final time, but is determined by an ODE dx/dt = f(x) and an initial value x(t_0) = x_0. In GeomDAE, then all the business of resetting initial values and initial times, e.g.,

prob.tspan = (prob.tspan(2),prob.tspan(2)+dt*ntimes)

which you'd also have to do for the initial value x_0, is taken care of in the solution, i.e., the solution is initialised from an ODE and than the ODE is not touched anymore. To run several integration cycles on the ODE, one just copies the last time/solution value in an SolutionODE the the initial values. From a user perspective I prefer this over fiddling with the ODE.

[ChrisRackauckas]

For reference, the output control is described here:

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/basics/common_solver_opts.html#Output-Control-1

Right now, GeomDAE.jl doesn't support adaptive time stepping. However, I have been thinking about how to implement this. In the adaptive case, I would not use ntimes to determine the number of time steps but the number of points to save the solution w.r.t. the initial time step. For the non-adaptive case I also have nsave to determine how often to save. That probably needs to be unified in a better way when adding the adaptive case, so I can also change it altogether. The biggest question here is how to prescribe when to save the solution for output. I will have another look at how you are doing this.

nsave is is like timeseries_steps, i.e. "save every n steps". (File an issue if you want to discuss naming. some of the names may be unintuitive since what they do morphed over time (but now stabilized)).

An interesting point in that respect is that you mentioned that you use interpolation for getting the output values at the prescribed times. For me, this is not an option as that would inherently not preserve structure (symplecticity, momenta, etc.). Instead, I was thinking about always using the actual integrator to hit the required times of the output array, i.e., if t < timeseries[i] and t+dt > timeseries[i], then set dt = timeseries[i] - t, where t is the current time and timeseries is the array of output times.

This is the tstops option. Interpolated points are through saveat. These are different because saveat is usually faster than taking more steps (interpolations are usually much cheaper), but tstops is also an option because interpolations aren't always possible (like in your case. My other example is solving stochastic differential equations: interpolations don't make too much sense in that case either).

While in principle it wouldn't be a problem to add the final time to the ODE problem, I was thinking of this in terms of an initial value problem, which is to say that it doesn't have a final time, but is determined by an ODE dx/dt = f(x) and an initial value x(t_0) = x_0.

What's a use case here? I'd like a concrete example. In the ones that I know of it's because the end time is determined by some event, but then a tspan like (0.0,inf) works, and you just use events to make the ODE exit when the event occurs. What's the use case where you want to solve an ODE indefinitely, but don't want to do it by an event? If it's after n iterations, then note that if you ever get above maxiters the ODE solvers will automatically exit at that point, so that's already part of the interface (but it will throw a warning).

the solution is initialised from an ODE and than the ODE is not touched anymore. To run several integration cycles on the ODE, one just copies the last time/solution value in an SolutionODE the the initial values. From a user perspective I prefer this over fiddling with the ODE.

This could be made easy by one more dispatch for ODEProblem:

ODEProblem(sol::AbstractODESolution,tend) = ODEProblem(sol.prob.f,sol[end],(sol.t[end],tend))

Then the solution + a new specified endpoint (can be inf again) makes a new problem starting from the end of the last solution. Is that what you're getting at?

[michakraus]

While in principle it wouldn't be a problem to add the final time to the ODE problem, I was thinking of this in terms of an initial value problem, which is to say that it doesn't have a final time, but is determined by an ODE dx/dt = f(x) and an initial value x(t_0) = x_0.

What's a use case here? I'd like a concrete example. In the ones that I know of it's because the end time is determined by some event, but then a tspan like (0.0,inf) works, and you just use events to make the ODE exit when the event occurs. What's the use case where you want to solve an ODE indefinitely, but don't want to do it by an event? If it's after n iterations, then note that if you ever get above maxiters the ODE solvers will automatically exit at that point, so that's already part of the interface (but it will throw a warning).

It's more a matter of how I see ODE. And that is as the definition of an initial value problem in the mathematical sense, which doesn't have a final time. But I am not overly dogmatic. I do not really care too much about that. I am fine adding a final time. Maybe we'll have shooting methods for solving boundary value problems at some point, and then this is actually needed.
Of course, we do not run our codes indefinitely, but we usually run for a prescribed time. That is, we do not stop e.g. when something happens in the solution.

This could be made easy by one more dispatch for ODEProblem:

ODEProblem(sol::AbstractODESolution,tend) = ODEProblem(sol.prob.f,sol[end],(sol.t[end],tend))

Then the solution + a new specified endpoint (can be inf again) makes a new problem starting from the end of the last solution. Is that what you're getting at?

I see that, but what I do not like so much, is doing the restart of the integrator by instantiating a new ODE. Right now this is dealt with in the solution and this is where I would prefer it to stay. If I have to create a new ODE every time I restart the integrator, in principle I also have to create a new solution all the time. After all, the integrator does not know that I am integrating the same equation for the same period of time again. I would like to avoid that. Therefore in GeomDAE, the runtime is an argument of the integrator and not part of the definition of the ODE.

Yet, that doesn't seem to prevent adding another interface like that of DifferentialEquations. I mean we can have an ODE with start and end time, and then the number of timesteps, output time series, etc., can be an argument of the solver call. When calling solve(), or rather integrate(), without pre-instatiated data structures, everything is created on the fly anyway, and that could be done in a way suitable for DifferentialEquations. Still, when one pre-instantiates everything by hand and then runs integrate on those data structures, the re-initialisation could still be done on the level of the solution instead of the ODE (or, actually just as well on the level of the ODE, if one wishes to do that).
So it seems like it is not necessary to change much of the current interface, but just to add a few more convenience functions.

[ChrisRackauckas]

Yet, that doesn't seem to prevent adding another interface like that of DifferentialEquations. I mean we can have an ODE with start and end time, and then the number of timesteps, output time series, etc., can be an argument of the solver call. When calling solve(), or rather integrate(), without pre-instatiated data structures, everything is created on the fly anyway, and that could be done in a way suitable for DifferentialEquations. Still, when one pre-instantiates everything by hand and then runs integrate on those data structures, the re-initialisation could still be done on the level of the solution instead of the ODE (or, actually just as well on the level of the ODE, if one wishes to do that).
So it seems like it is not necessary to change much of the current interface, but just to add a few more convenience functions.

Can you show what this would look like?

Note that:

and then the number of timesteps ... can be an argument of the solver call

already exists with maxiters. Normally you want that to warn though. So maybe a toggle to switch to another style?

output time series, etc., can be an argument of the solver call

Yeah, if you take a look at OrdinaryDiffEq's solve, you'll see this is already there but not documented:

function solve{uType,tType,isinplace,T<:OrdinaryDiffEqAlgorithm,F}(
  prob::AbstractODEProblem{uType,tType,isinplace,F},
  alg::T,timeseries=[],ts=[],ks=[]; # then kwargs

so output arrays can be passed in. However, the issue is that other than the timeseries for u, t, and possibly du, any other array like this is going to likely be package/algorithm specific (here, ks is used for the higher order interpolation information, or is the derivatives if the interpolation is the simple Hermite fallback).

I cannot think of how to make this less "package-specific", other than yes, the init infrastructure proposed by @mauro3 . Is that what you're thinking of?

[ChrisRackauckas]

Bumping this thread to make sure it doesn't get lost.

[ChrisRackauckas]

By the way, I updated the benchmarks.

pendulum_runtime.pdf

pendulum_memory_allocation.pdf

I'll submit a PR soon with an analysis of the profiles to explain the behaviors. The runtime in particular was interesting:

[ChrisRackauckas]

If you want to discuss some of the final details in chat, you can always find me in the JuliaDiffEq Gitter.

https://gitter.im/JuliaDiffEq/Lobby

This may be something that requires some back and forth. Once we get this worked out, I'd be willing to submit some PRs to get the common interface bindings working.

[michakraus]

Yes, once I got a little more time, we should do this. But right now I am a little too busy to wrap my head around this (delivering block lectures this and next week).

[ChrisRackauckas]

No worries, take your time.

[ChrisRackauckas]

Hey,
We have added a new integrator interface:

http://docs.juliadiffeq.org/latest/basics/integrator.html

Does this solve the problems you had before? If so, I'd be willing to write the JuliaDiffEq bindings for you and submit the PR, I just want to make sure it theoretically will work before taking the time to learn it all.

[Luapulu]

Hey, I'm new to the world of differential equations in Julia. I'm curious, what is the state of this currently? @michakraus

Am I right in thinking this package has retained its interface design and GeometricIntegratorsDiffEq.jl wraps this package to make it compatible with the DifferentialEquations.jl ecosystem?

Does GeometricIntegratorsDiffEq.jl wrapper work? The CI badge says failing.

Is there much interest in merging the interfaces and joining the DifferentialEquations.jl ecosystem or is it a lower priority issue or are there fundamental reasons for not wanting to merge, because of design differences or something like that?

[michakraus]

Thanks @Luapulu for bringing this up again. I never really found the time to look into this too deeply, but it seems like a good idea to spend some time here.

Your assumption is right. The interface wasn‘t changed, but GeometricIntegratorsDiffEq.jl can be used to run the integrators implemented in GeometricIntegrators.jl from the DifferentialEquations.jl ecosystem.

The wrapper used to work, but there seems to be some dependency incompatibility right now. Should be easy to fix, though. Besides, GeometricIntegratorsDiffEq.jl should be updated for GeometricIntegrators.jl v0.8.

Merging the interfaces is probably not necessary, as the wrapper works fine. There are a few points regarding interoperability that could be improved, though. We have a point on the todo list, that would allow to accept vector fields (e.g. of ODEs) defined with the DifferentialEquations.jl interface to be used in GeometricIntegrators.jl. Another useful and related feature would be a converter for the different equation structs.

Furthermore, we have a mid-level integrator interface since some time. It would be preferable to wrap this in GeometricIntegratorsDiffEq.jl instead of the high-level interface that is wrapped right now.

[ChrisRackauckas]

Tests passed, they just were run before the package it actually wanted to use was merged since it was tracking master. SciML/GeometricIntegratorsDiffEq.jl#19