Merge pull request #127 from SciML/pl/speciesreference_duplication

fix handling g of duplicate speciesReference

Author: paulflang

Project.toml

@@ -1,7 +1,7 @@
 name = "SBMLToolkit"
 uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 authors = ["paulflang", "anandijain"]
-version = "0.1.22"
+version = "0.1.23"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"

src/reactions.jl

@@ -80,8 +80,8 @@ Get reagents
 function get_reagents(reactant_references::Vector{SBML.SpeciesReference},
                       product_references::Vector{SBML.SpeciesReference},
                       model::SBML.Model)
-    reactants = Num[]
-    products = Num[]
+    reactants = String[]
+    products = String[]
     rstoich = Float64[]
     pstoich = Float64[]
 
@@ -93,11 +93,21 @@ function get_reagents(reactant_references::Vector{SBML.SpeciesReference},
             stoich = 1.0
         end
         iszero(stoich) && @error("Stoichiometry of $sn must be non-zero")
-        push!(reactants, create_var(sn, IV))
-        push!(rstoich, stoich)
+        if sn in reactants
+            idx = findfirst(isequal(sn), reactants)
+            rstoich[idx] += stoich
+        else
+            push!(reactants, sn)
+            push!(rstoich, stoich)
+        end
         if model.species[sn].boundary_condition == true
-            push!(products, create_var(sn, IV))
-            push!(pstoich, stoich)
+            if sn in products
+                idx = findfirst(isequal(sn), products)
+                pstoich[idx] += stoich
+            else
+                push!(products, sn)
+                push!(pstoich, stoich)
+            end
         end
     end
     for pr in product_references
@@ -109,10 +119,17 @@ function get_reagents(reactant_references::Vector{SBML.SpeciesReference},
         end
         iszero(stoich) && @error("Stoichiometry of $sn must be non-zero")
         if model.species[sn].boundary_condition != true
-            push!(products, create_var(sn, IV))
-            push!(pstoich, stoich)
+            if sn in products
+                idx = findfirst(isequal(sn), products)
+                pstoich[idx] += stoich
+            else
+                push!(products, sn)
+                push!(pstoich, stoich)
+            end
         end
     end
+    reactants = map(x -> Num(create_var(x, IV)), reactants)
+    products = map(x -> Num(create_var(x, IV)), products)
 
     if (length(reactants) == 0)
         reactants = nothing

test/reactions.jl

@@ -105,6 +105,13 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r))
               SBMLToolkit.get_reagents(r.reactants, r.products, m))
 @test isequal((nothing, nothing, nothing, nothing),
               SBMLToolkit.get_reagents(r.products, r.reactants, m))
+              
+r = SBML.Reaction(reactants = [SBML.SpeciesReference(species = "s", stoichiometry = 1.0),
+                               SBML.SpeciesReference(species = "s", stoichiometry = 1.0)],
+                  reversible = false)
+m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r))
+@test isequal(([var], [var], [2.0], [2.0]),
+              SBMLToolkit.get_reagents(r.reactants, r.products, m))
 
 # Test use_rate
 @test isequal(SBMLToolkit.use_rate(k1 * s1, [s1], [1]), (k1, false))  # Case hOSU=true
@@ -113,7 +120,7 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r))
 
 # Test get_massaction
 @test isequal(SBMLToolkit.get_massaction(k1 * s1, [s1], [1]), k1)  # Case hOSU=true
-@test_broken isequal(SBMLToolkit.get_massaction(k1 * s1 / c1, [s1], [1]), k1 / c1)  # Case hOSU=false
+@test isequal(SBMLToolkit.get_massaction(k1 * s1 / c1, [s1], [1]), k1 / c1)  # Case hOSU=false
 @test isequal(SBMLToolkit.get_massaction(k1 + c1, nothing, nothing), k1 + c1)  # Case zero order kineticLaw
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2 / (c1 + s2), [s1], [1]))  # Case Michaelis-Menten kinetics
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * IV, [s1], [1]))  # Case kineticLaw with time