Pl/constant specs as ps

Project.toml

@@ -1,7 +1,7 @@
 name = "SBMLToolkit"
 uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 authors = ["paulflang", "anandijain"]
-version = "0.1.14"
+version = "0.1.15"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
@@ -10,7 +10,7 @@ SBML = "e5567a89-2604-4b09-9718-f5f78e97c3bb"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 
 [compat]
-Catalyst = "11.0.2"
+Catalyst = "12"
 MathML = "0.1"
 SBML = "0.10.2"
 SymbolicUtils = "0.17, 0.18, 0.19"

src/reactionsystem.jl

@@ -28,20 +28,25 @@ function get_mappings(model::SBML.Model)
     u0map = Pair[]
     parammap = Pair[]
     for (k, v) in model.species
-        var = create_var(k, IV; constant=v.constant,
-                         isbc=has_rule_type(k, model, SBML.RateRule) ||
-                              has_rule_type(k, model, SBML.AssignmentRule) ||
-                              (has_rule_type(k, model, SBML.AlgebraicRule) &&
-                              (
-                                  all([netstoich(k, r) == 0 for r in values(model.reactions)]) ||
-                                  v.boundary_condition == true)
-                              )
-                        )  # To remove species that are otherwise defined
-        push!(u0map, var => inits[k])
+        if v.constant == true
+            var = create_param(k; isconstantspecies=true)
+            push!(parammap, var => inits[k])
+        else
+            var = create_var(k, IV;
+                            isbcspecies=has_rule_type(k, model, SBML.RateRule) ||
+                                has_rule_type(k, model, SBML.AssignmentRule) ||
+                                (has_rule_type(k, model, SBML.AlgebraicRule) &&
+                                (
+                                    all([netstoich(k, r) == 0 for r in values(model.reactions)]) ||
+                                    v.boundary_condition == true)
+                                )
+                            )  # To remove species that are otherwise defined
+            push!(u0map, var => inits[k])
+        end
     end
     for (k, v) in model.parameters
         if v.constant == false && SBML.seemsdefined(k, model)
-            var = create_var(k, IV; isbc=true)
+            var = create_var(k, IV; isbcspecies=true)
             push!(u0map, var => v.value)
         else
             var = create_param(k)
@@ -50,7 +55,7 @@ function get_mappings(model::SBML.Model)
     end
     for (k, v) in model.compartments
         if v.constant == false && SBML.seemsdefined(k, model)
-            var = create_var(k, IV; isbc=true)
+            var = create_var(k, IV; isbcspecies=true)
             push!(u0map, var => v.size)
         else
             var = create_param(k)
@@ -66,16 +71,25 @@ function netstoich(id, reaction)
     netstoich += get(reaction.products, id, 0)
 end
 
-""" Check if conversion to ReactionSystem is possible """
-function checksupport(file::String)
-    not_implemented = ["listOfConstraints", "</delay>", "<priority>", "spatialDimensions=\"0\""]
-    sbml = isfile(file) ? open(file) do file
+""" Check if conversion of file to ReactionSystem is possible """
+function checksupport_file(filename::String)
+    string = open(filename) do file
         read(file, String)
-    end : file
-    for item in not_implemented
-        occursin(item, sbml) && throw(ErrorException("SBML models with $item are not yet implemented."))
     end
-    occursin("<sbml xmlns:fbc=", sbml) && throw(ErrorException("This model was designed for constrained-based optimisation. Please use COBREXA.jl instead of SBMLToolkit."))
+    checksupport_string(string)
+end
+
+""" Check if conversion of xml-string to ReactionSystem is possible """
+function checksupport_string(xml::String)
+    not_implemented = ["listOfConstraints", "</delay>",
+                       "<priority>", "spatialDimensions=\"0\"",
+                       "factorial", "00387",  # Case 00387 requires event directionality
+                       "</eventAssignment>\n          <eventAssignment"]
+                       for item in not_implemented
+        occursin(item, xml) && throw(ErrorException("SBML models with $item are not yet implemented."))
+    end
+    occursin("<sbml xmlns:fbc=", xml) && throw(ErrorException("This model was designed for constrained-based optimisation. Please use COBREXA.jl instead of SBMLToolkit."))
+    !(occursin("<reaction", xml) || occursin("rateRule", xml)) && throw(ErrorException("Models that contain neither reactions or rateRules will fail in simulation."))
     true
 end
 
