You might want to use an online meanvar calculation. The ensemble analysis stuff has an implementation of Welford:
https://github.com/SciML/DiffEqBase.jl/blob/master/src/ensemble/ensemble_analysis.jl#L173-L214
You can define a reduction that does something similarly without ever holding the full solution in memory.
@isaacsas: I tried:
jprob = JumpProblem(js, dprob, NRM(), save_positions=(false, false))
sol = solve(jprob, SSAStepper(), saveat=[0.1, 0.2, 0.3])
And although I do get the solution at the times specified by saveat, I get all the other time levels as well. That would suggest that something is wrong with the saveat argument of JumpProblem.
Here is sol.t after the calls above:
sol.t
570-element Array{Float64,1}:
0.0
0.09203065154447765
0.22198751565820424
0.1
⋮
17.412720644069154
Yeah, looks like we aren’t passing the option along in the MT code. Sorry about that. The MT jump generation is only a few weeks old, so you’re hitting some cases that just haven’t been looked at yet… I’ll add that to the fix for the performance of JumpProblem.
OK, fixes for both your issues are in. I’ll let you know once new package releases are available. Thanks for letting us know about these issues, and please let us know if you encounter any other problems!
I will be happy to keep reporting problems. Thank you for fixing this. I assume that and “update” command in Pkg will not provide me the new version automatically?
Once the new package versions are out update should pull them down automatically. They just aren’t released yet.
If you update DiffEqJump to v6.7.7 and ModelingToolkit v3.7.1 you should get all the fixes. After you update you can confirm you have the both using the status command in the package manager.
Got them! Thanks. Now let us consider my benchmarks. Here are timings with 30 nodes of degree 5:
0.000212 seconds (339 allocations: 21.031 KiB), Completed: DiscreteProblem
0.043566 seconds (117.32 k allocations: 3.456 MiB), Completed: JumpProblem
0.000082 seconds (95 allocations: 73.281 KiB), Completed: solve
and with 300 nodes, degree 5:
(0.003, 0.389, 0.000915) (JumpProblem is now taking 0.73 sec 30 min later. Weird)
and with 3000 nodes, degree 5
(.29, 5.5, 0.02)
(scaling is sublinear in the number of nodes).
Keeping 3000 nodes, but with degree 10:
(.81, 10.0, 0.03)
JumpProblem scales with the degree (if fixed). I would expect that.
All in all, a great update!
Gordon
Here are the results on a larger network: 20,000 nodes, with 5 edges per nodes, so 100,000 edges.
DiscreteProblem: 22 sec
JumpProblem: 140 sec
Solve: 20 executions: 4 sec (very good timing)
So my question is: while the speed is much much faster than before, I am wondering: is that what you expect? Do you have a feel for what the speed should be? I guess I would have to try it out in Python with Networkx in order to compare?
Yeah, it’s still surprising to me that the assembly of the problem takes so long. Seeing a timing from another system could be helpful.
I checked DiscretizationProblem and found the following, where, somehow all the time is spent. Would you agree? For a system of 100,000 equations, where would the time be spent?
struct DiscreteProblem{uType,tType,isinplace,P,F,K} <: AbstractDiscreteProblem{uType,tType,isinplace}
"""The function in the map."""
f::F
"""The initial condition."""
u0::uType
"""The timespan for the problem."""
tspan::tType
"""The parameter values of the function."""
p::P
""" A callback to be applied to every solver which uses the problem."""
kwargs::K
@add_kwonly function DiscreteProblem{iip}(f::AbstractDiscreteFunction{iip},
u0,tspan::Tuple,p=NullParameters();
kwargs...) where {iip}
_tspan = promote_tspan(tspan)
new{typeof(u0),typeof(_tspan),isinplace(f,4),
typeof(p),
typeof(f),typeof(kwargs)}(f,u0,_tspan,p,kwargs)
end
Chris, I looked at the output of JumpProblem, and here it is for a graph of 5 vertices and 10 edges. Imagine how long it would be when the graph has 20,000 vertices and 200,000 edges. Of course, I should not print it. But here is the source of inefficiencies. The equations are stored as an array of equations and I would not be surprised if that led to cache problems:
MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(γ₂, Pair{Operation,Int64}[I₂ => 1], Pair{Operation,Int64}[I₂ => -1, R₂ => 1]),
Below is the output from JumpProblem, which must be processed by DiscreteEquations. Could you point me to the specific function called in DiffEqBase when calling DiscreteProblem? I understand that types are handled very efficiently, but processing this complex datastructure for a complex graph must take lots of time.
JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}
js= JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}(MassActionJump[MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₁, Pair{Operation,Int64}[I₁ => 1, S₂ => 1], Pair{Operation,Int64}[I₂ => 2, S₂ => -1, I₁ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₂, Pair{Operation,Int64}[I₁ => 1, S₄ => 1], Pair{Operation,Int64}[S₄ => -1, I₄ => 2, I₁ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₃, Pair{Operation,Int64}[I₂ => 1, S₃ => 1], Pair{Operation,Int64}[I₂ => -1, I₃ => 2, S₃ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₄, Pair{Operation,Int64}[I₂ => 1, S₅ => 1], Pair{Operation,Int64}[I₂ => -1, S₅ => -1, I₅ => 2]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₅, Pair{Operation,Int64}[I₃ => 1, S₂ => 1], Pair{Operation,Int64}[I₂ => 2, S₂ => -1, I₃ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₆, Pair{Operation,Int64}[I₃ => 1, S₅ => 1], Pair{Operation,Int64}[S₅ => -1, I₃ => -1, I₅ => 2]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₇, Pair{Operation,Int64}[I₄ => 1, S₂ => 1], Pair{Operation,Int64}[I₂ => 2, I₄ => -1, S₂ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₈, Pair{Operation,Int64}[I₄ => 1, S₅ => 1], Pair{Operation,Int64}[S₅ => -1, I₄ => -1, I₅ => 2]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₉, Pair{Operation,Int64}[I₅ => 1, S₂ => 1], Pair{Operation,Int64}[I₂ => 2, S₂ => -1, I₅ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(β₁₀, Pair{Operation,Int64}[I₅ => 1, S₄ => 1], Pair{Operation,Int64}[S₄ => -1, I₄ => 2, I₅ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(γ₁, Pair{Operation,Int64}[I₁ => 1], Pair{Operation,Int64}[R₁ => 1, I₁ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(γ₂, Pair{Operation,Int64}[I₂ => 1], Pair{Operation,Int64}[I₂ => -1, R₂ => 1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(γ₃, Pair{Operation,Int64}[I₃ => 1], Pair{Operation,Int64}[R₃ => 1, I₃ => -1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(γ₄, Pair{Operation,Int64}[I₄ => 1], Pair{Operation,Int64}[I₄ => -1, R₄ => 1]), MassActionJump{Operation,Array{Pair{Operation,Int64},1},Array{Pair{Operation,Int64},1}}(γ₅, Pair{Operation,Int64}[I₅ => 1], Pair{Operation,Int64}[R₅ => 1, I₅ => -1])]ConstantRateJump[]VariableRateJump[], t, Variable[S₁, S₂, S₃, S₄, S₅, I₁, I₂, I₃, I₄, I₅, R₁, R₂, R₃, R₄, R₅], Variable[β₁, β₂, β₃, β₄, β₅, β₆, β₇, β₈, β₉, β₁₀, γ₁, γ₂, γ₃, γ₄, γ₅], Symbol("##ReactionSystem#375"), JumpSystem[])
If you @which DiscreteProblem(...) you’ll see that it’s not hitting DiffEqBase here, it’s probably all in:
Is there a better place for this discussion, which is getting beyond “usage”? Just checking?
Thanks.
I found:
# DiscreteProblem from AbstractReactionNetwork
function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0::Union{AbstractArray, Number}, tspan::Tuple, p=nothing, args...; kwargs...)
u0 = typeof(u0) <: Array{<:Pair} ? u0 : Pair.(rs.states,u0)
p = typeof(p) <: Array{<:Pair} ? p : Pair.(rs.ps,p)
return DiscreteProblem(convert(JumpSystem,rs), u0,tspan,p, args...; kwargs...)
end
My guess is that the problem is in convert():
function Base.convert(::Type{<:JumpSystem},rs::ReactionSystem)
eqs = assemble_jumps(rs)
JumpSystem(eqs,rs.iv,rs.states,rs.ps,name=rs.name,
systems=convert.(JumpSystem,rs.systems))
end
and maybe in assemble_jumps.
and from there, assemble_jumps
So another question is: what is the best way to explore ReactionSystems? Since debugging is an art with Julia in Atom (how do you do it?), what approach should I take?
If you’re trying to trace to find out what’s taking all of the time, I recommend using the profiler like macOS Python with numpy faster than Julia in training neural network - Stack Overflow
This needs profiling to see what is rate limiting. You could look at documentation for the Julia Profile package, and read about how to profile in Juno. I can profile it when I have some time to see if there is further room for improvement. These jumps ultimately get merged into one vector within DiffEqJump, so perhaps we can rework this to vectorize everything at the ModelingToolkit level if that will help.
Agree with the profiling. In addition - are you doing anything with the graph itself, or is it just a data structure you’re passing on to some other package? What graph algorithms are you using?
The graph isn’t involved very much here: it’s a nice representation of the model but essentially no graph algorithms are involved (well, there is something in the solving)
Thanks. It would be interesting to see the SIR data cached per-node using MetaGraphs.
I have tried some profiling on a small graph. Here is an image.
GraphViz.jl module (tried two different ways via Pkg). which is required. The graph structure is only used to simply the creation of the Reaction equations. I doubt it has anything to do with the inefficiencies.
Here are timings on a network with 3000 nodes, 3 edges per node:
0.005074 seconds (20.00 k allocations: 734.812 KiB)
++++ Completed convert to ReactionSystem
0.576390 seconds (839.05 k allocations: 33.172 MiB, 8.06% gc time)
++++ Completed convert to JumpSystem
0.147043 seconds (35.77 k allocations: 2.004 MiB)
++++ Completed DiscreteProblem
3.863394 seconds (7.91 M allocations: 225.916 MiB, 2.37% gc time)
++++ Completed JumpProblem
I view DiscreteProblem as taking too long. And it has nothing to do with LightGraphs.jl, but rather with ModelingToolkit. What I should probably do is create my equations via some other means and look at the timings.
Note that I have parameters for each edge, but do not see why the timings would be affected by that.
Here are the timings on a network with 30,000 nodes and 10 edges per node.
So 10x more nodes, and 3x more edges. How does timing scale?
0.066596 seconds (279.00 k allocations: 10.209 MiB)
++++ Completed convert to ReactionSystem
78.693881 seconds (11.48 M allocations: 516.547 MiB, 0.93% gc time)
++++ Completed convert to JumpSystem
103.081751 seconds (239.78 k allocations: 14.486 MiB)
++++ Completed DiscreteProblem
299.844236 seconds (152.86 M allocations: 4.403 GiB, 0.90% gc time)
++++ Completed JumpProblem
DiscreteProblem took 300x longer, while JumpProblem took 300 sec, or about 100x longer than the smaller graph. Something not quite right.