Merge pull request #130 from SciML/pl/fix_map_symbolics_ident

return correct variable type/metadata in map_symbolics_ident

Author: paulflang

Project.toml

@@ -1,7 +1,7 @@
 name = "SBMLToolkit"
 uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 authors = ["paulflang", "anandijain"]
-version = "0.1.23"
+version = "0.1.24"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
@@ -12,7 +12,7 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 [compat]
 Catalyst = "13"
 MathML = "0.1"
-SBML = "1.4"
+SBML = "1.4.4"
 SymbolicUtils = "1"
 julia = "1.6"
 

README.md

@@ -31,9 +31,9 @@ using SBMLToolkit, ModelingToolkit, 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)
+    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)

docs/make.jl

@@ -7,27 +7,27 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
 DocMeta.setdocmeta!(SBMLToolkit, :DocTestSetup, :(using SBMLToolkit); recursive = true)
 
 makedocs(;
-         modules = [SBMLToolkit],
-         authors = "paulflang, anandijain",
-         sitename = "SBMLToolkit.jl",
-         format = Documenter.HTML(;
-                                  prettyurls = get(ENV, "CI", "false") == "true",
-                                  canonical = "https://docs.sciml.ai/SBMLToolkit/stable/",
-                                  assets = ["assets/favicon.ico"]),
-         clean = true, doctest = false, linkcheck = true,
-         linkcheck_ignore = ["https://www.linkedin.com/in/paul-lang-7b54a81a3/"],
-         strict = [
-             :doctest,
-             :linkcheck,
-             :parse_error,
-             :example_block,
-             # Other available options are
-             # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block
-         ],
-         pages = [
-             "Home" => "index.md",
-             "API documentation" => "api.md",
-         ])
+    modules = [SBMLToolkit],
+    authors = "paulflang, anandijain",
+    sitename = "SBMLToolkit.jl",
+    format = Documenter.HTML(;
+        prettyurls = get(ENV, "CI", "false") == "true",
+        canonical = "https://docs.sciml.ai/SBMLToolkit/stable/",
+        assets = ["assets/favicon.ico"]),
+    clean = true, doctest = false, linkcheck = true,
+    linkcheck_ignore = ["https://www.linkedin.com/in/paul-lang-7b54a81a3/"],
+    strict = [
+        :doctest,
+        :linkcheck,
+        :parse_error,
+        :example_block,
+        # Other available options are
+        # :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block
+    ],
+    pages = [
+        "Home" => "index.md",
+        "API documentation" => "api.md",
+    ])
 
 deploydocs(;
-           repo = "github.com/SciML/SBMLToolkit.jl")
+    repo = "github.com/SciML/SBMLToolkit.jl")

docs/src/index.md

@@ -42,9 +42,9 @@ using SBMLToolkit, ModelingToolkit, 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)
+    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)
 

src/SBMLToolkit.jl

@@ -15,7 +15,7 @@ include("utils.jl")
 
 export ReactionSystem, ODESystem
 export readSBML, readSBMLFromString, set_level_and_version, convert_simplify_math,
-       convert_promotelocals_expandfuns
+    convert_promotelocals_expandfuns
 export DefaultImporter, ReactionSystemImporter, ODESystemImporter
 
 end

src/events.jl

@@ -7,7 +7,7 @@ function get_events(model)  # Todo: implement up or downpass and parameters
     mtk_evs = Pair{Vector{Equation}, Vector{Equation}}[]
     for (_, e) in evs
         trigger = SBML.extensive_kinetic_math(model, e.trigger.math)
-        trigger = Symbolics.unwrap(interpret_as_num(trigger))
+        trigger = Symbolics.unwrap(interpret_as_num(trigger, model))
         lhs, rhs = map(x -> substitute(x, subsdict), trigger.arguments)
         trig = [lhs ~ rhs]
         mtk_evas = Equation[]
@@ -20,8 +20,8 @@ function get_events(model)  # Todo: implement up or downpass and parameters
                 end
             end
             var = create_var(eva.variable, IV; irreducible = true)
-            math = substitute(Symbolics.unwrap(interpret_as_num(math)),
-                              subsdict)
+            math = substitute(Symbolics.unwrap(interpret_as_num(math, model)),
+                subsdict)
             effect = var ~ math
             push!(mtk_evas, effect)
         end

src/reactions.jl

@@ -6,17 +6,17 @@ function get_reactions(model::SBML.Model)
     rxs = Reaction[]
     for reaction in values(model.reactions)
         extensive_math = SBML.extensive_kinetic_math(model, reaction.kinetic_math)
