ChrisRackauckas
commented
7 years ago Hey, thanks for sharing!
Reading the paper, the implementation of DiffEqJump is in your terms a version of the event-handling TJM. It doesn't group the rates to a total rate for the variable-rate jumps (it's no longer necessary when doing the integral via event-handling since the whole reason it was required before was due to changes in the rate function between jumps). For now, I like it because it keeps time unmodified which makes it easy to do things like mix other events in. Note that this implementation keeps fixed-rate equations separate, since they can be handled by a Gillespie algorithm separately and just place tstops values with a DiscreteCallback to be exact, so this sections off the only difficult things to handle as the variable-rate jumps.
But in the future, the CHV does look interesting. I see that you remove the need for event handling by making the "current time" a variable tau? That's interesting. And the integration before a jump is on an interval of length S_n (exponentially distributed), and so the integration of tau handles what the time to the jump will be and the other part solves the ODE and adjust the rates. That's a really nice way to solve the problem. In JuliaDiffEq terms, this would mean that it can be solved just by placing a tstop at S_n and use a DiscreteCallback at each end. Handling the time-change would be slightly difficult though, but it would just involve a closure over the user's function (and event functions) to automatically do the change back to "real-time" to make it invisible to the user.
I'll definitely have to implement that sometime in the future. There are a bunch of tradeoffs here. Increasing the stiffness or number of stiff components in a system will surely increase the error itself, if using a stiff solver (they are generally lower order and accumulate more error). They are also slower than nonstiff solvers. But nonstiff solvers on a stiff equation will be slow as well, but with very minimal error. Either way, adding stiffness will slow things down, so it's hard to know if this has a real speed benefit once you tune the event handling. I guess the only way to know would be to test it.
I think the true gains would only come from being "smart" with stiffness: not only stiffness switching, but sectioning out some parts as stiff and other parts as nonstiff. We can't do that yet, but we're almost there. We can do switching like LSODA (or you can just use LSODA on the common interface of course) using the CompositeAlgorithm, but you'd have to come up with the switching logic yourself until there's a default.
One thing I would check is whether MATLAB's event handling had any issues "missing events". After doing a ton of work on it, I am quite convinced that event handling cannot rely on defaults, and I know that both MATLAB and Sundials' event handling (Sundials' root finding) don't let you tweak all of the defaults for interpolation spacing and the like, and I know a bunch of examples of how to make the default values fail. So I think that if the full power of event handling is tweaked properly, then it should be able to handle those difficult problems. In DifferentialEquations.jl terms that would be changing the interp_points, but yes this does increase the computational cost quite a bit and requires a bit of user intervention and understanding (which is probably why it's usually not part of the API), so while it would likely handle it, avoiding event handling is a better option in the long term.
This stuff all carries over to SDEs as long as the time adaptivity / pullback is handled appropriately.
In general, since it's all in one cohesive API that works together, any advance in one area (stiff solvers, stiff switching, IMEX, adaptivity, event handling, the algorithms for encoding the jumps, etc.) seems to help all of the other areas. Hell, one thing I would like to point out is that Rosenbrock methods were generally not seen as competitive with BDF/Implicit RK methods, but with automatic and symbolic differentiation (instead of numerical differentiation) I see that in the majority of cases (i.e. not super large systems of ODEs derived from PDEs) that is not true: Rosenbrock methods seem to do well or better than other methods if you have the full scientific stack to make differentiation and nonlinear solving "better".
So yeah, I'll give CHV a try. Of course, I think a full "what does best" will depend on many factors, and it'll take a bunch of benchmarking. But in the end, my goal with JuliaDiffEq is to make everything an option and make it easy to switch between options and benchmark them. Only after it's all done will we know what is really best in each case. And I'm sure every algorithm has a case where it shines brightest.
CHV requires to have good stiff solvers and there were not available in JuliaDiffEq at the time of writing the package.
Well, I can see from the contribution timeline that the bigger issue is that JuliaDiffEq/DifferentialEquations.jl didn't exist when the bulk of the package was written!
Hi,
This is not really an issue (so you can change the tag maybe), but I have a package PDMP.jl that does similar things for jump processes without diffusion. Many exact algorithms are implemented like rejection algorithms but also a trick, called CHV, which I discovered and described in Arxiv.
CHV requires to have good stiff solvers and there were not available in JuliaDiffEq at the time of writing the package.
Best regards
Romain