omit conversion to mass action for stoichiometries > 100 to avoid Sta…

…ckOverflow in simplify_fractions

README.md

@@ -27,28 +27,43 @@ Pkg.add("SBMLToolkit")
 SBML models can be simulated with the following steps (note that `sol` is in absolute quantities rather than concentrations):
 
 ```julia
-using SBMLToolkit, ModelingToolkit, OrdinaryDiffEq
+using SBMLToolkit, OrdinaryDiffEq
 
-SBMLToolkit.checksupport_file("my_model.xml")
-mdl = readSBML("my_model.xml", doc -> begin
-    set_level_and_version(3, 2)(doc)
-    convert_promotelocals_expandfuns(doc)
-end)
-
-rs = ReactionSystem(mdl)  # If you want to create a reaction system
-odesys = convert(ODESystem, rs)  # Alternatively: ODESystem(mdl)
+odesys = readSBML("my_model.xml", ODESystemImporter())
 
 tspan = (0.0, 1.0)
 prob = ODEProblem(odesys, [], tspan, [])
 sol = solve(prob, Tsit5())
 ```
 
-Alternatively, SBMLToolkit also provides more concise methods to import `SBML.Models`, `Catalyst.ReactionSystems`, and `ModelingToolkit.ODESystems`.
+While this imports an `ODESystem` directly, you can also import a Catalyst.jl `ReactionSystem`:
 
 ```julia
-mdl = readSBML("my_model.xml", DefaultImporter())
+using SBMLToolkit
+
 rs = readSBML("my_model.xml", ReactionSystemImporter())
-odesys = readSBML("my_model.xml", ODESystemImporter())
+```
+
+One common case where this is useful is if you want to run stochastic instead of ODE simulations.
+
+In the very unlikely case that you need fine-grained control over the SBML import, you can create an SBML.jl `Model` (we strongly recommend manually running `SBMLToolkit.checksupport_file("my_model.xml")` before)
+
+```julia
+using SBML
+
+mdl = readSBML("my_model.xml", doc -> begin
+    set_level_and_version(3, 2)(doc)
+    convert_promotelocals_expandfuns(doc)
+end)
+```
+
+The conversion to SBML level 3 version 2 is necessary, because older versions are not well supported in SBMLToolkit. `convert_promotelocals_expandfuns` basically flattens the SBML before the import. Once you have obtained the `Model`, you can convert it to a `ReactionSystem` and `ODESystem`.
+
+```julia
+using SBMLToolkit
+
+rs = ReactionSystem(mdl)
+odesys = convert(ODESystem, rs)
 ```
 
 ## License

src/reactions.jl

@@ -19,6 +19,7 @@ function get_reactions(model::SBML.Model)
                 enforce_rate = enforce_rate)
         else
             kl = substitute(symbolic_math, subsdict)  # Todo: use SUBSDICT
+            println(reaction)
             add_reaction!(rxs, kl, reactant_references, product_references, model)
         end
     end
@@ -61,7 +62,9 @@ function add_reaction!(rxs::Vector{Reaction},
     reactants, products, rstoichvals, pstoichvals = get_reagents(reactant_references,
         product_references, model)
     isnothing(reactants) && isnothing(products) && return
+    println("rstoichvals: $rstoichvals")
     rstoichvals = stoich_convert_to_ints(rstoichvals)
+    println("rstoichvals: $rstoichvals")
     pstoichvals = stoich_convert_to_ints(pstoichvals)
     kl, our = use_rate(kl, reactants, rstoichvals)
     our = enforce_rate ? true : our
@@ -183,7 +186,17 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing},
         rate_const = kl
     elseif isnothing(reactants) | isnothing(stoich)
         throw(ErrorException("`reactants` and `stoich` are inconsistent: `reactants` are $(reactants) and `stoich` is $(stoich)."))
+    elseif max(stoich...) > 100  # simplify_fractions might StackOverflow
+        rate_const = kl
     else
+        println("kl")
+        println(kl)
+        println("reactants")
+        println(reactants)
+        println(stoich)
+        println(*((.^(reactants, stoich))...))
+        println("fraction")
+        println(kl / *((.^(reactants, stoich))...))
         rate_const = SymbolicUtils.simplify_fractions(kl / *((.^(reactants, stoich))...))
     end
 

test/reactions.jl

@@ -135,6 +135,7 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r))
 @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
+@test isequal(SBMLToolkit.get_massaction(k1 * (s1 + s1), [s1], [101]), k1 * (s1 + s1))  # Case no simplication performed due to large rstoich
 
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1]))
 @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1]))