fix events and checksupport (#60)

* fix events and checksupport * fix not_implemented list and bump version * prevent a warning in tests
SciML · Feb 23, 2022 · 08230bf · 08230bf
1 parent 097f2fe
commit 08230bf
Show file tree

Hide file tree

Showing 7 changed files with 217 additions and 4 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "SBMLToolkit"
 uuid = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 authors = ["paulflang", "anandijain"]
-version = "0.1.12"
+version = "0.1.13"
 
 [deps]
 Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"

diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
@@ -18,12 +18,12 @@ end
 """ ODESystem constructor """
 function ModelingToolkit.ODESystem(model::SBML.Model; include_zero_odes = false, kwargs...)
     rs = ReactionSystem(model; kwargs...)
-    convert(ODESystem, rs; include_zero_odes = include_zero_odes)
+    convert(ODESystem, rs; include_zero_odes = include_zero_odes, continuous_events=get_events(model, rs))
 end
 
 """ Check if conversion to ReactionSystem is possible """
 function checksupport(filename::String)
-    not_implemented = ["listOfConstraints", "listOfEvents"]
+    not_implemented = ["listOfConstraints", "</delay>", "<priority>"]
     sbml = open(filename) do file
         read(file, String)
     end
@@ -365,3 +365,29 @@ function raterule_to_diffeq(model, rule)
         error()
     end
 end
+
+"""
+    Creates ContinuousVectorCallbacks
+
+
+Note that one limitation of Event support is that ReactionSystems do not have a field for it yet.
+So in order for the system to have events, you must call `ODESystem(m::SBML.Model)` rather than `convert(ODESystem, ReactionSystem(m::SBML.Model))`
+"""
+function get_events(model, rs)
+    subsdict = _get_substitutions(model)
+    evs = model.events
+    mtk_evs = Pair{Vector{Equation},Vector{Equation}}[]
+    for (_, e) in evs
+        args = convert(Num, e.trigger; convert_time = (x::SBML.MathTime) -> Catalyst.DEFAULT_IV).val.arguments
+        lhs, rhs = map(x -> substitute(x, subsdict), args)
+        trig = [lhs ~ rhs]
+        mtk_evas = Equation[]
+        for eva in e.event_assignments
+            var = Symbol(eva.variable)
+            pair = ModelingToolkit.getvar(rs, var) ~ convert(Num, eva.math)
+            push!(mtk_evas, pair)
+        end
+        push!(mtk_evs, trig => mtk_evas)
+    end
+    mtk_evs
+end
diff --git a/test/data/00026-sbml-l3v2.xml b/test/data/00026-sbml-l3v2.xml
@@ -0,0 +1,72 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
+  <model metaid="_case00026" id="case00026" name="case00026" timeUnits="time">
+    <listOfUnitDefinitions>
+      <unitDefinition id="volume">
+        <listOfUnits>
+          <unit kind="litre" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="substance">
+        <listOfUnits>
+          <unit kind="mole" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="time">
+        <listOfUnits>
+          <unit kind="second" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+    <listOfCompartments>
+      <compartment id="compartment" name="compartment" spatialDimensions="3" size="1" units="volume" constant="true"/>
+    </listOfCompartments>
+    <listOfSpecies>
+      <species id="S1" name="S1" compartment="compartment" initialAmount="1" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="S2" name="S2" compartment="compartment" initialAmount="0" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+    </listOfSpecies>
+    <listOfParameters>
+      <parameter id="k1" name="k1" value="1" constant="true"/>
+    </listOfParameters>
+    <listOfReactions>
+      <reaction id="reaction1" name="reaction1" reversible="false">
+        <listOfReactants>
+          <speciesReference species="S1" stoichiometry="1" constant="true"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="S2" stoichiometry="1" constant="true"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> compartment </ci>
+              <ci> k1 </ci>
+              <ci> S1 </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+    </listOfReactions>
+    <listOfEvents>
+      <event id="event1" name="event1" useValuesFromTriggerTime="true">
+        <trigger initialValue="true" persistent="true">
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <lt/>
+              <ci> S1 </ci>
+              <cn> 0.1 </cn>
+            </apply>
+          </math>
+        </trigger>
+        <listOfEventAssignments>
+          <eventAssignment variable="S1">
+            <math xmlns="http://www.w3.org/1998/Math/MathML">
+              <cn type="integer"> 1 </cn>
+            </math>
+          </eventAssignment>
+        </listOfEventAssignments>
+      </event>
+    </listOfEvents>
+  </model>
+</sbml>
diff --git a/test/data/00041-sbml-l3v2.xml b/test/data/00041-sbml-l3v2.xml
@@ -0,0 +1,90 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<sbml xmlns="http://www.sbml.org/sbml/level3/version2/core" level="3" version="2">
+  <model metaid="_case00041" id="case00041" name="case00041" timeUnits="time">
+    <listOfUnitDefinitions>
+      <unitDefinition id="volume">
+        <listOfUnits>
+          <unit kind="litre" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="substance">
+        <listOfUnits>
+          <unit kind="mole" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+      <unitDefinition id="time">
+        <listOfUnits>
+          <unit kind="second" exponent="1" scale="0" multiplier="1"/>
+        </listOfUnits>
+      </unitDefinition>
+    </listOfUnitDefinitions>
+    <listOfCompartments>
+      <compartment id="compartment" name="compartment" spatialDimensions="3" size="1" units="volume" constant="true"/>
+    </listOfCompartments>
+    <listOfSpecies>
+      <species id="S1" name="S1" compartment="compartment" initialAmount="1" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+      <species id="S2" name="S2" compartment="compartment" initialAmount="0" substanceUnits="substance" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
+    </listOfSpecies>
+    <listOfParameters>
+      <parameter id="k1" name="k1" value="1" constant="true"/>
+    </listOfParameters>
+    <listOfReactions>
+      <reaction id="reaction1" name="reaction1" reversible="false">
+        <listOfReactants>
+          <speciesReference species="S1" stoichiometry="1" constant="true"/>
+        </listOfReactants>
+        <listOfProducts>
+          <speciesReference species="S2" stoichiometry="1" constant="true"/>
+        </listOfProducts>
+        <kineticLaw>
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <times/>
+              <ci> compartment </ci>
+              <ci> k1 </ci>
+              <ci> S1 </ci>
+            </apply>
+          </math>
+        </kineticLaw>
+      </reaction>
+    </listOfReactions>
+    <listOfEvents>
+      <event id="event1" name="event1" useValuesFromTriggerTime="true">
+        <trigger initialValue="true" persistent="true">
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <lt/>
+              <ci> S1 </ci>
+              <cn> 0.1 </cn>
+            </apply>
+          </math>
+        </trigger>
+        <listOfEventAssignments>
+          <eventAssignment variable="S1">
+            <math xmlns="http://www.w3.org/1998/Math/MathML">
+              <cn type="integer"> 1 </cn>
+            </math>
+          </eventAssignment>
+        </listOfEventAssignments>
+      </event>
+      <event id="event2" name="event2" useValuesFromTriggerTime="true">
+        <trigger initialValue="true" persistent="true">
+          <math xmlns="http://www.w3.org/1998/Math/MathML">
+            <apply>
+              <gt/>
+              <ci> S2 </ci>
+              <cn> 0.5 </cn>
+            </apply>
+          </math>
+        </trigger>
+        <listOfEventAssignments>
+          <eventAssignment variable="S2">
+            <math xmlns="http://www.w3.org/1998/Math/MathML">
+              <cn type="integer"> 0 </cn>
+            </math>
+          </eventAssignment>
+        </listOfEventAssignments>
+      </event>
+    </listOfEvents>
+  </model>
+</sbml>
diff --git a/test/data/unsupported.sbml b/test/data/unsupported.sbml
@@ -1,3 +1,3 @@
 a
-listOfEvents
+</delay>
 b
diff --git a/test/events.jl b/test/events.jl
@@ -0,0 +1,24 @@
+myread(fn) = readSBML(fn, doc -> begin
+    set_level_and_version(3, 2)(doc)
+    convert_simplify_math(doc)
+end)
+
+# 1 single trigger, single affect
+fn = "data/00026-sbml-l3v2.xml"
+m = myread(fn)
+@named sys = ODESystem(m)
+ssys = structural_simplify(sys)
+@test length(ModelingToolkit.get_continuous_events(ssys)) == 1
+prob = ODEProblem(ssys, [], (0, 5.0); saveat = 0:0.1:5.0)
+sol = solve(prob, Tsit5())
+@test sol.destats.ncondition > 0
+
+# multiple events 
+fn = "data/00041-sbml-l3v2.xml"
+m = myread(fn)
+@named sys = ODESystem(m)
+ssys = structural_simplify(sys)
+@test length(ModelingToolkit.get_continuous_events(ssys)) == 2
+prob = ODEProblem(ssys, [], (0, 5.0); saveat = 0:0.1:5.0)
+sol = solve(prob, Tsit5())
+@test sol.destats.ncondition > 0
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -8,4 +8,5 @@ using OrdinaryDiffEq
 @testset "SBMLToolkit.jl" begin
     @testset "Model to MTK conversions" begin include("reactionsystem.jl") end
     @testset "Rules" begin include("rules.jl") end
+    @testset "Events" begin include("events.jl") end
 end