 # Simulating MJPs many times with different parameters as inner loop

Hi,

I’m currently coding up a particle filter (a.k.a. Sequential Monte Carlo, among other names) for approximate inference in Chemical Reaction Networks. My dynamics for the propagation of particles will be based on the SSA algorithm, as implemented within DifferentialEquations.jl.

The plan is to eventually use this filter within a Particle MCMC scheme (with some extra details), and so I need to be able to run the SSA algorithm many times, with many different initial conditions, and with different rate constants. For the moment, let’s just focus on varying the initial conditions.

My main concern is how I should construct the `DiscreteProblem, JumpProblem` objects which will be `solve`d in the main loop of the filter. Below is a toy version of the problem in which I am interested:

``````using Catalyst, DifferentialEquations, Distributions

lotvol_model = @reaction_network begin
c1, x1 --> 2x1
c2, x1 + x2 --> 2x2
c3, x2 --> 0
end c1 c2 c3

# Generate synthetic data
x0 = floor.(100 * rand(10))
y0 = [rand(Poisson(1)) for i in 1:10]

# Advance MJP from x by time Δt, with rates c

function MJPStep_LotVol(x, c, Δt)
MJP = DiscreteProblem(lotvol_model, x, (0.0, Δt), c)
jump_prob = JumpProblem(lotvol_model, MJP, Direct())
sol = solve(jump_prob, SSAStepper())
return sol(Δt)
end

# Importance Sampling with
# prior (independently across time) x, ... , x[T],
# evolution by MJP with rates c, and
# data y, ..., y[T]

function ImportanceSampling(x, y, c)
T = length(y)
Z = zeros(T)
W = zeros(T)

for t = 1:T
Z[t] = MJPStep_LotVol(x[t], c, 1)       # propagate particle using SSA
W[t] = logpdf(Poisson(Z[t]), y[t])      # weight particle by observation likelihood
end
end
``````

The situations in which I am interested correspond (roughly) to `T` in the above being very large.

My main worry is that the way in which I define `MJPStep_LotVol(x, c, Δt)` is probably quite inefficient - I’m guessing that forming the `DiscreteProblem, JumpProblem` at each step in the inner loop could be quite costly, or might be in some way unnecessary?

On the other hand, the problems which the function is trying to solve at each step are very similar (in this example, they only differ in initial conditions) and so it seems reasonable that there’s a simpler solution which isn’t so costly.

I cannot go deeply in your code to see if this makes sense, but maybe this thread is useful to give you insights:

1 Like

Thank you, this looks useful and relevant - I will have a read through it.

In the mean time, I’ve run (a debugged) version of the above code, and the cost of re-forming the sub-problems is not as bad as I had worried. There are probably still ways of speeding it up, though, so it should be worthwhile to read through the above.

It is unnecessary. I assume you just want to tell the stepper to take new random numbers?

New random numbers, possibly new initial conditions / new rate constants, but the same set of reactions throughout.

Then you just need to call `DiffEqJump.reset_aggregated_jumps!`. It’s not documented yet, but it should work just fine.

Would it be possible to elaborate on how this should be implemented?

e.g. if I do something like

``````DiffEqJump.reset_aggregated_jumps!
sol = solve(jump_prob, SSAStepper())
return sol(Δt)
``````

then I don’t ever pass the new rate constants `c` to the solver.

Do I modify `jump_prob` in-place somehow to hand it new initial conditions + rate constants?

I’m assuming that it is worth avoiding having to re-define

``````MJP = DiscreteProblem(lotvol_model, x, (0.0, Δt), c)
jump_prob = JumpProblem(lotvol_model, MJP, Direct())
``````

for each new `c`.

Sorry, just found this question again. You just add `DiffEqJump.reset_aggregated_jumps!(integrator)` inside of your callback.