Merge remote-tracking branch 'origin/master' into add_no_param_tests

Project.toml

@@ -10,10 +10,10 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 
 [compat]
-Catlab = "0.7.2"
+Catlab = "0.7.3"
 Latexify = "0.13.5"
-MacroTools = "0.5"
-ModelingToolkit = "3.14"
+MacroTools = "0.5.5"
+ModelingToolkit = "3.14.2"
 Reexport = "0.2"
 julia = "1.3"
 

src/graphs.jl

@@ -17,10 +17,10 @@ function edgify(δ, i, reverse::Bool)
 end
 
 # make distinguished edge based on rate constant
-function edgifyrates(rn)    
+function edgifyrates(rxs, specs)    
     es = Edge[]
-    for (i,rx) in enumerate(reactions(rn))
-        deps = get_variables(rx.rate, states(rn))
+    for (i,rx) in enumerate(rxs)
+        deps = get_variables(rx.rate, specs)
         for dep in deps
             val = String(dep.op.name)
             attr = Attributes(:color => "#d91111", :style => "dashed")
@@ -48,20 +48,21 @@ Notes:
   red and black arrows from `A` to the reaction node.
 """
 function Graph(rn::ReactionSystem)
-    rxs = reactions(rn)
-    statenodes = [Node(string(s.name), Attributes(:shape=>"circle", :color=>"#6C9AC3")) for s in species(rn)]
+    rxs   = reactions(rn)
+    specs = species(rn)
+    statenodes = [Node(string(s.name), Attributes(:shape=>"circle", :color=>"#6C9AC3")) for s in specs]
     transnodes = [Node(string("rx_$i"), Attributes(:shape=>"point", :color=>"#E28F41", :width=>".1")) for (i,r) in enumerate(rxs)]
 
     stmts = vcat(statenodes, transnodes)
     edges = map(enumerate(rxs)) do (i,r)
       vcat(edgify(zip(r.substrates,r.substoich), i, false),
            edgify(zip(r.products,r.prodstoich), i, true))
     end
-    es = edgifyrates(rn)
-    (!isempty(es)) && push!(edges, edgifyrates(rn))
+    es = edgifyrates(rxs, specs)
+    (!isempty(es)) && push!(edges, es)
 
     stmts = vcat(stmts, collect(flatten(edges)))
-    g = Graphviz.Graph("G", true, stmts, graph_attrs, node_attrs,edge_attrs)
+    g = Graphviz.Digraph("G", stmts; graph_attrs=graph_attrs, node_attrs=node_attrs, edge_attrs=edge_attrs)
     return g
 end
 

src/networkapi.jl

@@ -6,27 +6,39 @@
     species(network)
 
 Given a [`ReactionSystem`](@ref), return a vector of species `Variable`s.
+
+Notes:
+- If `network.systems` is not empty may allocate. Otherwise returns
+  `network.states`.
 """
 function species(network)
-    states(network)
+    isempty(network.systems) ? network.states : states(network)
 end
 
 """
     params(network)
 
 Given a [`ReactionSystem`](@ref), return a vector of parameter `Variable`s.
+
+Notes:
+- If `network.systems` is not empty may allocate. Otherwise returns
+  `network.ps`.
 """
 function params(network)
-    parameters(network)
+    isempty(network.systems) ? network.ps : parameters(network) 
 end
 
 """
     reactions(network)
 
 Given a [`ReactionSystem`](@ref), return a vector of all `Reactions` in the system.
+
+Notes:
+- If `network.systems` is not empty may allocate. Otherwise returns
+  `network.eqs`.
 """
 function reactions(network)
-    equations(network)
+    isempty(network.systems) ? network.eqs : equations(network)
 end
 
 """
@@ -55,7 +67,11 @@ end
 Return the number of species within the given [`ReactionSystem`](@ref).
 """
 function numspecies(network)
-    length(species(network))
+    ns = length(network.states)
+    for sys in network.systems
+        ns += numspecies(ns)
+    end
+    ns
 end
 
 """
@@ -64,7 +80,11 @@ end
 Return the number of reactions within the given [`ReactionSystem`](@ref).
 """
 function numreactions(network)
-    length(equations(network))
+    nr = length(network.eqs)
+    for sys in network.systems
+        nr += numreactions(sys)
+    end
+    nr
 end
 
 """
@@ -73,7 +93,11 @@ end
 Return the number of parameters within the given [`ReactionSystem`](@ref).
 """
 function numparams(network)
