Merge branch 'master' of https://github.com/LalitChauhan56/Catalyst.jl

src/compound.jl

@@ -149,6 +149,80 @@ function get_stoich(reaction::Reaction)
 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
+
+### 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(size(A,1))
+    b[end] = 1
+
+    # Compute x= A^-1 *b
+    x = inv(A)*b
+
+    # Normalize the stoichiometric coefficients to be integers.
+    smallest_value = minimum(x[x .> 0])  
+    x = x / smallest_value
+
+    return round.(Int, x)   # rounding is due to potential minimal numerical inaccuracies
+end
+
+
 function balance(reaction::Reaction)
     # Calculate the stoichiometric coefficients for the balanced reaction.
     stoichiometries = get_stoich(reaction)
@@ -177,4 +251,3 @@ end
 #         end
 #     end
 #     ex
-# end

test/miscellaneous_tests/compound_macro.jl

@@ -51,7 +51,7 @@ end
 let
     @variables t
     @species C(t) H(t) O(t)
-    @compound OH(t) O H
+    @compound OH(t) 1O 1H
     @compound C3H5OH3(t) 3C 5H 3OH
 
     @test !iscompound(O)