From a3aafe3e9eb82565e464d2896108a763bc14a083 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Wed, 2 Aug 2023 20:31:17 +0530 Subject: [PATCH 01/37] Adds balancing code and tests --- src/compound.jl | 173 +++++++++++++++++++++ test/miscellaneous_tests/compound_macro.jl | 36 +++++ 2 files changed, 209 insertions(+) diff --git a/src/compound.jl b/src/compound.jl index 9f2ac5bca2..3e540393e9 100644 --- a/src/compound.jl +++ b/src/compound.jl @@ -87,3 +87,176 @@ 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 + + +### Balancing Code + +function create_matrix(reaction::Catalyst.Reaction) + compounds = [reaction.substrates; reaction.products] + atoms = [] + n_atoms = 0 + A = zeros(Int, 0, length(compounds)) + + for (j, compound) in enumerate(compounds) + if iscompound(compound) + pairs = component_coefficients(compound) + if pairs == Nothing + continue + end + else + pairs = [(compound, 1)] + end + + for pair in pairs + atom, coeff = pair + i = findfirst(x -> isequal(x, atom), atoms) + if i === nothing + push!(atoms, atom) + n_atoms += 1 + A = [A; zeros(Int, 1, length(compounds))] + i = n_atoms + end + coeff = any(map(p -> isequal(p, compounds[j]), reaction.products)) ? -coeff : coeff + A[i, j] = coeff + end + end + # Append a row with last element as 1 + new_row = zeros(Int, 1, size(A, 2)) + new_row[end] = 1 + A = vcat(A, new_row) + return A +end + +function get_stoich(reaction::Reaction) + # Create the matrix A using create_matrix function. + A = create_matrix(reaction) + + # Create the vector b. The last element is 1 and others are 0. + B = zeros(Int64, size(A,1)) + B[end] = 1 + + # Concatenate B to A to form an augmented matrix + AB = [A B] + + # Apply the Bareiss algorithm + ModelingToolkit.bareiss!(AB) + + # Extract the transformed A and B + A_transformed = AB[:, 1:end-1] + B_transformed = AB[:, end] + + # Convert A_transformed to rational numbers + A_transformed_rational = Rational.(A_transformed) + + # Convert B_transformed to rational numbers + B_transformed_rational = Rational.(B_transformed) + + # Perform the division + X = A_transformed_rational \ B_transformed_rational + + # Get the denominators of the rational numbers in X + denominators = denominator.(X) + + # Compute the LCM of the denominators + lcm_value = reduce(lcm, denominators) + + # Multiply each element in X by the LCM of the denominators + X_multiplied = X .* lcm_value + + # Convert the rational numbers to integers + X_integers = numerator.(X_multiplied) + + return X_integers +end + +function balance(reaction::Reaction) + # Calculate the stoichiometric coefficients for the balanced reaction. + stoichiometries = get_stoich(reaction) + + # Divide the stoichiometry vector into substrate and product stoichiometries. + substoich = stoichiometries[1:length(reaction.substrates)] + prodstoich = stoichiometries[(length(reaction.substrates)) + 1:end] + + # Create a new reaction with the balanced stoichiometries + balanced_reaction = Reaction(reaction.rate, reaction.substrates, reaction.products, substoich, prodstoich) + + # Return the balanced reaction + return balanced_reaction +end + +# function get_stoich(reaction::Reaction) ## Gaussian Elimination +# # Create the matrix A using create_matrix function. +# A = create_matrix(reaction) + +# # Create the vector b. The last element is 1 and others are 0. +# B = zeros(Int64, size(A,1)) +# B[end] = 1 + +# # Convert A and B to sparse matrices +# A = sparse(A) + +# # Apply the LU factorization (Gaussian elimination) +# lu_f = lu(A) + +# # Convert B to a dense vector +# B_dense = Vector(B) + +# # Solve for X using the LU factorization +# X = lu_f \ B_dense + +# # Get the smallest positive value in X +# smallest_value = minimum(X[X .> 0]) + +# # Normalize X to be integers +# X_normalized = round.(Int64, X / smallest_value) + +# return X_normalized +# end + +# function get_stoich(reaction::Reaction) ##Bareiss Algorithm (Does not avoid floating points) +# # Create the matrix A using create_matrix function. +# A = create_matrix(reaction) + +# # Create the vector b. The last element is 1 and others are 0. +# B = zeros(Int64, size(A,1)) +# B[end] = 1 + +# # Concatenate B to A to form an augmented matrix +# AB = [A B] + +# # Apply the Bareiss algorithm +# ModelingToolkit.bareiss!(AB) + +# # Extract the transformed A and B +# A_transformed = AB[:, 1:end-1] +# B_transformed = AB[:, end] + +# # Solve for X using back substitution +# X = A_transformed \ B_transformed + +# # Get the smallest positive value in X +# smallest_value = minimum(X[X .> 0]) + +# # Normalize X to be integers +# X_normalized = round.(Int64, X / smallest_value) + +# return X_normalized +# end + +# @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 +# rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) + +# using LinearAlgebra +# using SparseArrays +# using SuiteSparse.UMFPACK + +# gcd_value = gcd(X...) +# X_scaled = X ./ gcd_value + +# X = [1,1,0.2,1] \ No newline at end of file diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index 24d72b4019..c37f4ea642 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -96,3 +96,39 @@ end # @test isequal(coefficients(C10H12)[1], 2) # end +#Check that balancing works. +let + @variables t + @parameters k + @species H(t) O(t) + @compound H2(t) 2H + @compound O2(t) 2O + @compound H2O(t) 2H 1O + + rx = Reaction(k,[H2,O2],[H2O]) + + @test isequal(create_matrix(rx),[2 0 -2; 0 2 -1; 0 0 1;]) + @test isequal(get_stoich(rx),[2, 1, 2]) + + balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) + @test isequal(balanced_rx, balance(rx)) + +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 + + rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) + + @test isequal(create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0; 0 0 0 1;]) + @test isequal(get_stoich(rx),[6, 6, 1, 6]) + + balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1,6]) + @test isequal(balanced_rx, balance(rx)) +end \ No newline at end of file From ca53de44fbe879b3e9f7994ad82cfb8cd46b548b Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Fri, 4 Aug 2023 13:10:45 +0530 Subject: [PATCH 02/37] Update to master --- src/Catalyst.jl | 4 +- test/miscellaneous_tests/compound_macro.jl | 114 ++++++++++++++------- 2 files changed, 81 insertions(+), 37 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 68570b0a4b..995b0699f5 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -95,7 +95,7 @@ export Graph, savegraph, complexgraph # for creating compounds include("compound.jl") export @compound -export components, iscompound, coefficients, component_coefficients, balance, get_stoich -export create_matrix +export components, iscompound, coefficients +export get_stoich, create_matrix, balance end # module diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index c37f4ea642..c9a860b2a4 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -14,11 +14,8 @@ let @test isequal([C, H, O], components(C6H12O2)) @test isequal([6, 12, 2], coefficients(C6H12O2)) - @test isequal([C => 6, H => 12, O => 2], component_coefficients(C6H12O2)) # Maybe chaneg "component_coefficients" name soem time. + @test isequal([C => 6, H => 12, O => 2], Catalyst.component_coefficients(C6H12O2)) @test all(!iscompound(i) for i in components(C6H12O2)) - # Test threw exception - # Expression: all(.!(iscompound.components(C6H12O2))) - # type #iscompound has no field components end let @@ -32,7 +29,7 @@ let @test isequal([O], components(O2)) @test isequal([2], coefficients(O2)) - @test isequal([O => 2], component_coefficients(O2)) # Maybe chaneg "component_coefficients" name soem time. + @test isequal([O => 2], Catalyst.component_coefficients(O2)) @test all(!iscompound(i) for i in components(O2)) end @@ -40,11 +37,11 @@ end let @variables t @species C(t) H(t) - @test_throws Exception @compound C6H12O2(t) 6C 12H 2O # @test_throws might behave weirdly here, if you think this fails as it should but the test comes out starneg, just report and we will figure it out. + @test_throws Exception @compound C6H12O2(t) 6C 12H 2O end let @variables t - @test_throws Exception @compound O2(t) 2O # @test_throws might behave weirdly here, if you think this fails as it should but the test comes out starneg, just report and we will figure it out. + @test_throws Exception @compound O2(t) 2O end # Checks that nested components works as expected. @@ -69,32 +66,79 @@ let @test isequal([3, 5, 3], coefficients(C3H5OH3)) end -# # Checks that interpolation works. -# 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 - -# @test isequal(components(C6H12O2_1), components(C6H12O2_2)) -# @test isequal(coefficients(C6H12O2_1), coefficients(C6H12O2_2)) -# end - -# let -# @variables t -# @species C(t) H(t) -# @compound Cyclopentadiene(t) 5C 6H -# C5H6 = Cyclopentadiene -# @compound C10H12(t) 2C5H6 - -# @test iscompound(C10H12) -# @test iscompound(components(C10H12)[1]) - -# @test isequal(components(C10H12)[1], C5H6) -# @test isequal(components(C10H12)[1], Cyclopentadiene) -# @test isequal(coefficients(C10H12)[1], 2) -# end +# Checks that interpolation works. +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 + + @test iscompound(C6H12O2_1) + @test iscompound(C6H12O2_2) + + @test isequal(components(C6H12O2_1), components(C6H12O2_2)) + @test isequal(coefficients(C6H12O2_1), coefficients(C6H12O2_2)) +end + +let + @variables t + @species C(t) H(t) + @compound Cyclopentadiene(t) 5C 6H + C5H6 = Cyclopentadiene + @compound C10H12(t) 2C5H6 + + @test iscompound(C10H12) + @test iscompound(components(C10H12)[1]) + + @test isequal(components(C10H12)[1], C5H6) + @test isequal(components(C10H12)[1], Cyclopentadiene) + @test isequal(coefficients(C10H12)[1], 2) +end + +let + @variables t + @species H(t) + + alpha = 2 + @compound H2_1(t) alpha*H + @compound H2_2(t) 2H + + @test iscompound(H2_1) + @test iscompound(H2_2) + + @test isequal(components(H2_1),components(H2_2)) + @test isequal(coefficients(H2_1),coefficients(H2_2)) +end + +let + @variables t + @parameters alpha = 2 + @species H(t) + + @compound H2_1(t) alpha*H + @compound H2_2(t) 2H + + @test iscompound(H2_1) + @test iscompound(H2_2) + + @test isequal(components(H2_1),components(H2_2)) + @test isequal(coefficients(H2_1), @parameters alpha = 2) +end + +let + @variables t + @species A(t) + B = A + @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)) +end #Check that balancing works. let @@ -109,7 +153,7 @@ let @test isequal(create_matrix(rx),[2 0 -2; 0 2 -1; 0 0 1;]) @test isequal(get_stoich(rx),[2, 1, 2]) - + balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) @test isequal(balanced_rx, balance(rx)) @@ -128,7 +172,7 @@ let @test isequal(create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0; 0 0 0 1;]) @test isequal(get_stoich(rx),[6, 6, 1, 6]) - + balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1,6]) @test isequal(balanced_rx, balance(rx)) end \ No newline at end of file From e9ad3e7da78d07e144feb47e231b5f17621236ad Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Fri, 4 Aug 2023 13:31:23 +0530 Subject: [PATCH 03/37] Update tests --- test/miscellaneous_tests/compound_macro.jl | 85 +++++++++++++++++++++- 1 file changed, 84 insertions(+), 1 deletion(-) diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index a8ecebe28b..c9a860b2a4 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -14,7 +14,7 @@ let @test isequal([C, H, O], components(C6H12O2)) @test isequal([6, 12, 2], coefficients(C6H12O2)) - @test isequal([C => 6, H => 12, O => 2], Catalyst.component_coefficients(C6H12O2)) + @test isequal([C => 6, H => 12, O => 2], Catalyst.component_coefficients(C6H12O2)) @test all(!iscompound(i) for i in components(C6H12O2)) end @@ -74,6 +74,9 @@ let @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)) end @@ -93,3 +96,83 @@ let @test isequal(coefficients(C10H12)[1], 2) end +let + @variables t + @species H(t) + + alpha = 2 + @compound H2_1(t) alpha*H + @compound H2_2(t) 2H + + @test iscompound(H2_1) + @test iscompound(H2_2) + + @test isequal(components(H2_1),components(H2_2)) + @test isequal(coefficients(H2_1),coefficients(H2_2)) +end + +let + @variables t + @parameters alpha = 2 + @species H(t) + + @compound H2_1(t) alpha*H + @compound H2_2(t) 2H + + @test iscompound(H2_1) + @test iscompound(H2_2) + + @test isequal(components(H2_1),components(H2_2)) + @test isequal(coefficients(H2_1), @parameters alpha = 2) +end + +let + @variables t + @species A(t) + B = A + @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)) +end + +#Check that balancing works. +let + @variables t + @parameters k + @species H(t) O(t) + @compound H2(t) 2H + @compound O2(t) 2O + @compound H2O(t) 2H 1O + + rx = Reaction(k,[H2,O2],[H2O]) + + @test isequal(create_matrix(rx),[2 0 -2; 0 2 -1; 0 0 1;]) + @test isequal(get_stoich(rx),[2, 1, 2]) + + balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) + @test isequal(balanced_rx, balance(rx)) + +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 + + rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) + + @test isequal(create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0; 0 0 0 1;]) + @test isequal(get_stoich(rx),[6, 6, 1, 6]) + + balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1,6]) + @test isequal(balanced_rx, balance(rx)) +end \ No newline at end of file From cf6f7c20a081dc6d80f8e1416a922522bd4fd1da Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Tue, 8 Aug 2023 21:19:25 +0530 Subject: [PATCH 04/37] Change file name --- src/{compound.jl => chemistry_functionality.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{compound.jl => chemistry_functionality.jl} (100%) diff --git a/src/compound.jl b/src/chemistry_functionality.jl similarity index 100% rename from src/compound.jl rename to src/chemistry_functionality.jl From 7745709b69e6d684f22ec9e8d927609b43cbc375 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Tue, 8 Aug 2023 21:34:47 +0530 Subject: [PATCH 05/37] Update Catalyst.jl --- src/Catalyst.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 995b0699f5..ae862b5b26 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -93,7 +93,7 @@ include("graphs.jl") export Graph, savegraph, complexgraph # for creating compounds -include("compound.jl") +include("chemistry_functionality.jl") export @compound export components, iscompound, coefficients export get_stoich, create_matrix, balance From b3a89803ab57da4512a66f9c7745f8cfef1c9e3f Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Thu, 10 Aug 2023 02:04:37 +0530 Subject: [PATCH 06/37] Add tests, modify balancing logic --- src/Catalyst.jl | 2 +- src/chemistry_functionality.jl | 75 ++++++++++- test/miscellaneous_tests/compound_macro.jl | 148 ++++++++++++++++++++- 3 files changed, 217 insertions(+), 8 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index ae862b5b26..02266fe3a7 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -96,6 +96,6 @@ export Graph, savegraph, complexgraph include("chemistry_functionality.jl") export @compound export components, iscompound, coefficients -export get_stoich, create_matrix, balance +export get_stoich, create_matrix, balance_reaction end # module diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 3e540393e9..d1cb17be8c 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -93,29 +93,35 @@ end function create_matrix(reaction::Catalyst.Reaction) compounds = [reaction.substrates; reaction.products] - atoms = [] + atoms = [] # Array to store unique atoms n_atoms = 0 A = zeros(Int, 0, length(compounds)) for (j, compound) in enumerate(compounds) + # Check if the compound is a valid compound if iscompound(compound) + # Get component coefficients of the compound pairs = component_coefficients(compound) if pairs == Nothing continue end else + # If not a compound, assume coefficient of 1 pairs = [(compound, 1)] end for pair in pairs + # Extract atom and coefficient from the pair atom, coeff = pair i = findfirst(x -> isequal(x, atom), atoms) if i === nothing + # Add the atom to the atoms array if it's not already present push!(atoms, atom) n_atoms += 1 A = [A; zeros(Int, 1, length(compounds))] i = n_atoms end + # Adjust coefficient based on whether the compound is a product or substrate coeff = any(map(p -> isequal(p, compounds[j]), reaction.products)) ? -coeff : coeff A[i, j] = coeff end @@ -154,6 +160,26 @@ function get_stoich(reaction::Reaction) # Perform the division X = A_transformed_rational \ B_transformed_rational + # Check if X is of Float64 type + if eltype(X) == Float64 + # Check if X contains 0.5 + if any(abs.(X .- 0.5) .< 1e-8) + # If so, apply custom rounding logic + X = X .+ 1 + for i in 1:length(X) + if abs(X[i] - floor(X[i])) >= 0.6 + X[i] = ceil(X[i]) + else + X[i] = floor(X[i]) + end + end + X = round.(Int, X) + else + # If not, just round normally + X = round.(Int, X) + end + end + # Get the denominators of the rational numbers in X denominators = denominator.(X) @@ -166,10 +192,10 @@ function get_stoich(reaction::Reaction) # Convert the rational numbers to integers X_integers = numerator.(X_multiplied) - return X_integers + return abs.(X_integers) end -function balance(reaction::Reaction) +function balance_reaction(reaction::Reaction) # Calculate the stoichiometric coefficients for the balanced reaction. stoichiometries = get_stoich(reaction) @@ -184,6 +210,49 @@ function balance(reaction::Reaction) return balanced_reaction end +# function get_stoich(reaction::Reaction) +# # Create the matrix A using create_matrix function. +# A = create_matrix(reaction) + +# # Create the vector b. The last element is 1 and others are 0. +# B = zeros(Int64, size(A,1)) +# B[end] = 1 + +# # Concatenate B to A to form an augmented matrix +# AB = [A B] + +# # Apply the Bareiss algorithm +# ModelingToolkit.bareiss!(AB) + +# # Extract the transformed A and B +# A_transformed = AB[:, 1:end-1] +# B_transformed = AB[:, end] + +# # Convert A_transformed to rational numbers +# A_transformed_rational = Rational.(A_transformed) + +# # Convert B_transformed to rational numbers +# B_transformed_rational = Rational.(B_transformed) + +# # Perform the division +# X = A_transformed_rational \ B_transformed_rational +# X = round.(Int, X) + +# # Get the denominators of the rational numbers in X +# denominators = denominator.(X) + +# # Compute the LCM of the denominators +# lcm_value = reduce(lcm, denominators) + +# # Multiply each element in X by the LCM of the denominators +# X_multiplied = X .* lcm_value + +# # Convert the rational numbers to integers +# X_integers = numerator.(X_multiplied) + +# return X_integers +# end + # function get_stoich(reaction::Reaction) ## Gaussian Elimination # # Create the matrix A using create_matrix function. # A = create_matrix(reaction) diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index c9a860b2a4..7b2d53cd2e 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -155,7 +155,7 @@ let @test isequal(get_stoich(rx),[2, 1, 2]) balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) - @test isequal(balanced_rx, balance(rx)) + @test isequal(balanced_rx, balance_reaction(rx)) end @@ -173,6 +173,146 @@ let @test isequal(create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0; 0 0 0 1;]) @test isequal(get_stoich(rx),[6, 6, 1, 6]) - balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1,6]) - @test isequal(balanced_rx, balance(rx)) -end \ No newline at end of file + balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1, 6]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, H2O --> H2O +let + @variables t + @species H(t) O(t) + @compound H2O(t) 2H O + + rx = Reaction(1.0, [H2O], [H2O], [2], [2]) + + balanced_rx = Reaction(1.0, [H2O], [H2O], [1], [1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2H + 1O --> H2O +let + @variables t + @species H(t) O(t) + @compound H2O(t) 2H O + + rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7]) + + balanced_rx = Reaction(1.0, [H,O], [H2O], [2, 1], [1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 1CH4 + 2O2 --> 1CO2 + 2H2O +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 + + rx = Reaction(1.0, [CH4, O2], [CO2, H2O]) + + balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, N2 + 3H2 --> 2NH3 +let + @variables t + @species H(t) N(t) + @compound N2(t) 2N + @compound H2(t) 2H + @compound NH3(t) N 3H + + rx = Reaction(1.0, [N2, H2], [NH3]) + + balanced_rx = Reaction(1.0, [N2, H2], [NH3], [1 ,3], [2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, C2H5OH + CH3COOH --> C4H8O2 + H2O +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 + + rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O]) + + balanced_rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O], [1, 1], [1, 1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2Ca3PO42 --> 6CaO + 1P4O10 +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 + + rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10]) + + balanced_rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10], [2], [6, 1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 4Fe + 3O2 + 6H2O --> 4FeOH3 +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 + + rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3]) + + balanced_rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3], [4, 3, 6], [4]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2NaOH + H2SO4 --> Na2SO4 + 2H2O +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 + + rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O]) + + balanced_rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O], [2, 1], [1, 2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2NO2 --> 1N2O4 +let + @variables t + @species N(t) O(t) + @compound NO2(t) N 2O + @compound N2O4(t) 2N 4O + + rx = Reaction(1.0, [NO2], [N2O4]) + + balanced_rx = Reaction(1.0, [NO2], [N2O4], [2], [1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 1CaCO3 + 2HCl --> CaCl2 + H2O + CO2 +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 + + rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O]) + balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end From bd7f43b00f0683660e0ca5a55654aba152586038 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Thu, 17 Aug 2023 19:44:52 +0530 Subject: [PATCH 07/37] Rational approach from start to end Adds backward_substitution function, eliminates rounding completely and avoids floating point numbers --- src/Catalyst.jl | 2 +- src/chemistry_functionality.jl | 57 +++++++++++----------- test/miscellaneous_tests/compound_macro.jl | 8 +-- 3 files changed, 33 insertions(+), 34 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 02266fe3a7..65d950aa46 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -96,6 +96,6 @@ export Graph, savegraph, complexgraph include("chemistry_functionality.jl") export @compound export components, iscompound, coefficients -export get_stoich, create_matrix, balance_reaction +export balance_reaction end # module diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index d1cb17be8c..9313e56135 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -126,10 +126,7 @@ function create_matrix(reaction::Catalyst.Reaction) A[i, j] = coeff end end - # Append a row with last element as 1 - new_row = zeros(Int, 1, size(A, 2)) - new_row[end] = 1 - A = vcat(A, new_row) + return A end @@ -137,9 +134,18 @@ function get_stoich(reaction::Reaction) # Create the matrix A using create_matrix function. A = create_matrix(reaction) - # Create the vector b. The last element is 1 and others are 0. - B = zeros(Int64, size(A,1)) - B[end] = 1 + m, n = size(A) + if m == n + B = zeros(Int64, size(A,1)) + else + # Append a row with last element as 1 + new_row = zeros(Int, 1, size(A, 2)) + new_row[end] = 1 + A = vcat(A, new_row) + + B = zeros(Int64, size(A,1)) + B[end] = 1 + end # Concatenate B to A to form an augmented matrix AB = [A B] @@ -157,28 +163,8 @@ function get_stoich(reaction::Reaction) # Convert B_transformed to rational numbers B_transformed_rational = Rational.(B_transformed) - # Perform the division - X = A_transformed_rational \ B_transformed_rational - - # Check if X is of Float64 type - if eltype(X) == Float64 - # Check if X contains 0.5 - if any(abs.(X .- 0.5) .< 1e-8) - # If so, apply custom rounding logic - X = X .+ 1 - for i in 1:length(X) - if abs(X[i] - floor(X[i])) >= 0.6 - X[i] = ceil(X[i]) - else - X[i] = floor(X[i]) - end - end - X = round.(Int, X) - else - # If not, just round normally - X = round.(Int, X) - end - end + # Perform backward substitution + X = backward_substitution(A_transformed_rational,B_transformed_rational) # Get the denominators of the rational numbers in X denominators = denominator.(X) @@ -210,6 +196,19 @@ function balance_reaction(reaction::Reaction) return balanced_reaction end +function backward_substitution(A::AbstractMatrix{T}, B::AbstractVector{T}) where T <: Number + n = length(B) + x = zeros(Rational{Int}, n) + for i in n:-1:1 + if all(A[i, :] .== 0) + x[i] = 1 + else + x[i] = (B[i] - A[i, i+1:n]' * x[i+1:n]) / A[i, i] + end + end + return x +end + # function get_stoich(reaction::Reaction) # # Create the matrix A using create_matrix function. # A = create_matrix(reaction) diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index 7b2d53cd2e..660d84fa80 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -151,8 +151,8 @@ let rx = Reaction(k,[H2,O2],[H2O]) - @test isequal(create_matrix(rx),[2 0 -2; 0 2 -1; 0 0 1;]) - @test isequal(get_stoich(rx),[2, 1, 2]) + @test isequal(Catalyst.create_matrix(rx),[2 0 -2; 0 2 -1;]) + @test isequal(Catalyst.get_stoich(rx),[2, 1, 2]) balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) @test isequal(balanced_rx, balance_reaction(rx)) @@ -170,8 +170,8 @@ let rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) - @test isequal(create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0; 0 0 0 1;]) - @test isequal(get_stoich(rx),[6, 6, 1, 6]) + @test isequal(Catalyst.create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0;]) + @test isequal(Catalyst.get_stoich(rx),[6, 6, 1, 6]) balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1, 6]) @test isequal(balanced_rx, balance_reaction(rx)) From 58e671899c5b9d9366c7100d314806e5ff68ba85 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Thu, 17 Aug 2023 21:05:56 +0530 Subject: [PATCH 08/37] Separate balancing test file --- test/miscellaneous_tests/compound_macro.jl | 176 ----------------- .../miscellaneous_tests/reaction_balancing.jl | 178 ++++++++++++++++++ test/runtests.jl | 1 + 3 files changed, 179 insertions(+), 176 deletions(-) create mode 100644 test/miscellaneous_tests/reaction_balancing.jl diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index 660d84fa80..be3afbb035 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -140,179 +140,3 @@ let @test isequal(coefficients(A2), coefficients(B2)) end -#Check that balancing works. -let - @variables t - @parameters k - @species H(t) O(t) - @compound H2(t) 2H - @compound O2(t) 2O - @compound H2O(t) 2H 1O - - rx = Reaction(k,[H2,O2],[H2O]) - - @test isequal(Catalyst.create_matrix(rx),[2 0 -2; 0 2 -1;]) - @test isequal(Catalyst.get_stoich(rx),[2, 1, 2]) - - balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) - @test isequal(balanced_rx, balance_reaction(rx)) - -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 - - rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) - - @test isequal(Catalyst.create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0;]) - @test isequal(Catalyst.get_stoich(rx),[6, 6, 1, 6]) - - balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1, 6]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, H2O --> H2O -let - @variables t - @species H(t) O(t) - @compound H2O(t) 2H O - - rx = Reaction(1.0, [H2O], [H2O], [2], [2]) - - balanced_rx = Reaction(1.0, [H2O], [H2O], [1], [1]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, 2H + 1O --> H2O -let - @variables t - @species H(t) O(t) - @compound H2O(t) 2H O - - rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7]) - - balanced_rx = Reaction(1.0, [H,O], [H2O], [2, 1], [1]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, 1CH4 + 2O2 --> 1CO2 + 2H2O -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 - - rx = Reaction(1.0, [CH4, O2], [CO2, H2O]) - - balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, N2 + 3H2 --> 2NH3 -let - @variables t - @species H(t) N(t) - @compound N2(t) 2N - @compound H2(t) 2H - @compound NH3(t) N 3H - - rx = Reaction(1.0, [N2, H2], [NH3]) - - balanced_rx = Reaction(1.0, [N2, H2], [NH3], [1 ,3], [2]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, C2H5OH + CH3COOH --> C4H8O2 + H2O -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 - - rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O]) - - balanced_rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O], [1, 1], [1, 1]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, 2Ca3PO42 --> 6CaO + 1P4O10 -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 - - rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10]) - - balanced_rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10], [2], [6, 1]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, 4Fe + 3O2 + 6H2O --> 4FeOH3 -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 - - rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3]) - - balanced_rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3], [4, 3, 6], [4]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, 2NaOH + H2SO4 --> Na2SO4 + 2H2O -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 - - rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O]) - - balanced_rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O], [2, 1], [1, 2]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, 2NO2 --> 1N2O4 -let - @variables t - @species N(t) O(t) - @compound NO2(t) N 2O - @compound N2O4(t) 2N 4O - - rx = Reaction(1.0, [NO2], [N2O4]) - - balanced_rx = Reaction(1.0, [NO2], [N2O4], [2], [1]) - @test isequal(balanced_rx, balance_reaction(rx)) -end - -# @reaction k, 1CaCO3 + 2HCl --> CaCl2 + H2O + CO2 -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 - - rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O]) - balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1]) - @test isequal(balanced_rx, balance_reaction(rx)) -end diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl new file mode 100644 index 0000000000..d259bd7f57 --- /dev/null +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -0,0 +1,178 @@ +using Catalyst, Test + +#Check that balancing works. +let + @variables t + @parameters k + @species H(t) O(t) + @compound H2(t) 2H + @compound O2(t) 2O + @compound H2O(t) 2H 1O + + rx = Reaction(k,[H2,O2],[H2O]) + + @test isequal(Catalyst.create_matrix(rx),[2 0 -2; 0 2 -1;]) + @test isequal(Catalyst.get_stoich(rx),[2, 1, 2]) + + balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) + @test isequal(balanced_rx, balance_reaction(rx)) + +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 + + rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) + + @test isequal(Catalyst.create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0;]) + @test isequal(Catalyst.get_stoich(rx),[6, 6, 1, 6]) + + balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1, 6]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, H2O --> H2O +let + @variables t + @species H(t) O(t) + @compound H2O(t) 2H O + + rx = Reaction(1.0, [H2O], [H2O], [2], [2]) + + balanced_rx = Reaction(1.0, [H2O], [H2O], [1], [1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2H + 1O --> H2O +let + @variables t + @species H(t) O(t) + @compound H2O(t) 2H O + + rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7]) + + balanced_rx = Reaction(1.0, [H,O], [H2O], [2, 1], [1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 1CH4 + 2O2 --> 1CO2 + 2H2O +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 + + rx = Reaction(1.0, [CH4, O2], [CO2, H2O]) + + balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, N2 + 3H2 --> 2NH3 +let + @variables t + @species H(t) N(t) + @compound N2(t) 2N + @compound H2(t) 2H + @compound NH3(t) N 3H + + rx = Reaction(1.0, [N2, H2], [NH3]) + + balanced_rx = Reaction(1.0, [N2, H2], [NH3], [1 ,3], [2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, C2H5OH + CH3COOH --> C4H8O2 + H2O +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 + + rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O]) + + balanced_rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O], [1, 1], [1, 1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2Ca3PO42 --> 6CaO + 1P4O10 +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 + + rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10]) + + balanced_rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10], [2], [6, 1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 4Fe + 3O2 + 6H2O --> 4FeOH3 +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 + + rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3]) + + balanced_rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3], [4, 3, 6], [4]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2NaOH + H2SO4 --> Na2SO4 + 2H2O +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 + + rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O]) + + balanced_rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O], [2, 1], [1, 2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2NO2 --> 1N2O4 +let + @variables t + @species N(t) O(t) + @compound NO2(t) N 2O + @compound N2O4(t) 2N 4O + + rx = Reaction(1.0, [NO2], [N2O4]) + + balanced_rx = Reaction(1.0, [NO2], [N2O4], [2], [1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 1CaCO3 + 2HCl --> CaCl2 + H2O + CO2 +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 + + rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O]) + balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end diff --git a/test/runtests.jl b/test/runtests.jl index 2bbfa84495..5a772722f2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,6 +23,7 @@ using SafeTestsets @time @safetestset "Symbolic Stoichiometry" begin include("miscellaneous_tests/symbolic_stoichiometry.jl") end @time @safetestset "Events" begin include("miscellaneous_tests/events.jl") end @time @safetestset "Compound species" begin include("miscellaneous_tests/compound_macro.jl") end + @time @safetestset "Reaction balancing" begin include("miscellaneous_tests/reaction_balancing.jl") end @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end ### Reaction network analysis. ### From a3e6a3fb3c20951d72d9804d202cfddf5c1a054b Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Sat, 19 Aug 2023 00:02:28 +0530 Subject: [PATCH 09/37] Change filename back to previous name --- ...chemistry_functionality.jl => compound.jl} | 119 ------------------ 1 file changed, 119 deletions(-) rename src/{chemistry_functionality.jl => compound.jl} (67%) diff --git a/src/chemistry_functionality.jl b/src/compound.jl similarity index 67% rename from src/chemistry_functionality.jl rename to src/compound.jl index 9313e56135..e66372a5e6 100644 --- a/src/chemistry_functionality.jl +++ b/src/compound.jl @@ -209,122 +209,3 @@ function backward_substitution(A::AbstractMatrix{T}, B::AbstractVector{T}) where return x end -# function get_stoich(reaction::Reaction) -# # Create the matrix A using create_matrix function. -# A = create_matrix(reaction) - -# # Create the vector b. The last element is 1 and others are 0. -# B = zeros(Int64, size(A,1)) -# B[end] = 1 - -# # Concatenate B to A to form an augmented matrix -# AB = [A B] - -# # Apply the Bareiss algorithm -# ModelingToolkit.bareiss!(AB) - -# # Extract the transformed A and B -# A_transformed = AB[:, 1:end-1] -# B_transformed = AB[:, end] - -# # Convert A_transformed to rational numbers -# A_transformed_rational = Rational.(A_transformed) - -# # Convert B_transformed to rational numbers -# B_transformed_rational = Rational.(B_transformed) - -# # Perform the division -# X = A_transformed_rational \ B_transformed_rational -# X = round.(Int, X) - -# # Get the denominators of the rational numbers in X -# denominators = denominator.(X) - -# # Compute the LCM of the denominators -# lcm_value = reduce(lcm, denominators) - -# # Multiply each element in X by the LCM of the denominators -# X_multiplied = X .* lcm_value - -# # Convert the rational numbers to integers -# X_integers = numerator.(X_multiplied) - -# return X_integers -# end - -# function get_stoich(reaction::Reaction) ## Gaussian Elimination -# # Create the matrix A using create_matrix function. -# A = create_matrix(reaction) - -# # Create the vector b. The last element is 1 and others are 0. -# B = zeros(Int64, size(A,1)) -# B[end] = 1 - -# # Convert A and B to sparse matrices -# A = sparse(A) - -# # Apply the LU factorization (Gaussian elimination) -# lu_f = lu(A) - -# # Convert B to a dense vector -# B_dense = Vector(B) - -# # Solve for X using the LU factorization -# X = lu_f \ B_dense - -# # Get the smallest positive value in X -# smallest_value = minimum(X[X .> 0]) - -# # Normalize X to be integers -# X_normalized = round.(Int64, X / smallest_value) - -# return X_normalized -# end - -# function get_stoich(reaction::Reaction) ##Bareiss Algorithm (Does not avoid floating points) -# # Create the matrix A using create_matrix function. -# A = create_matrix(reaction) - -# # Create the vector b. The last element is 1 and others are 0. -# B = zeros(Int64, size(A,1)) -# B[end] = 1 - -# # Concatenate B to A to form an augmented matrix -# AB = [A B] - -# # Apply the Bareiss algorithm -# ModelingToolkit.bareiss!(AB) - -# # Extract the transformed A and B -# A_transformed = AB[:, 1:end-1] -# B_transformed = AB[:, end] - -# # Solve for X using back substitution -# X = A_transformed \ B_transformed - -# # Get the smallest positive value in X -# smallest_value = minimum(X[X .> 0]) - -# # Normalize X to be integers -# X_normalized = round.(Int64, X / smallest_value) - -# return X_normalized -# end - -# @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 -# rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) - -# using LinearAlgebra -# using SparseArrays -# using SuiteSparse.UMFPACK - -# gcd_value = gcd(X...) -# X_scaled = X ./ gcd_value - -# X = [1,1,0.2,1] \ No newline at end of file From 89c72f4799135ac433dfd475517a24f9326ebaa6 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Sat, 19 Aug 2023 00:02:54 +0530 Subject: [PATCH 10/37] Update Catalyst.jl --- src/Catalyst.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 65d950aa46..1ef2719827 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -93,7 +93,7 @@ include("graphs.jl") export Graph, savegraph, complexgraph # for creating compounds -include("chemistry_functionality.jl") +include("compound.jl") export @compound export components, iscompound, coefficients export balance_reaction From c01c87038b2b72beceb8a03f01fbfbcb2d7ca8aa Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Sat, 19 Aug 2023 00:09:41 +0530 Subject: [PATCH 11/37] Final renaming file --- src/{compound.jl => chemistry_functionality.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{compound.jl => chemistry_functionality.jl} (100%) diff --git a/src/compound.jl b/src/chemistry_functionality.jl similarity index 100% rename from src/compound.jl rename to src/chemistry_functionality.jl From 5475f9692dd6b16808bf4df7b3edc2b90ab5fed7 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Sat, 19 Aug 2023 00:20:40 +0530 Subject: [PATCH 12/37] Update Catalyst.jl --- src/Catalyst.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 1ef2719827..65d950aa46 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -93,7 +93,7 @@ include("graphs.jl") export Graph, savegraph, complexgraph # for creating compounds -include("compound.jl") +include("chemistry_functionality.jl") export @compound export components, iscompound, coefficients export balance_reaction From 05b4264a0a7b67f718eff5c9edcbeb890326cd19 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Tue, 22 Aug 2023 23:24:02 +0530 Subject: [PATCH 13/37] Add more tests, optimize logic --- src/chemistry_functionality.jl | 3 + .../miscellaneous_tests/reaction_balancing.jl | 126 ++++++++++++++++++ 2 files changed, 129 insertions(+) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index e66372a5e6..243a585beb 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -137,6 +137,9 @@ function get_stoich(reaction::Reaction) m, n = size(A) if m == n B = zeros(Int64, size(A,1)) + elseif m > n + A = A[1:n, :] + B = zeros(Int64, size(A,1)) else # Append a row with last element as 1 new_row = zeros(Int, 1, size(A, 2)) diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index d259bd7f57..f6bbcc5669 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -176,3 +176,129 @@ let balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1]) @test isequal(balanced_rx, balance_reaction(rx)) end + +# @reaction k, SiCl4 + 4H2O → H4SiO4 + 4HCl +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 + + rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl]) + balanced_rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl], [1,4], [1,4]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2Al + 6HCl → 2AlCl3 + 3H2 +let + @variables t + @species Al(t) Cl(t) H(t) + @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]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, Na2CO3 + 2HCl → 2NaCl + H2O + CO2 +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 + + rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2]) + balanced_rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2], [1,2], [2,1,1]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2C7H6O2 + 15O2 → 14CO2 + 6H2O +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 + + rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O]) + balanced_rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O], [2,15], [14,6]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k,Fe2(SO4)3 + 6KOH → 3K2SO4 + 2Fe(OH)3 +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 + + rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3]) #5x4 matrix + balanced_rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3], [1,6], [3,2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 2Ca3(PO4)2 + 6SiO2 → P4O10 + 6CaSiO3 +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 + + rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3]) #5x4 matrix + balanced_rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3], [2,6] , [1,6]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, 4KClO3 → 3KClO4 + KCl +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 + + rx = Reaction(1.0,[KClO3],[KClO4,KCl]) + balanced_rx = Reaction(1.0,[KClO3],[KClO4,KCl], [4], [3,1]) + @test isequal(balanced_rx, balance_reaction(rx)) +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) + @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]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# @reaction k, H2SO4 + 8HI → H2S + 4I2 + 4H2O +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 + + rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O]) + balanced_rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O], [1,8], [1,4,4]) + @test isequal(balanced_rx, balance_reaction(rx)) +end From 377698dd89d12d1b5f0145cb8cb51e6113941740 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Tue, 29 Aug 2023 22:59:29 +0530 Subject: [PATCH 14/37] Use nullspace to get stoichiometries --- src/chemistry_functionality.jl | 81 ++++++------------- .../miscellaneous_tests/reaction_balancing.jl | 38 +++++++++ 2 files changed, 62 insertions(+), 57 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 243a585beb..2409992147 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -134,54 +134,21 @@ function get_stoich(reaction::Reaction) # Create the matrix A using create_matrix function. A = create_matrix(reaction) - m, n = size(A) - if m == n - B = zeros(Int64, size(A,1)) - elseif m > n - A = A[1:n, :] - B = zeros(Int64, size(A,1)) - else - # Append a row with last element as 1 - new_row = zeros(Int, 1, size(A, 2)) - new_row[end] = 1 - A = vcat(A, new_row) - - B = zeros(Int64, size(A,1)) - B[end] = 1 - end - - # Concatenate B to A to form an augmented matrix - AB = [A B] - - # Apply the Bareiss algorithm - ModelingToolkit.bareiss!(AB) - - # Extract the transformed A and B - A_transformed = AB[:, 1:end-1] - B_transformed = AB[:, end] - - # Convert A_transformed to rational numbers - A_transformed_rational = Rational.(A_transformed) - - # Convert B_transformed to rational numbers - B_transformed_rational = Rational.(B_transformed) - - # Perform backward substitution - X = backward_substitution(A_transformed_rational,B_transformed_rational) + X = ModelingToolkit.nullspace(A) - # Get the denominators of the rational numbers in X - denominators = denominator.(X) + m, n = size(X) + if n == 1 + X = abs.(vec(X)) + common_divisor = reduce(gcd, X) + X = X ./ common_divisor - # Compute the LCM of the denominators - lcm_value = reduce(lcm, denominators) - - # Multiply each element in X by the LCM of the denominators - X_multiplied = X .* lcm_value - - # Convert the rational numbers to integers - X_integers = numerator.(X_multiplied) + elseif n > 1 + error("Chemical equation can be balanced in infinitely many ways") + else + error("Chemical equation cannot be balanced") + end - return abs.(X_integers) + return X end function balance_reaction(reaction::Reaction) @@ -199,16 +166,16 @@ function balance_reaction(reaction::Reaction) return balanced_reaction end -function backward_substitution(A::AbstractMatrix{T}, B::AbstractVector{T}) where T <: Number - n = length(B) - x = zeros(Rational{Int}, n) - for i in n:-1:1 - if all(A[i, :] .== 0) - x[i] = 1 - else - x[i] = (B[i] - A[i, i+1:n]' * x[i+1:n]) / A[i, i] - end - end - return x -end +# function backward_substitution(A::AbstractMatrix{T}, B::AbstractVector{T}) where T <: Number +# n = length(B) +# x = zeros(Rational{Int}, n) +# for i in n:-1:1 +# if all(A[i, :] .== 0) +# x[i] = 1 +# else +# x[i] = (B[i] - A[i, i+1:n]' * x[i+1:n]) / A[i, i] +# end +# end +# return x +# end diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index f6bbcc5669..7093c09b2c 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -302,3 +302,41 @@ let balanced_rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O], [1,8], [1,4,4]) @test isequal(balanced_rx, balance_reaction(rx)) end + +# @reaction k, C2H4 + 3O2 ↔ 2CO2 + 2H2O +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 + + rx = Reaction(1.0,[C2H4,O2],[CO2,H2O]) + balanced_rx = Reaction(1.0,[C2H4,O2],[CO2,H2O],[1,3],[2,2]) + @test isequal(balanced_rx, balance_reaction(rx)) +end + +# Infinite solutions +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 + + rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O]) + + @test_throws ErrorException("Chemical equation can be balanced in infinitely many ways") balance_reaction(rx) +end + +# No way to balance +let + @variables t + @species Fe(t) S(t) O(t) H(t) N(t) + + rx = Reaction(1.0,[FeS2,HNO3],[Fe2S3O12,NO,H2SO4]) + @test_throws ErrorException("Chemical equation cannot be balanced") balance_reaction(rx) +end From 0db24ed8aad3990b57dca573cb9a411cce3b939b Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Tue, 29 Aug 2023 23:20:36 +0530 Subject: [PATCH 15/37] Fix test --- test/miscellaneous_tests/reaction_balancing.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index 7093c09b2c..e0681b9ec0 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -337,6 +337,12 @@ 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 + rx = Reaction(1.0,[FeS2,HNO3],[Fe2S3O12,NO,H2SO4]) @test_throws ErrorException("Chemical equation cannot be balanced") balance_reaction(rx) end From 519b99ff2f4e23671d876de9ea540c9611ae25a4 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Wed, 30 Aug 2023 19:42:46 +0530 Subject: [PATCH 16/37] Return a vector of Reactions when the nullspace dimension is greater than 1 --- src/chemistry_functionality.jl | 52 +++++++++++++++---- .../miscellaneous_tests/reaction_balancing.jl | 4 +- 2 files changed, 44 insertions(+), 12 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 2409992147..61b90919de 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -141,31 +141,63 @@ function get_stoich(reaction::Reaction) X = abs.(vec(X)) common_divisor = reduce(gcd, X) X = X ./ common_divisor - + return X elseif n > 1 - error("Chemical equation can be balanced in infinitely many ways") + # initialize final_vector + final_vector = Vector{Vector{Int64}}() + for i = 1:min(m, n) + X_i = abs.(vec(X[:, i])) + common_divisor = reduce(gcd, X_i) + X_i = X_i ./ common_divisor + push!(final_vector, X_i) + end + return final_vector else error("Chemical equation cannot be balanced") end - return X end function balance_reaction(reaction::Reaction) # Calculate the stoichiometric coefficients for the balanced reaction. stoichiometries = get_stoich(reaction) - # Divide the stoichiometry vector into substrate and product stoichiometries. - substoich = stoichiometries[1:length(reaction.substrates)] - prodstoich = stoichiometries[(length(reaction.substrates)) + 1:end] + # Check if stoichiometries is a vector of vectors + if typeof(stoichiometries[1]) == Vector{Int64} + # Initialize an empty vector to store all the reactions + possible_reactions = Vector{Reaction}() + + # Iterate over each stoichiometry vector and create a reaction + for stoich in stoichiometries + # Divide the stoichiometry vector into substrate and product stoichiometries. + substoich = stoich[1:length(reaction.substrates)] + prodstoich = stoich[(length(reaction.substrates)) + 1:end] - # Create a new reaction with the balanced stoichiometries - balanced_reaction = Reaction(reaction.rate, reaction.substrates, reaction.products, substoich, prodstoich) + # Create a new reaction with the balanced stoichiometries + balanced_reaction = Reaction(reaction.rate, reaction.substrates, reaction.products, substoich, prodstoich) + + # Add the reaction to the vector of all reactions + push!(possible_reactions, balanced_reaction) + end - # Return the balanced reaction - return balanced_reaction + println("Chemical equation $reaction can be balanced in infinitely many ways") + # Return the vector of all reactions + return possible_reactions + + else + # Divide the stoichiometry vector into substrate and product stoichiometries. + substoich = stoichiometries[1:length(reaction.substrates)] + prodstoich = stoichiometries[(length(reaction.substrates)) + 1:end] + + # Create a new reaction with the balanced stoichiometries + balanced_reaction = Reaction(reaction.rate, reaction.substrates, reaction.products, substoich, prodstoich) + + # Return the balanced reaction + return balanced_reaction + end end + # function backward_substitution(A::AbstractMatrix{T}, B::AbstractVector{T}) where T <: Number # n = length(B) # x = zeros(Rational{Int}, n) diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index e0681b9ec0..c177c1b604 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -328,8 +328,8 @@ let @compound H2O(t) 2H O rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O]) - - @test_throws ErrorException("Chemical equation can be balanced in infinitely many ways") balance_reaction(rx) + @test typeof(balance_reaction(rx)) == Vector{Reaction} + @test length(balance_reaction(rx)) == 2 end # No way to balance From 5828a2cd491ee0ca2f090d918338c8bae5684e30 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Fri, 1 Sep 2023 13:41:05 -0400 Subject: [PATCH 17/37] Fix namespacing --- src/Catalyst.jl | 9 ++++++--- src/compound.jl | 16 ++++++++-------- src/reaction_network.jl | 3 ++- src/reactionsystem.jl | 6 ++++-- .../component_based_model_creation.jl | 4 ++-- 5 files changed, 22 insertions(+), 16 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 1757e0241d..4aa78ca7f2 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -6,7 +6,8 @@ module Catalyst using DocStringExtensions using SparseArrays, DiffEqBase, Reexport using LaTeXStrings, Latexify, Requires -using JumpProcesses: JumpProcesses, JumpProblem, MassActionJump, ConstantRateJump, +using JumpProcesses: JumpProcesses, + JumpProblem, MassActionJump, ConstantRateJump, VariableRateJump # ModelingToolkit imports and convenience functions we use @@ -26,7 +27,8 @@ import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_v namespace_parameters, rename, renamespace, getname, flatten # internal but needed ModelingToolkit functions -import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, check_units, +import ModelingToolkit: check_variables, + check_parameters, _iszero, _merge, check_units, get_unit, check_equations import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show @@ -55,7 +57,8 @@ end include("reactionsystem.jl") export isspecies export Reaction, ReactionSystem, ismassaction, oderatelaw, jumpratelaw, isspatial -export ODEProblem, SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem, +export ODEProblem, + SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem, SteadyStateProblem # reaction_network macro diff --git a/src/compound.jl b/src/compound.jl index 9f2ac5bca2..01c9081232 100644 --- a/src/compound.jl +++ b/src/compound.jl @@ -18,12 +18,12 @@ macro compound(species_expr, arr_expr...) # Construct the expressions that define the species species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0), - Expr(:call, species_name, species_arg)) + 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)) + Catalyst.CompoundSpecies, + true)) # Ensure the expressions are evaluated in the correct scope by escaping them escaped_species_expr = esc(species_expr) @@ -49,22 +49,22 @@ macro compound(species_expr, arr_expr...) # Construct the expression to set the components metadata setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundComponents, - $species_expr)) + Catalyst.CompoundComponents, + $species_expr)) # Ensure the expression is evaluated in the correct scope by escaping it escaped_setcomponents_expr = esc(setcomponents_expr) # Construct the expression to set the coefficients metadata setcoefficients_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), - Catalyst.CompoundCoefficients, - $coeffs_expr)) + Catalyst.CompoundCoefficients, + $coeffs_expr)) escaped_setcoefficients_expr = esc(setcoefficients_expr) # Return a block that contains the escaped expressions return Expr(:block, escaped_species_expr, escaped_setmetadata_expr, - escaped_setcomponents_expr, escaped_setcoefficients_expr) + escaped_setcomponents_expr, escaped_setcoefficients_expr) end # Check if a species is a compound diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 05812ac3c7..46c099cc0c 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -316,7 +316,8 @@ end function forbidden_symbol_check(v) !isempty(intersect(forbidden_symbols_error, v)) && error("The following symbol(s) are used as species or parameters: " * - ((map(s -> "'" * string(s) * "', ", intersect(forbidden_symbols_error, v))...)) * + ((map(s -> "'" * string(s) * "', ", + intersect(forbidden_symbols_error, v))...)) * "this is not permited.") nothing end diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 77dd98710b..adb89ea02a 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -228,7 +228,7 @@ function apply_if_nonempty(f, v) s end -function ModelingToolkit.namespace_equation(rx::Reaction, name) +function ModelingToolkit.namespace_equation(rx::Reaction, name; kw...) f = Base.Fix2(namespace_expr, name) rate = f(rx.rate) subs = apply_if_nonempty(f, rx.substrates) @@ -780,11 +780,13 @@ MT.get_continuous_events(sys::ReactionSystem) = getfield(sys, :continuous_events # need a custom equations since ReactionSystem.eqs are a mix of Reactions and Equations function MT.equations(sys::ReactionSystem) + ivs = independent_variables(sys) eqs = get_eqs(sys) systems = get_systems(sys) if !isempty(systems) eqs = CatalystEqType[eqs; - reduce(vcat, MT.namespace_equations.(systems); init = Any[])] + reduce(vcat, MT.namespace_equations.(systems, (ivs,)); + init = Any[])] return sort!(eqs; by = eqsortby) end return eqs diff --git a/test/programmatic_model_creation/component_based_model_creation.jl b/test/programmatic_model_creation/component_based_model_creation.jl index 0e66eb9e55..19382d8988 100644 --- a/test/programmatic_model_creation/component_based_model_creation.jl +++ b/test/programmatic_model_creation/component_based_model_creation.jl @@ -138,7 +138,7 @@ let # and after conversion to an AbstractSystem system = convert(NonlinearSystem, extended) @test isequal(system.a, ModelingToolkit.namespace_expr(a, system)) - @test isequal(system.x, ModelingToolkit.namespace_expr(x, system)) + @test isequal(system.x, ModelingToolkit.namespace_expr(x, system; ivs = independent_variables(extended))) @test length(equations(system)) == 1 @test equations(system) == equations(constraints) @@ -150,7 +150,7 @@ let # and after conversion to an AbstractSystem. system = convert(NonlinearSystem, extended) @test isequal(system.a, ModelingToolkit.namespace_expr(a, system)) - @test isequal(system.x, ModelingToolkit.namespace_expr(x, system)) + @test isequal(system.x, ModelingToolkit.namespace_expr(x, system; ivs = independent_variables(extended))) @test length(equations(system)) == 1 @test Set(equations(system)) == Set(equations(constraints)) From 062acad09abec969a202475af20c9968c09e0388 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Fri, 1 Sep 2023 17:19:52 -0400 Subject: [PATCH 18/37] Bound MTK --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 3c9f94dfd1..e756a1ae2d 100644 --- a/Project.toml +++ b/Project.toml @@ -30,7 +30,7 @@ JumpProcesses = "9.3.2" LaTeXStrings = "1.3.0" Latexify = "0.14, 0.15, 0.16" MacroTools = "0.5.5" -ModelingToolkit = "8.47" +ModelingToolkit = "8.66" Parameters = "0.12" Reexport = "0.2, 1.0" Requires = "1.0" From 9d74a24adbccff0b04d2ac0b0d757eacccc08bdd Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 1 Sep 2023 18:28:40 -0400 Subject: [PATCH 19/37] cleanup code and handle edge cases --- src/chemistry_functionality.jl | 170 ++++++++--------- .../miscellaneous_tests/reaction_balancing.jl | 176 ++++++++++++------ 2 files changed, 194 insertions(+), 152 deletions(-) diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index 61b90919de..e7bf9fbdcd 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -90,124 +90,108 @@ end ### Balancing Code +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 function create_matrix(reaction::Catalyst.Reaction) - compounds = [reaction.substrates; reaction.products] - atoms = [] # Array to store unique atoms + @unpack substrates, products = reaction + unique_atoms = [] # Array to store unique atoms n_atoms = 0 - A = zeros(Int, 0, length(compounds)) - - for (j, compound) in enumerate(compounds) - # Check if the compound is a valid compound - if iscompound(compound) - # Get component coefficients of the compound - pairs = component_coefficients(compound) - if pairs == Nothing - continue + ncompounds = length(substrates) + length(products) + A = zeros(Int, 0, ncompounds) + + coeffsign = 1 + jbase = 0 + for compounds in (substrates, products) + for (j2, compound) in enumerate(compounds) + j = jbase + j2 + + if iscompound(compound) + atoms = components(compound) + any(iscompound, atoms) && throw(COMPOUND_OF_COMPOUND_ERROR) + coeffs = coefficients(compound) + (atoms == nothing || coeffs == nothing) && continue + else + # If not a compound, assume coefficient of 1 + atoms = [compound] + coeffs = [1] end - else - # If not a compound, assume coefficient of 1 - pairs = [(compound, 1)] - end - for pair in pairs - # Extract atom and coefficient from the pair - atom, coeff = pair - i = findfirst(x -> isequal(x, atom), atoms) - if i === nothing - # Add the atom to the atoms array if it's not already present - push!(atoms, atom) - n_atoms += 1 - A = [A; zeros(Int, 1, length(compounds))] - i = n_atoms + for (atom,coeff) in zip(atoms, coeffs) + # Extract atom and coefficient from the pair + i = findfirst(x -> isequal(x, atom), unique_atoms) + if i === nothing + # Add the atom to the atoms array if it's not already present + push!(unique_atoms, atom) + n_atoms += 1 + A = [A; zeros(Int, 1, ncompounds)] + i = n_atoms + end + + # Adjust coefficient based on whether the compound is a product or substrate + coeff *= coeffsign + + A[i, j] = coeff end - # Adjust coefficient based on whether the compound is a product or substrate - coeff = any(map(p -> isequal(p, compounds[j]), reaction.products)) ? -coeff : coeff - A[i, j] = coeff end + + # update for iterating through products + coeffsign = -1 + jbase = length(substrates) end return A end -function get_stoich(reaction::Reaction) - # Create the matrix A using create_matrix function. +function get_balanced_stoich(reaction::Reaction) + # Create the reaction matrix A that is m atoms by n compounds A = create_matrix(reaction) - - X = ModelingToolkit.nullspace(A) - m, n = size(X) - if n == 1 - X = abs.(vec(X)) - common_divisor = reduce(gcd, X) - X = X ./ common_divisor - return X - elseif n > 1 - # initialize final_vector - final_vector = Vector{Vector{Int64}}() - for i = 1:min(m, n) - X_i = abs.(vec(X[:, i])) - common_divisor = reduce(gcd, X_i) - X_i = X_i ./ common_divisor - push!(final_vector, X_i) - end - return final_vector - else - error("Chemical equation cannot be balanced") + # get an integer nullspace basis + X = ModelingToolkit.nullspace(A) + nullity = size(X, 2) + + stoichvecs = Vector{Vector{Int64}}() + for (j, col) in pairs(eachcol(X)) + signs = unique(col) + + # If there is only one basis vector and the signs are not all the same this means + # we have a solution that would require moving at least one substrate to be a + # product (or vice-versa). We therefore do not return anything in this case. + # If there are multiple basis vectors we don't currently determine if we can + # construct a linear combination giving a solution, so we just return them. + if (nullity > 1) || (all(>=(0), signs) || all(<=(0), signs)) + coefs = abs.(col) + common_divisor = reduce(gcd, coefs) + coefs .= div.(coefs, common_divisor) + push!(stoichvecs, coefs) end + end + return stoichvecs end function balance_reaction(reaction::Reaction) # Calculate the stoichiometric coefficients for the balanced reaction. - stoichiometries = get_stoich(reaction) - - # Check if stoichiometries is a vector of vectors - if typeof(stoichiometries[1]) == Vector{Int64} - # Initialize an empty vector to store all the reactions - possible_reactions = Vector{Reaction}() + stoichiometries = get_balanced_stoich(reaction) - # Iterate over each stoichiometry vector and create a reaction - for stoich in stoichiometries - # Divide the stoichiometry vector into substrate and product stoichiometries. - substoich = stoich[1:length(reaction.substrates)] - prodstoich = stoich[(length(reaction.substrates)) + 1:end] + balancedrxs = Vector{Reaction}(undef, length(stoichiometries)) - # Create a new reaction with the balanced stoichiometries - balanced_reaction = Reaction(reaction.rate, reaction.substrates, reaction.products, substoich, prodstoich) - - # Add the reaction to the vector of all reactions - push!(possible_reactions, balanced_reaction) - end - - println("Chemical equation $reaction can be balanced in infinitely many ways") - # Return the vector of all reactions - return possible_reactions - - else + # Iterate over each stoichiometry vector and create a reaction + for (i,stoich) in enumerate(stoichiometries) # Divide the stoichiometry vector into substrate and product stoichiometries. - substoich = stoichiometries[1:length(reaction.substrates)] - prodstoich = stoichiometries[(length(reaction.substrates)) + 1:end] + substoich = stoich[1:length(reaction.substrates)] + prodstoich = stoich[(length(reaction.substrates) + 1):end] # Create a new reaction with the balanced stoichiometries - balanced_reaction = Reaction(reaction.rate, reaction.substrates, reaction.products, substoich, prodstoich) + balancedrx = Reaction(reaction.rate, reaction.substrates, + reaction.products, substoich, prodstoich) - # Return the balanced reaction - return balanced_reaction + # Add the reaction to the vector of all reactions + balancedrxs[i] = balancedrx end -end - - -# function backward_substitution(A::AbstractMatrix{T}, B::AbstractVector{T}) where T <: Number -# n = length(B) -# x = zeros(Rational{Int}, n) -# for i in n:-1:1 -# if all(A[i, :] .== 0) -# x[i] = 1 -# else -# x[i] = (B[i] - A[i, i+1:n]' * x[i+1:n]) / A[i, i] -# end -# end -# return x -# end + isempty(balancedrxs) && (@warn "Unable to balance reaction.") + (length(balancedrxs) > 1) && (@warn "Infinite balanced reactions from ($reaction) are possible, returning a basis for them. Note that we do not check if they preserve the set of substrates and products from the original reaction.") + return balancedrxs +end diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index c177c1b604..65cb0411f7 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -1,7 +1,7 @@ using Catalyst, Test #Check that balancing works. -let +let @variables t @parameters k @species H(t) O(t) @@ -12,17 +12,18 @@ let rx = Reaction(k,[H2,O2],[H2O]) @test isequal(Catalyst.create_matrix(rx),[2 0 -2; 0 2 -1;]) - @test isequal(Catalyst.get_stoich(rx),[2, 1, 2]) + @test isequal(Catalyst.get_balanced_stoich(rx),[[2, 1, 2]]) balanced_rx = Reaction(k,[H2,O2],[H2O],[2,1],[2]) - @test isequal(balanced_rx, balance_reaction(rx)) - + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end -let +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 @@ -31,14 +32,16 @@ let rx = Reaction(k,[CO2,H2O],[C6H12O6,O2]) @test isequal(Catalyst.create_matrix(rx),[ 1 0 -6 0; 2 1 -6 -2; 0 2 -12 0;]) - @test isequal(Catalyst.get_stoich(rx),[6, 6, 1, 6]) + @test isequal(Catalyst.get_balanced_stoich(rx),[[6, 6, 1, 6]]) balanced_rx = Reaction(k,[CO2,H2O],[C6H12O6,O2],[6, 6], [1, 6]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end -# @reaction k, H2O --> H2O -let +# @reaction k, H2O --> H2O +let @variables t @species H(t) O(t) @compound H2O(t) 2H O @@ -46,11 +49,13 @@ let rx = Reaction(1.0, [H2O], [H2O], [2], [2]) balanced_rx = Reaction(1.0, [H2O], [H2O], [1], [1]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end -# @reaction k, 2H + 1O --> H2O -let +# @reaction k, 2H + 1O --> H2O +let @variables t @species H(t) O(t) @compound H2O(t) 2H O @@ -58,11 +63,13 @@ let rx = Reaction(1.0, [H, O], [H2O], [23, 1], [7]) balanced_rx = Reaction(1.0, [H,O], [H2O], [2, 1], [1]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 1CH4 + 2O2 --> 1CO2 + 2H2O -let +let @variables t @species H(t) O(t) C(t) @compound CH4(t) C 4H @@ -73,7 +80,9 @@ let rx = Reaction(1.0, [CH4, O2], [CO2, H2O]) balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, N2 + 3H2 --> 2NH3 @@ -87,10 +96,12 @@ let rx = Reaction(1.0, [N2, H2], [NH3]) balanced_rx = Reaction(1.0, [N2, H2], [NH3], [1 ,3], [2]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end -# @reaction k, C2H5OH + CH3COOH --> C4H8O2 + H2O +# @reaction k, C2H5OH + CH3COOH --> C4H8O2 + H2O let @variables t @species C(t) H(t) O(t) @@ -102,13 +113,15 @@ let rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O]) balanced_rx = Reaction(1.0, [C2H5OH, CH3COOH], [C4H8O2, H2O], [1, 1], [1, 1]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end -# @reaction k, 2Ca3PO42 --> 6CaO + 1P4O10 -let +# @reaction k, 2Ca3PO42 --> 6CaO + 1P4O10 +let @variables t - @species Ca(t) P(t) O(t) + @species Ca(t) P(t) O(t) @compound Ca3PO42(t) 3Ca 2P 8O @compound CaO(t) Ca O @compound P4O10(t) 4P 10O @@ -116,21 +129,25 @@ let rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10]) balanced_rx = Reaction(1.0, [Ca3PO42], [CaO, P4O10], [2], [6, 1]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 4Fe + 3O2 + 6H2O --> 4FeOH3 let @variables t @species Fe(t) O(t) H(t) - @compound O2(t) 2O + @compound O2(t) 2O @compound H2O(t) 2H O @compound FeOH3(t) Fe 3H 3O rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3]) balanced_rx = Reaction(1.0, [Fe, O2, H2O], [FeOH3], [4, 3, 6], [4]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 2NaOH + H2SO4 --> Na2SO4 + 2H2O @@ -146,10 +163,12 @@ let rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O]) balanced_rx = Reaction(1.0, [NaOH,H2SO4], [Na2SO4,H2O], [2, 1], [1, 2]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end -# @reaction k, 2NO2 --> 1N2O4 +# @reaction k, 2NO2 --> 1N2O4 let @variables t @species N(t) O(t) @@ -159,11 +178,13 @@ let rx = Reaction(1.0, [NO2], [N2O4]) balanced_rx = Reaction(1.0, [NO2], [N2O4], [2], [1]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 1CaCO3 + 2HCl --> CaCl2 + H2O + CO2 -let +let @variables t @species C(t) H(t) O(t) Ca(t) Cl(t) @compound H2O(t) 2H 1O @@ -174,11 +195,13 @@ let rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O]) balanced_rx = Reaction(1.0,[CaCO3,HCl],[CaCl2,CO2,H2O], [1, 2], [1, 1, 1]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, SiCl4 + 4H2O → H4SiO4 + 4HCl -let +let @variables t @species Si(t) Cl(t) H(t) O(t) @compound SiCl4(t) 1Si 4Cl @@ -188,11 +211,13 @@ let rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl]) balanced_rx = Reaction(1.0,[SiCl4,H2O],[H4SiO4,HCl], [1,4], [1,4]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 2Al + 6HCl → 2AlCl3 + 3H2 -let +let @variables t @species Al(t) Cl(t) H(t) @compound HCl(t) H Cl @@ -201,11 +226,13 @@ let rx = Reaction(1.0,[Al,HCl],[AlCl3,H2]) balanced_rx = Reaction(1.0,[Al,HCl],[AlCl3,H2],[2,6], [2,3]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, Na2CO3 + 2HCl → 2NaCl + H2O + CO2 -let +let @variables t @species Na(t) C(t) O(t) H(t) Cl(t) @compound Na2CO3(t) 2Na C 3O @@ -216,11 +243,13 @@ let rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2]) balanced_rx = Reaction(1.0,[Na2CO3,HCl],[NaCl,H2O,CO2], [1,2], [2,1,1]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 2C7H6O2 + 15O2 → 14CO2 + 6H2O -let +let @variables t @species C(t) H(t) O(t) @compound C7H6O2(t) 7C 6H 2O @@ -230,11 +259,13 @@ let rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O]) balanced_rx = Reaction(1.0,[C7H6O2,O2],[CO2,H2O], [2,15], [14,6]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k,Fe2(SO4)3 + 6KOH → 3K2SO4 + 2Fe(OH)3 -let +let @variables t @species Fe(t) S(t) O(t) H(t) K(t) @compound Fe2S3O12(t) 2Fe 3S 12O @@ -244,11 +275,13 @@ let rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3]) #5x4 matrix balanced_rx = Reaction(1.0,[Fe2S3O12,KOH],[K2SO4,FeO3H3], [1,6], [3,2]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 2Ca3(PO4)2 + 6SiO2 → P4O10 + 6CaSiO3 -let +let @variables t @species Ca(t) P(t) O(t) Si(t) @compound Ca3P2O8(t) 3Ca 2P 8O @@ -258,11 +291,13 @@ let rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3]) #5x4 matrix balanced_rx = Reaction(1.0,[Ca3P2O8,SiO2],[P4O10,CaSiO3], [2,6] , [1,6]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, 4KClO3 → 3KClO4 + KCl -let +let @variables t @species K(t) Cl(t) O(t) @compound KClO3(t) K Cl 3O @@ -270,12 +305,14 @@ let @compound KCl(t) K Cl rx = Reaction(1.0,[KClO3],[KClO4,KCl]) - balanced_rx = Reaction(1.0,[KClO3],[KClO4,KCl], [4], [3,1]) - @test isequal(balanced_rx, balance_reaction(rx)) + balanced_rx = Reaction(1.0,[KClO3],[KClO4,KCl], [4], [3,1]) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, Al2(SO4)3 + 3Ca(OH)2 → 2Al(OH)3 + 3CaSO4 -let +let @variables t @species Al(t) S(t) O(t) Ca(t) O(t) (H) @compound Al2S3O12(t) 2Al 3S 12O @@ -283,13 +320,15 @@ let @compound AlO3H3(t) Al 3O 3H @compound CaSO4(t) Ca S 4O - rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4]) + rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4]) balanced_rx = Reaction(1.0,[Al2S3O12,CaO2H2],[AlO3H3,CaSO4], [1,3], [2,3]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, H2SO4 + 8HI → H2S + 4I2 + 4H2O -let +let @variables t @species H(t) S(t) O(t) I(t) @compound H2SO4(t) 2H S 4O @@ -298,13 +337,15 @@ let @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]) - @test isequal(balanced_rx, balance_reaction(rx)) + rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O]) + balanced_rx = Reaction(1.0,[H2SO4,HI],[H2S,I2,H2O], [1,8], [1,4,4]) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # @reaction k, C2H4 + 3O2 ↔ 2CO2 + 2H2O -let +let @variables t @species C(t) H(t) O(t) @compound C2H4(t) 2C 4H @@ -314,11 +355,13 @@ let rx = Reaction(1.0,[C2H4,O2],[CO2,H2O]) balanced_rx = Reaction(1.0,[C2H4,O2],[CO2,H2O],[1,3],[2,2]) - @test isequal(balanced_rx, balance_reaction(rx)) + brxs = balance_reaction(rx) + @test length(brxs) == 1 + @test isequal(balanced_rx, first(brxs)) end # Infinite solutions -let +let @variables t @species C(t) H(t) O(t) @compound CO(t) C O @@ -328,12 +371,13 @@ let @compound H2O(t) 2H O rx = Reaction(1.0,[CO,CO2,H2],[CH4,H2O]) - @test typeof(balance_reaction(rx)) == Vector{Reaction} - @test length(balance_reaction(rx)) == 2 + brxs = balance_reaction(rx) + @test_logs (:warn, r"Infinite balanced reactions*") match_mode=:any balance_reaction(rx) + @test length(brxs) == 2 end # No way to balance -let +let @variables t @species Fe(t) S(t) O(t) H(t) N(t) @@ -344,5 +388,19 @@ let @compound H2SO4(t) 2H S 4O rx = Reaction(1.0,[FeS2,HNO3],[Fe2S3O12,NO,H2SO4]) - @test_throws ErrorException("Chemical equation cannot be balanced") balance_reaction(rx) + brxs = balance_reaction(rx) + @test_logs (:warn, "Unable to balance reaction.") match_mode=:any balance_reaction(rx) + @test isempty(brxs) end + +# test errors on compounds of compounds +let + @variables t + @species C(t) H(t) O(t) + @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) +end \ No newline at end of file From 4a7d1116cd194b460e7e01f32104bee7950a91fa Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 1 Sep 2023 18:59:48 -0400 Subject: [PATCH 20/37] format --- test/miscellaneous_tests/reaction_balancing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/miscellaneous_tests/reaction_balancing.jl b/test/miscellaneous_tests/reaction_balancing.jl index 65cb0411f7..b858cf7aa9 100644 --- a/test/miscellaneous_tests/reaction_balancing.jl +++ b/test/miscellaneous_tests/reaction_balancing.jl @@ -403,4 +403,4 @@ let rx = Reaction(1.0, [CO, H2], [COH2]) @test_throws Catalyst.COMPOUND_OF_COMPOUND_ERROR balance_reaction(rx) -end \ No newline at end of file +end From 6fc73cf16ddf84da52aaec61bcb026e2fa26c0aa Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 1 Sep 2023 19:42:50 -0400 Subject: [PATCH 21/37] try using enumerate instead of pairs --- 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 3cc3001d28..96b0118f39 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -152,7 +152,7 @@ function get_balanced_stoich(reaction::Reaction) nullity = size(X, 2) stoichvecs = Vector{Vector{Int64}}() - for (j, col) in pairs(eachcol(X)) + for (j, col) in enumerate(eachcol(X)) signs = unique(col) # If there is only one basis vector and the signs are not all the same this means From 88342032931fccc7e793cf2af32d682d9b20a28d Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 1 Sep 2023 21:00:25 -0400 Subject: [PATCH 22/37] update HISTORY for compounds --- HISTORY.md | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ Project.toml | 2 +- 2 files changed, 50 insertions(+), 1 deletion(-) diff --git a/HISTORY.md b/HISTORY.md index 7bc2db909e..b3bc6d0366 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -2,6 +2,55 @@ ## Catalyst unreleased (master branch) +## Catalyst 13.4 +- Added the ability to create species that represent chemical compounds and know + their constituents. For example, water can be created and queried as + ```julia + @variables t + @species H(t) O(t) + @compound H2O(t) 2*H O + iscompound(H2O) == true + isspecies(H2O) == true # compounds are also species, so can be used in reactions + isequal(components(H2O), [H, O]) + coefficients(H2O) == [2, 1] + ``` +- Added reaction balancing via the `balance_reaction` command, which returns a + vector of balanced reaction versions, i.e. + ```julia + @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 + + # unbalanced reaction to balance + rx = Reaction(1.0, [CH4, O2], [CO2, H2O]) + + # calculate a balanced version, this returns a vector + # storing a single balanced version of the reaction in this case + brxs = balance_reaction(rx) + + # what one would calculate by hand + balanced_rx = Reaction(1.0, [CH4, O2], [CO2, H2O], [1, 2], [1, 2]) + + # testing equality + @test isequal(balanced_rx, first(brxs)) + ``` +- Note that balancing works via calculating the nullspace of an associated + integer matrix that stores in entry `(i,j)` a signed integer representing the + number of times the `i`'th atom appears within the `j`th compound. The entry + is positive for a substrate and negative for a product. One cannot balance a + reaction involving compounds of compounds currently. A non-empty solution + vector is returned if the reaction can be balanced in exactly one way with + minimal coefficients while preserving the set of substrates and products, i.e. + if the dimension of the nullspace is one. If the dimension is greater than one + we return a `Reaction` for each nullspace basis vector, but note that they may + currently interchange substrates and products (i.e. we do not solve for if + there is a linear combination of them that preserves the set of substrates and + products). An empty `Reaction` vector indicates it is not possible to balance + the reaction. + ## Catalyst 13.2 - Array parameters, species, and variables can be use in the DSL if explicitly declared with `@parameters`, `@species`, or `@variables` respectively, i.e. diff --git a/Project.toml b/Project.toml index e756a1ae2d..e9e4005608 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Catalyst" uuid = "479239e8-5488-4da2-87a7-35f2df7eef83" -version = "13.3.0" +version = "13.4.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" From 1aeb9dbabe3e2e31ad4ce09159460c0e22640b3a Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 1 Sep 2023 22:04:38 -0400 Subject: [PATCH 23/37] fix hh docs --- docs/Project.toml | 2 +- docs/src/example_networks/hodgkin_huxley_equation.md | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 0b30ebf4c9..bae4e788f6 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -28,7 +28,7 @@ HomotopyContinuation = "2.6" Latexify = "0.15, 0.16" ModelingToolkit = "8.47" NonlinearSolve = "1.6.0" -Optimization = "3.11" +Optimization = "~3.15" OptimizationOptimisers = "0.1.1" OrdinaryDiffEq = "6" Plots = "1.36" diff --git a/docs/src/example_networks/hodgkin_huxley_equation.md b/docs/src/example_networks/hodgkin_huxley_equation.md index 6b0e67fb88..25a4c2e0b0 100644 --- a/docs/src/example_networks/hodgkin_huxley_equation.md +++ b/docs/src/example_networks/hodgkin_huxley_equation.md @@ -14,7 +14,7 @@ We begin by importing some necessary packages. ```@example hh1 using ModelingToolkit, Catalyst, NonlinearSolve using DifferentialEquations, Symbolics -using Plots, IfElse +using Plots ``` We'll build a simple Hodgkin-Huxley model for a single neuron, with the voltage, @@ -32,7 +32,7 @@ The transition rate functions, which depend on the voltage, $V(t)$, are given by ```@example hh1 function αₘ(V) theta = (V + 45) / 10 - IfElse.ifelse(theta == 0.0, 1.0, theta/(1 - exp(-theta))) + ifelse(theta == 0.0, 1.0, theta/(1 - exp(-theta))) end βₘ(V) = 4*exp(-(V + 70)/18) @@ -41,14 +41,14 @@ end function αₙ(V) theta = (V + 60) / 10 - IfElse.ifelse(theta == 0.0, .1, .1*theta / (1 - exp(-theta))) + ifelse(theta == 0.0, .1, .1*theta / (1 - exp(-theta))) end βₙ(V) = .125 * exp(-(V + 70)/80) nothing # hide ``` - We now declare the symbolic variable, `V(t)`, that will represent the - transmembrane potential +We now declare the symbolic variable, `V(t)`, that will represent the +transmembrane potential ```@example hh1 @variables t V(t) From 139b9c5687b681f066238cfa2960b8db208e63ca Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 1 Sep 2023 22:10:08 -0400 Subject: [PATCH 24/37] fix simulation structure tutorial too --- .../simulation_structure_interfacing.md | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/docs/src/catalyst_applications/simulation_structure_interfacing.md b/docs/src/catalyst_applications/simulation_structure_interfacing.md index c2be900c8d..a5ec97571a 100644 --- a/docs/src/catalyst_applications/simulation_structure_interfacing.md +++ b/docs/src/catalyst_applications/simulation_structure_interfacing.md @@ -1,12 +1,12 @@ # [Interfacing problems, integrators, and solutions](@id simulation_structure_interfacing) -When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on a problem, during which the state of the simulation is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access, or modify the state, or parameter, values of problems, integrators, and solutions structures. +When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on a problem, during which the state of the simulation is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access, or modify the state, or parameter, values of problems, integrators, and solutions structures. Generally, when we have a structure `simulation_struct` and want to interface with the state (or parameter) `G`, we use `simulation_struct[:G]` to access the value, and `simulation_struct[:G] = 5.0` to set it to a new value. However, see the following examples for full details. ## Interfacing problem objects We begin by demonstrating how we can interface with problem objects. We will demonstrate using a `ODEProblem`, however, it works similarily for other problem types. -```example ex1 +```@example ex1 using Catalyst rn = @reaction_network begin (k1,k2), X1 <--> X2 @@ -19,23 +19,23 @@ nothing # hide ``` We can find the value of a state simply by interfacing with the corresponding symbol: -```example ex1 +```@example ex1 oprob[:X1] ``` with the notation being identical for parameters: -```example ex1 +```@example ex1 oprob[:k1] ``` If we want to change a state's initial condition value, we use the following notation -```example ex1 +```@example ex1 oprob[:X1] = 10.0 ``` with parameters using the same notation. #### Remaking problems using the `remake` function Typically, when modifying problems, it is recommended to use the `remake` function. Unlike when we do `oprob[:X1] = 10.0` (which modifies the problem in question), `remake` creates a new problem object. The `remake` function takes a problem as input, and any fields you wish to modify (and their new values) as optional inputs. Thus, we can do: -```example ex1 +```@example ex1 using DifferentialEquations @unpack X1, X2, k1, k2 = rn oprob1 = ODEProblem(rn, u0, (0.0,10.0), p) @@ -43,19 +43,19 @@ oprob2 = remake(oprob1; u0=[X1 => 10.0, X2 => 50.0], tspan=(0.0,100.0), p=[k1 => nothing # hide ``` and we can now check the fields of `oprob2` -```example ex1 +```@example ex1 oprob2.u0 ``` -```example ex1 +```@example ex1 oprob2.tspan ``` -```example ex1 +```@example ex1 oprob2.p ``` Please note that, currently, `remake` does not work while giving `Symbol`s as input (e.g `[:X1 => 10.0, :X2 => 50.0]`), but we need to unpack the symbolic variables and use them instead (please see the end of this tutorial for more information on using symbolic variables rather than `Symbol`s). When using `remake`, we only have to provide the fields that we actually wish to change, e.g. -```example ex1 +```@example ex1 oprob3 = remake(oprob1; u0=[X1 => 10.0, X2 => 50.0]) nothing # hide ``` @@ -65,19 +65,19 @@ will only update the initial conditions. ## Interfacing integrator objects During a simulation, the solution is stored in an integrator object, we will here describe how to interface with these. The primary circumstance under which a user may wish to do so is when using [callbacks](@ref advanced_simulations_callbacks). We can create an integrator by calling `init` on our problem ([while circumstances where the user might want to use `init` function exist](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Initialization-and-Stepping), since integrators are automatically created during simulations, these are rare): -```example ex1 +```@example ex1 integrator = init(oprob) ``` Using a similar syntax to problems, we can get the current values of a state within the integrator: -```example ex1 +```@example ex1 integrator[:X1] ``` or a parameter: -```example ex1 +```@example ex1 integrator[:k1] ``` Similarly, we can update their values using: -```example ex1 +```@example ex1 integrator[:X1] = 10.0 ``` Please read [this](@ref advanced_simulations_ssa_callbacks) with regards to updating integrators of `JumpProblem`s. @@ -86,15 +86,15 @@ Please read [this](@ref advanced_simulations_ssa_callbacks) with regards to upda ## Interfacing solution objects Finally, we consider solution objects. First, we simulate our problem: -```example ex1 +```@example ex1 sol = solve(oprob) ``` For solutions, when we access a state, we get its whole simulation vector: -```example ex1 +```@example ex1 sol[:X1] ``` while when we access a parameter we only get a single value: -```example ex1 +```@example ex1 sol[:k1] ``` Finally, we note that we cannot change the values of solution states or parameters (i.e. both `sol[:X1] = 0.0` and `sol[:k1] = 0.0` generate errors). @@ -102,14 +102,14 @@ Finally, we note that we cannot change the values of solution states or paramete ## Interfacing using symbolic representation Catalyst is built on an *intermediary representation* implemented by (ModelingToolkit.jl)[https://github.com/SciML/ModelingToolkit.jl]. ModelingToolkit is a modelling framework where one first declares a set of symbolic variables and parameters using e.g. -```example ex2 +```@example ex2 using ModelingToolkit @parameters σ ρ β @variables t x(t) y(t) z(t) nothing # hide ``` -and then uses these to build systems of equations. Here, these symbolic variables (`x`, `y`, and `z`) and parameters (`σ`, `ρ`, and `β`) can be used to interface a `problem`, `integrator`, and `solution` object (like we did previously, but using Symbols, e.g. `:X`). Since Catalyst models are built on ModelingToolkit, these models also contain similar symbolic variables and parameters. -```example ex2 +and then uses these to build systems of equations. Here, these symbolic variables (`x`, `y`, and `z`) and parameters (`σ`, `ρ`, and `β`) can be used to interface a `problem`, `integrator`, and `solution` object (like we did previously, but using Symbols, e.g. `:X`). Since Catalyst models are built on ModelingToolkit, these models also contain similar symbolic variables and parameters. +```@example ex2 using Catalyst rn = @reaction_network begin (k1,k2), X1 <--> X2 @@ -118,7 +118,7 @@ end @unpack k1,k2,X1,X2 = rn ``` Here, we first list the parameters and variables (for reaction systems the latter are typically species) we wish to import (in this case we select all, but we could select only a subset), next we denote from which model (here `rn`) from which we wish to import from. Next, these values can be used directly to interface with e.g. an `ODEProblem`: -```example ex2 +```@example ex2 u0 = [X1 => 1.0, X2 => 5.0] p = [:k1 => 5.0, :k2 => 2.0] oprob = ODEProblem(rn, u0, (0.0,10.0), p) @@ -128,6 +128,6 @@ oprob[k1] To interface with integrators and solutions we use a similar syntax. Finally, instead of using `@unpack` to access a symbolic variable or parameter, we can access it directly using `rn.X1`, and thus access a state of our `ODEProblem` using -```example ex2 +```@example ex2 oprob[rn.X1] ``` From 363e17107dd7f2ac34bb379a3efa5b0583d850e6 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 2 Sep 2023 02:55:44 -0400 Subject: [PATCH 25/37] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e9e4005608..e1198f1ef7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Catalyst" uuid = "479239e8-5488-4da2-87a7-35f2df7eef83" -version = "13.4.0" +version = "13.4.1" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" From 80ca5f08c7f01e955a86a877c1e3c52d70a5e36a Mon Sep 17 00:00:00 2001 From: Vasily Ilin Date: Sun, 3 Sep 2023 17:34:06 -0700 Subject: [PATCH 26/37] Update introduction_to_catalyst.md Fix two links to `JumpProcesses` documentation. --- docs/src/introduction_to_catalyst/introduction_to_catalyst.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md index 4a289f147a..ca4daae098 100644 --- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md +++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md @@ -194,10 +194,10 @@ We see that oscillations remain, but become much noisier. Note, in constructing the `JumpProblem` we could have used any of the SSAs that are part of JumpProcesses instead of the `Direct` method, see the list of SSAs (i.e., constant rate jump aggregators) in the -[documentation](https://docs.sciml.ai/stable/modules/JumpProcesses/jump_types/#Constant-Rate-Jump-Aggregators). +[documentation](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation). Common questions that arise in using the JumpProcesses SSAs (i.e. Gillespie methods) -are collated in the [JumpProcesses FAQ](https://docs.sciml.ai/stable/modules/JumpProcesses/faq/). +are collated in the [JumpProcesses FAQ](https://docs.sciml.ai/JumpProcesses/stable/faq/). --- ## Chemical Langevin equation (CLE) stochastic differential equation (SDE) models From eb9a7a506622d5163a64007c2694e52daec3796f Mon Sep 17 00:00:00 2001 From: Elisabeth Roesch Date: Thu, 7 Sep 2023 13:14:33 +1200 Subject: [PATCH 27/37] Update index.md add ".jl" to ModelingToolkit in first paragraph --- docs/src/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index 0688e32dea..d182e7a813 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,7 +4,7 @@ Catalyst.jl is a symbolic modeling package for analysis and high performance simulation of chemical reaction networks. Catalyst defines symbolic [`ReactionSystem`](@ref)s, which can be created programmatically or easily specified using Catalyst's domain specific language (DSL). Leveraging -[ModelingToolkit](https://docs.sciml.ai/ModelingToolkit/stable/) and +[ModelingToolkit.jl](https://docs.sciml.ai/ModelingToolkit/stable/) and [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/), Catalyst enables large-scale simulations through auto-vectorization and parallelism. Symbolic `ReactionSystem`s can be used to generate ModelingToolkit-based models, allowing From a6a2ccc7b3e78299beecb91b2acfc6c7cd7f8c0a Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 11 Sep 2023 01:02:14 +0000 Subject: [PATCH 28/37] Bump actions/checkout from 3 to 4 Bumps [actions/checkout](https://github.com/actions/checkout) from 3 to 4. - [Release notes](https://github.com/actions/checkout/releases) - [Changelog](https://github.com/actions/checkout/blob/main/CHANGELOG.md) - [Commits](https://github.com/actions/checkout/compare/v3...v4) --- updated-dependencies: - dependency-name: actions/checkout dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/CI.yml | 2 +- .github/workflows/Documentation.yml | 2 +- .github/workflows/FormatCheck.yml | 2 +- .github/workflows/Invalidations.yml | 4 ++-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 0b1927705d..699698f275 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -17,7 +17,7 @@ jobs: - '1' - '1.6' steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 27e4bcf301..4cdc1cf3ba 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -11,7 +11,7 @@ jobs: build: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@latest with: version: '1' diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 329ae6549c..9cca2b9011 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -21,7 +21,7 @@ jobs: with: version: ${{ matrix.julia-version }} - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Install JuliaFormatter and format # This will use the latest version by default but you can set the version like so: # diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml index 4d0004e831..28b9ce2fad 100644 --- a/.github/workflows/Invalidations.yml +++ b/.github/workflows/Invalidations.yml @@ -19,12 +19,12 @@ jobs: - uses: julia-actions/setup-julia@v1 with: version: '1' - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-invalidations@v1 id: invs_pr - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: ref: ${{ github.event.repository.default_branch }} - uses: julia-actions/julia-buildpkg@v1 From f8332aea1c8efcbf056d45fe5835a124c87cd378 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 19:31:42 -0400 Subject: [PATCH 29/37] Initiate --- .../catalyst_functionality/dsl_description.md | 20 +++++++++++++++++++ .../reactionsystem.jl | 11 ++++++++++ 2 files changed, 31 insertions(+) diff --git a/docs/src/catalyst_functionality/dsl_description.md b/docs/src/catalyst_functionality/dsl_description.md index 7e2a140113..12af0e9db2 100644 --- a/docs/src/catalyst_functionality/dsl_description.md +++ b/docs/src/catalyst_functionality/dsl_description.md @@ -412,6 +412,26 @@ sol = solve(oprob) plot(sol) ``` +## Setting initial conditions that depend on parameters +It is possible to set the initial condition of one (or several) species so that they depend on some system parameter. This is done in a similar way as default initial conditions, but giving the parameter instead of a value. When doing this, we also need to ensure that the initial condition parameter is a variable of the system: +```@example tut2 +rn = @reaction_network begin + @species X(t)=X0 + @parameters X0 + p, 0 --> X + d, X --> ∅ +end +``` +We can now simulate the network without providing any initial conditions: +```@example tut2 +u0 = [] +tspan = (0.0, 10.0) +p = [:p => 2.0, :d => .1, :X0 => 1.0] +oprob = ODEProblem(rn, u0, tspan, p) +sol = solve(oprob) +plot(sol) +``` + ## Naming the generated `ReactionSystem` ModelingToolkit uses system names to allow for compositional and hierarchical models. To specify a name for the generated `ReactionSystem` via the diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_structure/reactionsystem.jl index 21e95f54bd..8fbc5f8f07 100644 --- a/test/reactionsystem_structure/reactionsystem.jl +++ b/test/reactionsystem_structure/reactionsystem.jl @@ -631,6 +631,17 @@ let @test issetequal(parameters(rs), [k, b]) end +# Test parameteric initial conditions. +let + @parameters d X0 + @variables t + @species X(t)=X0 + rx = Reaction(d, [X], nothing, [1], nothing) + @named rs = ReactionSystem([rx], t) + prob = ODEProblem(rs, [], (0.0, 1.0), [d => 1.0, X0 => 7.6]) + @test prob[X] == 7.6 +end + # Test balanced_bc_check. let @variables t From ef8e7e03d7aa4c8c91daddb9e1b955ff87ba3bab Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 19:56:44 -0400 Subject: [PATCH 30/37] Update docs/src/catalyst_functionality/dsl_description.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_functionality/dsl_description.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_functionality/dsl_description.md b/docs/src/catalyst_functionality/dsl_description.md index 12af0e9db2..7a52176f11 100644 --- a/docs/src/catalyst_functionality/dsl_description.md +++ b/docs/src/catalyst_functionality/dsl_description.md @@ -416,8 +416,8 @@ plot(sol) It is possible to set the initial condition of one (or several) species so that they depend on some system parameter. This is done in a similar way as default initial conditions, but giving the parameter instead of a value. When doing this, we also need to ensure that the initial condition parameter is a variable of the system: ```@example tut2 rn = @reaction_network begin - @species X(t)=X0 @parameters X0 + @species X(t)=X0 p, 0 --> X d, X --> ∅ end From a9cbe5eebe4c01fe2ab8ac326430bbdcbcac7fc7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 08:06:14 -0400 Subject: [PATCH 31/37] add missing stable rng usage --- test/model_simulation/simulate_PDEs.jl | 6 +++- .../reactionsystem.jl | 28 +++++++++++-------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/test/model_simulation/simulate_PDEs.jl b/test/model_simulation/simulate_PDEs.jl index b31e01db3e..ed57fb12b1 100644 --- a/test/model_simulation/simulate_PDEs.jl +++ b/test/model_simulation/simulate_PDEs.jl @@ -5,9 +5,13 @@ using Catalyst, OrdinaryDiffEq, Test using ModelingToolkit, DomainSets const MT = ModelingToolkit +# Sets rnd number. +using StableRNGs +rng = StableRNG(12345) + # Test function. function icfun(n, x, y, A) - float(rand(Poisson(round(n * A * 10))) / A / 10) + float(rand(rng, Poisson(round(n * A * 10))) / A / 10) end ### Run Tests ### diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_structure/reactionsystem.jl index 8fbc5f8f07..0180fc1ebc 100644 --- a/test/reactionsystem_structure/reactionsystem.jl +++ b/test/reactionsystem_structure/reactionsystem.jl @@ -4,6 +4,10 @@ using Catalyst, LinearAlgebra, JumpProcesses, Test, OrdinaryDiffEq, StochasticDiffEq const MT = ModelingToolkit +# Sets rnd number. +using StableRNGs +rng = StableRNG(12345) + # Create test network. @parameters k[1:20] @variables t @@ -140,8 +144,8 @@ end # Don't ask me (Torkel) why the statement before/after is needed. t = 0.0 let - p = rand(length(k)) - u = rand(length(k)) + p = rand(rng, length(k)) + u = rand(rng, length(k)) du = oderhs(u, p, t) G = sdenoise(u, p, t) sdesys = convert(SDESystem, rs) @@ -162,8 +166,8 @@ end # Tests the noise_scaling argument. let - p = rand(length(k) + 1) - u = rand(length(k)) + p = rand(rng, length(k) + 1) + u = rand(rng, length(k)) t = 0.0 G = p[21] * sdenoise(u, p, t) @variables η @@ -176,8 +180,8 @@ end # Tests the noise_scaling vector argument. let - p = rand(length(k) + 3) - u = rand(length(k)) + p = rand(rng, length(k) + 3) + u = rand(rng, length(k)) t = 0.0 G = vcat(fill(p[21], 8), fill(p[22], 3), fill(p[23], 9))' .* sdenoise(u, p, t) @variables η[1:3] @@ -192,8 +196,8 @@ end # Tests using previous parameter for noise scaling let - p = rand(length(k)) - u = rand(length(k)) + p = rand(rng, length(k)) + u = rand(rng, length(k)) t = 0.0 G = [p p p p]' .* sdenoise(u, p, t) sdesys_noise_scaling = convert(SDESystem, rs; noise_scaling = k) @@ -205,7 +209,7 @@ end # Test with JumpSystem. let - p = rand(length(k)) + p = rand(rng, length(k)) @variables t @species A(t) B(t) C(t) D(t) E(t) F(t) rxs = [Reaction(k[1], nothing, [A]), # 0 -> A @@ -239,9 +243,9 @@ let @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.ConstantRateJump, cidxs)) @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.VariableRateJump, vidxs)) - pars = rand(length(k)) - u0 = rand(2:10, 6) - ttt = rand() + pars = rand(rng, length(k)) + u0 = rand(rng, 2:10, 6) + ttt = rand(rng) jumps = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}(undef, length(rxs)) From 7effc66e914b4da3d26c4b4093da3ceb72e59e17 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 08:08:23 -0400 Subject: [PATCH 32/37] save 10s (remove repeat in higher order reactions) --- test/reactionsystem_structure/higher_order_reactions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/reactionsystem_structure/higher_order_reactions.jl b/test/reactionsystem_structure/higher_order_reactions.jl index ea38027f42..1930cb8695 100644 --- a/test/reactionsystem_structure/higher_order_reactions.jl +++ b/test/reactionsystem_structure/higher_order_reactions.jl @@ -63,7 +63,7 @@ let d * binomial(X10, 2), 2X10 ⟾ ∅ end - for factor in [1e-1, 1e0], repeat in 1:5 + for factor in [1e-1, 1e0] u0 = rand(rng, 1:Int64(factor * 100), length(get_states(higher_order_network_1))) p = factor * rand(rng, length(get_ps(higher_order_network_3))) prob1 = JumpProblem(higher_order_network_1, From 3cb02e3625da947060bdad1777ac4fe70ea6fa5a Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 08:09:54 -0400 Subject: [PATCH 33/37] remove now superfluous print statements --- test/model_simulation/simulate_ODEs.jl | 6 ++---- test/model_simulation/simulate_jumps.jl | 2 -- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/test/model_simulation/simulate_ODEs.jl b/test/model_simulation/simulate_ODEs.jl index dc24463452..4820b935f1 100644 --- a/test/model_simulation/simulate_ODEs.jl +++ b/test/model_simulation/simulate_ODEs.jl @@ -9,7 +9,7 @@ using StableRNGs rng = StableRNG(12345) # Fetch test networks. -include("../test_networks.jl") +@time include("../test_networks.jl") ### Compares to Known Solution ### @@ -119,8 +119,6 @@ end let for (i, network) in enumerate(reaction_networks_all) - (i % 5 == 0) && - println("Iteration " * string(i) * " at line 104 in file solve_ODEs.jl") for factor in [1e-1, 1e0, 1e1] u0 = factor * rand(rng, length(get_states(network))) p = factor * rand(rng, length(get_ps(network))) @@ -157,7 +155,7 @@ let p = [:k => 1.0] oprob = ODEProblem(rn, u0, tspan, p; combinatoric_ratelaws = false) du1 = du2 = zeros(2) - u = rand(2) + u = rand(rng, 2) oprob.f(du1, u, [1.0], 0.0) oderhs(du2, u, [1.0], 0.0) @test isapprox(du1, du2, rtol = 1e3 * eps()) diff --git a/test/model_simulation/simulate_jumps.jl b/test/model_simulation/simulate_jumps.jl index c1dccbb797..b10727cc4b 100644 --- a/test/model_simulation/simulate_jumps.jl +++ b/test/model_simulation/simulate_jumps.jl @@ -119,8 +119,6 @@ end let for (i, network) in enumerate(reaction_networks_all) - (i % 5 == 0) && - println("Iteration " * string(i) * " at line 102 in file solve_jumps.jl") for factor in [1e-1, 1e0, 1e1] u0 = rand(rng, 1:Int64(factor * 100), length(get_states(network))) p = factor * rand(rng, length(get_ps(network))) From f86d9f0ec17a625259d6570163780bfcdaf69e13 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 08:13:39 -0400 Subject: [PATCH 34/37] remove unnecessary using ModelingToolkit --- test/dsl/dsl_basics.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/dsl/dsl_basics.jl b/test/dsl/dsl_basics.jl index bb8a00eb8f..78fadfcacc 100644 --- a/test/dsl/dsl_basics.jl +++ b/test/dsl/dsl_basics.jl @@ -2,7 +2,7 @@ ### Fetch Packages and Set Global Variables ### -using Catalyst, ModelingToolkit +using Catalyst @variables t ### Naming Tests ### From e17e4253696ad36c51cd68f0db5b55a6a26ac533 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 08:35:15 -0400 Subject: [PATCH 35/37] remove using ModelingToolkit and rename function overwrite. --- .../symbolic_stoichiometry.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/miscellaneous_tests/symbolic_stoichiometry.jl b/test/miscellaneous_tests/symbolic_stoichiometry.jl index 142b21cdac..d817c40654 100644 --- a/test/miscellaneous_tests/symbolic_stoichiometry.jl +++ b/test/miscellaneous_tests/symbolic_stoichiometry.jl @@ -1,4 +1,4 @@ -using Catalyst, ModelingToolkit, OrdinaryDiffEq, Test, LinearAlgebra, JumpProcesses +using Catalyst, OrdinaryDiffEq, Test, LinearAlgebra, JumpProcesses ### Base Tests ### @@ -34,7 +34,7 @@ let du1 = zeros(size(oprob.u0)) oprob.f(du1, oprob.u0, oprob.p, 1.5) - function oderhs(du, u, p, t) + function oderhs1(du, u, p, t) k = p[1] α = p[2] A = u[1] @@ -51,7 +51,7 @@ let end u0 = [uv[2] for uv in u0map] p = Tuple(pv[2] for pv in pmap) - oprob2 = ODEProblem(oderhs, u0, tspan, p) + oprob2 = ODEProblem(oderhs1, u0, tspan, p) du2 = copy(du1) oprob2.f(du2, oprob2.u0, oprob2.p, 1.5) @test norm(du1 .- du2) < 100 * eps() @@ -59,7 +59,7 @@ let # Test without rate law scalings. osys = convert(ODESystem, rs, combinatoric_ratelaws = false) oprob = ODEProblem(osys, u0map, tspan, pmap) - function oderhs(du, u, p, t) + function oderhs2(du, u, p, t) k = p[1] α = p[2] A = u[1] @@ -74,7 +74,7 @@ let du[3] = k * rl2 du[4] = α * rl2 end - oprob2 = ODEProblem(oderhs, [uv[2] for uv in u0map], tspan, oprob.p) + oprob2 = ODEProblem(oderhs2, [uv[2] for uv in u0map], tspan, oprob.p) du1 .= 0 du2 .= 0 oprob.f(du1, oprob.u0, oprob.p, 1.5) @@ -85,7 +85,7 @@ let ssys = convert(SDESystem, rs) sf = SDEFunction{false}(ssys, states(ssys), parameters(ssys)) G = sf.g(u0, p, 1.0) - function sdenoise(u, p, t) + function sdenoise1(u, p, t) k = p[1] α = p[2] A = u[1] @@ -100,14 +100,14 @@ let 0.0 k*rl2; 0.0 α*rl2] end - G2 = sdenoise(u0, p, 1.0) + G2 = sdenoise1(u0, p, 1.0) @test norm(G - G2) < 100 * eps() # SDESystem test with no combinatoric rate laws. ssys = convert(SDESystem, rs, combinatoric_ratelaws = false) sf = SDEFunction{false}(ssys, states(ssys), parameters(ssys)) G = sf.g(u0, p, 1.0) - function sdenoise(u, p, t) + function sdenoise2(u, p, t) k = p[1] α = p[2] A = u[1] @@ -122,7 +122,7 @@ let 0.0 k*rl2; 0.0 α*rl2] end - G2 = sdenoise(u0, p, 1.0) + G2 = sdenoise2(u0, p, 1.0) @test norm(G - G2) < 100 * eps() # JumpSystem test. From 5fdbf75176480ffbc4eb81a20a1de5205c489efd Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 08:49:35 -0400 Subject: [PATCH 36/37] Gain 1 minute (cut unnecessary repeats) --- test/model_simulation/simulate_jumps.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/model_simulation/simulate_jumps.jl b/test/model_simulation/simulate_jumps.jl index b10727cc4b..a680ef44cb 100644 --- a/test/model_simulation/simulate_jumps.jl +++ b/test/model_simulation/simulate_jumps.jl @@ -93,8 +93,8 @@ let jumps_3 = (jump_3_1, jump_3_2, jump_3_3, jump_3_4, jump_3_5, jump_3_6) push!(identical_networks, reaction_networks_constraint[5] => jumps_3) - for (i, networks) in enumerate(identical_networks) - for factor in [1e-2, 1e-1, 1e0, 1e1], repeat in 1:3 + for networks in identical_networks + for factor in [1e-2, 1e-1, 1e0, 1e1] (i == 3) && (factor > 1e-1) && continue # Large numbers seems to crash it. u0 = rand(rng, 1:Int64(factor * 100), length(get_states(networks[1]))) p = factor * rand(rng, length(get_ps(networks[1]))) @@ -119,7 +119,7 @@ end let for (i, network) in enumerate(reaction_networks_all) - for factor in [1e-1, 1e0, 1e1] + for factor in [1e0] u0 = rand(rng, 1:Int64(factor * 100), length(get_states(network))) p = factor * rand(rng, length(get_ps(network))) prob = JumpProblem(network, DiscreteProblem(network, u0, (0.0, 1.0), p), @@ -134,7 +134,7 @@ end # No parameter test. let no_param_network = @reaction_network begin (1.2, 5), X1 ↔ X2 end - for factor in [1e1, 1e2] + for factor in [1e1] u0 = rand(rng, 1:Int64(factor * 100), length(get_states(no_param_network))) prob = JumpProblem(no_param_network, DiscreteProblem(no_param_network, u0, (0.0, 1000.0)), Direct()) From d695aab6e0e7fbe6ee087bbef5f3cad275c469b5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 19 Sep 2023 09:21:37 -0400 Subject: [PATCH 37/37] fix remove unnecessary enumerate statement --- test/model_simulation/simulate_jumps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/model_simulation/simulate_jumps.jl b/test/model_simulation/simulate_jumps.jl index a680ef44cb..07f8ccbb98 100644 --- a/test/model_simulation/simulate_jumps.jl +++ b/test/model_simulation/simulate_jumps.jl @@ -93,7 +93,7 @@ let jumps_3 = (jump_3_1, jump_3_2, jump_3_3, jump_3_4, jump_3_5, jump_3_6) push!(identical_networks, reaction_networks_constraint[5] => jumps_3) - for networks in identical_networks + for (i, networks) in enumerate(identical_networks) for factor in [1e-2, 1e-1, 1e0, 1e1] (i == 3) && (factor > 1e-1) && continue # Large numbers seems to crash it. u0 = rand(rng, 1:Int64(factor * 100), length(get_states(networks[1]))) @@ -118,7 +118,7 @@ end ### Checks Simulations Don't Error ### let - for (i, network) in enumerate(reaction_networks_all) + for network in reaction_networks_all for factor in [1e0] u0 = rand(rng, 1:Int64(factor * 100), length(get_states(network))) p = factor * rand(rng, length(get_ps(network)))