-    length(params(network))
+    np = length(network.ps)
+    for sys in network.systems
+        np += numparams(ns)
+    end
+    np
 end
 
 """
@@ -89,10 +113,11 @@ the returned vector would be `[W(t),X(t),Y(t)]`.
 
 Notes:
 - Allocates
+- Does not check for dependents within any subsystems.
 """
 function dependents(rx, network)
     if rx.rate isa Operation
-        rvars = ModelingToolkit.get_variables(rx.rate, states(network))
+        rvars = ModelingToolkit.get_variables(rx.rate, species(network))
         return union!(rvars, rx.substrates)
     end
 
@@ -139,6 +164,7 @@ order network components were defined.
 Notes:
 - *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
     considered different than `(A+1)^2`.
+- Flattens subsystems, and hence may allocate, when checking equality.
 """
 function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
     issetequal(species(rn1), species(rn2)) || return false
@@ -148,11 +174,12 @@ function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
 
     # the following fails for some reason, so need to use issubset
     #issetequal(equations(rn1), equations(rn2)) || return false
-    (issubset(equations(rn1),equations(rn2)) && issubset(equations(rn2),equations(rn1))) || return false
+    (issubset(reactions(rn1),reactions(rn2)) && issubset(reactions(rn2),reactions(rn1))) || return false
 
+    # BELOW SHOULD NOT BE NEEDED as species, params and equations flatten
     #issetequal(rn1.systems, rn2.systems) || return false
-    sys1 = rn1.systems; sys2 = rn2.systems
-    (issubset(sys1,sys2) && issubset(sys2,sys1)) || return false
+    # sys1 = rn1.systems; sys2 = rn2.systems
+    # (issubset(sys1,sys2) && issubset(sys2,sys1)) || return false
     true
 end
 
@@ -165,7 +192,7 @@ end
 Construct an empty [`ReactionSystem`](@ref). `iv` is the independent variable, usually time.
 """
 function make_empty_network(; iv=Variable(:t))
-    ReactionSystem(Reaction[], iv, Operation[], Operation[], gensym(:ReactionSystem), ReactionSystem[])
+    ReactionSystem(Reaction[], iv, Operation[], Operation[], Variable[], Equation[], gensym(:ReactionSystem), ReactionSystem[])
 end
 
 """
@@ -182,10 +209,12 @@ Notes:
   variable, as this will potentially leave the system in an undefined state.
 """
 function addspecies!(network::ReactionSystem, s::Variable; disablechecks=false)
-    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, s), species(network))
+
+    # we don't check subsystems since we will add it to the top level system...
+    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, s), network.states)
     if curidx === nothing
         push!(network.states, s)
-        return length(species(network))
+        return length(network.states)
     else
         return curidx
     end
@@ -221,10 +250,12 @@ id of the parameter within the system.
   variable, as this will potentially leave the system in an undefined state.
 """
 function addparam!(network::ReactionSystem, p::Variable; disablechecks=false)
-    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, p), params(network))
+
+    # we don't check subsystems since we will add it to the top level system...
+    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, p), network.ps)
     if curidx === nothing
         push!(network.ps, p)
-        return length(params(network))
+        return length(network.ps)
     else
         return curidx
     end
@@ -259,7 +290,7 @@ Notes:
 """
 function addreaction!(network::ReactionSystem, rx::Reaction)
     push!(network.eqs, rx)
-    length(equations(network))
+    length(network.eqs)
 end
 
 
@@ -271,14 +302,15 @@ Merge `network2` into `network1`.
 Notes:
 - Duplicate reactions between the two networks are not filtered out.
 - [`Reaction`](@ref)s are not deepcopied to minimize allocations, so both networks will share underlying data arrays.
+- Subsystems are not deepcopied between the two networks and will hence be shared.
 - Returns `network1`
 """
 function merge!(network1::ReactionSystem, network2::ReactionSystem)
     isequal(network1.iv, network2.iv) || error("Reaction networks must have the same independent variable to be mergable.")
-    specs = species(network1)
-    foreach(spec -> !(spec in specs) && push!(specs, spec), species(network2))
-    ps = params(network1)
-    foreach(p -> !(p in ps) && push!(ps, p), params(network2))
+    specs = network1.states
+    foreach(spec -> !(spec in specs) && push!(specs, spec), network2.states)
+    ps = network1.ps
+    foreach(p -> !(p in ps) && push!(ps, p), network2.ps)
     append!(network1.eqs, network2.eqs)
     append!(network1.systems, network2.systems)
     network1
@@ -292,6 +324,7 @@ Create a new network merging `network1` and `network2`.
 Notes:
 - Duplicate reactions between the two networks are not filtered out.
 - [`Reaction`](@ref)s are not deepcopied to minimize allocations, so the new network will share underlying data arrays.
+- Subsystems are not deepcopied between the two networks and will hence be shared.
 - Returns the merged network.
 """
 function merge(network1::ReactionSystem, network2::ReactionSystem)

