ODE/DAE solver interface #5
Comments
|
I am not sure how you'd do this or if it's possible. DAE solvers just require more information: You could make them all |
|
Sorry, I should have been clearer, I updated above post. |
|
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 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 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. |
|
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. |
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.
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.
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:
And the real biggie:
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. |
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.
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.
That is fine, no worries. So all the ODE solvers you're writing would go into DiffEqODE.jl (or DAE), right?
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
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.
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?
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 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: Thus a integrator needs to provide the constructor 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 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.
Same here ;-) |
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
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.
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
(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
endIt'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).
No,
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
Trying. DDEs/explicit Jacobian usage needs to get hammered out before those are stable.
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).
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. |
Yep, code always gets better with 20-20 hindsight.
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)
endand 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.
Sounds good.
This is just about naming, using
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.
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.
Yes, it would only make sense to standardize the common options to keep them consistent between solvers, e.g. |
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.
Yeah, just naming.
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. Is this design making sense? Have any suggestions?
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 Specialized options to solvers are not an issue if you just always 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. |
@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.
Yes this sounds right. To bikeshed names: I'd prefer
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,...)
endAn 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(...)
endI think something along these lines could work, but maybe too complicated and not general enough...
Yes, that is good and what I put in above example.
:-) your clock ticks at a different rate to mine. In a month's time is plenty plenty early! |
|
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 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 As for 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. |
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 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 As for the solution, I think it's very clear to keep the By keeping these separate, the documentation is also more clear. We can document the So I think it's very clear that
Let me give a specific example.
This solves three issues:
There's no need to have a |
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. FlexibleIt 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
endwhat you're doing is making the algorithm have a callback function function callback(t,y)
@show (t,y)
if t > 1
break
end
endwhich 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)
endand 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. SimpleNo. I agree with you that
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 @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,tand have it work (replacing the preamble with explicit code of basically saying 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 GeneralNot 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 MisleadingThe 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
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). PerformancePacking 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! ConclusionThere 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 |
|
Also, thanks guys for being very responsive. It helps to be able to knock this out and get back to coding! |
|
To make it a little more clear why using the iterator interface ends up being a user trap, consider the standard event handling example: the bouncing ball with dampening. Let the first value be the position, and the second value be the velocity. If someone wants to use the iterator interface, they might think the following code makes sense: for (t,y) in ODE.iterate(ode)
if y[1] < 0 # Check if the position is negative
y[2] = -a*y[2] # Reverse the velocity and dampen
end
endHowever, this code will run into many issue for some very fundamental reasons.
Another case is the simple ODE for (t,y) in ODE.iterate(ode)
if y < 0 # Check if the position is negative
break
end
endHowever, if you use any high order method, the timesteps will become large enough that there is a very good chance it will skip over the negative timepoints. This is why event handling functions check at intermediate interpolated points for events. Those are two very basic examples which show why event handling in diffeq solvers have evolved the way they have. This is why once you have the iterators, it's still not even a good idea for most users to directly use the iterator interface. This is why MATLAB and all other "standard solvers" have specific functionality for dealing with events: it is hard to get it correct. Yet handling of events is the one thing that the iterator interface seemingly gives you, yet it's clear it misleads users. The only other features are:
Again, in every case except for 1, there is an easier solution (and for 1, callbacks work just fine). I hope that helps explain the point. |
We never added any storage because of the iterators. We only actually implemented two solvers for the iterator interface so the sample size might not be significant but the temp storage in the solvers was necessary for them to work anyway.
Well, this is an upstream issue, over which we have no influence. In that case you simply don't publish an iterator interface but a top-level interface (with callbacks functionality) instead.
And vice versa, I don't see how one is more flexible than the other. You even say it yourself: "The iterator interface does make it easy to add more advanced features (for example, the callbacks)". I don't think your arguments are sufficient to claim that callbacks are "strictly more capable".
@ode_preamble I get the idea, this particular example is simpler then what PR49 proposes and it has a lot of appeal to non-Julia users. The more I look at this example the more uncomfortably convinced I feel. But then you start adding the macros for loopfooter and loopheader and it gets more and more complicated to clearly see the whole picture (and you have to add them for the callback to work). Even with the effort you put in keeping these macros as simple as possible (they don't even take any arguments!) they are still quite extensive (on the order of hundreds of lines). And if you want to do something less standard and fiddle with parts of the code you would have to redefine the whole macro. While I certainly see the appeal of simplicity, comparing a code that uses them with a code that doesn't is not a fair game. While I'm not saying that PR49 approach is simpler on the user side, I want to highlight that the PR49 makes the definition of the solver almost completely self-contained. You can have everything clearly laid out in a single file. You still have to define a type, a constructor and a few methods, which, compared to the level of simplicity you suggest, seems like a fancy programming. Maybe there is a way to get the most of both worlds? It should be possible to define macros generating the PR49-like solvers. Looking at the code in
I don't see how this is related. BVP (I guess you mean boundary value problem) and SDE are completely different animals, we are talking about an interface for deterministic IVPs (see the topic, "ODE/DAE solver interface"). The classical IVP solvers seem to be conform very well to this decomposition into single forward steps. The stochastic case is out of my area of expertise, but maybe ODE/DAE and SDEs are simply not similar enough to warrant a common interface. The same goes for BVPs, that you mention earlier. You've probably given it much more though then I did, so could you please elaborate on why it makes sense to fit all of these problems into a single interface? Would it be possible to first build a lower level interface for IVP and only unify it with the interfaces for SDE/BVP one step higher in the interface hierarchy? This way the more specific functionality of IVP solvers would be still accessible to users working only with this class of problems.
I agree that this example might be misleading but it only proves that this is a bad example. This could be fixed by a better documentation. You are also right about about users becoming confused and trying to (ab)use iterators for events, but then again, this might be fixed with better documentation and with first introducing the user to a top-level interface (similar to what you propose), and the iterator interface later. And while I'm at it, we had root finding implemented at some point we just removed it to simplify and refactor the code. We already have generic/algorithm-specific interpolations though, along with the dense output. This proves that you can build more complex functionality on top of the iterator interface without changing the integration algorithms. We were thinking on event functions, which might be more challenging but still possible.
This is a valid concern. We are always passing a mutable type to iterator related functions and this might lead to some slowdown. That would require a benchmark comparison between ODE.jl (or other existing packages) and the optimized PR49. I used the form "the optimized PR49" to mark a difference between the "current PR49" which has never been completely optimized.
I believe that you are drawing these conclusions too hastily. I addressed most (all?) of your points and your most important concerns were about simplicity and performance. For the latter I haven't seen any actual numbers but I acknowledge that there might be some hard limitations. As for the simplicity, in DifferentialEquations.jl it is mostly achieved by a set of macros. In principle you could use similar macros to generate iterators, so I don't really see that as a fundamental issue (my distaste for macros aside). |
You seem to be confusing the iterator interface with event functions and root finding. |
|
IMO the iterator idiom seems just seems strange for ODEs. Chris' lays out the arguments better than I could. Is there a precedent or some sort of application you are thinking where this iterator interface is necessary/preffered? |
That's not addressing the statement. You already saved the current value of
This would be an issue if the iterator interface was the common interface. However, I already offered a solution which completely avoids this.
Then don't use macros? Those macros were just for my codes. With the
I'd argue the iterators are even less accessible. But still, this has nothing to do with the interface. You can still write a dispatch for
Because most people see these equations as extensions of ODEs. When you follow the logic of papers, people talk about an ODE system and ask "what happens if you add noise?". People simulate a gene regulatory network and then ask "what changes if you add the fact that the construction of proteins is controlled by a delayed signal?". In many problems ecological problems, people ask "what happens if I add a spatial element? Some kind of diffusion"? So for users, there is a natural transition from ODEs to SDEs/DDEs/PDEs. This is also a developer issue. Not only are delay equations related to ODEs, but the standard methods for solving DDEs is to use ODEs with event handling (for discontinuity propagation). Solving PDEs is essentially making and solving ODEs, and SDEs can be solved by solving PDEs or using MCMC methods. So for developers, there is a natural translation between the problems which is commonly exploited. The
Why not use the simpler interface for both?
What extra functionality are you gaining from the iterator interface? If you're not supposed to be putting code there because naively doing so handles many cases incorrectly, what specific functionality does it give you (that isn't done with callbacks)?
No, you didn't address the point. What you talked about instead was this:
Sure, you don't find my code accessible and I know it's not the most accessible because while my loops "are just the math", I attempt every hack under the sun to have both good performance and a huge array of features. I don't find your code accessible because the iterator interface "doesn't look like the math" for lots of its internal logic. Clearly, the most accessible code is ODE.jl which is why it still has a purpose. But that's not even what we're talking about. The discussion is: What should a common solver interface look like?I proposed the following:
where users define a problem, and The other possible proposal is an iterator interface. My argument against it is that:
That doesn't mean that you can't write your solvers using an iterator interface, and have it plug into the common interface. Rather, you can dispatch |
If that's not what it's for, then what does it do for a user? Is the manual supposed to say: for (t,y) in ODE.iterate(ode)
# Don't put code here. It won't do what you want it to do
end? The pro that you seem to suggest is summed up here:
The pro is: you like to code a few IVPs in an iterator format. But that's personal (and not something I agree with). Should it be imposed as the common interface for all DiffEq solvers? Should everyone have to code in an iterator style in order to be a part of the common interface? I don't see how you get there. What is clear is that, with the And why? What did would we actually gain by having an iterator common interface? It doesn't make code faster. It doesn't allow for more features. It doesn't support more equations. It doesn't support more implementations styles, modes, or languages. And it's not always simpler to code. Again, what about it makes it a compelling must-have common interface? |
|
Again, let me be clear, I do not have any issue at all with having a package in the JuliaDiffEq ecosystem which internally uses iterators. The "Package-Oriented Development, Distributed Maintenance" policy is about letting there be a plurality of implementation styles to encourage more developers to have their package or their method join the ecosystem, while helping them maintain the package for use in a common interface. We will want to show users a common interface: "look! Julia gives you fast and feature-filled methods for solving ODEs/SDEs/DAEs/DDEs/PDEs!", and this common interface will allow us to easily test between different implementations (very useful for research!) and recommend/default to different algorithms/packages in different situations at no cost. The The iterator interface does not accommodate all of this well, which is why I am against it as the common interface. That doesn't mean that you can't have an iterator interface for your package, and plug into the common interface by defining and solving an iterator. And on your package's README / extra documentation, you can detail this. |
|
I'm sold on the
What interface do you plan to be the standard for the callback argument? Is it the one defined by the macros in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve/ode_callbacks.jl? Is this going to be a part of the common interface? |
I did mention earlier that I don't know. It probably got lost in everything else:
The issue, since a callback is essentially injecting code, I don't see an easy way to make there be one function definition which works for all implementations: it needs to be tailored to the algorithm it's dealing with in order to get all of the functionality (the callbacks as implemented are very powerful, and can even do things like resize the internal caches, which clearly requires some specificity to the internal algorithm). One cheap solution is to have a well-defined DSL (macro for callback, event detection, and event apply), and let the different implementations add different macros which have the same user-interface but make different functions (with similar functionality). So for example, the current But it lends itself to a good long term solution. One those are all implemented, you can make it so that way writing a callback using the DSL: callback = @ode_callback begin
@ode_event event_f apply_event!
endwill generate specific codes for each algorithm. This is done by:
This idea is designed to match JuMP which uses a DSL for all of its user definitions to handle the problem that its solvers are very different backends with very different requirements on the definitions. JuMP has handled this very well, and so I think we can do the same. We also have a leg up because the implementation for our case is very easy: each package/implementation implements the callback DSL on its stuff, and one common find/replace macro will work (this common one won't need to be more than like 100 lines, and could be in DiffEqBase). |
|
I pruned down to make this common interface proposal be things which just about every implementation uses (or will be using in the near future). We can bikeshed names if you want. Problem Typetype ODEProblem{uType,uEltype,tType,fType<:Base.Callable} <: AbstractODEProblem
f::fType
u0::uType
analytic::Base.Callable
knownanalytic::Bool
numvars::Int
isinplace::Bool
tspan::Vector{tType}
endMost of this should make sense. Note that the common constructor is: ODEProblem(f::Function,u0,tspan=[0,1];analytic=nothing)and has two things to be aware of: Other Problem TypesOther problem types look very similar. It makes obvious small changes as necessary (for example, two functions for SDEs, the function in DAEs has a different signature and needs du0 as well, etc). I'm not sure on how to handle "extra functions". To you this probably just means Jacobians, but since we have the power of symbolic calculations (see ParameterizedFunctions), we can do things like symbolically pre-calculate the inverse of I suspect most people will let Common ODE OptionsHere's what I got in categories: dt
adaptive
gamma
abstol
reltol
qmax
qmin
dtmax
dtmin
qoldinit
beta1
beta2
internalnormThe standard adaptive timestepping options. dense
save_timeseries
timeseries_steps
saveat
calckThe next set of algorithms control the saving behavior. Then there's timeseries_errors
dense_errorsThis turns on/off the errors to calculate if an analytical solution is given. This is because in some cases it can be costly. progressbar
progress_stepsTurns on a progressbar, and designates how many steps to update it (since it can be costly). The Juno team got a new progressbar which has time estimates, I'm using it and adding it to Sundials and ODEInterface by default callbacks. I suggest everyone implement this because it's easy to do and good for users. Then there's some important stragglers: autodiffMany algorithms can choose to autodifferentiate or not. This should always be an option because it's not always possible (requires a bit of type-looseness), but always better. tTypeThe type for time. The type for the independent variables callbackSets a callback function. Common Options for Other ProblemsThe common options look pretty much the same. Anything missing for other problems? I think that this even covers everything in PDEs (that aren't algorithm options). ODE Algorithm Specification and OptionsAlgorithms will be specified by an immutable type in some type hierarchy. In DiffEqBase, there will be an immutable DP5
order
adaptive_order
beta1
beta2
endand have the constructor immutable CrankNicholson
order
linear_solver
endand have So the interface would normally look something like: solve(prob,DP5();kwargs...)to choose the DP5 algorithm, which could dispatch to my packages Passing OptionsI am backtracking on supporting the options in a type. Here's why:
As previously noted, the other way around is much rougher because types have a fixed number of spots.
Any objections? Default AlgorithmThe default algorithm will have immutable DefaultODEAlgorithm
alg_hints::Vector{Symbol}
endwith a constructor so that a user can do something like solve(prob,:stiff,:interpolant,:highprecision;kwargs...)and then DefaultAlgorithm(opargs...) will slurp up those into DocumentationThe DifferentialEquations docs will update to have the stuff about the common interface. For each equation, there will be a page on the problem type and the solver. The solver page will talk about the available algorithms and have a table for how kwargs are implemented. The table will be a checkbox yes/no if something is implement, and the defaults for things which are different amongst different packages ( SolutionsSolutions are TimingI should have most of this implemented on OrdinaryDiffEq by the end of the week, so let's bikeshed the names (I'll just regex change them as necessary). I hope to hammer this out and get a blog post together explaining the transition. |
|
I am sorry I don't have the time to read the whole thing so I will only comment on the problem type.
As for the in-place or non-in-place definitions, we decided to focus on the in-place approach (it leads to faster code) and have a simple mutable type with a single field, say a This leads me to another question. What is the minimal type interface (traits) that we require from a type for it to be supported. This should be fixed as a part of the interface. You are often mentioning supporting Unitful types, or DualNumbers and so on, so what is, in your view, the minimal type that an integrator should support? A minimal working type for PR49.jl is declared here https://github.com/JuliaODE/ODE.jl/blob/master/src/tests/minimal_types.jl. From the mathematical point of view And one more thing, I don't know what the standard Julia indentation is but we should pick one and stick with it. I, and I think @mauro3 too, am using 4 spaces indent while you are using 2 spaces. If there is going to be a common part of the code we should use a proper indentation. I looked it up in Julia Base and it seems that 4 spaces are standard. |
|
Just one more quick note If we add these options would it mean that dense output and saving options would have to be supported on a per solver basis? The old arguments aside, one benefit of the iterator interface is that you could easily add a generic dense output, save options and root finding algorithms on top of the integrators. The suggested approach to the interface is going to encourage implementing these on a per-integrator basis, which again adds boilerplate code. The other option would be to use some standard macros, as you do in your solvers, but I had an impression that we agreed to be able to avoid them. The same goes for other suggested options these would be cool to have but I don't image Jack or Jane implementing them. But still it would be nice to have these options for all the algorithms as they involve essentially the same functionality. That's why it helps to decompose the integrator into its loopy part and be able to do arbitrary stuff in between each iteration. We use iterators for that, you use macros, but both approaches address the same problems. |
I guess it isn't needed in the problem interface.
Same. That can be pruned (the default constructor indeed does this calculation).
It should be a common enough use case that every single algorithm that is developed uses it (I mean, you should have tests). It doesn't affect users at all if the default constructor makes it optional. As in the example: ODEProblem(f::Function,u0,tspan=[0,1];analytic=nothing)means that most people even if might never know you could give it an analytical solution or ever use it, but it'll still be there for tests.
First of all, it's fine if you don't want to explicitly support it. In the code where I wrapped your solvers, I just had if typeof(u) <: Number
u = [u]
end
if !isinplace && typeof(u)<:AbstractArray
f! = (t,u,du) -> (du[:] = prob.f(t,u))
else
f! = prob.f
endand it works. However, I want to support it because there are many interesting differential equations on non-vectors. Usually these arise in the cases of DDEs and SDEs (one famous one is the Black-Scholes equation SDE which is probably one of the most solved stochastic equations!). The performance difference between working with an immutable vs a vector of size 1 is pretty massive. It does require doubling up code so I understand if not everyone wants to fully support it, but I gave a quick way out above gives support at the cost of the not having the performance gains.
If it has |
|
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
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.) |
|
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. |
|
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:
not including the explicit parameter functions and parameter derivatives which aren't necessary for ODEs. I also plan on added explicit functions for |
|
ODEProblem has just one parameter for the function Internally, going down the |
|
But how do you check if a |
|
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 |
|
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]
endPutting 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). |
|
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 |
Yes, this is the simplest example of a 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]
endI guess call-overloading on types isn't too widespread yet (it's something that was introduced in v0.5), but it's the immutable HiddenType end
(f::HiddenType)(x,y) = x^2 + ywhere 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).
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), ...) = 2is 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 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
endTo solve |
|
Concerning #5 (comment): 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 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 |
Oh I see, that will work.
Yeah, but if
Yes. ParameterizedFunctions.jl just has macros to make a
Oh yeah oops. So yeah, put and ignore the method optional arg in the first. |
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 |
What do you mean "straight methods"? |
|
Say I implement a Alg: |
|
How do you do |
|
@ChrisRackauckas can you give me access rights to DiffEqBase.jl (and Roadmap whilst you're at it)? I have a branch to push. |
Done. |
|
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? |
|
Try it. |
|
Thanks, I made two PRs: SciML/DiffEqBase.jl#2 and SciML/DiffEqBase.jl#3, maybe discuss it further there. |
|
Alright, I moved my discussion over there. Is there any other high level topic to discuss that doesn't have its own thread? |
|
Closing this monster. I think we can start talking about specifics with @mauro3 's PRs. |
|
Thanks to all! I think this was very productive. |
|
Yeah, I learned quite a bit. I'm excited for traits now too! 👍 |
|
Thanks for a nice discussion, I'm very happy we agreed on some common ground. And, sorry for a limited participation. |
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 probablyDiffEqBase) and extend somesolvefunction 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:
solvesignature:solve(p::Problem, ::Type{Alg}; kwargs...)ProblemtypesSolutiontypesThe text was updated successfully, but these errors were encountered: