@compounds in the macro

Author: TorkelE

src/chemistry_functionality.jl

@@ -39,9 +39,14 @@ function make_compound(expr)
     species_expr = expr.args[1]                                         # E.g. :(CO2(t))
     species_name = expr.args[1].args[1]                                 # E.g. :CO2
     composition = Catalyst.recursive_find_reactants!(expr.args[2].args[1], 1, Vector{ReactantStruct}(undef, 0))
-    components = :([])                                                  # Extra step required here to get something like :([C, O]), rather than :([:C, :O])
-    foreach(comp -> push!(components.args, comp.reactant), composition) # E.g. [C, O]
-    coefficients = getfield.(composition, :stoichiometry)               # E.g. [1, 2]
+
+    # Loops through all components, add the component and the coefficients to the corresponding vectors (cannot extract directly using e.g. "getfield.(composition, :reactant)" because then we get something like :([:C, :O]), rather than :([C, O]))
+    components = :([])                                                  # Becomes something like :([C, O]).                                         
+    coefficients = :([])                                                # Becomes something like :([1, 2]). 
+    for comp in composition
+        push!(components.args, comp.reactant)
+        push!(coefficients.args, comp.stoichiometry)
+    end
 
     # Creates the found expressions that will create the compound species.
     # The `Expr(:escape, :(...))` is required so that teh expressions are evaluated in the scope the users use the macro in (to e.g. detect already exiting species).
@@ -50,6 +55,14 @@ function make_compound(expr)
     components_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundComponents, $components)))        # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])`
     coefficients_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundCoefficients, $coefficients)))  # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])`
 
+    println(
+        return quote
+            $species_declaration_expr
+            $compound_designation_expr
+            $components_designation_expr
+            $coefficients_designation_expr
+        end)
+
     # Returns the rephrased expression.
     return quote
         $species_declaration_expr
@@ -85,9 +98,22 @@ end
 
 # Function managing the @compound macro.
 function make_compounds(expr)
-    # Creates a separate compound call for each compound.
-    compound_calls = [make_compound(line) for line in expr.args]                # Crates a vector containing the quote for each compound.
-    return Expr(:block, vcat([c_call.args for c_call in compound_calls]...)...) # Combines the quotes to a single one. Don't want the double "...". But had problems getting past them for various metaprogramming reasons :(.)
+    # Creates an empty block containing the output call.
+    compound_declarations = Expr(:block)
+
+    # Creates a compound creation set of lines (4 in total) for each line. Loops through all 4x[Number of compounds] liens and add them to compound_declarations.
+    compound_calls = [Catalyst.make_compound(line) for line in expr.args] 
+    for compound_call in compound_calls, line in MacroTools.striplines(compound_call).args
+        push!(compound_declarations.args, line)
+    end
+
+    # The output of the macros should be a vector with the compounds (same as for e.g. "@species begin ... end", also require for things to work in the DSL).
+    # Creates an output vector, and loops through all compounds, adding them to it.
+    push!(compound_declarations.args, :($(Expr(:escape, :([])))))
+    for arg in expr.args
+        push!(compound_declarations.args[end].args[1].args, arg.args[1].args[1])
+    end
+    return compound_declarations
 end
 
 ## Compound Getters ###

src/reaction_network.jl

@@ -123,7 +123,7 @@ const forbidden_variables_error = let
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -358,9 +358,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
                        option_lines))
 
+    # Reads compounds options.
+    compound_expr, compound_species = read_compound_options(options)
+
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
-    species_declared = extract_syms(options, :species)
+    species_declared = [extract_syms(options, :species); compound_species]
     parameters_declared = extract_syms(options, :parameters)
     variables = extract_syms(options, :variables)
 
@@ -399,6 +402,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
     ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps")
     vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars")
     sps, spssym = scalarize_macro(!isempty(species), sexprs, "specs")
+    comps, compssym = scalarize_macro(!isempty(compound_species), compound_expr, "comps")
     rxexprs = :(Catalyst.CatalystEqType[])
     for reaction in reactions
         push!(rxexprs.args, get_rxexprs(reaction))
@@ -410,8 +414,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem)))
         $ivexpr
         $vars
         $sps
+        $comps
 
-        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym),
+        Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym, $compssym),
                                               $pssym; name = $name,
                                               spatial_ivs = $sivs)
     end
@@ -461,6 +466,19 @@ function esc_dollars!(ex)
     ex
 end
 
+# When compound species are declared using the "@compound begin ... end" option, get a list of the compound species, and also the expression that crates them.
+function read_compound_options(opts)
+    # If the compound option is used retrive a list of compound species, and the option that creeates them
+    if haskey(opts, :compounds)
+        compound_expr = opts[:compounds]
+        compound_species = [arg.args[1].args[1] for arg in compound_expr.args[3].args] # Loops through where in the "@compound begin ... end" the compound species names are.
+    else
+        compound_expr = :()
+        compound_species = Union{Symbol, Expr}[]
+    end
+    return compound_expr, compound_species
+end
+
 # When the user have used the @species (or @parameters) options, extract species (or
 # parameters) from its input.
 function extract_syms(opts, vartype::Symbol)

test/miscellaneous_tests/compound_macro.jl

@@ -1,9 +1,9 @@
 using Catalyst, Test
 
-# Test base funcationality in two cases.
+# Test base functionality in two cases.
 let
     @variables t
-    @species C(t) H(t) + O(t)
+    @species C(t) H(t) O(t)
     @compound C6H12O2(t) = 6C + 12H + 2O
 
     @test iscompound(C6H12O2)
@@ -102,14 +102,23 @@ let
     @species H(t)
 
     alpha = 2
