Skip to content

Commit f7dca0b

Browse files
committed
update
1 parent 1f4e149 commit f7dca0b

File tree

5 files changed

+97
-70
lines changed

5 files changed

+97
-70
lines changed

docs/src/catalyst_functionality/chemistry_related_functionality.md

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -82,16 +82,17 @@ as the components `C`, `H`, and `O` are not declared as a species anywhere. Plea
8282
#### Designating metadata and default values for compounds
8383
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 `,`:
8484
```@example chem1
85-
@compound CO2, [unit="mol"] ~ C + 2O
85+
@compound (CO2, [unit="mol"]) ~ C + 2O
8686
```
87-
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 `()`:
87+
Default values are designated using `=`, and provided directly after the compound name.:
8888
```@example chem1
8989
@compound (CO2 = 2.0) ~ C + 2O
9090
```
9191
If both default values and meta data are provided, the metadata is provided after the default value:
9292
```@example chem1
9393
@compound (CO2 = 2.0, [unit="mol"]) ~ C + 2O
9494
```
95+
In all of these cases, the side to the left of the `~` must be enclosed within `()`.
9596

9697
## Balancing chemical reactions
9798
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:

src/chemistry_functionality.jl

Lines changed: 42 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -30,48 +30,64 @@ end
3030

3131
# Declares compound error messages:
3232
const COMPOUND_CREATION_ERROR_BASE = "Malformed input to @compound. Should use form like e.g. \"@compound CO2 ~ C + 2O\"."
33-
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\")."
34-
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."
33+
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\")."
34+
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."
3535

3636
# Function managing the @compound macro.
3737
function make_compound(expr)
3838
# Error checks.
3939
(expr isa Expr) || error(COMPOUND_CREATION_ERROR_BASE)
4040
((expr.head == :call) && (expr.args[1] == :~) && (length(expr.args) == 3)) || error(COMPOUND_CREATION_ERROR_BAD_SEPARATOR)
41-
if !(expr.args[2] isa Symbol)
42-
# Checks if a dependent variable was given (e.g. CO2(t)).
43-
(expr.args[2].head == :call) && error(COMPOUND_CREATION_ERROR_DEPENDENT_VAR_GIVEN)
44-
((expr.args[2].args[1] isa Expr) && expr.args[2].args[1].head == :call) && error(COMPOUND_CREATION_ERROR_DEPENDENT_VAR_GIVEN)
45-
((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)
46-
end
4741

48-
# 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]))
42+
# Loops through all components, add the component and the coefficients to the corresponding vectors
43+
# Cannot extract directly using e.g. "getfield.(composition, :reactant)" because then
44+
# we get something like :([:C, :O]), rather than :([C, O]).
4945
composition = Catalyst.recursive_find_reactants!(expr.args[3], 1, Vector{ReactantStruct}(undef, 0))
50-
components = :([]) # Becomes something like :([C, O]).
51-
coefficients = :([]) # Becomes something like :([1, 2]).
46+
components = :([]) # Becomes something like :([C, O]).
47+
coefficients = :([]) # Becomes something like :([1, 2]).
5248
for comp in composition
5349
push!(components.args, comp.reactant)
5450
push!(coefficients.args, comp.stoichiometry)
5551
end
5652

57-
# Extracts the composite species name, as well as the expression which creates it (with e.g. meta data and default values included).
58-
species_expr = expr.args[2] # E.g. :(CO2 = 1.0, [metadata=true])
59-
species_expr = insert_independent_variable(species_expr, :(..)) # Prepare iv addition, i.e. turns CO to CO(..).
60-
species_name = find_varname_in_declaration(expr.args[2]) # E.g. :CO2
61-
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.
53+
# Extracts:
54+
# - The compound species name (species_name, e.g. `:CO2`).
55+
# - Any ivs attached to it (ivs, e.g. `[]` or `[t,x]`).
56+
# - The expression which creates the compound (species_expr, e.g. `CO2 = 1.0, [metadata=true]`).
57+
species_expr = expr.args[2]
58+
species_name, ivs, _, _ = find_varinfo_in_declaration(expr.args[2])
59+
60+
# If no ivs were given, inserts `(..)` (e.g. turning `CO` to `CO(..)`).
61+
isempty(ivs) && (species_expr = insert_independent_variable(species_expr, :(..)))
62+
63+
# Expression which when evaluated gives a vector with all the ivs of the components.
64+
ivs_get_expr = :(unique(reduce(vcat,[arguments(ModelingToolkit.unwrap(iv)) for iv in $components])))
6265

6366
# Creates the found expressions that will create the compound species.
64-
# 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).
65-
species_declaration_expr = Expr(:escape, :(@species $species_expr)) # E.g. `@species CO2(..)`
66-
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.
67-
compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)`
68-
components_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundComponents, $components))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])`
69-
coefficients_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundCoefficients, $coefficients))) # E.g. `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])`
67+
# The `Expr(:escape, :(...))` is required so that the expressions are evaluated in
68+
# the scope the users use the macro in (to e.g. detect already exiting species).
69+
# Creates something like (where `compound_ivs` and `component_ivs` evaluates to all the compound's and components' ivs):
70+
# `@species CO2(..)`
71+
# `isempty([])` && length(component_ivs) && error("When ...)
72+
# `CO2 = CO2(component_ivs..)`
73+
# `issetequal(compound_ivs, component_ivs) || error("The ...)`
74+
# `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, true)`
75+
# `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [C, O])`
76+
# `CO2 = ModelingToolkit.setmetadata(CO2, Catalyst.CompoundSpecies, [1, 2])`
77+
species_declaration_expr = Expr(:escape, :(@species $species_expr))
78+
multiple_ivs_error_check_expr = Expr(:escape, :($(isempty(ivs)) && (length($ivs_get_expr) > 1) && error($COMPOUND_CREATION_ERROR_DEPENDENT_VAR_REQUIRED)))
79+
iv_designation_expr = Expr(:escape, :($(isempty(ivs)) && ($species_name = $(species_name)($(ivs_get_expr)...))))
80+
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.")))
81+
compound_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundSpecies, true)))
82+
components_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundComponents, $components)))
83+
coefficients_designation_expr = Expr(:escape, :($species_name = ModelingToolkit.setmetadata($species_name, Catalyst.CompoundCoefficients, $coefficients)))
7084

