Algebraic equations are delayed one step too much in DiscreteSystems

[baggepinnen]

In the following example, the variables x and y are equal, but y is delayed one step compared to x in the solution.

# test discrete integrator
@variables t x(t)=0 y(t)=0 [output=true]
Dₜ = Difference(t; dt=dt)
eqs = [
    Dₜ(x) ~ x + 1
    y ~ x
]
@named sys = DiscreteSystem(eqs, t)
prob = DiscreteProblem(sys, Pair[], (0.0, 10.0))
sol = solve(prob, FunctionMap())
yind = findfirst(isequal(y),states(sys)) # due to https://github.com/SciML/ModelingToolkit.jl/issues/1329
xind = findfirst(isequal(x),states(sys))
@test Array(sol)[yind,:] == Array(sol)[xind,:] # fails
@test Array(sol)[yind,2:end] == Array(sol)[xind,1:end-1] # passes, but shouldn't

[baggepinnen]

There are multiple problems here

1. If the system is not simplified, algebraic equations are not moved into observed and they are all passed to DiscreteFunction which assumes all equations are on the form x(t+1) ~ f(x(t)), causing the delay.
2. Simplification of the discrete system does not do the right thing because initialize_system_structure assumes an ODESystem, and only looks for Differential variables. The call vars!(vars, eq.rhs) defaults to vars!(vars, eq.rhs; op=Differential), but adding another vars!(vars, eq.rhs; op=Difference) does not solve the problem.

@YingboMa could there be some modifications to algorithms that should work for both differential and difference equations such that they operate on some common supertype of Differential and Difference? And similarly calls to functions isdifferential, isdifference would call isdiff(x) = isdifferential(x) || isdifference(x)?

[ChrisRackauckas]

Old discrete was removed.