refactor: initial conditions represent values at t-1

Author: AayushSabharwal

src/systems/diffeqs/abstractodesystem.jl

@@ -1045,7 +1045,7 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
         discrete_cbs = map(affects, clocks, svs) do affect, clock, sv
             if clock isa Clock
                 PeriodicCallback(DiscreteSaveAffect(affect, sv), clock.dt;
-                    final_affect = true)
+                    final_affect = true, initial_affect = true)
             elseif clock isa SolverStepClock
                 affect = DiscreteSaveAffect(affect, sv)
                 DiscreteCallback(Returns(true), affect,
@@ -1080,7 +1080,7 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
     prob = ODEProblem{iip}(f, u0, tspan, p, pt; kwargs1..., kwargs...)
     if !isempty(inits)
         for init in inits
-            init(prob.u0, prob.p, tspan[1])
+            # init(prob.u0, prob.p, tspan[1])
         end
     end
     prob

src/systems/discrete_system/discrete_system.jl

@@ -253,6 +253,7 @@ function SciMLBase.DiscreteProblem(
 
     f, u0, p = process_DiscreteProblem(
         DiscreteFunction, sys, u0map, parammap; eval_expression, eval_module)
+    u0 = f(u0, p, tspan[1])
     DiscreteProblem(f, u0, tspan, p; kwargs...)
 end
 

test/clock.jl

@@ -112,14 +112,18 @@ prob = ODEProblem(ss, [x => 0.1], (0.0, Tf),
     [kp => 1.0; ud => 2.1; ud(k - 1) => 2.0])
 # create integrator so callback is evaluated at t=0 and we can test correct param values
 int = init(prob, Tsit5(); kwargshandle = KeywordArgSilent)
-@test sort(vcat(int.p...)) == [0.1, 1.0, 2.0, 2.1, 2.1] # yd, kp, ud(k-1), ud, Hold(ud)
+@test sort(vcat(int.p...)) == [0.1, 1.0, 2.1, 2.1, 2.1] # yd, kp, ud(k-1), ud, Hold(ud)
+prob = ODEProblem(ss, [x => 0.1], (0.0, Tf),
+    [kp => 1.0; ud => 2.1; ud(k - 1) => 2.0]) # recreate problem to empty saved values
 sol = solve(prob, Tsit5(), kwargshandle = KeywordArgSilent)
 
 ss_nosplit = structural_simplify(sys; split = false)
 prob_nosplit = ODEProblem(ss_nosplit, [x => 0.1], (0.0, Tf),
     [kp => 1.0; ud => 2.1; ud(k - 1) => 2.0])
 int = init(prob_nosplit, Tsit5(); kwargshandle = KeywordArgSilent)
-@test sort(int.p) == [0.1, 1.0, 2.0, 2.1, 2.1] # yd, kp, ud(k-1), ud, Hold(ud)
+@test sort(int.p) == [0.1, 1.0, 2.1, 2.1, 2.1] # yd, kp, ud(k-1), ud, Hold(ud)
+prob_nosplit = ODEProblem(ss_nosplit, [x => 0.1], (0.0, Tf),
+    [kp => 1.0; ud => 2.1; ud(k - 1) => 2.0]) # recreate problem to empty saved values
 sol_nosplit = solve(prob_nosplit, Tsit5(), kwargshandle = KeywordArgSilent)
 # For all inputs in parameters, just initialize them to 0.0, and then set them
 # in the callback.
@@ -145,7 +149,8 @@ function affect!(integrator, saved_values)
     nothing
 end
 saved_values = SavedValues(Float64, Vector{Float64})
-cb = PeriodicCallback(Base.Fix2(affect!, saved_values), 0.1; final_affect = true)
+cb = PeriodicCallback(
+    Base.Fix2(affect!, saved_values), 0.1; final_affect = true, initial_affect = true)
 #                                           kp   ud
 prob = ODEProblem(foo!, [0.1], (0.0, Tf), [1.0, 2.1, 2.0], callback = cb)
 sol2 = solve(prob, Tsit5())
@@ -308,8 +313,8 @@ if VERSION >= v"1.7"
         integrator.p[3] = ud2
         nothing
     end
-    cb1 = PeriodicCallback(affect1!, dt; final_affect = true)
-    cb2 = PeriodicCallback(affect2!, dt2; final_affect = true)
+    cb1 = PeriodicCallback(affect1!, dt; final_affect = true, initial_affect = true)
+    cb2 = PeriodicCallback(affect2!, dt2; final_affect = true, initial_affect = true)
     cb = CallbackSet(cb1, cb2)
     #                                           kp   ud1  ud2
     prob = ODEProblem(foo!, [0.0], (0.0, 1.0), [1.0, 1.0, 1.0], callback = cb)
@@ -438,10 +443,10 @@ y = res.y[:]
 prob = ODEProblem(ssys,
     [model.plant.x => 0.0; model.controller.kp => 2.0; model.controller.ki => 2.0],
     (0.0, Tf))
-
-@test prob.ps[Hold(ssys.holder.input.u)] == 1 # constant output * kp issue https://github.com/SciML/ModelingToolkit.jl/issues/2356
-@test prob.ps[ssys.controller.x] == 0 # c2d
-@test prob.ps[Sample(d)(ssys.sampler.input.u)] == 0 # disc state
+int = init(prob, Tsit5(); kwargshandle = KeywordArgSilent)
+@test int.ps[Hold(ssys.holder.input.u)] == 2 # constant output * kp issue https://github.com/SciML/ModelingToolkit.jl/issues/2356
+@test int.ps[ssys.controller.x] == 1 # c2d
+@test int.ps[Sample(d)(ssys.sampler.input.u)] == 0 # disc state
 sol = solve(prob,
     Tsit5(),
     kwargshandle = KeywordArgSilent,

test/discrete_system.jl

@@ -99,7 +99,7 @@ function sir_map!(u_diff, u, p, t)
     end
     nothing
 end;
-u0 = [990.0, 10.0, 0.0];
+u0 = prob_map2.u0;
 p = [0.05, 10.0, 0.25, 0.1];
 prob_map = DiscreteProblem(sir_map!, u0, tspan, p);
 sol_map2 = solve(prob_map, FunctionMap());
@@ -216,4 +216,4 @@ eqs = [u ~ 1
 prob = DiscreteProblem(de, [x => 0.0], (0, 10))
 sol = solve(prob, FunctionMap())
 
-@test reduce(vcat, sol.u) == 0:10
+@test reduce(vcat, sol.u) == 1:11