# Memory usage while solving \`JumpProblem\` increases with run time, even with \`save\_positions=(false,false)\`

**URL:** <https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656>\
**Category:** New to Julia\
**Tags:** differentialequation, catalyst\
**Created:** [March 27, 2023, 11:47am UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656 "2023-03-27T11:47:31Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [March 27, 2023, 11:47am UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/1 "2023-03-27T11:47:31Z")

</div>

I am running into problems regarding memory allocation that depends on the actual runtime of the solver of a stochastic jump process (i.e. it depends on `T` that defines `tspan(0.0,T)`). This gives notable performance drops when multi-threading a `JumpProblem()` using `EnsembleProblem(args, EnsembleThreads())`.

To illustrate the problem, I have a MWE that I present below. It highlights the memory issue without defining an ensemble problem. Note that the real problem I am solving is more complex and there the memory issues actually lead to decreased performance when multi-threading and/or using distributed solutions (using `EnsembleThreads()` and/or `EnsembleDistributed()`).

( **All benchmarks below are subsequent calls to the functions. My code has all this in a single module, but I made it into a MWE for clarity.** )

## Minimal working example

As an illustration/MWE of my problem, let us consider a very simple 2-species Lotka-Volterra system with a carrying capacity K. This is just the standard 2-species generalized Lotka-Volterra model with unit interaction a\_{xy} = a\_{yx} = 1 and K\_x = K\_y = K.

\frac{dx}{dt} = \eta x \left(1 - \frac{x+y}{K}\right) \\ \frac{dy}{dt} = \eta y \left(1 - \frac{x+y}{K}\right)

For this problem, we expect simple logistic growth with rate \eta until x+y = K, after which the system is in a stable state.  
Let us define the function using `Catalyst.jl`, and solve the ODE directly:

```julia
function lotkavolterra(;T::Float64=100.0, nsaves::Int64=1)
    function growth(η,X,Y,K)
        return η*(1-(X+Y)/K)
    end
    
    lv_model = @reaction_network begin
        growth(η,X,Y,K), X --> 2X
        growth(η,Y,X,K), Y --> 2Y
    end

    p = (:η => 0.1, :K => 1000)
    u0 = [:X => 50, :Y => 100]
    tspan = (0.0,T)
    Δt = (tspan[end]-tspan[begin])/nsaves
    tsave = tspan[begin]:Δt:tspan[end]

    odeprob = ODEProblem(lv_model, u0, tspan, p)
    odesol = solve(odeprob, Tsit5(); saveat=tsave)
    return odesol
end

```

```julia
julia> sol = lotkavolterra(T=100.0, nsaves=100);
julia> plot(sol)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/5/d509bee7818be0ec2696ec064436d944b197b4d8.png)  
So far so good.

Additionally, we expect `Catalyst.jl` to build non-allocating functions (see also [this answer](https://discourse.julialang.org/t/improving-performance-of-solving-odesystem-catalyst-jl/83568/9)), so the total memory allocation should only depend on `nsaves`, not on `T`.  
(I am only leaving in the `Memory estimate`, as the time is at the moment not of interest.)

```julia
julia> using BenchmarkTools
julia> @benchmark lotkavolterra(T=100.0, nsaves=10)
BenchmarkTools.Trial: 10000 samples with 6 evaluations.
 Memory estimate: 6.91 KiB, allocs estimate: 59.
julia> @benchmark lotkavolterra(T=100.0, nsaves=100)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Memory estimate: 17.50 KiB, allocs estimate: 150
julia> @benchmark lotkavolterra(T=1000.0, nsaves=100)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Memory estimate: 17.50 KiB, allocs estimate: 150.

```

So far so good.  
Now we simulate the Lotka-Volterra system by defining a `JumpProblem`, and we are now only interested in the final state, thus:

```julia
function lotkavolterra(;T::Float64=100.0)
    function growth(η,X,Y,K)
        return η*(1-(X+Y)/K)
    end
    
    lv_model = @reaction_network begin
        growth(η,X,Y,K), X --> 2X
        growth(η,Y,X,K), Y --> 2Y
    end

    p = (:η => 0.1, :K => 1000)
    u0 = [:X => 50, :Y => 100]
    tspan = (0.0,T)
    Δt = (tspan[end]-tspan[begin])/nsaves
    tsave = tspan[begin]:Δt:tspan[end]

    discrete_prob = DiscreteProblem(lv_model, u0, tspan, p)
    jump_prob = JumpProblem(lv_model, discrete_prob, Direct())
    jsol = solve(jump_prob, SSAStepper())
    return jsol
end

```

First to double check:

```julia
>julia sol = lotkavolterra(T=100.0); plot(sol)

```

 ![image](https://global.discourse-cdn.com/julialang/original/3X/1/a/1a6002befd627e0e1c3648b227f59ffe7e118a10.png)  
Seems good.

## Memory usage increases as total time increases

Now the problem. What I expect is to see `Catalyst.jl` build non-allocating functions, similar to solving it with an ODE, but now I actually observe a dependence on the total run time `T`. Again omitting runtime information and only saving the final state by adding `save_positions=(false,false)` to the definition of `JumpProblem(...)`, i.e.

```julia
    jump_prob = JumpProblem(lv_model, discrete_prob, Direct(); save_position=(false,false))

```

```julia
julia> @benchmark lotkavolterra(T=100.0)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Memory estimate: 41.64 KiB, allocs estimate: 1703.
julia> @benchmark lotkavolterra(T=10.0)
BenchmarkTools.Trial: 10000 samples with 5 evaluations.
 Memory estimate: 8.78 KiB, allocs estimate: 301.

```

Why does the memory estimate increase?

## Investigations

Interestingly, when I change the growth rate function to an in-place version:

```julia
function growth!(η,X,Y,K)
    return η*(1-(X+Y)/K)
end

```

nothing really changes – as in, I get the same increase in memory and allocations.  
What is more interesting is that when the specify carrying capacity term is removed (as in [the examples of `JumpProcesses`](https://docs.sciml.ai/JumpProcesses/stable/tutorials/discrete_stochastic_example/)), the memory does **not** increase! Other solves, such as `RSSA()` also have the same problem.

### Question

What is happening here? Why does the memory usage increase when the runtime increases when I define a `JumpProblem()` from a `DiscreteProblem()`, even though only the last state is saved by setting `save_positions=(false,false)`. Does this have to do with the interpolation that is done internally? Is there any way to have terms, such as the carrying-capacity term in this example, without increasing the memory usage?  
I ask because my _real_ problem is very similar – basically an n-species Lotka-Volterra model with different carrying capacity terms, and recall I am constantly running into memory-related issues, especially when turning the `JumpProblem()` into an `EnsembleProblem()` for multi-threading and/or distributed computing later down the line.

Any help or tips would be greatly appreciated, thanks!

## Full code

> **Module**
>
> ```julia
> module TestJumpProcesses
> 
> export lotkavolterra
> 
> using JumpProcesses
> using Catalyst
> using DifferentialEquations
> using BenchmarkTools
> 
> function lotkavolterra(;T::Float64=100.0, nsaves::Int64=1)
> function growth(η,X,Y,K)
> return η*(1-(X+Y)/K)
> end
>     
> lv_model = @reaction_network begin
> growth(η,X,Y,K), X --> 2X
> growth(η,Y,X,K), Y --> 2Y
> end
> 
> p = (:η => 0.1, :K => 1000)
> u0 = [:X => 50, :Y => 100]
> tspan = (0.0,T)
> Δt = (tspan[end]-tspan[begin])/nsaves
> tsave = tspan[begin]:Δt:tspan[end]
> 
> @info "Benchmarking ODE problem..."
> odeprob = ODEProblem(lv_model, u0, tspan, p)
> odesol = solve(odeprob, Tsit5(); saveat=tsave)
> bode = @benchmark solve($odeprob, Tsit5(); saveat=$tsave)
> 
> @info "Benchmarking jump problem..."
> discrete_prob = DiscreteProblem(lv_model, u0, tspan, p)
> jump_prob = JumpProblem(lv_model, discrete_prob, Direct(); save_positions=(false,false))
> sol = solve(jump_prob, SSAStepper())
> bjump = @benchmark solve($jump_prob, SSAStepper())
> 
> return bode, bjump
> end
> 
> ```

---

<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:** [March 27, 2023, 12:28pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/2 "2023-03-27T12:28:01Z")

</div>

Open an issue on JumpProcesses.jl

---

<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:** [March 27, 2023, 1:00pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/3 "2023-03-27T13:00:07Z")

</div>

The memory use with `EnsembleProblem`s is probably because of it deepcopying the `JumpProblem` for each simulation (when using multi-threading):

> <https://github.com/SciML/JumpProcesses.jl/blob/f029c45049e425e249abe6f59b88bb1460f67d6c/src/SSA_stepper.jl#L135-L161>

You are benchmarking a lot of stuff beyond just calling `solve`. Do you see a memory use difference if you just benchmark `solve` using problems that only differ in their final time?

---

<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:** [March 27, 2023, 1:13pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/4 "2023-03-27T13:13:28Z")

</div>

Here is what I see just benchmarking `solve`:

```julia
julia> function lotkavolterra(;T::Float64=100.0)
           function growth(η,X,Y,K)
               return η*(1-(X+Y)/K)
           end
           
           lv_model = @reaction_network begin
               growth(η,X,Y,K), X --> 2X
               growth(η,Y,X,K), Y --> 2Y
           end

           p = (:η => 0.1, :K => 1000)
           u0 = [:X => 50, :Y => 100]
           tspan = (0.0,T)

           discrete_prob = DiscreteProblem(lv_model, u0, tspan, p)
           jump_prob = JumpProblem(lv_model, discrete_prob, Direct(), save_positions = (false,false))
           return jump_prob
       end
lotkavolterra (generic function with 1 method)

julia> jprob = lotkavolterra(T = 10.0)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 2
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 0

julia> @btime solve($jprob, $(SSAStepper()))
  8.750 μs (243 allocations: 7.42 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
  0.0
 10.0
u: 2-element Vector{Vector{Int64}}:
 [50, 100]
 [111, 235]

julia> jprob2 = lotkavolterra(T = 100.0)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 2
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 0

julia> @btime solve($jprob2, $(SSAStepper()))
  60.208 μs (1703 allocations: 41.64 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
   0.0
 100.0
u: 2-element Vector{Vector{Int64}}:
 [50, 100]
 [287, 712]

julia> jprob3 = lotkavolterra(T = 1000.0)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 2
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 0

julia> @btime solve($jprob3, $(SSAStepper()))
  60.167 μs (1711 allocations: 41.83 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
    0.0
 1000.0
u: 2-element Vector{Vector{Int64}}:
 [50, 100]
 [395, 605]

julia> jprob4 = lotkavolterra(T = 10000.0)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 2
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 0

julia> @btime solve($jprob4, $(SSAStepper()))
  60.000 μs (1711 allocations: 41.83 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
     0.0
 10000.0
u: 2-element Vector{Vector{Int64}}:
 [50, 100]
 [412, 588]

```

but note that the system will eventually stop evolving when the growth rate hits zero so I guess it isn’t surprising there is no further change.

Finally, here is a simpler example that does not “turn off” (i.e. reactions keep occurring for all time). I see no change in allocations even though as expected the running time increases:

```julia
julia> rn = @reaction_network begin
           @species A(t) = 0.0
           1.0, 0 --> A
           end
Model ##ReactionSystem#484
States (1):
  A(t) [defaults to 0.0]
Parameters (0):

julia> dprob = DiscreteProblem(rn, [], (0.0, 10.0))
DiscreteProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 10.0)
u0: 1-element Vector{Float64}:
 0.0

julia> jprob = JumpProblem(rn, dprob, Direct(); save_positions = (false,false))
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 0
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 1

julia> sol = solve(jprob, SSAStepper())
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
  0.0
 10.0
u: 2-element Vector{Vector{Float64}}:
 [0.0]
 [7.0]

julia> @btime solve($jprob, SSAStepper())
  768.817 ns (11 allocations: 1.94 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
  0.0
 10.0
u: 2-element Vector{Vector{Float64}}:
 [0.0]
 [7.0]

julia> dprob2 = DiscreteProblem(rn, [], (0.0, 100.0))
DiscreteProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 1-element Vector{Float64}:
 0.0

julia> jprob2 = JumpProblem(rn, dprob2, Direct(); save_positions = (false,false))
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 0
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 1

julia> sol = solve(jprob2, SSAStepper())
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
   0.0
 100.0
u: 2-element Vector{Vector{Float64}}:
 [0.0]
 [105.0]

julia> @btime solve($jprob2, SSAStepper())
  3.923 μs (11 allocations: 1.94 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
   0.0
 100.0
u: 2-element Vector{Vector{Float64}}:
 [0.0]
 [129.0]

```

I also see no memory use difference on this example if I change `1.0` to `1.0 + (1 - A/100)` which forces usage of a `ConstantRateJump`.

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [March 27, 2023, 2:51pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/5 "2023-03-27T14:51:40Z")

</div>

Yes sorry, I had included the [full code](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656#full-code-5) in order to avoid this, but it was easy to miss. In the full code I only benchmarked the `solve()` and got similar results to you. Note that in both my case and in yours memory usage increases from `T=10.0` versus `T=100.0`. Is this to be expected?

> Here is what I see just benchmarking `solve`

It is indeed as you say, probably because my MWE was too simple it stagnates at some point.  
Interestingly, the examples you provided (and some others I tried) do seem to ensure work…

For example

```julia
function loggrowth(;T::Float64=100.0)
    @parameters η, K
    @variables t
    @species A(t) B(t)
    rn = @reaction_network begin
        η*(1 - (A+B)/K), A --> 2A
        μ, B --> 0
    end
    p = (:η => 1, :μ => 1e-6, :K => 1e7)
    u0 = [:A => 1, :B => 1e5]
    tspan = (0.0, T)

    dprob = DiscreteProblem(rn, u0, tspan, p)
    jprob = JumpProblem(rn, dprob, Direct(); save_positions=(false,false))
    return jprob
end

```

```julia
julia> long_prob = loggrowth(T=1e7)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 1
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 1

julia> long_sol = solve(long_prob, SSAStepper())
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
 0.0
 1.0e7
u: 2-element Vector{Vector{Float64}}:
 [1.0, 100000.0]
 [9.999998e6, 2.0]

julia> @btime solve($long_prob, SSAStepper())
  157.063 ms (11 allocations: 2.09 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
 0.0
 1.0e7
u: 2-element Vector{Vector{Float64}}:
 [1.0, 100000.0]
 [9.999995e6, 5.0]
julia> short_prob = loggrowth(T=1e4)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 1
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 1

julia> short_sol = solve(short_prob, SSAStepper())
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
     0.0
 10000.0
u: 2-element Vector{Vector{Float64}}:
 [1.0, 100000.0]
 [9.900929e6, 99071.0]

julia> @btime solve($short_prob, SSAStepper())
  152.330 ms (11 allocations: 2.09 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
     0.0
 10000.0
u: 2-element Vector{Vector{Float64}}:
 [1.0, 100000.0]
 [9.901011e6, 98989.0]

```

### Question

So, what is different in the code that I provided earlier, and why does this one show memory increases? I will play around a bit more and come back to this, as the difference between the MWEs and my real code is very odd.

---

<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:** [March 27, 2023, 2:55pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/6 "2023-03-27T14:55:40Z")

</div>

Yeah, I have no idea yet why your original example shows allocation differences. It is indeed weird.

---

<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:** [March 27, 2023, 3:04pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/7 "2023-03-27T15:04:37Z")

</div>

Just so you are aware:

```julia
    @parameters η, K
    @variables t
    @species A(t) B(t)

```

in `loggrowth` is not doing anything in your latest example. What Catalyst version are you on? In the latest version the DSL would interpret any substrate/product as a species (`A` and `B`), and everything else as a parameter automatically. That is why your `loggrowth` function works ok. If you put those commands inside the `@reaction_network` macro then they would be used.

If you want the `@reaction_network` macro to see an external variable you have to interpolate it, like `$A`, when using it.

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [March 27, 2023, 3:07pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/8 "2023-03-27T15:07:32Z")

</div>

> Just so you are aware:

```julia
    @parameters η, K
    @variables t
    @species A(t) B(t)

```

> in `loggrowth` is not doing anything in your latest example. What Catalyst version are you on? In the latest version the DSL would interpret any substrate/product as a species (`A` and `B`), and everything else as a parameter automatically. That is why your `loggrowth` function works ok. If you put those commands inside the `@reaction_network` macro then they would be used.

I did not know this, I was just following the latest tutorials that the documentation `Catalyst.jl` provided (or perhaps it was `JumpProcesses.jl`) – but it appears this notation is indeed deprecated?

> What version of `Catalyst.jl`?

```julia
(examples) pkg> status Catalyst
Status `~/<path/to/project>/Project.toml`
  [479239e8] Catalyst v13.2.0

```

---

<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:** [March 27, 2023, 3:10pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/9 "2023-03-27T15:10:10Z")

</div>

Check out:

[https://docs.sciml.ai/Catalyst/stable/catalyst\_functionality/dsl\_description/#Explicit-specification-of-network-species-and-parameters](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/#Explicit-specification-of-network-species-and-parameters)

[https://docs.sciml.ai/Catalyst/stable/catalyst\_functionality/dsl\_description/#Specifying-alternative-time-variables-and/or-extra-independent-variables](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/#Specifying-alternative-time-variables-and/or-extra-independent-variables)

[https://docs.sciml.ai/Catalyst/stable/catalyst\_functionality/dsl\_description/#dsl\_description\_interpolation\_of\_variables](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/dsl_description/#dsl_description_interpolation_of_variables)

---

<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:** [March 27, 2023, 3:14pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/10 "2023-03-27T15:14:32Z")

</div>

Can you point me to the tutorial you are referencing? It sounds like something didn’t get updated (since you are on the latest version…).

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [March 27, 2023, 3:18pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/11 "2023-03-27T15:18:12Z")

</div>

I wish I could but it seems all of the links in my history now point to the updated version. Do you know when the docs was updated? I learned about `Catalyst.jl` about ~2 weeks ago and have not went back to check on any updated tutorials after following a few around that time. I do remember it being like this (perhaps it was in a blog post?). I will dig into my history a bit more and edit or reply if I find it.

But my guess it was this one, as I was generating reaction networks programmatically:  
[https://docs.sciml.ai/Catalyst/stable/catalyst\_functionality/programmatic\_CRN\_construction/](https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/programmatic_CRN_construction/)

---

<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:** [March 27, 2023, 3:21pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/12 "2023-03-27T15:21:34Z")

</div>

I don’t think anything has been updated in several weeks, and nothing in the newer notation for longer than that.

In any case, it sounds like you are now getting the correct versions. If you find something wrong somewhere just let us know.

I’ll try to play with your example with the extra allocations more, but if you make any more progress on distinguishing what features induce the allocations please do also let me know!

---

<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:** [March 27, 2023, 3:54pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/13 "2023-03-27T15:54:33Z")

</div>

@ChrisRackauckas could this be the issue? BenchmarkTools does report allocations in calling this function:

> <https://github.com/SciML/JumpProcesses.jl/blob/f029c45049e425e249abe6f59b88bb1460f67d6c/src/aggregators/ssajump.jl#L184>

I’m noticing the `getindex` here generates code that is a union type (since `p.affects!` is a tuple of functions). I’m not sure how to get around this though, unless we want to just drop non-`FunctionWrapper` versions of `Direct`. (There is also FunctionWranglers.jl, but it seems unmaintained…)

Switching to `DirectFW` reduces this but doesn’t eliminate the issue.

---

<div class="post-metadata">

**Author:** ![johannesnauta](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johannesnauta/32/47434_2.png) [@johannesnauta](https://discourse.julialang.org/u/johannesnauta)\
**Post date:** [March 27, 2023, 4:01pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/14 "2023-03-27T16:01:20Z")

</div>

Not sure if it’ll help, but I have found a problem that highlights the issue very clearly (if I am not mistaken…). It concerns a problem that is very close (basically a simplification) to the one presented in equation (1) of [this paper](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4453012/). That is, a problem of plasmid maintenance in bacteria:

\begin{align} \frac{dF}{dt} &= \eta F \left( 1 - \frac{F+P}{K} \right) - \gamma F P + \rho P - \mu F, \\ \frac{dP}{dt} &= \eta P \left( 1 - \frac{F+P}{K} \right) + \gamma F P - \rho P - \mu P \end{align}

where \eta the intrinsic growth rates for the plasmid-free (F) and plasmid-carrying (P) populations respectively. \gamma is the plasmid infection rate, \rho the recovery rate (through segregation error) and \mu the background mortality rate.

The code for this system using `Catalyst.jl`:

```julia
function catalyst_plasmid(;T::Float64=250.0)    
    function rate(η,X,Y,K)
        return η*(1-(X+Y)/K)
    end
    
    fp_model = @reaction_network begin
        rate(η,F,P,K), F --> 2F # Background reproduction F
        rate(η,P,F,K), P --> 2P # Background reproduction P !NOTE: order of F & P are switched
        μ, F --> 0 # Background mortality F
        μ, P --> 0 # Background mortality P
        γ, F + P --> 2P # Infection (conjugation)
        ρ, P --> F # Recovery (segregation error)
    end    
    p = (:η => 1.0, :μ => 0.1, :γ => 1e-5, :ρ => 0.01, :K => 1e4)
    u0 = [:F => 100, :P => 10]
    tspan = (0.0, T)

    dprob = DiscreteProblem(fp_model, u0, tspan, p)
    jump_prob = JumpProblem(fp_model, dprob, Direct(); save_positions=(false,false))
    return jump_prob
end

```

and the benchmarks

```julia
julia> using Catalyst, DifferentialEquations

julia> stepper = SSAStepper()
SSAStepper()

julia> short_prob = catalyst_plasmid(T=1e3);

julia> short_sol = solve(short_prob, stepper); 

julia> @btime solve($short_prob, $stepper)
  87.697 ms (1805709 allocations: 41.33 MiB)

julia> long_prob = catalyst_plasmid(T=1e5); # 100 times as long

julia> long_sol = solve(long_prob, stepper);

julia> @btime solve($long_prob, $stepper)
  10.803 s (180003301 allocations: 4.02 GiB)

```

So increasing runtime by 100 also seems to increase memory usage by a factor of 100. There must be wrong somewhere – let’s hope it is in my implementation!  
Sadly, I cannot remove the carrying-capacity term as this allows for infinite growth and, rightfully, crashes the REPL.

---

<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:** [March 27, 2023, 4:09pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/15 "2023-03-27T16:09:38Z")

</div>

> [@isaacsas](#):
>
> I’m noticing the `getindex` here generates code that is a union type (since `p.affects!` is a tuple of functions). I’m not sure how to get around this though, unless we want to just drop non-`FunctionWrapper` versions of `Direct`. (There is also FunctionWranglers.jl, but it seems unmaintained…)

For the non-functionwrapper versions we should codegen to manual split this to runtime branching.

```julia
if idx === 1
  @inbounds p.affects![1](integrator)
elseif idx === 2
  @inbounds p.affects![2](integrator)
...

```

would make it go away.

---

<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:** [March 27, 2023, 4:11pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/16 "2023-03-27T16:11:18Z")

</div>

Yeah, though even using `FunctionWrappers` I’m still seeing small allocations.

Could somehow the generated function from Symbolics be allocating a small amount here?

---

<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:** [March 27, 2023, 4:18pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/17 "2023-03-27T16:18:13Z")

</div>

First time called compilation? Not all affects will always be called.

---

<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:** [March 27, 2023, 4:23pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/18 "2023-03-27T16:23:27Z")

</div>

Yeah, when calling via the `FunctionWrapper` to call the `affect!` function in `DirectFW` I’m still getting a reported single 32 byte allocation. It occurs consistently, not just on the first call.

---

<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:** [March 27, 2023, 4:30pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/19 "2023-03-27T16:30:21Z")

</div>

@ChrisRackauckas can confirm this pure JumpProcesses version of the first network does not exhibit this issue. So this seems to be coming from the Symbolics codegen or RuntimeGeneratedFunctions.

```julia
julia> function lotkavolterradirect(; T = 100.0)
           r1(u,p,t) = growth(p[1], u[1], u[2], p[2]) * u[1]
           r2(u,p,t) = growth(p[1], u[2], u[1], p[2]) * u[2]
           aff1!(int) = int.u[1] += 1
           aff2!(int) = int.u[2] += 1
           p = (.1, 1000)
           u0 = [50, 100]
           tspan = (0.0, T)
           dprob = DiscreteProblem(u0, tspan, p)
           jprob = JumpProblem(dprob, Direct(), ConstantRateJump(r1,aff1!), ConstantRateJump(r2, aff2!); save_positions = (false,false))
           jprob
       end
lotkavolterradirect (generic function with 1 method)

julia> jprob = lotkavolterradirect(T = 10.0)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 2
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 0

julia> sol = solve(jprob, SSAStepper())
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
  0.0
 10.0
u: 2-element Vector{Vector{Int64}}:
 [50, 100]
 [102, 216]

julia> @btime solve($jprob, SSAStepper())
  5.858 μs (11 allocations: 1.03 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
  0.0
 10.0
u: 2-element Vector{Vector{Int64}}:
 [50, 100]
 [109, 228]

julia> jprob2 = lotkavolterradirect(T = 100.0)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 2
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 0

julia> @btime solve($jprob2, SSAStepper())
  30.625 μs (11 allocations: 1.03 KiB)
retcode: Success
Interpolation: Piecewise constant interpolation
t: 2-element Vector{Float64}:
   0.0
 100.0
u: 2-element Vector{Vector{Int64}}:
 [50, 100]
 [402, 598]

```

Example:

```julia
julia> jprob = lotkavolterra(T = 10.0)
JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 2
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 0

julia> aff1 = jprob.discrete_jump_aggregation.affects![1]
RuntimeGeneratedFunction(#=in ModelingToolkit=#, #=using ModelingToolkit=#, :((var"##MTKIntegrator#1026",)->begin
          #= /Users/isaacsas/.julia/packages/SymbolicUtils/1JRDc/src/code.jl:350 =#
          #= /Users/isaacsas/.julia/packages/SymbolicUtils/1JRDc/src/code.jl:351 =#
          #= /Users/isaacsas/.julia/packages/SymbolicUtils/1JRDc/src/code.jl:352 =#
          begin
              ˍ₋out = (var"##MTKIntegrator#1026").u
              ˍ₋arg1 = (var"##MTKIntegrator#1026").u
              ˍ₋arg2 = (var"##MTKIntegrator#1026").p
              t = (var"##MTKIntegrator#1026").t
              begin
                  begin
                      #= /Users/isaacsas/.julia/packages/Symbolics/VIBnK/src/build_function.jl:520 =#
                      #= /Users/isaacsas/.julia/packages/SymbolicUtils/1JRDc/src/code.jl:399 =# @inbounds begin
                              #= /Users/isaacsas/.julia/packages/SymbolicUtils/1JRDc/src/code.jl:395 =#
                              ˍ₋out[1] = (+)(1, ˍ₋arg1[1])
                              #= /Users/isaacsas/.julia/packages/SymbolicUtils/1JRDc/src/code.jl:397 =#
                              nothing
                          end
                  end
              end
          end
      end))

julia> int = (u = [1.0], p = (1.0, 1.0, 1.0), t = 1.0)
(u = [1.0], p = (1.0, 1.0, 1.0), t = 1.0)

julia> @btime aff1($int)
  14.194 ns (1 allocation: 48 bytes)

julia> @btime aff1($int)
  14.194 ns (1 allocation: 48 bytes)

```

---

<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:** [March 27, 2023, 4:35pm UTC](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656/20 "2023-03-27T16:35:58Z")

</div>

If you take the generated function and just eval it then does it go away?

[Next page](https://discourse.julialang.org/t/memory-usage-while-solving-jumpproblem-increases-with-run-time-even-with-save-positions-false-false/96656.md?page=2)