update after Sam's comments

Author: TorkelE

src/dsl.jl

@@ -868,4 +868,31 @@ function get_reaction(line)
         error("Malformed reaction. @reaction macro only creates a single reaction. E.g. double arrows, such as `<-->` are not supported.")
     end
     return reaction[1]
+end
+
+
+### Generic Expression Manipulation ###
+
+# Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
+function recursive_expand_functions!(expr::ExprValues)
+    (typeof(expr) != Expr) && (return expr)
+    foreach(i -> expr.args[i] = recursive_expand_functions!(expr.args[i]),
+            1:length(expr.args))
+    if expr.head == :call
+        !isdefined(Catalyst, expr.args[1]) && (expr.args[1] = esc(expr.args[1]))
+    end
+    expr
+end
+
+# Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical).
+function tup_leng(ex::ExprValues)
+    (typeof(ex) == Expr && ex.head == :tuple) && (return length(ex.args))
+    return 1
+end
+
+# Gets the ith element in a expression tuple, or returns the input itself if it is not an expression tuple
+# (probably a  Symbol/Numerical).
+function get_tup_arg(ex::ExprValues, i::Int)
+    (tup_leng(ex) == 1) && (return ex)
+    return ex.args[i]
 end

src/expression_utils.jl

@@ -20,33 +20,6 @@ function is_escaped_expr(expr)
 end
 
 
-### Generic Expression Manipulation ###
-
-# Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
-function recursive_expand_functions!(expr::ExprValues)
-    (typeof(expr) != Expr) && (return expr)
-    foreach(i -> expr.args[i] = recursive_expand_functions!(expr.args[i]),
-            1:length(expr.args))
-    if expr.head == :call
-        !isdefined(Catalyst, expr.args[1]) && (expr.args[1] = esc(expr.args[1]))
-    end
-    expr
-end
-
-# Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical).
-function tup_leng(ex::ExprValues)
-    (typeof(ex) == Expr && ex.head == :tuple) && (return length(ex.args))
-    return 1
-end
-
-# Gets the ith element in a expression tuple, or returns the input itself if it is not an expression tuple
-# (probably a  Symbol/Numerical).
-function get_tup_arg(ex::ExprValues, i::Int)
-    (tup_leng(ex) == 1) && (return ex)
-    return ex.args[i]
-end
-
-
 ### Parameters/Species/Variables Symbols Correctness Checking ###
 
 # Throws an error when a forbidden symbol is used.

src/reaction_structure.jl renamed to src/reaction.jl

@@ -242,17 +242,9 @@ end
 # Union type for `Reaction`s and `Equation`s.
 const CatalystEqType = Union{Reaction, Equation}
 
-### Base Functions ###
-
-# Show function for `Reaction`s.
-function Base.show(io::IO, rx::Reaction)
-    print(io, rx.rate, ", ")
-    print_rxside(io, rx.substrates, rx.substoich)
-    arrow = rx.only_use_rate ? "⇒" : "-->"
-    print(io, " ", arrow, " ")
-    print_rxside(io, rx.products, rx.prodstoich)
-end
+### Base Function Dispatches ###
 
+# Used by `Base.show`.
 function print_rxside(io::IO, specs, stoich)
     # reactants/substrates
     if isempty(specs)
@@ -275,6 +267,15 @@ function print_rxside(io::IO, specs, stoich)
     nothing
 end
 
