up

Author: TorkelE

src/reactionsystem_conversions.jl

@@ -42,19 +42,6 @@ function oderatelaw(rx; combinatoric_ratelaw = true)
     rl
 end
 
-function assemble_drift(rs, ispcs; combinatoric_ratelaws = true, as_odes = true,
-                        include_zero_odes = true, remove_conserved = false)
-    rhsvec = assemble_oderhs(rs, ispcs; combinatoric_ratelaws, remove_conserved)
-    if as_odes
-        D = Differential(get_iv(rs))
-        eqs = [Equation(D(x), rhs)
-               for (x, rhs) in zip(ispcs, rhsvec) if (include_zero_odes || (!_iszero(rhs)))]
-    else
-        eqs = [Equation(0, rhs) for rhs in rhsvec if (include_zero_odes || (!_iszero(rhs)))]
-    end
-    eqs
-end
-
 function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserved = false)
     nps = get_networkproperties(rs)
     species_to_idx = Dict(x => i for (i, x) in enumerate(ispcs))
@@ -98,6 +85,19 @@ function assemble_oderhs(rs, ispcs; combinatoric_ratelaws = true, remove_conserv
     rhsvec
 end
 
+function assemble_drift(rs, ispcs; combinatoric_ratelaws = true, as_odes = true,
+                        include_zero_odes = true, remove_conserved = false)
+    rhsvec = assemble_oderhs(rs, ispcs; combinatoric_ratelaws, remove_conserved)
+    if as_odes
+        D = Differential(get_iv(rs))
+        eqs = [Equation(D(x), rhs)
+               for (x, rhs) in zip(ispcs, rhsvec) if (include_zero_odes || (!_iszero(rhs)))]
+    else
+        eqs = [Equation(0, rhs) for rhs in rhsvec if (include_zero_odes || (!_iszero(rhs)))]
+    end
+    eqs
+end
+
 # this doesn't work with constraint equations currently
 function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true,
                             remove_conserved = false)
@@ -258,6 +258,57 @@ function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
     return true
 end
 
+@inline function makemajump(rx; combinatoric_ratelaw = true)
+    @unpack rate, substrates, substoich, netstoich = rx
+    zeroorder = (length(substoich) == 0)
+    reactant_stoch = Vector{Pair{Any, eltype(substoich)}}()
+    @inbounds for (i, spec) in enumerate(substrates)
+        # move constant species into the rate
+        if isconstant(spec)
+            rate *= spec
+            isone(substoich[i]) && continue
+            for i in 1:(substoich[i] - 1)
+                rate *= spec - i
+            end
+        else
+            push!(reactant_stoch, substrates[i] => substoich[i])
+        end
+    end
+
+    if (!zeroorder) && combinatoric_ratelaw
+        coef = prod(factorial, substoich)
+        (!isone(coef)) && (rate /= coef)
+    end
+
+    net_stoch = filter(p -> !drop_dynamics(p[1]), netstoich)
+    isempty(net_stoch) &&
+        error("$rx has no net stoichiometry change once accounting for constant and boundary condition species. This is not supported.")
+
+    MassActionJump(Num(rate), reactant_stoch, net_stoch, scale_rates = false,
+                   useiszero = false)
+end
+
+# recursively visit each neighbor's rooted tree and mark everything in it as vrj
+function dfs_mark!(isvrjvec, visited, depgraph, i)
+    visited[i] = true
+    nhbrs = depgraph[i]
+    for nhbr in nhbrs
+        if !visited[nhbr]
+            isvrjvec[nhbr] = true
+            dfs_mark!(isvrjvec, visited, depgraph, nhbr)
+        end
+    end
+    nothing
+end
+
+# get_depgraph(rs)[i] is the list of reactions with rates depending on species changed by
+# i'th reaction.
+function get_depgraph(rs)
+    jdeps = asgraph(rs)
+    vdeps = variable_dependencies(rs)
+    eqeq_dependencies(jdeps, vdeps).fadjlist
+end
+
 function assemble_jumps(rs; combinatoric_ratelaws = true)
     meqs = MassActionJump[]
     ceqs = ConstantRateJump[]
