Merge pull request #117 from isaacsas/better-ismassaction-test

Better ismassaction test and == test for networks

Author: isaacsas

src/DiffEqBiological.jl

@@ -12,6 +12,8 @@ using Parameters
 @reexport using DiffEqBase, DiffEqJump
 using Compat
 
+import Base: (==)
+
 const ExprValues = Union{Expr,Symbol,Float64,Int}                   
 
 include("reaction_network.jl")

src/maketype.jl

@@ -133,6 +133,35 @@ function gentypefun_exprs(name; esc_exprs=true, gen_inplace=true, gen_outofplace
     exprs
 end
 
+######################## reaction network operators #######################
+
+# note, this is allocating, and could potentially be made more efficient by iterating
+# through reaction components directly...
+"""
+    ==(rn1::DiffEqBase.AbstractReactionNetwork, rn2::DiffEqBase.AbstractReactionNetwork)
+
+Tests whether the underlying species symbols, parameter symbols and reactions
+are the same in the two networks. Ignores order network components were defined,
+so the integer id of any individual species/parameters/reactions may be
+different between the two networks. *Does not* currently account for different
+reaction definitions, so "k, X+Y --> Y + Z" will not be the same as "k, Y+X -->
+Y + Z"
+"""
+function (==)(rn1::DiffEqBase.AbstractReactionNetwork, rn2::DiffEqBase.AbstractReactionNetwork)
+    issetequal(species(rn1), species(rn2)) || return false
+    issetequal(params(rn1), params(rn2)) || return false
+    nr1 = numreactions(rn1)
+    nr2 = numreactions(rn2)
+    (nr1 == nr2) || return false
+    idx1 = 1:nr1
+    idx2 = 1:nr2
+    issetequal(substrates.(rn1, idx1), substrates.(rn2, idx2)) || return false
+    issetequal(products.(rn1, idx1), products.(rn2, idx2)) || return false
+    issetequal(dependants.(rn1, idx1), dependants.(rn2, idx2)) || return false
+    issetequal(rateexpr.(rn1, idx1), rateexpr.(rn2, idx2)) || return false
+    issetequal(ismassaction.(rn1, idx1), ismassaction.(rn2, idx2)) || return false
+end
+
 ######################## functions to extend a network ####################
 
 """
@@ -209,88 +238,88 @@ the noise scaling coefficient.
 add_scale_noise_param!(rn::DiffEqBase.AbstractReactionNetwork, scale_noise_name::String) = add_scale_noise_param!(rn, Symbol(scale_noise_name))
 
 """
-    addreaction!(network, rateexpr::Union{Expr,Symbol,Int,Float64}, rxexpr::Expr)
+    addreaction!(network, rateex::Union{Expr,Symbol,Int,Float64}, rxexpr::Expr)
 
 Given an AbstractReaction network, add a reaction with the passed in rate and
 reaction expressions. i.e. a reaction of the form
 ```julia
 k*X, 2X + Y --> 2W
 ```
-would have `rateexpr=:(k*X)` and `rxexpr=:(2X + Y --> W)`, 
+would have `rateex=:(k*X)` and `rxexpr=:(2X + Y --> W)`, 
 ```julia
 10.5, 0 --> X
 ```
-would have `rateexpr=10.5` and `rxexpr=:(0 --> X)`, and
+would have `rateex=10.5` and `rxexpr=:(0 --> X)`, and
 ```julia
 k, X+X --> Z
 ```
-would have `rateexpr=:k` and `rxexpr=:(X+X --> Z)`.
+would have `rateex=:k` and `rxexpr=:(X+X --> Z)`.
 
 All normal DSL reaction definition notation should be supported.
 """
-function addreaction!(rn::DiffEqBase.AbstractReactionNetwork, rateexpr::ExprValues, rxexpr::Expr)
-    ex = Expr(:block, :(($rateexpr, $rxexpr)))
+function addreaction!(rn::DiffEqBase.AbstractReactionNetwork, rateex::ExprValues, rxexpr::Expr)
+    ex = Expr(:block, :(($rateex, $rxexpr)))
     newrxs = get_reactions(ex)
     foreach(rx -> push!(rn.reactions,ReactionStruct(rx, species(rn))), newrxs)
     nothing
 end
 
 """
-    addreaction!(network, rateexpr::Union{Expr,Symbol,Int,Float64}, substrates, products)
+    addreaction!(network, rateex::Union{Expr,Symbol,Int,Float64}, substrates, products)
 
 Given an AbstractReaction network, add a reaction with the passed in rate,
-`rateexpr`, substrate stoichiometry, and product stoichiometry. Stoichiometries
+`rateex`, substrate stoichiometry, and product stoichiometry. Stoichiometries
 are represented as tuples of `Pair{Symbol,Int}`. i.e. a reaction of the form
 ```julia
 k*X, 2X + Y --> 2W
 ```
-would have `rateexpr=:(k*X)`, `substrates=(:X=>2, :Y=>2)`` and
+would have `rateex=:(k*X)`, `substrates=(:X=>2, :Y=>2)`` and
 `products=(W=>2,)`, 
 ```julia
 10.5, 0 --> X
 ```
-would have `rateexpr=10.5`, `substrates=()` and `products=(:X=>1,)`, and
+would have `rateex=10.5`, `substrates=()` and `products=(:X=>1,)`, and
 ```julia
 k, X+X --> Z
 ```
-would have `rateexpr=:k`, `substrates=(:X=>2,)` and `products=(:Z=>2,)`.
+would have `rateex=:k`, `substrates=(:X=>2,)` and `products=(:Z=>2,)`.
 
 All normal DSL reaction definition notation should be supported for the
-`rateexpr`.
+`rateex`.
 """
-function addreaction!(rn::DiffEqBase.AbstractReactionNetwork, rateexpr::ExprValues, 
+function addreaction!(rn::DiffEqBase.AbstractReactionNetwork, rateex::ExprValues, 
                                         subs::Tuple{Vararg{Pair{Symbol,Int}}}, 
                                         prods::Tuple{Vararg{Pair{Symbol,Int}}}) where {T <: Number}
 
     substrates = ReactantStruct[ReactantStruct(p[1],p[2]) for p in subs]
     dependents = Symbol[p[1] for p in subs]
     products = ReactantStruct[ReactantStruct(p[1],p[2]) for p in prods]
-    rate_DE = mass_rate_DE(substrates, true, rateexpr)
-    rate_SSA = mass_rate_SSA(substrates, true, rateexpr)   
+    rate_DE = mass_rate_DE(substrates, true, rateex)
+    rate_SSA = mass_rate_SSA(substrates, true, rateex)   
 
-    # resolve dependents from rateexpr
-    if rateexpr isa Number
+    # resolve dependents from rateex
+    if rateex isa Number
         ismassaction = true        
-    elseif rateexpr isa Symbol
-        if haskey(speciesmap(rn), rateexpr)
+    elseif rateex isa Symbol
+        if haskey(speciesmap(rn), rateex)
             ismassaction = false
-            if rateexpr ∉ dependents
-                push!(dependents, rateexpr)
+            if rateex ∉ dependents
+                push!(dependents, rateex)
             end
-        elseif haskey(paramsmap(rn), rateexpr)
+        elseif haskey(paramsmap(rn), rateex)
             ismassaction = true
         else
-            error("rateexpr is a symbol that is neither a species or parameter.")
+            error("rateex is a symbol that is neither a species or parameter.")
         end
     else # isa Expr
 
         # mimicing ReactionStruct constructor for now, but this should be optimized...
         newdeps = unique!(recursive_content(rate_DE, speciesmap(rn), Vector{Symbol}()))
-        ismassaction = length(newdeps)==length(intersect!(dependents, newdeps)) ? true : false
+        ismassaction = issetequal(dependents,newdeps) 
         dependents = newdeps
     end
 
-    push!(rn.reactions, ReactionStruct(substrates, products, rateexpr, rate_DE, rate_SSA, dependents, ismassaction))
+    push!(rn.reactions, ReactionStruct(substrates, products, rateex, rate_DE, rate_SSA, dependents, ismassaction))
     nothing
 end
 

src/network_properties.jl

@@ -313,11 +313,25 @@ end
 Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
 vector of pairs, mapping ids of species that serve as substrates in the reaction
 to the corresponding stoichiometric coefficient as a substrate. 
+
+Allocates a new vector to store the pairs.
 """ 
 function substratestoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
     substratestoich(rn.reactions[rxidx], speciesmap(rn))
 end
 
+"""
+    substratesymstoich(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+`ReactantStruct`, mapping the symbols of species that serve as substrates in the
+reaction to the corresponding stoichiometric coefficient as a substrate. 
+
+Non-allocating, returns underlying field within the reaction_network.
+""" 
+function substratesymstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    rn.reactions[rxidx].substrates
+end
 
 function productstoich(rs::ReactionStruct, specmap)
     sort( [specmap[p.reactant] => p.stoichiometry for p in rs.products] )    
@@ -329,11 +343,26 @@ end
 Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
 vector of pairs, mapping ids of species that are products in the reaction to the
 corresponding stoichiometric coefficient as a product.
+
+Allocates a new vector to store the pairs.
 """
 function productstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
     productstoich(rn.reactions[rxidx], speciesmap(rn))
 end
 
+"""
+    productsymstoich(network, rxidx)
+
+Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
+`ReactantStruct`, mapping the symbols of species that are products in the
+reaction to the corresponding stoichiometric coefficient as a product. 
+
+Non-allocating, returns underlying field within the reaction_network.
+""" 
+function productsymstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
+    rn.reactions[rxidx].products
+end
+
 function netstoich(rs::ReactionStruct, specmap)
     nsdict = Dict{Int,Int}(specmap[s.reactant] => -s.stoichiometry for s in rs.substrates)
 
@@ -357,6 +386,8 @@ Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
 vector of pairs, mapping ids of species that participate in the reaction to the
 net stoichiometric coefficient of the species (i.e. net change in the species
 due to the reaction).
+
+Allocates a new vector to store the pairs.
 """
 function netstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
     netstoich(rn.reactions[rxidx], speciesmap(rn))
@@ -376,6 +407,8 @@ would return true, while reactions with state-dependent rates like
 `k*X, X + Y --> Z`
 
 would return false.
+
+Non-allocating, returns underlying field within the reaction_network.
 """
 function ismassaction(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
     rn.reactions[rxidx].is_pure_mass_action
@@ -391,6 +424,8 @@ vector of symbols of species the *reaction rate law* depends on. i.e. for
 `k*W, 2X + 3Y --> 5Z + W`
 
 the returned vector would be `[:W,:X,:Y]`.
+
+Non-allocating, returns underlying field within the reaction_network.
 """
 function dependents(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
     rn.reactions[rxidx].dependants
@@ -415,9 +450,11 @@ i.e. for
 `k*W, X + 3Y --> X + W`
 
 the returned vector would be `[:X,:Y]`.
+
+Allocates a new vector to store the symbols.
 """
 function substrates(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
-    rn.reactions[rxidx].substrates
+    [sub.reactant for sub in rn.reactions[rxidx].substrates]
 end
 
 """
@@ -430,9 +467,11 @@ i.e. for
 `k*W, X + 3Y --> X + W`
 
 the returned vector would be `[:X,:W]`.
+
+Allocates a new vector to store the symbols.
 """
 function products(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
-    rn.reactions[rxidx].products
+    [prod.reactant for prod in rn.reactions[rxidx].products]
 end
 
 ######### Network Properties: #########