@@ -310,17 +324,17 @@ function get_substitutions(model)
     subsdict
 end
 
-function create_var(x; constant=false, isbc=false)
+function create_var(x; isbcspecies=false)
     sym = Symbol(x)
-    Symbolics.unwrap(first(@variables $sym [isconstant=constant, isbc=isbc]))
+    Symbolics.unwrap(first(@variables $sym [isbcspecies=isbcspecies]))
 end
-function create_var(x, iv; constant=false, isbc=false)
+function create_var(x, iv; isbcspecies=false)
     sym = Symbol(x)
-    Symbolics.unwrap(first(@variables $sym(iv) [isconstant=constant, isbc=isbc]))
+    Symbolics.unwrap(first(@variables $sym(iv) [isbcspecies=isbcspecies]))
 end
-function create_param(x)
+function create_param(x; isconstantspecies=false)
     sym = Symbol(x)
-    Symbolics.unwrap(first(@parameters $sym))
+    Symbolics.unwrap(first(@parameters $sym [isconstantspecies=isconstantspecies]))
 end
 
 function has_rule_type(id::String, m::SBML.Model, T::Type{<:SBML.Rule})
@@ -354,8 +368,10 @@ function get_events(model, rs)  # Todo: implement up or downpass and parameters
                     math = SBML.MathApply("*", [SBML.MathIdent(vc), math])
                 end
             end
-            var = Symbol(eva.variable)  # Todo: try create_var(eva.variable, IV)
-            effect = ModelingToolkit.getvar(rs, var) ~ Symbolics.unwrap(interpret_as_num(math))
+            var = create_var(eva.variable, IV)
+            math = substitute(Symbolics.unwrap(interpret_as_num(math)),
+                              subsdict)
+            effect = var ~ math
             push!(mtk_evas, effect)
         end
         push!(mtk_evs, trig => mtk_evas)

test/misc.jl

@@ -1,11 +1,12 @@
 const case_ids = [7,  # boundary_condition
                   22,  # non-integer stoichiometry
+                  23,  # species with constant=boundaryCondition="true"
+                  140,  # compartment size overridden with assignmentRule
                   170,  # Model using parameters and rules only
                   325,  # One reactions and two rate rules with four species in a 2D compartment
                   679  # Initial value calculated by assignmentRule in compartment of non-unit size
                  ]
 
