From b369300158ea33bcb5ffff80de2a61dc5d18d735 Mon Sep 17 00:00:00 2001 From: Lalit Chauhan Date: Thu, 27 Jul 2023 23:16:19 +0530 Subject: [PATCH] Fix commit --- Project.toml | 3 +- src/compound.jl | 220 ++++++++++++++++-------------------------------- 2 files changed, 73 insertions(+), 150 deletions(-) diff --git a/Project.toml b/Project.toml index e38955b5f0..3c9f94dfd1 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5" LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" @@ -56,4 +57,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["DomainSets", "Graphviz_jll", "LinearAlgebra", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"] +test = ["DomainSets", "Graphviz_jll", "NonlinearSolve", "OrdinaryDiffEq", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"] diff --git a/src/compound.jl b/src/compound.jl index 3bbfa953c7..b6711643a1 100644 --- a/src/compound.jl +++ b/src/compound.jl @@ -1,4 +1,3 @@ -const MT = ModelingToolkit struct CompoundSpecies end struct CompoundComponents end struct CompoundCoefficients end @@ -91,168 +90,91 @@ function component_coefficients(s) return [c => co for (c, co) in zip(components(s), coefficients(s))] end -# function esc_dollars!(ex) -# if ex isa Expr -# if ex.head == :call && ex.args[1] == :* && ex.args[3] isa Expr && ex.args[3].head == :$ -# coeff = ex.args[2] -# value = eval(ex.args[3].args[1]) -# return coeff * value -# else -# for i in 1:length(ex.args) -# ex.args[i] = esc_dollars!(ex.args[i]) -# end -# end -# end -# ex -# end - -# macro compound(species_expr, arr_expr...) -# # Ensure the species name is a valid expression -# if !(species_expr isa Expr && species_expr.head == :call) -# error("Invalid species name in @compound macro") -# end +### 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 -# # Parse the species name to extract the species name and argument -# species_name = species_expr.args[1] -# species_arg = species_expr.args[2] + 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 -# # Construct the expressions that define the species -# species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0), -# Expr(:call, species_name, species_arg)) + # 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) -# # Construct the expression to set the iscompound metadata -# setmetadata_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), -# Catalyst.CompoundSpecies, -# true)) + return A +end -# # Ensure the expressions are evaluated in the correct scope by escaping them -# escaped_species_expr = esc(species_expr) -# escaped_setmetadata_expr = esc(setmetadata_expr) +function get_stoich(reaction::Reaction) + # Create the matrix A using create_matrix function. + A = create_matrix(reaction) -# # Construct the array from the remaining arguments -# coeffs = [] -# species = [] - -# for expr in arr_expr -# if isa(expr, Expr) && expr.head == :call && expr.args[1] == :* -# coeff = expr.args[2] -# specie = expr.args[3] -# # check if the coefficient is a variable and get its value -# if isa(coeff, Expr) && coeff.head == :$ -# coeff = getfield(Main, coeff.args[1]) -# end -# push!(coeffs, coeff) -# push!(species, specie) -# else -# push!(coeffs, 1) -# push!(species, expr) -# end -# end - -# coeffs_expr = Expr(:vect, coeffs...) -# species_expr = Expr(:vect, species...) - -# # Construct the expression to set the components metadata -# setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), -# 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)) - -# 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) -# end - + # 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 -# macro compound(species_expr, arr_expr...) -# # Ensure the species name is a valid expression -# if !(species_expr isa Expr && species_expr.head == :call) -# error("Invalid species name in @compound macro") -# end + # Normalize the stoichiometric coefficients to be integers. + smallest_value = minimum(x[x .> 0]) + x = x / smallest_value -# # Parse the species name to extract the species name and argument -# species_name = species_expr.args[1] -# species_arg = species_expr.args[2] + return round.(Int, x) # rounding is due to potential minimal numerical inaccuracies +end -# # Construct the expressions that define the species -# species_expr = Expr(:macrocall, Symbol("@species"), LineNumberNode(0), -# Expr(:call, species_name, species_arg)) -# # Construct the expression to set the iscompound metadata -# setmetadata_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), -# Catalyst.CompoundSpecies, -# true)) +function balance(reaction::Reaction) + # Calculate the stoichiometric coefficients for the balanced reaction. + stoichiometries = get_stoich(reaction) -# # Ensure the expressions are evaluated in the correct scope by escaping them -# escaped_species_expr = esc(species_expr) -# escaped_setmetadata_expr = esc(setmetadata_expr) + # Divide the stoichiometry vector into substrate and product stoichiometries. + substoich = stoichiometries[1:length(reaction.substrates)] + prodstoich = stoichiometries[(length(reaction.substrates)) + 1:end] -# # Construct the array from the remaining arguments -# arr = Expr(:vect, (arr_expr)...) + # Create a new reaction with the balanced stoichiometries + balanced_reaction = Reaction(reaction.rate, reaction.substrates, reaction.products, substoich, prodstoich) -# coeffs = [] -# species = [] + # Return the balanced reaction + return balanced_reaction +end -# # Add process_expr! function within the macro -# function process_expr!(ex) -# if ex isa Expr && ex.head == :call && ex.args[1] == :* -# # Extract the coefficient and the species +# function esc_dollars!(ex) +# if ex isa Expr +# if ex.head == :call && ex.args[1] == :* && ex.args[3] isa Expr && ex.args[3].head == :$ # coeff = ex.args[2] -# species = ex.args[3] - -# # Check if the species is an interpolated variable -# if species isa Expr && species.head == :$ -# # Replace the interpolated variable with its value -# species = getfield(Main, species.args[1]) - -# # Multiply the coefficient with the species -# return Expr(:call, :*, coeff, species) -# end -# end -# return ex -# end - -# for expr in arr_expr -# # Call process_expr! on the expression to process it -# processed_expr = process_expr!(expr) - -# if isa(processed_expr, Expr) && processed_expr.head == :call && processed_expr.args[1] == :* -# push!(coeffs, processed_expr.args[2]) -# push!(species, processed_expr.args[3]) +# value = eval(ex.args[3].args[1]) +# return coeff * value # else -# push!(coeffs, 1) -# push!(species, processed_expr) +# for i in 1:length(ex.args) +# ex.args[i] = esc_dollars!(ex.args[i]) +# end # end # end - -# coeffs_expr = Expr(:vect, coeffs...) -# species_expr = Expr(:vect, species...) - -# # Construct the expression to set the components metadata -# setcomponents_expr = :($(species_name) = ModelingToolkit.setmetadata($(species_name), -# 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)) - -# 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) -# end +# ex +# end \ No newline at end of file