# Simulating MJPs many times with different parameters as inner loop

**URL:** <https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [October 26, 2020, 5:37pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061 "2020-10-26T17:37:33Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 26, 2020, 5:37pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/1 "2020-10-26T17:37:33Z")

</div>

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:

```julia
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[1], ... , x[T], 
# evolution by MJP with rates c, and
# data y[1], ..., 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.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 26, 2020, 6:25pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/2 "2020-10-26T18:25:51Z")

</div>

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

> [@Define data-dependent function, various ways](https://discourse.julialang.org/t/define-data-dependent-function-various-ways/48988):
>
> This is a topic that is frequently raised here (I have asked this before). The problem is that of defining a function that will be input as a parameter for a solver, but the function depends on data while the function call inside the solver does not have a field to accept the data. Usually, the first try of a new programmer is to set the data on global variables, which is very detrimental to performance. Soon one finds out that if one declares the data as a constant global the performance is goo…

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 26, 2020, 8:20pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/3 "2020-10-26T20:20:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 26, 2020, 11:00pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/4 "2020-10-26T23:00:42Z")

</div>

> [@Sam\_Power](#):
>
> 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?

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

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [October 26, 2020, 11:14pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/5 "2020-10-26T23:14:15Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 31, 2020, 7:15am UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/6 "2020-10-31T07:15:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [November 18, 2020, 4:13pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/7 "2020-11-18T16:13:41Z")

</div>

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

e.g. if I do something like

```julia
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

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

```

for each new `c`.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 22, 2020, 4:21pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/8 "2020-11-22T16:21:09Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [November 27, 2020, 11:13am UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/9 "2020-11-27T11:13:52Z")

</div>

No problem, and thank you. I apologise, but I’m still not totally following. I’ll try to simplify my example, and be a bit more specific.

I’ll say that I have a few different initial conditions from which I want to solve the MJP, and call them `x[1], x[2], x[3]`.

```julia
using Catalyst, DifferentialEquations, Distributions

# specify parametrised MJP system
birthdeath_model = @reaction_network begin
    c1, x1 --> 2x1
    c2, x1 -> 0
end c1 c2

X_Initial = [10, 40, 70]
C_Global = [1, 2]

X_Final = zeros(3)

for i = 1:3
   MJP = DiscreteProblem(birthdeath_model, X_Initial[i], (0.0, 1.0), C_global)
   jump_prob = JumpProblem(birthdeath_model, MJP, Direct())
   X_Final[i] = solve(jump_prob, SSAStepper())

```

Where in the `for` loop does the `DiffEqJump.reset_aggregated_jumps!` call go?

I can sort of see how if it were the same `JumpProblem` for each `i`, then I could do something like

```julia
for i = 1:3
   DiffEqJump.reset_aggregated_jumps!(SSAStepper())
   X_Final[i] = solve(jump_prob, SSAStepper())

```

- would this be correct?

For my problem, I need to change the parameters of the problem at each step as well, and I’m assuming that `reset_aggregated_jumps` doesn’t handle this issue (?). So if I’m going to get rid of the lines of code which start with `MJP =`, `jump_prob = `, then I will need to specify the new initial conditions somewhere in the loop - otherwise I will end up re-solving the original problem with the same initial conditions. Maybe I’m missing something simple.

Thanks again for your help and patience.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 28, 2020, 9:46am UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/10 "2020-11-28T09:46:37Z")

</div>

> [@Sam\_Power](#):
>
> For my problem, I need to change the parameters of the problem at each step as well, and I’m assuming that `reset_aggregated_jumps` doesn’t handle this issue (?). So if I’m going to get rid of the lines of code which start with `MJP =` , `jump_prob = ` , then I will need to specify the new initial conditions somewhere in the loop - otherwise I will end up re-solving the original problem with the same initial conditions. Maybe I’m missing something simple.

It does handle it. Whether nothing changes or everything changes, it’s correct because of the no-memory property of the exponential distribution. But if anything changes you have to have it of course.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [November 28, 2020, 5:16pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/11 "2020-11-28T17:16:43Z")

</div>

@Sam_Power For your example, a `MassActionJump` should be getting created. Changing the parameter vector is unfortunately not currently propagated to the `MassActionJump`, which internally stores the rates for each jump, unless you create an entirely new `JumpProblem`. That said, it should be possible to get the `MassActionJump` within an existing `JumpProblem` and then change entries in the `scaled_rates` field directly. After that one would call `reset_aggregated_jumps!` as @ChrisRackauckas mentions.

I’ll try to put an example together of doing this for you in the next day or two, and verify it does indeed propagate the changed rates ok. (Sorry that I had missed your initial question. Feel free to ping me on Slack through the DiffEq-Bridged or SciML channels with future questions, or post to the DiffEqJump issues on Github – sometimes it is easy to miss stuff on Discourse.)

More generally we probably need some API functions that make this an easier process for users.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [November 28, 2020, 9:11pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/12 "2020-11-28T21:11:28Z")

</div>

OK, so on a _quick_ glance, if you only need to change parameters between solves this seems to work (but note the caveats I mention at the end):

```julia
using Catalyst, DiffEqBase, DiffEqJump
 rn = @reaction_network begin
           α, S + I --> 2I
           β, I --> R
         end α β
p = [.1/1000, .01]
tspan = (0.0,250.0)
u0 = [999,1,0] 
dprob = DiscreteProblem(rn,u0,tspan,p)
jprob = JumpProblem(rs, dprob, Direct())
sol = solve(jprob, SSAStepper())

# change the rate of the first reaction
jprob.massaction_jump.scaled_rates[1] = .001
sol = solve(jprob, SSAStepper())

```

I think this should work ok for changes between simulations, especially for `Direct`, but _please_ test with some examples of your own before trusting! I don’t think `reset_aggregated_jumps!` is actually needed here since the jump data structures are reinitialized with each call to solve, but again please double check on your own examples and let us know if you have any issues.

Now the caveats: if any of your reactions are not mass action, i.e. the rate is time-dependent or something like

```julia
k1*S, S+I -> R

```

where the rate expression involves a species, then you need to figure out the internal reaction order. In the generated `JumpSystem`, we internally order `MassActionJump`s in the order they appear within the `reaction_network` macro, so that will tell you the ordering of rate constants in the `scaled_rates` field. If you want to update the rate for a reaction involving a time-dependent component or as above, i.e. a reaction that is internally stored as a `ConstantRateJump` or `VariableRateJump`, you need to update the `p` vector and call `reset_aggregated_jumps!` as Chris mentioned.

The other thing to be aware of is that for a reaction like

```julia
k1, A+A -> B

```

the corresponding entry in `scaled_rates` will be `k1/2`, so if you update such an entry you should make sure to divide by two to preserve the correct scaling. (See [Model Simulation · Catalyst.jl](https://catalyst.sciml.ai/dev/tutorials/models/#Reaction-rate-laws-used-in-simulations))

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [December 2, 2020, 10:39am UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/13 "2020-12-02T10:39:39Z")

</div>

Thanks for this - this is very close to what I’m looking for! I’ll try to get in touch via Slack going forward.

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [December 2, 2020, 1:37pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/14 "2020-12-02T13:37:18Z")

</div>

Posting here is fine too! If you drop a link in the Slack channels though that will help with noticing your post if it seems we missed it.

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [December 2, 2020, 3:24pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/15 "2020-12-02T15:24:23Z")

</div>

Great, that sounds good. I’ve had a look through the API, and I think that I’m a bit closer to a solution now. I follow that

`jprob.massaction_jump.scaled_rates[1] = .001`

is being used to change the parameters, and so I’m guessing that something broadly similar should work to change the initial conditions.

At the moment, I’ve tried a couple of things lke

`jprob.states = [499,1,0]`  
`jprob.prob.states = [499,1,0]`  
`jprob.prob.u0 = [499,1,0]`

but these don’t seem to be working quite right. It seems like `reinit` is the way to handle this for ODEs, but doesn’t apply to `jprob` since it’s not an integrator per se (?). Is there an analog for jump problems?

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [December 2, 2020, 3:30pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/16 "2020-12-02T15:30:55Z")

</div>

What about `jprob.prob.u0 .= [499,1,0]`?

---

<div class="post-metadata">

**Author:** ![Sam\_Power](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sam_power/32/18972_2.png) [@Sam\_Power](https://discourse.julialang.org/u/Sam_Power)\
**Post date:** [December 2, 2020, 3:45pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/17 "2020-12-02T15:45:05Z")

</div>

There we go - it seems to be working, thank you!

---

<div class="post-metadata">

**Author:** ![isaacsas](https://avatars.discourse-cdn.com/v4/letter/i/f6c823/32.png) [@isaacsas](https://discourse.julialang.org/u/isaacsas)\
**Post date:** [August 19, 2021, 7:50pm UTC](https://discourse.julialang.org/t/simulating-mjps-many-times-with-different-parameters-as-inner-loop/49061/18 "2021-08-19T19:50:25Z")

</div>

Just to followup on this; `remake` should now properly change the initial condition and parameters for `JumpProblem`s, and we’ve added examples for this to the docs.