Efficient memory allocation for large system of ODEs



I have a system of differential equations that mimic the clonal population in an immune system. For whoever is curious, the problem definition can be found in the post:

At present, I am simulating with a small number of clones and thus, am able to deal with the memory required. But, for higher number of clones, the performance of my present code is pretty detrimental. Say, for 2000 clones (variable n), the code runs for about 15 hours and stops prematurely because of the shortage of memory. As advised by @ChrisRackauckas in his tutorials, I tried saving only the last point of the simulation, but again faced the same situation.

The code is clearly memory inefficient. I tried using minimizing global variables by defining the function main. But, the memory and the time required increases for these cases. What else can be done for managing the memory?


Code: https://www.pastiebin.com/5b0d2ad13e206

PositiveDomain callback and isoutofdomain implementation

Have you read (and maybe re-read) the “Performance tips” section of the manual? Have you profiled your code: https://docs.julialang.org/en/stable/manual/profile/#Profiling-1?

p = [n, g, p2, p3, p4, p5, p6, h1, h2, z];

n is an integer which makes this array not of the same types, which makes it typeof(p) = Vector{Any}, which makes every parameter in your inner loop not inferred and thus your whole code not type-stable. I would instead recommend using a tuple or something like Parameters.jl. Just pull that inner function aside to see:

n = 4; # number of clones
h1 = [841.0, 292.0, 860.0, 882.0];
h2 = h1 .+ [564.0, 374.0, 727.0, 813.0];
g = [0.791508, 0.116203, 0.702549, 0.832428];
p2 = [2.0, 1.0, 3.0, 7.0];
p3 = 0.01+(1-0.01)*[0.739588, 0.546912, 0.725891, 0.0332758];
p4 = 0.001+(1-0.001)*[0.040059, 0.804641, 0.639816, 0.513092];
p5 = 0.001+(1-0.001)*[0.9515, 0.2969, 0.624892, 0.34171];
p6 = 0.01+(1-0.01)*[0.455292, 0.339854, 0.0306784, 0.437118];
z = zeros(n);

y0 = [1.0; vec(hcat([65.0, 14.0, 56.0, 65.0],[0.492182, 0.431219, 0.538526, 0.517137])')];
f = similar(y0)
y = y0
t = 20.0

p = [n, g, p2, p3, p4, p5, p6, h1, h2, z]
@code_warntype Scaled_Null(f,y,p,t)

p = (n, g, p2, p3, p4, p5, p6, h1, h2, z)
@code_warntype Scaled_Null(f,y,p,t)

So if you make that a tuple it should be much better.