Merge pull request #2436 from AayushSabharwal/as/require-completed

refactor!: require systems to be `complete`d before creating `XProblem`/`XProblemExpr`

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "9.0.0"
+version = "8.76.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

docs/src/basics/Events.md

@@ -73,7 +73,7 @@ function UnitMassWithFriction(k; name)
     ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity
 end
 @named m = UnitMassWithFriction(0.7)
-prob = ODEProblem(m, Pair[], (0, 10pi))
+prob = ODEProblem(complete(m), Pair[], (0, 10pi))
 sol = solve(prob, Tsit5())
 plot(sol)
 ```
@@ -244,7 +244,7 @@ u0 = [N => 0.0]
 tspan = (0.0, 20.0)
 p = [α => 100.0, tinject => 10.0, M => 50]
 @named osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection)
-oprob = ODEProblem(osys, u0, tspan, p)
+oprob = ODEProblem(complete(osys), u0, tspan, p)
 sol = solve(oprob, Tsit5(); tstops = 10.0)
 plot(sol)
 ```

docs/src/examples/parsing.md

@@ -27,6 +27,6 @@ using ModelingToolkit, NonlinearSolve
 vars = union(ModelingToolkit.vars.(eqs)...)
 @named ns = NonlinearSystem(eqs, vars, [])
 
-prob = NonlinearProblem(ns, [1.0, 1.0, 1.0])
+prob = NonlinearProblem(complete(ns), [1.0, 1.0, 1.0])
 sol = solve(prob, NewtonRaphson())
 ```

docs/src/examples/sparse_jacobians.md

@@ -55,6 +55,7 @@ Now let's use `modelingtoolkitize` to generate the symbolic version:
 
 ```@example sparsejac
 sys = modelingtoolkitize(prob);
+sys = complete(sys);
 nothing # hide
 ```
 

docs/src/tutorials/bifurcation_diagram_computation.md

@@ -13,6 +13,7 @@ using ModelingToolkit
 eqs = [0 ~ μ * x - x^3 + α * y,
     0 ~ -y]
 @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+nsys = complete(nsys)
 ```
 
 we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information:
@@ -97,6 +98,7 @@ D = Differential(t)
 eqs = [D(x) ~ μ * x - y - x * (x^2 + y^2),
     D(y) ~ x + μ * y - y * (x^2 + y^2)]
 @named osys = ODESystem(eqs, t)
+osys = complete(osys)
 
 bif_par = μ
 plot_var = x

docs/src/tutorials/modelingtoolkitize.md

@@ -58,5 +58,5 @@ sys = modelingtoolkitize(prob)
 Using this, we can symbolically build the Jacobian and then rebuild the ODEProblem:
 
 ```@example mtkize
-prob_jac = ODEProblem(sys, [], (0.0, 1e5), jac = true)
+prob_jac = ODEProblem(complete(sys), [], (0.0, 1e5), jac = true)
 ```

docs/src/tutorials/nonlinear.md

@@ -17,6 +17,7 @@ eqs = [0 ~ σ * (y - x),
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z]
 @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+ns = complete(ns)
 
 guess = [x => 1.0,
     y => 0.0,

docs/src/tutorials/optimization.md

@@ -50,7 +50,7 @@ u0 = [x => 1.0
 p = [a => 1.0
     b => 100.0]
 
-prob = OptimizationProblem(sys, u0, p, grad = true, hess = true)
+prob = OptimizationProblem(complete(sys), u0, p, grad = true, hess = true)
 solve(prob, GradientDescent())
 ```
 
@@ -71,7 +71,7 @@ cons = [
 @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
 u0 = [x => 0.14
     y => 0.14]
-prob = OptimizationProblem(sys, u0, grad = true, hess = true, cons_j = true, cons_h = true)
+prob = OptimizationProblem(complete(sys), u0, grad = true, hess = true, cons_j = true, cons_h = true)
 solve(prob, IPNewton())
 ```
 

docs/src/tutorials/stochastic_diffeq.md

@@ -24,6 +24,7 @@ noiseeqs = [0.1 * x,
     0.1 * z]
 
 @named de = SDESystem(eqs, noiseeqs, t, [x, y, z], [σ, ρ, β])
+de = complete(de)
 
 u0map = [
     x => 1.0,

ext/MTKBifurcationKitExt.jl

@@ -88,6 +88,9 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
         record_from_solution = BifurcationKit.record_sol_default,
         jac = true,
         kwargs...)
+    if !ModelingToolkit.iscomplete(nsys)
+        error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `BifurcationProblem`")
+    end
     # Creates F and J functions.
     ofun = NonlinearFunction(nsys; jac = jac)
     F = ofun.f
@@ -133,11 +136,14 @@ end
 
 # When input is a ODESystem.
 function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
+    if !ModelingToolkit.iscomplete(osys)
+        error("A completed `ODESystem` is required. Call `complete` or `structural_simplify` on the system before creating a `BifurcationProblem`")
+    end
     nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)],
         unknowns(osys),
         parameters(osys);
         name = nameof(osys))
