Set known defaults for a CellModels ODESystem

CellML files contain `initial_value` fields specifying initial values
for parameters and states in the system. Retrieving these for the
`ODESystem` was already possible, but had to be done separately or as
part of the `ODEProblem` constructor defined for `CellModel`s.

Defaults are only set for simplified `ODESystem`s. Unsimplified
`ODESystem`s may contain parameters and states for which no defaults are
available resulting in errors trying to construct the `ODESystem`. These
errors will still occur when using the result of calling
`process_components` with `simplify = false`, but as noted in the code
that should only be done for debugging purposes and this functionality
is not exposed at the higher API levels.

Fixes #56.

Author: bauglir

src/components.jl

@@ -234,6 +234,11 @@ function process_components(doc::Document; simplify = true)
     if simplify
         sys = structural_simplify(sys)
         @set! sys.eqs = substitute_eqs(equations(sys), post_sub)
+
+        # Defaults need to be set after simplifying as otherwise parameters and
+        # states for which no defaults are available may still be present in
+        # the system
+        @set! sys.defaults = Dict(find_list_value(doc, vcat(parameters(sys), states(sys))))
     end
 
     return sys

test/beeler.jl

@@ -17,6 +17,15 @@ V2 = map(x -> x[1], sol2.u)
 err = sum(abs.(V1 .- V2)) / length(V1)
 @test err < 0.1
 
+# 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))
+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 err2 ≈ 0
+
 # prob = ODEProblem(ml, (0,10000.0); jac=true)
 # sol3 = solve(prob, TRBDF2(), dtmax=0.5, saveat=1.0)
 # V3 = map(x -> x[1], sol2.u)

test/noble_1962.jl

@@ -6,5 +6,14 @@ sol1 = solve(prob, Euler(), dt = 0.01, saveat = 1.0)
 sol2 = solve(prob, TRBDF2(), dtmax = 0.5, saveat = 1.0)
 V1 = map(x -> x[2], sol1.u)
 V2 = map(x -> x[2], sol2.u)
-err = sum(abs.(V1 .- V2)) / length(V1)
-@test err < 1.0
+err1 = sum(abs.(V1 .- V2)) / length(V1)
+@test err1 < 1.0
+
+# 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))
+sol3 = solve(prob, Euler(), dt = 0.01, saveat = 1.0)
+V3 = map(x -> x[2], sol3.u)
+err2 = sum(abs.(V1 .- V3)) / length(V1)
+@test err2 ≈ 0