Simplify out explicit integrals from continuous states

[ChrisRackauckas]

In the case where u' = f(t)dt, you can simply solve it via an integral, u = Int(f(t)dt).

In the case where you have a system of equations [u1',u2'] = [f1(u1), f2(u1)], you can remove u2 from the state and solve u1 = f1(u1) and then tag along u2 = Int(f2(u1)dt). This is what's enabled by the IntegratingCallback https://docs.sciml.ai/DiffEqCallbacks/stable/integrating/ and is a trick heavily used in things like the adjoint implementation. This can be very important for stiff systems since this can reduce the overall size of the Jacobian and its factorizations. Also, the integrated system can use Gauss-Legendre points, and so it can actually achieve very good convergence as long as the interpolation in time for the system is sufficiently good.

[isaacsas]

If we modeled the variable rate jumps symbolically and MTK did this optimization, we'd actually be quite a bit faster in many cases. But we'd have to move the rootfinding to the callback as well. But, because the interpolation of the ODE solution is a polynomial, that rootfinding would be a polynomial rootfind on the integral of the interpolating polynomial, and so we could use HomotopyContinuation to guarantee that we find all roots w.r.t. t. That would make it much more robust and quite a bit faster.

It's a good chunk of work but that is probably what we really need to make the hybrid modeling shine.

The nice part of doing it explicitly through MTK though is that this optimization can apply to many systems. A lot of systems people write down can simplify out part of it, so you'd get some good speedups in lots of places.

ExtendedJumpArray only had its weird behavior because SII didn't exist. Without SII and symbolic indexing, changing someone's u0 behind the scenes as part of jump "magic" can screw up the user's ability to understand the array, so we faked it as "an array that didn't change". In the context of MTK, we don't really need to do that though because there's already interfaces for handling the fact that the state can change and the user needs to respect that.

So with all of the SII things, and then also the inline SCC nonlinear solving stuff (which is an even bigger change to most Catalyst integration, especially with homotopy integration), we'll really need to reconsider the variable rate jump implementation since I would not be surprised to see it get 100x faster through connecting with the other optimizations going on.

Copying this from the other issue since this seems like a more appropriate place to discuss.

One thing we might want to consider is that propensities are non-negative, so the integrals should always be monotonic over time (except when they are reset because of an event firing). That should be exploitable by the integration scheme and the rootfind.

[isaacsas]

Another thing I was thinking is that we currently test each reaction separately for when it fires (essentially a first reaction method approach), but one could also presumably build a direct method equivalent of integrating and root finding to find the next reaction time and then sampling the reaction that occurs at that time. I wonder if that could actually be faster for many systems.

[ChrisRackauckas]

How would you sample though? You could single rootfind but you'd still have to integrate all of them.

[isaacsas]

I haven't worked with PDMPs much, so perhaps I'm missing something, but why do you need to integrate them individually instead of integrating their sum? And given the next jump time, wouldn't the jump that occurs just be a categorical sample using the values of the propensities/intensities at the jump time (not their integrals)?

[ChrisRackauckas]

Integrating the sum only can be used for the rootfinding. That would improve the rootfinding because then you'd have a single callback. But

wouldn't the jump that occurs just be a categorical sample using the values of the propensities/intensities at the jump time (not their integrals)?

The instantaneous rate doesn't really mean anything. Think about what happens if you have a rate that is zero, but in the epsilon time before a jump, that rate spikes to almost infinity. Should that one have an almost 100% chance of being selected? No, because the rate was zero for essentially the whole time the probability that jump fired is practically zero. So the probability of that one having fired is based on the proportional of the integrals, not the instantaneous rates. So you still have to calculate integrals for each one. This can be done via a Gauss formula though, so it's still better than the current form, but requires a modified integratingsumcallback. It would be a good student project that could get good results pretty quickly.

[isaacsas]

But that reaction would never be selected as its instantaneous rate is zero at the sampled reaction time.

I looked around and Anderson describes the method I mention above in Section V of https://arxiv.org/pdf/0708.0370

edit: and if you mean that in the small time before the next jump it transitions to very large and stays large, then it should be very likely to be the sampled reaction at the sampled time — it’s large values right before that time presumably determined why the next reaction time was at that point and not earlier.

[isaacsas]

This seems to use similar sampling: https://minds.wisconsin.edu/bitstream/handle/1793/10884/file_1.pdf?sequence=1&isAllowed=y

[ChrisRackauckas]

It's surprising that would work? I guess it must balance out because the probability you have a wild shift in the propensity at the reaction time must average out. Nobody actually is proving that though in the papers? But yeah if that is true, then that's a much simpler implementation numerically.

[isaacsas]

The way I was thinking of it the joint jump type ($k$) and time increment ($\tau$) density given the last event was at $t$ should be:

$$\rho(\tau,k | t) = a_k(t + \tau) e^{-\sum_k \int_t^{t + \tau} a_k(s) \, ds}$$

Where $a_k(t)$ is the propensity / intensity of the $k$’th jump type and the exponential is the probability no jump occurs until $t + \tau$. The time increment satisfies

$$\rho(\tau | t) = \left( \sum_{k} a_k(t + \tau) \right) e^{-\sum_k \int_t^{t + \tau} a_k(s) \, ds}$$

and the next jump type distribution is therefore

$$\rho(k | \tau, t) = \frac{a_k(t + \tau)}{\sum_j a_j(t + \tau)}$$

Maybe you are thinking about the probability the next jump is of type $k$, and not the probability it is type $k$ given the time it occurs.

[ChrisRackauckas]

Maybe you are thinking about the probability the next jump is of type k

indeed, that would need the integral. It's not clear to me though why you can just drop the integral when you condition on t and tau though.

[isaacsas]

If the conditional involved an integral then couldn’t one sample a jump for which the propensity is currently zero (but wasn’t zero in the recent past)?