-    return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
+    return BifurcationKit.BifurcationProblem(complete(nsys), args...; kwargs...)
 end
 
 end # module

src/systems/abstractsystem.jl

@@ -1330,6 +1330,8 @@ information. E.g.
 ```julia-repl
 julia> sys = debug_system(sys);
 
+julia> sys = complete(sys);
+
 julia> prob = ODEProblem(sys, [], (0, 1.0));
 
 julia> du = zero(prob.u0);

src/systems/diffeqs/abstractodesystem.jl

@@ -308,6 +308,9 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = u
         analytic = nothing,
         split_idxs = nothing,
         kwargs...) where {iip, specialize}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression},
         expression_module = eval_module, checkbounds = checkbounds,
         kwargs...)
@@ -504,6 +507,9 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys)
         eval_module = @__MODULE__,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; implicit_dae = true,
         expression = Val{eval_expression},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -579,6 +585,9 @@ function DiffEqBase.DDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys)
         eval_module = @__MODULE__,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `DDEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; isdde = true,
         expression = Val{true},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -603,6 +612,9 @@ function DiffEqBase.SDDEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys
         eval_module = @__MODULE__,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `SDDEFunction`")
+    end
     f_gen = generate_function(sys, dvs, ps; isdde = true,
         expression = Val{true},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -656,6 +668,9 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = unknowns(sys),
         sparsity = false,
         observedfun_exp = nothing,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEFunctionExpr`")
+    end
     f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)
 
     dict = Dict()
@@ -830,6 +845,9 @@ function DAEFunctionExpr{iip}(sys::AbstractODESystem, dvs = unknowns(sys),
         linenumbers = false,
         sparse = false, simplify = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `DAEFunctionExpr`")
+    end
     f_oop, f_iip = generate_function(sys, dvs, ps; expression = Val{true},
         implicit_dae = true, kwargs...)
     fsym = gensym(:f)
@@ -891,6 +909,9 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map =
         callback = nothing,
         check_length = true,
         kwargs...) where {iip, specialize}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
+    end
     f, u0, p = process_DEProblem(ODEFunction{iip, specialize}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
         check_length, kwargs...)
@@ -959,6 +980,9 @@ end
 function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
         parammap = DiffEqBase.NullParameters();
         check_length = true, kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblem`")
+    end
     f, du0, u0, p = process_DEProblem(DAEFunction{iip}, sys, u0map, parammap;
         implicit_dae = true, du0map = du0map, check_length,
         kwargs...)
@@ -984,6 +1008,9 @@ function DiffEqBase.DDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         callback = nothing,
         check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DDEProblem`")
+    end
     f, u0, p = process_DEProblem(DDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
         symbolic_u0 = true,
@@ -1042,6 +1069,9 @@ function DiffEqBase.SDDEProblem{iip}(sys::AbstractODESystem, u0map = [],
         check_length = true,
         sparsenoise = nothing,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `SDDEProblem`")
+    end
     f, u0, p = process_DEProblem(SDDEFunction{iip}, sys, u0map, parammap;
         t = tspan !== nothing ? tspan[1] : tspan,
         symbolic_u0 = true,
@@ -1126,6 +1156,9 @@ struct ODEProblemExpr{iip} end
 function ODEProblemExpr{iip}(sys::AbstractODESystem, u0map, tspan,
         parammap = DiffEqBase.NullParameters(); check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `ODEProblemExpr`")
+    end
     f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap; check_length,
         kwargs...)
     linenumbers = get(kwargs, :linenumbers, true)
@@ -1168,6 +1201,9 @@ struct DAEProblemExpr{iip} end
 function DAEProblemExpr{iip}(sys::AbstractODESystem, du0map, u0map, tspan,
         parammap = DiffEqBase.NullParameters(); check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblemExpr`")
+    end
     f, du0, u0, p = process_DEProblem(DAEFunctionExpr{iip}, sys, u0map, parammap;
         implicit_dae = true, du0map = du0map, check_length,
         kwargs...)
@@ -1216,6 +1252,9 @@ end
 function DiffEqBase.SteadyStateProblem{iip}(sys::AbstractODESystem, u0map,
         parammap = SciMLBase.NullParameters();
         check_length = true, kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `SteadyStateProblem`")
+    end
     f, u0, p = process_DEProblem(ODEFunction{iip}, sys, u0map, parammap;
         steady_state = true,
         check_length, kwargs...)
@@ -1245,6 +1284,9 @@ function SteadyStateProblemExpr{iip}(sys::AbstractODESystem, u0map,
         parammap = SciMLBase.NullParameters();
         check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `SteadyStateProblemExpr`")
+    end
     f, u0, p = process_DEProblem(ODEFunctionExpr{iip}, sys, u0map, parammap;
         steady_state = true,
         check_length, kwargs...)

src/systems/diffeqs/basic_transformations.jl

@@ -31,7 +31,7 @@ u0 = [x => 1.0,
       y => 1.0,
       trJ => 1.0]
 
-prob = ODEProblem(sys2,u0,tspan,p)
+prob = ODEProblem(complete(sys2),u0,tspan,p)
 sol = solve(prob,Tsit5())
 ```
 

src/systems/diffeqs/sdesystem.jl

@@ -322,7 +322,7 @@ parammap = [
     β => 1.0
 ]
 
-probmod = SDEProblem(demod,u0modmap,(0.0,1.0),parammap)
+probmod = SDEProblem(complete(demod),u0modmap,(0.0,1.0),parammap)
 ensemble_probmod = EnsembleProblem(probmod;
           output_func = (sol,i) -> (g(sol[x,end])*sol[demod.weight,end],false),
           )
@@ -393,6 +393,9 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = unknowns(sys),
         jac = false, Wfact = false, eval_expression = true,
         checkbounds = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEFunction`")
+    end
     dvs = scalarize.(dvs)
     ps = scalarize.(ps)
 
@@ -515,6 +518,9 @@ function SDEFunctionExpr{iip}(sys::SDESystem, dvs = unknowns(sys),
         jac = false, Wfact = false,
         sparse = false, linenumbers = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEFunctionExpr`")
+    end
     idx = iip ? 2 : 1
     f = generate_function(sys, dvs, ps; expression = Val{true}, kwargs...)[idx]
     g = generate_diffusion_function(sys, dvs, ps; expression = Val{true}, kwargs...)[idx]
@@ -573,6 +579,9 @@ function DiffEqBase.SDEProblem{iip}(sys::SDESystem, u0map = [], tspan = get_tspa
         parammap = DiffEqBase.NullParameters();
         sparsenoise = nothing, check_length = true,
         callback = nothing, kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEProblem`")
+    end
     f, u0, p = process_DEProblem(SDEFunction{iip}, sys, u0map, parammap; check_length,
         kwargs...)
     cbs = process_events(sys; callback)
@@ -632,6 +641,9 @@ function SDEProblemExpr{iip}(sys::SDESystem, u0map, tspan,
         parammap = DiffEqBase.NullParameters();
         sparsenoise = nothing, check_length = true,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `SDESystem` is required. Call `complete` or `structural_simplify` on the system before creating an `SDEProblemExpr`")
+    end
     f, u0, p = process_DEProblem(SDEFunctionExpr{iip}, sys, u0map, parammap; check_length,
         kwargs...)
     linenumbers = get(kwargs, :linenumbers, true)

src/systems/jumps/jumpsystem.jl

@@ -119,6 +119,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
             connector_type, devents, metadata, gui_metadata, complete)
     end
 end
+JumpSystem(tag, ap, iv, states, ps, var_to_name, args...; kwargs...) = JumpSystem{typeof(ap)}(tag, ap, iv, states, ps, var_to_name, args...; kwargs...)
 
 function JumpSystem(eqs, iv, unknowns, ps;
         observed = Equation[],
@@ -289,14 +290,17 @@ using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => 0.1 / 1000, γ => 0.01]
 tspan = (0.0, 250.0)
-dprob = DiscreteProblem(js, u₀map, tspan, parammap)
+dprob = DiscreteProblem(complete(js), u₀map, tspan, parammap)
 ```
 """
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothing},
         parammap = DiffEqBase.NullParameters();
         checkbounds = false,
         use_union = true,
         kwargs...)
+    if !iscomplete(sys)
+        error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`")
+    end
     dvs = unknowns(sys)
     ps = parameters(sys)
 
@@ -344,7 +348,7 @@ using DiffEqBase, JumpProcesses
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => 0.1 / 1000, γ => 0.01]
 tspan = (0.0, 250.0)
-dprob = DiscreteProblem(js, u₀map, tspan, parammap)
+dprob = DiscreteProblem(complete(js), u₀map, tspan, parammap)
 ```
 """
 struct DiscreteProblemExpr{iip} end
@@ -353,6 +357,9 @@ function DiscreteProblemExpr{iip}(sys::JumpSystem, u0map, tspan::Union{Tuple, No
         parammap = DiffEqBase.NullParameters();
         use_union = false,
         kwargs...) where {iip}
+    if !iscomplete(sys)
+        error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblemExpr`")
+    end
     dvs = unknowns(sys)
     ps = parameters(sys)
     defs = defaults(sys)
@@ -382,12 +389,15 @@ Generates a JumpProblem from a JumpSystem.
 Continuing the example from the [`DiscreteProblem`](@ref) definition:
 
 ```julia
-jprob = JumpProblem(js, dprob, Direct())
+jprob = JumpProblem(complete(js), dprob, Direct())
 sol = solve(jprob, SSAStepper())
 ```
 """
 function JumpProcesses.JumpProblem(js::JumpSystem, prob, aggregator; callback = nothing,
         kwargs...)
+    if !iscomplete(js)
+        error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `JumpProblem`")
+    end
     unknowntoid = Dict(value(unknown) => i for (i, unknown) in enumerate(unknowns(js)))
     eqs = equations(js)
     invttype = prob.tspan[1] === nothing ? Float64 : typeof(1 / prob.tspan[2])