init

Author: TorkelE

docs/src/api/catalyst_api.md

@@ -178,6 +178,9 @@ substoichmat
 prodstoichmat
 netstoichmat
 reactionrates
+get_metadata_dict
+has_metadata
+get_metadata
 ```
 
 ## [Functions to extend or modify a network](@id api_network_extension_and_modification)

src/Catalyst.jl

@@ -67,7 +67,7 @@ export ODEProblem,
        SteadyStateProblem
 
 # reaction_network macro
-const ExprValues = Union{Expr, Symbol, Float64, Int}
+const ExprValues = Union{Expr, Symbol, Float64, Int, Bool}
 include("expression_utils.jl")
 include("reaction_network.jl")
 export @reaction_network, @add_reactions, @reaction, @species
@@ -80,6 +80,7 @@ export mm, mmr, hill, hillr, hillar
 include("networkapi.jl")
 export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap
 export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap
+export get_metadata_dict, has_metadata, get_metadata
 export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants

src/expression_utils.jl

@@ -4,7 +4,8 @@ function tup_leng(ex::ExprValues)
     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).
+# 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]

src/reaction_network.jl

@@ -293,16 +293,17 @@ struct ReactionStruct
     substrates::Vector{ReactantStruct}
     products::Vector{ReactantStruct}
     rate::ExprValues
-    only_use_rate::Bool
+    metadata::Expr
 
-    function ReactionStruct(sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues,
-                            only_use_rate::Bool)
+    function ReactionStruct(sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, 
+                            metadata_line::ExprValues)
         sub = recursive_find_reactants!(sub_line, 1, Vector{ReactantStruct}(undef, 0))
         prod = recursive_find_reactants!(prod_line, 1, Vector{ReactantStruct}(undef, 0))
-        new(sub, prod, rate, only_use_rate)
+        metadata = extract_metadata(metadata_line)
+        new(sub, prod, rate, metadata)
     end
 end
-
+ 
 ### Functions rephrasing the macro input as a ReactionSystem structure. ###
 
 function forbidden_variable_check(v)
@@ -565,8 +566,8 @@ function get_rxexprs(rxstruct)
     prod_init = isempty(rxstruct.products) ? nothing : :([])
     prod_stoich_init = deepcopy(prod_init)
     reaction_func = :(Reaction($(recursive_expand_functions!(rxstruct.rate)), $subs_init,
-                               $prod_init, $subs_stoich_init, $prod_stoich_init,
-                               only_use_rate = $(rxstruct.only_use_rate)))
+                               $prod_init, $subs_stoich_init, $prod_stoich_init, 
+                               metadata = Dict($(rxstruct.metadata)),))
     for sub in rxstruct.substrates
         push!(reaction_func.args[3].args, sub.reactant)
         push!(reaction_func.args[5].args, sub.stoichiometry)
@@ -580,70 +581,82 @@ end
 
 ### Functions for extracting the reactions from a DSL expression, and putting them ReactionStruct vector. ###
 
-# Reads a line and creates the corresponding ReactionStruct.
+# Reads a single line and creates the corresponding ReactionStruct.
 function get_reaction(line)
-    (rate, r_line) = line.args
-    (r_line.head == :-->) && (r_line = Expr(:call, :→, r_line.args[1], r_line.args[2]))
-
-    arrow = r_line.args[1]
-    in(arrow, double_arrows) && error("Double arrows not allowed for single reactions.")
-
-    only_use_rate = in(arrow, pure_rate_arrows)
-    if in(arrow, fwd_arrows)
-        rs = create_ReactionStruct(r_line.args[2], r_line.args[3], rate, only_use_rate)
-    elseif in(arrow, bwd_arrows)
-        rs = create_ReactionStruct(r_line.args[3], r_line.args[2], rate, only_use_rate)
-    else
-        throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))")
+    reaction = get_reactions([line])
+    if (length(reaction) != 1)
+        error("Malformed reaction. @reaction macro only creates a single reaction. E.g. double arrows, such as `<-->` are not supported.")
     end
-
-    rs
+    return reaction[1]
 end
+
 # Generates a vector containing a number of reaction structures, each containing the information about one reaction.
 function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0))
     for line in exprs
-        (rate, r_line) = line.args
-        (r_line.head == :-->) && (r_line = Expr(:call, :→, r_line.args[1], r_line.args[2]))
+        # Reads core reaction information.
+        arrow, rate, reaction, metadata = read_reaction_line(line)
 
-        arrow = r_line.args[1]
-        only_use_rate = in(arrow, pure_rate_arrows)
+        # Checks the type of arrow used, and creates the corresponding reaction(s). Returns them in an array.
         if in(arrow, double_arrows)
-            (typeof(rate) == Expr && rate.head == :tuple) ||
+            if typeof(rate) != Expr || rate.head != :tuple
                 error("Error: Must provide a tuple of reaction rates when declaring a bi-directional reaction.")
-            push_reactions!(reactions, r_line.args[2], r_line.args[3], rate.args[1],
-                            only_use_rate)
-            push_reactions!(reactions, r_line.args[3], r_line.args[2], rate.args[2],
-                            only_use_rate)
+            end
+            push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow)
+            push_reactions!(reactions, reaction.args[3], reaction.args[2], rate.args[2], metadata.args[2], arrow)
         elseif in(arrow, fwd_arrows)
-            push_reactions!(reactions, r_line.args[2], r_line.args[3], rate, only_use_rate)
+            push_reactions!(reactions, reaction.args[2], reaction.args[3], rate, metadata, arrow)
         elseif in(arrow, bwd_arrows)
-            push_reactions!(reactions, r_line.args[3], r_line.args[2], rate, only_use_rate)
+            push_reactions!(reactions, reaction.args[3], reaction.args[2], rate, metadata, arrow)
         else
             throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))")
         end
     end
-    reactions
-end
+    return reactions
+end
+
+# Extracts the rate, reaction, and metadata fields (the last one optional) from a reaction line.
+function read_reaction_line(line::Expr)
+    # Handles rate, reaction, and arrow.
+    # Special routine required for  the`-->` case, which creates different expression from all other cases.
+    rate = line.args[1]
+    reaction = line.args[2]
+    (reaction.head == :-->) && (reaction = Expr(:call, :→, reaction.args[1], reaction.args[2]))
+    arrow = reaction.args[1]
+
+    # Handles metadata. If not provided, empty metadata is created.
+    if length(line.args) == 2
+        metadata = in(arrow, double_arrows) ? :(([], [])) : :([])
+    elseif length(line.args) == 3
+        metadata = line.args[3]
+    else
+        error("The following reaction line: \"$line\" was malformed. It should have a form \"rate, reaction\" or a form \"rate, reaction, metadata\". It has neither.")
+    end
 
-# Creates a ReactionStruct from the information in a single line.
-function create_ReactionStruct(sub_line::ExprValues, prod_line::ExprValues,
-                               rate::ExprValues, only_use_rate::Bool)
-    all(==(1), (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate))) ||
-        error("Malformed reaction, line appears to be defining multiple reactions incorrectly: rate=$rate, subs=$sub_line, prods=$prod_line.")
-    ReactionStruct(get_tup_arg(sub_line, 1), get_tup_arg(prod_line, 1),
-                   get_tup_arg(rate, 1), only_use_rate)
+    return arrow, rate, reaction, metadata
 end
 
-#Takes a reaction line and creates reactions from it and pushes those to the reaction array. Used to create multiple reactions from, for instance, 1.0, (X,Y) --> 0.
-function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues,
-                         prod_line::ExprValues, rate::ExprValues, only_use_rate::Bool)
-    lengs = (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate))
-    for i in 1:maximum(lengs)
-        (count(lengs .== 1) + count(lengs .== maximum(lengs)) < 3) &&
-            (throw("Malformed reaction, rate=$rate, subs=$sub_line, prods=$prod_line."))
-        push!(reactions,
-              ReactionStruct(get_tup_arg(sub_line, i), get_tup_arg(prod_line, i),
-                             get_tup_arg(rate, i), only_use_rate))
+# Takes a reaction line and creates reaction(s) from it and pushes those to the reaction array.
+# Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`.
+# Handles metadata, e.g. `1.0, Z --> 0, [noisescaling=η]`.
+function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, 
+                         rate::ExprValues, metadata::ExprValues, arrow::Symbol)  
+    # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`).
+    # This finds the length of these tuples (or 1 if not in tuple forms). Errors if lengs inconsistent.
+    lengs = (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate), tup_leng(metadata))
+    if any(!(leng in [1, maximum(lengs)]) for leng in lengs)
+        throw("Malformed reaction, rate=$rate, subs=$sub_line, prods=$prod_line, metadata=$metadata.")
+    end
+
+    # Loops through each reaction encoded by the reaction composites. Adds the reaction to the reactions vector.
+    for i in 1:maximum(lengs)                       
+        # If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used.
+        metadata_i = get_tup_arg(metadata, i)
+        if all([arg.args[1] != :only_use_rate for arg in metadata_i.args])
+            push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows))))
+        end
+
+        push!(reactions, ReactionStruct(get_tup_arg(sub_line, i), get_tup_arg(prod_line, i),
+                                        get_tup_arg(rate, i), metadata_i))
     end
 end
 
