support for observables

Author: Brad Carman

src/systems/abstractsystem.jl

@@ -767,25 +767,32 @@ function toexpr(sys::AbstractSystem)
     psname = gensym(:ps)
     ps = parameters(sys)
     push_vars!(stmt, psname, Symbol("@parameters"), ps)
+    obs = observed(sys)
+    obsvars = [o.lhs for o in obs]
+    obsvarsname = gensym(:obs)
+    push_vars!(stmt, obsvarsname, Symbol("@variables"), obsvars)
 
     var2name = Dict{Any, Symbol}()
-    for v in Iterators.flatten((sts, ps))
+    for v in Iterators.flatten((sts, ps, obsvars))
         var2name[v] = getname(v)
     end
 
-    eqs_name = push_eqs!(stmt, equations(sys), var2name)
+    eqs_name = push_eqs!(stmt, full_equations(sys), var2name)
     defs_name = push_defaults!(stmt, defaults(sys), var2name)
+    obs_name = push_eqs!(stmt, obs, var2name)
 
     if sys isa ODESystem
         iv = get_iv(sys)
         ivname = gensym(:iv)
         push!(stmt, :($ivname = (@variables $(getname(iv)))[1]))
         push!(stmt,
             :($ODESystem($eqs_name, $ivname, $stsname, $psname; defaults = $defs_name,
+                observed = $obs_name,
                 name = $name, checks = false)))
     elseif sys isa NonlinearSystem
         push!(stmt,
             :($NonlinearSystem($eqs_name, $stsname, $psname; defaults = $defs_name,
+                observed = $obs_name,
                 name = $name, checks = false)))
     end
 

src/systems/diffeqs/abstractodesystem.jl

@@ -516,6 +516,7 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
     sparse = false, simplify = false,
     steady_state = false,
     sparsity = false,
+    observedfun_exp = nothing,
     kwargs...) where {iip}
     f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)
 
@@ -567,7 +568,8 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
             syms = $(Symbol.(states(sys))),
             indepsym = $(QuoteNode(Symbol(get_iv(sys)))),
             paramsyms = $(Symbol.(parameters(sys))),
-            sparsity = $(sparsity ? jacobian_sparsity(sys) : nothing))
+            sparsity = $(sparsity ? jacobian_sparsity(sys) : nothing),
+            observed = $observedfun_exp)
     end
     !linenumbers ? striplines(ex) : ex
 end

test/serialization.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, SciMLBase, Serialization
+using ModelingToolkit, SciMLBase, Serialization, OrdinaryDiffEq
 
 @parameters t
 @variables x(t)
@@ -28,3 +28,41 @@ write(io, rc_model)
 str = String(take!(io))
 sys = include_string(@__MODULE__, str)
 @test sys == flatten(rc_model) # this actually kind of works, but the variables would have different identities.
+
+# check answer
+ss = structural_simplify(rc_model)
+all_obs = [o.lhs for o in observed(ss)]
+prob = ODEProblem(ss, [], (0, 0.1))
+sol = solve(prob, ImplicitEuler())
+
+## Check ODESystem with Observables ----------
+ss_exp = ModelingToolkit.toexpr(ss)
+ss_ = eval(ss_exp)
+prob_ = ODEProblem(ss_, [], (0, 0.1))
+sol_ = solve(prob_, ImplicitEuler())
+@test sol[all_obs] == sol_[all_obs]
+
+## Check ODEProblemExpr with Observables -----------
+
+# build the observable function expression
+obs_exps = []
+for var in all_obs
+    f = ModelingToolkit.build_explicit_observed_function(ss, var; expression = true)
+    sym = ModelingToolkit.getname(var) |> string
+    ex = :(if name == Symbol($sym)
+        return $f(u0, p, t)
+    end)
+    push!(obs_exps, ex)
+end
+# observedfun expression for ODEFunctionExpr
+observedfun_exp = :(function (var, u0, p, t)
+    name = ModelingToolkit.getname(var)
+    $(obs_exps...)
+end)
+
+# ODEProblemExpr with observedfun_exp included
+probexpr = ODEProblemExpr{true}(ss, [], (0, 0.1); observedfun_exp);
+prob_obs = eval(probexpr)
+sol_obs = solve(prob_obs, ImplicitEuler())
+
+@test sol_obs[all_obs] == sol[all_obs]