Skip to content

Commit

Permalink
Fix commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@LalitChauhan56
LalitChauhan56 committed Jul 27, 2023
1 parent 6324db1 commit b369300
Show file tree
Hide file tree
Showing 2 changed files with 73 additions and 150 deletions.
3 changes: 2 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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"]
220 changes: 71 additions & 149 deletions src/compound.jl
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
const MT = ModelingToolkit
struct CompoundSpecies end
struct CompoundComponents end
struct CompoundCoefficients end
Expand Down Expand Up @@ -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

0 comments on commit b369300

Please sign in to comment.