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:

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

[ChrisRackauckas]

I'm not quite sure what you mean with below statement, could you post an example:

You should be able to do something like const f = prob.f, though it may need a const f::typeof(prob.f) = prob.f. Then it should know how to compile the rest of the code knowing f is strictly typed while prob.f is not. But this is not ideal or developer friendly. (But even this may not work, though it's essentially simulating using a function barrier so I would suggest always using the function barrier instead).

For the functions which have super long compile times, would there be a trick to let them precompile but not specialize. Probably involving ANY.

The trick is to put them in a type and use a "reverse function barrier", i.e. put it in a type which is not strictly typed on the function, pass this through the function barrier, and then pull the function out from that type. Then only the top-level interface has to recompile for each function (which is a trivial cost). But as I said earlier, I'm doubtful that this has any benefits because the places where this would really cause a large difference in compile times are for very specialized functions which only make sense in long running high-precision applications (like the Feagin methods).

Making all of this strictly typed could help @fastmath too.

[mauro3]

I tried here writing const ff::typeof(prob.f) = prob.f but it did not change the timings, whilst removing the T from here still got me back to the original timings (edit: i.e. the slowest times). I guess this may have to do with https://github.com/JuliaLang/julia/issues/5148.

[ChrisRackauckas]

Ahh okay, even more of a trap than I thought.

[mauro3]

My experience is that type annotations generally help a bit but not 100%.

I'll try and look through your other posts later today.

[ChrisRackauckas]

Okay cool. Maybe this should've been split into different posts on the different parts, but I think we're nearing the end.

[ChrisRackauckas]

My experience is that type annotations generally help a bit but not 100%.

Annotations are very different from declarations. Declarations help a lot more at compile time in these edge cases. But yes, the easiest is just to strictly type the function.

[mauro3]

I should be more careful: the T here is a type assertion and that seems to help a bit. Adding a type declaration like local ff::typeof(prob.f) = prob.f here does not help. Is that what you mean with "declarations"?

[ChrisRackauckas]

Adding a type declaration like local ff::typeof(prob.f) = prob.f here does not help.

Surprised that doesn't help with functions because it will help with other types that I know of. That must be because the type of the function can only be gotten at runtime.

[pwl]

Sorry for the delay, but I finally got to read up on the Solution type. I like most of the ideas, so I will only comment on the things that I think we could improve on.

1. This looks a little inconsistent to me

sol[i] # Gets the ith value of u, the saved timeseries
sol.t[i] # Gets the timepoint of the ith save

A single step consists of a pair (t[i],u[i]), yet you access u in a different way then t (u[i]=sol[i] and t[i]=sol.t[i]). What do you think about the following interface?

sol[i] # Gets the ith pair (t[i],u[i])
sol.u[i] # Gets the ith value of u
sol.t[i] # Gets the ith timepoint

This way you can access u and t as arrays (by directly referring to sol.t) and as a bonus you can easily get a pair (t[i],u[i]) that you could iterate over. I wrote the next points to be separate but they could be viewed as arguments for supporting this interface.

1. For some applications you need to save the derivative du, is there any room for having it in the Solution? It might sound trivial to get du for ODEs (you just use du=F(t,u)) but it's an additional computational effort, which was performend by the integrator anyway. This issue becomes nontrivial if we consider DAEs, for which obtaining du requires solving a nonlinear equation 0=F(t,u,du). I'm not saying we should save du every time but there should be a standardized optional way of saving and accessing it. If we add du as a type field it would be reasonable to return it as a part of the tuple sol[i]=(t[i],u[i],du[i]) (this might potentially be a problem if du was not saved, in that case we could simply return (t[i],u[i],nothing) or something like that). The advantage of having sol[i] returning a tuple is that one can simply ignore additional terms and write [(...) for (t,u) in sol] and it would work independently of whether the solver saves du or not.
2. (related to 2.) is there any way to solve additional per-step data for the stepper? I guess this could be covered by a callback, but we could potentially have an integrator_data::Vector{MyIntegratorData}, which would be optionally allocated by the integrator. This could be useful for testing integrators or debugging and we could also return it with sol[i]=(t[i],u[i],du[i],data[i]).
3. As with the ODEProblem type I would prefer if we kept the analytic solution and the errors in a separate Solution type (TestSolution?).

[ChrisRackauckas]

As with the ODEProblem type I would prefer if we kept the analytic solution and the errors in a separate Solution type (TestSolution?).

I saw that coming and split the types already. It's no build deal, and I made a build_ode_solution which dispatches off of prob to build the right one.

(related to 2.) is there any way to solve additional per-step data for the stepper? I guess this could be covered by a callback, but we could potentially have an integrator_data::Vector{MyIntegratorData}, which would be optionally allocated by the integrator. This could be useful for testing integrators or debugging and we could also return it with sol[i]=(t[i],u[i],du[i],data[i])

If the solution type is non-concrete and this is just a template simple solution type, then you can add concrete instances however you want. All of the general tools (plotting, parameter estimation, etc.) only assume the Solution interface as discussed, but if you make a concrete subtype of AbstractODESolution you can add extra fields as necessary and still get these features (which is what I plan to do with Sundials. Note that actually AbstactSDESolution <:AbstactODESolution and AbstactDAESolution <: AbstactODESolution and, when we get there, I think it would make sense to have AbstactDDESolution <: AbstactODESolution.

For some applications you need to save the derivative du, is there any room for having it in the Solution? It might sound trivial to get du for ODEs (you just use du=F(t,u)) but it's an additional computational effort, which was performend by the integrator anyway. This issue becomes nontrivial if we consider DAEs, for which obtaining du requires solving a nonlinear equation 0=F(t,u,du). I'm not saying we should save du every time but there should be a standardized optional way of saving and accessing it. If we add du as a type field it would be reasonable to return it as a part of the tuple sol[i]=(t[i],u[i],du[i]) (this might potentially be a problem if du was not saved, in that case we could simply return (t[i],u[i],nothing) or something like that). The advantage of having sol[i] returning a tuple is that one can simply ignore additional terms and write [(...) for (t,u) in sol] and it would work independently of whether the solver saves du or not.

Sure, we can add a du field. We will need to add an attribute to turn on/off its saving: in many PDEs things get memory-bound and I think many might not care about the derivative. save_derivative? I actually kind of had this in there as k (or k[i][1] depending on the algorithm), which should probably be given the more informative name interp_data.

What do you think about the following interface?

The advantage of having sol[i] returning a tuple is that one can simply ignore additional terms and write [(...) for (t,u) in sol] and it would work independently of whether the solver saves du or not.

I understand the consistency, and I understand that this looks kinda nice, but I don't know of a case where this is actually convenient. For plotting you want to keep the vectors. For parameter estimation you want to compare the interpolated solution. For sensitivity analysis, you want to break apart the u to the component systems. For PDEs you want to do extensive analysis on just the u. And when working with PR49, I had to specifically ignore that interface because it always got in the way (just about everything wants the vectors!). This makes the abstract array interface on sol an array of tuples, not an array of the solution values, which adds the cost of repeatedly building thousands of tuples and adding dereferencing costs, and I'm not sure of what you'd get.

On top of that, I am pretty sure that if you tried to do all of the nice Julia vectorized things over this, such as broadcast, exp.(sol) wouldn't work anymore? Also, returning these tuples would be inconsistent with the interpolations, unless you recommend the interpolation building a tuple (t,u,du,data), so now do we also have to include du as a standard part of the interpolation? While it would be a nice option to add to the interp interface (interp(t,k) for the kth derivative), making it standard wouldn't make sense (since the most used part of this feature is plotting and checking against data!), but now we have a tuple which sometimes has a bunch of blank spots, and other times not. If your tuple seems to have blank spots more than half of the time, then maybe the tuple isn't the right idea and working on the index is. The comprehension could still be done with indices, and that would be more flexible to including / excluding parts of the type.

Mind you, I am not against it because I would have to re-write some things. It wouldn't be hard to change every instance of sol[...] to sol.u[...] (though the change in the array interface might break a few things and need to turn them into a loop). However, it would change sol[...] from something that is widely used throughout all of the diffeq things to something that isn't useful in any of the tools, which is an indication that it might not be a good idea.

The major point of the interface was to make working with Vector{Vector{T}} stuff easier (that's one big thing gained by the plot recipes). The proposal of making this iterator essentially a Vector{Tuple{Vector{tType},Vector{uType}}} makes the iterator so tough to work with that there's almost no point in doing so, except for the one case where you want to do a comprehension and don't want to touch indices.

Another major point of the interface is that it makes the variable naming agnostic. I use u, but I know others who stay only in ODEs tend to use y. Notice that the entire interface, from defining the problem object to analyzing the solution, allows the user to not ever be forced to u. It's all naming agnostic. But if you want to get a vector out and the solution make a vector of tuples, you will need to sol.u.

[ChrisRackauckas]

One last thing, I the comprehensions on the solution will not work well anyways.

[(...) for s in sol]

returns the solution with tslocation mutated to match the state. This is the only way that I know of to make the animations on the iterator interface still use the plot recipe (otherwise we'd have to special case all of the plotting, and that seemed like a huge problem to solve with an easier solution: just use the plot recipe with a tslocation in the solution!). We would need a good solution for this in order to drop the "iterator for animation interface", because otherwise we get a Plots dependency which doesn't play nice (due to conditional deps, you've probably read about this everywhere already).

Meanwhile,

[@show (sol.t[i],sol[i]) for i in eachindex(sol)]

always works not matter how you want to change it around, add things, etc.

[ChrisRackauckas]

I went around and asked my lab and had a bunch of Gitter discussions and it seems agreed upon that sol[i] should stay sol.u[i]. Another thing I forgot to mention is that it iterates like sol[i,j] == sol.u[i][j] and so it "acts like the matrix output" people are familiar with from other DE solvers (like Sundials, Fortran codes, etc.).

That doesn't fix your iterating issue. Thankfully, we found a solution.

[@show (t,y) for (t,y) in tuple.(sol.t,sol.u)]

Works really nicely and is extendable for whatever you want to put in there. I defined tuples(sol::DESolution = tuple.(sol.t,sol.u) as the standard case, so

[@show (t,y) for (t,y) in tuples(sol)]

works. That fixes both the comprehensions, keeps the interface for those who like the matrix style, allows for animate(sol) to work, keeps it variable-name agnostic (notice this uses y the whole time, you can do that if you want. Just the interface for developers needs u as the field name), and still allows for performant plotting. I think that puts this issue to rest.

As for du, someone rightfully pointed out that in many cases it doesn't make sense. For example, most SDEs are not differentiable anywhere, and many PDEs are not differentiable most of the time. So instead du should be specific for certain types of equations. What we should do is document specific extra parts of the solution which can be available along with the problem objects. So DAEs and ODEs can have a du as given by save_derivative, but it's not part of the standard solver interface. This also makes sense for many PDEs where we might want sol.x, sol.y, etc. to be meaningful.

[mauro3]

Thanks @ChrisRackauckas for the big write ups. I haven't thought that much about solution types yet so my input is very limited. Whilst @pwl's suggestion is "purer", I'm firmly on @ChrisRackauckas side that sol[i] should just index u. I love that sol[:,j] just gives the time series of the j-th component. The rest looks good (without having spent a great deal pondering it). Making dense a type parameter seems good, maybe the two :dense and :points could work. A bool is a bit cryptic.

Jacobians & Parameterized Functions

Ideally, we should just go with overloaded functions as we don't need any of the bells and whistles. This is what every julia user knows. So tell them if they want to provide a Jacobian, then define f(::Val{Jacobian}, ...). However, if would we want to internally add the Jacobian, say with ForwardDiff, then that doesn't work:

julia> function fn(f::Function, x)
       f(::Val{:a}, uu) = uu.^2
       end
ERROR: syntax: cannot add method to function argument f

Is this an issue or not? Does this work with the call overload feature?

DAE problem types

I think we could just have separate types for implicitly defined ODE/DAEs and ones with a mass matrix. It should then be possible to set up some promotion system to convert the less general ones to the more general ones.

Misc

The isoutofdomain I introduced for adaptive time stepping. Say you got a ODE for which y=>0, then that can be encoded in isoutofdomain and force the stepper to take a smaller step if y<0 occurs.

[ChrisRackauckas]

Making dense a type parameter seems good, maybe the two :dense and :points could work. A bool is a bit cryptic.

This isn't type stable for ODEs (for equations for which interpolations don't exist it is, but if the user can turn density on/off (which is necessary for memory management in PDEs), then it is not). I think that amount of over-typing can lead to headaches. (What would :points be for? We don't need a boolean for if there are intermediate steps saved in the timeseries: you can just check the size. The interpolation is different since to check it, you have to try interpolating and catch the error).

I think we could just have separate types for implicitly defined ODE/DAEs and ones with a mass matrix. It should then be possible to set up some promotion system to convert the less general ones to the more general ones.

I didn't even think about a promotion system. That makes perfect sense and would automatically make the DAE methods implicit solvers for ODEs. Since the solution of a DAE subtypes and AbstractODESolution, the plots and everything would still be the same.

The isoutofdomain I introduced for adaptive time stepping. Say you got a ODE for which y=>0, then that can be encoded in isoutofdomain and force the stepper to take a smaller step if y<0 occurs.

Huh, I've never heard of this. I assume if it's out of the domain you just count it as a standard rejection? I.e. if q>1 || isoutofdomain, reject? Do you have a reference on how this affects solution properties?

[ChrisRackauckas]

I love that sol[:,j] just gives the time series of the j-th component.

I actually was missing that one but just added it. It splats so if your independent variable has shape you can get the timeseries by that shape too. So if it's a matrix, sol[:,4,2] is the timeseries of the [4,2] coordinate. If you just do sol[:,3], it will get the timeseries of the 3rd component by linear indexing. Don't know how I missed this dispatch for getindex because now after you mention it, this is super useful!

[ChrisRackauckas]

I changed ParameterizedFunctions.jl to have the following interface:

f.a # accesses the parameter a
f(t,u,du) # Call the function
f(t,u,du,params) # Call the function to calculate with parameters params (vector)
f(Val{:a},t,u,2.0,du) # Call the explicit parameter function with a=2.0
f(Val{:a},Val{:Deriv},t,u,2.0,df) # Call the explicit parameter derivative function with a=2.0
f(Val{:param_Jac},t,u,J,params) # Call the explicit parameter Jacobian function
f(Val{:Jac},t,u,J) # Call the explicit Jacobian function
f(Val{:InvJac},t,u,iJ) # Call the explicit Inverse Jacobian function
f(Val{:Hes},t,u,H) # Call the explicit Hessian function
f(Val{:InvHes},t,u,iH) # Call the explicit Inverse Hessian function

with overload checkers:

jac_exists(f)
invjac_exists(f)
hes_exists(f)
invhes_exists(f)
paramjac_exists(f)
pfunc_exists(f)
pderiv_exists(f)

using a method similar to what @mauro3 described (all those booleans are gone now, that's really nice!). Obviously, the solver interface can use any overloaded function, but this is the most complete interface (so far). The booleans require that a user types the first argument of the function if you make an explicit parameter dispatch, since f(t,u,du,params) comes up as all Any and so the dispatches for Jacobians and the like catch it. That's such a rare non-issue though.

On that note:

Ideally, we should just go with overloaded functions as we don't need any of the bells and whistles. This is what every julia user knows. So tell them if they want to provide a Jacobian, then define f(::Val{Jacobian}, ...). However, if would we want to internally add the Jacobian, say with ForwardDiff, then that doesn't work:

Is this an issue or not? Does this work with the call overload feature?

immutable Tester end
function add_dispatches()
  (t::Tester)(u) = 2
end
add_dispatches()
t = Tester()
t(3)

That's why I say overloaded types. It's literally the same thing (it's how anonymous functions are internally made) except it works in more crazy cases. Of course, overloaded functions are the same exact thing, so if you say that you think it should be overloaded functions, then we are agreeing. The solver would interact with the variable the same way, and the user can make it however they choose.

[mauro3]

A few comments:

• yes, let's not make dense a type parameter. However, instead of having a field dense, how about just leaving interp undefined when there is no dense output?
• isoutofdomain: yes, it just rejects the step, just like finding a NaN also rejects the step. I thought I saw it in Hairer & Wanner but couldn't find it. Maybe I just made it up. Intuitively, I would have thought that it could only help, but I don't have proof.
• above I should have told you about potential problems with method_exists and parameterized functions. Check out https://github.com/JuliaLang/julia/issues/8959. I hope this does not impact you, otherwise you'll have to revert to the old system for now.
• Is this right: we could let people just use normal overloading for providing the jacobian. If they want extra features then they could use the parameterized function interface. if so I'm for going with that option, barring any show stoppers from https://github.com/JuliaLang/julia/issues/8959.

[pwl]

@ChrisRackauckas +1 for asking the users. After giving it some thought, I like your proposed interface for Solution. It's slightly less consistent but it offers more functionality and is more user friendly.

I think we could just have separate types for implicitly defined ODE/DAEs and ones with a mass matrix. It should then be possible to set up some promotion system to convert the less general ones to the more general ones.

I totally agree, it seems like a better option then having everything crammed into a single type (like what we have done for PR49). It might require adding some redundant code but this shouldn't be a big issue.

Is this right: we could let people just use normal overloading for providing the jacobian. If they want extra features then they could use the parameterized function interface. if so I'm for going with that option, barring any show stoppers from JuliaLang/julia#8959.

I got a little bit lost in the ParameterizedFunctions.jl discussion but this question suggests that ParameterizedFunction is now the only option to define an ODEProblem. Is that so? Will we still be able to define the right hand side as a plain function?

[ChrisRackauckas]

yes, let's not make dense a type parameter. However, instead of having a field dense, how about just leaving interp undefined when there is no dense output?

How do you leave a field as undefined? It looks like you have to do it with an inner constructor which chops off the last few fields? So you'd have to have an ordering on "if this is undefined, so are the ones further down", you'd have to set the field orderings to that, and you'd have to have different inner constructors for different cases? I don't think this will scale well, and a boolean is a pretty simple option that does scale. Both using a boolean and an undefined reference require runtime checks too, so there shouldn't be a performance difference.

above I should have told you about potential problems with method_exists and parameterized functions. Check out JuliaLang/julia#8959. I hope this does not impact you, otherwise you'll have to revert to the old system for now.

Nope this doesn't affect this setup. There is only one "bad interaction" in this setup and it's because if you define the "explicit parameter function" f(t,u,params,du) without any type assertions, then that's (Any,Any,Any,Any) so the other method_exists which are looking at a 4-argument type signature come up as true because of the full-Any dispatch. That seems like such a rare case (how many people will make that function by hand?) to simply ignore and leave in documentation.

I got a little bit lost in the ParameterizedFunctions.jl discussion but this question suggests that ParameterizedFunction is now the only option to define an ODEProblem. Is that so? Will we still be able to define the right hand side as a plain function?

Of course! You just define a function f(t,u) or f(t,u,du) and put that in the ODE and it will work. However, there's now an extended interface for adding / checking for extra functions via dispatches. So for example the interface for Jacobians is f(Val{:Jac},t,u,J). If you make a dispatch on that signature, for example:

function f(::Type{Val{:Jac}},t,u,J)
  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

Then f contains both the function and the Jacobian, and the developers can check if the user defined the Jacobian with jac_exists(f). However, by using dispatch like this it allows for a much richer set of functions than just Jacobians like traditional solvers were limited to. ParameterizedFunctions.jl uses SymEngine to symbolically calculate and build all of the available "extra functions", but adding as many or as few dispatches as you can/want by hand will have the same effect.

It's just dispatch, Julia's core feature, to give an interface which doesn't include a field for the 10's of special "extra functions" which could exist to speedup calculations. I am building ParameterizedFunctions.jl to make this easy for users to have it (and things like the inverse Jacobian) symbolically calculated so all of these performance improvements are automatic (and so it's something which is much more important than in other languages: now in our ecosystem defining a function with a Jacobian or Inverse Jacobian is rather trivial for the standard system of ODEs and so you really should make sure you have a way of using user Jacobians!). However, this direction has limitations, for example it's not clear how to extend it to most types of PDEs (I say most because I have some tricks up my sleeve for a few PDEs which are important to me), but since this method still allows you to add Jacobians by hand (in the natural Julian way that asserts their relation) that isn't a problem.

isoutofdomain: yes, it just rejects the step, just like finding a NaN also rejects the step. I thought I saw it in Hairer & Wanner but couldn't find it. Maybe I just made it up. Intuitively, I would have thought that it could only help, but I don't have proof.

Okay. I would explain it as rejecting the step: otherwise saying that it "reduces the timestep" makes it sound vague and leaves open the question of "how exactly?". I would not use this to reject a NaN because in any case where you've unstably exploded to a NaN, you've already take like 5 steps in that direction, meaning that rejecting a timestep will never help pull the solution back in. Instead, it would just keep going off in the wrong direction and you'd be delaying the inevitable failure. You should just throw the failure immediately in this situation.

[In fact, in a lot of the benchmarks on difficult problems PR49 fails due to hitting the maximum iterations even when it's set very large. Do you have isoutofdomain false on NaNs by default? Because that would explain a lot.]

However, the case of "keeping it positive" is interesting and intuitive. Since it doesn't hurt, I'll add this to the common interface and I'll default to just using @inline isoutofdomain(t,u)=false.

[ChrisRackauckas]

Question: if we aren't instantiating the type, how are we switching compile-time information?

const T = Rosenbrock23{true}
T.parameters[1] == true

function get_type_parameter{bool}(T::Type{Rosenbrock23{bool}})
  bool == true
end
@code_llvm get_type_parameter(T)

T2 = Rosenbrock23
T2.parameters[1] == true
@code_llvm get_type_parameter(T2)

The second fails. Would I have to require that the user always give a type parameter? That seems not very good. It seems that passing Rosenbrock23() which has a default constructor set a type parameter, and then changing it with Rosenbrock23(false) is easier? This is needed to switch on/off different parts of algorithms such as ForwardDiff and leaving the branches in the code (and not compiling them out) can lead to huge type-instabilities (so autodiff needs to not be a keyword argument and instead be part of the algorithm choice, which makes sense because you want it to compile a different algorithm!)

Also, the jac_exists from before isn't known at compile time because the methods table is a global:

using DiffEqBase
f = prob.f
function testing(f)
  if jac_exists(f)
    a = 2
  else
    a = 3.0
  end
  a
end
@code_llvm testing(f)

So I instead suggest that we should have boolean parameters in ParameterizedFunctions, and for convenience make typealias FunctionWithJacobian ParameterizedFunctions{true,false,false,...,false} which then means that you just have to use overloading and specify the subtype:

immutable MyFuncType <: FunctionWithJacobian end
function (f::MyFuncType)(t,u,du)
  du[1] = u[1] - u[1]*u[2]
  du[2] = -u[2] + u[1]*u[2]
end
function (f::MyFuncType)(::Type{Val{:Jac}},t,u,J)
  J[1,1] = 1 - u[2]
  J[1,2] = -u[1]
  J[2,1] = u[2]
  J[2,2] = u[1] - 1
  nothing
end

Still very standard Julia overloading, just now this is able to fix all the branches at compile-time and make things type-stable.

Thoughts?

[mauro3]

Concerning "dense": as the interp field is a function, would it benefit performance to get a type parameter? If so, we could set interp=nothing and get a different type then when interp=some_function. Otherwise, let's go with the bool-field.

[In fact, in a lot of the benchmarks on difficult problems PR49 fails due to hitting the maximum iterations even when it's set very large. Do you have isoutofdomain false on NaNs by default? Because that would explain a lot.]

I'll check!

Concerning your last post, if I understand correctly, I'd just use typealias. Is the type-para for dense output? In PR49 we have the RK methods all being the same type with different first two parameters (see), e.g. RKIntegratorAdaptive{:dopri5}. But one could/should probably make a full alias, at least for the common methods, e.g. RKdopri5.

[ChrisRackauckas]

Concerning "dense": as the interp field is a function, would it benefit performance to get a type parameter? If so, we could set interp=nothing and get a different type then when interp=some_function. Otherwise, let's go with the bool-field.

That would make the type always be uninferrable as well (you'd have to know the value of dense to know the type parameters) and so there wouldn't be much of a performance gain (while checking types for parameters in abstract cases is a hassle: typeof(sol).parameters[i] and search? Ouch). If you then wanted to dispatch on the solution's type inside of a function, you'd have to pay the price of dynamic dispatch (along with a new compile) each time. I would expect this to only decrease performance (of course, you could benchmark it).

Concerning your last post, if I understand correctly, I'd just use typealias. Is the type-para for dense output? In PR49 we have the RK methods all being the same type with different first two parameters (see), e.g. RKIntegratorAdaptive{:dopri5}. But one could/should probably make a full alias, at least for the common methods, e.g. RKdopri5.

How would you handle the cases for Rosenbrock effectively? Autonomous and not autonomous, autodifferentiation or not, choice of linear solver, etc.? I guess make the most common choices type aliases, and then document the full parameter set? Can you alias no-parameters to something? Rosenbrock23 -> Rosenbrock23{true,false,:cg}? Type-aliasing sounds like it's almost there, and it's probably more powerful than I know how to use, so maybe this does have a solution.

[ChrisRackauckas]

And any thoughts on the Jacobians / type-stability?

I think it's very necessary for the Jacobians and algorithms to carry these things as type information because it allows the branches to compile away (which in many cases results in type-stability whereas having a bunch of existing branches could cause types to be un-inferrable even if it's clear what they should be!). Since the methods table is a global and not inferrable, I can't think of anything else other than requiring the Jacobians to be done by call-overloading and aliasing the FunctionWithJacobian to a ParameterizedFunction subtype, or something like that.

[mauro3]

How would you handle the cases for Rosenbrock effectively? Autonomous and not autonomous, autodifferentiation or not, choice of linear solver, etc.? I guess make the most common choices type aliases, and then document the full parameter set? Can you alias no-parameters to something?

I think those things should be keyword options (unless you think one of those options makes a completely different solver). Inside solve you can then do whatever you want with those kwargs, make them fields of an options-type, turn them into type-parameters of the solver type, ignore them, etc. I don't think a user should need to keep track of what "option" goes where, just kwarg all of them.

[mauro3]

while checking types for parameters in abstract cases is a hassle: typeof(sol).parameters[i] and search

The proper way of doing this just got documented: https://github.com/JuliaLang/julia/pull/19306/files#diff-78f2e86ca66ebcee42787a87d298c710R451

[ChrisRackauckas]

Won't treating those things as keyword args cause dynamic dispatch?

[mauro3]

Yes. But it was my understanding is that it is ok as using kwargs puts a spanner into efficiency anyway. The idea would be to put a function barrier into solve. The stuff above about init shows a possible way to do that (https://github.com/JuliaDiffEq/Roadmap/issues/5#issuecomment-258791153) and also giving an option to call solve type-stabelly with a instantiated solver type (instantiated via init).

[ChrisRackauckas]

The stuff above about init shows a possible way to do that (#5 (comment)) and also giving an option to call solve type-stabelly with a instantiated solver type (instantiated via init).

Wouldn't that lead to a lot of redundant code? That's basically asking for that large ODEIntegrator type I have to part a part of every type? It also makes it so that way one can't use parametric types to distinguish between related algorithms anyways because there will be so many more parameters, right?

[mauro3]

I will have to look at your code to see (tonight). But it should be fine in general I would have thought.

[ChrisRackauckas]

Okay thanks! Also, reading (a lot of your stuff) on traits: overloaded types having Jacobians should be a trait, etc. I just really need to find a way to make code which gets rid of the other branches when this is the case, and it looks like using has_jacobian(f) with the methods table doesn't cut it. I think the two standing issues are:

1. How to make "Functions with Jacobians" something that's easy to setup and has type information. I think just using the standard call overloading syntax with an alias to a ParameterizedFunction (which would have a lot of booleans) isn't a bad solution.
2. How to get algorithm sub-choices like "autodiff" to be type-information which doesn't cause dynamic dispatch. Pre-instantiating the types would work, using types parameters on algorithms would work, and initializing could work if there's a good clean way to do it.

I think when those two issues are solved we can close this and the remaining problems could be handled in their own threads (callbacks, DAE types, etc.)

[mauro3]

Regarding 1., rather than all the of bools, could you use an approach like https://github.com/JuliaODE/ODE.jl/blob/e47748ba852ebf4f9b1acf4ce0771e577bbe92e5/src/base.jl#L51: i.e. each (possible) function-field has a type parameter. If it is used then its the function type, if not then it's Void.

[ChrisRackauckas]

We could do that. And a separate dispatch to make an ODEProblem from a ParameterizedFunction to put the right functions in the right places. The only problem is, are you sure you want that many fields? The kinds of functions I have already for this is:

1. Jacobian
2. Inverse Jacobian
3. Rosenbrock Inverse-W
4. Rosenbrock Transformed Inverse-W
5. Hessian
6. Inverse Hessian

not including the explicit parameter functions and parameter derivatives which aren't necessary for ODEs. I also plan on added explicit functions for exp(dt*J) (for very fast exponential integrator methods) and who knows what else. Since I have the power to symbolically compute and use macro programming to build these as explicit in-place updating functions (and the domain I work in, Systems Biology, has most of its equations compatible with what I am doing, even in the PDE form after certain transformations), I plan to be researching this area quite a bit since I think it can result in "stupid fast" methods (in the specific cases where it's possible). We could put them all in the ODEProblem and make the ODEProblem have some keyword arguments for these things, but maybe there's a better way, and I'm sure a few months from now you'll be asking for the ODEProblem to be cleaned up with the amount of functions put in there.

[mauro3]

ODEProblem has just one parameter for the function f. That parameter can then be the type of just a function, say typeof(sin), or the complex beast of a ParameterizedFunction. If you add them later, say an analytic jacobian, you will get a different problem instance. I think, this can (ans should) be done internally without the user ever knowing.

Internally, going down the init path, this would mean that init also potentially re-creates the problem instance: p0_, alg = init(p0, Alg; kwargs...). Then, say, typeof(p0.f)==typeof(sin) but typeof(p0_.f)==ParameterizedFunction{...}.

[ChrisRackauckas]

But how do you check if a ParameterizedFunction has a specific dispatch in a type-stable way?

[mauro3]

Could this work:

type ParaF{..., tJ, ...}
...
J::tJ
...
end

# typealias not needed but nice to avoid juggling function-parameters:
typealias ParaFwithJ{...,tJ,...} ParaF{...,tJ,...}
typealias ParaFwithoutJ{...} ParaF{...,Void,...}

f(::ParaFwithJ) = 1
f(::ParaFwithoutJ) = 2

[ChrisRackauckas]

No, the whole reason for using overloading was to have the parameter implicit forms:

(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

Putting functions as fields won't allow this. Also, having the functions as fields isn't any easier to part for existence than putting bools in the type parameter (since that's what you'd essentially be doing here: checking the type parameters).

[mauro3]

I see. Is each ParameterizedFunction its own type? (I think so and your above example suggests it)

If yes, then you could go with traits because then you can add functions later without having to re-do the type because one of the parameters changes.

using SimpleTraits # you can hand code it as well
@traitdef HasJ{X}
@traitdef HasH{X}

immutable ParaF # no parameters
end

...
(p::LotkaVolterra)(Val{:J}, t, u, J) = ...
@traitimpl HasJ{LotkaVolterra}

@traitfn fn(p::ParaF::HasJ, ...) = 1
@traitfn fn(p::ParaF::(!HasJ), ...) = 2

[ChrisRackauckas]

I see. Is each ParameterizedFunction its own type? (I think so and your above example suggests it)

Yes, this is the simplest example of a ParameterizedFunction:

type  LotkaVolterra <: ParameterizedFunction
         a::Float64
         b::Float64
end
f = LotkaVolterra(0.0,0.0)
(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

I guess call-overloading on types isn't too widespread yet (it's something that was introduced in v0.5), but it's the (p::LotkaVolterra)(t,u,du) which allows you to make your type have a call. In v0.5, anonymous functions are defined like:

immutable HiddenType end
(f::HiddenType)(x,y)  = x^2 + y

where it generates a type with some weird unique name. So this is really just a more general way of defining functions (since you now have the whole type, so you can put fields in there).

If yes, then you could go with traits because then you can add functions later without having to re-do the type because one of the parameters changes.

Thanks for showing how to use traits. Should we just embrace them before a Base release? I mean, then the code for defining a Jacobian manually would be:

type  LotkaVolterra <: ParameterizedFunction
         a::Float64
         b::Float64
end
f = LotkaVolterra(0.0,0.0)
(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
(p::LotkaVolterra)(Val{:J}, t, u, J) = ...
@traitimpl HasK{LotkaVolterra}

which I don't think looks bad at all (users who aren't familiar with such overloading would probably just need one example like this, right? Maybe it's one of those things where since I already know it that it looks easy).

What I really like is then:

@traitfn fn(p::ParaF::HasJ, ...) = 1
@traitfn fn(p::ParaF::(!HasJ), ...) = 2

is extremely nice for writing the solvers. This would make it extremely easy to write something nested. For example, in Rosenbrock methods you really want to calculate W^(-1) = (I - gamma*J), and so if I am understanding the traits correctly the Jacobian code could be:

@traitfn calc_jac(p::ParaF::HasJ,t,u,J)    = p(Val{:J}, t, u,J)
@traitfn calc_jac(p::ParaF::(!HasJ),t,u,J,method=Val{:FD}) = # Calculate the Jacobian using ForwardDiff or Calculus.jl?

and then calc_jac could always be used and it will essentially auto-optimize? That's really freaking cool! I would propose a separate interface on this then, having calc_jac in DiffEqBase, and the DiffEq universe could just always use this code agnostic of whether the it's using a user-defined Jacobian or not.

But, if I'm understanding it correctly, it gives really complex things like

# Define Win-calculations
@traitfn function Wsolve(p::ParaF::HasWinv,t,u,out,b,Jcache)
    p(Val{:Winv}, t, u, Jcache)
    @into out = Jcache*b # uses InPlaceOps.jl to make `*` in-place
    Jcache # return Winv
end
@traitfn function Wsolve(p::ParaF::(!HasWinv),t,u,out,b,Jcache)
   calc_jac(p::ParaF::HasJ,t,u,Jcache)
   fac = lufact((I - gamma*Jcache))
   @into out = fac \ b # uses InPlaceOps.jl to make `\` in-place
    fac # return the factorization
end

To solve Wx=b for W = (I - gamma*J) in the most efficient manner given what functions have been defined (without overhead). Is this on the right track. Of course, there's still something to cleanup, (like how we need Wsolve_cached(...) which solves the next steps using either the cached Winv or factorization, etc.), but does traits actually give us this design? Because if that's the case, I am 10000% for traits.

[mauro3]

Concerning https://github.com/JuliaDiffEq/Roadmap/issues/5#issuecomment-260288332: The way we handle it in PR49 is to have the options part of the Alg type, whereas you (if I understand correctly), have the Alg type part of the options type (ODEIntegrator). I don't think that makes that much difference. The advantage of our approach is that each Alg only has the options it can deal with, whereas your approach probably makes unifying the options a bit easier. The init function in your case would then be p0_, odeint = init(p0, Alg; kwargs...) and the kwarg-less solve invokation solve(p0_, odeint). Anyway, I think it should be fine but I may well be missing something.

Concerning your last post: yep, I suspect that this should work. I'm a bit hesitant to make SimpleTraits a requirement, but maybe that would be ok. If someone doesn't want to do trait-dispatch they could always use a provided check has_jac(f) and do hand (or dynamic) dispatch. Am I correct that the common interface wouldn't need to depend on ParameterizedFunctions? It would just specify to provide a Jacobian by overloading f(::Val{:jac},...). (but that would be compatible with ParameterizedFunctions)

Note that this example of yours

@traitfn calc_jac(p::ParaF::HasJ,t,u,J)    = p(Val{:J}, t, u,J)
@traitfn calc_jac(p::ParaF::(!HasJ),t,u,J,method=Val{:FD}) = # Calculate the Jacobian using ForwardDiff or Calculus.jl?

doesn't quite work as dispatch first happens on types as usual. Only then, if a traitmethod is dispatched to, the trait dispatch kicks in. So it would have to be kwargs ...,J;method=Val{:FD}) or probably let the first definition also have the method argument (but just ignore it).

[ChrisRackauckas]

The init function in your case would then be p0, odeint = init(p0, Alg; kwargs...) and the kwarg-less solve invokation solve(p0, odeint). Anyway, I think it should be fine but I may well be missing something.

Oh I see, that will work.

I'm a bit hesitant to make SimpleTraits a requirement, but maybe that would be ok. If someone doesn't want to do trait-dispatch they could always use a provided check has_jac(f) and do hand (or dynamic) dispatch.

Yeah, but if has_jac(f) is defined in DiffEqBase they might as well use it. It would require people who define their Jacobians by hand to use SimpleTraits though.

Am I correct that the common interface wouldn't need to depend on ParameterizedFunctions? It would just specify to provide a Jacobian by overloading f(::Val{:jac},...). (but that would be compatible with ParameterizedFunctions)

Yes. ParameterizedFunctions.jl just has macros to make a ParameterizedFunction, but the abstract type is in DiffEqBase. So all code can be written without ParameterizedFunctions.jl and just DiffEqBase.

doesn't quite work as dispatch first happens on types as usual. Only then, if a traitmethod is dispatched to, the trait dispatch kicks in. So it would have to be kwargs ...,J;method=Val{:FD}) or probably let the first definition also have the method argument (but just ignore it).

Oh yeah oops. So yeah, put and ignore the method optional arg in the first.

[mauro3]

Yes. ParameterizedFunctions.jl just has macros to make a ParameterizedFunction, but the abstract type is in DiffEqBase. So all code can be written without ParameterizedFunctions.jl and just DiffEqBase.

Ok, but what about just using straight methods? That should work just fine too even if they are not call-overloaded. Or would someone implementing the DiffEqBase interface need to subtype from ParameterizedFunction?

[ChrisRackauckas]

Ok, but what about just using straight methods?

What do you mean "straight methods"?

[mauro3]

Say I implement a Alg: solve(p::Prob, ::Type{Euler}; kwargs...) = .... Inside the solve method, I do not need to care whether p.f<:Function or f<:ParameterizedFunction. I think for "normal" algorithms, f could be either and I wouldn't need to know which it is. Thus, I think we don't need ParameterizedFunction type part of DiffEqBase. But I might be wrong, if so, I don't understand how. Maybe you can point me to the places where it is needed?

[ChrisRackauckas]

How do you do has_jac(f) without it? The traits example dispatches off of ParameterizedFunction. Using a type-alias like FunctionWithJacobian is just using type parameters (in a hidden and easier style), but still requires subtyping. The methods table is not inferrable since it's using a global. Is there a way to write a function which can check for the dispatches on a generic function at compile time?

[mauro3]

@ChrisRackauckas can you give me access rights to DiffEqBase.jl (and Roadmap whilst you're at it)? I have a branch to push.

[ChrisRackauckas]

@ChrisRackauckas can you give me access rights to DiffEqBase.jl (and Roadmap whilst you're at it)? I have a branch to push.

Done.

[mauro3]

No, that didn't work. I can't push to either nor can I access the settings on the web-page. Can you try again?

[ChrisRackauckas]

Try it.

[mauro3]

Tnx!

Here a function interface draft.

Here an example implementing the interface.

[mauro3]

Thanks, I made two PRs: https://github.com/JuliaDiffEq/DiffEqBase.jl/pull/2 and https://github.com/JuliaDiffEq/DiffEqBase.jl/pull/3, maybe discuss it further there.