Merge pull request #163 from SciML/pl/support_catalystv14

Support Catalyst v14

Author: ChrisRackauckas

Project.toml

@@ -11,14 +11,14 @@ SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 [compat]
 Aqua = "0.8"
 CSV = "0.10.5"
-Catalyst = "13"
+Catalyst = "14"
 DataFrames = "1"
 Downloads = "1"
-ModelingToolkit = "8.51"
+ModelingToolkit = "9"
 OrdinaryDiffEq = "6.42"
 Plots = "1.11"
 SBML = "1.4.4"
-SBMLToolkitTestSuite = "0.0.3"
+SBMLToolkitTestSuite = "0.0.4"
 SafeTestsets = "0.1"
 Sundials = "4.14"
 SymbolicUtils = "1, 2"

README.md

@@ -11,7 +11,7 @@
 
 SBMLToolkit.jl is a lightweight tool to import models specified in the Systems Biology Markup Language (SBML) into the Julia SciML ecosystem. There are multiple ways to specify the same model in SBML. Please help us improving SBMLToolkit.jl by creating a GitHub issue if you experience errors when converting your SBML model.
 
-SBMLToolkit uses the [SBML.jl](https://github.com/LCSB-BioCore/SBML.jl) wrapper of the [libSBML](https://model.caltech.edu/software/libsbml/) library to lower dynamical SBML models into dynamical systems. If you would like to import BioNetGen models use the `writeSBML()` export function or import the `.net` file with [ReactionNetworkImporters.jl](https://github.com/SciML/ReactionNetworkImporters.jl). For constrained-based modeling, please have a look at [COBREXA.jl](https://github.com/LCSB-BioCore/COBREXA.jl). We also recommend trying [SBMLImporter.jl](https://github.com/sebapersson/SBMLImporter.jl) (Pros: respects directionality of events, outputs concentrations instead of amounts. Cons: imports ODESystems, i.e. does not interface with Catalyst).
+SBMLToolkit uses the [SBML.jl](https://github.com/LCSB-BioCore/SBML.jl) wrapper of the [libSBML](https://model.caltech.edu/software/libsbml/) library to lower dynamical SBML models into completed dynamical systems. If you would like to import BioNetGen models use the `writeSBML()` export function or import the `.net` file with [ReactionNetworkImporters.jl](https://github.com/SciML/ReactionNetworkImporters.jl). For constrained-based modeling, please have a look at [COBREXA.jl](https://github.com/LCSB-BioCore/COBREXA.jl). We also recommend trying [SBMLImporter.jl](https://github.com/sebapersson/SBMLImporter.jl). While SBMLToolkit.jl has a slightly cleaner interface, SBMLImporter.jl respects directionality of events, can output concentrations or amounts, and provides better simulation performance for models including time-triggered events and SBML `piecewise` expressions.
 
 ## Installation
 
@@ -46,7 +46,7 @@ rs = readSBML("my_model.xml", ReactionSystemImporter())
 
 One common case where this is useful is if you want to run stochastic instead of ODE simulations.
 
-In the very unlikely case that you need fine-grained control over the SBML import, you can create an SBML.jl `Model` (we strongly recommend manually running `SBMLToolkit.checksupport_file("my_model.xml")` before)
+In the very unlikely case that you need fine-grained control over the SBML import, you can create an SBML.jl `Model` (we strongly recommend manually running `checksupport_file("my_model.xml")` before)
 
 ```julia
 using SBML

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, checksupport_file
 export DefaultImporter, ReactionSystemImporter, ODESystemImporter
 
 end

src/reactions.jl

@@ -30,7 +30,7 @@ Infer forward and reverse components of bidirectional kineticLaw
 """
 function get_unidirectional_components(bidirectional_math)
     bm = ModelingToolkit.value(bidirectional_math)  #  Symbolics.tosymbol(bidirectional_math)
-    bm = simplify(bm; expand = true)
+    bm = expand(simplify(bm))
     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))
@@ -168,7 +168,7 @@ function get_massaction(kl::Num, reactants::Union{Vector{Num}, Nothing},
 
     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)
+            if isnan(check_args(arg)) || isequal(arg, default_t())
                 return NaN
             end
         end

src/systems.jl

@@ -39,8 +39,9 @@ See also [`Model`](@ref) and [`ODESystemImporter`](@ref).
 """
 function SBML.readSBML(sbmlfile::String, ::ODESystemImporter;
         include_zero_odes::Bool = true, kwargs...)  # Returns an MTK.ODESystem
-    convert(ODESystem, readSBML(sbmlfile, ReactionSystemImporter(), kwargs...),
+    odesys = convert(ODESystem, readSBML(sbmlfile, ReactionSystemImporter(), kwargs...),
         include_zero_odes = include_zero_odes)
+    complete(odesys)
 end
 
 """
@@ -90,11 +91,12 @@ function Catalyst.ReactionSystem(model::SBML.Model; kwargs...)  # Todo: requires
         defs = ModelingToolkit._merge(defs, kwargs[:defaults])
         kwargs = filter(x -> !isequal(first(x), :defaults), kwargs)
     end
-    ReactionSystem([rxs..., algrules..., raterules_subs..., obsrules_rearranged...],
+    rs = 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...)
+    return complete(rs)  # Todo: maybe add a `complete=True` kwarg
 end
 
 """
@@ -107,7 +109,8 @@ See also [`ReactionSystem`](@ref).
 function ModelingToolkit.ODESystem(model::SBML.Model; include_zero_odes::Bool = true,
         kwargs...)
     rs = ReactionSystem(model; kwargs...)
-    convert(ODESystem, rs; include_zero_odes = include_zero_odes)
+    odesys = convert(ODESystem, rs; include_zero_odes = include_zero_odes)
+    complete(odesys)
 end
 
 function get_mappings(model::SBML.Model)

src/utils.jl

@@ -1,6 +1,6 @@
 # Conversion to symbolics
-const IV = Catalyst.DEFAULT_IV
-const D = Differential(IV)
+const IV = default_t()
+const D = default_time_deriv()
 symbolicsRateOf(x) = D(x)
 const symbolics_mapping = Dict(SBML.default_function_mapping...,
     "rateOf" => symbolicsRateOf)

test/events.jl

@@ -2,7 +2,7 @@ using SBMLToolkit
 using Catalyst, SBMLToolkitTestSuite
 using Test
 
-const IV = Catalyst.DEFAULT_IV
+const IV = default_t()
 @parameters compartment
 @species S1(IV) S2(IV)
 

test/reactions.jl

@@ -4,7 +4,7 @@ using Test
 
 cd(@__DIR__)
 sbmlfile = joinpath("data", "reactionsystem_01.xml")
-const IV = Catalyst.DEFAULT_IV
+const IV = default_t()
 @parameters k1, c1
 @species s1(IV), s2(IV), s1s2(IV)
 

test/rules.jl

@@ -3,7 +3,7 @@ using SBML, SBMLToolkitTestSuite
 using Catalyst, ModelingToolkit, OrdinaryDiffEq
 using Test
 
-const IV = Catalyst.DEFAULT_IV
+const IV = default_t()
 @parameters k1, compartment
 @species S1(IV), S2(IV)
 
@@ -24,7 +24,7 @@ o_true = [S1 ~ 7 * compartment]
 sbml, _, _ = SBMLToolkitTestSuite.read_case("00031")  # rateRule
 m = readmodel(sbml)
 a, o, r = SBMLToolkit.get_rules(m)
-r_true = [Differential(IV)(S1) ~ 7 * compartment]
+r_true = [default_time_deriv()(S1) ~ 7 * compartment]
 @test isequal(r, r_true)
 
 sbml, _, _ = SBMLToolkitTestSuite.read_case("00039")  # algebraicRule
@@ -59,16 +59,16 @@ m = readmodel(sbml)
 vc = SBMLToolkit.get_volume_correction(m, "S1")
 @test isnothing(vc)
 
-# tests that non-constant parameters become states
+# tests that non-constant parameters become variables
 sbml, _, _ = SBMLToolkitTestSuite.read_case("00033")
 m = readmodel(sbml)
 @named sys = ODESystem(m)
 @species k1(IV)
-@test isequal(k1, states(sys)[end])
+@test isequal(k1, unknowns(sys)[end])
 
-# tests that non-constant compartments become states
+# tests that non-constant compartments become variables
 sbml, _, _ = SBMLToolkitTestSuite.read_case("00051")  # hOSU="true" species
 m = readmodel(sbml)
 @named sys = ODESystem(m)
 @species C(IV)
-@test isequal(C, states(sys)[end])
+@test isequal(C, unknowns(sys)[end])

test/systems.jl

@@ -4,7 +4,7 @@ using Test
 
 cd(@__DIR__)
 sbmlfile = joinpath("data", "reactionsystem_01.xml")
-const IV = Catalyst.DEFAULT_IV
+const IV = default_t()
 @parameters k1, c1
 @species s1(IV), s2(IV), s1s2(IV)
 
@@ -49,7 +49,7 @@ rs = ReactionSystem(MODEL1)
 @test isequal(Catalyst.get_eqs(rs),
     Catalyst.Reaction[Catalyst.Reaction(k1, nothing, [s1], nothing, [1.0])])
 @test isequal(Catalyst.get_iv(rs), IV)
-@test isequal(Catalyst.get_states(rs), [s1])
+@test isequal(Catalyst.get_species(rs), [s1])
 @test isequal(Catalyst.get_ps(rs), [k1, c1])
 @named rs = ReactionSystem(MODEL1)
 isequal(nameof(rs), :rs)
@@ -59,7 +59,7 @@ rs = ReactionSystem(readSBML(sbmlfile))
     Catalyst.Reaction[Catalyst.Reaction(k1 / c1, [s1, s2], [s1s2], [1.0, 1.0],
         [1.0])])
 @test isequal(Catalyst.get_iv(rs), IV)
-@test isequal(Catalyst.get_states(rs), [s1, s1s2, s2])
+@test isequal(Catalyst.get_species(rs), [s1, s1s2, s2])
 @test isequal(Catalyst.get_ps(rs), [k1, c1])
 @named rs = ReactionSystem(MODEL1)
 isequal(nameof(rs), :rs)
@@ -71,18 +71,18 @@ rs = ReactionSystem(MODEL2)  # Contains reversible reaction
         Catalyst.Reaction(k1, nothing, [s1],
             nothing, [1])])
 @test isequal(Catalyst.get_iv(rs), IV)
-@test isequal(Catalyst.get_states(rs), [s1])
+@test isequal(Catalyst.get_species(rs), [s1])
 @test isequal(Catalyst.get_ps(rs), [k1, c1])
 
 @test convert(ModelingToolkit.ODESystem, rs) isa ODESystem
 @test structural_simplify(convert(ModelingToolkit.ODESystem, rs)) isa ODESystem
 
 # Test ODESystem constructor
 odesys = ODESystem(MODEL1)
-trueeqs = Equation[Differential(IV)(s1) ~ k1]
+trueeqs = Equation[default_time_deriv()(s1) ~ k1]
 @test isequal(Catalyst.get_eqs(odesys), trueeqs)
 @test isequal(Catalyst.get_iv(odesys), IV)
-@test isequal(Catalyst.get_states(odesys), [s1])
+@test isequal(Catalyst.get_unknowns(odesys), [s1])
 @test isequal(Catalyst.get_ps(odesys), [k1, c1])
 u0 = [s1 => 1.0]
 par = [k1 => 1.0, c1 => 2.0]
@@ -93,12 +93,12 @@ isequal(nameof(odesys), :odesys)
 
 odesys = ODESystem(readSBML(sbmlfile))
 m = readSBML(sbmlfile)
-trueeqs = Equation[Differential(IV)(s1) ~ -(k1 * s1 * s2) / c1,
-    Differential(IV)(s1s2) ~ (k1 * s1 * s2) / c1,
-    Differential(IV)(s2) ~ -(k1 * s1 * s2) / c1]
+trueeqs = Equation[default_time_deriv()(s1) ~ -((k1 * s1 * s2) / c1),
+    default_time_deriv()(s1s2) ~ (k1 * s1 * s2) / c1,
+    default_time_deriv()(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 isequal(Catalyst.get_unknowns(odesys), [s1, s1s2, s2])
 @test isequal(Catalyst.get_ps(odesys), [k1, c1])
 u0 = [s1 => 2 * 1.0, s2 => 2 * 1.0, s1s2 => 2 * 1.0]
 par = [k1 => 1.0, c1 => 2.0]