I apologize in advance for this length thread, but the problem setting is relatively complicated and this was the smallest MWE that I could think of, while still remaining clear.

### Problem setting

I am encountering some issues with memory usage that is so big that my `Julia`

processes get killed, yet I do not understand why so much memory is being allocated.

Essentially, I am trying to solve (generalized) Lotka-Volterra systems. We can consider the very simple form;

I am using `ModelingToolkit.jl`

in conjunction with `DifferentialEquations.jl`

to create an `EnsembleProblem`

and solve many instances, with random A_{ij} in parallel. As the system at hand is known for being chaotic, for now let us assume it suffices to assume that for the given parameters there is a unique fixed point to which all systems converge. Note that this fixed point obviously depends on all A_{ij}, but not on on the initial conditions. I am interested in the state of the system at some (large) time t, so I only save the final state of the solution β this should also save a lot of memory.

The code I am using to investigate this problem is comprised of some functions. The main two are to generate an `ODEProblem`

and the `EnsembleProblem`

, and there are some βhelperβ functions that update an existing problem and a callback that puts βspeciesβ with low abundance to zero, i.e. x_i = 0 when x_i < \varepsilon.

### The code

The function to generate the `ODEProblem`

:

```
function define_odeproblem(S, a; x0=rand(S), k=1.0, tspan=(0.0, 1_000.0))
@variables t
@variables (x(t))[1:S]
@parameters A[1:S,1:S]
D = Differential(t)
eqns = [D(x[i]) ~ x[i] * (1.0 - sum([A[i,j]*x[j] for j in 1:S])) for i in 1:S]
@named odesys = ODESystem(eqns, t)
params = vec([A[i,j] => a[i,j] for i in 1:S, j in 1:S])
u0 = [x[i] => x0[i] for i in 1:S]
odeprob = ODEProblem(complete(odesys), u0, tspan, params)
return odeprob
end
```

The function to generate the `EnsembleProblem`

:

```
function define_ensembleproblem(S, ΞΌ, Ο; nseeds=8, tspan=(0.0, 1_000.0))
#/ Initialize interactions
#~ (does not really matter how, as they'll be changed in set_interaction(..) anyways
function generate_interactions(S, ΞΌ, Ο; seed=42)
a = ΞΌ/S .+ Ο*randn(Random.Xoshiro(1234*seed),S,S)/sqrt(S)
a[CartesianIndex.(1:S, 1:S)] .= 1.0
return a
end
amatrices = [generate_interactions(S, ΞΌ, Ο) for _ in 1:nseeds]
#/ Define ODEProblem
prob = define_odeproblem(S, amatrices[begin]; tspan=tspan)
function set_interactions(prob, k, nrepeats)
#~ Set interactions randomly and update existing ODEProblem
@unpack A = prob.f.sys
Amap = vec([A[i,j] => amatrices[k][i,j] for i in 1:S, j in 1:S])
return update_prob(prob, Amap)
end
#/ Define EnsembleProblem
#~ Save only final state
output_func(sol, i) = (last(sol), false)
eprob = EnsembleProblem(prob, output_func=output_func, prob_func=set_interactions)
return eprob, amatrices
end
```

and the βhelperβ functions

```
"Update ODEProblem"
function update_prob(prob, pmap)
#~ Updates given ODEProblem with the new parameters given in pmap
#~ Old parameters (not in pmap) remain the same
params = ModelingToolkit.parameters(prob.f.sys)
pdict = Dict(params[n] => prob.p[n] for n in 1:length(prob.p))
p = ModelingToolkit.varmap_to_vars(pmap, params, defaults=pdict)
updated_prob = remake(prob; p=p, u0=rand()*prob.u0)
return updated_prob
end
"Callback to put extinct species to 0.0"
function cb_extinction(threshold::Float64)
#/ Create callback to let species with low abundances go extinct
function condition(u, t, integrator)
any(u .< threshold)
end
function affect!(integrator)
integrator.u[integrator.u .< threshold] .= 0.0
end
return DiscreteCallback(condition, affect!)
end
```

The function that updates the problem is written such that, in principle, not all parameters need to be given for the problem to be updated. One can, for example, give only `A_{11}`

and update the problem in that way. If there is a better way to write this, please let me know.

### The problem

I am encountering issues when increasing the no. of species and/or the number of seeds. The problems are not with the time, but there are many memory allocations of which I do not know the origin. Intuitively, I would expect that perhaps some more temporary memory is used, however I have observed crashes when there is no more RAM available.

Consider the following example (I omitted all `include`

βs for brevity)

```
julia> S = 32; ΞΌ = 4.0; Ο = 0.5;
julia> eprob, amatrices = define_ensembleproblem(S, ΞΌ, Ο; nseeds=32);
julia> esol = solve_ensembleproblem(eprob, length(amatrices))
julia> @benchmark solve_ensembleproblem(eprob, length(amatrices))
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min β¦ max): 516.055 ms β¦ 989.800 ms β GC (min β¦ max): 0.00% β¦ 46.69%
Time (median): 539.820 ms β GC (median): 0.00%
Time (mean Β± Ο): 582.574 ms Β± 143.338 ms β GC (mean Β± Ο): 7.93% Β± 14.77%
ββ
β
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ
β
516 ms Histogram: frequency by time 990 ms <
Memory estimate: 528.11 MiB, allocs estimate: 7988942.
julia> eprob, amatrices = define_ensembleproblem(S, ΞΌ, Ο; nseeds=64);
julia> esol = solve_ensembleproblem(eprob, length(amatrices));
julia> @benchmark solve_ensembleproblem(eprob, length(amatrices))
BenchmarkTools.Trial: 5 samples with 1 evaluation.
Range (min β¦ max): 1.096 s β¦ 1.626 s β GC (min β¦ max): 0.00% β¦ 30.79%
Time (median): 1.130 s β GC (median): 0.00%
Time (mean Β± Ο): 1.222 s Β± 226.789 ms β GC (mean Β± Ο): 8.20% Β± 13.77%
β βββ β
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
1.1 s Histogram: frequency by time 1.63 s <
Memory estimate: 1.03 GiB, allocs estimate: 15977622.
```

First, the total memory that is used is, in my opinion, really large. I suspect that the `ODEProblem`

that is defined using `ModelingToolkit`

is *not* in-place (i.e., it is allocating)? Could this be the underlying issue? If so, how could this be fixed?

This would kind of makes sense as the memory estimate doubles when the no. of seeds doubles. However I explicitly tell the solver to not store anything but the final state, so I would guess that the memory should only increase slightly to accommodate the size of the solution (which is double the size).

Are there further ways I can improve the (memory) efficiency of this code? I am looking to model systems where both S and `nseeds`

are in the order of 10^2 to 10^3, so I need a way to run these large models without crashing, at least.