-    @compound H2_1(t) = alpha*H
-    @compound H2_2(t) = 2H
+    h = H
+    @compound H2_1(t) = 2*H
+    @compound H2_2(t) = alpha*H
+    @compound H2_3(t) = 2*h
+    @compound H2_4(t) = alpha*2H
 
     @test iscompound(H2_1)
     @test iscompound(H2_2)
+    @test iscompound(H2_2)
+    @test iscompound(H2_4)
 
     @test isequal(components(H2_1),components(H2_2))
+    @test isequal(components(H2_2),components(H2_3))
+    @test isequal(components(H2_3),components(H2_4))
     @test isequal(coefficients(H2_1),coefficients(H2_2))
+    @test isequal(coefficients(H2_2),coefficients(H2_3))
+    @test isequal(coefficients(H2_3),coefficients(H2_4))
 end
 
 let
@@ -185,4 +194,79 @@ let
     @test isequal(components(comp),components(comp_alt))
     @test isequal(coefficients(comp), coefficients(comp_alt))
     @test isequal(component_coefficients(comp), component_coefficients(comp_alt))
+end
+
+### Compounds in DSL ###
+
+# Checks with a single compound.
+# Checks using @unpack.
+# Check where compounds and components does not occur in reactions.
+let 
+    rn = @reaction_network begin
+        @species C(t) O(t)
+        @compounds begin
+            CO2(t) = C + 2O
+        end
+    end
+    @unpack C, O, CO2 = rn
+    
+    @test length(species(rn)) == 3
+    @test iscompound(CO2)
+    @test isequal([C, O], components(CO2))
+    @test isequal([1, 2], coefficients(CO2))
+    @test isequal([C => 1, O => 2], component_coefficients(CO2)) 
+end
+
+# Test using multiple compounds.
+# Test using rn. notation to fetch species.
+let 
+    rn = complete(@reaction_network begin
+        @species C(t) O(t) H(t)
+        @compounds begin
+            CH4(t) = C + 4H
+            O2(t) = 2O
+            CO2(t) = C + 2O
+            H2O(t) = 2H + O
+        end
+        k, CH4 + O2 --> CO2 + H2O
+    end)
+    species(rn)
+    
+    @test length(species(rn)) == 7
+    @test isequal([rn.C, rn.H], components(rn.CH4))
+    @test isequal([1, 4], coefficients(rn.CH4))
+    @test isequal([rn.C => 1, rn.H => 4], component_coefficients(rn.CH4)) 
+    @test isequal([rn.O], components(rn.O2))
+    @test isequal([2], coefficients(rn.O2))
+    @test isequal([rn.O => 2], component_coefficients(rn.O2)) 
+    @test isequal([rn.C, rn.O], components(rn.CO2))
+    @test isequal([1, 2], coefficients(rn.CO2))
+    @test isequal([rn.C => 1, rn.O => 2], component_coefficients(rn.CO2)) 
+    @test isequal([rn.H, rn.O], components(rn.H2O))
+    @test isequal([2, 1], coefficients(rn.H2O))
+    @test isequal([rn.H => 2, rn.O => 1], component_coefficients(rn.H2O)) 
+end
+
+# Tests using compounds of compounds.
+# Tests where species are part of reactions and not declared using "@species".
+let
+    rn = @reaction_network begin
+        @compounds begin
+            SO2(t) = S + 2O
+            S2O4(t) = 2SO2
+        end
+        dS, S --> 0
+        dO, O --> 0
+    end
+    species(rn)
+    @unpack S, O, SO2, S2O4 = rn
+    
+    @test length(species(rn)) == 4
+    
+    @test isequal([S, O], components(SO2))
+    @test isequal([1, 2], coefficients(SO2))
+    @test isequal([S => 1, O => 2], component_coefficients(SO2)) 
+    @test isequal([SO2], components(S2O4))
+    @test isequal([2], coefficients(S2O4))
+    @test isequal([SO2 => 2], component_coefficients(S2O4)) 
 end

test/miscellaneous_tests/reaction_balancing.jl

@@ -8,7 +8,7 @@ let
     @compounds begin
         H2(t) = 2H
         O2(t) = 2O
-        H2O(t) = 2H + 1O2
+        H2O(t) = 2H + 1O
     end
 
     rx = Reaction(k,[H2,O2],[H2O])
@@ -29,7 +29,7 @@ end
 let
     @variables t
     @parameters k
-    @species C(t) H(t) = O(t)
+    @species C(t) H(t) O(t)
     @compound O2(t) = 2O
     @compound CO2(t) = 1C + 2O
     @compound H2O(t) = 2H + 1O
@@ -416,3 +416,25 @@ let
     rx = Reaction(1.0, [CO, H2], [COH2])
     @test_throws Catalyst.COMPOUND_OF_COMPOUND_ERROR balance_reaction(rx)
 end
+
+# Checks that balancing work for a reaction from a reaction_network.
+let
+    rn = complete(@reaction_network begin
+        @species C(t) H(t) O(t)
+        @compounds begin
+            O2(t) = 2O
+            CO2(t) = 1C + 2O
+            H2O(t) = 2H + 1O
+            C6H12O6(t) = 6C + 12H + 6O
+        end
+        k, CO2 + H2O --> C6H12O6 + O2
+    end)
+    
+    brxs = balance_reaction(reactions(rn)[1])[1]
+    
+    @test isequal(rn.k, brxs.rate)
+    @test isequal([rn.CO2,  rn.H2O], brxs.substrates)
+    @test isequal([rn.C6H12O6,  rn.O2], brxs.products)
+    @test isequal([6, 6], brxs.substoich)
+    @test isequal([1,  6], brxs.prodstoich)    
+end