Support initialAssignments with the DefaultImporter

Project.toml

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

src/systems.jl

@@ -53,9 +53,9 @@
 function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires unique parameters (i.e. SBML must have been imported with localParameter promotion in libSBML)
     # length(model.events) > 0 ? error("Model contains events. Please import with `ODESystem(model)`") : nothing  @Anand: how to suppress this when called from ODESystem
     rxs = get_reactions(model)
-    u0map, parammap = get_mappings(model)
+    u0map, parammap, initial_assignment_map = get_mappings(model)
     defs = Dict{Num, Any}()
-    for (k, v) in vcat(u0map, parammap)
+    for (k, v) in vcat(u0map, parammap, initial_assignment_map)  # initial_assignments override u0map and parammap
         defs[k] = v
     end
     # defs = ModelingToolkit._merge(Dict(u0map), Dict(parammap))
@@ -114,45 +114,44 @@
     inits = Dict(SBML.initial_amounts(model, convert_concentrations = true))
     u0map = Pair[]
     parammap = Pair[]
+    initial_assignment_map = Pair[]
+
     for (k, v) in model.species
+        var = create_symbol(k, model)
         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
+        var = create_symbol(k, model)
         if v.constant == false &&
            (SBML.seemsdefined(k, model) || is_event_assignment(k, model))
-            var = create_var(k, IV; isbcspecies = true)
             push!(u0map, var => v.value)
         elseif v.constant == true && isnothing(v.value)  # Todo: maybe add this branch also to model.compartments
-            var = create_param(k)
             val = model.initial_assignments[k]
             push!(parammap, var => interpret_as_num(val, model))
         else
-            var = create_param(k)
             push!(parammap, var => v.value)
         end
     end
+
     for (k, v) in model.compartments
+        var = create_symbol(k, model)
         if v.constant == false && SBML.seemsdefined(k, model)
-            var = create_var(k, IV; isbcspecies = true)
             push!(u0map, var => v.size)
         else
-            var = create_param(k)
             push!(parammap, var => v.size)
         end
     end
-    u0map, parammap
+
+    for (k, v) in model.initial_assignments
+        var = create_symbol(k, model)
+        push!(initial_assignment_map, var => interpret_as_num(v, model))
+    end
+    u0map, parammap, initial_assignment_map
 end
 
 function netstoich(id, reaction)

src/utils.jl

@@ -101,3 +101,36 @@ const importdefaults = doc -> begin
     set_level_and_version(3, 2)(doc)
     convert_promotelocals_expandfuns(doc)
 end
+
+function create_symbol(k::String, model::SBML.Model)
+    if k in keys(model.species)
+        v = model.species[k]
+        if v.constant == true
+            sym = create_param(k; isconstantspecies = true)
+        else
+            sym = 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
+        end
+    elseif k in keys(model.parameters)
+        v = model.parameters[k]
+        if v.constant == false &&
+           (SBML.seemsdefined(k, model) || is_event_assignment(k, model))
+            sym = create_var(k, IV; isbcspecies = true)
+        else
+            sym = create_param(k)
+        end
+    elseif k in keys(model.compartments)
+        v = model.compartments[k]
+        if v.constant == false && SBML.seemsdefined(k, model)
+            sym = create_var(k, IV; isbcspecies = true)
+        else
+            sym = create_param(k)
+        end
+    end
+    sym
+end

test/systems.jl

@@ -114,16 +114,18 @@ isequal(nameof(odesys), :odesys)
 # @test_nowarn ODEProblem(ODESystem(readSBML(sbmlfile)), [], [0.0, 1.0], [])
 
 # Test get_mappings
-u0map, parammap = SBMLToolkit.get_mappings(MODEL1)
+u0map, parammap, initial_assignment_map = SBMLToolkit.get_mappings(MODEL1)
 u0map_true = [s1 => 1.0]
 parammap_true = [k1 => 1.0, c1 => 2.0]
+initial_assignment_map_true = Pair[]
 @test isequal(u0map, u0map_true)
 @test isequal(parammap, parammap_true)
+@test isequal(initial_assignment_map, initial_assignment_map_true)
 
 p = SBML.Parameter(name = "k2", value = 1.0, constant = false)
 m = SBML.Model(parameters = Dict("k2" => p),
     rules = SBML.Rule[SBML.RateRule("k2", KINETICMATH2)])
-u0map, parammap = SBMLToolkit.get_mappings(m)
+u0map, parammap, initial_assignment_map = SBMLToolkit.get_mappings(m)
 u0map_true = [SBMLToolkit.create_var("k2", IV; isbcspecies = true) => 1.0]
 @test isequal(u0map, u0map_true)
 @test Catalyst.isbc(first(u0map[1]))
@@ -132,13 +134,15 @@ p = SBML.Parameter(name = "k2", value = nothing, constant = true)
 ia = Dict("k2" => KINETICMATH1)
 m = SBML.Model(parameters = Dict("k1" => PARAM1, "k2" => p),
     initial_assignments = ia)
-u0map, parammap = SBMLToolkit.get_mappings(m)
+u0map, parammap, initial_assignment_map = SBMLToolkit.get_mappings(m)
 parammap_true = [k1 => 1.0, SBMLToolkit.create_var("k2") => k1]
+initial_assignment_map_true = [SBMLToolkit.create_var("k2") => k1]
 @test isequal(parammap, parammap_true)
+@test isequal(initial_assignment_map, initial_assignment_map_true)
 
 m = SBML.Model(species = Dict("s2" => SPECIES2),
     rules = SBML.Rule[SBML.AlgebraicRule(KINETICMATH2)])
-u0map, parammap = SBMLToolkit.get_mappings(m)
+u0map, parammap, initial_assignment_map = SBMLToolkit.get_mappings(m)
 Catalyst.isbc(first(u0map[1]))
 
 # Test get_netstoich

test/utils.jl

@@ -85,6 +85,21 @@ par = SBMLToolkit.create_param("k")
 par = SBMLToolkit.create_param("k", isconstantspecies = true)
 @test Catalyst.isconstant(par)
 
+# Test create_symbol
+# Comment in when https://github.com/SciML/ModelingToolkit.jl/issues/2228 is fixed
+# @species B(IV) Dv(IV)
+# @parameters C D Bc
+# sym = SBMLToolkit.create_symbol("B", MODEL1)
+# @test isequal(B, sym)
+# sym = SBMLToolkit.create_symbol("Dv", MODEL1)
+# @test isequal(Dv, sym)
+# sym = SBMLToolkit.create_symbol("C", MODEL1)
+# @test isequal(C, sym)
+# sym = SBMLToolkit.create_symbol("D", MODEL1)
+# @test isequal(D, sym)
+# sym = SBMLToolkit.create_symbol("Bc", MODEL1)
+# @test isequal(Bc, sym)
+
 # Test has_rule_type
 sbml, _, _ = SBMLToolkitTestSuite.read_case("00031")  # rateRule
 m = readmodel(sbml)