Extemely slow creation of ODEProblem from ReactionSystem
#1051
Comments
|
Thanks for the example; we'll have to look into this. ModelingToolkit development moves fast these days, and it looks like they have changed around some of the input processing from how it was done in the past. We may need to make updates here to mimic how they now internally process inputs for systems, and how array type symbolics are handled. I'll try to play around with this when I have some time, unfortunately it may not be until next week due to other deadlines. |
|
No worries. I just found it very odd as it was not something I expected and, I believe, it worked for me before without issue. I am not super familiar with the |
|
ModelingToolkit 9 has made a lot of under the hood changes, including how array symbolics are handled I think. It may be that we are not using the right approach with how to handle them any more. You could look at the ODESystem constructors in ModelingToolkit and compare them to the Catalyst ReactionSystem constructor to see differences. It might also be good to actually profile the code and see what step is taking so long. Finally, it might be worth trying the four-argument constructor to |
|
Could you share a StatProfilerHTMlL.jl zip file report? That will make it easy to see what the culprit is. |
|
@johannesnauta just to make sure we understand. There is no issue with the time it takes to make the systems then, it is with creating an If so, can you give us your code for creating the problems that you actually benchmarked. I don't see it above. |
|
Your flamegraph seems to show that all the slowdown is actually in ModelingToolkit building an |
|
Yes exactly. The only issue is when creating the function create_ODEProblem(
rs::ReactionSystem;
S=length(species(rs)),
rv=ones(S), Av=ones(S, S), xv=ones(S),
tspan=(0.0,100.0)
)
pmap = [rs.r => rv, rs.A => Av]
xmap = [rs.x => xv]
prob = ODEProblem(rs, xmap, tspan, pmap)
return prob
end
function create_ODEProblem(
osys::ODESystem;
S=length(equations(osys)),
rv=ones(S), Av=ones(S, S), xv=ones(S),
tspan=(0.0,100.0)
)
pmap = [osys.r => rv, osys.A => Av]
xmap = [osys.x => xv]
prob = ODEProblem(osys, xmap, tspan, pmap)
return prob
end |
|
@isaacsas's guess that it might have to do with the vectorised parameters/species sounds like a good guess, there have been quite a lot of movements there lately. There was quite a lot of discussions a while ago, have forgotten some of the turns, but I did test building rather large systems using vector species/parameters, so it did work fine a couple of months ago. Probably requires some digging to sort out properly. The first step is probably to check whether the issue persist when only parameters/only species are on vector form (I remember these were treated differently), and maybe also with neither on that form. |
|
Seems that it might be connected to symbolic arrays indeed, here's a crude test (it's cumbersome to time out all interaction parameters by hand, obviously, so just a small system of 2 species for now). function create_scalarrs()
t = default_t()
@parameters r1 r2
@parameters A11 A12 A21 A22
@species x1(t) x2(t)
rxs = [
Reaction(r1, [x1], [x1], [1], [2]),
Reaction(r2, [x2], [x2], [1], [2]),
Reaction(r1*A11, [x1], [x1], [2], [1]),
Reaction(r1*A12, [x1,x2], [x2], [1,1], [1]),
Reaction(r2*A21, [x1,x2], [x1], [1,1], [1]),
Reaction(r2*A22, [x2], [x1], [2], [1])
]
x = [x1, x2]
p = [r1, r2, A11, A12, A21, A22]
@named glvrs = ReactionSystem(rxs, t, x, p, combinatoric_ratelaws=false)
return complete(glvrs)
endwhich gives julia> scalrs = create_scalarrs()
julia> scalprob = create_scalarODEProblem(scalrs);
julia> @btime scalprob = create_scalarODEProblem(scalrs);
718.716 μs (7103 allocations: 481.45 KiB)
julia> @btime scalprob = create_scalarODEProblem(scalrs);
720.363 μs (7103 allocations: 481.45 KiB)while using symbolic arrays with the julia> vecrs = create_rs(2);
julia> vecode = create_ODEProblem(vecrs);
julia> @btime vecode = create_ODEProblem(vecrs);
1.309 ms (10633 allocations: 647.30 KiB)
julia> @btime vecode = create_ODEProblem(vecrs);
1.309 ms (10633 allocations: 647.30 KiB) |
|
From what I can tell ModelingToolkit no longer scalarizes, however, unknowns and parameters are handled differently: julia> function create_ODESystem(S::Int)
#/ Define variables and parameters
t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits
@variables (x(t))[1:S]
@parameters r[1:S] A[1:S,1:S]
eqns = [D(x[i]) ~ r[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)
return ModelingToolkit.complete(odesys)
end
create_ODESystem (generic function with 1 method)
julia> osys = create_ODESystem(3)
Model odesys with 3 equations
Unknowns (3):
(x(t))[1]
(x(t))[2]
(x(t))[3]
Parameters (2):
A
rNotice that while it isn't scalarizing it is detecting and explicitly storing each component of the vector unknown We only scalarize in a few places, but I have no idea how this might mess things up in Catalyst since now the number of parameters is not the same as the length of the parameter vector (which we might assume), and we'd now have to handle in places expanding a vector |
|
@johannesnauta I just wanted to point out that it looks like your systems aren't actually exactly the same? When I take the difference of the generated ODEs rhs and simplify them I get some leftover terms: 3-element Vector{Any}:
0
A[2, 1]*r[1]*(x(t))[1]*(x(t))[2] - A[2, 1]*r[2]*(x(t))[1]*(x(t))[2]
A[3, 1]*r[1]*(x(t))[1]*(x(t))[3] - A[3, 1]*r[3]*(x(t))[1]*(x(t))[3] + A[3, 2]*r[2]*(x(t))[2]*(x(t))[3] - A[3, 2]*r[3]*(x(t))[2]*(x(t))[3](This is the ODE model minus the converted reaction system model when I guess you aren't exactly setting the rate consistently in the two models? |
|
So good news; I think it should be relatively easy to make Catalyst match MTK's new conventions, see #1052. Bad news; it appears our tests will need lots of modifications if we make such a change. Unfortunately I don't think this is something I have the time to handle in the next few weeks due to other deadlines, but that PR seems to get things fixed in terms of the generated ODESystem if you want to test it out on your example. Now someone has to go through all the tests and update them accordingly to no longer assume all inputs are scalars. |
|
This will probably also require a breaking release of Catalyst once updated since it does change how parameters are represented internally and stored in what is passed to users (only returning the vector/matrix symbolics instead of their components). |
I was indeed sloppy with the coding, and in this case there as additionally a typo. The original snippets only gave equal equations when the growth rates are function create_rs(S::Int)
t = default_t()
@parameters r[1:S] A[1:S,1:S]
@species (x(t))[1:S]
#/ Allocate and fill
growthrxs = Array{Reaction}(undef, S)
interactionrxs = Array{Reaction}(undef, S, S)
for i in 1:S
#~ Growth
growthrxs[i] = Reaction(r[i], [x[i]], [x[i]], [1], [2])
#~ Diagonal terms
interactionrxs[i,i] = Reaction(r[i]*A[i,i], [x[i]], [x[i]], [2], [1])
for j in (i+1):S
#~ Off-diagonal terms
interactionrxs[i,j] = Reaction(r[i]*A[i,j], [x[i], x[j]], [x[j]], [1, 1], [1])
interactionrxs[j,i] = Reaction(r[j]*A[j,i], [x[j], x[i]], [x[i]], [1, 1], [1])
end
end
#/ Gather reactions
rxs = [growthrxs..., interactionrxs...]
@named glvrs = ReactionSystem(rxs, t, combinatoric_ratelaws=false)
return Catalyst.complete(glvrs)
end
function create_ODESystem(S::Int)
#/ Define variables and parameters
@variables (x(t))[1:S]
@parameters r[1:S] A[1:S,1:S]
eqns = [D(x[i]) ~ x[i]*(r[i] - sum(r[i]*[A[i,j]*x[j] for j in 1:S])) for i in 1:S]
@named odesys = ODESystem(eqns, t)
return ModelingToolkit.complete(odesys)
endGiving, julia> S = 3;
julia> rs = create_rs(S);
julia> rsode = convert(ODESystem, rs);
julia> rseqs = [eq.rhs for eq in equations(rsode)];
julia> mtkode = create_ODESystem(S);
julia> mtkeqs = [eq.rhs for eq in equations(mtkode)];
julia> expand.(mtkeqs) .- rseqs
3-element Vector{Int64}:
0
0
0 |
I will test these out and report back to you when I find some time later today. Thanks a lot in any case! |
|
Please do let us know what you find. If there is still a big difference in calling ODEProblem with a directly built ODESystem vs. one you have converted from Catalyst I'll have to dig deeper. |
|
Yes it does seem to be that there is still a big difference, sadly. julia> S = 8;
julia> rs = create_rs(S);
julia> rsprob = create_ODEProblem(rs);
julia> sys = create_ODESystem(S);
julia> sysprob = create_ODEProblem(sys);
julia> @btime create_ODEProblem(rs);
65.161 ms (360819 allocations: 19.75 MiB)
julia> @btime create_ODEProblem(sys);
4.152 ms (50613 allocations: 2.87 MiB)also note that the parameters are still scalarized julia> converted_rs = convert(ODESystem, rs)
Model glvrs with 8 equations
Unknowns (8):
(x(t))[1]
(x(t))[2]
(x(t))[3]
(x(t))[4]
(x(t))[5]
(x(t))[6]
(x(t))[7]
(x(t))[8]
Parameters (72):
r[1]
r[2]
r[3]
r[4]
r[5]
r[6]
r[7]
r[8]
A[1, 1]
⋮
julia> sys
Model odesys with 8 equations
Unknowns (8):
(x(t))[1]
(x(t))[2]
(x(t))[3]
(x(t))[4]
(x(t))[5]
(x(t))[6]
(x(t))[7]
(x(t))[8]
Parameters (2):
A
rIn this case, I guess that the fact that the interaction matrix growth as |
|
But where is the slowdown? Is it in |
|
We add terms one by one into the ODEs. So if the slow part is the |
|
The slowdown is in creating the |
|
Wait, parameters are still scalarized using the PR? OK, that seems like I missed something then. When I had tested it that wasn't the case. |
|
@johannesnauta I don't think you used the PR branch (it isn't merged into master). I ran your example with it and the parameters are not expanded for me. For julia> @btime create_ODEProblem($osysrs)
3.581 ms (56764 allocations: 3.22 MiB)
julia> @btime create_ODEProblem($osys)
3.379 ms (50579 allocations: 2.87 MiB)So a much smaller difference. |
|
Wait, then how should one use a specific version/commit of a package? What I did was just pkg> add Catalyst#6dde4d6fbc1c61d011d744eda2256d0ca420323fwhich I believe to be related to #1052, correct? Also when I looked through the source code, the |
|
That is the correct PR but I think you maybe only pulled the first commit on it instead of the latest commit? I usually just checkout PRs via the Github Desktop app. |
|
The latest commit is 871d3da |
|
Ah I assumed commits on Github were listed in reverse chronological order, but apparently they are chronological instead. Your fix gives me similar results indeed, with the |
Describe the bug 🐞
Me again. I am simulating Lotka-Volterra systems, and have recently jumped back on
Catalyst.jlfor the flexibility and comparison with other solvers (SDEs, SSAs, etc.) other than solving ODEs. However, I find that creating (not solving) theODEProblemfrom aReactionSystemthat is sufficiently large takes an incredible amount of time. The time it takes reaches the point where trying to create a problem of medium sizes (~50 species) basically does not finish. Additionally, when using anODESystemusingModelingToolkit.jldirectly, and creating theODEProblemdoes not seem to have this issue.Minimal Reproducible Example 👇
To illustrate, consider the following functions, one using
Catalyst.jl, and one usingModelingToolkit.jl. They both create a simple Lotka-Volterra system, that isUsing
Catalyst.jl(note: if there are other/better ways to generate the
ReactionSystem, let me know!)Using
ModelingToolkit.jlThen, running the following:
From this the problem is already apparent; about two orders of magnitude longer time is needed to create the
ODEProblemfrom theReactionSystem. This issue is even more pressing the larger the no. of speciesS.Now, I know I have read somewhere that it is expected to take a bit longer. But even converting the
ReactionSystemto anODESystem, i.e. usingconvert(ODESystem, rs), and then making theODEProblemdoes not change anything for me.Can someone explain to me what is going on? Am I missing something? Why does creating the
ODEProblemtakes so much time (and memory)? Why is this different from theModelingToolkit.jlimplementation?I would like to learn how to do this properly, as I really would like to use the
ReactionSystemto create other types of problems, such asSDEProblem, andJumpProblem-- which I think is one of the main selling points for usingCatalyst.jl.Any help or comments are greatly appreciated!
Environment (please complete the following information):
Details
The text was updated successfully, but these errors were encountered: