Merge pull request #3698 from AayushSabharwal/as/concrete-getu

[v9] hotfix: fix new implementation of `concrete_getu`

Author: AayushSabharwal

.github/workflows/Downstream.yml

@@ -1,7 +1,7 @@
 name: IntegrationTest
 on:
   push:
-    branches: [master]
+    branches: [master, 'backport-v9']
     tags: [v*]
   pull_request:
     paths-ignore:

.github/workflows/ReleaseTest.yml

@@ -1,7 +1,7 @@
 name: ReleaseTest
 on:
   push:
-    branches: [master]
+    branches: [master, 'backport-v9']
     tags: [v*]
   pull_request:
     paths-ignore:

.github/workflows/Tests.yml

@@ -5,11 +5,13 @@ on:
     branches:
       - master
       - 'release-'
+      - 'backport-v9'
     paths-ignore:
       - 'docs/**'
   push:
     branches:
       - master
+      - 'backport-v9'
     paths-ignore:
       - 'docs/**'
 

src/systems/diffeqs/odesystem.jl

@@ -476,6 +476,8 @@ Generates a function that computes the observed value(s) `ts` in the system `sys
 - `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist.
 - `mkarray`: only used if the output is an array (that is, `!isscalar(ts)`  and `ts` is not a tuple, in which case the result will always be a tuple). Called as `mkarray(ts, output_type)` where `ts` are the expressions to put in the array and `output_type` is the argument of the same name passed to build_explicit_observed_function.
 - `cse = true`: Whether to use Common Subexpression Elimination (CSE) to generate a more efficient function.
+- `wrap_delays = is_dde(sys)`: Whether to add an argument for the history function and use
+  it to calculate all delayed variables.
 
 ## Returns
 
@@ -514,7 +516,8 @@ function build_explicit_observed_function(sys, ts;
         op = Operator,
         throw = true,
         cse = true,
-        mkarray = nothing)
+        mkarray = nothing,
+        wrap_delays = is_dde(sys))
     is_tuple = ts isa Tuple
     if is_tuple
         ts = collect(ts)
@@ -600,14 +603,15 @@ function build_explicit_observed_function(sys, ts;
     p_end = length(dvs) + length(inputs) + length(ps)
     fns = build_function_wrapper(
         sys, ts, args...; p_start, p_end, filter_observed = obsfilter,
-        output_type, mkarray, try_namespaced = true, expression = Val{true}, cse)
+        output_type, mkarray, try_namespaced = true, expression = Val{true}, cse,
+        wrap_delays)
     if fns isa Tuple
         if expression
             return return_inplace ? fns : fns[1]
         end
         oop, iip = eval_or_rgf.(fns; eval_expression, eval_module)
         f = GeneratedFunctionWrapper{(
-            p_start + is_dde(sys), length(args) - length(ps) + 1 + is_dde(sys), is_split(sys))}(
+            p_start + wrap_delays, length(args) - length(ps) + 1 + wrap_delays, is_split(sys))}(
             oop, iip)
         return return_inplace ? (f, f) : f
     else
@@ -616,7 +620,7 @@ function build_explicit_observed_function(sys, ts;
         end
         f = eval_or_rgf(fns; eval_expression, eval_module)
         f = GeneratedFunctionWrapper{(
-            p_start + is_dde(sys), length(args) - length(ps) + 1 + is_dde(sys), is_split(sys))}(
+            p_start + wrap_delays, length(args) - length(ps) + 1 + wrap_delays, is_split(sys))}(
             f, nothing)
         return f
     end

src/systems/problem_utils.jl

@@ -631,7 +631,9 @@ end
 ObservedWrapper{TD}(f::F) where {TD, F} = ObservedWrapper{TD, F}(f)
 
 function (ow::ObservedWrapper{true})(prob)
-    ow.f(state_values(prob), parameter_values(prob), current_time(prob))
+    # Edge case for steady state problems
+    t = applicable(current_time, prob) ? current_time(prob) : Inf
+    ow.f(state_values(prob), parameter_values(prob), t)
 end
 
 function (ow::ObservedWrapper{false})(prob)
@@ -649,7 +651,7 @@ function. It does NOT work for solutions.
 """
 Base.@nospecializeinfer function concrete_getu(indp, syms::AbstractVector)
     @nospecialize
-    obsfn = SymbolicIndexingInterface.observed(indp, syms)
+    obsfn = build_explicit_observed_function(indp, syms; wrap_delays = false)
     return ObservedWrapper{is_time_dependent(indp)}(obsfn)
 end
 

src/utils.jl

@@ -1083,7 +1083,7 @@ Keyword arguments:
   `available_vars` will not be searched for in the observed equations.
 """
 function observed_equations_used_by(sys::AbstractSystem, exprs;
-        involved_vars = vars(exprs; op = Union{Shift, Differential}), obs = observed(sys), available_vars = [])
+        involved_vars = vars(exprs; op = Union{Shift, Differential, Initial}), obs = observed(sys), available_vars = [])
     obsvars = getproperty.(obs, :lhs)
     graph = observed_dependency_graph(obs)
     if !(available_vars isa Set)

test/odesystem.jl

@@ -5,10 +5,14 @@ using OrdinaryDiffEq, Sundials
 using DiffEqBase, SparseArrays
 using StaticArrays
 using Test
-using SymbolicUtils: issym
+using SymbolicUtils.Code
+using SymbolicUtils: Sym, issym
 using ForwardDiff
 using ModelingToolkit: value
 using ModelingToolkit: t_nounits as t, D_nounits as D
+using Symbolics
+using Symbolics: unwrap
+using DiffEqBase: isinplace
 
 # Define some variables
 @parameters σ ρ β
@@ -607,13 +611,6 @@ sys = complete(sys)
 @test_throws Any ODEFunction(sys)
 
 @testset "Preface tests" begin
-    using OrdinaryDiffEq
-    using Symbolics
-    using DiffEqBase: isinplace
-    using ModelingToolkit
-    using SymbolicUtils.Code
-    using SymbolicUtils: Sym
-
     c = [0]
     function f(c, du::AbstractVector{Float64}, u::AbstractVector{Float64}, p, t::Float64)
         c .= [c[1] + 1]
@@ -656,7 +653,9 @@ sys = complete(sys)
 
     @named sys = ODESystem(eqs, t, us, ps; defaults = defs, preface = preface)
     sys = complete(sys)
-    prob = ODEProblem(sys, [], (0.0, 1.0))
+    # don't build initializeprob because it will use preface in other functions and
+    # affect `c`
+    prob = ODEProblem(sys, [], (0.0, 1.0); build_initializeprob = false)
     sol = solve(prob, Euler(); dt = 0.1)
 
     @test c[1] == length(sol)