Merge pull request #1832 from SciML/myb/obsfix

Fix observed function generation

Author: YingboMa

src/structural_transformation/StructuralTransformations.jl

@@ -22,7 +22,7 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
                        get_postprocess_fbody, vars!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
                        invalidate_cache!, Substitutions, get_or_construct_tearing_state,
-                       AliasGraph
+                       AliasGraph, filter_kwargs
 
 using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview

src/structural_transformation/codegen.jl

@@ -384,7 +384,6 @@ function build_observed_function(state, ts, var_eq_matching, var_sccs,
 
     fullvars = state.fullvars
     s = state.structure
-    graph = s.graph
     solver_states = fullvars[is_solver_state_idxs]
     algvars = fullvars[.!is_solver_state_idxs]
 
@@ -525,6 +524,7 @@ function ODAEProblem{iip}(sys,
     has_difference = any(isdifferenceeq, eqs)
     cbs = process_events(sys; callback, has_difference, kwargs...)
 
+    kwargs = filter_kwargs(kwargs)
     if cbs === nothing
         ODEProblem{iip}(fun, u0, tspan, p; kwargs...)
     else

src/structural_transformation/symbolics_tearing.jl

@@ -241,6 +241,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
     #removed_vars = Int[]
     removed_obs = Int[]
     diff_to_var = invview(var_to_diff)
+    dummy_sub = Dict()
     for var in 1:length(fullvars)
         dv = var_to_diff[var]
         dv === nothing && continue
@@ -253,7 +254,8 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
                 push!(removed_obs, idx)
             end
             for eq in 𝑑neighbors(graph, dv)
-                neweqs[eq] = substitute(neweqs[eq], fullvars[dv] => v_t)
+                dummy_sub[dd] = v_t
+                neweqs[eq] = substitute(neweqs[eq], dd => v_t)
             end
             fullvars[dv] = v_t
             # update the structural information
@@ -590,11 +592,17 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
         oldobs[idx] = (lhs ~ neweqs[eqidx].rhs)
     end
     deleteat!(oldobs, sort!(removed_obs))
-    @set! sys.observed = [oldobs; subeqs]
     @set! sys.substitutions = Substitutions(subeqs, deps)
+
+    obs_sub = dummy_sub
+    for eq in equations(sys)
+        isdiffeq(eq) || continue
+        obs_sub[eq.lhs] = eq.rhs
+    end
+    obs = substitute.([oldobs; subeqs], (obs_sub,))
+    @set! sys.observed = obs
     @set! state.sys = sys
     @set! sys.tearing_state = state
-
     return invalidate_cache!(sys)
 end
 

src/systems/diffeqs/odesystem.jl

@@ -323,25 +323,10 @@ function build_explicit_observed_function(sys, ts;
     end
     ts = map(t -> substitute(t, subs), ts)
     obsexprs = []
-    eqs_cache = Ref{Any}(nothing)
     for i in 1:maxidx
         eq = obs[i]
         lhs = eq.lhs
         rhs = eq.rhs
-        vars!(vars, rhs)
-        for v in vars
-            isdifferential(v) || continue
-            if eqs_cache[] === nothing
-                eqs_cache[] = Dict(eq.lhs => eq.rhs for eq in equations(sys))
-            end
-            eqs_dict = eqs_cache[]
-            rhs_diffeq = get(eqs_dict, v, nothing)
-            push!(obsexprs, v ← rhs_diffeq)
-            if rhs === nothing
-                error("The observed variable $(eq.lhs) depends on the differentiated variable $v, but it's not explicit solved. Fix file an issue if you are sure that the system is valid.")
-            end
-        end
-        empty!(vars)
         push!(obsexprs, lhs ← rhs)
     end