up

Author: TorkelE

docs/src/catalyst_functionality/chemistry_related_functionality.md

@@ -15,9 +15,9 @@ using Catalyst
 ```
 Next, we create the `CO2` compound species:
 ```@example chem1
-@compound CO2(t) = C + 2O
+@compound CO2 ~ C + 2O
 ```
-Here, the compound is the first argument to the macro, followed by its component. The `(t)` indicates that `CO2` is a time dependant species. Components with non-unitary stoichiometries have this value written before the component (generally, the rules for designating the components of a compound are identical to those of designating the substrates or products of a reaction). The created compound, `CO2`, is a species in every sense, and can be used wherever e.g. `C` can be used:
+Here, the compound is the first argument to the macro, followed by its component (with the left-hand and right-hand sides separated by a `~` sign). While non-compound species (such as `C` and `O`) have their independent variable (in this case `t`) designated, independent variables as not designated for compounds (but instead directly inferred from their components). Components with non-unitary stoichiometries have this value written before the component (generally, the rules for designating the components of a compound are identical to those of designating the substrates or products of a reaction). The created compound, `CO2`, is a species in every sense, and can be used wherever e.g. `C` can be used:
 ```@example chem1
 isspecies(CO2)
 ```
@@ -40,16 +40,16 @@ iscompound(CO2)
 Compound components that are also compounds are allowed, e.g. we can create a carbonic acid compound (H₂CO₃) that consists of CO₂ and H₂O:
 ```@example chem1
 @species H(t)
-@compound H2O(t) = 2H + O
-@compound H2CO3(t) = CO2 + H2O
+@compound H2O ~ 2H + O
+@compound H2CO3 ~ CO2 + H2O
 ```
 
 When multiple compounds are created, they can be created simultaneously using the `@compounds` macro, e.g. the previous code-block can be re-written as:
 ```@example chem1
 @species H(t)
 @compounds begin
-    H2O(t) = 2H + O
-    H2CO3(t) = CO2 + H2O
+    H2O ~ 2H + O
+    H2CO3 ~ CO2 + H2O
 end
 ```
 
@@ -59,9 +59,9 @@ It is also possible to declare species as compound species within the `@reaction
 rn = @reaction_network begin
     @species C(t) H(t) O(t)
     @compounds begin
-        C2O(t) = C + 2O
-        H2O(t) = 2H + O
-        H2CO3(t) = CO2 + H2O
+        C2O ~ C + 2O
+        H2O ~ 2H + O
+        H2CO3 ~ CO2 + H2O
     end
     (k1,k2), H2O+ CO2 <--> H2CO3
 end
@@ -70,15 +70,29 @@ When creating compound species using the DSL, it is important to note that *ever
 ```julia 
 rn = @reaction_network begin
     @compounds begin
-        C2O(t) = C + 2O
-        H2O(t) = 2H + O
-        H2CO3(t) = CO2 + H2O
+        C2O ~ C + 2O
+        H2O ~ 2H + O
+        H2CO3 ~ CO2 + H2O
     end
     (k1,k2), H2O+ CO2 <--> H2CO3
 end
 ```
 as the components `C`, `H`, and `O` are not declared as a species anywhere. Please also note that only `@compounds` can be used as an option in the DSL, not `@compound`.
 
+#### Designating metadata and default values for compounds
+Just like for normal species, it is possible to designate metadata and default values for compounds. Metadata is provided after the compound name, but separated from it by a `,`:
+```@example chem1
+@compound CO2, [unit="mol"] ~ C + 2O
+```
+Default values are designated using `=`, and provided directly after the compound name. If default values are given, the left-hand side must be grouped using `()`:
+```@example chem1
+@compound (CO2 = 2.0) ~ C + 2O
+```
+If both default values and meta data are provided, the metadata is provided after teh default value:
+```@example chem1
+@compound (CO2 = 2.0, [unit="mol"]) ~ C + 2O
+```
+
 ## Balancing chemical reactions
 One use of defining a species as a compound is that they can be used to balance reactions to that the number of components are the same on both sides. Catalyst provides the `balance_reaction` function, which takes a reaction, and returns a balanced version. E.g. let us consider a reaction when carbon dioxide is formed from carbon and oxide `C + O --> CO2`. Here, `balance_reaction` enables us to find coefficients creating a balanced reaction (in this case, where the number of carbon and oxygen atoms are the same on both sides). To demonstrate, we first created the unbalanced reactions:
 ```@example chem1
@@ -96,10 +110,10 @@ using Catalyst # hide
 @variables t
 @species N(t) H(t) O(t) 
 @compounds begin
-    NH3(t) = N + 3H
-    O2(t) = 2O
-    NO(t) = N + O
-    H2O(t) = 2H + O
+    NH3 ~ N + 3H
+    O2 ~ 2O
+    NO ~ N + O
+    H2O ~ 2H + O
 end
 unbalanced_reaction = @reaction k, $NH3 + $O2 --> $NO + $H2O
 ```

