Merge pull request #135 from AayushSabharwal/as/mtkv10

chore: update to ModelingToolkit@10

Author: ChrisRackauckas

Project.toml

@@ -16,7 +16,7 @@ Aqua = "0.8"
 EzXML = "1.1"
 MathML = "0.1.14"
 Memoize = "0.4.2"
-ModelingToolkit = "9"
+ModelingToolkit = "10"
 OrdinaryDiffEq = "6.56"
 SafeTestsets = "0.1"
 Setfield = "1"

README.md

@@ -69,7 +69,7 @@ The models directory contains a few CellML model examples. Let's start with a si
   ml = CellModel("models/lorenz.cellml.xml")
 ```
 
-Now, `ml` is a `CellModel` structure that contains both a list of the loaded XML files and their components (accessible as `ml.doc`) and a ModelingToolkit `ODESystem` that defines variables and dynamics and can be accessed as `getsys(ml)`.
+Now, `ml` is a `CellModel` structure that contains both a list of the loaded XML files and their components (accessible as `ml.doc`) and a ModelingToolkit `System` that defines variables and dynamics and can be accessed as `getsys(ml)`.
 
 The next step is to convert `ml` into an `ODEProblem`, ready to be solved.
 

src/CellMLToolkit.jl

@@ -5,8 +5,8 @@ using MathML: extract_mathml, parse_node
 using Memoize: @memoize
 using SymbolicUtils: operation
 using ModelingToolkit: ModelingToolkit, @parameters, @variables, Differential,
-                       Equation, ODEProblem, ODESystem,
-                       Symbolics, equations, parameters, structural_simplify,
+                       Equation, ODEProblem, System,
+                       Symbolics, equations, parameters, mtkcompile,
                        substitute, unknowns
 using Setfield: @set!
 
@@ -50,7 +50,7 @@ import ModelingToolkit.ODEProblem
 function ODEProblem(ml::CellModel, tspan;
         jac = false, level = 1, p = last.(list_params(ml)),
         u0 = last.(list_states(ml)))
-    ODEProblem(ml.sys, u0, tspan, p; jac = jac)
+    ODEProblem(ml.sys, Pair[unknowns(ml.sys) .=> u0; parameters(ml.sys) .=> p], tspan; jac = jac)
 end
 
 function update_list!(l, sym, val)

src/analysis.jl

@@ -7,7 +7,7 @@ function simplify_systems(systems)
         sys = last(x)
         print("simplifying $(nameof(sys))")
         try
-            sys = structural_simplify(sys)
+            sys = mtkcompile(sys)
             k = max(1, 50 - length(string(nameof(sys))))
             printstyled(repeat(" ", k) * "OK!"; color = :green)
             println()

src/components.jl

@@ -168,7 +168,7 @@ end
 
 """
     pre_substitution generates the substitution rules to be applied to
-    individual systems before structural_simplify merges them together
+    individual systems before mtkcompile merges them together
 
     it morphes variables and parameters to their correct form for ModelingToolkit
     based on the global CellML doc information
@@ -193,7 +193,7 @@ end
 
 """
     post_substitution generates the substitution rules to be applied to
-    the merged system after structural_simplify is applied
+    the merged system after mtkcompile is applied
 
     if changes the names of the independent variable (iv) in each system
     to the global iv name
@@ -212,7 +212,7 @@ substitute_eqs(eqs, s) = [substitute(eq.lhs, s) ~ substitute(eq.rhs, s) for eq i
 
 """
     process_components is the main entry point
-    it processes an doc document and returns a merged ODESystem
+    it processes an doc document and returns a merged System
 
     use simplify=false only for debugging purposes!
 """
@@ -224,13 +224,13 @@ function process_components(doc::Document; simplify = true)
     systems = subsystems(doc, class)
     post_sub = post_substitution(doc, systems)
 
-    sys = ODESystem(translate_connections(doc, systems, class),
+    sys = System(Vector{Equation}(translate_connections(doc, systems, class)),
         get_ivₚ(doc),
         systems = collect(values(systems)),
         name = gensym(:cellml))
 
     if simplify
-        sys = structural_simplify(sys)
+        sys = mtkcompile(sys)
         @set! sys.eqs = substitute_eqs(equations(sys), post_sub)
 
         # Defaults need to be set after simplifying as otherwise parameters and
@@ -250,12 +250,12 @@ end
     class is the output of classify_variables
 """
 function subsystems(doc::Document, class)
-    Dict{Symbol, ODESystem}(to_symbol(comp) => process_component(doc, comp, class)
+    Dict{Symbol, System}(to_symbol(comp) => process_component(doc, comp, class)
     for comp in components(doc))
 end
 
 """
-    process_component converts a single CellML component to an ODESystem
+    process_component converts a single CellML component to an System
 
     comp in the name of the component
     class is the output of classify_variables
@@ -270,7 +270,7 @@ function process_component(doc::Document, comp, class)
         eqs = vcat(parse_node.(math)...)
         eqs = substitute_eqs(eqs, pre_sub)
         sub_rhs = remove_rhs_diff(eqs)
-        eqs = [eq.lhs ~ substitute(eq.rhs, sub_rhs) for eq in eqs]
+        eqs = Equation[eq.lhs ~ substitute(eq.rhs, sub_rhs) for eq in eqs]
     end
 
     ivₚ = get_ivₚ(doc)
@@ -280,7 +280,7 @@ function process_component(doc::Document, comp, class)
     states = [last(x)
               for x in values(pre_sub)
               if !ModelingToolkit.isparameter(last(x)) && !isequal(last(x), ivₚ)]
-    ODESystem(eqs, ivₚ, states, ps; name = to_symbol(comp))
+    System(eqs, ivₚ, states, ps; name = to_symbol(comp))
 end
 
 ###############################################################################

src/structures.jl

@@ -32,5 +32,5 @@ connections(doc::Document) = doc.conns
 
 struct CellModel
     doc::Document
-    sys::ODESystem
+    sys::System
 end

test/beeler.jl

@@ -24,7 +24,7 @@ err = sum(abs.(V1 .- V2)) / length(V1)
 # Ensure defaults are set and that generating an `ODEProblem` directly from the
 # `ODESystem` is equivalent to doing so from a `CellModel`
 @test length(ModelingToolkit.defaults(ml.sys)) > 0
-sys_prob = ODEProblem(ml.sys; tspan = (0, 10000.0))
+sys_prob = ODEProblem(ml.sys, nothing, (0, 10000.0))
 sol3 = solve(prob, Euler(), dt = 0.01, saveat = 1.0)
 V3 = map(x -> x[1], sol3.u)
 err2 = sum(abs.(V1 .- V3)) / length(V1)

test/noble_1962.jl

@@ -15,7 +15,7 @@ err1 = sum(abs.(V1 .- V2)) / length(V1)
 # Ensure defaults are set and that generating an `ODEProblem` directly from the
 # `ODESystem` is equivalent to doing so from a `CellModel`
 @test length(ModelingToolkit.defaults(ml.sys)) > 0
-sys_prob = ODEProblem(ml.sys; tspan = (0, 10000.0))
+sys_prob = ODEProblem(ml.sys, nothing, (0, 10000.0))
 sol3 = solve(prob, Euler(), dt = 0.01, saveat = 1.0)
 V3 = map(x -> x[2], sol3.u)
 err2 = sum(abs.(V1 .- V3)) / length(V1)