Merge pull request #828 from SciML/myb/obs

Populate the `observed` field in ODEFunction

Author: YingboMa

src/systems/diffeqs/abstractodesystem.jl

@@ -162,6 +162,15 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
 
     _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)
 
+    observedfun = let sys = sys, dict = Dict()
+        function generated_observed(obsvar, u, p, t)
+            obs = get!(dict, value(obsvar)) do
+                build_explicit_observed_function(sys, obsvar)
+            end
+            obs(u, p, t)
+        end
+    end
+
     ODEFunction{iip}(
                      f,
                      jac = _jac === nothing ? nothing : _jac,
@@ -170,6 +179,7 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                      jac_prototype = sparse ? similar(get_jac(sys)[],Float64) : nothing,
                      syms = Symbol.(states(sys)),
                      indepsym = Symbol(independent_variable(sys)),
+                     observed = observedfun,
                     )
 end
 

src/systems/diffeqs/odesystem.jl

@@ -230,3 +230,36 @@ function flatten(sys::ODESystem)
 end
 
 ODESystem(eq::Equation, args...; kwargs...) = ODESystem([eq], args...; kwargs...)
+
+"""
+$(SIGNATURES)
+
+Build the observed function assuming the observed equations are all explicit,
+i.e. there are no cycles or dependencies.
+"""
+function build_explicit_observed_function(
+        sys, syms;
+        expression=false,
+        output_type=Array)
+
+    if (isscalar = !(syms isa Vector))
+        syms = [syms]
+    end
+    syms = value.(syms)
+
+    obs = observed(sys)
+    observed_idx = Dict(map(x->x.lhs, obs) .=> 1:length(obs))
+    output = map(sym->obs[observed_idx[sym]].rhs, syms)
+
+    ex = Func(
+        [
+         DestructuredArgs(states(sys))
+         DestructuredArgs(parameters(sys))
+         independent_variable(sys)
+        ],
+        [],
+        isscalar ? output[1] : MakeArray(output, output_type)
+    ) |> toexpr
+
+    expression ? ex : @RuntimeGeneratedFunction(ex)
+end

test/components.jl

@@ -93,6 +93,12 @@ u0 = [
       resistor.v => 0.0
      ]
 prob = ODEProblem(sys, u0, (0, 10.0))
-sol = solve(prob, Rodas5())
+sol = solve(prob, Rodas4())
+
+@test sol[resistor.p.i] == sol[capacitor.p.i]
+@test sol[resistor.n.i] == -sol[capacitor.p.i]
+@test sol[capacitor.n.i] == -sol[capacitor.p.i]
+@test iszero(sol[ground.g.i])
+@test iszero(sol[ground.g.v])
 #using Plots
 #plot(sol)