@@ -323,57 +374,6 @@ function assemble_jumps(rs; combinatoric_ratelaws = true)
     vcat(meqs, ceqs, veqs)
 end
 
-# get_depgraph(rs)[i] is the list of reactions with rates depending on species changed by
-# i'th reaction.
-function get_depgraph(rs)
-    jdeps = asgraph(rs)
-    vdeps = variable_dependencies(rs)
-    eqeq_dependencies(jdeps, vdeps).fadjlist
-end
-
-# recursively visit each neighbor's rooted tree and mark everything in it as vrj
-function dfs_mark!(isvrjvec, visited, depgraph, i)
-    visited[i] = true
-    nhbrs = depgraph[i]
-    for nhbr in nhbrs
-        if !visited[nhbr]
-            isvrjvec[nhbr] = true
-            dfs_mark!(isvrjvec, visited, depgraph, nhbr)
-        end
-    end
-    nothing
-end
-
-@inline function makemajump(rx; combinatoric_ratelaw = true)
-    @unpack rate, substrates, substoich, netstoich = rx
-    zeroorder = (length(substoich) == 0)
-    reactant_stoch = Vector{Pair{Any, eltype(substoich)}}()
-    @inbounds for (i, spec) in enumerate(substrates)
-        # move constant species into the rate
-        if isconstant(spec)
-            rate *= spec
-            isone(substoich[i]) && continue
-            for i in 1:(substoich[i] - 1)
-                rate *= spec - i
-            end
-        else
-            push!(reactant_stoch, substrates[i] => substoich[i])
-        end
-    end
-
-    if (!zeroorder) && combinatoric_ratelaw
-        coef = prod(factorial, substoich)
-        (!isone(coef)) && (rate /= coef)
-    end
-
-    net_stoch = filter(p -> !drop_dynamics(p[1]), netstoich)
-    isempty(net_stoch) &&
-        error("$rx has no net stoichiometry change once accounting for constant and boundary condition species. This is not supported.")
-
-    MassActionJump(Num(rate), reactant_stoch, net_stoch, scale_rates = false,
-                   useiszero = false)
-end
-
 
 ### Equation Coupling ###
 

src/reactionsystem_structure.jl

@@ -794,6 +794,13 @@ function reactionparamsmap(network)
     Dict(p => i for (i, p) in enumerate(reactionparams(network)))
 end
 
+# used in the next function (`reactions(network)`).
+function namespace_reactions(network::ReactionSystem)
+    rxs = reactions(network)
+    isempty(rxs) && return Reaction[]
+    map(rx -> namespace_equation(rx, network), rxs)
+end
+
 """
     reactions(network)
 
@@ -809,11 +816,6 @@ function reactions(network)
     [rxs; reduce(vcat, namespace_reactions.(systems); init = Reaction[])]
 end
 
-function namespace_reactions(network::ReactionSystem)
-    rxs = reactions(network)
-    isempty(rxs) && return Reaction[]
-    map(rx -> namespace_equation(rx, network), rxs)
-end
 
 """
     numreactions(network)
@@ -1267,12 +1269,6 @@ function make_empty_network(; iv = DEFAULT_IV, name = gensym(:ReactionSystem))
 end
 
 # A helper function used in `flatten`.
-function getsubsystypes(sys)
-    typeset = Set{Type}()
-    getsubsystypes!(typeset, sys)
-    typeset
-end
-
 function getsubsystypes!(typeset::Set{Type}, sys::T) where {T <: MT.AbstractSystem}
     push!(typeset, T)
     for subsys in get_systems(sys)
@@ -1281,6 +1277,12 @@ function getsubsystypes!(typeset::Set{Type}, sys::T) where {T <: MT.AbstractSyst
     typeset
 end
 
+function getsubsystypes(sys)
+    typeset = Set{Type}()
+    getsubsystypes!(typeset, sys)
+    typeset
+end
+
 """
     Catalyst.flatten(rs::ReactionSystem)
 
@@ -1433,4 +1435,4 @@ function validate(rs::ReactionSystem, info::String = "")
 end
 
 # Checks if a unit consist of exponents with base 1 (and is this unitless).
-unitless_exp(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)
+(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)