DiscreteProblem without and with `ifelse`

[bolognam]

Describe the bug 🐞

I have two models that I believe are identical but one gives a missing variable error message while the other works as expected.

Expected behavior

The isapprox line at the bottom of each script should return true values across all timesteps.

Minimal Reproducible Example 👇

Working script

using ModelingToolkit
using ModelingToolkit: t_nounits as t
using OrdinaryDiffEq: solve, FunctionMap
using OrdinaryDiffEq
using DelimitedFiles
# Equation:
# var(k) ~ var(k - 1) + var(k + 1) + shk(k)
# But positive shifts aren't allowed, so shift down by the max lead:
# var(k - 1) ~ var(k - 2) + var(k) + shk(k - 1)
# Equivalent to:
# var(k) ~ var(k - 1) - var(k - 2) - shk(k - 1)
maxlag = 1
maxlead = 1
shk_data = [0.0, 1.0, 0.0, 2.0, 0.0, 3.0, 0.0, 4.0, 0.0, 0.0]
var_res = [1.0, 1.0, -1.0, -2.0, -3.0, -1.0, -1.0, 0.0, -3.0, -3.0]
@assert length(var_res) == length(shk_data)
tstart = 0 # Time at `var_res[1+maxlag]` and `shk_data[1+maxlag]`
tend = length(var_res) - 1 - maxlag - maxlead # Time at `var_res[end-maxlead]` and `shk_data[end-maxlead]`
k = ShiftIndex(t)
getdata(buffer, t) = buffer[Int(t)+maxlag+1]
@register_symbolic getdata(buffer::Vector, t)
function System(; name, buffer)
    pars = @parameters begin
        shk::Vector{Float64} = buffer
    end
    vars = @variables begin
        var(t)
        time(t)
    end
    eqs = [
        time(k) ~ time(k-1) + 1
        var(k - 1) ~ var(k - 2) + var(k) + getdata(shk, time - 1)
    ]
    return DiscreteSystem(eqs, t, vars, pars; name)
end
@named de = System(; buffer = shk_data)
simp = structural_simplify(de)
s = complete(simp)
prob = DiscreteProblem(s, [s.time(k-1) => tstart - 1, s.var(k-2) => 0.0, s.var(k-1) => var_res[1]], (tstart, tend + maxlead))
sol = solve(prob, FunctionMap())
res = [var_res[1:maxlag]; sol[s.var]]
@show isapprox.(res, var_res; atol = 1e-9)
res

Non-working script (produces error below)

using ModelingToolkit
using ModelingToolkit: t_nounits as t
using OrdinaryDiffEq: solve, FunctionMap
using OrdinaryDiffEq
using DelimitedFiles

# Equation:
# var(k) ~ var(k - 1) + var(k + 1) + shk(k)
# But positive shifts aren't allowed, so shift down by the max lead:
# var(k - 1) ~ var(k - 2) + var(k) + shk(k - 1)
# Equivalent to:
# var(k) ~ var(k - 1) - var(k - 2) - shk(k - 1)

maxlag = 1
maxlead = 1

shk_data = [0.0, 1.0, 0.0, 2.0, 0.0, 3.0, 0.0, 4.0, 0.0, 0.0]
var_res = [1.0, 1.0, -1.0, -2.0, -3.0, -1.0, -1.0, 0.0, -3.0, -3.0]
@assert length(var_res) == length(shk_data)
#var_res = [0; var_res]

tstart = 0 # Time at `var_res[1+maxlag]` and `shk_data[1+maxlag]`
tend = length(var_res) - 1 - maxlag - maxlead # Time at `var_res[end-maxlead]` and `shk_data[end-maxlead]`
k = ShiftIndex(t)

getdata(buffer, t) = buffer[@show(Int(t))+maxlag+1]
@register_symbolic getdata(buffer::Vector, t)

function System(; name, buffer)

    pars = @parameters begin
        shk::Vector{Float64} = buffer
    end

    vars = @variables begin
        var(t)
        time(t)
    end

    eqs = [
        time(k) ~ time(k-1) + 1
        var(k - 1) ~ ifelse(time - 1 < tstart,
            getdata(var_res, time - 1),
            ifelse(time - 1 > tend,
                getdata(var_res, time - 1),
                var(k - 2) + var(k) + getdata(shk, time - 1),
            ),
        )
    ]

    return DiscreteSystem(eqs, t, vars, pars; name)

end

@named de = System(; buffer = shk_data)
simp = structural_simplify(de)
s = complete(simp)
prob = DiscreteProblem(s, [s.time(k-1) => tstart - 1, s.var(k-1) => var_res[1]], (tstart, tend + maxlead))

sol = solve(prob, FunctionMap())

res = [var_res[(1:maxlag) .+ 0]; sol[s.var]]
@show isapprox.(res, var_res[1:end]; atol = 1e-9)
res

Error & Stacktrace ⚠️

ERROR: LoadError: Initial condition for Shift(t, -1)(var(t)) not provided.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] process_DiscreteProblem(constructor::Type{…}, sys::DiscreteSystem, u0map::Vector{…}, parammap::SciMLBase.NullParameters; linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, kwargs::@Kwargs{…})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/w72sG/src/systems/discrete_system/discrete_system.jl:246
 [3] process_DiscreteProblem
   @ ~/.julia/packages/ModelingToolkit/w72sG/src/systems/discrete_system/discrete_system.jl:223 [inlined]
 [4] DiscreteProblem(sys::DiscreteSystem, u0map::Vector{…}, tspan::Tuple{…}, parammap::SciMLBase.NullParameters; eval_module::Module, eval_expression::Bool, use_union::Bool, kwargs::@Kwargs{})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/w72sG/src/systems/discrete_system/discrete_system.jl:286
 [5] DiscreteProblem
   @ ~/.julia/packages/ModelingToolkit/w72sG/src/systems/discrete_system/discrete_system.jl:270 [inlined]
 [6] DiscreteProblem(sys::DiscreteSystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Int64, Int64})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/w72sG/src/systems/discrete_system/discrete_system.jl:270
 [7] top-level scope
   @ ~/Documents/glcs/Bank of Canada/src/BankOfCanada.jl/scripts/test.jl:58
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:489
 [9] top-level scope
   @ REPL[4]:1
in expression starting at /Users/mbologna/Documents/glcs/Bank of Canada/src/BankOfCanada.jl/scripts/test.jl:58
Some type information was truncated. Use `show(err)` to see complete types.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [8bb1440f] DelimitedFiles v1.9.1
  [961ee093] ModelingToolkit v9.13.0
  [1dea7af3] OrdinaryDiffEq v6.76.0

• Output of versioninfo()

Julia Version 1.10.2
Commit bd47eca2c8a (2024-03-01 10:14 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
Threads: 8 default, 0 interactive, 4 GC (on 8 virtual cores)
Environment:
  JULIA_NUM_THREADS = auto
  JULIA_COPY_STACKS = yes
  JULIA_PKG_SERVER = juliahub.com
  JULIA_PKG_USE_CLI_GIT = true