Pl/catalystupdate

Project.toml

@@ -1,7 +1,7 @@
 name = "SBMLToolkit"
 uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 authors = ["paulflang", "anandijain"]
-version = "0.1.20"
+version = "0.1.21"
 
 [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 = "12, 13"
+Catalyst = "13"
 MathML = "0.1"
 SBML = "1.4"
 SymbolicUtils = "0.17, 0.18, 0.19, 1"

src/reactions.jl

@@ -29,24 +29,24 @@ end
 Infer forward and reverse components of bidirectional kineticLaw
 """
 function get_unidirectional_components(bidirectional_math)
-    bm = Symbolics.tosymbol(bidirectional_math)
+    bm = ModelingToolkit.value(bidirectional_math)  #  Symbolics.tosymbol(bidirectional_math)
     bm = simplify(bm; expand = true)
-    if (bm isa Union{Real, Symbol}) || (SymbolicUtils.operation(bm) != +)
-        @warn "Cannot separate bidirectional kineticLaw `$bidirectional_math` to forward and reverse part. Setting forward to `$bidirectional_math` and reverse to `0`."
+    if !SymbolicUtils.isadd(bm)
+        @warn "Cannot separate bidirectional kineticLaw `$bidirectional_math` to forward and reverse part. Setting forward to `$bidirectional_math` and reverse to `0`. Stochastic simulations will be inexact."
         return (bidirectional_math, Num(0))
     end
-    terms = SymbolicUtils.arguments(bm)
+    terms = SymbolicUtils.arguments(ModelingToolkit.value(bm))
     fw_terms = []
     rv_terms = []
     for term in terms
-        if (term isa SymbolicUtils.Mul) && (term.coeff < 0)
+        if SymbolicUtils.ismul(ModelingToolkit.value(term)) && (term.coeff < 0)
             push!(rv_terms, Num(-term))  # PL: @Anand: Perhaps we should to create_var(term) or so?
         else
             push!(fw_terms, Num(term))  # PL: @Anand: Perhaps we should to create_var(term) or so?
         end
     end
     if (length(fw_terms) != 1) || (length(rv_terms) != 1)
-        @warn "Cannot separate bidirectional kineticLaw `$bidirectional_math` to forward and reverse part. Setting forward to `$bidirectional_math` and reverse to `0`."
+        @warn "Cannot separate bidirectional kineticLaw `$bidirectional_math` to forward and reverse part. Setting forward to `$bidirectional_math` and reverse to `0`. Stochastic simulations will be inexact."
         return (bidirectional_math, Num(0))
     end
     return (fw_terms[1], rv_terms[1])
@@ -145,17 +145,19 @@ Get rate constant of mass action kineticLaws
 """
 function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing},
                         stoich::Union{Vector{<:Real}, Nothing})
-    function check_args(x::SymbolicUtils.Symbolic{Real})
+    function check_args(x::SymbolicUtils.BasicSymbolic{Real})
+        check_args(Val(SymbolicUtils.istree(x)), x)
+    end
+    function check_args(::Val{true}, x::SymbolicUtils.BasicSymbolic{Real})
         for arg in SymbolicUtils.arguments(x)
             if isnan(check_args(arg)) || isequal(arg, Catalyst.DEFAULT_IV)
                 return NaN
             end
         end
         return 0
     end
-    check_args(_::Term{Real, Nothing}) = NaN  # Variable leaf node
-    check_args(_::Sym{Real, Base.ImmutableDict{DataType, Any}}) = 0  # Parameter leaf node
-    check_args(_::Real) = 0  # Real leaf node
+    check_args(::Val{false}, x::SymbolicUtils.BasicSymbolic{Real}) = isspecies(x) ? NaN : 0  # Species or Parameter leaf node
+    check_args(::Real) = 0  # Real leaf node
     check_args(x) = throw(ErrorException("Cannot handle $(typeof(x)) types."))  # Unknow leaf node
     if isnothing(reactants) && isnothing(stoich)
         rate_const = kl
