Repeatedly creating JumpProblem using same system gradually increase creation runtime (requires non-constant rates)

[TorkelE]

The following code:

using Catalyst, JumpProcesses, BenchmarkTools
rn = @reaction_network begin
    d*X, X --> 0
end
dprob = DiscreteProblem(rn, [rand()], (0.0,10.0), [rand()])

@btime jprob = JumpProblem(rn, dprob, Direct())

gives

  114.131 ms (95237 allocations: 6.32 MiB)

One would expect the time for repeats of the same call. However,

@btime jprob = JumpProblem(rn, dprob, Direct())
@btime jprob = JumpProblem(rn, dprob, Direct())
@btime jprob = JumpProblem(rn, dprob, Direct())
@btime jprob = JumpProblem(rn, dprob, Direct())
@btime jprob = JumpProblem(rn, dprob, Direct())
@btime jprob = JumpProblem(rn, dprob, Direct())

gives:

  121.393 ms (95237 allocations: 6.32 MiB)
  125.955 ms (95237 allocations: 6.32 MiB)
  130.440 ms (95237 allocations: 6.32 MiB)
  136.026 ms (95237 allocations: 6.32 MiB)
  140.068 ms (95237 allocations: 6.32 MiB)
  145.800 ms (95237 allocations: 6.32 MiB)

And it just seems to continue rising. One should note that rerunning

rn = @reaction_network begin
    d*X, X --> 0
end
dprob = DiscreteProblem(rn, [rand()], (0.0,10.0), [rand()])

does not reset the time. Indeed, even running

rn2 = @reaction_network begin
    d*X, X --> 0
end
dprob = DiscreteProblem(rn2, [rand()], (0.0,10.0), [rand()])
@btime jprob = JumpProblem(rn2, dprob, Direct())

I still see the same increased runtimes.

It should be noted that if instead have

rn = @reaction_network begin
    d, 2X --> X
end

then everything is fine, suggesting that if you only have mass action jumps this problem does not materialise.

[isaacsas]

Very weird, especially given that the allocations are staying constant.

Note you should interpolate in calling JumpProblem, i.e.

@btime JumpProblem($rn, $dprob, Direct())

but that doesn't seem to change this.

What led you to repeatedly create a JumpProblem from the same network like this?

[isaacsas]

I'm not seeing this when directly building a ConstantRateJump based JumpProblem, but am seeing the same issue if you convert rn to a JumpSystem and build a system using the latter. So it seems to be MTK related.

[isaacsas]

@ChrisRackauckas I think this is arising from https://github.com/SciML/ModelingToolkit.jl/blob/01aa16633be9fcb17a34c95e94c08ffa6d9d38a7/src/systems/jumps/jumpsystem.jl#L219 in MTK. Could this be due to repeatedly creating new RuntimeGeneratedFunctions?

[isaacsas]

This line https://github.com/SciML/ModelingToolkit.jl/blob/01aa16633be9fcb17a34c95e94c08ffa6d9d38a7/src/systems/jumps/jumpsystem.jl#L223 seems to be one of the culprits. I don't see any increase if I just call generate_affect_function, but I do see a consistent increase in running time if I wrap that in @RuntimeGeneratedFunction.

[isaacsas]

Transferred to MTK as this isn't a Catalyst issue from what I can see.

[TorkelE]

What led you to repeatedly create a JumpProblem from the same network like this? I had a system where I wanted to simulate it for a large number of parameter sets and check its behaviour (and for each parameter set for a couple of different initial conditions). So I just created a simulate_model function which took parameters and initial conditions etc., created the corresponding DiscreteProblem and then JumpProblem and simulated it.

I figured that the time to create the problem would be negligible compared to longer simulations. Now I rewrote it so that I just call remake on an initial JumpaProblem, so should be fine.

Thanks for having such a quick look at this!