ODE/DAE solver interface

[mauro3]

I think it would be great if all the ODE/DAE solvers could be called through the same interface. This would essentially mirror the approach of https://github.com/JuliaOpt/MathProgBase.jl . Would it be worth trying to hammer out this interface in this road-map issue?

Edit: By "same interface" I don't mean that ODE and DAE solvers have the same interface. But that all ODE solvers are wrapped to conform to the same interface (similarly all DAE, and maybe all IMEX, etc). So, Sundials.jl, ODE.jl, etc. would then depend on ODEProgBase (or probably DiffEqBase) and extend some solve function there, taking as input some problem specification type. Likewise all higher level Diff-Eq stuff (currently DifferentialEquations.jl and IVPTestSuite.jl) would program against this interface.

TODO

Specify:

• solve signature: solve(p::Problem, ::Type{Alg}; kwargs...)
• Problem types
• Solution types
• default options
• for callback options see Callback Interface #10

[ChrisRackauckas]

I am not sure how you'd do this or if it's possible. DAE solvers just require more information: (t,u,du,out) or whatever you want to name them. There's no reason for the fourth argument for an ODE solver, but the fourth argument is what the DAE solver is "zeroing". You can't make the DAE solvers in-place either because Sundials requires the in-place style (we wrap it inside the simple interface so that f! = (t,u,du,out) -> (out[:] = f(t,u,du,out)), but now the benefits of in-place updates wouldn't be available). Or we could have ODE functions take 4 arguments and ignore 1: I don't see that as a good idea because it would be slightly confusing.

You could make them all f(t,u,du) and have the DAE solver be a call overloaded type which write the output in-place to f.out, and use that to build the function for Sundials. However, then using that would require users to write the function as both in-place and as a call overloaded type which is really asking a lot unless you essentially have the main option being a DSL with macros. Also this would put heavy limitations on the usage: you wouldn't be able to use something like ForwardDiff because the restrictions on the type parameters for f.out (which would have to be strict for performance) would be too strict to allow DualNumbers (this is why ParameterizedFunctions.jl ends up generating extra "parameter explicit" functions to allow for explicit parameter derivatives in sensitivity analysis).

[mauro3]

Sorry, I should have been clearer, I updated above post.

[ChrisRackauckas]

There is already a unified interface in DifferentialEquations.jl. The only difference is that the problem object needs to be specified slightly differently (the function specification difference noted above, and DAEs need an initial condition for the derivative). You can compare this two pages (note the DAE docs here a little bit behind because of the modularization changes):

http://juliadiffeq.github.io/DifferentialEquations.jl/latest/tutorials/ode_example.html
http://juliadiffeq.github.io/DifferentialEquations.jl/latest/tutorials/dae_example.html

All of the other analysis and plotting functions work the same (I don't have the overloaded interpolations on the DAEs because I need to make a special output for Sundials for that though) since this is all handled (in the newest modularization change) in DiffEqBase. The same interface/API and functionality is then provided to SDEs, PDEs, etc, and will be provided to DDEs by the end of next month. I think I've shown that it's a good and extendable solution, though it's likely in need of dev docs (which I'm adding soon now that the modularization split happened).

If someone wants to take a different approach, that's fine. I don't know what exactly you would gain (or what a different approach would look like) but offshoots of experimental approaches aren't harmful.

[mauro3]

Yes, I know, the interface of DifferentialEquations.jl resides in https://github.com/JuliaDiffEq/DiffEqBase.jl. Here I'm wondering whether it would make sense to declare it, or some variation thereof as the "community standard". Packages like Sundials.jl, OrdinaryDiffEq.jl, ODEInterface.jl, ODE.jl, etc. could then implement this interface (or make a wrapper package) and this would open up their use elsewhere. Also making this a "standard" would make them more accessible from outside the DifferentialEquations.jl-universe. For me as a user this would be useful, I'd depend on DiffEqBase.jl and could call into all solvers. I know that currently I can do this with DifferentialEquations.jl, but usually I prefer hooking into a lower level package.

Then, on a personal note, for me as a contributor, this would also help to alleviate my feeling that DifferentialEquations.jl is steamrolling all ode related stuff, that any contributions of mine would be greeted by a "I'll wrap this in DifferentialEquations.jl" and just end up being an insignificantly tiny cog in its machinery. If I could instead program against such an interface, I would probably feel a bit more independent and motivated to contribute. (Whilst it would still be immediately usable to DifferentialEquations.jl) Of course, I realize that my personal feelings should not be an indication of how to take JuliaDiffEq forward, even more so as the numbers and quality of my contributions will stay modest at best. At the end the best & best-maintained software should and will be the "winner".

Last, maybe it's still too early to do this. DifferentialEquations.jl is young, so is the iterator interface in ODE.jl. Thus it may not be clear yet what a useful interface would be.

[ChrisRackauckas]

Here I'm wondering whether it would make sense to declare it, or some variation thereof as the "community standard".

The thing is, I don't know what right I, or anyone else, would have to declare a "community standard". What I can do is build everything I do off of this standard, document it, and personally recommend it. I am at least building it as a standard type hierarchy and interface for everything that I am doing. If there are ways to to make this a more useful interface in a broader sense, let's discuss it.

On that note Sundials.jl is undergoing an API change. I was trying to keep it "as Sundials-y" as possible (i.e. not introduce problem/solution types) but am realizing that to take full advantage of Sundials I am going to want at least a special solution type. To handle this I was going to write an alternative Sundials API using DiffEqBase.jl things, and use this through OrdinaryDiffEq.jl so as to not affect the current Sundials API too much.

Then, on a personal note, for me as a contributor, this would also help to alleviate my feeling that DifferentialEquations.jl is steamrolling all ode related stuff, that any contributions of mine would be greeted by a "I'll wrap this in DifferentialEquations.jl" and just end up being an insignificantly tiny cog in its machinery.

I hope that a small contribution isn't seen as insignificant. If it's the right algorithm for the right problem then it's a significant contribution for that purpose. I plan on making sure Julia has (optimized) implementations available in DifferentialEquations.jl for every algorithm that I know of or can find, so in that sense any algorithm is a small but meaningful contribution.

The thing I really need to think about is how to deal with existing implementations which I think are not optimized and may need to be re-written. For example, the GSoC project is written in a variable-coefficient format instead of a Nystrom format which is known to not be as well-optimized (and it uses a lot of mutable data structures in the inner loops). There's a double-edged sword here that I am going to be implementing a Nystrom version of ABM methods in order to have a well-optimized version available, but yes it would be in some sense stepping on another's toes (and the change would pretty much require a re-design of the current ABM anyways). I do see how this could build a culture against contributions and a "ehh just let Chris do it attitude" which is not good for the me or the community because, even though I have been pretty productive, this has a huge scope and I do need help.

Last, maybe it's still too early to do this. DifferentialEquations.jl is young, so is the iterator interface in ODE.jl. Thus it may not be clear yet what a useful interface would be.

Indeed, that's why I'm taking a kind of brute-force tactic of "build it all and see what sticks". However, maybe we can hammer out a type-structure which works for all of us. Maybe this should then move over to DiffEqBase.jl, but here's some thoughts on that:

• It may be better in many cases to define an interface over abstract types. For example, the plot recipes work on AbstractODESolution types, and AbstractSDESolution <: AbstractODESolution and contain fields t and timeseries so they plot the same. A clear definition of the interface would allow people to write tools over them, or get all of this stuff "for free".
• The problem types I use are pretty simple, see this. Maybe that could go over to DiffEqBase.jl with some modifications (add places for Jacobians and Hessians? I think that the Jacobian and Hessians are properties of the function, and should thus be interfaced as a dispatch on the function itself like in ParameterizedFunctions.jl. Because we are using such a powerful language, we should also allow for symbolically calculated inverse/pre-factorized Jacobians and Hessians which again ParameterizedFunctions.jl already does, and so I think it's another case where an interface through overloading the function makes it easier for people to use this without the extra machinery if it doesn't exist)
• The Solution objects definitely have to be tailored to the solver output. Think about how for Sundials the mem type has to be stored and all interpolations need to be built for that, which is completely different than native algorithm dense output. Another case where this crops up is that with Hermite interpolations you want to store the derivative at each timepoint, whereas for higher order interpolations you want a vector of derivatives at evaluation points for each timepoint, so both what is stored and the algorithm on it are very different. What we could do is just define an abstract interface here, where for an AbstractODESolution, dense==true means that the object has its call overloaded with the interpolating function, that is sol(t) is the continuous solution at time t. This is already pretty much standardized, and other interfaces can be documented like that. Then, even if the implementation are very different for solutions of different algorithms, and algorithm written on an AbstractODESolution (like parameter estimation or sensitivity analysis) can know to just use sol(t), no matter the algorithm.

And the real biggie:

• ode_solve.jl in OrdinaryDiffEq.jl is a mess. I am not sure using alg=:DP5 and what not for choosing the algorithm was the right idea. If this was instead based off of types and dispatch (with a smart tolerance and specifier-based default dispatch using a plan_ode method, or just a dispatch from a DefaultODEAlgorithm), then conditional dependencies wouldn't be necessary, and would could just add new algorithms to OrdinaryDiffEq.jl in their own separate package by adding a new dispatch (you can already do this, but you would currently need to add a lot of stuff to internal dictionaries). If fixing up this algorithm selection mechanism would help bring in contributors, then this will get bumped up the priority list.

All in all, the modularization changes were really to help pull in the possibility for others to contribute. DifferentialEquations.jl was undergoing some pretty large changes daily which would have made it a little bit tough to contribute, but with it being modular now, at least some of the packages during some timeframes will be more stable (I'll actually be more productive now that this is transitioning from being a side project I work on at night. Julia and JuliaDiffEq is now part of my research instead of a distraction, and so I can easily double my amount of development time towards this stuff now. So it's moreso that certain packages will be stable since I am concentrating more time elsewere, mostly with stochastic ODEs/PDEs and inverse problems. I never intended to work on ODEs, but I just needed something which was fast enough for the PDEs I was solving). I hope the stability makes it easier for others to contribute, and I hope to have a load of dev docs explaining a lot of what I've done.

[mauro3]

Here I'm wondering whether it would make sense to declare it, or some variation thereof as the "community standard".

The thing is, I don't know what right I, or anyone else, would have to declare a "community standard". What I can do is build everything I do off of this standard, document it, and personally recommend it. I am at least building it as a standard type hierarchy and interface for everything that I am doing. If there are ways to to make this a more useful interface in a broader sense, let's discuss it.

I think that is fine as a community standard. I guess one thing is that the interface ought to be stable-ish, which sounds like you're heading towards.

On that note Sundials.jl is undergoing an API change. I was trying to keep it "as Sundials-y" as possible (i.e. not introduce problem/solution types) but am realizing that to take full advantage of Sundials I am going to want at least a special solution type. To handle this I was going to write an alternative Sundials API using DiffEqBase.jl things, and use this through OrdinaryDiffEq.jl so as to not affect the current Sundials API too much.

Yes, it would be just fine for solvers to have their own interface too. Although, one indication of a good the standard interface would be that most adapt it as the only interface eventually.

The thing I really need to think about is how to deal with existing implementations which I think are not optimized and may need to be re-written. For example, the GSoC project is written in a variable-coefficient format instead of a Nystrom format which is known to not be as well-optimized (and it uses a lot of mutable data structures in the inner loops). There's a double-edged sword here that I am going to be implementing a Nystrom version of ABM methods in order to have a well-optimized version available, but yes it would be in some sense stepping on another's toes (and the change would pretty much require a re-design of the current ABM anyways). I do see how this could build a culture against contributions and a "ehh just let Chris do it attitude" which is not good for the me or the community because, even though I have been pretty productive, this has a huge scope and I do need help.

That is fine, no worries. So all the ODE solvers you're writing would go into DiffEqODE.jl (or DAE), right?

Last, maybe it's still too early to do this. DifferentialEquations.jl is young, so is the iterator interface in ODE.jl. Thus it may not be clear yet what a useful interface would be.

Indeed, that's why I'm taking a kind of brute-force tactic of "build it all and see what sticks". However, maybe we can hammer out a type-structure which works for all of us. Maybe this should then move over to DiffEqBase.jl, but here's some thoughts on that:

• It may be better in many cases to define an interface over abstract types. For example, the plot recipes work on AbstractODESolution types, and AbstractSDESolution <: AbstractODESolution and contain fields t and timeseries so they plot the same. A clear definition of the interface would allow people to write tools over them, or get all of this stuff "for free".

Yes, that sounds right. Maybe it would make sense to have a standard type which has one field which can be used freely:

type ODESolution{T,Y,O}
    t::T # times
    y::Y # solution
    ...
    other::O # put whatever else is needed for a particular solver in-here
end

• The problem types I use are pretty simple, see this. Maybe that could go over to DiffEqBase.jl with some modifications (add places for Jacobians and Hessians? I think that the Jacobian and Hessians are properties of the function, and should thus be interfaced as a dispatch on the function itself like in ParameterizedFunctions.jl. Because we are using such a powerful language, we should also allow for symbolically calculated inverse/pre-factorized Jacobians and Hessians which again ParameterizedFunctions.jl already does, and so I think it's another case where an interface through overloading the function makes it easier for people to use this without the extra machinery if it doesn't exist)

I like simple. Julia's powers make it easy to come up with too complicated solutions. Having said this, I'll need to check ParameterizedFunctions.jl again as that did seem quite complicated to me.

• The Solution objects definitely have to be tailored to the solver output. Think about how for Sundials the mem type has to be stored and all interpolations need to be built for that, which is completely different than native algorithm dense output. Another case where this crops up is that with Hermite interpolations you want to store the derivative at each timepoint, whereas for higher order interpolations you want a vector of derivatives at evaluation points for each timepoint, so both what is stored and the algorithm on it are very different. What we could do is just define an abstract interface here, where for an AbstractODESolution, dense==true means that the object has its call overloaded with the interpolating function, that is sol(t) is the continuous solution at time t. This is already pretty much standardized, and other interfaces can be documented like that. Then, even if the implementation are very different for solutions of different algorithms, and algorithm written on an AbstractODESolution (like parameter estimation or sensitivity analysis) can know to just use sol(t), no matter the algorithm.

In all types there needs to be room for custom fields. Letting t use to get the solution at time t looks good. Maybe use [] to get it though to mirror Interpolations.jl?

And the real biggie:

• ode_solve.jl in OrdinaryDiffEq.jl is a mess. I am not sure using alg=:DP5 and what not for choosing the algorithm was the right idea. If this was instead based off of types and dispatch (with a smart tolerance and specifier-based default dispatch using a plan_ode method, or just a dispatch from a DefaultODEAlgorithm), then conditional dependencies wouldn't be necessary, and would could just add new algorithms to OrdinaryDiffEq.jl in their own separate package by adding a new dispatch (you can already do this, but you would currently need to add a lot of stuff to internal dictionaries). If fixing up this algorithm selection mechanism would help bring in contributors, then this will get bumped up the priority list.

We were/are looking at this over in PR49 and we do dispatch on the solver-type because of the concern you expressed above. Our strategy is the following: a problem is a pair: (IVP-specification, integrator), which corresponds to your (ODEProblem, ODEIntegrator). In our case the types are called (IVP, AbstractSolver) and wrapped in a Problem type. (Note, a bit confusing: what you call "Problem" is our IVP-specs only, whereas we call the pair IVP+integrator the "Problem"). See here.

The thing we noted then is that the integrator instance (<:AbstractSolver) can only be made once the IVP is specified as its type-parameters depend on the parameters of IVP (say time is in Float32). Thus to create an instance of a Solver you need to pass it the IVP-instance: Solver(ivp; opts...) and thus Problem(ivp, Solver(ivp; opts...)). (And the same goes for the structure holding the options as its type parameters also depend on the IVP-parameters and solver parameters. This also goes in direction of addressing your plan_ode concern.)

Thus a integrator needs to provide the constructor Solver(ivp; opts...) and then it can be used in PR49. A simple example is here and more complicated one is here

TLDR Maybe a suitable interface would thus be:

abstract DEProblem
abstract AbstractIntegrator

solve{I<:AbstractIntegrator}(de::DEProblem, ::Type{I}; options...) = solve(de, I(de, options))

If for some reason a package cannot implement the I(de, options) constructor, then it can just overload solve directly solve(de::DEProblem, ::Type{MyIntegrator}; options...) = ...

We should provide standard subtypes of DEProblem for say ODE, DAE, etc. But if a solver requires a special subtype of DEProblem then it can define one.

Not sure whether we should try to standardize the options too. In PR49 we try this: https://github.com/JuliaODE/ODE.jl/blob/master/src/options.jl.

Also, it maybe good to have several levels to the interface. Above would be the top-level which would make most functionality of down-stream packages available. But then if something extra is needed, say sensitivity analysis requires some extra bits, then a deeper layer of the interface also needs implementing to unlock that feature.

I never intended to work on ODEs, but I just needed something which was fast enough for the PDEs I was solving.

Same here ;-)

[ChrisRackauckas]

I think that is fine as a community standard. I guess one thing is that the interface ought to be stable-ish, which sounds like you're heading towards.

Indeed, I am aiming for this to be stable and be pretty much untouched. In order to get there, it's pretty bare and mostly includes abstract types with no concrete types. I need to power through the DDEs, some FDM implementations, and usage of explicit Jacobians/Hessians, to really have the Problem types nailed down. It's easy to think it's setup alright for these things, but there's probably something that I missed.

So all the ODE solvers you're writing would go into DiffEqODE.jl (or DAE), right?

Yes, all directly into the ode_integrators.jl file. There's a top level above it, ode_solve.jl, which chooses to call the algorithms out of there, or interfaces with Sundials, ODEInterface, ODE. ode_integrators.jl has a bunch of associated tools so that writing algorithms is fast and easy (and they end up with top-notch performance), but because of that there is quite a lot of complexity in there. But given this structure, that's why I am fine with having solvers in different packages because I can just bring them all together anyways as long as the API is "close enough". However, given the DSL in ode_integrators.jl, maybe I should just dev doc that more because you can essentially copy/paste algorithms out of a paper to get it running in there, so the only reason why others wouldn't use it is because it's not documented enough.

Yes, that sounds right. Maybe it would make sense to have a standard type which has one field which can be used freely

At that point though we're doing it because we want an interface. I think it should just go the route of a documented interface like the rest of Julia.

As with the use of y, I don't think you're going to get me to adopt that bad idea. If you use y, you have to do a bunch of things:

• You have to drop the compatibility of the basic parts of the interface between most PDEs and ODEs. If you have an interface with t as time and y as the final solution, but then y doesn't make sense for the solution since it has space? There's a reason why PDE books use u.
• y doesn't even make sense for many of the test ODEs, like Lorenz and Pleides. Anything physical which has a 2D location has a y as a variable, and so you already have name clashing.
• y(t) isn't even standard amongst ODE books. x(t) is another standard, or for books which will do PDEs later, many just do u(t) from the beginning.

I like simple. Julia's powers make it easy to come up with too complicated solutions. Having said this, I'll need to check ParameterizedFunctions.jl again as that did seem quite complicated to me.

(p::LotkaVolterra)(t,u,du) = begin
         du[1] = p.a * u[1] - p.b * u[1]*u[2]
         du[2] = -3 * u[2] + u[1]*u[2]
end
function (p::LotkaVolterra)(t,u,du,J,::Type{Val{:Jac}})
  J[1,1] = p.a - p.b * u[2]
  J[1,2] = -(p.b) * u[1]
  J[2,1] = 1 * u[2]
  J[2,2] = -3 + u[1]
  nothing
end

It's just an interface on call overloaded, which is now just a standard and widely used feature in Julia. Most users shouldn't have to use it directly anyways because the macro is powerful enough (or if it's not, make it more powerful).

In all types there needs to be room for custom fields. Letting t use to get the solution at time t looks good. Maybe use [] to get it though to mirror Interpolations.jl?

No, [] clashes with the array interface. sol[i] is the ith numerical solution. This makes the solution object work in the array interface and thus makes it easy to plug into already duck-typed algorithms. sol(t) does based off of time. If you make getindex dispatch off of floats as time and ints and the solution point, then sol[2.0] != sol[2] which is bad behavior. You can't guarantee that everyone will want continuous output either because it takes more memory and for some equations/problems a continuous output is neither easy or always possible (always have SPDEs in mind, or even better, RPDDAEs).

Note, a bit confusing: what you call "Problem" is our IVP-specs only, whereas we call the pair IVP+integrator the "Problem"

That is highly confusing, because when I tell people about an ODE problem, it's always independent of numerical methods. There's a semantic issue with having the two together: the integrator is how you solved it, not the specification of the mathematical problem (it is part of the computational problem specification I guess). Optim.jl's dispatch method is better than both of ours since it lets you choose the algorithm in the solve command (which the alg=:DP5 already lets you do), but then does the choice based off of non-dynamic dispatch (which the PR49 method does). It's the same as solve(de::DEProblem, ::Type{MyIntegrator}; options...) = ....

We should provide standard subtypes of DEProblem for say ODE, DAE, etc. But if a solver requires a special subtype of DEProblem then it can define one.

Trying. DDEs/explicit Jacobian usage needs to get hammered out before those are stable.

Not sure whether we should try to standardize the options too.

Yes, but that set of options is far too restrictive. It doesn't allow control even over the true options available in ODE.jl, which is a small subset of anything OrdinaryDiffEq or ODEInterface does. This is the set of options I've needed so far: https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve/ode_constants.jl#L53 . It doesn't include norm choice because a devectorized and non-allocating form for that needs to be hammered out. Output controls, verbosity (progressbar) control, stabilization controls, and callbacks are all part of standard ODE interfaces which are not present in the ODE.jl list of options but are essential for usage. It gets complicated because I have to alias between the naming differences of the different packages (you can see all the aliases on that page).

Also, it maybe good to have several levels to the interface. Above would be the top-level which would make most functionality of down-stream packages available. But then if something extra is needed, say sensitivity analysis requires some extra bits, then a deeper layer of the interface also needs implementing to unlock that feature.

Of course. Call overloading makes this seamless on the problem specification side. Surprisingly it's the solution handling which is much more complex, but that's developer and not user details so it just needs docs.

[mauro3]

It's easy to think it's setup alright for these things, but there's probably something that I missed.

Yep, code always gets better with 20-20 hindsight.

So all the ODE solvers you're writing would go into DiffEqODE.jl (or DAE), right?

Yes, all directly into the ode_integrators.jl file. There's a top level above it, ode_solve.jl, which chooses to call the algorithms out of there, or interfaces with Sundials, ODEInterface, ODE. ode_integrators.jl has a bunch of associated tools so that writing algorithms is fast and easy (and they end up with top-notch performance), but because of that there is quite a lot of complexity in there. But given this structure, that's why I am fine with having solvers in different packages because I can just bring them all together anyways as long as the API is "close enough". However, given the DSL in ode_integrators.jl, maybe I should just dev doc that more because you can essentially copy/paste algorithms out of a paper to get it running in there, so the only reason why others wouldn't use it is because it's not documented enough.

Maybe I'm missing something but how about a slightly different organization: DiffEqBase.jl has all in it to make a callable, standard wrap/interface for a solver-package using it. So, for instance, Sundials.jl would get an additional file where that interface is implemented (as opposed to in OrdinaryDiffEq.jl). For packages outside of JuliaDiffEq we could have a wrapper package implementing the interface (or submit a PR to the package in question). OrdinaryDiffEq.jl would then be the place where you have the solvers you coded and not where all the wraps are too. Example:

A package needing an ODE integrator could then only depend on DiffEqBase:

module SuperDuperModel
using DiffEqBase

...

run_it{I<:AbstractIntegrator}(model::SuperDuperType, odesolver::Type{I}) = solve(model.eq, odesolver)

end

and then a user of that package could pick their solver of choice to use with it:

using SuperDuperModel
import Sundials
model = ...
run_it(model, Sundials.Ida)

I think this decoupling would make the whole more modular. As an added bonus it would also clear my personal discomfort mentioned above.

Yes, that sounds right. Maybe it would make sense to have a standard type which has one field which can be used freely

At that point though we're doing it because we want an interface. I think it should just go the route of a documented interface like the rest of Julia.

Sounds good.

As with the use of y, I don't think you're going to get me to adopt that bad idea. If you use y, you have to do a bunch of things:

This is just about naming, using y vs u, right? I'm easy. But preferably consistent. Looks like, you use u, t and (internally) for step size \Delta t. I'll go with that in PR49 too.

Note, a bit confusing: what you call "Problem" is our IVP-specs only, whereas we call the pair IVP+integrator the "Problem"

That is highly confusing, because when I tell people about an ODE problem, it's always independent of numerical methods. There's a semantic issue with having the two together: the integrator is how you solved it, not the specification of the mathematical problem (it is part of the computational problem specification I guess).

I'm not advocating any particular names here, my above statement was just to not get you confused. But yes, we thought of the "Problem" as all that was needed to actually produce an output, i.e. maths+numerics.

Optim.jl's dispatch method is better than both of ours since it lets you choose the algorithm in the solve command (which the alg=:DP5 already lets you do), but then does the choice based off of non-dynamic dispatch (which the PR49 method does). It's the same as solve(de::DEProblem, ::Type{MyIntegrator}; options...) = ....

I guess the difference is that Optim uses an instance and PR49 uses a type for dispatch. The reason we went for the latter is that the type parameter of the created instance could be set to match those of the IVP-specs.

Not sure whether we should try to standardize the options too.

Yes, but that set of options is far too restrictive.

Yes, it would only make sense to standardize the common options to keep them consistent between solvers, e.g. abstol, reltol... But many/most solvers will have their specialized options too.

[ChrisRackauckas]

Maybe I'm missing something but how about a slightly different organization: DiffEqBase.jl has all in it to make a callable, standard wrap/interface for a solver-package using it. So, for instance, Sundials.jl would get an additional file where that interface is implemented (as opposed to in OrdinaryDiffEq.jl). For packages outside of JuliaDiffEq we could have a wrapper package implementing the interface (or submit a PR to the package in question). OrdinaryDiffEq.jl would then be the place where you have the solvers you coded and not where all the wraps are too.

Yes, using dispatch will make this better and work almost automatically. Then it will make OrdinaryDiffEq.jl simpler, but also get rid of the conditional dependencies that it uses.

This is just about naming, using y vs u, right? I'm easy. But preferably consistent. Looks like, you use u, t and (internally) for step size \Delta t. I'll go with that in PR49 too.

Yeah, just naming. t and u. \Delta t was a bad idea which bit me when I had to use a cluster which only supported ASCII. After that experience, I want to make sure anything that users have to interface with doesn't require Unicode. The answer to that is to have dt alias with \Delta t so you can set the option however you wish. I think that's the only user-facing Unicode that I had.

I guess the difference is that Optim uses an instance and PR49 uses a type for dispatch. The reason we went for the latter is that the type parameter of the created instance could be set to match those of the IVP-specs.

It probably will be types instead of instances of immutable types.

solve(de::AbstractODEProblem, tspan, alg=DefaultODEAlgorithm; options...)

(Should tspan be moved to be in the problem type?)

Dispatches to different packages could then be made off of abstract types

function solve{alg<:OrdinaryDiffEqAlgs}(...)

with OrdinaryDiffEqAlgs<:AbstractODEAlgs. This would allow packages to have "common setup functions", like a function which sets up output to any Sundials algorithm, or the like. The setup is necessary because this outer function naturally puts a function-barrier between any type-instability. Then instead of having dictionaries of algorithms to choose solvers, this would naturally make dispatch for things like CVODE_BDF <: SundialsODEAlgs <: AbstractODEAlgs. DefaultODEAlgorithm would then have an interface before a solver in order to choose an algorithm, based on solver options.

Is this design making sense? Have any suggestions?

Yes, it would only make sense to standardize the common options to keep them consistent between solvers, e.g. abstol, reltol... But many/most solvers will have their specialized options too.

Even more difficult is when default options are different. And the biggest difficulty, the one which lead to the most code bloat for me, is algorithm-specific defaults. For example, good default stabilization parameters for the Verner methods are different from what should be used with DP5, so even kwargs don't work well because there needs to be a good way to specialize the options to the algorithm choices.

Specialized options to solvers are not an issue if you just always kwargs... since that will simply pass things on if you don't know what to do with them. Writing the solve command like that would allow different implementations to choose which options to use, but ignore others (without error). Documenting it well require a Plots.jl-like attribute table.

---

This is planning a large refactor which I won't be implementing right now (I'm refactor tired out... this modularization has been a lot of work for "no gains"... well, organizational gains, but I didn't get to see anything go faster...). Let's try to hammer this design out and I probably won't get to it until the end of the month, but will handle getting both OrdinaryDiffEq and Sundials together. ODEInterface kind of does its own thing... so what I have in OrdinaryDiffEq might need to essentially just turn into a wrapper package which adds this interface to ODEInterface. This will get rid of the conditional dependencies and make Tony happy too.

[mauro3]

solve(de::AbstractODEProblem, tspan, alg=DefaultODEAlgorithm; options...)

(Should tspan be moved to be in the problem type?)

@pwl and I discussed this at length but did not come to a conclusive conclusion. Good to see that you're struggling too ;-) I think that it should go into the AbstractODEProblem: t0 definitely belongs into AbstractODEProblem. So, I think it would be good the have the whole of tspan there.

Dispatches to different packages could then be made off of abstract types

function solve{alg<:OrdinaryDiffEqAlgs}(...)

with OrdinaryDiffEqAlgs<:AbstractODEAlgs. This would allow packages to have "common setup functions", like a function which sets up output to any Sundials algorithm, or the like. The setup is necessary because this outer function naturally puts a function-barrier between any type-instability. Then instead of having dictionaries of algorithms to choose solvers, this would naturally make dispatch for things like CVODE_BDF <: SundialsODEAlgs <: AbstractODEAlgs. DefaultODEAlgorithm would then have an interface before a solver in order to choose an algorithm, based on solver options.

Is this design making sense? Have any suggestions?

Yes this sounds right. To bikeshed names: I'd prefer AbstractIntegrator instead of AbstractODEAlg, as Alg is a bit non-specific. But maybe Integrator is not specific enough either when considering the whole DiffEq-universe.

Yes, it would only make sense to standardize the common options to keep them consistent between solvers, e.g. abstol, reltol... But many/most solvers will have their specialized options too.

Even more difficult is when default options are different. And the biggest difficulty, the one which lead to the most code bloat for me, is algorithm-specific defaults. For example, good default stabilization parameters for the Verner methods are different from what should be used with DP5, so even kwargs don't work well because there needs to be a good way to specialize the options to the algorithm choices.

How about something like this for standard options (adapted form PR49):

immutable AdaptiveOptions{T,N<:Function,O<:Function} <: Options{T}
    reltol::T
    abstol::T
    ...
end
# default constructor defines defaults, dispatches on AbstractIntegrator:
function (::Type{AdaptiveOptions{T}}){I<:AbstractIntegrator, T,N,O}(::Type{I};
                                                     reltol   = eps(T)^T(1//3)/T(10),
                                                     abstol   = eps(T)^T(1//2)/T(10),
                                                     ...,
                                                     kargs...)
    AdaptiveOptions{T,N,O}(reltol,abstol,...)
end

...

## For a solver which has different defaults overload constructor:
function (::Type{AdaptiveOptions{T}}){I<:SunidalsIntegrator, T,N,O}(::Type{I};
                                                     reltol   = 1e-4,
                                                     abstol   = 1e-6,
                                                     ...,
                                                     kargs...)
    AdaptiveOptions{T,N,O}(reltol,abstol,...)
end

An integrator type would then look like:

abstract SundialsIntegrator{T} <: AbstractIntegrator{T}
immutable SunidalsIda{T} <: SundialsIntegrator{T}
    opts::AdaptiveOptions{T}
    other_opts::IdaOpts{T} # solver specific options
    ...
end

function SunidalsIda{T}(::AbstractODEProblem{T}; opts...)
    SunidalsIda{T}( AdaptiveOptions(SunidalsIda; opts...), IdaOpts(opts...) )
end

# and finally in solve it would look
function solve{I<:SunidalsIntegrators}(de::AbstractODEProblem, alg::Type{I}; opts...)
    solver = I(de, opts...)
    sundials_solve(...)
end

I think something along these lines could work, but maybe too complicated and not general enough...

Specialized options to solvers are not an issue if you just always kwargs... since that will simply pass things on if you don't know what to do with them. Writing the solve command like that would allow different implementations to choose which options to use, but ignore others (without error). Documenting it well require a Plots.jl-like attribute table.

Yes, that is good and what I put in above example.

This is planning a large refactor which I won't be implementing right now (I'm refactor tired out... this modularization has been a lot of work for "no gains"... well, organizational gains, but I didn't get to see anything go faster...). Let's try to hammer this design out and I probably won't get to it until the end of the month, but will handle getting both OrdinaryDiffEq and Sundials together. ODEInterface kind of does its own thing... so what I have in OrdinaryDiffEq might need to essentially just turn into a wrapper package which adds this interface to ODEInterface. This will get rid of the conditional dependencies and make Tony happy too.

:-) your clock ticks at a different rate to mine. In a month's time is plenty plenty early!

[pwl]

I've just caught up with the whole discussion. All of this looks quite right to me, and seems to be intersecting PR49 at surprisingly many places.

As for the naming I would rather have Problem for IVP and Solver for the combination of IVP+algorithm, but I wouldn't mind the current PR49 naming.

The example https://github.com/JuliaODE/ODE.jl/blob/master/examples/custom_integrator.jl shows how one can add a solver to the current PR49 interface and also shows how to handle options ask kargs. The Options type that we implemented in PR49 was merely for convenience and it's use is not obligatory, it just implements the common parameters used by most algorithms. I think that implementing a number of things to do when implementing a new solver should be kept down to a minimum, as @mauro3 suggested, we could have a simple dependency (using DiffEqBase or whatnot) and then define a few methods to define a solver. In this example we are only defining init, onestep! and the constructor for EulerIntegrator{T} and the integrator is ready to go with minimal dependencies.

As for tspan=[t0,t1] being a part of IVP, we only kept it separate because of the iterator interface. In principle the iterators don't need t1 as they don't need to know the stop time. If you don't plan on adding iterators to your interface there is no need to keep tspan out of IVP. That said, I still think that the iterators are the way to go with ODE solvers as they give you more control over the interrupting, saving and plotting of the solution without having to implement these inside the integrators as callbacks or options. At the same time you can build the higher-level interface on top of them so the end user doesn't have to use or understand them, you also don't have to care much about them when implementing a new solver, as illustrated in the example of EulerIntegrator{T} I posted earlier.

I will try to follow this discussion as much as time permits, but @mauro3 seems to be doing a good job in conveying the problems and ideas that we encountered while developing PR49. Congrats on becoming a full time JuliaDiffEq developer @ChrisRackauckas! I basically only have weekends to think about this as my daily schedule is now completely covered from 8am to 8pm.

[ChrisRackauckas]

As for the naming I would rather have Problem for IVP and Solver for the combination of IVP+algorithm

That does work in any case where one wants to define problems independently of the solver. It incorrectly adds a difficulty to anyone who wants to "change the method". In something like DiffEqProblemLibrary.jl, do we just set default methods? The current work is to move the problem definitions out of OrdinaryDiffEq.jl et. al. and into DiffEqBase so that way DiffEqProblemLibrary can have no dependents, but if you have to define the problem with an algorithm it will still need the dependency on the solver algorithm packages. If the solution to this problem is to put dummy BlankSolver into each problem that is defined (or if it gets DefaultSolver, then it would still require the package which defines that!), then that's admitting there's a problem with putting the solver method with the algorithm selection. This problem pops up in all of the *Models packages too which are for creating models and domain-specific problems which will be solver independent (and not have a dependency on the solvers) so that way the tools can be easily used outside of DiffEq (MultiScaleModels.jl is already picking up other uses).

So because it's semantically weird and only causes problems for actual implementation, I don't think they should be put together. I know that you want to put them together because in the iterator interface you want to pass one object around to everything. That's a developer issue, and not a user issue, so it shouldn't propagate beyond the solver. Anyways, it has an easy fix: the first thing you do in your methods is you put the two in an immutable which holds both types like you currently have it. Done. As noted above, the opposite way does not have an easy fix, and imposes a structure on everything even when it doesn't make sense.

As for the names, we are pretty much settled on Problem, except you abbreviate it in IVP. Julia types usually are not abbreviations, see the Style Guide and check out Base conventions (in Base source). I think that is a good reason to not use IVP and instead write out Problem. I think ODEProblem is clearly the IVP, and we can make ODEBoundaryProblem for the BVP. This keeps going as HeatProblem is clearly the problem for the Heat Equation, SDEProblem is clearly for SDEs whereas SIVP is not an abbreviation which is used. These are the justifications for the Problem naming scheme.

As for the solution, I think it's very clear to keep the Problem from the Solution. Again, it's a developer-side issue that for iterator schemes this ends up being the same thing because you're just mutating the Problem-type. But by keeping Problem and Solution different, we can have different interfaces. This is important because you can require parts of the Solution interface for certain functionality, for example a full timeseries for plot recipes and an interpolation scheme for parameter estimation. The other problem is that different solvers need different solution objects. Sundials saves its memory quite differently, which makes overloading the interpolant very different than RK methods where you have full access to the k vector. Making different Problem types with different fields to initialize for different solvers again adds unnecessary user burden when there is a simple solution: keep Problem separate from Solution.

By keeping these separate, the documentation is also more clear. We can document the Problem interface and the Solution interface. People who only write new parameter estimation techniques need to know the Solution interface well, but not necessarily the Problem interface.

So I think it's very clear that Problem and Solution is here to stay. It's easy to document, extends to other problem domains well, and from all of the feedback I've had, is very intuitive. However, that doesn't mean I'm not conceding to any changes, there are plenty of changes in the PR49 direction that I will be making. Here are some changes and justifications:

• Put tspan in the Problem. This should've been done in the first place. In DiffEqProblemLibrary.jl you can see that in the docstrings I usually mention the timespan for the problems... that makes it very clear I made a mistake.
• Move the options over to some type. It's very clear that this dictionary method was not a good idea. Your usage of types here is much better. However, the true solution is the JuliaML/Optim way of using instances of the types. The reason is because that way types can be immutables which hold their default values (stabilization betas, etc.) and override the defaults. Each standard option will have a boolean for if the default is overwritten. It's the solver's job the handle kwargs seamlessly. This gives algorithm-specific defaults which are general enough for the ecosystem, but also easy on the user.

Let me give a specific example. q is the stepsize modifier, where the next step you adapt to q*dt (this is Hairer's notation?). It's common to have a qmax. It's also common for that qmax to be dependent on the algorithm: DP5 uses 10.0 while DP8 uses 6.0 following the Hairer convention / tests. 10.0 is a standard default for most algorithms. So the logic goes as follows:

• If the user hasn't set a qmax, then the option type is made with qmax=10.0 and qmax_set=false. DP8 has a field qmax=6.0, and has a function set_alg_defaults!(opts,alg) which sees qmax_set=false and replaces (or makes a new option type, depending on if it's mutable) the value with qmax=10.0.
• If the user passes qmax=10.0 as a kwarg, then the solver does a generic dispatch where it reads the kwargs and makes an option type from then, then calls solver with the options. Thus algorithm developers don't need to worry about kwarg support (which can get complicated due to aliasing).
• For the case where the user passed the kwarg, qmax=10.0 and qmax_set=true. Inside set_alg_defaults(opts,alg), it sees qmax_set=true and thus does not override the value.

This solves three issues:

1. It makes it so that way algorithm defaults are optional. Less issues for the developers to handle.
2. It handles aliasing at the highest level. This keeps things simple for those who aren't aware with how the terms are used throughout the ecosystem, both developers and users.
3. It gives a way for algorithm defaults to not override user input, even if the user input is the global default.

set_alg_defaults will just have to be part of the advanced solver interface, but many people can ignore it. Sane defaults will be set in the option definition (i.e. matching the MATLAB/SciPy/ODEInterface/Sundials/LSODE/DifferentialEquations conventions most people know and which are well tested). The last thing is that the options type will have a dictionary at the end for nonstandard / algorithm-specific options.

There's no need to have a FixedOptions and AdaptiveOptions. There's an adaptive option to turn on/off adaptive timestepping, and if you don't need some options, just ignore them. There's not really a cost to ignoring options, but there's a big conceptual/implementation cost to splitting apart and having multiple options types. Instead just ODEOptions, SDEOptions, etc. since between problem domains these are very different. I think it's much easier to have a pretty encompassing set of options which "most algorithms use" and have a chart mentioning what's available to each algorithm

[ChrisRackauckas]

PR49 was never about implementing newer or faster solvers but more about designing a flexible and simple interface, so if the plan is to wrap it up with another layer of interfaces anyway, then maybe there no point in having PR49 at all.

While I think it was an interesting experiment, I don't think that it's flexible or simple. In fact, I think it's main problem is that it can be very misleading. Here's my full take on it.

Flexible

It is flexible, but not anymore flexible than callback functions. Essentially when you write:

for (t,y) in ODE.iterate(ode)
    @show (t,y)
    if t > 1
        break
    end
end

what you're doing is making the algorithm have a callback function

function callback(t,y)
    @show (t,y)
    if t > 1
        break
    end
end

which occurs after each iteration. Not only is that strictly less flexible in ability than having the saving behavior in the callback (in the iterator interface, you need the temps part of the object so you have some form of saving already occurring, and this causes problems when you have immutables unless you make the type mutable which then reduces performance), but it's not as flexible in algorithm-choice: only iterator style algorithms can do this. This excludes many algorithms. Sure, you can put an iterator interface on Sundials, but you can't on ODEInterface. And if someone writes GPU algorithms? If you want to port more C algorithms? They will all allow for callbacks, but not iterator interfaces.

But while callbacks are strictly more capable than the iterator interface, you can still use a callback function with the iterator interface. Because you have the iterator interface, I can define the solver to use

for (t,y) in ODE.iterate(ode)
   callback(t,y)
end

and then if someone passes a callback function, it will work with PR49.

So while it's more flexible than what ODE had before and makes it easy to setup ODE with callbacks, it's still less flexible than callbacks and so this is not a point that's pro iterators.

Simple

No. I agree with you that

I think that implementing a number of things to do when implementing a new solver should be kept down to a minimum

That means an iterator interface shouldn't be required. It takes a good amount of Julia knowledge to really make the iterator interface work for something as complex as solving differential equations. Indeed, I still don't understand some of the choices you had to make in full detail. The choice of an iterator interface that people have to write towards would exclude the simplest way to implement methods: a for loop. If we want to make the interface accessible to those who might have very little development/Julia knowledge, then they should just be able to write a loop which uses the input to compute the formula from the paper. Someone should be able to write an algorithm as:

  @ode_preamble
  halfΔt::tType = Δt/2
  local du::rateType
  while t < T
      k = f(t+halfΔt,u+halfΔt*f(t,u))
      u = u + Δt*k
      t += Δt
  end
  return u,t

and have it work (replacing the preamble with explicit code of basically saying f=prob.f etc.). That's clearly the midpoint method. Does this version have callbacks, adaptive timestepping, etc.? No, but it should work. This will allow people to contribute algorithms for things like PDEs without knowing much Julia development. Someone can come later and add to it, but I think allowing the simplest algorithm should be allowed, and that's not an iterator interface on a custom type, it's a loop.

The iterator interface does make it easy to add more advanced features (for example, the callbacks), but that's an advanced developer issue, not something to complicte both users and not as advanced developers.

Not General

Not only is the iterator interface not flexible, but it's also not general. While all algorithms iterate in some form, the iterations can mean something very different in different algorithms. The iterator interface you have works well in time, but not all algorithms do that. Some algorithms don't necessary step foward in time, and instead iterates a vector of steps torwards a solution (some parareal methods). BVP tends to take this form as well. Some SDE methods don't necessarily treat time always in a foward motion. So restricting the interface to iterating (t,y) is not general, and you'd have to say "well, for some algorithms you do this... and for others...". But they can all use solve. And whether you have to do callbacks differently (or what callbacks even mean for some algorithms!) is a higher order issue with extra documentation that should only need to be addressed if the user has specific requirements.

Misleading

The iterator interface allows for misleadingly simple callbacks. By misleading I mean, it lets you easily do bad things. In your manual you already show one:

for (t,y) in ODE.iterate(ode)
    @show (t,y)
    if t > 1
        break
    end
end

if t>1 isn't good because you step over t=1. So then you backtrack and say "nonono, don't do that!". And this is precisely the problem. Anytime you want to do event handling like this, you most likely want to use interpolations, root finding, etc. to get it right. A simple check like x<0 will hit the wrong time and worse, it will not even find the events if there's oscillations (every timestep can be positive, with it negative inbetween!). It's even worse for stochasticity. I have a research project in event handling for stochasticity. Think about this: if you're very close to x=0, then you "probably hit it already" in the stochastic case. So what's the distribution you should use to pull back your time?

This is not easy stuff to do correctly, and by telling them "oh you define events right here!" is leaving the correctness as user-burden. A user has to know proper event handling algorithms to do this correctly. I think that most people who develop differential equation algorithms don't even know this stuff! Yet this interface leaves it ripe for abuse. If the discontinuties are not handled appropriately, you lose all of your accuracy (if it's handled with a conditional as you show, then every algorithm is first order O(dt)!). You need to find the the point of the event, backtrack the solver to that time, discard whatever was after and save the new value, and decide how to continue. This is not shown in your example, and it won't be easy to do in general (different algorithms have differerent interpolation schemes. Just showing the way to do it for Hermite interpolation is really bad because it's 3rd order with a large constant for the error!)

This is why event handling really needs a DSL. There is a lot of functionality that you want to make simple, but is actually much more complex. This still isn't very stable, but check out the docs for it. How to handle the DSL for different packages is still an issue that needs to be worked out (if other packages get as far as event handling) since currently the macros make callbacks specific for OrdinaryDiffEq (likely, the callbacks/macros will have to be package-specific. Making it work for all OrdinaryDiffEq algorithms was already a feat which pushed the limits of what's doable).

Performance

Packing everything into a type (as is required by the iterator interface) makes you use more mutable structures for when things change. It's much easier to write an algorithm using only stack-allocated immutable, limited function calls (and testing that each function inlines), etc. when it's not in an iterator interface. Julia optimizes the standard loops very well, and so we should allow someone to just write the loop!

Conclusion

There are a lot of downsides to the iterator interface. It's a Julia-specific interface which only makes sense in some cases (though those cases are more normal than not), but also is much more complex, can hinder performance, and makes it extremely easy to destroy the accuracy of the methods. It also adds more user-burden since the iterator interface is undeniably more complex than "solve the problem", and there's developer-burden if it's made into a requirement (some algorithms I'm not quite sure how I'd even do it!). If there wasn't a simple solution which doesn't have these problems, I would go with the iterator interface. But there is a simple solution which doesn't have these problems and is perfectly compatible with non-Julia codes: callback functions.

That doesn't mean there aren't any pros to the iterator interface. It's clear that, if you design your set of algorithms with an iterator interface, it's easy to put a common footer on each algorithm and uniformly implement callbacks. Also, by folding things into types, "algorithm-switching" (detecting stiffness and switching to a different algorithm) can be easier to implement (indeed, this research path is why I looked so deeply into iterators myself). However, there are too many pervasive and hard to avoid negatives that I cannot recommend that it be the common interface. I think it was a good experiment, but the conclusion is that the iterator interfaces work much better on machine learning / optimization problems, and iterative solvers of Ax=b, than they do for differential equations.