src/reaction_network.jl

@@ -83,7 +83,7 @@ end
 #Returns a empty network (with, or without, parameters declared)
 macro reaction_network(parameters...)
     !isempty(intersect(forbidden_symbols,parameters)) && error("The following symbol(s) are used as reactants or parameters: "*((map(s -> "'"*string(s)*"', ",intersect(forbidden_symbols,reactants,parameters))...))*"this is not permited.")
-    return Expr(:block,:(@parameters $((:t,parameters...)...)), :(ReactionSystem(Reaction[], t, Operation[], [$(parameters...)] , gensym(:ReactionSystem), ReactionSystem[])))
+    return Expr(:block,:(@parameters $((:t,parameters...)...)), :(ReactionSystem(Reaction[], t, Operation[], [$(parameters...)] , Variable[], Equation[], gensym(:ReactionSystem), ReactionSystem[])))
 end
 
 """

test/api.jl

@@ -37,6 +37,19 @@ addreaction!(rs2, Reaction(k3, [S], [D]))
 addreaction!(rs2, Reaction(k4, [S,I], [D]))
 @test rs2 == rs
 
+rxs = [Reaction(k1, [S,I], [I], [1,1], [2]),
+       Reaction(k2, [I], [R]) ]
+rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
+rs3 = make_empty_network()
+addspecies!(rs3, S)
+addspecies!(rs3, D.op)
+addparam!(rs3, k3)
+addparam!(rs3, k4.op)
+addreaction!(rs3, Reaction(k3, [S], [D]))
+addreaction!(rs3, Reaction(k4, [S,I], [D]))
+rs4 = merge(rs, rs3)
+@test rs2 == rs4
+
 rxs = [Reaction(k1*S, [S,I], [I], [2,3], [2]),
        Reaction(k2*R, [I], [R]) ]
 rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
@@ -49,4 +62,4 @@ addspecies!(rs, Variable(:S), disablechecks=true)
 addparam!(rs, Variable(:k1))
 @test numparams(rs) == 2
 addparam!(rs, Variable(:k1), disablechecks=true)
-@test numparams(rs) == 3
+@test numparams(rs) == 3

test/graphs.jl

@@ -0,0 +1,14 @@
+using Catalyst
+rn = @reaction_network begin
+    α, S + I --> 2I
+    β, I --> R
+    S^2, R --> 0
+end α β
+
+# check can make a graph
+gr = Graph(rn)
+
+# check can save a graph
+fname = Base.Filesystem.tempname()
+savegraph(gr, fname, "png")
+

test/make_model.jl

@@ -338,12 +338,12 @@ g1 = SDEFunction(convert(SDESystem,reaction_networks_constraint[1]))
 g2 = SDEFunction(convert(SDESystem,time_network))
 for factor in [1e-2, 1e-1, 1e0, 1e1, 1e2, 1e3]
     u0 = factor*rand(length(time_network.states))
-    k2 = factor*rand(); k3 = factor*rand(); k6 = factor*rand();
-    t = rand()
-    p1 = [t, k2, k3, t, t, k6]; p2 = [k2, k3, k6];
-    @test all(abs.(f1(u0,p1,t) .- f2(u0,p2,t)) .< 100*eps())
-    @test all(abs.(f1.jac(u0,p1,t) .- f2.jac(u0,p2,t)) .< 100*eps())
-    @test all(abs.(g1(u0,p1,t) .- g2(u0,p2,t)) .< 100*eps())
+    κ2 = factor*rand(); κ3 = factor*rand(); κ6 = factor*rand();
+    τ = rand()
+    p1 = [τ, κ2, κ3, τ, τ, κ6]; p2 = [κ2, κ3, κ6];
+    @test all(abs.(f1(u0,p1,τ) .- f2(u0,p2,τ)) .< 100*eps())
+    @test all(abs.(f1.jac(u0,p1,τ) .- f2.jac(u0,p2,τ)) .< 100*eps())
+    @test all(abs.(g1(u0,p1,τ) .- g2(u0,p2,τ)) .< 100*eps())
 end
 
 

test/runtests.jl

@@ -31,4 +31,7 @@ using SafeTestsets
 #@time @safetestset "Basic Plotting" begin include("plotting.jl") end
 @time @safetestset "Latexify" begin include("latexify.jl") end
 
+# the following can't really be run until there is an artifact for Graphviz
+#@time @safetestset "Graphs" begin include("graphs.jl") end
+
 end # @time