@@ -164,5 +166,5 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing},
     else
         rate_const = SymbolicUtils.simplify_fractions(kl / *((.^(reactants, stoich))...))
     end
-    isnan(check_args(rate_const.val)) ? NaN : rate_const
+    isnan(check_args(ModelingToolkit.value(rate_const))) ? NaN : rate_const
 end

src/systems.jl

@@ -41,12 +41,10 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires
         defs = ModelingToolkit._merge(defs, kwargs[:defaults])
         kwargs = filter(x -> !isequal(first(x), :defaults), kwargs)
     end
-    constraints_sys = ODESystem(vcat(algrules, raterules_subs, obsrules_rearranged),
-                                IV; name = gensym(:CONSTRAINTS),
-                                continuous_events = get_events(model))
-    ReactionSystem(rxs, IV, first.(u0map), first.(parammap);
+    ReactionSystem([rxs..., algrules..., raterules_subs..., obsrules_rearranged...],
+                   IV, first.(u0map), first.(parammap);
                    defaults = defs, name = gensym(:SBML),
-                   constraints = constraints_sys,
+                   continuous_events = get_events(model),
                    combinatoric_ratelaws = false, kwargs...)
 end
 

src/utils.jl

@@ -9,7 +9,7 @@ const symbolics_mapping = Dict(SBML.default_function_mapping...,
 
 map_symbolics_ident(x) = begin
     sym = Symbol(x.id)
-    first(@variables $sym)
+    first(@species $sym)
 end
 
 function interpret_as_num(x::SBML.Math)
@@ -42,11 +42,11 @@ end
 
 function create_var(x; isbcspecies = false)
     sym = Symbol(x)
-    Symbolics.unwrap(first(@variables $sym [isbcspecies = isbcspecies]))
+    Symbolics.unwrap(first(@species $sym [isbcspecies = isbcspecies]))
 end
 function create_var(x, iv; isbcspecies = false, irreducible = false)
     sym = Symbol(x)
-    Symbolics.unwrap(first(@variables $sym(iv) [
+    Symbolics.unwrap(first(@species $sym(iv) [
                                isbcspecies = isbcspecies,
                                irreducible = irreducible,
                            ]))

test/events.jl

@@ -4,7 +4,7 @@ using Test
 
 const IV = Catalyst.DEFAULT_IV
 @parameters compartment
-@variables S1(IV) S2(IV)
+@species S1(IV) S2(IV)
 
 function readmodel(sbml)
     SBMLToolkit.readSBMLFromString(sbml, doc -> begin

test/reactions.jl

@@ -6,7 +6,7 @@ cd(@__DIR__)
 sbmlfile = joinpath("data", "reactionsystem_01.xml")
 const IV = Catalyst.DEFAULT_IV
 @parameters k1, c1
-@variables s1(IV), s2(IV), s1s2(IV)
+@species s1(IV), s2(IV), s1s2(IV)
 
 COMP1 = SBML.Compartment("c1", true, 3, 2.0, "nl", nothing, nothing, nothing, nothing,
                          SBML.CVTerm[])
@@ -59,7 +59,7 @@ km = SBML.MathIdent("s1s2")
 sm1 = SBMLToolkit.interpret_as_num(km)
 sm2 = sm - sm1
 @test isequal(SBMLToolkit.get_unidirectional_components(sm2), (sm2, Num(0)))
-@test isequal(SBMLToolkit.get_unidirectional_components(:k1), (:k1, Num(0)))
+@test isequal(SBMLToolkit.get_unidirectional_components(k1), (k1, Num(0)))
 
 # Test add_reaction!
 rxs = Catalyst.Reaction[]

test/rules.jl

@@ -5,7 +5,7 @@ using Test
 
 const IV = Catalyst.DEFAULT_IV
 @parameters k1, compartment
-@variables S1(IV), S2(IV)
+@species S1(IV), S2(IV)
 
 function readmodel(sbml)
     SBMLToolkit.readSBMLFromString(sbml, doc -> begin
@@ -62,12 +62,12 @@ vc = SBMLToolkit.get_volume_correction(m, "S1")
 sbml, _, _ = SBMLToolkitTestSuite.read_case("00033")
 m = readmodel(sbml)
 @named sys = ODESystem(m)
-@variables k1(IV)
+@species k1(IV)
 @test isequal(k1, states(sys)[end])
 
 # tests that non-constant compartments become states
 sbml, _, _ = SBMLToolkitTestSuite.read_case("00051")  # hOSU="true" species
 m = readmodel(sbml)
 @named sys = ODESystem(m)
-@variables C(IV)
+@species C(IV)
 @test isequal(C, states(sys)[end])

test/systems.jl

@@ -6,7 +6,7 @@ cd(@__DIR__)
 sbmlfile = joinpath("data", "reactionsystem_01.xml")
 const IV = Catalyst.DEFAULT_IV
 @parameters k1, c1
-@variables s1(IV), s2(IV), s1s2(IV)
+@species s1(IV), s2(IV), s1s2(IV)
 
 COMP1 = SBML.Compartment("c1", true, 3, 2.0, "nl", nothing, nothing, nothing, nothing,
                          SBML.CVTerm[])
@@ -91,9 +91,9 @@ isequal(nameof(odesys), :odesys)
 
 odesys = ODESystem(readSBML(sbmlfile))
 m = readSBML(sbmlfile)
-trueeqs = Equation[Differential(IV)(s1) ~ -((k1 * s1 * s2) / c1),
+trueeqs = Equation[Differential(IV)(s1) ~ -(k1 * s1 * s2) / c1,
                    Differential(IV)(s1s2) ~ (k1 * s1 * s2) / c1,
-                   Differential(IV)(s2) ~ -((k1 * s1 * s2) / c1)]
+                   Differential(IV)(s2) ~ -(k1 * s1 * s2) / c1]
 @test isequal(Catalyst.get_eqs(odesys), trueeqs)
 @test isequal(Catalyst.get_iv(odesys), IV)
 @test isequal(Catalyst.get_states(odesys), [s1, s1s2, s2])

test/utils.jl

@@ -10,19 +10,19 @@ function readmodel(sbml)
 end
 
 const IV = Catalyst.DEFAULT_IV
-@variables s1(IV)
+@species s1(IV)
 # Test symbolicsRateOf
 rate = SBMLToolkit.symbolicsRateOf(s1)
 rate_true = SBMLToolkit.D(s1)
 @test isequal(rate, rate_true)
 
 # Test map_sumbolics_ident
 sym = SBMLToolkit.map_symbolics_ident(SBML.MathIdent("s1"))
-@variables s1
+@species s1
 @test isequal(sym, s1)
 
 # Test interpret_as_num
-@variables B C D Time
+@species B C D Time
 
 test = SBML.MathApply("*",
                       SBML.Math[SBML.MathApply("+",
@@ -42,7 +42,7 @@ sbml, _, _ = SBMLToolkitTestSuite.read_case("00001")
 m = readmodel(sbml)
 subsdict = SBMLToolkit.get_substitutions(m)
 @parameters k1, compartment
-@variables S1(IV), S2(IV)
+@species S1(IV), S2(IV)
 subsdict_true = Dict(Num(Symbolics.variable(:S1; T = Real)) => S1,
                      Num(Symbolics.variable(:S2; T = Real)) => S2,
                      Num(Symbolics.variable(:k1; T = Real)) => k1,
@@ -51,7 +51,7 @@ subsdict_true = Dict(Num(Symbolics.variable(:S1; T = Real)) => S1,
 
 # Test create_var
 var = SBMLToolkit.create_var("s2")
-@variables s2
+@species s2
 @test isequal(var, s2)
 
 var = SBMLToolkit.create_var("s2", IV, isbcspecies = true)