Merge pull request #2440 from AayushSabharwal/rmdiff

refactor!: remove old discrete system

Author: ChrisRackauckas

docs/pages.jl

@@ -32,7 +32,6 @@ pages = [
         "systems/JumpSystem.md",
         "systems/NonlinearSystem.md",
         "systems/OptimizationSystem.md",
-        "systems/DiscreteSystem.md",
         "systems/PDESystem.md"],
     "comparison.md",
     "internals.md",

docs/src/systems/DiscreteSystem.md

@@ -24,7 +24,7 @@ structural_simplify
 ## Problem Constructors
 
 ```@docs; canonical=false
-DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, u0map, tspan)
+DiscreteProblem(sys::JumpSystem, u0map, tspan)
 ```
 
 ```@docs

docs/src/systems/JumpSystem.md

@@ -150,7 +150,7 @@ include("systems/optimization/modelingtoolkitize.jl")
 include("systems/pde/pdesystem.jl")
 
 include("systems/sparsematrixclil.jl")
-include("systems/discrete_system/discrete_system.jl")
+
 include("systems/unit_check.jl")
 include("systems/validation.jl")
 include("systems/dependency_graphs.jl")
@@ -196,7 +196,7 @@ export NonlinearFunction, NonlinearFunctionExpr
 export NonlinearProblem, BlockNonlinearProblem, NonlinearProblemExpr
 export OptimizationProblem, OptimizationProblemExpr, constraints
 export SteadyStateProblem, SteadyStateProblemExpr
-export JumpProblem, DiscreteProblem
+export JumpProblem
 export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, Flow, Stream, instream
@@ -214,9 +214,6 @@ export SymScope, LocalScope, ParentScope, DelayParentScope, GlobalScope
 export independent_variable, equations, controls,
     observed, structure, full_equations
 export structural_simplify, expand_connections, linearize, linearization_function
-export DiscreteSystem,
-    DiscreteProblem, DiscreteProblemExpr, DiscreteFunction,
-    DiscreteFunctionExpr
 
 export calculate_jacobian, generate_jacobian, generate_function
 export calculate_control_jacobian, generate_control_jacobian
@@ -227,7 +224,6 @@ export calculate_hessian, generate_hessian
 export calculate_massmatrix, generate_diffusion_function
 export stochastic_integral_transform
 export TearingState, StateSelectionState
-export generate_difference_cb
 
 export BipartiteGraph, equation_dependencies, variable_dependencies
 export eqeq_dependencies, varvar_dependencies

src/ModelingToolkit.jl

@@ -42,7 +42,7 @@ end
 has_time_domain(x::Num) = has_time_domain(value(x))
 has_time_domain(x) = false
 
-for op in [Differential, Difference]
+for op in [Differential]
     @eval input_timedomain(::$op, arg = nothing) = Continuous()
     @eval output_timedomain(::$op, arg = nothing) = Continuous()
 end

src/clock.jl

@@ -215,7 +215,7 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     end
     inputs = map(x -> time_varying_as_func(value(x), sys), inputs)
 
-    eqs = [eq for eq in full_equations(sys) if !isdifferenceeq(eq)]
+    eqs = [eq for eq in full_equations(sys)]
     if disturbance_inputs !== nothing
         # Set all disturbance *inputs* to zero (we just want to keep the disturbance state)
         subs = Dict(disturbance_inputs .=> 0)

src/inputoutput.jl

@@ -1,6 +1,6 @@
 using LinearAlgebra
 
-using ModelingToolkit: isdifferenceeq, process_events, get_preprocess_constants
+using ModelingToolkit: process_events, get_preprocess_constants
 
 const MAX_INLINE_NLSOLVE_SIZE = 8
 
@@ -543,8 +543,7 @@ function ODAEProblem{iip}(sys,
     u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults = defs, tofloat = true)
     p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults = defs, tofloat, use_union)
 
-    has_difference = any(isdifferenceeq, eqs)
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
 
     kwargs = filter_kwargs(kwargs)
     if cbs === nothing

src/structural_transformation/codegen.jl

@@ -498,7 +498,7 @@ merge_cb(::Nothing, x) = merge_cb(x, nothing)
 merge_cb(x, ::Nothing) = x
 merge_cb(x, y) = CallbackSet(x, y)
 
-function process_events(sys; callback = nothing, has_difference = false, kwargs...)
+function process_events(sys; callback = nothing, kwargs...)
     if has_continuous_events(sys)
         contin_cb = generate_rootfinding_callback(sys; kwargs...)
     else
@@ -509,9 +509,7 @@ function process_events(sys; callback = nothing, has_difference = false, kwargs.
     else
         discrete_cb = nothing
     end
-    difference_cb = has_difference ? generate_difference_cb(sys; kwargs...) : nothing
 
-    cb = merge_cb(contin_cb, difference_cb)
-    cb = merge_cb(cb, callback)
-    (discrete_cb === nothing) ? cb : CallbackSet(cb, discrete_cb...)
+    cb = merge_cb(contin_cb, callback)
+    (discrete_cb === nothing) ? cb : CallbackSet(contin_cb, discrete_cb...)
 end

src/systems/callbacks.jl

@@ -410,7 +410,7 @@ end
 
 function generate_connection_equations_and_stream_connections(csets::AbstractVector{
         <:ConnectionSet,
-})
+    })
     eqs = Equation[]
     stream_connections = ConnectionSet[]
 

src/systems/connectors.jl

@@ -144,12 +144,11 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys), ps = par
         ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) :
                nothing,
         isdde = false,
-        has_difference = false,
         kwargs...)
     if isdde
         eqs = delay_to_function(sys)
     else
-        eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
+        eqs = [eq for eq in equations(sys)]
     end
     if !implicit_dae
         check_operator_variables(eqs, Differential)
@@ -167,8 +166,7 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys), ps = par
     if isdde
         build_function(rhss, u, DDE_HISTORY_FUN, p, t; kwargs...)
     else
-        pre, sol_states = get_substitutions_and_solved_unknowns(sys,
-            no_postprocess = has_difference)
+        pre, sol_states = get_substitutions_and_solved_unknowns(sys)
 
         if implicit_dae
             build_function(rhss, ddvs, u, p, t; postprocess_fbody = pre,
@@ -226,54 +224,8 @@ function delay_to_function(expr, iv, sts, ps, h)
     end
 end
 
-function generate_difference_cb(sys::ODESystem, dvs = unknowns(sys), ps = parameters(sys);
-        kwargs...)
-    eqs = equations(sys)
-    check_operator_variables(eqs, Difference)
-
-    var2eq = Dict(arguments(eq.lhs)[1] => eq for eq in eqs if isdifference(eq.lhs))
-
-    u = map(x -> time_varying_as_func(value(x), sys), dvs)
-    p = map(x -> time_varying_as_func(value(x), sys), ps)
-    t = get_iv(sys)
-
-    body = map(dvs) do v
-        eq = get(var2eq, v, nothing)
-        eq === nothing && return v
-        d = operation(eq.lhs)
-        d.update ? eq.rhs : eq.rhs + v
-    end
-
-    pre = get_postprocess_fbody(sys)
-    cpre = get_preprocess_constants(body)
-    pre2 = x -> pre(cpre(x))
-    f_oop, f_iip = build_function(body, u, p, t; expression = Val{false},
-        postprocess_fbody = pre2, kwargs...)
-
-    cb_affect! = let f_oop = f_oop, f_iip = f_iip
-        function cb_affect!(integ)
-            if DiffEqBase.isinplace(integ.sol.prob)
-                tmp, = DiffEqBase.get_tmp_cache(integ)
-                f_iip(tmp, integ.u, integ.p, integ.t) # aliasing `integ.u` would be bad.
-                copyto!(integ.u, tmp)
-            else
-                integ.u = f_oop(integ.u, integ.p, integ.t)
-            end
-            return nothing
-        end
-    end
-
-    getdt(eq) = operation(eq.lhs).dt
-    deqs = values(var2eq)
-    dt = getdt(first(deqs))
-    all(dt == getdt(eq) for eq in deqs) ||
-        error("All difference variables should have same time steps.")
-
-    PeriodicCallback(cb_affect!, first(dt))
-end
-
 function calculate_massmatrix(sys::AbstractODESystem; simplify = false)
-    eqs = [eq for eq in equations(sys) if !isdifferenceeq(eq)]
+    eqs = [eq for eq in equations(sys)]
     dvs = unknowns(sys)
     M = zeros(length(eqs), length(eqs))
     unknown2idx = Dict(s => i for (i, s) in enumerate(dvs))
@@ -939,12 +891,10 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
         callback = nothing,
         check_length = true,
         kwargs...) where {iip, specialize}
-    has_difference = any(isdifferenceeq, equations(sys))
     f, u0, p = process_DEProblem(ODEFunction{iip, specialize}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
-        has_difference = has_difference,
         check_length, kwargs...)
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
     inits = []
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
         affects, inits, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...)
@@ -1009,24 +959,16 @@ end
 function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
         parammap = DiffEqBase.NullParameters();
         check_length = true, kwargs...) where {iip}
-    has_difference = any(isdifferenceeq, equations(sys))
     f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap;
-        implicit_dae = true, du0map = du0map,
-        has_difference = has_difference, check_length,
+        implicit_dae = true, du0map = du0map, check_length,
         kwargs...)
     diffvars = collect_differential_variables(sys)
     sts = unknowns(sys)
     differential_vars = map(Base.Fix2(in, diffvars), sts)
     kwargs = filter_kwargs(kwargs)
 
-    if has_difference
-        DAEProblem{iip}(f, du0, u0, tspan, p;
-            difference_cb = generate_difference_cb(sys; kwargs...),
-            differential_vars = differential_vars, kwargs...)
-    else
-        DAEProblem{iip}(f, du0, u0, tspan, p; differential_vars = differential_vars,
-            kwargs...)
-    end
+    DAEProblem{iip}(f, du0, u0, tspan, p; differential_vars = differential_vars,
+        kwargs...)
 end
 
 function generate_history(sys::AbstractODESystem, u0; kwargs...)
@@ -1042,16 +984,14 @@ function DiffEqBase.DDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         callback = nothing,
         check_length = true,
         kwargs...) where {iip}
-    has_difference = any(isdifferenceeq, equations(sys))
     f, u0, p = process_DEProblem(DDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
-        has_difference = has_difference,
         symbolic_u0 = true,
         check_length, kwargs...)
     h_oop, h_iip = generate_history(sys, u0)
     h = h_oop
     u0 = h(p, tspan[1])
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
     inits = []
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
         affects, inits, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...)
@@ -1102,17 +1042,15 @@ function DiffEqBase.SDDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         check_length = true,
         sparsenoise = nothing,
         kwargs...) where {iip}
-    has_difference = any(isdifferenceeq, equations(sys))
     f, u0, p = process_DEProblem(SDDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
-        has_difference = has_difference,
         symbolic_u0 = true,
         check_length, kwargs...)
     h_oop, h_iip = generate_history(sys, u0)
     h(out, p, t) = h_iip(out, p, t)
     h(p, t) = h_oop(p, t)
     u0 = h(p, tspan[1])
-    cbs = process_events(sys; callback, has_difference, kwargs...)
+    cbs = process_events(sys; callback, kwargs...)
     inits = []
     if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing
         affects, inits, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...)