7185
# Returns the rephrased expression.
7286
return quote
7387
$species_declaration_expr
74-
$iv_designation_expression
88+
$multiple_ivs_error_check_expr
89+
$iv_designation_expr
90+
$iv_check_expr
7591
$compound_designation_expr
7692
$components_designation_expr
7793
$coefficients_designation_expr
@@ -107,7 +123,8 @@ function make_compounds(expr)
107123
# Creates an empty block containing the output call.
108124
compound_declarations = Expr(:block)
109125

110-
# 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.
126+
# Creates a compound creation set of lines (4 in total) for each compound line in expr.
127+
# Loops through all 4x[Number of compounds] lines and add them to compound_declarations.
111128
compound_calls = [Catalyst.make_compound(line) for line in expr.args]
112129
for compound_call in compound_calls, line in MacroTools.striplines(compound_call).args
113130
push!(compound_declarations.args, line)
@@ -118,7 +135,7 @@ function make_compounds(expr)
118135
push!(compound_declarations.args, :($(Expr(:escape, :([])))))
119136
compound_syms = :([])
120137
for arg in expr.args
121-
push!(compound_syms.args, find_varname_in_declaration(arg.args[2]))
138+
push!(compound_syms.args, find_varinfo_in_declaration(arg.args[2])[1])
122139
end
123140
push!(compound_declarations.args, :($(Expr(:escape, :($(compound_syms))))))
124141

src/expression_utils.jl

