Pl/symbolic bool

Project.toml

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

src/drafts.jl

@@ -1,5 +1,5 @@
 function is_event_assignment(k, model)
-    for ev in values(model.events)
+    for ev in last.(model.events)
         for as in ev.event_assignments
             if as.variable == k
                 return true

src/reactions.jl

@@ -177,6 +177,7 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing},
 
     check_args(::Val{false}, x::SymbolicUtils.BasicSymbolic{Real}) = isspecies(x) ? NaN : 0  # Species vs Parameter leaf node
     check_args(::Real) = 0  # Real leaf node
+    check_args(::SymbolicUtils.BasicSymbolic{Bool}) = NaN
     check_args(x) = throw(ErrorException("Cannot handle $(typeof(x)) types."))  # Unknow leaf node
 
     if isnothing(reactants) && isnothing(stoich)

src/rules.jl

@@ -21,12 +21,56 @@ function get_rules(model)
     algeqs, obseqs, raterules
 end
 
+function extensive_kinetic_math(m::SBML.Model, formula::SBML.Math)
+    SBML.interpret_math(  # TODO: move this to SBML.jl
+        formula,
+        map_apply = (x, rec) -> SBML.MathApply(x.fn, rec.(x.args)),
+        map_const = identity,
+        map_ident = (x::SBML.MathIdent) -> begin
+            haskey(m.reactions, x.id) && return m.reactions[x.id].kinetic_math
+            haskey(m.species, x.id) || return x
+            sp = m.species[x.id]
+            sp.only_substance_units && return x
+            if isnothing(m.compartments[sp.compartment].size) &&
+               !seemsdefined(sp.compartment, m)
+                if m.compartments[sp.compartment].spatial_dimensions == 0
+                    # If the compartment ID doesn't seem directly defined anywhere
+                    # and it is a zero-dimensional unsized compartment, just avoid
+                    # any sizing questions.
+                    return x
+                else
+                    # In case the compartment is expected to be defined, complain.
+                    throw(
+                        DomainError(
+                        sp.compartment,
+                        "compartment size is insufficiently defined"
+                    ),
+                    )
+                end
+            else
+                # Now we are sure that the model either has the compartment with
+                # constant size, or the definition is easily reachable. So just use
+                # the compartment ID as a variable to compute the concentration (or
+                # area-centration etc, with different dimensionalities) by dividing
+                # it.
+                return SBML.MathApply("/", [x, SBML.MathIdent(sp.compartment)])
+            end
+        end,
+        map_lambda = (x, _) -> error(
+            ErrorException("converting lambdas to extensive kinetic math is not supported"),
+        ),
+        map_time = identity,
+        map_avogadro = identity,
+        map_value = identity
+    )
+end
+
 function get_var_and_assignment(model, rule)
     if !haskey(merge(model.species, model.compartments, model.parameters), rule.variable)
         error("Cannot find target for rule with ID `$(rule.variable)`")
     end
     var = create_var(rule.variable, IV)
-    math = SBML.extensive_kinetic_math(model, rule.math)
+    math = extensive_kinetic_math(model, rule.math)
     vc = get_volume_correction(model, rule.variable)
     if !isnothing(vc)
         math = SBML.MathApply("*", [SBML.MathIdent(vc), math])

src/systems.jl

@@ -87,11 +87,21 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires
         #                         Dict(Catalyst.DEFAULT_IV.val => 0)))
         push!(raterules_subs, rhs ~ o.lhs)
     end
+    zero_rates = []
+    for (k, v) in merge(model.parameters, model.compartments)
+        if is_event_assignment(k, model)
+            if v.constant
+                ErrorException("$k appears in event assignment but should be constant.")
+            end
+            push!(zero_rates, D(create_var(k, IV; isbcspecies = true)) ~ 0)
+        end
+    end
     if haskey(kwargs, :defaults)
         defs = ModelingToolkit._merge(defs, kwargs[:defaults])
         kwargs = filter(x -> !isequal(first(x), :defaults), kwargs)
     end
-    rs = ReactionSystem([rxs..., algrules..., raterules_subs..., obsrules_rearranged...],
+    rs = ReactionSystem(
+        [rxs..., algrules..., raterules_subs..., obsrules_rearranged..., zero_rates...],
         IV, first.(u0map), first.(parammap);
         defaults = defs, name = gensym(:SBML),
         continuous_events = get_events(model),

test/reactions.jl

@@ -139,4 +139,5 @@ m = SBML.Model(species = Dict("s" => s), reactions = Dict("r1" => r))
 
 @test isnan(SBMLToolkit.get_massaction(k1 * s1 * s2, [s1], [1]))
 @test isnan(SBMLToolkit.get_massaction(k1 + c1, [s1], [1]))
+@test isnan(SBMLToolkit.get_massaction(s1 + ifelse(s1 < k1, s1, k1), [s1], [1]))
 @test_throws ErrorException SBMLToolkit.get_massaction(k1, nothing, [1])

test/simresults.jl

@@ -8,7 +8,8 @@ const case_ids = [7,  # boundary_condition
     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    # 1208, # Non-hOSU species with rateRule in variable compartment -> require MTK fix.
+    679,  # Initial value calculated by assignmentRule in compartment of non-unit size    # 1208, # Non-hOSU species with rateRule in variable compartment -> require MTK fix.
+    944  # Non-constant parameter as target of event_assignment
 ]
 
 const logdir = joinpath(@__DIR__, "logs")