+# Show function for `Reaction`s.
+function Base.show(io::IO, rx::Reaction)
+    print(io, rx.rate, ", ")
+    print_rxside(io, rx.substrates, rx.substoich)
+    arrow = rx.only_use_rate ? "⇒" : "-->"
+    print(io, " ", arrow, " ")
+    print_rxside(io, rx.products, rx.prodstoich)
+end
+
 """
     ==(rx1::Reaction, rx2::Reaction)
 
@@ -311,10 +312,18 @@ function hash(rx::Reaction, h::UInt)
 end
 
 
-### ModelingToolkit-inherited Functions ###
+### ModelingToolkit Function Dispatches ###
+
+# Used by ModelingToolkit.namespace_equation.
+function apply_if_nonempty(f, v)
+    isempty(v) && return v
+    s = similar(v)
+    map!(f, s, v)
+    s
+end
 
 # Returns a name-spaced version of a reaction.
-function ModelingToolkit.namespace_equation(rx::Reaction, name; kw...)
+function MT.namespace_equation(rx::Reaction, name; kw...)
     f = Base.Fix2(namespace_expr, name)
     rate = f(rx.rate)
     subs = apply_if_nonempty(f, rx.substrates)
@@ -330,22 +339,15 @@ function ModelingToolkit.namespace_equation(rx::Reaction, name; kw...)
     Reaction(rate, subs, prods, substoich, prodstoich, netstoich, rx.only_use_rate, rx.metadata)
 end
 
-function apply_if_nonempty(f, v)
-    isempty(v) && return v
-    s = similar(v)
-    map!(f, s, v)
-    s
-end
-
 # Overwrites equation-type functions to give the correct input for `Reaction`s.
-ModelingToolkit.is_diff_equation(rx::Reaction) = false
-ModelingToolkit.is_alg_equation(rx::Reaction) = false
+MT.is_diff_equation(rx::Reaction) = false
+MT.is_alg_equation(rx::Reaction) = false
 
 
 ### Dependency-related Functions ###
 
 # determine which unknowns a reaction depends on
-function ModelingToolkit.get_variables!(deps::Set, rx::Reaction, variables)
+function MT.get_variables!(deps::Set, rx::Reaction, variables)
     (rx.rate isa Symbolic) && get_variables!(deps, rx.rate, variables)
     for s in rx.substrates
         # parametric stoichiometry means may have a parameter as a substrate
@@ -355,14 +357,14 @@ function ModelingToolkit.get_variables!(deps::Set, rx::Reaction, variables)
 end
 
 # determine which species a reaction modifies
-function ModelingToolkit.modified_unknowns!(munknowns, rx::Reaction, sts::Set)
+function MT.modified_unknowns!(munknowns, rx::Reaction, sts::Set)
     for (species, stoich) in rx.netstoich
         (species in sts) && push!(munknowns, species)
     end
     munknowns
 end
 
-function ModelingToolkit.modified_unknowns!(munknowns, rx::Reaction, sts::AbstractVector)
+function MT.modified_unknowns!(munknowns, rx::Reaction, sts::AbstractVector)
     for (species, stoich) in rx.netstoich
         any(isequal(species), sts) && push!(munknowns, species)
     end

src/reactionsystem_structure.jl renamed to src/reactionsystem.jl

@@ -523,7 +523,7 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; spa
 end
 
 
-### Base Functions ###
+### Base Function Dispatches ###
 
 """
     ==(rn1::ReactionSystem, rn2::ReactionSystem)
@@ -574,38 +574,6 @@ function isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = tr
 end
 
 
-### ModelingToolkit-inherited Functions ###
-
-# Retrives events.
-MT.get_continuous_events(sys::ReactionSystem) = getfield(sys, :continuous_events)
-# `MT.get_discrete_events(sys::ReactionSystem) = getfield(sys, :get_discrete_events)` should be added here.
-
-# need a custom equations since ReactionSystem.eqs are a mix of Reactions and Equations
-function MT.equations(sys::ReactionSystem)
-    ivs = independent_variables(sys)
-    eqs = get_eqs(sys)
-    systems = get_systems(sys)
-    if !isempty(systems)
-        eqs = CatalystEqType[eqs;
-                             reduce(vcat, MT.namespace_equations.(systems, (ivs,));
-                                    init = Any[])]
-        return sort!(eqs; by = eqsortby)
-    end
-    return eqs
-end
-
-function MT.unknowns(sys::ReactionSystem)
-    sts = get_unknowns(sys)
-    systems = get_systems(sys)
-    if !isempty(systems)
-        sts = unique!([sts; reduce(vcat, namespace_variables.(systems))])
-        sort!(sts; by = !isspecies)
-        return sts
-    end
-    return sts
-end
-
-
 ### Basic `ReactionSystem`-specific Accessors ###
 
 """
@@ -663,6 +631,24 @@ function combinatoric_ratelaws(sys::ReactionSystem)
     mapreduce(combinatoric_ratelaws, |, subsys; init = crl)
 end
 
+# Gets sub systems that are also reaction systems.
+# Used by several subsequent API functions.
+function filter_nonrxsys(network)
+    systems = get_systems(network)
+    rxsystems = ReactionSystem[]
+    for sys in systems
+        (sys isa ReactionSystem) && push!(rxsystems, sys)
+    end
+    rxsystems
+end
+
+# Special species but which take a set of states and names spaces them according to another
+# `ReactionSystem`.
+# Used by `species(network)`.
+function species(network::ReactionSystem, sts)
+    [MT.renamespace(network, st) for st in sts]
+end
+
 """
     species(network)
 
@@ -681,12 +667,6 @@ function species(network)
     unique([sts; reduce(vcat, map(sys -> species(sys, species(sys)), systems))])
 end
 