Lines changed: 32 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
#Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical).
1+
# Returns the length of a expression tuple, or 1 if it is not an expression tuple (probably a Symbol/Numerical).
22
function tup_leng(ex::ExprValues)
33
(typeof(ex) == Expr && ex.head == :tuple) && (return length(ex.args))
44
return 1
@@ -19,27 +19,37 @@ end
1919
# X(t) = 1.0
2020
# X(t), [metadata=true]
2121
# X(t) = 1.0, [metadata=true]
22-
# Finds the variable name (here, X).
23-
# Currently does not support e.g. "X, [metadata=true]" (when metadata does not have a comma before).
24-
function find_varname_in_declaration(expr)
25-
(expr isa Symbol) && (return expr) # Case: X
26-
(expr.head == :call) && (return ex5.args[1]) # Case: X(t)
22+
# Finds the: Variable name (X), Independent variable name(s) ([t]), default value (2.0), and metadata (:([metadata=true])).
23+
# If a field does not exist (e.g. independent variable in `X, [metadata=true]`), gives nothing.
24+
# The independent variables are given as a vector (empty if none given).
25+
# Does not support e.g. "X [metadata=true]" (when metadata does not have a comma before).
26+
function find_varinfo_in_declaration(expr)
27+
# Case: X
28+
(expr isa Symbol) && (return expr, [], nothing, nothing)
29+
# Case: X(t)
30+
(expr.head == :call) && (return expr.args[1], expr.args[2:end], nothing, nothing)
2731
if expr.head == :(=)
28-
(expr.args[1] isa Symbol) && (return expr.args[1]) # Case: X = 1.0
29-
(expr.args[1].head == :call) && (return expr.args[1].args[1]) # Case: X(t) = 1.0
32+
# Case: X = 1.0
33+
(expr.args[1] isa Symbol) && (return expr.args[1], [], expr.args[2], nothing)
34+
# Case: X(t) = 1.0
35+
(expr.args[1].head == :call) && (return expr.args[1].args[1], expr.args[1].args[2:end], expr.args[2].args[1], nothing)
3036
end
3137
if expr.head == :tuple
32-
(expr.args[1] isa Symbol) && (return expr.args[1]) # Case: X, [metadata=true]
33-
(expr.args[1].head == :call) && (return expr.args[1].args[1]) # Case: X(t), [metadata=true]
38+
# Case: X, [metadata=true]
39+
(expr.args[1] isa Symbol) && (return expr.args[1], [], nothing, expr.args[2])
40+
# Case: X(t), [metadata=true]
41+
(expr.args[1].head == :call) && (return expr.args[1].args[1], expr.args[1].args[2:end], nothing, expr.args[2])
3442
if (expr.args[1].head == :(=))
35-
(expr.args[1].args[1] isa Symbol) && (return expr.args[1].args[1]) # Case: X = 1.0, [metadata=true]
36-
(expr.args[1].args[1].head == :call) && (return expr.args[1].args[1].args[1]) # Case: X(t) = 1.0, [metadata=true]
43+
# Case: X = 1.0, [metadata=true]
44+
(expr.args[1].args[1] isa Symbol) && (return expr.args[1].args[1], [], expr.args[1].args[2], expr.args[2])
45+
# Case: X(t) = 1.0, [metadata=true]
46+
(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])
3747
end
3848
end
3949
error("Unable to detect the variable declared in expression: $expr.")
4050
end
4151

42-
# Converts an expression of the form:
52+
# Converts an expression of the forms:
4353
# X
4454
# X = 1.0
4555
# X, [metadata=true]
@@ -53,29 +63,28 @@ end
5363
# Here, the iv is a iv_expr, which can be anything, which is inserted
5464
function insert_independent_variable(expr_in, iv_expr)
5565
# If expr is a symbol, just attach the iv. If not we have to create a new expr and mutate it.
56-
# Because Symbols (a possible input) cannot be mutated, this function cannot mutate the input (would have been easier if Expr input was guaranteed).
66+
# Because Symbols (a possible input) cannot be mutated, this function cannot mutate the input
67+
# (would have been easier if Expr input was guaranteed).
5768
(expr_in isa Symbol) && (return Expr(:call, expr_in, iv_expr))
5869
expr = deepcopy(expr_in)
5970

60-
if expr.head == :(=) # Case: :(X = 1.0)
71+
# Loops through possible cases.
72+
if expr.head == :(=)
73+
# Case: :(X = 1.0)
6174
expr.args[1] = Expr(:call, expr.args[1], iv_expr)
6275
elseif expr.head == :tuple
63-
if expr.args[1] isa Symbol # Case: :(X, [metadata=true])
76+
if expr.args[1] isa Symbol
77+
# Case: :(X, [metadata=true])
6478
expr.args[1] = Expr(:call, expr.args[1], iv_expr)
65-
elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol) # Case: :(X = 1.0, [metadata=true])
79+
elseif (expr.args[1].head == :(=)) && (expr.args[1].args[1] isa Symbol)
80+
# Case: :(X = 1.0, [metadata=true])
6681
expr.args[1].args[1] = Expr(:call, expr.args[1].args[1], iv_expr)
6782
end
6883
end
69-
(expr == expr_in) && error("Failed to add independent variable $(iv) to expression: $expr_in")
7084
return expr
7185
end
7286

7387

74-
75-
76-
77-
78-
7988
### Old Stuff ###
8089

8190
#This will be called whenever a function stored in funcdict is called.

src/reaction_network.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -472,7 +472,7 @@ function read_compound_options(opts)
472472
if haskey(opts, :compounds)
473473
compound_expr = opts[:compounds]
474474
# Find compound species names, and append the independent variable.
475-
compound_species = [find_varname_in_declaration(arg.args[2]) for arg in compound_expr.args[3].args]
475+
compound_species = [find_varinfo_in_declaration(arg.args[2])[1] for arg in compound_expr.args[3].args]
476476
else # If option is not used, return empty vectors and expressions.
477477
compound_expr = :()
478478
compound_species = Union{Symbol, Expr}[]