src/chemistry_functionality.jl

@@ -28,17 +28,28 @@ macro compound(expr)
     make_compound(MacroTools.striplines(expr))
 end
 
+# Declares compound error messages:
+const COMPOUND_CREATION_ERROR_BASE = "Malformed input to @compound. Should use form like e.g. \"@compound CO2 ~ C + 2O\"."
+const COMPOUND_CREATION_ERROR_BAD_SEPARATOR = "Malformed input to @compound. Left-hand side (the compound) and the right-hand side (the components) should be separated by a \"~\" (e.g. \"@compound CO2 ~ C + 2O\"). If the left hand side contains metadata or default values, this should be enclosed by \"()\" (e.g. \"@compound (CO2 = 1.0, [output=true]) ~ C + 2O\")."
+const COMPOUND_CREATION_ERROR_DEPENDENT_VAR_GIVEN = "Left hand side of @compound is malformed. Does the compound depend on a independent variable (e.g. \"CO2(t)\")? If so, that should not be the case, these are inferred from the compounds."
+
 # Function managing the @compound macro.
 function make_compound(expr)
     # Error checks.
-    !(expr isa Expr) || (expr.head != :(=)) && error("Malformed expression. Compounds should be declared using a \"=\".")
-    (length(expr.args) != 2) && error("Malformed expression. Compounds should be consists of two expression, separated by a \"=\".")
-    ((expr.args[1] isa Symbol) || (expr.args[1].head != :call)) && error("Malformed expression. There must be a single compound which depend on an independent variable, e.g. \"CO2(t)\".")
+    (expr isa Expr) || error(COMPOUND_CREATION_ERROR_BASE)
+    ((expr.head == :call) && (expr.args[1] == :~) && (length(expr.args) == 3)) || error(COMPOUND_CREATION_ERROR_BAD_SEPARATOR)
+    if !(expr.args[2] isa Symbol) 
+        # Checks if a dependent variable was given (e.g. CO2(t)).
+        (expr.args[2].head == :call) && error(COMPOUND_CREATION_ERROR_DEPENDENT_VAR_GIVEN)
+        ((expr.args[2].args[1] isa Expr) && expr.args[2].args[1].head == :call) && error(COMPOUND_CREATION_ERROR_DEPENDENT_VAR_GIVEN)
+        ((expr.args[2].args[1] isa Expr) && (expr.args[2].args[1].args[1] isa Expr) && expr.args[2].args[1].args[1].head == :call) && error(COMPOUND_CREATION_ERROR_DEPENDENT_VAR_GIVEN) 
+    end
 
     # Extracts the composite species name, and a Vector{ReactantStruct} of its components.
-    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))
+    species_expr = expr.args[2]                                         # E.g. :(CO2(t))
+    species_expr = insert_independent_variable(species_expr, :t)        # Add independent variable (e.g. goes from `CO2` to `CO2(t)`).
+    species_name = find_varname_in_declaration(expr.args[2])            # E.g. :CO2
+    composition = Catalyst.recursive_find_reactants!(expr.args[3], 1, Vector{ReactantStruct}(undef, 0))
 
     # 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]).                                         
@@ -102,9 +113,13 @@ function make_compounds(expr)
     # 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, :([])))))
+    compound_syms = :([])
     for arg in expr.args