-# const case_ids = [325]
 const cases = map(x -> x[end-4:end], .*("0000", string.(case_ids)))
 
 const algomap = Dict("00177" => Rodas4(),
@@ -17,25 +18,22 @@ const algomap = Dict("00177" => Rodas4(),
                   "00882" => Rodas4()
                 )
 
-const special_tolerances = Dict("00201" => 100)
+const special_tolerances = Dict("00172" => 100,
+                                "00201" => 100,
+                                "00358" => 100,
+                                "00387" => 100)
 
-const ss_fail = ["00023", "00024"]
 const logdir = joinpath(@__DIR__, "logs")
 ispath(logdir) && rm(logdir,recursive=true)
 mkdir(logdir)
 
 const expected_errs = 
-    ["Model contains no reactions.",
-    "are not yet implemented.",
+    ["are not yet implemented.",
     "Please make reaction irreversible or rearrange kineticLaw to the form `term1 - term2`.",
-    "BoundsError(String[], (1,))",  # Occurs where no V3L2 file is available
     "COBREXA.jl",  # Occurs when model requires fbc package
-    "no method matching length(::Nothing)", "MethodError(iterate, (nothing,),", # Occurs for insance in case 00029, where S1(t) = 7 is the only eqn.
     "Stoichiometry must be a non-negative integer.",
-    "NaN result for non-NaN input.",  # Todo: remove this once you can handle factorials
     "RequestError(",
-    "structural_simplify"  # Todo: remove once structural_simplify works with `constant` and `isbc`.
-    ]
+    "neither reactions or rateRules"]
 
 function setup_settings_txt(text)
     ls = split(text, "\n")
@@ -48,11 +46,12 @@ function to_concentrations(sol, ml, res_df, ia)
     volumes = [1.]
     sol_df = DataFrame(sol)
     for sn in names(sol_df)[2:end]
-        if haskey(ml.species, sn[1:3-end])
+        if haskey(ml.species, sn[1:end-3])
             spec = ml.species[sn[1:end-3]]
             comp = ml.compartments[spec.compartment]
             ic = spec.initial_concentration
-            isnothing(ic) || haskey(ia, sn[1:end-3]) ? push!(volumes, 1.) : push!(volumes, comp.size)
+            # isnothing(ic) || haskey(ia, sn[1:end-3]) ? push!(volumes, 1.) : push!(volumes, comp.size)
+            isnothing(spec.initial_amount) ? push!(volumes, comp.size) : push!(volumes, 1.)  # Todo: see if this line works better than the above
         else
             push!(volumes, 1.)
         end
@@ -104,17 +103,18 @@ function verify_case(case; verbose=true)
     try
     # Read case
         sbml, settings, res_df = read_case(case)
-
         # Read SBML
-        SBMLToolkit.checksupport(sbml)
+        SBMLToolkit.checksupport_string(sbml)
         ml = readSBMLFromString(sbml, doc -> begin
             set_level_and_version(3, 2)(doc)
             convert_simplify_math(doc)
         end)
-        ia = readSBMLFromString(sbml, doc -> begin
-            set_level_and_version(3, 2)(doc)
-        end)
-        ia = ia.initial_assignments
+        # ia = readSBMLFromString(sbml, doc -> begin
+        #     set_level_and_version(3, 2)(doc)
+        # end)
+
+        # ia = ia.initial_assignments
+        ia = Dict()  # Todo: figure out if ia are divided by volume or not (688 isn't)
         k = 1
 
         rs = ReactionSystem(ml)
@@ -126,11 +126,11 @@ function verify_case(case; verbose=true)
         end
         k = 3
 
-        # ssys = structural_simplify(sys)
+        ssys = structural_simplify(sys)
         k = 4
 
         ts = res_df[:, 1]  # LinRange(settings["start"], settings["duration"], settings["steps"]+1)
-        prob = ODEProblem(sys, Pair[], (settings["start"], Float64(settings["duration"])); saveat=ts)
+        prob = ODEProblem(ssys, Pair[], (settings["start"], Float64(settings["duration"])); saveat=ts)
         k = 5
 
         algorithm = get(algomap, case, Sundials.CVODE_BDF())

test/reactionsystem.jl

@@ -125,10 +125,10 @@ sol = solve(oprob, Tsit5())
 @test_nowarn ODEProblem(ODESystem(readSBML(sbmlfile)), [], [0.0, 1.0], [])
 
 # Test checksupport
-@test_nowarn SBMLToolkit.checksupport(sbmlfile)
-@test_throws ErrorException SBMLToolkit.checksupport(joinpath("data", "unsupported.sbml"))
-@test_nowarn SBMLToolkit.checksupport("all good")
-@test_throws ErrorException SBMLToolkit.checksupport("contains </delay>")
+@test_nowarn SBMLToolkit.checksupport_file(sbmlfile)
+@test_throws ErrorException SBMLToolkit.checksupport_file(joinpath("data", "unsupported.sbml"))
+@test_nowarn SBMLToolkit.checksupport_string("all good <reaction")
+@test_throws ErrorException SBMLToolkit.checksupport_string("contains </delay>")
 
 # Test get_substitutions
 truesubs = Dict(Num(Symbolics.variable(:c1; T = Real)) => c1,
@@ -205,7 +205,7 @@ m = SBML.Model(
     parameters = Dict("k2" => p),
     rules = SBML.Rule[SBML.RateRule("k2", KINETICMATH2)])
 u0map, parammap = SBMLToolkit.get_mappings(m)
-u0map_true = [SBMLToolkit.create_var("k2", IV; isbc=true) => 1.0]
+u0map_true = [SBMLToolkit.create_var("k2", IV; isbcspecies=true) => 1.0]
 @test isequal(u0map, u0map_true)
 @test Catalyst.isbc(first(u0map[1]))
 
@@ -225,6 +225,3 @@ Catalyst.isbc(first(u0map[1]))
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1]))
 @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1]))
 @test_throws ErrorException SBMLToolkit.get_massaction(k1, nothing, [1])
-
-# default test
-@test_broken ModelingToolkit.defaults(m) == ModelingToolkit.defaults(ReactionSystem(m))