-        symbolic_math = interpret_as_num(extensive_math)
+        symbolic_math = interpret_as_num(extensive_math, model)
         reactant_references = reaction.reactants
         product_references = reaction.products
         if reaction.reversible
             symbolic_math = get_unidirectional_components(symbolic_math)
             kl_fw, kl_rv = [substitute(x, subsdict) for x in symbolic_math]
             enforce_rate = isequal(kl_rv, 0)
             add_reaction!(rxs, kl_fw, reactant_references, product_references, model;
-                          enforce_rate = enforce_rate)
+                enforce_rate = enforce_rate)
             add_reaction!(rxs, kl_rv, product_references, reactant_references, model;
-                          enforce_rate = enforce_rate)
+                enforce_rate = enforce_rate)
         else
             kl = substitute(symbolic_math, subsdict)  # Todo: use SUBSDICT
             add_reaction!(rxs, kl, reactant_references, product_references, model)
@@ -53,21 +53,21 @@ function get_unidirectional_components(bidirectional_math)
 end
 
 function add_reaction!(rxs::Vector{Reaction},
-                       kl::Num,
-                       reactant_references::Vector{SBML.SpeciesReference},
-                       product_references::Vector{SBML.SpeciesReference},
-                       model::SBML.Model;
-                       enforce_rate = false)
+    kl::Num,
+    reactant_references::Vector{SBML.SpeciesReference},
+    product_references::Vector{SBML.SpeciesReference},
+    model::SBML.Model;
+    enforce_rate = false)
     reactants, products, rstoichvals, pstoichvals = get_reagents(reactant_references,
-                                                                 product_references, model)
+        product_references, model)
     isnothing(reactants) && isnothing(products) && return
     rstoichvals = stoich_convert_to_ints(rstoichvals)
     pstoichvals = stoich_convert_to_ints(pstoichvals)
     kl, our = use_rate(kl, reactants, rstoichvals)
     our = enforce_rate ? true : our
     push!(rxs,
-          Catalyst.Reaction(kl, reactants, products, rstoichvals, pstoichvals;
-                            only_use_rate = our))
+        Catalyst.Reaction(kl, reactants, products, rstoichvals, pstoichvals;
+            only_use_rate = our))
 end
 
 function stoich_convert_to_ints(xs)
@@ -78,8 +78,8 @@ end
 Get reagents
 """
 function get_reagents(reactant_references::Vector{SBML.SpeciesReference},
-                      product_references::Vector{SBML.SpeciesReference},
-                      model::SBML.Model)
+    product_references::Vector{SBML.SpeciesReference},
+    model::SBML.Model)
     reactants = String[]
     products = String[]
     rstoich = Float64[]
@@ -146,7 +146,7 @@ end
 Get kineticLaw for use in MTK.Reaction
 """
 function use_rate(kl::Num, react::Union{Vector{Num}, Nothing},
-                  stoich::Union{Vector{<:Real}, Nothing})
+    stoich::Union{Vector{<:Real}, Nothing})
     rate_const = get_massaction(kl, react, stoich)
     if !isnan(rate_const)
         kl = rate_const