@@ -687,6 +700,17 @@ function recursive_find_reactants!(ex::ExprValues, mult::ExprValues,
     reactants
 end
 
+# Finds the metadata from a metadata expresion (`[key=val, ...]`) and returns as a Vector{Pair{Symbol,ExprValues}}.
+function extract_metadata(metadata_line::Expr)
+    metadata = :([])
+    for arg in metadata_line.args
+        (arg.head != :(=)) && error("Malformatted metadata line: $metadata_line. Each entry in the vector should contain a `=`.")
+        (arg.args[1] isa Symbol) || error("Malformatted metadata entry: $arg. Entries left-hand-side should be a single symbol.")
+        push!(metadata.args, :($(QuoteNode(arg.args[1])) => $(arg.args[2])))
+    end
+    return metadata
+end
+
 ### DSL Options Handling ###
 # Most options handled in previous sections, when code re-organised, these should ideally be moved to the same place.
 
@@ -888,3 +912,13 @@ end
 #     end
 #     return defaults
 # end
+
+
+# # Creates a ReactionStruct from the information in a single line.
+# function create_ReactionStruct(sub_line::ExprValues, prod_line::ExprValues,
+#     rate::ExprValues, only_use_rate::Bool)
+#     all(==(1), (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate))) ||
+#     error("Malformed reaction, line appears to be defining multiple reactions incorrectly: rate=$rate, subs=$sub_line, prods=$prod_line.")
+#     ReactionStruct(get_tup_arg(sub_line, 1), get_tup_arg(prod_line, 1),
+#     get_tup_arg(rate, 1), only_use_rate)
+# end

src/reactionsystem.jl

@@ -98,7 +98,7 @@ Notes:
 - The three-argument form assumes all reactant and product stoichiometric coefficients
   are one.
 """
-struct Reaction{S, T}
+struct Reaction{R, S, T}
     """The rate function (excluding mass action terms)."""
     rate::Any
     """Reaction substrates."""
@@ -116,6 +116,11 @@ struct Reaction{S, T}
     `true` if `rate` represents the full reaction rate law.
     """
     only_use_rate::Bool
+    """
+    Contain additional data, such whenever the reaction have a specific noise-scaling expression for
+    the chemical Langevin equation.
+    """
+    metadata::Dict{Symbol, R}
 end
 
 """
@@ -126,7 +131,8 @@ Test if a species is valid as a reactant (i.e. a species variable or a constant
 isvalidreactant(s) = MT.isparameter(s) ? isconstant(s) : (isspecies(s) && !isconstant(s))
 
 function Reaction(rate, subs, prods, substoich, prodstoich;
-                  netstoich = nothing, only_use_rate = false,
+                  netstoich = nothing, metadata = Dict{Symbol,Any}(), 
+                  only_use_rate = (haskey(metadata, :only_use_rate) && Bool(metadata[:only_use_rate])),
                   kwargs...)
     (isnothing(prods) && isnothing(subs)) &&
         throw(ArgumentError("A reaction requires a non-nothing substrate or product vector."))
@@ -181,7 +187,9 @@ function Reaction(rate, subs, prods, substoich, prodstoich;
         convert.(stoich_type, netstoich) : netstoich
     end
 
-    Reaction(value(rate), subs, prods, substoich′, prodstoich′, ns, only_use_rate)
+    haskey(metadata, :only_use_rate) && delete!(metadata, :only_use_rate)
+
+    Reaction(value(rate), subs, prods, substoich′, prodstoich′, ns, only_use_rate, metadata)
 end
 
 # three argument constructor assumes stoichiometric coefs are one and integers
@@ -241,7 +249,7 @@ function ModelingToolkit.namespace_equation(rx::Reaction, name; kw...)
         ns = similar(rx.netstoich)
         map!(n -> f(n[1]) => f(n[2]), ns, rx.netstoich)
     end
-    Reaction(rate, subs, prods, substoich, prodstoich, netstoich, rx.only_use_rate)
+    Reaction(rate, subs, prods, substoich, prodstoich, netstoich, rx.only_use_rate, rx.metadata)
 end
 
 netstoich_stoichtype(::Vector{Pair{S, T}}) where {S, T} = T
@@ -830,6 +838,67 @@ function get_indep_sts(rs::ReactionSystem, remove_conserved = false)
     indepsts, filter(isspecies, indepsts)
 end
 
+######################## Other accessors ##############################
+
+"""
+    get_metadata_dict(reaction::Reaction)
+
+Retrives the `ImmutableDict` containing all of the metadata associated with a specific reaction.
+
+Arguments:
+- `reaction`: The reaction for which we wish to retrive all metadata.
+
+Example:
+```julia
+reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
+get_metadata_dict(reaction)
+```
+"""
+function get_metadata_dict(reaction::Reaction)
+    return reaction.metadata
+end
+
+"""
+    has_metadata(reaction::Reaction, md_key::Symbol)
+
+Checks if a `Reaction` have a certain metadata field. If it does, returns `true` (else returns `false`).
+
+Arguments:
+- `reaction`: The reaction for which we wish to check if a specific metadata field exist.
+- `md_key`: The metadata for which we wish to check existence of.
+
+Example:
+```julia
+reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
+has_metadata(reaction, :noise_scaling)
+```
+"""
+function has_metadata(reaction::Reaction, md_key::Symbol)
+    return haskey(get_metadata_dict(reaction), md_key)
+end
+
+"""
+    get_metadata(reaction::Reaction, md_key::Symbol)
+
+Retrives a certain metadata value from a `Reaction`. If the metadata does not exists, throws an error.
+
+Arguments:
+- `reaction`: The reaction for which we wish to retrive a specific metadata value.
+- `md_key`: The metadata for which we wish to retrive.
+
+Example:
+```julia
+reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
+get_metadata(reaction, :noise_scaling)
+```
+"""
+function get_metadata(reaction::Reaction, md_key::Symbol)
+    if !has_metadata(reaction, md_key) 
+        error("The reaction does not have a metadata field $md_key. It does have the following metadata fields: $(keys(get_metadata_dict(reaction))).")
+    end
+    return get_metadata_dict(reaction)[md_key]
+end
+
 ######################## Conversion to ODEs/SDEs/jump, etc ##############################
 
 """

src/registered_functions.jl

@@ -134,7 +134,8 @@ function expand_registered_functions(expr)
 end
 # If applied to a Reaction, return a reaction with its rate modified.
 function expand_registered_functions(rx::Reaction)
-    Reaction(expand_registered_functions(rx.rate), rx.substrates, rx.products, rx.substoich, rx.prodstoich, rx.netstoich, rx.only_use_rate)
+    Reaction(expand_registered_functions(rx.rate), rx.substrates, rx.products, rx.substoich, 
+             rx.prodstoich, rx.netstoich, rx.only_use_rate, rx.metadata)
 end
 # If applied to a Equation, returns it with it applied to lhs and rhs
 function expand_registered_functions(eq::Equation)