feat: initial implementation of new DiscreteSystem

Author: AayushSabharwal

src/ModelingToolkit.jl

@@ -147,6 +147,8 @@ include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
 include("systems/diffeqs/basic_transformations.jl")
 
+include("systems/discrete_system/discrete_system.jl")
+
 include("systems/jumps/jumpsystem.jl")
 
 include("systems/optimization/constraints_system.jl")
@@ -209,6 +211,7 @@ export ODESystem,
 export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDEProblemExpr
 export SystemStructure
+export DiscreteSystem, DiscreteProblem, DiscreteFunction, DiscreteFunctionExpr
 export JumpSystem
 export ODEProblem, SDEProblem
 export NonlinearFunction, NonlinearFunctionExpr

src/clock.jl

@@ -146,3 +146,8 @@ Base.hash(c::SolverStepClock, seed::UInt) = seed ⊻ 0x953d7b9a18874b91
 function Base.:(==)(c1::SolverStepClock, c2::SolverStepClock)
     ((c1.t === nothing || c2.t === nothing) || isequal(c1.t, c2.t))
 end
+
+struct IntegerSequence <: AbstractClock
+    t::Union{Nothing, Symbolic}
+    IntegerSequence(t::Union{Num, Symbolic}) = new(value(t))
+end

src/discretedomain.jl

@@ -28,7 +28,7 @@ struct Shift <: Operator
     Shift(t, steps = 1) = new(value(t), steps)
 end
 Shift(steps::Int) = new(nothing, steps)
-normalize_to_differential(s::Shift) = Differential(s.t)^s.steps
+normalize_to_differential(s::Shift) = Differential(s.t)^abs(s.steps)
 function (D::Shift)(x, allow_zero = false)
     !allow_zero && D.steps == 0 && return x
     Term{symtype(x)}(D, Any[x])
@@ -168,6 +168,7 @@ struct ShiftIndex
     steps::Int
     ShiftIndex(clock::TimeDomain = Inferred(), steps::Int = 0) = new(clock, steps)
     ShiftIndex(t::Num, dt::Real, steps::Int = 0) = new(Clock(t, dt), steps)
+    ShiftIndex(t::Num, steps::Int = 0) = new(IntegerSequence(t), steps)
 end
 
 function (xn::Num)(k::ShiftIndex)

src/systems/diffeqs/abstractodesystem.jl

@@ -812,7 +812,8 @@ function get_u0_p(sys,
     u0, p, defs
 end
 
-function get_u0(sys, u0map, parammap = nothing; symbolic_u0 = false)
+function get_u0(
+        sys, u0map, parammap = nothing; symbolic_u0 = false, toterm = default_toterm)
     dvs = unknowns(sys)
     ps = parameters(sys)
     defs = defaults(sys)
@@ -821,9 +822,10 @@ function get_u0(sys, u0map, parammap = nothing; symbolic_u0 = false)
     end
     defs = mergedefaults(defs, u0map, dvs)
     if symbolic_u0
-        u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = false, use_union = false)
+        u0 = varmap_to_vars(
+            u0map, dvs; defaults = defs, tofloat = false, use_union = false, toterm)
     else
-        u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true)
+        u0 = varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true, toterm)
     end
     return u0, defs
 end