@@ -161,7 +161,7 @@ end
 Get rate constant of mass action kineticLaws
 """
 function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing},
-                        stoich::Union{Vector{<:Real}, Nothing})
+    stoich::Union{Vector{<:Real}, Nothing})
     function check_args(x::SymbolicUtils.BasicSymbolic{Real})
         check_args(Val(SymbolicUtils.istree(x)), x)
     end

src/rules.jl

@@ -5,7 +5,7 @@ function get_rules(model)
     raterules = Equation[]
     for r in model.rules
         if r isa SBML.AlgebraicRule
-            push!(algeqs, 0 ~ interpret_as_num(r.math))
+            push!(algeqs, 0 ~ interpret_as_num(r.math, model))
         elseif r isa SBML.AssignmentRule
             var, ass = get_var_and_assignment(model, r)
             push!(obseqs, var ~ ass)
@@ -17,7 +17,7 @@ function get_rules(model)
         end
     end
     algeqs, obseqs, raterules = map(x -> substitute(x, subsdict),
-                                    (algeqs, obseqs, raterules))
+        (algeqs, obseqs, raterules))
     algeqs, obseqs, raterules
 end
 
@@ -31,7 +31,7 @@ function get_var_and_assignment(model, rule)
     if !isnothing(vc)
         math = SBML.MathApply("*", [SBML.MathIdent(vc), math])
     end
-    assignment = interpret_as_num(math)
+    assignment = interpret_as_num(math, model)
     if rule isa SBML.RateRule && haskey(model.species, rule.variable)
         sp = model.species[rule.variable]
         comp = model.compartments[sp.compartment]
@@ -52,6 +52,6 @@ function get_volume_correction(model, s_id)
     isnothing(comp.size) && !SBML.seemsdefined(sp.compartment, model) &&
         comp.spatial_dimensions != 0 &&  # remove this line when SBML test suite is fixed
         throw(DomainError(sp.compartment,
-                          "compartment size is insufficiently defined"))
+            "compartment size is insufficiently defined"))
     sp.compartment
 end

src/systems.jl

@@ -38,9 +38,9 @@ Create a `ModelingToolkit.ODESystem` from an SBML file, using default import set
 See also [`Model`](@ref) and [`ODESystemImporter`](@ref).
 """
 function SBML.readSBML(sbmlfile::String, ::ODESystemImporter;
-                       include_zero_odes::Bool = true, kwargs...)  # Returns an MTK.ODESystem
+    include_zero_odes::Bool = true, kwargs...)  # Returns an MTK.ODESystem
     convert(ODESystem, readSBML(sbmlfile, ReactionSystemImporter(), kwargs...),
-            include_zero_odes = include_zero_odes)
+        include_zero_odes = include_zero_odes)
 end
 
 """
@@ -54,7 +54,11 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires
     # 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)
-    defs = ModelingToolkit._merge(Dict(u0map), Dict(parammap))
+    defs = Dict{Num, Any}()
+    for (k, v) in vcat(u0map, parammap)
+        defs[k] = v
+    end
+    # defs = ModelingToolkit._merge(Dict(u0map), Dict(parammap))
     algrules, obsrules, raterules = get_rules(model)
     obsrules_rearranged = Equation[]
     for o in obsrules
@@ -64,9 +68,9 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires
                 rhs = r.rhs
             end
         end
-        defs[o.lhs] = substitute(rhs,
-                                 ModelingToolkit._merge(defs,
-                                                        Dict(Catalyst.DEFAULT_IV.val => 0)))
+        defs[o.lhs] = ModelingToolkit.fixpoint_sub(rhs, defs)
+        #  ModelingToolkit._merge(defs,
+        #                         Dict(Catalyst.DEFAULT_IV.val => 0)))
         push!(obsrules_rearranged, 0 ~ rhs - o.lhs)
     end
     raterules_subs = []
@@ -77,20 +81,20 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires
                 rhs = r.rhs
             end
         end
-        defs[o.lhs] = substitute(rhs,
-                                 ModelingToolkit._merge(defs,
-                                                        Dict(Catalyst.DEFAULT_IV.val => 0)))
+        defs[o.lhs] = ModelingToolkit.fixpoint_sub(rhs, defs)
+        #  ModelingToolkit._merge(defs,
+        #                         Dict(Catalyst.DEFAULT_IV.val => 0)))
         push!(raterules_subs, rhs ~ o.lhs)
     end
     if haskey(kwargs, :defaults)
         defs = ModelingToolkit._merge(defs, kwargs[:defaults])
         kwargs = filter(x -> !isequal(first(x), :defaults), kwargs)
     end
     ReactionSystem([rxs..., algrules..., raterules_subs..., obsrules_rearranged...],
-                   IV, first.(u0map), first.(parammap);
-                   defaults = defs, name = gensym(:SBML),
-                   continuous_events = get_events(model),
-                   combinatoric_ratelaws = false, kwargs...)
+        IV, first.(u0map), first.(parammap);
+        defaults = defs, name = gensym(:SBML),
+        continuous_events = get_events(model),
+        combinatoric_ratelaws = false, kwargs...)
 end
 
 """
@@ -101,7 +105,7 @@ Create an `ODESystem` from an `SBML.Model`.
 See also [`ReactionSystem`](@ref).
 """
 function ModelingToolkit.ODESystem(model::SBML.Model; include_zero_odes::Bool = true,
-                                   kwargs...)
+    kwargs...)
     rs = ReactionSystem(model; kwargs...)
     convert(ODESystem, rs; include_zero_odes = include_zero_odes)
 end
@@ -116,12 +120,12 @@ function get_mappings(model::SBML.Model)
             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
+                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
@@ -133,7 +137,7 @@ function get_mappings(model::SBML.Model)
         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))
+            push!(parammap, var => interpret_as_num(val, model))
         else
             var = create_param(k)
             push!(parammap, var => v.value)
@@ -154,9 +158,9 @@ end
 function netstoich(id, reaction)
     netstoich = 0
     rdict = Dict(getproperty.(reaction.reactants, :species) .=>
-                     getproperty.(reaction.reactants, :stoichiometry))
+        getproperty.(reaction.reactants, :stoichiometry))
     pdict = Dict(getproperty.(reaction.products, :species) .=>
-                     getproperty.(reaction.products, :stoichiometry))
+        getproperty.(reaction.products, :stoichiometry))
     netstoich -= get(rdict, id, 0)
     netstoich += get(pdict, id, 0)
 end