test/miscellaneous_tests/compound_macro.jl

Lines changed: 19 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ let
88

99
# Basic cases that should pass:
1010
@compound H2O_1 ~ 2H + O
11-
@compound H2O_2, [output=true] ~ 2H + O
11+
@compound (H2O_2, [output=true]) ~ 2H + O
1212
@compound (H2O_3 = 1.5) ~ 2H + O
1313
@compound (H2O_4 = 4, [output=true]) ~ 2H + O
1414
@compound (H2O_5 = p1, [output=true]) ~ 2H + p2*O
@@ -18,26 +18,19 @@ let
1818
@test iscompound(H2O_4)
1919
@test iscompound(H2O_5)
2020

21-
# Independent variable given.
22-
@test_throws LoadError @eval @compound H2O(t) ~ 2H + O
23-
@test_throws LoadError @eval @compound H2O(t), [output=true] ~ 2H + O
24-
@test_throws LoadError @eval @compound (H2O(t) = 1.5) ~ 2H + O
25-
@test_throws LoadError @eval @compound (H2O(t) = 4, [output=true]) ~ 2H + O
26-
@test_throws LoadError @eval @compound (H2O(t) = p1, [output=true]) ~ 2H + p2*O
27-
2821
# Other errors.
29-
@test_throws LoadError @eval @compound H2O(t) = 2H + O
30-
@test_throws LoadError @eval @compound H2O(t), [output=true] = 2H + O
31-
@test_throws LoadError @eval @compound H2O(t) = 1.5 ~ 2H + O
32-
@test_throws LoadError @eval @compound H2O(t) = 4, [output=true] ~ 2H + O
33-
@test_throws LoadError @eval @compound H2O(t) = p1, [output=true] ~ 2H + p2*O
22+
@test_throws LoadError @eval @compound H2O = 2H + O
23+
@test_throws LoadError @eval @compound (H2O, [output=true]) = 2H + O
24+
@test_throws LoadError @eval @compound H2O = 1.5 ~ 2H + O
25+
@test_throws LoadError @eval @compound (H2O = 4, [output=true]) ~ 2H + O
26+
@test_throws LoadError @eval @compound (H2O = p1, [output=true]) ~ 2H + p2*O
3427

3528
# Compounds created in block notation.
3629
@compounds begin
3730
CO2_1 ~ 2H + O
3831
end
3932
@compounds begin
40-
CO2_2, [output=true] ~ 2H + O
33+
(CO2_2, [output=true]) ~ 2H + O
4134
(CO2_3 = 1.5) ~ 2H + O
4235
(CO2_4 = 4, [output=true]) ~ 2H + O
4336
(CO2_5 = p1, [output=true]) ~ 2H + p2*O
@@ -54,7 +47,7 @@ let
5447
@parameters p1 p2
5548
@compounds begin
5649
NH3_1 ~ N + 3H
57-
NH3_2, [output=true] ~ N + 3H
50+
(NH3_2, [output=true]) ~ N + 3H
5851
(NH3_3 = 1.5) ~ N + 3H
5952
(NH3_4 = 4, [output=true]) ~ N + 3H
6053
(NH3_5 = p1, [output=true]) ~ N + p2*H
@@ -71,13 +64,20 @@ end
7164
let
7265
@variables t x y z
7366
@species C(t) H(x) N(x) O(t) P(t,x) S(x,y)
67+
68+
# Checks that wrong (or absent) independent variable produces errors.
69+
@test_throws Exception @eval @compound CO2(t,x) ~ C + 2O
70+
@test_throws Exception @eval @compound (NH4(s), [output=true]) ~ N + 4H
71+
@test_throws Exception @eval @compound (H2O = 2.0) ~ 2H + O
72+
@test_throws Exception @eval @compound PH4(x) ~ P + 4H
73+
@test_throws Exception @eval @compound SO2(t,y) ~ S + 2O
7474

7575
# Creates compounds.
7676
@compound CO2 ~ C + 2O
77-
@compound NH4, [output=true] ~ N + 4H
78-
@compound (H2O = 2.0) ~ 2H + O
79-
@compound PH4 ~ P + 4H
80-
@compound SO2 ~ S + 2O
77+
@compound (NH4, [output=true]) ~ N + 4H
78+
@compound (H2O(t,x) = 2.0) ~ 2H + O
79+
@compound PH4(t,x) ~ P + 4H
80+
@compound SO2(t,x,y) ~ S + 2O
8181

8282
# Checks they have the correct independent variables.
8383
@test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t])

0 commit comments

Comments
 (0)