Merge pull request #319 from isaacsas/allow-sys-names-and-more

Allow sys names and other API fixes

Author: isaacsas

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "6.9.1"
+version = "6.10.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -17,7 +17,7 @@ DataStructures = "0.18"
 DocStringExtensions = "0.8"
 Latexify = "0.14"
 MacroTools = "0.5.5"
-ModelingToolkit = "5.12.0"
+ModelingToolkit = "5.13.0"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"

docs/src/api/catalyst_api.md

@@ -11,8 +11,9 @@ be converted to other `ModelingToolkit.AbstractSystem`s, including a
 `ModelingToolkit.ODESystem`, `ModelingToolkit.SDESystem`, or
 `ModelingToolkit.JumpSystem`.
 
-An empty network can be generated using [`@reaction_network`](@ref) with no arguments or
-the [`make_empty_network`](@ref) function. These can then be extended
+An empty network can be generated using [`@reaction_network`](@ref) with no
+arguments (or one argument to name the system), or the
+[`make_empty_network`](@ref) function. These can then be extended
 programmatically using [`addspecies!`](@ref), [`addparam!`](@ref), and
 [`addreaction!`](@ref).
 

docs/src/tutorials/generated_systems.md

@@ -33,7 +33,7 @@ Each `Reaction` within `reactions(rn)` has a number of subfields. For `rx` a
   non-filled arrows and should ignore mass action kinetics. `false` by default.
 
 Empty `ReactionSystem`s can be generated via [`make_empty_network`](@ref) or
-[`@reaction_network`](@ref) with no arguments. `ReactionSystem`s can be
-programmatically extended using [`addspecies!`](@ref), [`addparam!`](@ref),
-[`addreaction!`](@ref), [`@add_reactions`](@ref), or composed using `merge` and
-`merge!`.
+[`@reaction_network`](@ref) with no arguments (giving one argument to the latter
+will specify a system name). `ReactionSystem`s can be programmatically extended
+using [`addspecies!`](@ref), [`addparam!`](@ref), [`addreaction!`](@ref),
+[`@add_reactions`](@ref), or composed using `merge` and `merge!`.

src/networkapi.jl

@@ -69,7 +69,7 @@ Return the number of species within the given [`ReactionSystem`](@ref).
 function numspecies(network)
     ns = length(ModelingToolkit.ModelingToolkit.get_states(network))
     for sys in ModelingToolkit.get_systems(network)
-        ns += numspecies(ns)
+        ns += numspecies(sys)
     end
     ns
 end
@@ -95,7 +95,7 @@ Return the number of parameters within the given [`ReactionSystem`](@ref).
 function numparams(network)
     np = length(ModelingToolkit.get_ps(network))
     for sys in ModelingToolkit.get_systems(network)
-        np += numparams(ns)
+        np += numparams(sys)
     end
     np
 end
@@ -219,13 +219,14 @@ end
 ######################## functions to extend a network ####################
 
 """
-    make_empty_network(; iv=Sym{ModelingToolkit.Parameter{Real}}(:t))
+    make_empty_network(; iv=Sym{ModelingToolkit.Parameter{Real}}(:t), 
+                         name=gensym(:ReactionSystem))
 
 Construct an empty [`ReactionSystem`](@ref). `iv` is the independent variable,
-usually time.
+usually time, and `name` is the name to give the `ReactionSystem`.
 """
-function make_empty_network(; iv=Sym{ModelingToolkit.Parameter{Real}}(:t))
-    ReactionSystem(Reaction[], iv, [], [], Equation[], gensym(:ReactionSystem), ReactionSystem[])
+function make_empty_network(; iv=Sym{ModelingToolkit.Parameter{Real}}(:t), name=gensym(:ReactionSystem))
+    ReactionSystem(Reaction[], iv, [], [], Equation[], name, ReactionSystem[])
 end
 
 """
@@ -244,10 +245,10 @@ Notes:
 function addspecies!(network::ReactionSystem, s::Symbolic; disablechecks=false)
 
     # we don't check subsystems since we will add it to the top-level system...
-    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, s), ModelingToolkit.ModelingToolkit.get_states(network))
+    curidx = disablechecks ? nothing : findfirst(S -> isequal(S, s), ModelingToolkit.get_states(network))
     if curidx === nothing
-        push!(ModelingToolkit.ModelingToolkit.get_states(network), s)
-        return length(ModelingToolkit.ModelingToolkit.get_states(network))
+        push!(ModelingToolkit.get_states(network), s)
+        return length(ModelingToolkit.get_states(network))
     else
         return curidx
     end

src/reaction_network.jl

@@ -65,6 +65,14 @@ bwd_arrows = Set{Symbol}([:<, :←, :↢, :↤, :⇽, :⟵, :⟻, :⥚, :⥞, :
 double_arrows = Set{Symbol}([:↔, :⟷, :⇄, :⇆, :⇌, :⇋, :⇔, :⟺])
 pure_rate_arrows = Set{Symbol}([:⇐, :⟽, :⇒, :⟾, :⇔, :⟺])
 
+# unfortunately the following doesn't seem to work on 1.5, so supporting these 
+# operators seems to require us to only support 1.6 and up...
+# @static if VERSION >= v"1.6.0"
+#     push!(bwd_arrows, :<--)
+#     push!(double_arrows, :<-->)
+# end
+
+
 # Declares symbols which may neither be used as parameters not varriables.
 forbidden_symbols = [:t, :π, :pi, :ℯ, :im, :nothing, :∅]
 
@@ -78,26 +86,50 @@ network.
 
 See the [Catalyst.jl for Reaction Models](@ref) documentation for details on
 parameters to the macro.
+
+Examples:
+```julia
+# a basic SIR model, with name SIR
+sir_model = @reaction_network SIR begin
+    c1, s + i --> 2i
+    c2, i --> r
+end c1 c2
+
+# a basic SIR model, with random generated name
+sir_model = @reaction_network begin
+    c1, s + i --> 2i
+    c2, i --> r
+end c1 c2
+
+# an empty network with name empty
+emptyrn = @reaction_network empty
+
+# an empty network with random generated name
+emptyrn = @reaction_network 
+```
 """
+macro reaction_network(name::Symbol, ex::Expr, parameters...)
+    make_reaction_system(MacroTools.striplines(ex), parameters; name=name)
+end
+
 macro reaction_network(ex::Expr, parameters...)
     make_reaction_system(MacroTools.striplines(ex), parameters)
 end
 
-### Macros used for manipulating, and successively builing up, reaction systems. ###
-
-#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 parameters: "*((map(s -> "'"*string(s)*"', ",intersect(forbidden_symbols,parameters))...))*"this is not permited.")
-    return Expr(:block,:(@parameters $((:t,parameters...)...)),
+#Returns a empty network (with, or without, a declared name)
+macro reaction_network(name::Symbol=gensym(:ReactionSystem))
+    return Expr(:block,:(@parameters t),
                 :(ReactionSystem(Reaction[],
                                  t,
                                  [],
-                                 [$(parameters...)], 
+                                 [], 
                                  Equation[],
-                                 gensym(:ReactionSystem),
+                                 $(QuoteNode(name)),
                                  ReactionSystem[])))
 end
 
+### Macros used for manipulating, and successively builing up, reaction systems. ###
+
 """
     @add_reactions
 
@@ -139,16 +171,15 @@ end
 ### Functions that process the input and rephrase it as a reaction system ###
 
 # Takes the reactions, and rephrases it as a "ReactionSystem" call, as designated by the ModelingToolkit IR.
-function make_reaction_system(ex::Expr, parameters)
+function make_reaction_system(ex::Expr, parameters; name=gensym(:ReactionSystem))
     reactions = get_reactions(ex)
     reactants = get_reactants(reactions)
-    !isempty(intersect(forbidden_symbols,union(reactants,parameters))) && error("The following symbol(s) are used as reactants or parameters: "*((map(s -> "'"*string(s)*"', ",intersect(forbidden_symbols,union(reactants,parameters)))...))*"this is not permited.")
-
-    network_code = Expr(:block,:(@parameters t),:(@variables), :(ReactionSystem([],t,[],[])))
+    !isempty(intersect(forbidden_symbols,union(reactants,parameters))) && error("The following symbol(s) are used as reactants or parameters: "*((map(s -> "'"*string(s)*"', ",intersect(forbidden_symbols,union(reactants,parameters)))...))*"this is not permited.")    
+    network_code = Expr(:block,:(@parameters t),:(@variables), :(ReactionSystem([],t,[],[]; name=$(QuoteNode(name)))))
     foreach(parameter-> push!(network_code.args[1].args, parameter), parameters)
     foreach(reactant -> push!(network_code.args[2].args, Expr(:call,reactant,:t)), reactants)
-    foreach(parameter-> push!(network_code.args[3].args[5].args, parameter), parameters)
-    foreach(reactant -> push!(network_code.args[3].args[4].args, reactant), reactants)
+    foreach(parameter-> push!(network_code.args[3].args[6].args, parameter), parameters)
+    foreach(reactant -> push!(network_code.args[3].args[5].args, reactant), reactants)
     for reaction in reactions
         subs_init = isempty(reaction.substrates) ? nothing : :([]); subs_stoich_init = deepcopy(subs_init)
         prod_init = isempty(reaction.products) ? nothing : :([]); prod_stoich_init = deepcopy(prod_init)
@@ -161,7 +192,7 @@ function make_reaction_system(ex::Expr, parameters)
             push!(reaction_func.args[4].args, prod.reactant)
             push!(reaction_func.args[6].args, prod.stoichiometry)
         end
-        push!(network_code.args[3].args[2].args,reaction_func)
+        push!(network_code.args[3].args[3].args,reaction_func)
     end
     return network_code
 end

test/make_model.jl

@@ -8,10 +8,8 @@ rng = StableRNG(12345)
 include("test_networks.jl")
 
 
-### Debugg functions ###
-
+### Debug functions ###
 opname(x) = istree(x) ? nameof(operation(x)) : nameof(x)
-
 alleq(xs,ys) = all(isequal(x,y) for (x, y) in zip(xs, ys))
 
 # Gets all the reactants in a set of equations.
@@ -405,3 +403,13 @@ rn = @reaction_network begin
     k2, I --> R
 end k1 k2
 @test isequal(opname(species(rn)[2]),:I)
+
+
+# test names work
+rn = @reaction_network SIR1 begin
+    k1, S + I --> 2I
+    k2, I --> R
+end k1 k2
+@test nameof(rn) == :SIR1
+
+

test/model_modification.jl

@@ -13,7 +13,13 @@ empty_network_1 = @reaction_network
 @test length(empty_network_1.states) == 0
 @test length(ps) == 0
 
-empty_network_2 = @reaction_network p1 p2 p3 p4 p5
+empty_network_2 = @reaction_network 
+@variables p1 p2 p3 p4 p5
+addparam!(empty_network_2, p1)
+addparam!(empty_network_2, p2)
+addparam!(empty_network_2, p3)
+addparam!(empty_network_2, p4)
+addparam!(empty_network_2, p5)
 @unpack eqs,iv,ps,name,systems = empty_network_2
 @test length(eqs) == 0
 @test iv.name == :t
@@ -129,7 +135,9 @@ end k1
 end k4 m J4
 push!(identical_networks, reaction_networks_real[3] => step_by_step_network_7)
 
-step_by_step_network_8 = @reaction_network k1
+step_by_step_network_8 = @reaction_network 
+@parameters k1
+addparam!(step_by_step_network_8,k1)
 @add_reactions step_by_step_network_8 begin
     k1, X1 → X2
     0, X2 → X3