-# Special species but which take a set of states and names spaces them according to another
-# `ReactionSystem`.
-function species(network::ReactionSystem, sts)
-    [MT.renamespace(network, st) for st in sts]
-end
-
 """
     numspecies(network)
 
@@ -863,14 +843,36 @@ Returns whether `rn` has any spatial independent variables (i.e. is a spatial ne
 """
 isspatial(rn::ReactionSystem) = !isempty(get_sivs(rn))
 
-# Gets sub systems that are also reaction systems.
-function filter_nonrxsys(network)
-    systems = get_systems(network)
-    rxsystems = ReactionSystem[]
-    for sys in systems
-        (sys isa ReactionSystem) && push!(rxsystems, sys)
+
+### ModelingToolkit Function Dispatches ###
+
+# Retrives events.
+MT.get_continuous_events(sys::ReactionSystem) = getfield(sys, :continuous_events)
+# `MT.get_discrete_events(sys::ReactionSystem) = getfield(sys, :get_discrete_events)` should be added here.
+
+# need a custom equations since ReactionSystem.eqs are a mix of Reactions and Equations
+function MT.equations(sys::ReactionSystem)
+    ivs = independent_variables(sys)
+    eqs = get_eqs(sys)
+    systems = get_systems(sys)
+    if !isempty(systems)
+        eqs = CatalystEqType[eqs;
+                             reduce(vcat, MT.namespace_equations.(systems, (ivs,));
+                                    init = Any[])]
+        return sort!(eqs; by = eqsortby)
     end
-    rxsystems
+    return eqs
+end
+
+function MT.unknowns(sys::ReactionSystem)
+    sts = get_unknowns(sys)
+    systems = get_systems(sys)
+    if !isempty(systems)
+        sts = unique!([sts; reduce(vcat, namespace_variables.(systems))])
+        sort!(sts; by = !isspecies)
+        return sts
+    end
+    return sts
 end
 
 
@@ -1435,4 +1437,4 @@ function validate(rs::ReactionSystem, info::String = "")
 end
 
 # Checks if a unit consist of exponents with base 1 (and is this unitless).
-(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)
+unitless_exp(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)

src/reactionsystem_conversions.jl

@@ -42,6 +42,10 @@ function oderatelaw(rx; combinatoric_ratelaw = true)
     rl
 end
 
+# Function returning `true` for species which shouldn't change from the reactions, 
+# including non-species variables.
+drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s))
+
 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))
@@ -432,10 +436,6 @@ end
 
 ### Utility ###
 
-# Function returning `true` for species which shouldn't change from the reactions, 
-# including non-species variables.
-drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s))
-
 # Throws an error when attempting to convert a spatial system to an unssuported type.
 function spatial_convert_err(rs::ReactionSystem, systype)
     isspatial(rs) && error("Conversion to $systype is not supported for spatial networks.")
@@ -789,6 +789,24 @@ end
 
 ### Symbolic Variable/Symbol Conversions ###
 
+# convert symbol of the form :sys.a.b.c to a symbolic a.b.c
+function _symbol_to_var(sys, sym)
+    if hasproperty(sys, sym)
+        var = getproperty(sys, sym, namespace = false)
+    else
+        strs = split(String(sym), "₊")   # need to check if this should be split of not!!!
+        if length(strs) > 1
+            var = getproperty(sys, Symbol(strs[1]), namespace = false)
+            for str in view(strs, 2:length(strs))
+                var = getproperty(var, Symbol(str), namespace = true)
+            end
+        else
+            throw(ArgumentError("System $(nameof(sys)): variable $sym does not exist"))
+        end
+    end
+    var
+end
+
 """
     symmap_to_varmap(sys, symmap)
 
@@ -853,24 +871,6 @@ end
 symmap_to_varmap(sys, symmap) = symmap
 #error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.")
 
-# convert symbol of the form :sys.a.b.c to a symbolic a.b.c
-function _symbol_to_var(sys, sym)
-    if hasproperty(sys, sym)
-        var = getproperty(sys, sym, namespace = false)
-    else
-        strs = split(String(sym), "₊")   # need to check if this should be split of not!!!
-        if length(strs) > 1
-            var = getproperty(sys, Symbol(strs[1]), namespace = false)
-            for str in view(strs, 2:length(strs))
-                var = getproperty(var, Symbol(str), namespace = true)
-            end
-        else
-            throw(ArgumentError("System $(nameof(sys)): variable $sym does not exist"))
-        end
-    end
-    var
-end
-
 
 ### Other Conversion-related Functions ###