From f7039d1993fa64f063033b5c9b425bf85c65a5a5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 7 Nov 2023 16:32:20 -0500 Subject: [PATCH 01/30] init --- .../homotopy_continuation_extension.jl | 25 ++++++++++++++++++- src/reactionsystem.jl | 10 ++++++++ 2 files changed, 34 insertions(+), 1 deletion(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index f9c1a5c19b..844cedfc19 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -104,4 +104,27 @@ function filter_negative_f(sols; neg_thres=-1e-20) (neg_thres < sol[idx] < 0) && (sol[idx] = 0) end return filter(sol -> all(>=(0), sol), sols) -end \ No newline at end of file +end + +### Archived ### + +# # Unfolds a function (like mm or hill). +# function deregister(fs::Vector{T}, expr) where T +# for f in fs +# expr = deregister(f, expr) +# end +# return expr +# end +# # Provided by Shashi Gowda. +# deregister(f, expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___deregister(f)))(unwrap(expr))) +# function ___deregister(f) +# (expr) -> +# if istree(expr) && operation(expr) == f +# args = arguments(expr) +# invoke_with = map(args) do a +# t = typeof(a) +# issym(a) || istree(a) ? wrap(a) => symtype(a) : a => typeof(a) +# end +# invoke(f, Tuple{last.(invoke_with)...}, first.(invoke_with)...) +# end +# end \ No newline at end of file diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 4fdb99f3ee..dad8562c3a 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1296,6 +1296,11 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs) eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved, include_zero_odes) eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved) + + # Converts expressions like mm(X,v,K) to v*X/(X+K). + expand_functions && for eq in eqs + eq.rhs = expand_registered_functions!(eq.rhs) + end ODESystem(eqs, get_iv(fullrs), sts, ps; observed = obs, @@ -1404,6 +1409,11 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; eqs, sts, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved) ps = (noise_scaling === nothing) ? ps : vcat(ps, toparam(noise_scaling)) + # Converts expressions like mm(X,v,K) to v*X/(X+K). + expand_functions && for eq in eqs + eq.rhs = expand_registered_functions!(eq.rhs) + end + if any(isbc, get_states(flatrs)) @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead." end From e77abcf91c3adef67b485c42374bb3a8abfd76c7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 7 Nov 2023 17:04:57 -0500 Subject: [PATCH 02/30] add tests --- src/reactionsystem.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index dad8562c3a..19525ff953 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1298,9 +1298,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs) eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved) # Converts expressions like mm(X,v,K) to v*X/(X+K). - expand_functions && for eq in eqs - eq.rhs = expand_registered_functions!(eq.rhs) - end + expand_functions && (eqs = [eq.lhs ~ expand_registered_functions!(eq.rhs) for eq in eqs]) ODESystem(eqs, get_iv(fullrs), sts, ps; observed = obs, @@ -1410,8 +1408,9 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; ps = (noise_scaling === nothing) ? ps : vcat(ps, toparam(noise_scaling)) # Converts expressions like mm(X,v,K) to v*X/(X+K). - expand_functions && for eq in eqs - eq.rhs = expand_registered_functions!(eq.rhs) + if expand_functions + eqs = [eq.lhs ~ expand_registered_functions!(eq.rhs) for eq in eqs] + noiseeqs = [expand_registered_functions!(neq) for neq in noiseeqs] end if any(isbc, get_states(flatrs)) From a921ed4d22d51f8e4831cd6f6255560abdc203d3 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 13 Nov 2023 22:22:12 -0500 Subject: [PATCH 03/30] partial progress --- src/reactionsystem.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 19525ff953..a1168c98a7 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -538,6 +538,11 @@ struct ReactionSystem{V <: NetworkProperties} <: checks && validate(rs) rs end + + # Copies a reaction system, but with the option of having some fields replaced + function ReactionSystem(rs::ReactionSystem; eqs = rs.eqs, rxs = rs.rxs, iv = rs.iv, sivs = rs.sivs, states = rs.states, species = rs.species, ps = rs.ps, var_to_name = rs.var_to_name, observed = rs.observed, name = rs.name, systems = rs.systems, defaults = rs.defaults, connection_type = rs.connection_type, networkproperties = rs.networkproperties, combinatoric_ratelaws = rs.combinatoric_ratelaws, continuous_events = rs.continuous_events, discrete_events = rs.discrete_events, complete = rs.complete) + new{typeof(networkproperties)}(eqs, rxs, ModelingToolkit.unwrap(iv), ModelingToolkit.unwrap.(sivs), ModelingToolkit.unwrap.(states), ModelingToolkit.unwrap.(species), ModelingToolkit.unwrap.(ps), var_to_name, observed, name, systems, defaults, connection_type, networkproperties, combinatoric_ratelaws, continuous_events, discrete_events, complete) + end end function get_speciestype(iv, states, systems) From bad52efbcb165f581d6d048862e01d6388b57ad4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 14 Nov 2023 08:43:33 -0500 Subject: [PATCH 04/30] up --- .../homotopy_continuation_extension.jl | 25 +------------------ src/reactionsystem.jl | 14 ----------- 2 files changed, 1 insertion(+), 38 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 844cedfc19..f9c1a5c19b 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -104,27 +104,4 @@ function filter_negative_f(sols; neg_thres=-1e-20) (neg_thres < sol[idx] < 0) && (sol[idx] = 0) end return filter(sol -> all(>=(0), sol), sols) -end - -### Archived ### - -# # Unfolds a function (like mm or hill). -# function deregister(fs::Vector{T}, expr) where T -# for f in fs -# expr = deregister(f, expr) -# end -# return expr -# end -# # Provided by Shashi Gowda. -# deregister(f, expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___deregister(f)))(unwrap(expr))) -# function ___deregister(f) -# (expr) -> -# if istree(expr) && operation(expr) == f -# args = arguments(expr) -# invoke_with = map(args) do a -# t = typeof(a) -# issym(a) || istree(a) ? wrap(a) => symtype(a) : a => typeof(a) -# end -# invoke(f, Tuple{last.(invoke_with)...}, first.(invoke_with)...) -# end -# end \ No newline at end of file +end \ No newline at end of file diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index a1168c98a7..4fdb99f3ee 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -538,11 +538,6 @@ struct ReactionSystem{V <: NetworkProperties} <: checks && validate(rs) rs end - - # Copies a reaction system, but with the option of having some fields replaced - function ReactionSystem(rs::ReactionSystem; eqs = rs.eqs, rxs = rs.rxs, iv = rs.iv, sivs = rs.sivs, states = rs.states, species = rs.species, ps = rs.ps, var_to_name = rs.var_to_name, observed = rs.observed, name = rs.name, systems = rs.systems, defaults = rs.defaults, connection_type = rs.connection_type, networkproperties = rs.networkproperties, combinatoric_ratelaws = rs.combinatoric_ratelaws, continuous_events = rs.continuous_events, discrete_events = rs.discrete_events, complete = rs.complete) - new{typeof(networkproperties)}(eqs, rxs, ModelingToolkit.unwrap(iv), ModelingToolkit.unwrap.(sivs), ModelingToolkit.unwrap.(states), ModelingToolkit.unwrap.(species), ModelingToolkit.unwrap.(ps), var_to_name, observed, name, systems, defaults, connection_type, networkproperties, combinatoric_ratelaws, continuous_events, discrete_events, complete) - end end function get_speciestype(iv, states, systems) @@ -1301,9 +1296,6 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs) eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved, include_zero_odes) eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved) - - # Converts expressions like mm(X,v,K) to v*X/(X+K). - expand_functions && (eqs = [eq.lhs ~ expand_registered_functions!(eq.rhs) for eq in eqs]) ODESystem(eqs, get_iv(fullrs), sts, ps; observed = obs, @@ -1412,12 +1404,6 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; eqs, sts, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved) ps = (noise_scaling === nothing) ? ps : vcat(ps, toparam(noise_scaling)) - # Converts expressions like mm(X,v,K) to v*X/(X+K). - if expand_functions - eqs = [eq.lhs ~ expand_registered_functions!(eq.rhs) for eq in eqs] - noiseeqs = [expand_registered_functions!(neq) for neq in noiseeqs] - end - if any(isbc, get_states(flatrs)) @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead." end From f9cab9e76b7ad28b4c59a4bae8d8f66f34b23c3d Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 10 Nov 2023 17:42:35 -0500 Subject: [PATCH 05/30] Init --- docs/pages.jl | 1 + .../chemistry_related_functionality.md | 85 +++++++ .../constraint_equations.md | 2 +- src/Catalyst.jl | 4 +- src/chemistry_functionality.jl | 211 +++++++++++++----- test/miscellaneous_tests/compound_macro.jl | 80 +++++-- 6 files changed, 304 insertions(+), 79 deletions(-) create mode 100644 docs/src/catalyst_functionality/chemistry_related_functionality.md diff --git a/docs/pages.jl b/docs/pages.jl index 9b21e12b10..51f3c02863 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -7,6 +7,7 @@ pages = Any["Home" => "index.md", "catalyst_functionality/constraint_equations.md", "catalyst_functionality/parametric_stoichiometry.md", "catalyst_functionality/network_analysis.md", + "catalyst_functionality/chemistry_related_functionality.md", "Model creation examples" => Any["catalyst_functionality/example_networks/basic_CRN_examples.md", "catalyst_functionality/example_networks/hodgkin_huxley_equation.md", "catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md"]], diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md new file mode 100644 index 0000000000..d3e9eee86f --- /dev/null +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -0,0 +1,85 @@ +# [Chemistry-related functionality](@id chemistry_functionality) + +While Catalyst has primarily been designed around the modelling of biological systems, reaction network models are also common across chemistry. This section describes two types of functionality, that while of general interest, should be especially useful in the modelling of chemical systems. +- The `@compound` option, allowing the user to designate that a specific species is composed of certain subspecies. +- The `balance_reaction` function enabling the user to balance a reactions so the same number of compounds occur on both sides. + +## Modelling with compound species +Defining compound species is currently only supported for [programmatic construction](@ref programmatic_CRN_construction) of reactions and reaction network models. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ compound species, consisting of 1 C species and 2 O species. First we create species corresponding to the components: +```@example chem1 +@variables t +@species C(t) O(t) +``` +Next, we create the `CO2` compound: +```@example chem1 +@compound CO2(t) = 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 variable. Components with non-unitary stoichiometry have this value written before the component (generally, the rules for designating the components of a compounds are identical to hose 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` could be used: +```@example chem1 +isspecies(CO2) +``` +However, in its metadata is stored the information of its components, which can be retrieved using the `components` (returning a vector of its component species) and `coefficients` (returning a vector with each component's stoichiometry) functions: +```@example chem1 +components(CO2) +``` +```@example chem1 +coefficients(CO2) +``` +Alternatively, we can retrieve the components and their stoichiometric coefficients as a single vector using the `component_coefficients` function: +```@example chem1 +component_coefficients(CO2) +``` +Finally, it is possible to check whether a species is a compound or not using the `iscompound` function: +```@example chem1 +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 +``` + +When multiple compounds are created, they can be created simultaneously using the `@compounds` macro, e.g. the previous code-block could have been executed using: +```@example chem1 +@species H(t) +@compounds begin + H2O(t) = 2H + O + H2CO3(t) = CO2 + H2O +end +``` + +One use defining a species as a compound is that they can be used to balance reactions to that the number of compounds are the same on both side. + +## Balancing chemical reactions +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` enable us to find coefficients creating a balanced reaction (in this case, where the number of carbon and oxygen atoms are teh same on both sides). To demonstrate, we first created the unbalanced reactions: +```@example chem1 +rx = @reaction k, C + O --> $CO2 +``` +Here, we set a reaction rate `k` (which is not involved in the reaction balancing). We also use interpolation of `CO2`, ensuring that the `CO2` used in the reaction is the same one we previously defined as a compound of `C` and `O`. Next, we call the `balance_reaction` function +```@example chem1 +balance_reaction(rx) +``` +which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). + +Let us consider a more elaborate example, the reaction between ammonia (NH₃) and oxygen (O₂) to form nitrogen monoxide (NO) and water (H₂O). Let us first create the components and the unbalanced reaction: +```@example chem2 +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 +end +unbalanced_reaction = @reaction k, $NH3 + $O2 --> $NO + $H2O +``` +We can now created a balanced version (where the amount of H, N, and O is the same on both sides): +```@example chem2 +balanced_reaction = balance_reaction(unbalanced_reaction)[1] +``` + +!!! note + Reaction balancing is currently not supported for reactions involving compounds of compounds. \ No newline at end of file diff --git a/docs/src/catalyst_functionality/constraint_equations.md b/docs/src/catalyst_functionality/constraint_equations.md index 29dcbb2fc1..c0fd74bf99 100644 --- a/docs/src/catalyst_functionality/constraint_equations.md +++ b/docs/src/catalyst_functionality/constraint_equations.md @@ -17,7 +17,7 @@ $\lambda$. For now we'll assume the cell can grow indefinitely. We'll also keep track of one protein $P(t)$, which is produced at a rate proportional $V$ and can be degraded. -## Coupling ODE constraints via extending a system +## [Coupling ODE constraints via extending a system](@id constraint_equations_coupling_constratins) There are several ways we can create our Catalyst model with the two reactions and ODE for $V(t)$. One approach is to use compositional modeling, create diff --git a/src/Catalyst.jl b/src/Catalyst.jl index fa2b3e51c2..42398a7c95 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -102,8 +102,8 @@ export Graph, savegraph, complexgraph # for creating compounds include("chemistry_functionality.jl") -export @compound -export components, iscompound, coefficients +export @compound, @compounds +export iscompound, components, coefficients, component_coefficients export balance_reaction ### Extensions ### diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 96b0118f39..ee1f809f48 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -1,3 +1,4 @@ +### Declares Compound Related Metadata ### struct CompoundSpecies end struct CompoundComponents end struct CompoundCoefficients end @@ -6,93 +7,138 @@ Symbolics.option_to_metadata_type(::Val{:iscompound}) = CompoundSpecies Symbolics.option_to_metadata_type(::Val{:components}) = CompoundComponents Symbolics.option_to_metadata_type(::Val{:coefficients}) = CompoundCoefficients -macro compound(species_expr, arr_expr...) - # Ensure the species name is a valid expression - if !(species_expr isa Expr && species_expr.head == :call) - error("Invalid species name in @compound macro") - end +### Create @compound Macro(s) ### - # Parse the species name to extract the species name and argument - species_name = species_expr.args[1] - species_arg = species_expr.args[2] - - # Construct the expressions that define the species - species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0), - Expr(:call, species_name, species_arg)) - - # Construct the expression to set the iscompound metadata - setmetadata_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundSpecies, - true)) - - # Ensure the expressions are evaluated in the correct scope by escaping them - escaped_species_expr = esc(species_expr) - escaped_setmetadata_expr = esc(setmetadata_expr) - - # Construct the array from the remaining arguments - arr = Expr(:vect, (arr_expr)...) - coeffs = [] - species = [] - - for expr in arr_expr - if isa(expr, Expr) && expr.head == :call && expr.args[1] == :* - push!(coeffs, expr.args[2]) - push!(species, expr.args[3]) - else - push!(coeffs, 1) - push!(species, expr) - end - end +""" + @compound + +Macro that creates a compound species, which is composed of smaller constituent species. + +Example: +```julia +@variables t +@species C(t) O(t) +@compound CO2(t) = C + 2O +``` + +Notes: +- The constituent species must be defined before using the `@compound` macro. +""" +macro compound(expr) + make_compound(MacroTools.striplines(expr)) +end - coeffs_expr = Expr(:vect, coeffs...) - species_expr = Expr(:vect, species...) +# 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)\".") + + # 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)) + 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] + + # 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). + species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)` + compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)` + 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])` + + # Returns the rephrased expression. + return quote + $species_declaration_expr + $compound_designation_expr + $components_designation_expr + $coefficients_designation_expr + end +end - # Construct the expression to set the components metadata - setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundComponents, - $species_expr)) +""" + @compounds - # Ensure the expression is evaluated in the correct scope by escaping it - escaped_setcomponents_expr = esc(setcomponents_expr) +Macro that creates several compound species, which each is composed of smaller constituent species. Uses the same syntax as `@compound`, but with one compound species one each line. - # Construct the expression to set the coefficients metadata - setcoefficients_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundCoefficients, - $coeffs_expr)) +Example: +```julia +@variables t +@species C(t) H(t) O(t) +@compounds + CH4(t) = C + 4H + O2(t) = 2O + CO2(t) = C + 2O + H2O(t) = 2H + O +end +``` - escaped_setcoefficients_expr = esc(setcoefficients_expr) +Notes: +- The constituent species must be defined before using the `@compound` macro. +""" +macro compounds(expr) + make_compounds(MacroTools.striplines(expr)) +end - # Return a block that contains the escaped expressions - return Expr(:block, escaped_species_expr, escaped_setmetadata_expr, - escaped_setcomponents_expr, escaped_setcoefficients_expr) +# 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 :(.) end -# Check if a species is a compound +## Compound Getters ### + +""" + iscompound(s) + +Returns `true` if the input is a compound species (else false). +""" iscompound(s::Num) = iscompound(MT.value(s)) function iscompound(s) MT.getmetadata(s, CompoundSpecies, false) end -coefficients(s::Num) = coefficients(MT.value(s)) -function coefficients(s) - MT.getmetadata(s, CompoundCoefficients) -end +""" + components(s) +Returns a vector with a list of all the components of a compound species (created using e.g. the @compound macro). +""" components(s::Num) = components(MT.value(s)) function components(s) MT.getmetadata(s, CompoundComponents) end +""" + coefficients(s) + +Returns a vector with a list of all the stoichiometric coefficients of the components of a compound species (created using e.g. the @compound macro). +""" +coefficients(s::Num) = coefficients(MT.value(s)) +function coefficients(s) + MT.getmetadata(s, CompoundCoefficients) +end + +""" + component_coefficients(s) + +Returns a Vector{Pari{Symbol,Int64}}, listing a compounds species (created using e.g. the @compound macro) all the coefficients and their stoichiometric coefficients. +""" component_coefficients(s::Num) = component_coefficients(MT.value(s)) function component_coefficients(s) return [c => co for (c, co) in zip(components(s), coefficients(s))] end +### Reaction Balancing Functionality ### -### Balancing Code +# Reaction balancing error. const COMPOUND_OF_COMPOUND_ERROR = ErrorException("Reaction balancing does not currently work for reactions involving compounds of compounds.") -# note this does not correctly handle compounds of compounds currently +# Note this does not correctly handle compounds of compounds currently. +# Internal function used by "balance_reaction" (via "get_balanced_stoich"). function create_matrix(reaction::Catalyst.Reaction) @unpack substrates, products = reaction unique_atoms = [] # Array to store unique atoms @@ -143,6 +189,7 @@ function create_matrix(reaction::Catalyst.Reaction) return A end +# Internal function used by "balance_reaction". function get_balanced_stoich(reaction::Reaction) # Create the reaction matrix A that is m atoms by n compounds A = create_matrix(reaction) @@ -171,6 +218,52 @@ function get_balanced_stoich(reaction::Reaction) return stoichvecs end +""" + balance_reaction(reaction::Reaction) + +Returns a vector of all possible stoichiometrically balanced `Reaction` objects for the given `Reaction`. + +Example: +```julia +@variables t +@species Si(t) Cl(t) H(t) O(t) +@compound SiCl4(t) = Si + 4Cl +@compound H2O(t) = 2H + O +@compound H4SiO4(t) = 4H + Si + 4O +@compound HCl(t) = H + Cl +rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl]) +balance_reaction(rx) # Exactly one solution. +``` + +```julia +@variables t +@species C(t) H(t) O(t) +@compound CO(t) = C + O +@compound CO2(t) = C + 2O +@compound H2(t) = 2H +@compound CH4(t) = C + 4H +@compound H2O(t) = 2H + O +rx = Reaction(1.0, [CO, CO2, H2], [CH4, H2O]) +balance_reaction(rx) # Multiple solutions. +``` + +```julia +@variables t +@species Fe(t) S(t) O(t) H(t) N(t) +@compound FeS2(t) = Fe + 2S +@compound HNO3(t) = H + N + 3O +@compound Fe2S3O12(t) = 2Fe + 3S + 12O +@compound NO(t) = N + O +@compound H2SO4(t) = 2H + S + 4O +rx = Reaction(1.0, [FeS2, HNO3], [Fe2S3O12, NO, H2SO4]) +brxs = balance_reaction(rx) # No solution. +``` + +Notes: +- Balancing reactions that contain compounds of compounds is currently not supported. +- A reaction may not always yield a single solution; it could have an infinite number of solutions or none at all. When there are multiple solutions, a vector of all possible `Reaction` objects is returned. However, substrates and products may be interchanged as we currently do not solve for a linear combination that maintains the set of substrates and products. +- If the reaction cannot be balanced, an empty `Reaction` vector is returned. +""" function balance_reaction(reaction::Reaction) # Calculate the stoichiometric coefficients for the balanced reaction. stoichiometries = get_balanced_stoich(reaction) diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index be3afbb035..d13c72cd7b 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -3,8 +3,8 @@ using Catalyst, Test # Test base funcationality in two cases. let @variables t - @species C(t) H(t) O(t) - @compound C6H12O2(t) 6C 12H 2O + @species C(t) H(t) + O(t) + @compound C6H12O2(t) = 6C + 12H + 2O @test iscompound(C6H12O2) @test isspecies(C6H12O2) @@ -21,7 +21,7 @@ end let @variables t @species O(t) - @compound O2(t) 2O + @compound O2(t) = 2O @test iscompound(O2) @test isspecies(O2) @@ -37,19 +37,19 @@ end let @variables t @species C(t) H(t) - @test_throws Exception @compound C6H12O2(t) 6C 12H 2O + @test_throws Exception @compound C6H12O2(t) = 6C + 12H + 2O end let @variables t - @test_throws Exception @compound O2(t) 2O + @test_throws Exception @compound O2(t) = 2O end # Checks that nested components works as expected. let @variables t @species C(t) H(t) O(t) - @compound OH(t) 1O 1H - @compound C3H5OH3(t) 3C 5H 3OH + @compound OH(t) = 1O + 1H + @compound C3H5OH3(t) = 3C + 5H + 3OH @test !iscompound(O) @test !iscompound(H) @@ -71,22 +71,23 @@ let @variables t @species C(t) H(t) O(t) s = C - @compound C6H12O2_1(t) 6s 12H 2O - @compound C6H12O2_2(t) 6C 12H 2O + @compound C6H12O2_1(t) = 6s + 12H + 2O + @compound C6H12O2_2(t) = 6C + 12H + 2O @test iscompound(C6H12O2_1) @test iscompound(C6H12O2_2) @test isequal(components(C6H12O2_1), components(C6H12O2_2)) @test isequal(coefficients(C6H12O2_1), coefficients(C6H12O2_2)) + @test isequal(component_coefficients(C6H12O2_1), component_coefficients(C6H12O2_2)) end let @variables t @species C(t) H(t) - @compound Cyclopentadiene(t) 5C 6H + @compound Cyclopentadiene(t) = 5C + 6H C5H6 = Cyclopentadiene - @compound C10H12(t) 2C5H6 + @compound C10H12(t) = 2C5H6 @test iscompound(C10H12) @test iscompound(components(C10H12)[1]) @@ -101,8 +102,8 @@ let @species H(t) alpha = 2 - @compound H2_1(t) alpha*H - @compound H2_2(t) 2H + @compound H2_1(t) = alpha*H + @compound H2_2(t) = 2H @test iscompound(H2_1) @test iscompound(H2_2) @@ -116,8 +117,8 @@ let @parameters alpha = 2 @species H(t) - @compound H2_1(t) alpha*H - @compound H2_2(t) 2H + @compound H2_1(t) = alpha*H + @compound H2_2(t) = 2H @test iscompound(H2_1) @test iscompound(H2_2) @@ -130,13 +131,58 @@ let @variables t @species A(t) B = A - @compound A2(t) 2A - @compound B2(t) 2B + @compound A2(t) = 2A + @compound B2(t) = 2B @test iscompound(A2) @test iscompound(B2) @test isequal(components(A2),components(B2)) @test isequal(coefficients(A2), coefficients(B2)) + @test isequal(component_coefficients(A2), component_coefficients(B2)) end +### Check @compounds Macro ### + +# Basic syntax. +let + @variables t + @species C(t) H(t) O(t) + @compound OH(t) = 1O + 1H + @compound C3H5OH3(t) = 3C + 5H + 3OH + + @compounds begin + OH_alt(t) = 1O + 1H + C3H5OH3_alt(t) = 3C + 5H + 3OH + end + + @test iscompound(OH_alt) + @test iscompound(C3H5OH3_alt) + + @test isequal(components(OH),components(OH_alt)) + @test isequal(coefficients(OH), coefficients(OH_alt)) + @test isequal(component_coefficients(OH), component_coefficients(OH_alt)) + @test isequal(components(C3H5OH3),components(C3H5OH3_alt)) + @test isequal(coefficients(C3H5OH3), coefficients(C3H5OH3_alt)) + @test isequal(component_coefficients(C3H5OH3), component_coefficients(C3H5OH3_alt)) +end + +# Interpolation +let + @variables t + @species s1(t) s2(t) s3(t) + s2_alt = s2 + s3_alt = s3 + + @compounds begin + comp(t) = s1 + s2 + 4s3 + comp_alt(t) = s1 + s2_alt + 4s3_alt + end + + @test iscompound(comp) + @test iscompound(comp_alt) + + @test isequal(components(comp),components(comp_alt)) + @test isequal(coefficients(comp), coefficients(comp_alt)) + @test isequal(component_coefficients(comp), component_coefficients(comp_alt)) +end \ No newline at end of file From 99237cf2a9640748fd63fdf3f819b5ad576e3ef6 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 10 Nov 2023 18:26:01 -0500 Subject: [PATCH 06/30] misc improvements --- .../chemistry_related_functionality.md | 20 +- src/chemistry_functionality.jl | 9 +- .../miscellaneous_tests/reaction_balancing.jl | 196 ++++++++++-------- 3 files changed, 119 insertions(+), 106 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index d3e9eee86f..4e21d83ae4 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -2,10 +2,10 @@ While Catalyst has primarily been designed around the modelling of biological systems, reaction network models are also common across chemistry. This section describes two types of functionality, that while of general interest, should be especially useful in the modelling of chemical systems. - The `@compound` option, allowing the user to designate that a specific species is composed of certain subspecies. -- The `balance_reaction` function enabling the user to balance a reactions so the same number of compounds occur on both sides. +- The `balance_reaction` function enabling the user to balance a reaction so the same number of components occur on both sides. ## Modelling with compound species -Defining compound species is currently only supported for [programmatic construction](@ref programmatic_CRN_construction) of reactions and reaction network models. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ compound species, consisting of 1 C species and 2 O species. First we create species corresponding to the components: +Defining compound species is currently only supported for [programmatic construction](@ref programmatic_CRN_construction) of reactions and reaction network models. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ compound species, consisting of one C species and two O species. First, we create species corresponding to the components: ```@example chem1 @variables t @species C(t) O(t) @@ -14,11 +14,11 @@ Next, we create the `CO2` compound: ```@example chem1 @compound CO2(t) = 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 variable. Components with non-unitary stoichiometry have this value written before the component (generally, the rules for designating the components of a compounds are identical to hose 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` could be used: +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: ```@example chem1 isspecies(CO2) ``` -However, in its metadata is stored the information of its components, which can be retrieved using the `components` (returning a vector of its component species) and `coefficients` (returning a vector with each component's stoichiometry) functions: +In its metadata, however, is stored information of its components, which can be retrieved using the `components` (returning a vector of its component species) and `coefficients` (returning a vector with each component's stoichiometry) functions: ```@example chem1 components(CO2) ``` @@ -41,7 +41,7 @@ Compound components that are also compounds are allowed, e.g. we can create a ca @compound H2CO3(t) = CO2 + H2O ``` -When multiple compounds are created, they can be created simultaneously using the `@compounds` macro, e.g. the previous code-block could have been executed using: +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 @@ -50,18 +50,18 @@ When multiple compounds are created, they can be created simultaneously using th end ``` -One use defining a species as a compound is that they can be used to balance reactions to that the number of compounds are the same on both side. +One use of defining a species as a compound is that they can be used to balance reactions to that the number of compounds are the same on both sides. ## Balancing chemical reactions -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` enable us to find coefficients creating a balanced reaction (in this case, where the number of carbon and oxygen atoms are teh same on both sides). To demonstrate, we first created the unbalanced reactions: +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 rx = @reaction k, C + O --> $CO2 ``` -Here, we set a reaction rate `k` (which is not involved in the reaction balancing). We also use interpolation of `CO2`, ensuring that the `CO2` used in the reaction is the same one we previously defined as a compound of `C` and `O`. Next, we call the `balance_reaction` function +Here, the reaction rate (`k`) is not involved in the reaction balancing. We use interpolation for `CO2`, ensuring that the `CO2` used in the reaction is the same one we previously defined as a compound of `C` and `O`. Next, we call the `balance_reaction` function ```@example chem1 balance_reaction(rx) ``` -which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). +which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). No, or an infinite number of, solutions is also possible. Let us consider a more elaborate example, the reaction between ammonia (NH₃) and oxygen (O₂) to form nitrogen monoxide (NO) and water (H₂O). Let us first create the components and the unbalanced reaction: ```@example chem2 @@ -76,7 +76,7 @@ using Catalyst # hide end unbalanced_reaction = @reaction k, $NH3 + $O2 --> $NO + $H2O ``` -We can now created a balanced version (where the amount of H, N, and O is the same on both sides): +We can now create a balanced version (where the amount of H, N, and O is the same on both sides): ```@example chem2 balanced_reaction = balance_reaction(unbalanced_reaction)[1] ``` diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index ee1f809f48..afd12397c8 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -12,7 +12,7 @@ Symbolics.option_to_metadata_type(::Val{:coefficients}) = CompoundCoefficients """ @compound -Macro that creates a compound species, which is composed of smaller constituent species. +Macro that creates a compound species, which is composed of smaller component species. Example: ```julia @@ -22,7 +22,7 @@ Example: ``` Notes: -- The constituent species must be defined before using the `@compound` macro. +- The component species must be defined before using the `@compound` macro. """ macro compound(expr) make_compound(MacroTools.striplines(expr)) @@ -62,7 +62,7 @@ end """ @compounds -Macro that creates several compound species, which each is composed of smaller constituent species. Uses the same syntax as `@compound`, but with one compound species one each line. +Macro that creates several compound species, which each is composed of smaller component species. Uses the same syntax as `@compound`, but with one compound species one each line. Example: ```julia @@ -77,7 +77,7 @@ end ``` Notes: -- The constituent species must be defined before using the `@compound` macro. +- The component species must be defined before using the `@compound` macro. """ macro compounds(expr) make_compounds(MacroTools.striplines(expr)) @@ -132,6 +132,7 @@ function component_coefficients(s) return [c => co for (c, co) in zip(components(s), coefficients(s))] end + ### Reaction Balancing Functionality ### # Reaction balancing error. diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index b858cf7aa9..41e7ebd201 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -5,9 +5,11 @@ let @variables t @parameters k @species H(t) O(t) - @compound H2(t) 2H - @compound O2(t) 2O - @compound H2O(t) 2H 1O + @compounds begin + H2(t) = 2H + O2(t) = 2O + H2O(t) = 2H + 1O2 + end rx = Reaction(k,[H2,O2],[H2O]) @@ -18,16 +20,20 @@ let brxs = balance_reaction(rx) @test length(brxs) == 1 @test isequal(balanced_rx, first(brxs)) + + rx_macro = @reaction k, $H2 + $O2 --> $H2O + brxs_macro = balance_reaction(rx_macro) + @test isequal(brxs, brxs_macro) end let @variables t @parameters k - @species C(t) H(t) O(t) - @compound O2(t) 2O - @compound CO2(t) 1C 2O - @compound H2O(t) 2H 1O - @compound C6H12O6(t) 6C 12H 6O + @species C(t) H(t) = O(t) + @compound O2(t) = 2O + @compound CO2(t) = 1C + 2O + @compound H2O(t) = 2H + 1O + @compound C6H12O6(t) = 6C + 12H + 6O rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) @@ -38,13 +44,17 @@ let brxs = balance_reaction(rx) @test length(brxs) == 1 @test isequal(balanced_rx, first(brxs)) + + rx_macro = @reaction k, $CO2 + $H2O --> $C6H12O6 + $O2 + brxs_macro = balance_reaction(rx_macro) + @test isequal(brxs, brxs_macro) end # @reaction k, H2O --> H2O let @variables t @species H(t) O(t) - @compound H2O(t) 2H O + @compound H2O(t) = 2H + O rx = Reaction(1.0, [H2O], [H2O], [2], [2]) @@ -58,7 +68,7 @@ end let @variables t @species H(t) O(t) - @compound H2O(t) 2H O + @compound H2O(t) = 2H + O rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7]) @@ -72,10 +82,12 @@ end let @variables t @species H(t) O(t) C(t) - @compound CH4(t) C 4H - @compound O2(t) 2O - @compound CO2(t) C 2O - @compound H2O(t) 2H O + @compounds begin + CH4(t) = C + 4H + O2(t) = 2O + CO2(t) = C + 2O + H2O(t) = 2H + O + end rx = Reaction(1.0, [CH4, O2], [CO2, H2O]) @@ -89,9 +101,9 @@ end let @variables t @species H(t) N(t) - @compound N2(t) 2N - @compound H2(t) 2H - @compound NH3(t) N 3H + @compound N2(t) = 2N + @compound H2(t) = 2H + @compound NH3(t) = N + 3H rx = Reaction(1.0, [N2, H2], [NH3]) @@ -105,10 +117,10 @@ end let @variables t @species C(t) H(t) O(t) - @compound C2H5OH(t) 2C 6H O - @compound CH3COOH(t) 2C 4H 2O - @compound C4H8O2(t) 4C 8H 2O - @compound H2O(t) 2H O + @compound C2H5OH(t) = 2C + 6H + O + @compound CH3COOH(t) = 2C + 4H + 2O + @compound C4H8O2(t) = 4C + 8H + 2O + @compound H2O(t) = 2H + O rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O]) @@ -122,9 +134,9 @@ end let @variables t @species Ca(t) P(t) O(t) - @compound Ca3PO42(t) 3Ca 2P 8O - @compound CaO(t) Ca O - @compound P4O10(t) 4P 10O + @compound Ca3PO42(t) = 3Ca + 2P + 8O + @compound CaO(t) = Ca + O + @compound P4O10(t) = 4P + 10O rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10]) @@ -138,9 +150,9 @@ end let @variables t @species Fe(t) O(t) H(t) - @compound O2(t) 2O - @compound H2O(t) 2H O - @compound FeOH3(t) Fe 3H 3O + @compound O2(t) = 2O + @compound H2O(t) = 2H + O + @compound FeOH3(t) = Fe + 3H + 3O rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3]) @@ -154,11 +166,11 @@ end let @variables t @species Na(t) O(t) H(t) S(t) - @compound SO4(t) S 4O - @compound NaOH(t) Na O H - @compound H2SO4(t) 2H 1S 4O - @compound Na2SO4(t) 2Na 1S 4O - @compound H2O(t) 2H O + @compound SO4(t) = S + 4O + @compound NaOH(t) = Na + O + H + @compound H2SO4(t) = 2H + 1S + 4O + @compound Na2SO4(t) = 2Na + 1S + 4O + @compound H2O(t) = 2H + O rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O]) @@ -172,8 +184,8 @@ end let @variables t @species N(t) O(t) - @compound NO2(t) N 2O - @compound N2O4(t) 2N 4O + @compound NO2(t) = N + 2O + @compound N2O4(t) = 2N + 4O rx = Reaction(1.0, [NO2], [N2O4]) @@ -187,11 +199,11 @@ end let @variables t @species C(t) H(t) O(t) Ca(t) Cl(t) - @compound H2O(t) 2H 1O - @compound CO2(t) 1C 2O - @compound CaCO3(t) 1Ca 1C 3O - @compound HCl(t) 1H 1Cl - @compound CaCl2(t) 1Ca 2Cl + @compound H2O(t) = 2H + 1O + @compound CO2(t) = 1C + 2O + @compound CaCO3(t) = 1Ca + 1C + 3O + @compound HCl(t) = 1H + 1Cl + @compound CaCl2(t) = 1Ca + 2Cl rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O]) balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1]) @@ -204,10 +216,10 @@ end let @variables t @species Si(t) Cl(t) H(t) O(t) - @compound SiCl4(t) 1Si 4Cl - @compound H2O(t) 2H O - @compound H4SiO4(t) 4H Si 4O - @compound HCl(t) H Cl + @compound SiCl4(t) = 1Si + 4Cl + @compound H2O(t) = 2H + O + @compound H4SiO4(t) = 4H + Si + 4O + @compound HCl(t) = H + Cl rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl]) balanced_rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl], [1,4], [1,4]) @@ -220,9 +232,9 @@ end let @variables t @species Al(t) Cl(t) H(t) - @compound HCl(t) H Cl - @compound AlCl3(t) Al 3Cl - @compound H2(t) 2H + @compound HCl(t) = H + Cl + @compound AlCl3(t) = Al + 3Cl + @compound H2(t) = 2H rx = Reaction(1.0,[Al,HCl],[AlCl3,H2]) balanced_rx = Reaction(1.0,[Al,HCl],[AlCl3,H2],[2,6], [2,3]) @@ -235,11 +247,11 @@ end let @variables t @species Na(t) C(t) O(t) H(t) Cl(t) - @compound Na2CO3(t) 2Na C 3O - @compound HCl(t) H Cl - @compound NaCl(t) Na Cl - @compound H2O(t) 2H O - @compound CO2(t) C 2O + @compound Na2CO3(t) = 2Na + C + 3O + @compound HCl(t) = H + Cl + @compound NaCl(t) = Na + Cl + @compound H2O(t) = 2H + O + @compound CO2(t) = C + 2O rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2]) balanced_rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2], [1,2], [2,1,1]) @@ -252,10 +264,10 @@ end let @variables t @species C(t) H(t) O(t) - @compound C7H6O2(t) 7C 6H 2O - @compound O2(t) 2O - @compound CO2(t) C 2O - @compound H2O(t) 2H O + @compound C7H6O2(t) = 7C + 6H + 2O + @compound O2(t) = 2O + @compound CO2(t) = C + 2O + @compound H2O(t) = 2H + O rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O]) balanced_rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O], [2,15], [14,6]) @@ -268,10 +280,10 @@ end let @variables t @species Fe(t) S(t) O(t) H(t) K(t) - @compound Fe2S3O12(t) 2Fe 3S 12O - @compound KOH(t) K O H - @compound K2SO4(t) 2K S 4O - @compound FeO3H3(t) Fe 3O 3H + @compound Fe2S3O12(t) = 2Fe + 3S + 12O + @compound KOH(t) = K + O + H + @compound K2SO4(t) = 2K + S + 4O + @compound FeO3H3(t) = Fe + 3O + 3H rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3]) #5x4 matrix balanced_rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3], [1,6], [3,2]) @@ -284,10 +296,10 @@ end let @variables t @species Ca(t) P(t) O(t) Si(t) - @compound Ca3P2O8(t) 3Ca 2P 8O - @compound SiO2(t) Si 2O - @compound P4O10(t) 4P 10O - @compound CaSiO3(t) Ca Si 3O + @compound Ca3P2O8(t) = 3Ca + 2P + 8O + @compound SiO2(t) = Si + 2O + @compound P4O10(t) = 4P + 10O + @compound CaSiO3(t) = Ca + Si + 3O rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3]) #5x4 matrix balanced_rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3], [2,6] , [1,6]) @@ -300,9 +312,9 @@ end let @variables t @species K(t) Cl(t) O(t) - @compound KClO3(t) K Cl 3O - @compound KClO4(t) K Cl 4O - @compound KCl(t) K Cl + @compound KClO3(t) = K + Cl + 3O + @compound KClO4(t) = K + Cl + 4O + @compound KCl(t) = K + Cl rx = Reaction(1.0,[KClO3],[KClO4,KCl]) balanced_rx = Reaction(1.0,[KClO3],[KClO4,KCl], [4], [3,1]) @@ -315,10 +327,10 @@ end let @variables t @species Al(t) S(t) O(t) Ca(t) O(t) (H) - @compound Al2S3O12(t) 2Al 3S 12O - @compound CaO2H2(t) Ca 2O 2H - @compound AlO3H3(t) Al 3O 3H - @compound CaSO4(t) Ca S 4O + @compound Al2S3O12(t) = 2Al + 3S + 12O + @compound CaO2H2(t) = Ca + 2O + 2H + @compound AlO3H3(t) = Al + 3O + 3H + @compound CaSO4(t) = Ca + S + 4O rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4]) balanced_rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4], [1,3], [2,3]) @@ -331,11 +343,11 @@ end let @variables t @species H(t) S(t) O(t) I(t) - @compound H2SO4(t) 2H S 4O - @compound HI(t) H I - @compound H2S(t) 2H S - @compound I2(t) 2I - @compound H2O(t) 2H O + @compound H2SO4(t) = 2H + S + 4O + @compound HI(t) = H + I + @compound H2S(t) = 2H + S + @compound I2(t) = 2I + @compound H2O(t) = 2H + O rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O]) balanced_rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O], [1,8], [1,4,4]) @@ -348,10 +360,10 @@ end let @variables t @species C(t) H(t) O(t) - @compound C2H4(t) 2C 4H - @compound O2(t) 2O - @compound CO2(t) C 2O - @compound H2O(t) 2H O + @compound C2H4(t) = 2C + 4H + @compound O2(t) = 2O + @compound CO2(t) = C + 2O + @compound H2O(t) = 2H + O rx = Reaction(1.0,[C2H4,O2],[CO2,H2O]) balanced_rx = Reaction(1.0,[C2H4,O2],[CO2,H2O],[1,3],[2,2]) @@ -364,11 +376,11 @@ end let @variables t @species C(t) H(t) O(t) - @compound CO(t) C O - @compound CO2(t) C 2O - @compound H2(t) 2H - @compound CH4(t) C 4H - @compound H2O(t) 2H O + @compound CO(t) = C + O + @compound CO2(t) = C + 2O + @compound H2(t) = 2H + @compound CH4(t) = C + 4H + @compound H2O(t) = 2H + O rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O]) brxs = balance_reaction(rx) @@ -381,11 +393,11 @@ let @variables t @species Fe(t) S(t) O(t) H(t) N(t) - @compound FeS2(t) Fe 2S - @compound HNO3(t) H N 3O - @compound Fe2S3O12(t) 2Fe 3S 12O - @compound NO(t) N O - @compound H2SO4(t) 2H S 4O + @compound FeS2(t) = Fe + 2S + @compound HNO3(t) = H + N + 3O + @compound Fe2S3O12(t) = 2Fe + 3S + 12O + @compound NO(t) = N + O + @compound H2SO4(t) = 2H + S + 4O rx = Reaction(1.0,[FeS2,HNO3],[Fe2S3O12,NO,H2SO4]) brxs = balance_reaction(rx) @@ -397,9 +409,9 @@ end let @variables t @species C(t) H(t) O(t) - @compound CO(t) C O - @compound H2(t) 2H - @compound COH2(t) CO H2 + @compound CO(t) = C + O + @compound H2(t) = 2H + @compound COH2(t) = CO + H2 rx = Reaction(1.0, [CO, H2], [COH2]) @test_throws Catalyst.COMPOUND_OF_COMPOUND_ERROR balance_reaction(rx) From b0c027d1447c1a829bf57684544ee49584b68b4e Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 10 Nov 2023 20:03:02 -0500 Subject: [PATCH 07/30] @compounds in the macro --- src/chemistry_functionality.jl | 38 ++++++-- src/reaction_network.jl | 24 ++++- test/miscellaneous_tests/compound_macro.jl | 92 ++++++++++++++++++- .../miscellaneous_tests/reaction_balancing.jl | 26 +++++- 4 files changed, 165 insertions(+), 15 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index afd12397c8..3b2de11fe1 100644 --- a/src/chemistry_functionality.jl +++ b/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 ### diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 6a5fb6871c..1b40f66a93 100644 --- a/src/reaction_network.jl +++ b/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) diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index d13c72cd7b..b06279c7fe 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/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 \ No newline at end of file diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index 41e7ebd201..0cce80e40d 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/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 \ No newline at end of file From 8f29ced18011bf34b250a792f8883e710219f33a Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 10 Nov 2023 20:11:44 -0500 Subject: [PATCH 08/30] update doc with DSL compounds option. --- .../chemistry_related_functionality.md | 34 +++++++++++++++++-- 1 file changed, 31 insertions(+), 3 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 4e21d83ae4..f55c174cd0 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -5,7 +5,9 @@ While Catalyst has primarily been designed around the modelling of biological sy - The `balance_reaction` function enabling the user to balance a reaction so the same number of components occur on both sides. ## Modelling with compound species -Defining compound species is currently only supported for [programmatic construction](@ref programmatic_CRN_construction) of reactions and reaction network models. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ compound species, consisting of one C species and two O species. First, we create species corresponding to the components: + +#### Creating compound species programmatically +We will first show how to create compound species through [programmatic construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ compound species, consisting of one C species and two O species. First, we create species corresponding to the components: ```@example chem1 @variables t @species C(t) O(t) @@ -50,10 +52,34 @@ When multiple compounds are created, they can be created simultaneously using th end ``` -One use of defining a species as a compound is that they can be used to balance reactions to that the number of compounds are the same on both sides. +#### Creating compound species programmatically +It is also possible to declare species as compound species within the `@reaction_network` DSL, using the `@compounds` options: +```@example chem1 +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 + end + (k1,k2), H2O+ CO2 <--> H2CO3 +end +``` +When creating compound species using the DSL, it is important to note that *every component must be known to the system as a species, either by being declared using the `@species` option, or by appearing in a reaction*. E.g. the following is not valid +```julia +rn = @reaction_network begin + @compounds begin + C2O(t) = C + 2O + H2O(t) = 2H + O + H2CO3(t) = CO2 + H2O + end + (k1,k2), H2O+ CO2 <--> H2CO3 +end +``` +as the components `C`, `H`, and `O` are not explicitly declared as a species anywhere. Please also note that only `@compounds` can be used as an option in the DSL, not `@compound`. ## Balancing chemical reactions -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: +One use of defining a species as a compound is that they can be used to balance reactions to that the number of compounds 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 rx = @reaction k, C + O --> $CO2 ``` @@ -81,5 +107,7 @@ We can now create a balanced version (where the amount of H, N, and O is the sam balanced_reaction = balance_reaction(unbalanced_reaction)[1] ``` +Reactions declared as a part of a `ReactionSystem` (e.g. using the DSL) can be retrieved for balancing using the `reaction` function. Please note that balancing these will not mutate the `ReactionSystem`, but a new reaction system will need to be created using the balanced reactions. + !!! note Reaction balancing is currently not supported for reactions involving compounds of compounds. \ No newline at end of file From 55a9fa3392c3fd329f5ec8e4cf4539bbd9046484 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 10 Nov 2023 20:25:04 -0500 Subject: [PATCH 09/30] update --- .../chemistry_related_functionality.md | 8 ++++---- src/chemistry_functionality.jl | 8 -------- src/reaction_network.jl | 4 ++-- test/miscellaneous_tests/compound_macro.jl | 2 +- 4 files changed, 7 insertions(+), 15 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index f55c174cd0..5f2366d7bd 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -7,12 +7,12 @@ While Catalyst has primarily been designed around the modelling of biological sy ## Modelling with compound species #### Creating compound species programmatically -We will first show how to create compound species through [programmatic construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ compound species, consisting of one C species and two O species. First, we create species corresponding to the components: +We will first show how to create compound species through [programmatic construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ molecule, consisting of one C atom and two O atoms. First, we create species corresponding to the components: ```@example chem1 @variables t @species C(t) O(t) ``` -Next, we create the `CO2` compound: +Next, we create the `CO2` compound species: ```@example chem1 @compound CO2(t) = C + 2O ``` @@ -76,10 +76,10 @@ rn = @reaction_network begin (k1,k2), H2O+ CO2 <--> H2CO3 end ``` -as the components `C`, `H`, and `O` are not explicitly declared as a species anywhere. Please also note that only `@compounds` can be used as an option in the DSL, not `@compound`. +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`. ## 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 compounds 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: +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 rx = @reaction k, C + O --> $CO2 ``` diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 3b2de11fe1..3c6fc20de9 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -55,14 +55,6 @@ 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 diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 1b40f66a93..322ee5e68d 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -468,11 +468,11 @@ 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 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. - else + else # If option is not used, return empty vectors and expressions. compound_expr = :() compound_species = Union{Symbol, Expr}[] end diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index b06279c7fe..50e9677721 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -106,7 +106,7 @@ let @compound H2_1(t) = 2*H @compound H2_2(t) = alpha*H @compound H2_3(t) = 2*h - @compound H2_4(t) = alpha*2H + @compound H2_4(t) = alpha*H @test iscompound(H2_1) @test iscompound(H2_2) From 28757242560f2f0a988fcf643feba8da1df39cf9 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 10 Nov 2023 20:27:13 -0500 Subject: [PATCH 10/30] add API in docs --- docs/src/api/catalyst_api.md | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md index 2aa153c3e7..a6c48d03cc 100644 --- a/docs/src/api/catalyst_api.md +++ b/docs/src/api/catalyst_api.md @@ -280,6 +280,17 @@ Base.convert ModelingToolkit.structural_simplify ``` +## Chemistry-related functionalities +Various functionalities primarily relevant to modelling of chemical systems (but potentially also in biology). +```@docs +@compound +@compounds +iscompound +components +coefficients +component_coefficients +``` + ## Unit validation ```@docs validate(rx::Reaction; info::String = "") From f6906b6b06de299c5a39f235f266a17b6cb58c64 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 11 Nov 2023 14:02:22 -0500 Subject: [PATCH 11/30] up --- .../catalyst_functionality/chemistry_related_functionality.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 5f2366d7bd..b40f29eea1 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -9,6 +9,7 @@ While Catalyst has primarily been designed around the modelling of biological sy #### Creating compound species programmatically We will first show how to create compound species through [programmatic construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ molecule, consisting of one C atom and two O atoms. First, we create species corresponding to the components: ```@example chem1 +using Catalyst @variables t @species C(t) O(t) ``` From 713962576d8a766e8d5934133410b7dad6151322 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 13 Nov 2023 21:36:32 -0500 Subject: [PATCH 12/30] up --- .../chemistry_related_functionality.md | 46 +++--- src/chemistry_functionality.jl | 29 +++- src/expression_utils.jl | 67 +++++++++ src/reaction_network.jl | 3 +- test/miscellaneous_tests/compound_macro.jl | 131 +++++++++++++----- 5 files changed, 221 insertions(+), 55 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index b40f29eea1..672f0e8b3d 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/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 ``` diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 3c6fc20de9..bfc9257f42 100644 --- a/src/chemistry_functionality.jl +++ b/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 diff --git a/src/expression_utils.jl b/src/expression_utils.jl index 2778772e33..c5731becda 100644 --- a/src/expression_utils.jl +++ b/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)) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 322ee5e68d..c214f634b3 100644 --- a/src/reaction_network.jl +++ b/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}[] diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index 50e9677721..b6ee23069e 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -1,10 +1,79 @@ using Catalyst, Test +### Tests Main Macro Creation Forms ### +let + @variables t + @species C(t) H(t) O(t) + @parameters p1 p2 + + # Basic cases that should pass: + @compound H2O_1 ~ 2H + O + @compound H2O_2, [output=true] ~ 2H + O + @compound (H2O_3 = 1.5) ~ 2H + O + @compound (H2O_4 = 4, [output=true]) ~ 2H + O + @compound (H2O_5 = p1, [output=true]) ~ 2H + p2*O + @test iscompound(H2O_1) + @test iscompound(H2O_2) + @test iscompound(H2O_3) + @test iscompound(H2O_4) + @test iscompound(H2O_5) + + # Independent variable given. + @test_throws LoadError @eval @compound H2O(t) ~ 2H + O + @test_throws LoadError @eval @compound H2O(t), [output=true] ~ 2H + O + @test_throws LoadError @eval @compound (H2O(t) = 1.5) ~ 2H + O + @test_throws LoadError @eval @compound (H2O(t) = 4, [output=true]) ~ 2H + O + @test_throws LoadError @eval @compound (H2O(t) = p1, [output=true]) ~ 2H + p2*O + + # Other errors. + @test_throws LoadError @eval @compound H2O(t) = 2H + O + @test_throws LoadError @eval @compound H2O(t), [output=true] = 2H + O + @test_throws LoadError @eval @compound H2O(t) = 1.5 ~ 2H + O + @test_throws LoadError @eval @compound H2O(t) = 4, [output=true] ~ 2H + O + @test_throws LoadError @eval @compound H2O(t) = p1, [output=true] ~ 2H + p2*O + + # Compounds created in block notation. + @compounds begin + CO2_1 ~ 2H + O + end + @compounds begin + CO2_2, [output=true] ~ 2H + O + (CO2_3 = 1.5) ~ 2H + O + (CO2_4 = 4, [output=true]) ~ 2H + O + (CO2_5 = p1, [output=true]) ~ 2H + p2*O + end + @test iscompound(CO2_1) + @test iscompound(CO2_2) + @test iscompound(CO2_3) + @test iscompound(CO2_4) + @test iscompound(CO2_5) + + # Declares stuff in the DSL. + rn = @reaction_network begin + @species N(t) H(t) + @parameters p1 p2 + @compounds begin + NH3_1 ~ N + 3H + NH3_2, [output=true] ~ N + 3H + (NH3_3 = 1.5) ~ N + 3H + (NH3_4 = 4, [output=true]) ~ N + 3H + (NH3_5 = p1, [output=true]) ~ N + p2*H + end + end + @test iscompound(rn.NH3_1) + @test iscompound(rn.NH3_2) + @test iscompound(rn.NH3_3) + @test iscompound(rn.NH3_4) + @test iscompound(rn.NH3_5) +end + +### Other Minor Tests ### + # Test base functionality in two cases. let @variables t @species C(t) H(t) O(t) - @compound C6H12O2(t) = 6C + 12H + 2O + @compound C6H12O2 ~ 6C + 12H + 2O @test iscompound(C6H12O2) @test isspecies(C6H12O2) @@ -21,7 +90,7 @@ end let @variables t @species O(t) - @compound O2(t) = 2O + @compound O2 ~ 2O @test iscompound(O2) @test isspecies(O2) @@ -37,19 +106,19 @@ end let @variables t @species C(t) H(t) - @test_throws Exception @compound C6H12O2(t) = 6C + 12H + 2O + @test_throws Exception @compound C6H12O2 ~ 6C + 12H + 2O end let @variables t - @test_throws Exception @compound O2(t) = 2O + @test_throws Exception @compound O2 ~ 2O end # Checks that nested components works as expected. let @variables t @species C(t) H(t) O(t) - @compound OH(t) = 1O + 1H - @compound C3H5OH3(t) = 3C + 5H + 3OH + @compound OH ~ 1O + 1H + @compound C3H5OH3 ~ 3C + 5H + 3OH @test !iscompound(O) @test !iscompound(H) @@ -71,8 +140,8 @@ let @variables t @species C(t) H(t) O(t) s = C - @compound C6H12O2_1(t) = 6s + 12H + 2O - @compound C6H12O2_2(t) = 6C + 12H + 2O + @compound C6H12O2_1 ~ 6s + 12H + 2O + @compound C6H12O2_2 ~ 6C + 12H + 2O @test iscompound(C6H12O2_1) @test iscompound(C6H12O2_2) @@ -85,9 +154,9 @@ end let @variables t @species C(t) H(t) - @compound Cyclopentadiene(t) = 5C + 6H + @compound Cyclopentadiene ~ 5C + 6H C5H6 = Cyclopentadiene - @compound C10H12(t) = 2C5H6 + @compound C10H12 ~ 2C5H6 @test iscompound(C10H12) @test iscompound(components(C10H12)[1]) @@ -103,10 +172,10 @@ let alpha = 2 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*H + @compound H2_1 ~ 2*H + @compound H2_2 ~ alpha*H + @compound H2_3 ~ 2*h + @compound H2_4 ~ alpha*H @test iscompound(H2_1) @test iscompound(H2_2) @@ -126,8 +195,8 @@ let @parameters alpha = 2 @species H(t) - @compound H2_1(t) = alpha*H - @compound H2_2(t) = 2H + @compound H2_1 ~ alpha*H + @compound H2_2 ~ 2H @test iscompound(H2_1) @test iscompound(H2_2) @@ -140,8 +209,8 @@ let @variables t @species A(t) B = A - @compound A2(t) = 2A - @compound B2(t) = 2B + @compound A2 ~ 2A + @compound B2 ~ 2B @test iscompound(A2) @test iscompound(B2) @@ -157,12 +226,12 @@ end let @variables t @species C(t) H(t) O(t) - @compound OH(t) = 1O + 1H - @compound C3H5OH3(t) = 3C + 5H + 3OH + @compound OH ~ 1O + 1H + @compound C3H5OH3 ~ 3C + 5H + 3OH @compounds begin - OH_alt(t) = 1O + 1H - C3H5OH3_alt(t) = 3C + 5H + 3OH + OH_alt ~ 1O + 1H + C3H5OH3_alt ~ 3C + 5H + 3OH end @test iscompound(OH_alt) @@ -184,8 +253,8 @@ let s3_alt = s3 @compounds begin - comp(t) = s1 + s2 + 4s3 - comp_alt(t) = s1 + s2_alt + 4s3_alt + comp ~ s1 + s2 + 4s3 + comp_alt ~ s1 + s2_alt + 4s3_alt end @test iscompound(comp) @@ -205,7 +274,7 @@ let rn = @reaction_network begin @species C(t) O(t) @compounds begin - CO2(t) = C + 2O + CO2 ~ C + 2O end end @unpack C, O, CO2 = rn @@ -223,10 +292,10 @@ 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 + CH4 ~ C + 4H + O2 ~ 2O + CO2 ~ C + 2O + H2O ~ 2H + O end k, CH4 + O2 --> CO2 + H2O end) @@ -252,8 +321,8 @@ end let rn = @reaction_network begin @compounds begin - SO2(t) = S + 2O - S2O4(t) = 2SO2 + SO2 ~ S + 2O + S2O4 ~ 2SO2 end dS, S --> 0 dO, O --> 0 From cc47e15e7f16b84e8a20eb9103726cf5fb3611da Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 13 Nov 2023 21:38:44 -0500 Subject: [PATCH 13/30] reaction balancing test new compound style --- .../miscellaneous_tests/reaction_balancing.jl | 190 +++++++++--------- 1 file changed, 95 insertions(+), 95 deletions(-) diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index 0cce80e40d..110b5318a6 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -6,9 +6,9 @@ let @parameters k @species H(t) O(t) @compounds begin - H2(t) = 2H - O2(t) = 2O - H2O(t) = 2H + 1O + H2 ~ 2H + O2 ~ 2O + H2O ~ 2H + 1O end rx = Reaction(k,[H2,O2],[H2O]) @@ -30,10 +30,10 @@ let @variables t @parameters k @species C(t) H(t) O(t) - @compound O2(t) = 2O - @compound CO2(t) = 1C + 2O - @compound H2O(t) = 2H + 1O - @compound C6H12O6(t) = 6C + 12H + 6O + @compound O2 ~ 2O + @compound CO2 ~ 1C + 2O + @compound H2O ~ 2H + 1O + @compound C6H12O6 ~ 6C + 12H + 6O rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) @@ -54,7 +54,7 @@ end let @variables t @species H(t) O(t) - @compound H2O(t) = 2H + O + @compound H2O ~ 2H + O rx = Reaction(1.0, [H2O], [H2O], [2], [2]) @@ -68,7 +68,7 @@ end let @variables t @species H(t) O(t) - @compound H2O(t) = 2H + O + @compound H2O ~ 2H + O rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7]) @@ -83,10 +83,10 @@ let @variables t @species H(t) O(t) C(t) @compounds begin - CH4(t) = C + 4H - O2(t) = 2O - CO2(t) = C + 2O - H2O(t) = 2H + O + CH4 ~ C + 4H + O2 ~ 2O + CO2 ~ C + 2O + H2O ~ 2H + O end rx = Reaction(1.0, [CH4, O2], [CO2, H2O]) @@ -101,9 +101,9 @@ end let @variables t @species H(t) N(t) - @compound N2(t) = 2N - @compound H2(t) = 2H - @compound NH3(t) = N + 3H + @compound N2 ~ 2N + @compound H2 ~ 2H + @compound NH3 ~ N + 3H rx = Reaction(1.0, [N2, H2], [NH3]) @@ -117,10 +117,10 @@ end let @variables t @species C(t) H(t) O(t) - @compound C2H5OH(t) = 2C + 6H + O - @compound CH3COOH(t) = 2C + 4H + 2O - @compound C4H8O2(t) = 4C + 8H + 2O - @compound H2O(t) = 2H + O + @compound C2H5OH ~ 2C + 6H + O + @compound CH3COOH ~ 2C + 4H + 2O + @compound C4H8O2 ~ 4C + 8H + 2O + @compound H2O ~ 2H + O rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O]) @@ -134,9 +134,9 @@ end let @variables t @species Ca(t) P(t) O(t) - @compound Ca3PO42(t) = 3Ca + 2P + 8O - @compound CaO(t) = Ca + O - @compound P4O10(t) = 4P + 10O + @compound Ca3PO42 ~ 3Ca + 2P + 8O + @compound CaO ~ Ca + O + @compound P4O10 ~ 4P + 10O rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10]) @@ -150,9 +150,9 @@ end let @variables t @species Fe(t) O(t) H(t) - @compound O2(t) = 2O - @compound H2O(t) = 2H + O - @compound FeOH3(t) = Fe + 3H + 3O + @compound O2 ~ 2O + @compound H2O ~ 2H + O + @compound FeOH3 ~ Fe + 3H + 3O rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3]) @@ -166,11 +166,11 @@ end let @variables t @species Na(t) O(t) H(t) S(t) - @compound SO4(t) = S + 4O - @compound NaOH(t) = Na + O + H - @compound H2SO4(t) = 2H + 1S + 4O - @compound Na2SO4(t) = 2Na + 1S + 4O - @compound H2O(t) = 2H + O + @compound SO4 ~ S + 4O + @compound NaOH ~ Na + O + H + @compound H2SO4 ~ 2H + 1S + 4O + @compound Na2SO4 ~ 2Na + 1S + 4O + @compound H2O ~ 2H + O rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O]) @@ -184,8 +184,8 @@ end let @variables t @species N(t) O(t) - @compound NO2(t) = N + 2O - @compound N2O4(t) = 2N + 4O + @compound NO2 ~ N + 2O + @compound N2O4 ~ 2N + 4O rx = Reaction(1.0, [NO2], [N2O4]) @@ -199,11 +199,11 @@ end let @variables t @species C(t) H(t) O(t) Ca(t) Cl(t) - @compound H2O(t) = 2H + 1O - @compound CO2(t) = 1C + 2O - @compound CaCO3(t) = 1Ca + 1C + 3O - @compound HCl(t) = 1H + 1Cl - @compound CaCl2(t) = 1Ca + 2Cl + @compound H2O ~ 2H + 1O + @compound CO2 ~ 1C + 2O + @compound CaCO3 ~ 1Ca + 1C + 3O + @compound HCl ~ 1H + 1Cl + @compound CaCl2 ~ 1Ca + 2Cl rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O]) balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1]) @@ -216,10 +216,10 @@ end let @variables t @species Si(t) Cl(t) H(t) O(t) - @compound SiCl4(t) = 1Si + 4Cl - @compound H2O(t) = 2H + O - @compound H4SiO4(t) = 4H + Si + 4O - @compound HCl(t) = H + Cl + @compound SiCl4 ~ 1Si + 4Cl + @compound H2O ~ 2H + O + @compound H4SiO4 ~ 4H + Si + 4O + @compound HCl ~ H + Cl rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl]) balanced_rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl], [1,4], [1,4]) @@ -232,9 +232,9 @@ end let @variables t @species Al(t) Cl(t) H(t) - @compound HCl(t) = H + Cl - @compound AlCl3(t) = Al + 3Cl - @compound H2(t) = 2H + @compound HCl ~ H + Cl + @compound AlCl3 ~ Al + 3Cl + @compound H2 ~ 2H rx = Reaction(1.0,[Al,HCl],[AlCl3,H2]) balanced_rx = Reaction(1.0,[Al,HCl],[AlCl3,H2],[2,6], [2,3]) @@ -247,11 +247,11 @@ end let @variables t @species Na(t) C(t) O(t) H(t) Cl(t) - @compound Na2CO3(t) = 2Na + C + 3O - @compound HCl(t) = H + Cl - @compound NaCl(t) = Na + Cl - @compound H2O(t) = 2H + O - @compound CO2(t) = C + 2O + @compound Na2CO3 ~ 2Na + C + 3O + @compound HCl ~ H + Cl + @compound NaCl ~ Na + Cl + @compound H2O ~ 2H + O + @compound CO2 ~ C + 2O rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2]) balanced_rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2], [1,2], [2,1,1]) @@ -264,10 +264,10 @@ end let @variables t @species C(t) H(t) O(t) - @compound C7H6O2(t) = 7C + 6H + 2O - @compound O2(t) = 2O - @compound CO2(t) = C + 2O - @compound H2O(t) = 2H + O + @compound C7H6O2 ~ 7C + 6H + 2O + @compound O2 ~ 2O + @compound CO2 ~ C + 2O + @compound H2O ~ 2H + O rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O]) balanced_rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O], [2,15], [14,6]) @@ -280,10 +280,10 @@ end let @variables t @species Fe(t) S(t) O(t) H(t) K(t) - @compound Fe2S3O12(t) = 2Fe + 3S + 12O - @compound KOH(t) = K + O + H - @compound K2SO4(t) = 2K + S + 4O - @compound FeO3H3(t) = Fe + 3O + 3H + @compound Fe2S3O12 ~ 2Fe + 3S + 12O + @compound KOH ~ K + O + H + @compound K2SO4 ~ 2K + S + 4O + @compound FeO3H3 ~ Fe + 3O + 3H rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3]) #5x4 matrix balanced_rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3], [1,6], [3,2]) @@ -296,10 +296,10 @@ end let @variables t @species Ca(t) P(t) O(t) Si(t) - @compound Ca3P2O8(t) = 3Ca + 2P + 8O - @compound SiO2(t) = Si + 2O - @compound P4O10(t) = 4P + 10O - @compound CaSiO3(t) = Ca + Si + 3O + @compound Ca3P2O8 ~ 3Ca + 2P + 8O + @compound SiO2 ~ Si + 2O + @compound P4O10 ~ 4P + 10O + @compound CaSiO3 ~ Ca + Si + 3O rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3]) #5x4 matrix balanced_rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3], [2,6] , [1,6]) @@ -312,9 +312,9 @@ end let @variables t @species K(t) Cl(t) O(t) - @compound KClO3(t) = K + Cl + 3O - @compound KClO4(t) = K + Cl + 4O - @compound KCl(t) = K + Cl + @compound KClO3 ~ K + Cl + 3O + @compound KClO4 ~ K + Cl + 4O + @compound KCl ~ K + Cl rx = Reaction(1.0,[KClO3],[KClO4,KCl]) balanced_rx = Reaction(1.0,[KClO3],[KClO4,KCl], [4], [3,1]) @@ -327,10 +327,10 @@ end let @variables t @species Al(t) S(t) O(t) Ca(t) O(t) (H) - @compound Al2S3O12(t) = 2Al + 3S + 12O - @compound CaO2H2(t) = Ca + 2O + 2H - @compound AlO3H3(t) = Al + 3O + 3H - @compound CaSO4(t) = Ca + S + 4O + @compound Al2S3O12 ~ 2Al + 3S + 12O + @compound CaO2H2 ~ Ca + 2O + 2H + @compound AlO3H3 ~ Al + 3O + 3H + @compound CaSO4 ~ Ca + S + 4O rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4]) balanced_rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4], [1,3], [2,3]) @@ -343,11 +343,11 @@ end let @variables t @species H(t) S(t) O(t) I(t) - @compound H2SO4(t) = 2H + S + 4O - @compound HI(t) = H + I - @compound H2S(t) = 2H + S - @compound I2(t) = 2I - @compound H2O(t) = 2H + O + @compound H2SO4 ~ 2H + S + 4O + @compound HI ~ H + I + @compound H2S ~ 2H + S + @compound I2 ~ 2I + @compound H2O ~ 2H + O rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O]) balanced_rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O], [1,8], [1,4,4]) @@ -360,10 +360,10 @@ end let @variables t @species C(t) H(t) O(t) - @compound C2H4(t) = 2C + 4H - @compound O2(t) = 2O - @compound CO2(t) = C + 2O - @compound H2O(t) = 2H + O + @compound C2H4 ~ 2C + 4H + @compound O2 ~ 2O + @compound CO2 ~ C + 2O + @compound H2O ~ 2H + O rx = Reaction(1.0,[C2H4,O2],[CO2,H2O]) balanced_rx = Reaction(1.0,[C2H4,O2],[CO2,H2O],[1,3],[2,2]) @@ -376,11 +376,11 @@ end let @variables t @species C(t) H(t) O(t) - @compound CO(t) = C + O - @compound CO2(t) = C + 2O - @compound H2(t) = 2H - @compound CH4(t) = C + 4H - @compound H2O(t) = 2H + O + @compound CO ~ C + O + @compound CO2 ~ C + 2O + @compound H2 ~ 2H + @compound CH4 ~ C + 4H + @compound H2O ~ 2H + O rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O]) brxs = balance_reaction(rx) @@ -393,11 +393,11 @@ let @variables t @species Fe(t) S(t) O(t) H(t) N(t) - @compound FeS2(t) = Fe + 2S - @compound HNO3(t) = H + N + 3O - @compound Fe2S3O12(t) = 2Fe + 3S + 12O - @compound NO(t) = N + O - @compound H2SO4(t) = 2H + S + 4O + @compound FeS2 ~ Fe + 2S + @compound HNO3 ~ H + N + 3O + @compound Fe2S3O12 ~ 2Fe + 3S + 12O + @compound NO ~ N + O + @compound H2SO4 ~ 2H + S + 4O rx = Reaction(1.0,[FeS2,HNO3],[Fe2S3O12,NO,H2SO4]) brxs = balance_reaction(rx) @@ -409,9 +409,9 @@ end let @variables t @species C(t) H(t) O(t) - @compound CO(t) = C + O - @compound H2(t) = 2H - @compound COH2(t) = CO + H2 + @compound CO ~ C + O + @compound H2 ~ 2H + @compound COH2 ~ CO + H2 rx = Reaction(1.0, [CO, H2], [COH2]) @test_throws Catalyst.COMPOUND_OF_COMPOUND_ERROR balance_reaction(rx) @@ -422,10 +422,10 @@ 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 + O2 ~ 2O + CO2 ~ 1C + 2O + H2O ~ 2H + 1O + C6H12O6 ~ 6C + 12H + 6O end k, CO2 + H2O --> C6H12O6 + O2 end) From c1e8e45c4439200bab75664a62662ce31c601efc Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 14 Nov 2023 16:54:01 -0500 Subject: [PATCH 14/30] update --- .../chemistry_related_functionality.md | 10 ++++++++ src/chemistry_functionality.jl | 9 +++++++- src/expression_utils.jl | 12 ++++++---- test/miscellaneous_tests/compound_macro.jl | 23 +++++++++++++++++++ 4 files changed, 49 insertions(+), 5 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 672f0e8b3d..00d2f3407a 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -37,6 +37,16 @@ Finally, it is possible to check whether a species is a compound or not using th iscompound(CO2) ``` +!!! note + Currently, compounds with components that depend on variables that are not `t` are not supported. E.g. + ```julia + @variables x y + @species O(x, y) + @compound O2 ~ 2O + ``` + will currently throw an error. + + 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) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index bfc9257f42..b13aac9853 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -47,7 +47,7 @@ function make_compound(expr) # Extracts the composite species name, and a Vector{ReactantStruct} of its components. 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_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)) @@ -61,13 +61,20 @@ function make_compound(expr) # 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). + non_t_iv_error_check_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)` species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)` compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)` 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])` + # Currently, non-t independent variables are not supported for compounds. If there are any like these, we throw an error: + non_t_iv_error_check_expr = Expr(:escape, :(issetequal(unique(reduce(vcat, arguments.(ModelingToolkit.unwrap.($components)))), [t]) || error("Currently, compounds depending on components that are not \"t\" are not supported."))) + +println(non_t_iv_error_check_expr) + # Returns the rephrased expression. return quote + $non_t_iv_error_check_expr $species_declaration_expr $compound_designation_expr $components_designation_expr diff --git a/src/expression_utils.jl b/src/expression_utils.jl index c5731becda..379fd305b1 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -50,19 +50,23 @@ end # X(t), [metadata=true] # X(t) = 1.0, [metadata=true] # (In this example the independent variable t was inserted). +# Extra dispatch is to ensure so that iv is a vector (so that we can handle both single or multiple iv insertion in one function). +# - This way we don't have to make a check for how to create the final expression (which is different for symbol or vector of symbols). function insert_independent_variable(expr_in, iv) + # If iv is a symbol (e.g. :t), we convert to vector (e.g. [:t]). This way we can handle single and multiple ivs using the same code, without needing to check at the end. + (iv isa Symbol) && (iv = [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_in isa Symbol) && (return Expr(:call, expr_in, iv...)) expr = deepcopy(expr_in) if expr.head == :(=) # Case: :(X = 1.0) - expr.args[1] = :($(expr.args[1])($iv)) + expr.args[1] = Expr(:call, 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)) + expr.args[1] = Expr(:call, 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)) + expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], iv...) end end (expr == expr_in) && error("Failed to add independent variable $(iv) to expression: $expr_in") diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index b6ee23069e..a0df95f423 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -67,6 +67,29 @@ let @test iscompound(rn.NH3_5) end +### Test Various Independent Variables ### +let + @variables t x y z + @species C(t) H(x) N(x) O(t) P(t,x) S(x,y) + + @compound CO2 ~ C + 2O + @test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t]) + + # These 4 tests checks currently non-supported feature. + @test_broken false + @test_broken false + @test_broken false + @test_broken false + # @compound NH4, [output=true] ~ N + 4H + # @compound (H2O = 2.0) ~ 2H + N + # @compound PH4 ~ P + 4H + # @compound SO2 ~ S + 2O + # @test issetequal(arguments(ModelingToolkit.unwrap(NH4)), [x]) + # @test issetequal(arguments(ModelingToolkit.unwrap(H2O)), [t, x]) + # @test issetequal(arguments(ModelingToolkit.unwrap(PH4)), [t, x]) + # @test issetequal(arguments(ModelingToolkit.unwrap(SO2)), [t, x, y]) +end + ### Other Minor Tests ### # Test base functionality in two cases. From 6e6f43ba7a5d896868d7322a75f27fdee36c5540 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 14 Nov 2023 18:40:34 -0500 Subject: [PATCH 15/30] update --- src/chemistry_functionality.jl | 2 -- test/miscellaneous_tests/reaction_balancing.jl | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index b13aac9853..64d9ee7b56 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -70,8 +70,6 @@ function make_compound(expr) # Currently, non-t independent variables are not supported for compounds. If there are any like these, we throw an error: non_t_iv_error_check_expr = Expr(:escape, :(issetequal(unique(reduce(vcat, arguments.(ModelingToolkit.unwrap.($components)))), [t]) || error("Currently, compounds depending on components that are not \"t\" are not supported."))) -println(non_t_iv_error_check_expr) - # Returns the rephrased expression. return quote $non_t_iv_error_check_expr diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index 110b5318a6..2a5194a080 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -326,7 +326,7 @@ end # @reaction k, Al2(SO4)3 + 3Ca(OH)2 → 2Al(OH)3 + 3CaSO4 let @variables t - @species Al(t) S(t) O(t) Ca(t) O(t) (H) + @species Al(t) S(t) O(t) Ca(t) O(t) H(t) @compound Al2S3O12 ~ 2Al + 3S + 12O @compound CaO2H2 ~ Ca + 2O + 2H @compound AlO3H3 ~ Al + 3O + 3H From ed705d4165fd6d9683097d5315f3553b7123703f Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 15 Nov 2023 10:30:21 -0500 Subject: [PATCH 16/30] saved progress on ivs --- src/chemistry_functionality.jl | 26 +++++++++++++++----------- src/expression_utils.jl | 12 +++++------- 2 files changed, 20 insertions(+), 18 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 64d9ee7b56..ca30e79927 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -45,13 +45,8 @@ function make_compound(expr) ((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[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])) + composition = Catalyst.recursive_find_reactants!(expr.args[3], 1, Vector{ReactantStruct}(undef, 0)) components = :([]) # Becomes something like :([C, O]). coefficients = :([]) # Becomes something like :([1, 2]). for comp in composition @@ -59,20 +54,29 @@ function make_compound(expr) push!(coefficients.args, comp.stoichiometry) end + # Extracts the composite species name, as well as the expression which creates it (with e.g. meta data and default values included). + species_expr = expr.args[2] # E.g. :(CO2 = 1.0, [metadata=true]) + iv_expr = :(unique(reduce(vcat, arguments.(ModelingToolkit.unwrap.($(components)))))) # Expression which, when evaluated, becoems a list of all the independent variables the compound should depend on. + species_expr = insert_independent_variable(species_expr, iv_expr) # Add independent variable (e.g. goes from `CO2` to `CO2(t)`). + species_name = find_varname_in_declaration(expr.args[2]) # E.g. :CO2 + + # 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). - non_t_iv_error_check_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)` + # 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). # E.g. `@species CO2(t)` species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)` compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)` 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])` - # Currently, non-t independent variables are not supported for compounds. If there are any like these, we throw an error: - non_t_iv_error_check_expr = Expr(:escape, :(issetequal(unique(reduce(vcat, arguments.(ModelingToolkit.unwrap.($components)))), [t]) || error("Currently, compounds depending on components that are not \"t\" are not supported."))) + println(quote + $species_declaration_expr + $compound_designation_expr + $components_designation_expr + $coefficients_designation_expr + end) # Returns the rephrased expression. return quote - $non_t_iv_error_check_expr $species_declaration_expr $compound_designation_expr $components_designation_expr diff --git a/src/expression_utils.jl b/src/expression_utils.jl index 379fd305b1..e3df736e41 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -52,21 +52,19 @@ end # (In this example the independent variable t was inserted). # Extra dispatch is to ensure so that iv is a vector (so that we can handle both single or multiple iv insertion in one function). # - This way we don't have to make a check for how to create the final expression (which is different for symbol or vector of symbols). -function insert_independent_variable(expr_in, iv) - # If iv is a symbol (e.g. :t), we convert to vector (e.g. [:t]). This way we can handle single and multiple ivs using the same code, without needing to check at the end. - (iv isa Symbol) && (iv = [iv]) +function insert_independent_variable(expr_in, ivs) # 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(:call, expr_in, iv...)) + (expr_in isa Symbol) && (return Expr(:call, expr_in, :($ivs...))) expr = deepcopy(expr_in) if expr.head == :(=) # Case: :(X = 1.0) - expr.args[1] = Expr(:call, expr.args[1], iv...) + expr.args[1] = Expr(:call, expr.args[1], :($ivs...)) elseif expr.head == :tuple if expr.args[1] isa Symbol # Case: :(X, [metadata=true]) - expr.args[1] = Expr(:call, expr.args[1], iv...) + expr.args[1] = Expr(:call, expr.args[1], :($ivs...)) elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol) # Case: :(X = 1.0, [metadata=true]) - expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], iv...) + expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], :($ivs...)) end end (expr == expr_in) && error("Failed to add independent variable $(iv) to expression: $expr_in") From d6d767b61eb36b78414f8fb8d0daf777c5c0ed21 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 15 Nov 2023 14:56:11 -0500 Subject: [PATCH 17/30] fix iv inference --- src/chemistry_functionality.jl | 21 ++++++++--------- src/expression_utils.jl | 15 ++++++------ test/miscellaneous_tests/compound_macro.jl | 27 ++++++++++------------ 3 files changed, 28 insertions(+), 35 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index ca30e79927..738523f4f6 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -56,28 +56,22 @@ function make_compound(expr) # Extracts the composite species name, as well as the expression which creates it (with e.g. meta data and default values included). species_expr = expr.args[2] # E.g. :(CO2 = 1.0, [metadata=true]) - iv_expr = :(unique(reduce(vcat, arguments.(ModelingToolkit.unwrap.($(components)))))) # Expression which, when evaluated, becoems a list of all the independent variables the compound should depend on. - species_expr = insert_independent_variable(species_expr, iv_expr) # Add independent variable (e.g. goes from `CO2` to `CO2(t)`). + species_expr = insert_independent_variable(species_expr, :(..)) # Prepare iv addition, i.e. turns CO to CO(..). species_name = find_varname_in_declaration(expr.args[2]) # E.g. :CO2 + ivs_get_expr = :(unique(reduce(vcat,[arguments(ModelingToolkit.unwrap(iv)) for iv in $components]))) # Expression which when evaluated gives a vector with all the ivs of the components. - # 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). # E.g. `@species CO2(t)` - species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(t)` + # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in the scope the users use the macro in (to e.g. detect already exiting species). + species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(..)` + iv_designation_expression = Expr(:escape, :($species_name = $(species_name)($(ivs_get_expr)...))) # E.g. `CO2 = CO2([...]..)` where [...] is something which evaluates to a vector with all the ivs of the components. compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)` 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(quote - $species_declaration_expr - $compound_designation_expr - $components_designation_expr - $coefficients_designation_expr - end) - # Returns the rephrased expression. return quote $species_declaration_expr + $iv_designation_expression $compound_designation_expr $components_designation_expr $coefficients_designation_expr @@ -128,6 +122,9 @@ function make_compounds(expr) end push!(compound_declarations.args, :($(Expr(:escape, :($(compound_syms)))))) + # The output needs to be converted to Vector{Num} (from Vector{SymbolicUtils.BasicSymbolic{Real}}) to be consistent with e.g. @variables. + compound_declarations.args[end] = :([ModelingToolkit.wrap(cmp) for cmp in $(compound_declarations.args[end])]) + # Returns output that. return compound_declarations end diff --git a/src/expression_utils.jl b/src/expression_utils.jl index e3df736e41..32cbf0f62a 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -49,22 +49,21 @@ end # X(t) = 1.0 # X(t), [metadata=true] # X(t) = 1.0, [metadata=true] -# (In this example the independent variable t was inserted). -# Extra dispatch is to ensure so that iv is a vector (so that we can handle both single or multiple iv insertion in one function). -# - This way we don't have to make a check for how to create the final expression (which is different for symbol or vector of symbols). -function insert_independent_variable(expr_in, ivs) +# (In this example the independent variable :t was inserted). +# Here, the iv is a iv_expr, which can be anything, which is inserted +function insert_independent_variable(expr_in, iv_expr) # 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(:call, expr_in, :($ivs...))) + (expr_in isa Symbol) && (return Expr(:call, expr_in, iv_expr)) expr = deepcopy(expr_in) if expr.head == :(=) # Case: :(X = 1.0) - expr.args[1] = Expr(:call, expr.args[1], :($ivs...)) + expr.args[1] = Expr(:call, expr.args[1], iv_expr) elseif expr.head == :tuple if expr.args[1] isa Symbol # Case: :(X, [metadata=true]) - expr.args[1] = Expr(:call, expr.args[1], :($ivs...)) + expr.args[1] = Expr(:call, expr.args[1], iv_expr) elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol) # Case: :(X = 1.0, [metadata=true]) - expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], :($ivs...)) + expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], iv_expr) end end (expr == expr_in) && error("Failed to add independent variable $(iv) to expression: $expr_in") diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index a0df95f423..892b53d41b 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -71,23 +71,20 @@ end let @variables t x y z @species C(t) H(x) N(x) O(t) P(t,x) S(x,y) - + + # Creates compounds. @compound CO2 ~ C + 2O - @test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t]) + @compound NH4, [output=true] ~ N + 4H + @compound (H2O = 2.0) ~ 2H + O + @compound PH4 ~ P + 4H + @compound SO2 ~ S + 2O - # These 4 tests checks currently non-supported feature. - @test_broken false - @test_broken false - @test_broken false - @test_broken false - # @compound NH4, [output=true] ~ N + 4H - # @compound (H2O = 2.0) ~ 2H + N - # @compound PH4 ~ P + 4H - # @compound SO2 ~ S + 2O - # @test issetequal(arguments(ModelingToolkit.unwrap(NH4)), [x]) - # @test issetequal(arguments(ModelingToolkit.unwrap(H2O)), [t, x]) - # @test issetequal(arguments(ModelingToolkit.unwrap(PH4)), [t, x]) - # @test issetequal(arguments(ModelingToolkit.unwrap(SO2)), [t, x, y]) + # Checks they have the correct independent variables. + @test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t]) + @test issetequal(arguments(ModelingToolkit.unwrap(NH4)), [x]) + @test issetequal(arguments(ModelingToolkit.unwrap(H2O)), [t, x]) + @test issetequal(arguments(ModelingToolkit.unwrap(PH4)), [t, x]) + @test issetequal(arguments(ModelingToolkit.unwrap(SO2)), [t, x, y]) end ### Other Minor Tests ### From 025e07d5316160b837eb77e0873c291395fa9b00 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 15 Nov 2023 14:58:24 -0500 Subject: [PATCH 18/30] iv doc update --- .../chemistry_related_functionality.md | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 00d2f3407a..36591d13a5 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -17,7 +17,7 @@ Next, we create the `CO2` compound species: ```@example chem1 @compound CO2 ~ C + 2O ``` -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: +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 are not designated for compounds (these are 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) ``` @@ -37,16 +37,6 @@ Finally, it is possible to check whether a species is a compound or not using th iscompound(CO2) ``` -!!! note - Currently, compounds with components that depend on variables that are not `t` are not supported. E.g. - ```julia - @variables x y - @species O(x, y) - @compound O2 ~ 2O - ``` - will currently throw an error. - - 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) From f957caa43565bc56388f82df46edc4594b12962c Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 15 Nov 2023 16:19:42 -0500 Subject: [PATCH 19/30] history file update --- HISTORY.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/HISTORY.md b/HISTORY.md index d786a47bdd..02ac9a266e 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -24,6 +24,13 @@ ``` X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.) +- Update how compounds are created. E.g. use +```julia +@variables t C(t) O(t) +@compound CO2 ~ C + 2O +``` +to create a compound species `CO2` that consists of `C` and 2 `O`. +- Added documentation for chemistry related functionality (compound creation and reaction balancing). - Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage: ```julia using Catalyst From 1f4e149525f7624bf825347a1def44bfee568186 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 6 Dec 2023 16:43:38 -0500 Subject: [PATCH 20/30] doc update --- .../chemistry_related_functionality.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 36591d13a5..6186779f9b 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -1,13 +1,13 @@ # [Chemistry-related functionality](@id chemistry_functionality) While Catalyst has primarily been designed around the modelling of biological systems, reaction network models are also common across chemistry. This section describes two types of functionality, that while of general interest, should be especially useful in the modelling of chemical systems. -- The `@compound` option, allowing the user to designate that a specific species is composed of certain subspecies. -- The `balance_reaction` function enabling the user to balance a reaction so the same number of components occur on both sides. +- The `@compound` option, which enables the user to designate that a specific species is composed of certain subspecies. +- The `balance_reaction` function, which enables the user to balance a reaction so the same number of components occur on both sides. ## Modelling with compound species #### Creating compound species programmatically -We will first show how to create compound species through [programmatic construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and stoichiometries). In this example, we will create a CO₂ molecule, consisting of one C atom and two O atoms. First, we create species corresponding to the components: +We will first show how to create compound species through [programmatic model construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and their stoichiometries). In this example, we will create a CO₂ molecule, consisting of one C atom and two O atoms. First, we create species corresponding to the components: ```@example chem1 using Catalyst @variables t @@ -53,7 +53,7 @@ When multiple compounds are created, they can be created simultaneously using th end ``` -#### Creating compound species programmatically +#### Creating compound species within the DSL It is also possible to declare species as compound species within the `@reaction_network` DSL, using the `@compounds` options: ```@example chem1 rn = @reaction_network begin @@ -88,7 +88,7 @@ Default values are designated using `=`, and provided directly after the compoun ```@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: +If both default values and meta data are provided, the metadata is provided after the default value: ```@example chem1 @compound (CO2 = 2.0, [unit="mol"]) ~ C + 2O ``` From f7dca0b7232d248c16204104222db015aa6fe4e2 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 20 Dec 2023 15:48:41 +0100 Subject: [PATCH 21/30] update --- .../chemistry_related_functionality.md | 5 +- src/chemistry_functionality.jl | 67 ++++++++++++------- src/expression_utils.jl | 55 ++++++++------- src/reaction_network.jl | 2 +- test/miscellaneous_tests/compound_macro.jl | 38 +++++------ 5 files changed, 97 insertions(+), 70 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 6186779f9b..002f4812f4 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -82,9 +82,9 @@ as the components `C`, `H`, and `O` are not declared as a species anywhere. Plea #### 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 +@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 `()`: +Default values are designated using `=`, and provided directly after the compound name.: ```@example chem1 @compound (CO2 = 2.0) ~ C + 2O ``` @@ -92,6 +92,7 @@ If both default values and meta data are provided, the metadata is provided afte ```@example chem1 @compound (CO2 = 2.0, [unit="mol"]) ~ C + 2O ``` +In all of these cases, the side to the left of the `~` must be enclosed within `()`. ## 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: diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 738523f4f6..4157f1749d 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -30,48 +30,64 @@ 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." +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 and/or default values, this should be enclosed by \"()\" (e.g. \"@compound (CO2 = 1.0, [output=true]) ~ C + 2O\")." +const COMPOUND_CREATION_ERROR_DEPENDENT_VAR_REQUIRED = "When the components (collectively) have more than 1 independent variable, independent variables have to be specified for the compound (e.g. `@compound CO2(t,s) ~ C + 2O`). This is required since the @compound macro cannot infer the correct order of the independent variables." # Function managing the @compound macro. function make_compound(expr) # Error checks. (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 - # 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])) + # 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]). composition = Catalyst.recursive_find_reactants!(expr.args[3], 1, Vector{ReactantStruct}(undef, 0)) - components = :([]) # Becomes something like :([C, O]). - coefficients = :([]) # Becomes something like :([1, 2]). + 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 - # Extracts the composite species name, as well as the expression which creates it (with e.g. meta data and default values included). - species_expr = expr.args[2] # E.g. :(CO2 = 1.0, [metadata=true]) - species_expr = insert_independent_variable(species_expr, :(..)) # Prepare iv addition, i.e. turns CO to CO(..). - species_name = find_varname_in_declaration(expr.args[2]) # E.g. :CO2 - ivs_get_expr = :(unique(reduce(vcat,[arguments(ModelingToolkit.unwrap(iv)) for iv in $components]))) # Expression which when evaluated gives a vector with all the ivs of the components. + # Extracts: + # - The compound species name (species_name, e.g. `:CO2`). + # - Any ivs attached to it (ivs, e.g. `[]` or `[t,x]`). + # - The expression which creates the compound (species_expr, e.g. `CO2 = 1.0, [metadata=true]`). + species_expr = expr.args[2] + species_name, ivs, _, _ = find_varinfo_in_declaration(expr.args[2]) + + # If no ivs were given, inserts `(..)` (e.g. turning `CO` to `CO(..)`). + isempty(ivs) && (species_expr = insert_independent_variable(species_expr, :(..))) + + # Expression which when evaluated gives a vector with all the ivs of the components. + ivs_get_expr = :(unique(reduce(vcat,[arguments(ModelingToolkit.unwrap(iv)) for iv in $components]))) # Creates the found expressions that will create the compound species. - # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in the scope the users use the macro in (to e.g. detect already exiting species). - species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(..)` - iv_designation_expression = Expr(:escape, :($species_name = $(species_name)($(ivs_get_expr)...))) # E.g. `CO2 = CO2([...]..)` where [...] is something which evaluates to a vector with all the ivs of the components. - compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)` - 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])` + # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in + # the scope the users use the macro in (to e.g. detect already exiting species). + # Creates something like (where `compound_ivs` and `component_ivs` evaluates to all the compound's and components' ivs): + # `@species CO2(..)` + # `isempty([])` && length(component_ivs) && error("When ...) + # `CO2 = CO2(component_ivs..)` + # `issetequal(compound_ivs, component_ivs) || error("The ...)` + # `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)` + # `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])` + # `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])` + species_declaration_expr = Expr(:escape, :(@species $species_expr)) + multiple_ivs_error_check_expr = Expr(:escape, :($(isempty(ivs)) && (length($ivs_get_expr) > 1) && error($COMPOUND_CREATION_ERROR_DEPENDENT_VAR_REQUIRED))) + iv_designation_expr = Expr(:escape, :($(isempty(ivs)) && ($species_name = $(species_name)($(ivs_get_expr)...)))) + iv_check_expr = Expr(:escape, :(issetequal(arguments(ModelingToolkit.unwrap($species_name)), $ivs_get_expr) || error("The independent variable(S) provided to the compound ($(arguments(ModelingToolkit.unwrap($species_name)))), and those of its components ($($ivs_get_expr)))), are not identical."))) + compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) + components_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundComponents, $components))) + coefficients_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundCoefficients, $coefficients))) # Returns the rephrased expression. return quote $species_declaration_expr - $iv_designation_expression + $multiple_ivs_error_check_expr + $iv_designation_expr + $iv_check_expr $compound_designation_expr $components_designation_expr $coefficients_designation_expr @@ -107,7 +123,8 @@ function make_compounds(expr) # 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. + # Creates a compound creation set of lines (4 in total) for each compound line in expr. + # Loops through all 4x[Number of compounds] lines 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) @@ -118,7 +135,7 @@ function make_compounds(expr) push!(compound_declarations.args, :($(Expr(:escape, :([]))))) compound_syms = :([]) for arg in expr.args - push!(compound_syms.args, find_varname_in_declaration(arg.args[2])) + push!(compound_syms.args, find_varinfo_in_declaration(arg.args[2])[1]) end push!(compound_declarations.args, :($(Expr(:escape, :($(compound_syms)))))) diff --git a/src/expression_utils.jl b/src/expression_utils.jl index 32cbf0f62a..48c0cfb5ed 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -1,4 +1,4 @@ -#Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical). +# 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 @@ -19,27 +19,37 @@ end # 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) +# Finds the: Variable name (X), Independent variable name(s) ([t]), default value (2.0), and metadata (:([metadata=true])). +# If a field does not exist (e.g. independent variable in `X, [metadata=true]`), gives nothing. +# The independent variables are given as a vector (empty if none given). +# Does not support e.g. "X [metadata=true]" (when metadata does not have a comma before). +function find_varinfo_in_declaration(expr) + # Case: X + (expr isa Symbol) && (return expr, [], nothing, nothing) + # Case: X(t) + (expr.head == :call) && (return expr.args[1], expr.args[2:end], nothing, nothing) 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 + # Case: X = 1.0 + (expr.args[1] isa Symbol) && (return expr.args[1], [], expr.args[2], nothing) + # Case: X(t) = 1.0 + (expr.args[1].head == :call) && (return expr.args[1].args[1], expr.args[1].args[2:end], expr.args[2].args[1], nothing) 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] + # Case: X, [metadata=true] + (expr.args[1] isa Symbol) && (return expr.args[1], [], nothing, expr.args[2]) + # Case: X(t), [metadata=true] + (expr.args[1].head == :call) && (return expr.args[1].args[1], expr.args[1].args[2:end], nothing, expr.args[2]) 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] + # Case: X = 1.0, [metadata=true] + (expr.args[1].args[1] isa Symbol) && (return expr.args[1].args[1], [], expr.args[1].args[2], expr.args[2]) + # Case: X(t) = 1.0, [metadata=true] + (expr.args[1].args[1].head == :call) && (return expr.args[1].args[1].args[1], expr.args[1].args[1].args[2:end], expr.args[1].args[2].args[1], expr.args[2]) end end error("Unable to detect the variable declared in expression: $expr.") end -# Converts an expression of the form: +# Converts an expression of the forms: # X # X = 1.0 # X, [metadata=true] @@ -53,29 +63,28 @@ end # Here, the iv is a iv_expr, which can be anything, which is inserted function insert_independent_variable(expr_in, iv_expr) # 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). + # 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(:call, expr_in, iv_expr)) expr = deepcopy(expr_in) - if expr.head == :(=) # Case: :(X = 1.0) + # Loops through possible cases. + if expr.head == :(=) + # Case: :(X = 1.0) expr.args[1] = Expr(:call, expr.args[1], iv_expr) elseif expr.head == :tuple - if expr.args[1] isa Symbol # Case: :(X, [metadata=true]) + if expr.args[1] isa Symbol + # Case: :(X, [metadata=true]) expr.args[1] = Expr(:call, expr.args[1], iv_expr) - elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol) # Case: :(X = 1.0, [metadata=true]) + elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol) + # Case: :(X = 1.0, [metadata=true]) expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], iv_expr) 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. diff --git a/src/reaction_network.jl b/src/reaction_network.jl index c214f634b3..68d62796ba 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -472,7 +472,7 @@ function read_compound_options(opts) if haskey(opts, :compounds) compound_expr = opts[:compounds] # 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] + compound_species = [find_varinfo_in_declaration(arg.args[2])[1] 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}[] diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index 892b53d41b..8943316937 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -8,7 +8,7 @@ let # Basic cases that should pass: @compound H2O_1 ~ 2H + O - @compound H2O_2, [output=true] ~ 2H + O + @compound (H2O_2, [output=true]) ~ 2H + O @compound (H2O_3 = 1.5) ~ 2H + O @compound (H2O_4 = 4, [output=true]) ~ 2H + O @compound (H2O_5 = p1, [output=true]) ~ 2H + p2*O @@ -18,26 +18,19 @@ let @test iscompound(H2O_4) @test iscompound(H2O_5) - # Independent variable given. - @test_throws LoadError @eval @compound H2O(t) ~ 2H + O - @test_throws LoadError @eval @compound H2O(t), [output=true] ~ 2H + O - @test_throws LoadError @eval @compound (H2O(t) = 1.5) ~ 2H + O - @test_throws LoadError @eval @compound (H2O(t) = 4, [output=true]) ~ 2H + O - @test_throws LoadError @eval @compound (H2O(t) = p1, [output=true]) ~ 2H + p2*O - # Other errors. - @test_throws LoadError @eval @compound H2O(t) = 2H + O - @test_throws LoadError @eval @compound H2O(t), [output=true] = 2H + O - @test_throws LoadError @eval @compound H2O(t) = 1.5 ~ 2H + O - @test_throws LoadError @eval @compound H2O(t) = 4, [output=true] ~ 2H + O - @test_throws LoadError @eval @compound H2O(t) = p1, [output=true] ~ 2H + p2*O + @test_throws LoadError @eval @compound H2O = 2H + O + @test_throws LoadError @eval @compound (H2O, [output=true]) = 2H + O + @test_throws LoadError @eval @compound H2O = 1.5 ~ 2H + O + @test_throws LoadError @eval @compound (H2O = 4, [output=true]) ~ 2H + O + @test_throws LoadError @eval @compound (H2O = p1, [output=true]) ~ 2H + p2*O # Compounds created in block notation. @compounds begin CO2_1 ~ 2H + O end @compounds begin - CO2_2, [output=true] ~ 2H + O + (CO2_2, [output=true]) ~ 2H + O (CO2_3 = 1.5) ~ 2H + O (CO2_4 = 4, [output=true]) ~ 2H + O (CO2_5 = p1, [output=true]) ~ 2H + p2*O @@ -54,7 +47,7 @@ let @parameters p1 p2 @compounds begin NH3_1 ~ N + 3H - NH3_2, [output=true] ~ N + 3H + (NH3_2, [output=true]) ~ N + 3H (NH3_3 = 1.5) ~ N + 3H (NH3_4 = 4, [output=true]) ~ N + 3H (NH3_5 = p1, [output=true]) ~ N + p2*H @@ -71,13 +64,20 @@ end let @variables t x y z @species C(t) H(x) N(x) O(t) P(t,x) S(x,y) + + # Checks that wrong (or absent) independent variable produces errors. + @test_throws Exception @eval @compound CO2(t,x) ~ C + 2O + @test_throws Exception @eval @compound (NH4(s), [output=true]) ~ N + 4H + @test_throws Exception @eval @compound (H2O = 2.0) ~ 2H + O + @test_throws Exception @eval @compound PH4(x) ~ P + 4H + @test_throws Exception @eval @compound SO2(t,y) ~ S + 2O # Creates compounds. @compound CO2 ~ C + 2O - @compound NH4, [output=true] ~ N + 4H - @compound (H2O = 2.0) ~ 2H + O - @compound PH4 ~ P + 4H - @compound SO2 ~ S + 2O + @compound (NH4, [output=true]) ~ N + 4H + @compound (H2O(t,x) = 2.0) ~ 2H + O + @compound PH4(t,x) ~ P + 4H + @compound SO2(t,x,y) ~ S + 2O # Checks they have the correct independent variables. @test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t]) From e5c2070fda0a1ab431c4a4aff00256a4efbc74cd Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 20 Dec 2023 16:08:34 +0100 Subject: [PATCH 22/30] test fix --- test/miscellaneous_tests/compound_macro.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index 8943316937..a93d6dfcbb 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -19,11 +19,11 @@ let @test iscompound(H2O_5) # Other errors. - @test_throws LoadError @eval @compound H2O = 2H + O - @test_throws LoadError @eval @compound (H2O, [output=true]) = 2H + O - @test_throws LoadError @eval @compound H2O = 1.5 ~ 2H + O - @test_throws LoadError @eval @compound (H2O = 4, [output=true]) ~ 2H + O - @test_throws LoadError @eval @compound (H2O = p1, [output=true]) ~ 2H + p2*O + @test_throws Exception @eval @compound H2O = 2H + O + @test_throws Exception @eval @compound (H2O, [output=true]) = 2H + O + @test_throws Exception @eval @compound (H2O = 1.5) ~ 2H + O + @test_throws Exception @eval @compound (H2O = 4, [output=true]) ~ 2H + O + @test_throws Exception @eval @compound (H2O = p1, [output=true]) ~ 2H + p2*O # Compounds created in block notation. @compounds begin From 850316a96212f0173694010fc79ca3c3fadc87c2 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 20 Dec 2023 16:29:40 +0100 Subject: [PATCH 23/30] improve a code comment. --- src/chemistry_functionality.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 4157f1749d..b01ed8dcf8 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -123,8 +123,8 @@ function make_compounds(expr) # Creates an empty block containing the output call. compound_declarations = Expr(:block) - # Creates a compound creation set of lines (4 in total) for each compound line in expr. - # Loops through all 4x[Number of compounds] lines and add them to compound_declarations. + # For each compound in `expr`, creates the set of 7 compound creation lines (using `make_compound`). + # Next, loops through all 7*[Number of compounds] lines 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) From 71aec3a640c0eca871652dc4a8b8c2ff76074935 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 20 Dec 2023 17:21:51 +0100 Subject: [PATCH 24/30] text improvement --- src/chemistry_functionality.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index b01ed8dcf8..85f969bbfe 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -61,7 +61,7 @@ function make_compound(expr) isempty(ivs) && (species_expr = insert_independent_variable(species_expr, :(..))) # Expression which when evaluated gives a vector with all the ivs of the components. - ivs_get_expr = :(unique(reduce(vcat,[arguments(ModelingToolkit.unwrap(iv)) for iv in $components]))) + ivs_get_expr = :(unique(reduce(vcat,[arguments(ModelingToolkit.unwrap(comp)) for comp in $components]))) # Creates the found expressions that will create the compound species. # The `Expr(:escape, :(...))` is required so that the expressions are evaluated in From 8772169183c6b8e88e8f66363821c1b6a9d8236d Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Sat, 30 Dec 2023 05:23:03 -0500 Subject: [PATCH 25/30] Update docs/src/catalyst_functionality/chemistry_related_functionality.md Co-authored-by: Sam Isaacson --- .../catalyst_functionality/chemistry_related_functionality.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 002f4812f4..5d4e0d7329 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -17,7 +17,7 @@ Next, we create the `CO2` compound species: ```@example chem1 @compound CO2 ~ C + 2O ``` -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 are not designated for compounds (these are 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: +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 are not designated for compounds (these are instead directly inferred from their components). Components with non-unit 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 also a species, and can be used wherever e.g. `C` can be used: ```@example chem1 isspecies(CO2) ``` From a6e17426353d7743aefa32763cfcc96e8a55736c Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Sat, 30 Dec 2023 05:24:04 -0500 Subject: [PATCH 26/30] Update docs/src/catalyst_functionality/chemistry_related_functionality.md Co-authored-by: Sam Isaacson --- .../catalyst_functionality/chemistry_related_functionality.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 5d4e0d7329..5b55f6d4dd 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -95,7 +95,7 @@ If both default values and meta data are provided, the metadata is provided afte In all of these cases, the side to the left of the `~` must be enclosed within `()`. ## 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: +One use of defining a species as a compound is that they can be used to balance reactions so 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 rx = @reaction k, C + O --> $CO2 ``` From e710efa06cb46164aa3770e8cdf852b5e0c57159 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Sat, 30 Dec 2023 05:24:17 -0500 Subject: [PATCH 27/30] Update docs/src/catalyst_functionality/chemistry_related_functionality.md Co-authored-by: Sam Isaacson --- .../catalyst_functionality/chemistry_related_functionality.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 5b55f6d4dd..96bfd0fb11 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -103,7 +103,7 @@ Here, the reaction rate (`k`) is not involved in the reaction balancing. We use ```@example chem1 balance_reaction(rx) ``` -which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). No, or an infinite number of, solutions is also possible. +which correctly finds the (rather trivial) solution `C + 2O --> CO2`. Here we note that `balance_reaction` actually returns a vector. The reason is that the reaction balancing problem may have several solutions. Typically, there is only a single solution (in which case this is the vector's only element). No, or an infinite number of, solutions is also possible depending on the given reaction. Let us consider a more elaborate example, the reaction between ammonia (NH₃) and oxygen (O₂) to form nitrogen monoxide (NO) and water (H₂O). Let us first create the components and the unbalanced reaction: ```@example chem2 From ab54e3c27e215f6f67420b471719ae4d5acc7b95 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Sat, 30 Dec 2023 05:24:46 -0500 Subject: [PATCH 28/30] Update src/chemistry_functionality.jl Co-authored-by: Sam Isaacson --- src/chemistry_functionality.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 85f969bbfe..f51c0b5940 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -18,7 +18,7 @@ Example: ```julia @variables t @species C(t) O(t) -@compound CO2(t) = C + 2O +@compound CO2(t) ~ C + 2O ``` Notes: From cb6c12d3e6668af995d176097000a1c82b95d7f3 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Sat, 30 Dec 2023 05:35:23 -0500 Subject: [PATCH 29/30] Update docs/src/catalyst_functionality/chemistry_related_functionality.md Co-authored-by: Sam Isaacson --- .../catalyst_functionality/chemistry_related_functionality.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 96bfd0fb11..83f6042827 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -123,7 +123,7 @@ We can now create a balanced version (where the amount of H, N, and O is the sam balanced_reaction = balance_reaction(unbalanced_reaction)[1] ``` -Reactions declared as a part of a `ReactionSystem` (e.g. using the DSL) can be retrieved for balancing using the `reaction` function. Please note that balancing these will not mutate the `ReactionSystem`, but a new reaction system will need to be created using the balanced reactions. +Reactions declared as a part of a `ReactionSystem` (e.g. using the DSL) can be retrieved for balancing using the [`reactions`](@ref) function. Please note that balancing these will not mutate the `ReactionSystem`, but a new reaction system will need to be created using the balanced reactions. !!! note Reaction balancing is currently not supported for reactions involving compounds of compounds. \ No newline at end of file From 50a96d76c2224c95479903c976efc276c319d023 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 30 Dec 2023 11:49:12 +0100 Subject: [PATCH 30/30] up --- .../chemistry_related_functionality.md | 23 +++++++++---- src/chemistry_functionality.jl | 34 +++++++++---------- 2 files changed, 33 insertions(+), 24 deletions(-) diff --git a/docs/src/catalyst_functionality/chemistry_related_functionality.md b/docs/src/catalyst_functionality/chemistry_related_functionality.md index 83f6042827..01c44094ef 100644 --- a/docs/src/catalyst_functionality/chemistry_related_functionality.md +++ b/docs/src/catalyst_functionality/chemistry_related_functionality.md @@ -6,7 +6,7 @@ While Catalyst has primarily been designed around the modelling of biological sy ## Modelling with compound species -#### Creating compound species programmatically +### Creating compound species programmatically We will first show how to create compound species through [programmatic model construction](@ref programmatic_CRN_construction), and then demonstrate using the DSL. To create a compound species, use the `@compound` macro, first designating the compound, followed by its components (and their stoichiometries). In this example, we will create a CO₂ molecule, consisting of one C atom and two O atoms. First, we create species corresponding to the components: ```@example chem1 using Catalyst @@ -17,7 +17,7 @@ Next, we create the `CO2` compound species: ```@example chem1 @compound CO2 ~ C + 2O ``` -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 are not designated for compounds (these are instead directly inferred from their components). Components with non-unit 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 also a species, 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 are generally not designated for compounds (these are instead directly inferred from their components). Components with non-unit 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 also a species, and can be used wherever e.g. `C` can be used: ```@example chem1 isspecies(CO2) ``` @@ -53,7 +53,7 @@ When multiple compounds are created, they can be created simultaneously using th end ``` -#### Creating compound species within the DSL +### Creating compound species within the DSL It is also possible to declare species as compound species within the `@reaction_network` DSL, using the `@compounds` options: ```@example chem1 rn = @reaction_network begin @@ -63,12 +63,12 @@ rn = @reaction_network begin H2O ~ 2H + O H2CO3 ~ CO2 + H2O end - (k1,k2), H2O+ CO2 <--> H2CO3 + (k1,k2), H2O + CO2 <--> H2CO3 end ``` -When creating compound species using the DSL, it is important to note that *every component must be known to the system as a species, either by being declared using the `@species` option, or by appearing in a reaction*. E.g. the following is not valid +When creating compound species using the DSL, it is important to note that *every component must be known to the system as a species, either by being declared using the `@species` or `@compound` options, or by appearing in a reaction*. E.g. the following is not valid ```julia -rn = @reaction_network begin +rn = @reaction_network begin`` @compounds begin C2O ~ C + 2O H2O ~ 2H + O @@ -79,7 +79,7 @@ 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 +### 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 @@ -94,6 +94,15 @@ If both default values and meta data are provided, the metadata is provided afte ``` In all of these cases, the side to the left of the `~` must be enclosed within `()`. +### Compounds with multiple independent variables +While we generally do not need to specify independent variables for compound, if the components (together) have more than one independent variable, this have to be done: +```@example chem1 +@variables t s +@species N(s) O(t) +@compound NO2(t,s) ~ N + 2O +``` +Here, `NO2` depend both on a spatial independent variable (`s`) and a time one (`t`). This is required since, while multiple independent variables can be inferred, their internal order cannot (and must hence be provided by the user). + ## Balancing chemical reactions One use of defining a species as a compound is that they can be used to balance reactions so 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 diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index f51c0b5940..a839385ab3 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -284,35 +284,35 @@ Example: ```julia @variables t @species Si(t) Cl(t) H(t) O(t) -@compound SiCl4(t) = Si + 4Cl -@compound H2O(t) = 2H + O -@compound H4SiO4(t) = 4H + Si + 4O -@compound HCl(t) = H + Cl -rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl]) +@compound SiCl4 ~ Si + 4Cl +@compound H2O ~ 2H + O +@compound H4SiO4 ~ 4H + Si + 4O +@compound HCl ~ H + Cl +rx = @reaction 1.0, SiCl4 + H2O --> H4SiO4 HCl balance_reaction(rx) # Exactly one solution. ``` ```julia @variables t @species C(t) H(t) O(t) -@compound CO(t) = C + O -@compound CO2(t) = C + 2O -@compound H2(t) = 2H -@compound CH4(t) = C + 4H -@compound H2O(t) = 2H + O -rx = Reaction(1.0, [CO, CO2, H2], [CH4, H2O]) +@compound CO ~ C + O +@compound CO2 ~ C + 2O +@compound H2 ~ 2H +@compound CH4 ~ C + 4H +@compound H2O ~ 2H + O +rx = @reaction 1.0, CO + CO2 + H2--> CH4 H2O balance_reaction(rx) # Multiple solutions. ``` ```julia @variables t @species Fe(t) S(t) O(t) H(t) N(t) -@compound FeS2(t) = Fe + 2S -@compound HNO3(t) = H + N + 3O -@compound Fe2S3O12(t) = 2Fe + 3S + 12O -@compound NO(t) = N + O -@compound H2SO4(t) = 2H + S + 4O -rx = Reaction(1.0, [FeS2, HNO3], [Fe2S3O12, NO, H2SO4]) +@compound FeS2 ~ Fe + 2S +@compound HNO3 ~ H + N + 3O +@compound Fe2S3O12 ~ 2Fe + 3S + 12O +@compound NO ~ N + O +@compound H2SO4 ~ 2H + S + 4O +rx = @reaction 1.0, FeS2 + HNO3 --> Fe2S3O12 NO + H2SO4 brxs = balance_reaction(rx) # No solution. ```