start depgraphs for AbstractSys

Author: isaacsas

Project.toml

@@ -12,6 +12,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
+LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"

src/ModelingToolkit.jl

@@ -101,6 +101,7 @@ include("systems/optimization/optimizationsystem.jl")
 include("systems/pde/pdesystem.jl")
 
 include("systems/reaction/reactionsystem.jl")
+include("systems/dependency_graphs.jl")
 
 include("latexify_recipes.jl")
 include("build_function.jl")

src/systems/dependency_graphs.jl

@@ -0,0 +1,33 @@
+
+# each system type should define extract_variables! for a single equation
+function extract_variables(eqs, vars)
+    deps = Set{Variable}()
+    depeqs_to_vars = Vector{Vector{Variable}}(undef,length(eqs))
+
+    for (i,eq) in enumerate(eqs)        
+        depeqs_to_vars[i] = collect(extract_variables!(deps, eq, vars))
+        empty!(deps)
+    end
+end
+
+equation_dependencies(sys) = extract_variables(equations(sys), states(sys))
+
+# map each variable to the eqs depending on it 
+function variables_to_depeqs(sys::AbstractSystem; equationdeps = nothing, statestoids = nothing)
+    sts   = states(sys)
+    stoi  = isnothing(statestoids) ? Dict(convert(Variable,state) => i for (i,state) in enumerate(sts)) : statestoids
+
+    # map from eqs to states they depend on
+    eqdeps = isnothing(equationdeps) ? equation_dependencies(equations(sys), sts) : equationdeps
+
+    # reverse map and switch to integer indices of states
+    dg = [Vector{Int}() for i = 1:length(sts)]
+    for (k,dep) in enumerate(eqdeps)
+        for state in dep
+            push!(dg[stoi[state]],k)
+        end
+    end
+    foreach(dep -> sort!(dep), dg)
+
+    dg
+end

src/systems/jumps/jumpsystem.jl

@@ -103,6 +103,7 @@ function DiffEqBase.DiscreteProblem(sys::AbstractSystem, u0map, tspan::Tuple,
     DiscreteProblem(df, u0, tspan, p; kwargs...)
 end
 
+
 """
 ```julia
 function DiffEqBase.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
@@ -111,59 +112,34 @@ function DiffEqBase.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
 Generates a JumpProblem from a JumpSystem.
 """
 function DiffEqJump.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
-    vrjs = Vector{VariableRateJump}()
-    crjs = Vector{ConstantRateJump}()
-    majs = Vector{MassActionJump}()
-    pvars = parameters(js)
-    statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
-    parammap = map((x,y)->Pair(x(),y),pvars,prob.p)
-    eqs      = equations(js)
 
-    foreach(j -> push!(majs, assemble_maj(js, j, statetoid, parammap)), eqs.x[1]) 
-    foreach(j -> push!(crjs, assemble_crj(js, j, statetoid)), eqs.x[2]) 
-    foreach(j -> push!(vrjs, assemble_vrj(js, j, statetoid)), eqs.x[3]) 
+    statetoid = Dict(convert(Variable,state) => i for (i,state) in enumerate(states(js)))
+    parammap  = map((x,y)->Pair(x(),y), parameters(js), prob.p)
+    eqs       = equations(js)
 
+    majs = MassActionJump[assemble_maj(js, j, statetoid, parammap) for j in eqs.x[1]]
+    crjs = ConstantRateJump[assemble_crj(js, j, statetoid) for j in eqs.x[2]]
+    vrjs = VariableRateJump[assemble_vrj(js, j, statetoid) for j in eqs.x[3]]
     ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps")
+    
     jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, isempty(majs) ? nothing : majs)
     JumpProblem(prob, aggregator, jset)
 end
 
 
-### Functions to determine which jumps depend on a given state
+### Functions to determine which states a jump depends on
 
-function ratevars!(dep, jump::Union{ConstantRateJump,VariableRateJump}, sts)
+function extract_variables!(dep, jump::Union{ConstantRateJump,VariableRateJump}, sts)
     foreach(st -> (convert(Variable,st) in sts) && push!(dep,st), vars(jump.rate))
+    dep
 end
 
-function ratevars!(dep, jump::MassActionJump, args...)    
-    foreach(st -> push!(dep,convert(Variable,st[1])), jump.reactant_stoch)
-end
-
-function ratevars(jump, sts)
-    dep = Set{Variable}()
-    ratevars!(dep,jump,sts)
+function extract_variables!(dep, jump::MassActionJump, args...)    
+    foreach(st -> push!(dep,convert(Variable,st[1])), jump.reactant_stoch)    
     dep
 end
 
-# map each state (as an Int index) to the jumps with a rate that depend on it 
-function statetojump_depgraph(js::JumpSystem, statestoids = nothing)
-    sts   = states(js)
-    stoi  = isnothing(statestoids) ? Dict(convert(Variable,state) => i for (i,state) in enumerate(sts)) : statestoids
 
-    # map from jump to states (as Variables) it depends on through the rate
-    deps  = [ratevars(eq,sts) for eq in equations(js)]    
-
-    # reverse map and switch to integer indices of states
-    dg = [Vector{Int}() for i = 1:length(sts)]
-    for (k,dep) in enumerate(deps)
-        for state in dep
-            push!(dg[stoi[state]],k)
-        end
-    end
-    foreach(dep -> sort!(dep), dg)
-
-    dg
-end
 
 ### Functions to determine which states are modified by a given jump