Merge pull request #494 from isaacsas/dsl_updates

Parametric stoich via DSL

Author: isaacsas

src/reaction_network.jl

@@ -192,7 +192,7 @@ end
 #Structure containing information about one reactant in one reaction.
 struct ReactantStruct
     reactant::Union{Symbol,Expr}
-    stoichiometry::Number
+    stoichiometry::ExprValues
 end
 #Structure containing information about one Reaction. Contain all its substrates and products as well as its rate. Contains a specialized constructor.
 struct ReactionStruct
@@ -261,7 +261,7 @@ function make_reaction_system(ex::Expr, parameters; name=:(gensym(:ReactionSyste
     # parse DSL lines
     reactions = get_reactions(ex)
     reactants = get_reactants(reactions)
-    allspecies = union(reactants, get_rate_species(reactions,parameters))
+    allspecies = union(reactants, get_rx_species(reactions,parameters))
     !isempty(intersect(forbidden_symbols,union(allspecies,parameters))) &&
         error("The following symbol(s) are used as species or parameters: "*((map(s -> "'"*string(s)*"', ",intersect(forbidden_symbols,union(species,parameters)))...))*"this is not permited.")
 
@@ -290,7 +290,7 @@ function make_reaction(ex::Expr)
     # parse DSL lines
     reaction   = get_reaction(ex)
     allspecies = get_reactants(reaction)   # species defined by stoich
-    parameters = get_rate_species([reaction],Symbol[])  # anything in a rate is a parameter
+    parameters = get_rx_species([reaction],Symbol[])  # anything in a rate is a parameter
     !isempty(intersect(forbidden_symbols,union(allspecies,parameters))) &&
         error("The following symbol(s) are used as species or parameters: "*((map(s -> "'"*string(s)*"', ",intersect(forbidden_symbols,union(species,parameters)))...))*"this is not permited.")
 
@@ -305,11 +305,17 @@ function make_reaction(ex::Expr)
     end
 end
 
-function get_rate_species(rxs, ps)
+function get_rx_species(rxs, ps)
     pset = Set(ps)
     species_set = Set{Symbol}()
     for rx in rxs
         find_species_in_rate!(species_set, rx.rate, pset)
+        for sub in rx.substrates
+            find_species_in_rate!(species_set, sub.stoichiometry, pset)
+        end
+        for prod in rx.products
+            find_species_in_rate!(species_set, prod.stoichiometry, pset)
+        end
     end
     collect(species_set)
 end
@@ -387,18 +393,31 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues
     end
 end
 
+function processmult(op, mult, stoich)
+    if (mult isa Number) && (stoich isa Number)
+        op(mult, stoich)
+    else
+        :($op($mult,$stoich))
+    end
+end
+
 #Recursive function that loops through the reaction line and finds the reactants and their stoichiometry. Recursion makes it able to handle weird cases like 2(X+Y+3(Z+XY)).
-function recursive_find_reactants!(ex::ExprValues, mult::Number, reactants::Vector{ReactantStruct})
+function recursive_find_reactants!(ex::ExprValues, mult::ExprValues, reactants::Vector{ReactantStruct})
     if typeof(ex)!=Expr || (ex.head == :escape)
         (ex == 0 || in(ex,empty_set)) && (return reactants)
-        if in(ex, getfield.(reactants,:reactant))
+        if any(ex==reactant.reactant for reactant in reactants)
             idx = findall(x -> x==ex, getfield.(reactants,:reactant))[1]
-            reactants[idx] = ReactantStruct(ex,mult+reactants[idx].stoichiometry)
+            reactants[idx] = ReactantStruct(ex,processmult(+,mult,reactants[idx].stoichiometry))
         else
             push!(reactants, ReactantStruct(ex,mult))
         end
     elseif ex.args[1] == :*
-        recursive_find_reactants!(ex.args[3],mult*ex.args[2],reactants)
+        if length(ex.args) == 3
+            recursive_find_reactants!(ex.args[3],processmult(*,mult,ex.args[2]),reactants)
+        else
+            newmult = processmult(*, mult, Expr(:call,ex.args[1:end-1]...))
+            recursive_find_reactants!(ex.args[end],newmult,reactants)
+        end
     elseif ex.args[1] == :+
         for i = 2:length(ex.args)
             recursive_find_reactants!(ex.args[i],mult,reactants)
@@ -409,6 +428,7 @@ function recursive_find_reactants!(ex::ExprValues, mult::Number, reactants::Vect
     return reactants
 end
 
+
 function get_reactants(reaction::ReactionStruct, reactants=Vector{Union{Symbol,Expr}}())
     for reactant in Iterators.flatten((reaction.substrates,reaction.products))
         !in(reactant.reactant,reactants) && push!(reactants,reactant.reactant)

test/symbolic_stoich.jl

@@ -1,6 +1,6 @@
 using Catalyst, ModelingToolkit, OrdinaryDiffEq, Test, LinearAlgebra, DiffEqJump
 
-@parameters t k α
+@parameters k α
 @variables t, A(t), B(t), C(t), D(t)
 rxs = [Reaction(t*k, [A], [B], [2*α^2], [k+α*C])
        Reaction(1.0, [A,B], [C,D], [α,2], [k,α])
@@ -10,6 +10,19 @@ rxs = [Reaction(t*k, [A], [B], [2*α^2], [k+α*C])
 @test issetequal(parameters(rs), [k,α])
 osys = convert(ODESystem, rs)
 
+g = (k+α*C)
+rs2 = @reaction_network rs begin
+    t*k, 2*α^2*A --> $g*B    
+    1.0, α*A + 2*B --> k*C + α*D
+end k α
+@test rs2 == rs
+
+
+rxs2 = [(@reaction t*k, 2*α^2*A --> $g*B),
+        (@reaction 1.0, α*A + 2*B --> k*C + α*D)]
+rs3  = ReactionSystem(rxs2, t; name=:rs)
+@test rs3 == rs
+
 u0map = [A => 3.0, B => 2.0, C => 3.0, D => 1.5]
 pmap  = (k => 2.5, α => 2)
 tspan = (0.0,5.0)
@@ -37,6 +50,8 @@ du2 = copy(du1)
 oprob2.f(du2,oprob2.u0,oprob2.p,1.5)
 @test norm(du1 .- du2) < 100*eps()
 
+
+
 # test without rate law scalings
 osys = convert(ODESystem, rs, combinatoric_ratelaws=false)
 oprob = ODEProblem(osys, u0map, tspan, pmap)