feat: support observed function generation for DDEs

Author: AayushSabharwal

src/systems/abstractsystem.jl

@@ -805,6 +805,10 @@ function SymbolicIndexingInterface.observed(
         return let _fn = _fn
             fn1(u, p, t) = _fn(u, p, t)
             fn1(u, p::MTKParameters, t) = _fn(u, p..., t)
+
+            # DDEs
+            fn1(u, histfn, p, t) = _fn(u, histfn, p, t)
+            fn1(u, histfn, p::MTKParameters, t) = _fn(u, histfn, p..., t)
             fn1
         end
     else

src/systems/diffeqs/odesystem.jl

@@ -475,13 +475,22 @@ function build_explicit_observed_function(sys, ts;
     end
     ts = map(t -> substitute(t, subs), ts)
     obsexprs = []
+
+    histfn_name = gensym(:hist)
+    histfn = only(@variables ($histfn_name)(..)::Vector{Real})
+
+    delays = Set()
+    _vars_util = Set()
     for i in 1:maxidx
         eq = obs[i]
         lhs = eq.lhs
         rhs = eq.rhs
+        populate_delays(delays, obsexprs, histfn, sys, rhs)
         push!(obsexprs, lhs ← rhs)
     end
 
+    populate_delays(delays, obsexprs, histfn, sys, ts)
+
     if inputs !== nothing
         ps = setdiff(ps, inputs) # Inputs have been converted to parameters by io_preprocessing, remove those from the parameter list
     end
@@ -497,12 +506,17 @@ function build_explicit_observed_function(sys, ts;
         ps = (DestructuredArgs(unwrap.(ps), inbounds = !checkbounds),)
     end
     dvs = DestructuredArgs(unknowns(sys), inbounds = !checkbounds)
+    if !isempty(delays)
+        dvs = (dvs, histfn)
+    else
+        dvs = (dvs,)
+    end
     if inputs === nothing
-        args = param_only ? [ps..., ivs...] : [dvs, ps..., ivs...]
+        args = param_only ? [ps..., ivs...] : [dvs..., ps..., ivs...]
     else
         inputs = unwrap.(inputs)
         ipts = DestructuredArgs(inputs, inbounds = !checkbounds)
-        args = param_only ? [ipts, ps..., ivs...] : [dvs, ipts, ps..., ivs...]
+        args = param_only ? [ipts, ps..., ivs...] : [dvs..., ipts, ps..., ivs...]
     end
     pre = get_postprocess_fbody(sys)
 
@@ -538,6 +552,20 @@ function build_explicit_observed_function(sys, ts;
     end
 end
 
+function populate_delays(delays::Set, obsexprs, histfn, sys, sym)
+    _vars_util = vars(sym)
+    for v in _vars_util
+        v in delays && continue
+        iscall(v) && issym(operation(v)) && (args = arguments(v); length(args) == 1) &&
+            iscall(only(args)) || continue
+
+        idx = variable_index(sys, operation(v)(get_iv(sys)))
+        idx === nothing && error("Delay term $v is not an unknown in the system")
+        push!(delays, v)
+        push!(obsexprs, v ← histfn(only(args))[idx])
+    end
+end
+
 function _eq_unordered(a, b)
     length(a) === length(b) || return false
     n = length(a)