Add discrete callbacks and tstops to SSAStepper about jumpprocesses.jl HOT 20 CLOSED

That's fine, it's not major for Pumas because pharmacometrics hasn't started using the these methods yet, but given their utility in biological models, my plan is to introduce it into the community with some nice examples.

from jumpprocesses.jl.

We'd have to time it. Tstops is one of the things it's avoiding so that it doesn't have to check a heap. At least a "fastest version of a stepper that is a bit less restrictive" can be made.

from jumpprocesses.jl.

Why do we need a heap? Is that to allow callbacks to dynamically add additional Tstops during a simulation? If we don't allow injecting additional Tstops, couldn't we just sort the initial Tstop vector and progress through it iteratively? Such a check would be one additional comparison per step, which for anything with a few jumps shouldn't be the dominant cost. I'd think the bigger issue is checking all the callbacks each iteration to see if they fired off (with the current callback setup).

Maybe we want to have a type of callback that is set to (one or more) fixed, user given times times, and a general termination condition call back? This would be a reduced set of functionality, but could cover many use cases I imagine.

from jumpprocesses.jl.

Actually, I’m not really sure why checking the top of a heap would be an issue either. That’s essentially what the NRM does every step, and certainly isn’t the dominant cost for it that makes it generally slower than Direct for simple systems.

from jumpprocesses.jl.

Actually, I’m not really sure why checking the top of a heap would be an issue either. That’s essentially what the NRM does every step, and certainly isn’t the dominant cost for it that makes it generally slower than Direct for simple systems.

It used to be a bigger deal, but it doesn't show up in as many profiles these days. It's worth checking. I would start from

https://github.com/JuliaDiffEq/SimpleDiffEq.jl/blob/master/src/functionmap.jl

and add in the pieces for SSA, and see how it compares to SSAStepper. If it's just as fast, then that could be the new backing to the SSAStepper, and then adding callbacks is simple.

from jumpprocesses.jl.

Ok, thanks for the pointer. I have to familiarize myself better with the integrator interface as I haven’t used it before. This is definitely something I’d like to take a look at once I get past those mid-November deadlines I mentioned.

from jumpprocesses.jl.

There's also another issue we really should address longer term -- we need an interface to update the SSAs internally so their state is consistent after a callback runs. With the exception of Direct, most have internal state that would need to be reset after changes to u.

from jumpprocesses.jl.

u_modified! is how it's handled in the ODE solvers. We might need a hook for that.

from jumpprocesses.jl.

If I understand how that works, we could provide an SSAIntegrator dispatch of reeval_internals_due_to_modification! for u_modified! to call into, which in turn then calls into some (as of yet undefined) interface to reinitialize an aggregator.

It might be good to get this setup more generally though, so that other SSAs than Direct could be used with existing integrators that already support callbacks and such. The lack of state reinitialization means pretty much all the non-Direct SSAs can only be used in SSAStepper currently.

from jumpprocesses.jl.

That's about it. We do need hooks into Direct anyways, since if you change u, the rates change. How to mathematically handle that correctly in general is 🤷‍♂ , but even if the starting rates are 0 you have an issue, since you need to change the next event from Inf to finite, which is the issue I am having in

https://github.com/PumasAI/Pumas.jl/issues/136
https://github.com/PumasAI/Pumas.jl/pull/472

We're going to need to brainstorm this more.

from jumpprocesses.jl.

You're completely right. I totally forgot the code flow is to execute a reaction and then calculate new rates and the next reaction time. (For some reason I was thinking we calculate the time and then execute the reaction... It's been a while since I actually looked at that code.)

The easiest thing to do would just be to reinitialize all rates, resort any associated data structures, and recalculate the next event time. If we know what components of u changed it would be reasonable to only update rates that depend on species that changed (but this would only be possible for dependency graph methods). This might actually be a place where RSSA shines. It would only require us to update rates that depend on changed u components that also moved outside their bracketing interval.

This will probably clobber performance if the callback happens often, but shouldn't be a big deal for ones that only occur a few times a given simulation.

from jumpprocesses.jl.

Correctness is what I'm worried about. If it's Inf, then you just replace it with the timing from the next rate. That's clear. But if it's not Inf, then you need to take new rates, but how do you discount the time that already occurred? If callbacks keep happening, then you can get a bad case where the chance of a reaction happening is essentially zero? Think of a periodic callback hitting and doing integrator.u[1] += 0 every small dt.

At least that's what I was thinking, but double check me on this though. Because the exponential distribution is memory-less, we can get away with a form of "chopping". What we instead need to calculate is the time to the next reaction given a total rate lambda conditional on an amount of time t already having happened. But that conditional distribution, because of the memory-less property, is just the same exponential distribution (!?), so we actually don't need to do anything special.

Another way to interpret this would be back to the definition of the Poisson distribution as having a constant time of occurring in any interval. Given that, the fact that you don't change how you're choosing the rate has to be the same with each ``integrator.u[1] += 0callback call in adt`, so then the chance of a reaction ever occurring is the same between every callback.

So, for awhile I was scared that if I do something too simple that I would ruin the statistics, but it seems like the simple thing is the correct thing. Please double check this though, and if it's right, then essentially what we need to do is re-initialize with every u_modified!. Implementing that might need a hook of some sorts.

from jumpprocesses.jl.

Let's take the (equivalent to Direct method) first reaction-method approach, and think of all individual molecules or pairs of molecules as having their own individual reaction times (instead of thinking about total rates among all molecules of each species). I would (naively) think that for the individual molecules that are unaffected by the callback, the times for their individual reactions would not change. So if we have A+B -> C and lose some A's and B's from a callback, the times for the remaining A and B pairs to individually react shouldn't change. By the Markovian property we can always resample their individual pair times to react, or equivalently switch to sampling from the total rate for some pair in the whole population to react. This is hand wavy though.

The right way to work this out is to look at the master equation or the Kurtz process representation with the callback dynamics added into it. Can a fixed callback time be represented as a counting process with a delta-function intensity centered at that time? (edit: don't think so). Maybe it makes more sense to use the more general Poisson random measure representation for the dynamics here.

edit: I have to include it as a discontinuous deterministic term in the equation.

from jumpprocesses.jl.

I guess you can view it as a PDMP too. If you consider the process Y(t) = (t, X(t)), then there are two ways to jump:

1. a stochastic one associated with some transition rates
2. a deterministic one when the deterministic flow (here (x,t)->t) hits a boundary.

Davis' book is a reference about this.

from jumpprocesses.jl.

So yeah, this is cool because it means that handling events is fundamentally "easy", so we just computationally need to handle it in some way that it knows to re-sample the rates when another callback fires. Instead of doing a pure callback set, I wonder if we can do some kind of wrapper to make it general. But at least SSAStepper could directly handle this in a simple manner.

from jumpprocesses.jl.

As long as we're planning this out, do we want to have a single reinit function that works if either u or a parameter has changed?

Am I right that we can rely on the uprev field in the integrator to have the old u value before the callback was called? In that case we can figure out which values changed and only update the relevant rates (for dependency graph methods).

from jumpprocesses.jl.

from jumpprocesses.jl.

Can that function assume the passed in integrator uprev value will have the value of u right before the callback was applied? If we don't want to reinitialize all propensity functions we would need to be able to determine which components of u had changed.

from jumpprocesses.jl.

Yes it can make that assumption

from jumpprocesses.jl.

OK, then this should be completely doable, and near trivial for some of the methods.

Once my November 8 deadline passes I'll have more time and put this near the top of my priority list. I hadn't realized this was something that was blocking any downstream users/packages... (I don't follow the Pumas issues.)

from jumpprocesses.jl.