-        push!(compound_declarations.args[end].args[1].args, arg.args[1].args[1])
+        push!(compound_syms.args, find_varname_in_declaration(arg.args[2]))
     end
+    push!(compound_declarations.args, :($(Expr(:escape, :($(compound_syms))))))
+
+    # Returns output that.
     return compound_declarations
 end
 

src/expression_utils.jl

@@ -10,6 +10,73 @@ function get_tup_arg(ex::ExprValues, i::Int)
     return ex.args[i]
 end
 
+# In variable/species/parameters on the forms like:
+# X
+# X = 1.0
+# X, [metadata=true]
+# X = 1.0, [metadata=true]
+# X(t)
+# X(t) = 1.0
+# X(t), [metadata=true]
+# X(t) = 1.0, [metadata=true]
+# Finds the variable name (here, X).
+# Currently does not support e.g. "X, [metadata=true]" (when metadata does not have a comma before).
+function find_varname_in_declaration(expr)
+    (expr isa Symbol) && (return expr)           # Case: X
+    (expr.head == :call) && (return ex5.args[1]) # Case: X(t)
+    if expr.head == :(=)
+        (expr.args[1] isa Symbol) && (return expr.args[1])            # Case: X = 1.0
+        (expr.args[1].head == :call) && (return expr.args[1].args[1]) # Case: X(t) = 1.0
+    end
+    if expr.head == :tuple
+        (expr.args[1] isa Symbol) && (return expr.args[1])            # Case: X, [metadata=true]
+        (expr.args[1].head == :call) && (return expr.args[1].args[1]) # Case: X(t), [metadata=true]
+        if (expr.args[1].head == :(=)) 
+            (expr.args[1].args[1] isa Symbol) && (return expr.args[1].args[1])            # Case: X = 1.0, [metadata=true]
+            (expr.args[1].args[1].head == :call) && (return expr.args[1].args[1].args[1]) # Case: X(t) = 1.0, [metadata=true]
+        end
+    end
+    error("Unable to detect the variable declared in expression: $expr.")
+end
+
+# Converts an expression of the form:
+# X
+# X = 1.0
+# X, [metadata=true]
+# X = 1.0, [metadata=true]
+# To the form:
+# X(t)
+# X(t) = 1.0
+# X(t), [metadata=true]
+# X(t) = 1.0, [metadata=true]
+# (In this example the independent variable t was inserted).
+function insert_independent_variable(expr_in, iv)
+    # If expr is a symbol, just attach the iv. If not we have to create a new expr and mutate it. 
+    # Because Symbols (a possible input) cannot be mutated, this function cannot mutate the input (would have been easier if Expr input was guaranteed).
+    (expr_in isa Symbol) && (return :($(expr_in)($iv)))
+    expr = deepcopy(expr_in)
+
+    if expr.head == :(=) # Case: :(X = 1.0)
+        expr.args[1] = :($(expr.args[1])($iv))
+    elseif expr.head == :tuple
+        if expr.args[1] isa Symbol # Case: :(X, [metadata=true])
+            expr.args[1] = :($(expr.args[1])($iv))
+        elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol) # Case: :(X = 1.0, [metadata=true])
+            expr.args[1].args[1] = :($(expr.args[1].args[1])($iv))
+        end
+    end
+    (expr == expr_in) && error("Failed to add independent variable $(iv) to expression: $expr_in")
+    return expr
+end
+
+
+
+
+
+
+
+### Old Stuff ###
+
 #This will be called whenever a function stored in funcdict is called.
 # function replace_names(expr, old_names, new_names)
 #     mapping = Dict(zip(old_names, new_names))

src/reaction_network.jl

@@ -471,7 +471,8 @@ function read_compound_options(opts)
     # If the compound option is used retrive a list of compound species (need to be added to teh reaction system's species), and the option that creates them (used to declare them as compounds at the end).
     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.
+        # Find compound species names, and append the independent variable.
+        compound_species = [find_varname_in_declaration(arg.args[2]) for arg in compound_expr.args[3].args]
     else  # If option is not used, return empty vectors and expressions.
         compound_expr = :()
         compound_species = Union{Symbol, Expr}[]