diff --git a/ext/MTKCasADiDynamicOptExt.jl b/ext/MTKCasADiDynamicOptExt.jl index a5400e3c8d..8f6b2c345e 100644 --- a/ext/MTKCasADiDynamicOptExt.jl +++ b/ext/MTKCasADiDynamicOptExt.jl @@ -76,8 +76,10 @@ function MTK.CasADiDynamicOptProblem(sys::System, u0map, tspan, pmap; steps = nothing, guesses = Dict(), kwargs...) MTK.warn_overdetermined(sys, u0map) - _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) - f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, _u0map, pmap; + _u0map = has_alg_eqs(sys) ? MTK.to_varmap(u0map, unknowns(sys)) : + merge(Dict(u0map), Dict(guesses)) + pmap = MTK.to_varmap(pmap, parameters(sys)) + f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, merge(_u0map, pmap); t = tspan !== nothing ? tspan[1] : tspan, output_type = MX, kwargs...) pmap = MTK.recursive_unwrap(MTK.AnyDict(pmap)) diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl index f4fff61dff..40eb6f5264 100644 --- a/ext/MTKInfiniteOptExt.jl +++ b/ext/MTKInfiniteOptExt.jl @@ -59,8 +59,10 @@ function MTK.JuMPDynamicOptProblem(sys::System, u0map, tspan, pmap; steps = nothing, guesses = Dict(), kwargs...) MTK.warn_overdetermined(sys, u0map) - _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) - f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, _u0map, pmap; + _u0map = has_alg_eqs(sys) ? MTK.to_varmap(u0map, unknowns(sys)) : + merge(Dict(u0map), Dict(guesses)) + pmap = MTK.to_varmap(pmap, parameters(sys)) + f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, merge(_u0map, pmap); t = tspan !== nothing ? tspan[1] : tspan, kwargs...) pmap = MTK.recursive_unwrap(MTK.AnyDict(pmap)) @@ -86,8 +88,10 @@ function MTK.InfiniteOptDynamicOptProblem(sys::System, u0map, tspan, pmap; steps = nothing, guesses = Dict(), kwargs...) MTK.warn_overdetermined(sys, u0map) - _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) - f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, _u0map, pmap; + _u0map = has_alg_eqs(sys) ? MTK.to_varmap(u0map, unknowns(sys)) : + merge(Dict(u0map), Dict(guesses)) + pmap = MTK.to_varmap(pmap, parameters(sys)) + f, u0, p = MTK.process_SciMLProblem(ODEInputFunction, sys, merge(_u0map, pmap); t = tspan !== nothing ? tspan[1] : tspan, kwargs...) pmap = MTK.recursive_unwrap(MTK.AnyDict(pmap)) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index e94c8bb5af..63468917d1 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -234,7 +234,7 @@ PrecompileTools.@compile_workload begin using ModelingToolkit @variables x(ModelingToolkit.t_nounits) @named sys = System([ModelingToolkit.D_nounits(x) ~ -x], ModelingToolkit.t_nounits) - prob = ODEProblem(mtkcompile(sys), [x => 30.0], (0, 100), [], jac = true) + prob = ODEProblem(mtkcompile(sys), [x => 30.0], (0, 100), jac = true) @mtkmodel __testmod__ begin @constants begin c = 1.0 diff --git a/src/deprecations.jl b/src/deprecations.jl index c8aee32253..73df50ecc4 100644 --- a/src/deprecations.jl +++ b/src/deprecations.jl @@ -8,4 +8,3 @@ macro mtkbuild(exprs...) @mtkcompile $(exprs...) end |> esc end - diff --git a/src/linearization.jl b/src/linearization.jl index 757d6eece5..0c174e4cc3 100644 --- a/src/linearization.jl +++ b/src/linearization.jl @@ -73,7 +73,7 @@ function linearization_function(sys::AbstractSystem, inputs, end prob = ODEProblem{true, SciMLBase.FullSpecialize}( - sys, op, (nothing, nothing), p; allow_incomplete = true, + sys, merge(op, anydict(p)), (nothing, nothing); allow_incomplete = true, algebraic_only = true, guesses) u0 = state_values(prob) diff --git a/src/problems/bvproblem.jl b/src/problems/bvproblem.jl index 82783077c1..30bd7e61c3 100644 --- a/src/problems/bvproblem.jl +++ b/src/problems/bvproblem.jl @@ -44,7 +44,7 @@ If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting `BVProblem` must be solved using BVDAE solvers, such as Ascher. """ @fallback_iip_specialize function SciMLBase.BVProblem{iip, spec}( - sys::System, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; check_compatibility = true, cse = true, checkbounds = false, eval_expression = false, eval_module = @__MODULE__, expression = Val{false}, guesses = Dict(), callback = nothing, @@ -55,22 +55,23 @@ If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting # Systems without algebraic equations should use both fixed values + guesses # for initialization. - _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) + _op = has_alg_eqs(sys) ? op : merge(Dict(op), Dict(guesses)) fode, u0, p = process_SciMLProblem( - ODEFunction{iip, spec}, sys, _u0map, parammap; guesses, + ODEFunction{iip, spec}, sys, _op; guesses, t = tspan !== nothing ? tspan[1] : tspan, check_compatibility = false, cse, checkbounds, time_dependent_init = false, expression, kwargs...) dvs = unknowns(sys) stidxmap = Dict([v => i for (i, v) in enumerate(dvs)]) - u0_idxs = has_alg_eqs(sys) ? collect(1:length(dvs)) : [stidxmap[k] for (k, v) in u0map] + u0_idxs = has_alg_eqs(sys) ? collect(1:length(dvs)) : + [stidxmap[k] for (k, v) in op if haskey(stidxmap, k)] fbc = generate_boundary_conditions( sys, u0, u0_idxs, tspan[1]; expression = Val{false}, wrap_gfw = Val{true}, cse, checkbounds) - if (length(constraints(sys)) + length(u0map) > length(dvs)) - @warn "The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The BVP solvers will default to doing a nonlinear least-squares optimization." + if (length(constraints(sys)) + length(op) > length(dvs)) + @warn "The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by op) exceeds the total number of states. The BVP solvers will default to doing a nonlinear least-squares optimization." end kwargs = process_kwargs(sys; expression, kwargs...) diff --git a/src/problems/daeproblem.jl b/src/problems/daeproblem.jl index d1f8893cd5..860f2d2da2 100644 --- a/src/problems/daeproblem.jl +++ b/src/problems/daeproblem.jl @@ -62,15 +62,15 @@ end @fallback_iip_specialize function SciMLBase.DAEProblem{iip, spec}( - sys::System, du0map, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; callback = nothing, check_length = true, eval_expression = false, eval_module = @__MODULE__, check_compatibility = true, expression = Val{false}, kwargs...) where {iip, spec} check_complete(sys, DAEProblem) check_compatibility && check_compatible_system(DAEProblem, sys) - f, du0, u0, p = process_SciMLProblem(DAEFunction{iip, spec}, sys, u0map, parammap; - du0map, t = tspan !== nothing ? tspan[1] : tspan, check_length, eval_expression, + f, du0, u0, p = process_SciMLProblem(DAEFunction{iip, spec}, sys, op; + t = tspan !== nothing ? tspan[1] : tspan, check_length, eval_expression, eval_module, check_compatibility, implicit_dae = true, expression, kwargs...) kwargs = process_kwargs(sys; expression, callback, eval_expression, eval_module, diff --git a/src/problems/ddeproblem.jl b/src/problems/ddeproblem.jl index dfc660effd..6f072b3d51 100644 --- a/src/problems/ddeproblem.jl +++ b/src/problems/ddeproblem.jl @@ -42,14 +42,14 @@ end @fallback_iip_specialize function SciMLBase.DDEProblem{iip, spec}( - sys::System, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; callback = nothing, check_length = true, cse = true, checkbounds = false, eval_expression = false, eval_module = @__MODULE__, check_compatibility = true, u0_constructor = identity, expression = Val{false}, kwargs...) where {iip, spec} check_complete(sys, DDEProblem) check_compatibility && check_compatible_system(DDEProblem, sys) - f, u0, p = process_SciMLProblem(DDEFunction{iip, spec}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(DDEFunction{iip, spec}, sys, op; t = tspan !== nothing ? tspan[1] : tspan, check_length, cse, checkbounds, eval_expression, eval_module, check_compatibility, symbolic_u0 = true, expression, u0_constructor, kwargs...) diff --git a/src/problems/discreteproblem.jl b/src/problems/discreteproblem.jl index d77efed225..f640ee9ff5 100644 --- a/src/problems/discreteproblem.jl +++ b/src/problems/discreteproblem.jl @@ -39,7 +39,7 @@ end @fallback_iip_specialize function SciMLBase.DiscreteProblem{iip, spec}( - sys::System, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; check_compatibility = true, expression = Val{false}, kwargs...) where {iip, spec} check_complete(sys, DiscreteProblem) check_compatibility && check_compatible_system(DiscreteProblem, sys) @@ -47,7 +47,7 @@ end dvs = unknowns(sys) u0map = to_varmap(u0map, dvs) add_toterms!(u0map; replace = true) - f, u0, p = process_SciMLProblem(DiscreteFunction{iip, spec}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(DiscreteFunction{iip, spec}, sys, op; t = tspan !== nothing ? tspan[1] : tspan, check_compatibility, expression, kwargs...) diff --git a/src/problems/implicitdiscreteproblem.jl b/src/problems/implicitdiscreteproblem.jl index 476265e049..5b61e34df5 100644 --- a/src/problems/implicitdiscreteproblem.jl +++ b/src/problems/implicitdiscreteproblem.jl @@ -43,16 +43,16 @@ end @fallback_iip_specialize function SciMLBase.ImplicitDiscreteProblem{iip, spec}( - sys::System, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; check_compatibility = true, expression = Val{false}, kwargs...) where {iip, spec} check_complete(sys, ImplicitDiscreteProblem) check_compatibility && check_compatible_system(ImplicitDiscreteProblem, sys) dvs = unknowns(sys) - u0map = to_varmap(u0map, dvs) - add_toterms!(u0map; replace = true) + op = to_varmap(op, dvs) + add_toterms!(op; replace = true) f, u0, p = process_SciMLProblem( - ImplicitDiscreteFunction{iip, spec}, sys, u0map, parammap; + ImplicitDiscreteFunction{iip, spec}, sys, op; t = tspan !== nothing ? tspan[1] : tspan, check_compatibility, expression, kwargs...) diff --git a/src/problems/initializationproblem.jl b/src/problems/initializationproblem.jl index be535da3e6..e957fa7652 100644 --- a/src/problems/initializationproblem.jl +++ b/src/problems/initializationproblem.jl @@ -2,8 +2,7 @@ struct InitializationProblem{iip, specialization} end """ ```julia -InitializationProblem{iip}(sys::AbstractSystem, t, u0map, - parammap = DiffEqBase.NullParameters(); +InitializationProblem{iip}(sys::AbstractSystem, t, op; version = nothing, tgrad = false, jac = false, checkbounds = false, sparse = false, @@ -20,8 +19,7 @@ initial conditions for the given DAE. """ @fallback_iip_specialize function InitializationProblem{iip, specialize}( sys::AbstractSystem, - t, u0map = [], - parammap = DiffEqBase.NullParameters(); + t, op = Dict(); guesses = [], check_length = true, warn_initialize_determined = true, @@ -37,18 +35,24 @@ initial conditions for the given DAE. if !iscomplete(sys) error("A completed system is required. Call `complete` or `mtkcompile` on the system before creating an `ODEProblem`") end - if isempty(u0map) && get_initializesystem(sys) !== nothing + has_u0_ics = false + op = copy(anydict(op)) + for k in keys(op) + has_u0_ics |= is_variable(sys, k) || isdifferential(k) || + symbolic_type(k) == ArraySymbolic() && + is_sized_array_symbolic(k) && is_variable(sys, first(collect(k))) + end + if !has_u0_ics && get_initializesystem(sys) !== nothing isys = get_initializesystem(sys; initialization_eqs, check_units) simplify_system = false - elseif isempty(u0map) && get_initializesystem(sys) === nothing + elseif !has_u0_ics && get_initializesystem(sys) === nothing isys = generate_initializesystem( - sys; initialization_eqs, check_units, pmap = parammap, - guesses, algebraic_only) + sys; initialization_eqs, check_units, op, guesses, algebraic_only) simplify_system = true else isys = generate_initializesystem( - sys; u0map, initialization_eqs, check_units, time_dependent_init, - pmap = parammap, guesses, algebraic_only) + sys; op, initialization_eqs, check_units, time_dependent_init, + guesses, algebraic_only) simplify_system = true end @@ -106,20 +110,17 @@ initial conditions for the given DAE. @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. $(scc_message)To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true" end - parammap = recursive_unwrap(anydict(parammap)) if t !== nothing - parammap[get_iv(sys)] = t + op[get_iv(sys)] = t end - filter!(kvp -> kvp[2] !== missing, parammap) + filter!(kvp -> kvp[2] !== missing, op) - u0map = to_varmap(u0map, unknowns(sys)) if isempty(guesses) guesses = Dict() end - filter_missing_values!(u0map) - filter_missing_values!(parammap) - u0map = merge(ModelingToolkit.guesses(sys), todict(guesses), u0map) + filter_missing_values!(op) + op = merge(ModelingToolkit.guesses(sys), todict(guesses), op) TProb = if neqs == nunknown && isempty(unassigned_vars) if use_scc && neqs > 0 @@ -135,8 +136,7 @@ initial conditions for the given DAE. else NonlinearLeastSquaresProblem end - TProb{iip}(isys, u0map, parammap; kwargs..., - build_initializeprob = false, is_initializeprob = true) + TProb{iip}(isys, op; kwargs..., build_initializeprob = false, is_initializeprob = true) end const INCOMPLETE_INITIALIZATION_MESSAGE = """ diff --git a/src/problems/intervalnonlinearproblem.jl b/src/problems/intervalnonlinearproblem.jl index 9c9e5d8688..65e7f7cb62 100644 --- a/src/problems/intervalnonlinearproblem.jl +++ b/src/problems/intervalnonlinearproblem.jl @@ -33,7 +33,9 @@ function SciMLBase.IntervalNonlinearProblem( check_compatibility && check_compatible_system(IntervalNonlinearProblem, sys) u0map = unknowns(sys) .=> uspan[1] - f, u0, p = process_SciMLProblem(IntervalNonlinearFunction, sys, u0map, parammap; + op = anydict([unknowns(sys)[1] => uspan[1]]) + merge!(op, to_varmap(parammap, parameters(sys))) + f, u0, p = process_SciMLProblem(IntervalNonlinearFunction, sys, op; check_compatibility, expression, kwargs...) kwargs = process_kwargs(sys; kwargs...) diff --git a/src/problems/jumpproblem.jl b/src/problems/jumpproblem.jl index 242b9549fa..28de46c35f 100644 --- a/src/problems/jumpproblem.jl +++ b/src/problems/jumpproblem.jl @@ -1,5 +1,5 @@ @fallback_iip_specialize function JumpProcesses.JumpProblem{iip, spec}( - sys::System, u0map, tspan::Union{Tuple, Nothing}, pmap = SciMLBase.NullParameters(); + sys::System, op, tspan::Union{Tuple, Nothing}; check_compatibility = true, eval_expression = false, eval_module = @__MODULE__, checkbounds = false, cse = true, aggregator = JumpProcesses.NullAggregator(), callback = nothing, rng = nothing, kwargs...) where {iip, spec} @@ -13,18 +13,18 @@ if (has_vrjs || has_eqs) if has_eqs && has_noise prob = SDEProblem{iip, spec}( - sys, u0map, tspan, pmap; check_compatibility = false, + sys, op, tspan; check_compatibility = false, build_initializeprob = false, checkbounds, cse, check_length = false, kwargs...) elseif has_eqs prob = ODEProblem{iip, spec}( - sys, u0map, tspan, pmap; check_compatibility = false, + sys, op, tspan; check_compatibility = false, build_initializeprob = false, checkbounds, cse, check_length = false, kwargs...) else - _, u0, p = process_SciMLProblem(EmptySciMLFunction{iip}, sys, u0map, pmap; - t = tspan === nothing ? nothing : tspan[1], tofloat = false, - check_length = false, build_initializeprob = false) + _, u0, p = process_SciMLProblem(EmptySciMLFunction{iip}, sys, op; + t = tspan === nothing ? nothing : tspan[1], + check_length = false, build_initializeprob = false, kwargs...) observedfun = ObservedFunctionCache(sys; eval_expression, eval_module, checkbounds, cse) f = (du, u, p, t) -> (du .= 0; nothing) @@ -32,8 +32,8 @@ prob = ODEProblem{true}(df, u0, tspan, p; kwargs...) end else - _f, u0, p = process_SciMLProblem(EmptySciMLFunction{iip}, sys, u0map, pmap; - t = tspan === nothing ? nothing : tspan[1], tofloat = false, check_length = false, build_initializeprob = false, cse) + _f, u0, p = process_SciMLProblem(EmptySciMLFunction{iip}, sys, op; + t = tspan === nothing ? nothing : tspan[1], check_length = false, build_initializeprob = false, cse, kwargs...) f = DiffEqBase.DISCRETE_INPLACE_DEFAULT observedfun = ObservedFunctionCache( diff --git a/src/problems/nonlinearproblem.jl b/src/problems/nonlinearproblem.jl index 0d3ffad17f..ef36179e06 100644 --- a/src/problems/nonlinearproblem.jl +++ b/src/problems/nonlinearproblem.jl @@ -57,7 +57,7 @@ end @fallback_iip_specialize function SciMLBase.NonlinearProblem{iip, spec}( - sys::System, u0map, parammap = SciMLBase.NullParameters(); expression = Val{false}, + sys::System, op; expression = Val{false}, check_length = true, check_compatibility = true, kwargs...) where {iip, spec} check_complete(sys, NonlinearProblem) if is_time_dependent(sys) @@ -65,7 +65,7 @@ end end check_compatibility && check_compatible_system(NonlinearProblem, sys) - f, u0, p = process_SciMLProblem(NonlinearFunction{iip, spec}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(NonlinearFunction{iip, spec}, sys, op; check_length, check_compatibility, expression, kwargs...) kwargs = process_kwargs(sys; kwargs...) @@ -75,12 +75,12 @@ end end @fallback_iip_specialize function SciMLBase.NonlinearLeastSquaresProblem{iip, spec}( - sys::System, u0map, parammap = DiffEqBase.NullParameters(); check_length = false, + sys::System, op; check_length = false, check_compatibility = true, expression = Val{false}, kwargs...) where {iip, spec} check_complete(sys, NonlinearLeastSquaresProblem) check_compatibility && check_compatible_system(NonlinearLeastSquaresProblem, sys) - f, u0, p = process_SciMLProblem(NonlinearFunction{iip}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(NonlinearFunction{iip}, sys, op; check_length, expression, kwargs...) kwargs = process_kwargs(sys; kwargs...) diff --git a/src/problems/odeproblem.jl b/src/problems/odeproblem.jl index 99cff129e0..e2901292ee 100644 --- a/src/problems/odeproblem.jl +++ b/src/problems/odeproblem.jl @@ -65,14 +65,14 @@ end @fallback_iip_specialize function SciMLBase.ODEProblem{iip, spec}( - sys::System, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; callback = nothing, check_length = true, eval_expression = false, expression = Val{false}, eval_module = @__MODULE__, check_compatibility = true, kwargs...) where {iip, spec} check_complete(sys, ODEProblem) check_compatibility && check_compatible_system(ODEProblem, sys) - f, u0, p = process_SciMLProblem(ODEFunction{iip, spec}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(ODEFunction{iip, spec}, sys, op; t = tspan !== nothing ? tspan[1] : tspan, check_length, eval_expression, eval_module, expression, check_compatibility, kwargs...) @@ -98,12 +98,12 @@ Generates an SteadyStateProblem from a `System` of ODEs and allows for automatic symbolically calculating numerical enhancements. """ @fallback_iip_specialize function DiffEqBase.SteadyStateProblem{iip, spec}( - sys::System, u0map, parammap; check_length = true, check_compatibility = true, + sys::System, op; check_length = true, check_compatibility = true, expression = Val{false}, kwargs...) where {iip, spec} check_complete(sys, SteadyStateProblem) check_compatibility && check_compatible_system(SteadyStateProblem, sys) - f, u0, p = process_SciMLProblem(ODEFunction{iip}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(ODEFunction{iip}, sys, op; steady_state = true, check_length, check_compatibility, expression, force_initialization_time_independent = true, kwargs...) diff --git a/src/problems/optimizationproblem.jl b/src/problems/optimizationproblem.jl index e0de2f78ff..3f121e3d2a 100644 --- a/src/problems/optimizationproblem.jl +++ b/src/problems/optimizationproblem.jl @@ -87,13 +87,13 @@ function SciMLBase.OptimizationProblem(sys::System, args...; kwargs...) end function SciMLBase.OptimizationProblem{iip}( - sys::System, u0map, parammap = SciMLBase.NullParameters(); lb = nothing, + sys::System, op; lb = nothing, ub = nothing, check_compatibility = true, expression = Val{false}, kwargs...) where {iip} check_complete(sys, OptimizationProblem) check_compatibility && check_compatible_system(OptimizationProblem, sys) - f, u0, p = process_SciMLProblem(OptimizationFunction{iip}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(OptimizationFunction{iip}, sys, op; check_compatibility, tofloat = false, check_length = false, expression, kwargs...) dvs = unknowns(sys) @@ -114,7 +114,7 @@ function SciMLBase.OptimizationProblem{iip}( end ps = parameters(sys) - defs = merge(defaults(sys), to_varmap(parammap, ps), to_varmap(u0map, dvs)) + defs = merge(defaults(sys), to_varmap(op, dvs)) lb = varmap_to_vars(dvs .=> lb, dvs; defaults = defs, tofloat = false) ub = varmap_to_vars(dvs .=> ub, dvs; defaults = defs, tofloat = false) diff --git a/src/problems/sccnonlinearproblem.jl b/src/problems/sccnonlinearproblem.jl index 63058a08f3..3e00170ad9 100644 --- a/src/problems/sccnonlinearproblem.jl +++ b/src/problems/sccnonlinearproblem.jl @@ -69,9 +69,8 @@ function SciMLBase.SCCNonlinearProblem(sys::System, args...; kwargs...) SCCNonlinearProblem{true}(sys, args...; kwargs...) end -function SciMLBase.SCCNonlinearProblem{iip}(sys::System, u0map, - parammap = SciMLBase.NullParameters(); eval_expression = false, eval_module = @__MODULE__, - cse = true, kwargs...) where {iip} +function SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; eval_expression = false, + eval_module = @__MODULE__, cse = true, kwargs...) where {iip} if !iscomplete(sys) || get_tearing_state(sys) === nothing error("A simplified `System` is required. Call `mtkcompile` on the system before creating an `SCCNonlinearProblem`.") end @@ -85,7 +84,7 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, u0map, if length(var_sccs) == 1 return NonlinearProblem{iip}( - sys, u0map, parammap; eval_expression, eval_module, kwargs...) + sys, op; eval_expression, eval_module, kwargs...) end condensed_graph = MatchedCondensationGraph( @@ -102,7 +101,7 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::System, u0map, obs = observed(sys) _, u0, p = process_SciMLProblem( - EmptySciMLFunction{iip}, sys, u0map, parammap; eval_expression, eval_module, kwargs...) + EmptySciMLFunction{iip}, sys, op; eval_expression, eval_module, kwargs...) explicitfuns = [] nlfuns = [] diff --git a/src/problems/sddeproblem.jl b/src/problems/sddeproblem.jl index 3bc20c0412..0daa04fcf7 100644 --- a/src/problems/sddeproblem.jl +++ b/src/problems/sddeproblem.jl @@ -43,7 +43,7 @@ end @fallback_iip_specialize function SciMLBase.SDDEProblem{iip, spec}( - sys::System, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; callback = nothing, check_length = true, cse = true, checkbounds = false, eval_expression = false, eval_module = @__MODULE__, check_compatibility = true, u0_constructor = identity, sparse = false, sparsenoise = sparse, @@ -51,7 +51,7 @@ end check_complete(sys, SDDEProblem) check_compatibility && check_compatible_system(SDDEProblem, sys) - f, u0, p = process_SciMLProblem(SDDEFunction{iip, spec}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(SDDEFunction{iip, spec}, sys, op; t = tspan !== nothing ? tspan[1] : tspan, check_length, cse, checkbounds, eval_expression, eval_module, check_compatibility, sparse, symbolic_u0 = true, expression, u0_constructor, kwargs...) diff --git a/src/problems/sdeproblem.jl b/src/problems/sdeproblem.jl index a6270973f2..57ca120641 100644 --- a/src/problems/sdeproblem.jl +++ b/src/problems/sdeproblem.jl @@ -67,14 +67,14 @@ end @fallback_iip_specialize function SciMLBase.SDEProblem{iip, spec}( - sys::System, u0map, tspan, parammap = SciMLBase.NullParameters(); + sys::System, op, tspan; callback = nothing, check_length = true, eval_expression = false, eval_module = @__MODULE__, check_compatibility = true, sparse = false, sparsenoise = sparse, expression = Val{false}, kwargs...) where {iip, spec} check_complete(sys, SDEProblem) check_compatibility && check_compatible_system(SDEProblem, sys) - f, u0, p = process_SciMLProblem(SDEFunction{iip, spec}, sys, u0map, parammap; + f, u0, p = process_SciMLProblem(SDEFunction{iip, spec}, sys, op; t = tspan !== nothing ? tspan[1] : tspan, check_length, eval_expression, eval_module, check_compatibility, sparse, expression, kwargs...) diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index f70669eca3..db6e6b9648 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -870,8 +870,9 @@ function compile_equational_affect( u_getter = getsym(affsys, [sys_map[u] for u in dvs_to_update]) p_getter = getsym(affsys, [sys_map[p] for p in ps_to_update]) - affprob = ImplicitDiscreteProblem(affsys, [dv => 0 for dv in unknowns(affsys)], - (0, 0), [p => 0.0 for p in parameters(affsys)]; + affprob = ImplicitDiscreteProblem( + affsys, Pair[unknowns(affsys) .=> 0; parameters(affsys) .=> 0], + (0, 0); build_initializeprob = false, check_length = false, eval_expression, eval_module, check_compatibility = false) diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl index ddc0b72aed..d74e4bd226 100644 --- a/src/systems/nonlinear/homotopy_continuation.jl +++ b/src/systems/nonlinear/homotopy_continuation.jl @@ -549,13 +549,13 @@ function HomotopyContinuationProblem{false}(sys::System, args...; kwargs...) end function HomotopyContinuationProblem{iip, spec}( - sys::System, u0map, pmap = SciMLBase.NullParameters(); + sys::System, op; kwargs...) where {iip, spec} if !iscomplete(sys) error("A completed `System` is required. Call `complete` or `mtkcompile` on the system before creating a `HomotopyContinuationProblem`") end f, u0, p = process_SciMLProblem( - HomotopyNonlinearFunction{iip, spec}, sys, u0map, pmap; kwargs...) + HomotopyNonlinearFunction{iip, spec}, sys, op; kwargs...) kwargs = filter_kwargs(kwargs) return NonlinearProblem{iip}(f, u0, p; kwargs...) diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index a364caaa18..ab09c71280 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -13,8 +13,7 @@ $(TYPEDSIGNATURES) Generate `System` of nonlinear equations which initializes a problem from specified initial conditions of an `AbstractTimeDependentSystem`. """ function generate_initializesystem_timevarying(sys::AbstractSystem; - u0map = Dict(), - pmap = Dict(), + op = Dict(), initialization_eqs = [], guesses = Dict(), default_dd_guess = Bool(0), @@ -40,8 +39,16 @@ function generate_initializesystem_timevarying(sys::AbstractSystem; idxs_diff = isdiffeq.(eqs) # PREPROCESSING - u0map = copy(anydict(u0map)) - pmap = anydict(pmap) + op = anydict(op) + u0map = anydict() + pmap = anydict() + build_operating_point!(sys, op, u0map, pmap, defs, unknowns(sys), + parameters(sys; initial_parameters = true)) + for (k, v) in op + if has_parameter_dependency_with_lhs(sys, k) && is_variable_floatingpoint(k) + pmap[k] = v + end + end initsys_preprocessing!(u0map, defs) # 1) Use algebraic equations of system as initialization constraints @@ -65,6 +72,13 @@ function generate_initializesystem_timevarying(sys::AbstractSystem; function process_u0map_with_dummysubs(y, x) y = get(schedule.dummy_sub, y, y) y = fixpoint_sub(y, diffmap) + # FIXME: DAEs provide initial conditions that require reducing the system + # to index zero. If `isdifferential(y)`, an initial condition was given for an + # algebraic variable, so ignore it. Otherwise, the initialization system + # gets a `D(y) ~ ...` equation and errors. This is the same behavior as v9. + if isdifferential(y) + return + end # If we have `D(x) ~ x` and provide [D(x) => x, x => 1.0] to `u0map`, then # without this condition `defs` would get `x => x` instead of retaining # `x => 1.0`. @@ -170,8 +184,7 @@ $(TYPEDSIGNATURES) Generate `System` of nonlinear equations which initializes a problem from specified initial conditions of an `AbstractTimeDependentSystem`. """ function generate_initializesystem_timeindependent(sys::AbstractSystem; - u0map = Dict(), - pmap = Dict(), + op = Dict(), initialization_eqs = [], guesses = Dict(), algebraic_only = false, @@ -189,8 +202,16 @@ function generate_initializesystem_timeindependent(sys::AbstractSystem; guesses = merge(get_guesses(sys), additional_guesses) # PREPROCESSING - u0map = copy(anydict(u0map)) - pmap = anydict(pmap) + op = anydict(op) + u0map = anydict() + pmap = anydict() + build_operating_point!(sys, op, u0map, pmap, defs, unknowns(sys), + parameters(sys; initial_parameters = true)) + for (k, v) in op + if has_parameter_dependency_with_lhs(sys, k) && is_variable_floatingpoint(k) + pmap[k] = v + end + end initsys_preprocessing!(u0map, defs) # Calculate valid `Initial` parameters. These are unknowns for @@ -543,11 +564,11 @@ function SciMLBase.remake_initialization_data( defs = defaults(sys) use_scc = true initialization_eqs = Equation[] + op = anydict() if oldinitdata !== nothing && oldinitdata.metadata isa InitializationMetadata meta = oldinitdata.metadata - u0map = merge(meta.u0map, u0map) - pmap = merge(meta.pmap, pmap) + op = copy(meta.op) merge!(guesses, meta.guesses) use_scc = meta.use_scc initialization_eqs = meta.additional_initialization_eqs @@ -591,10 +612,12 @@ function SciMLBase.remake_initialization_data( if t0 === nothing && is_time_dependent(sys) t0 = 0.0 end - filter_missing_values!(u0map) - filter_missing_values!(pmap) + merge!(op, u0map, pmap) + filter_missing_values!(op) - op, missing_unknowns, missing_pars = build_operating_point!(sys, + u0map = anydict() + pmap = anydict() + missing_unknowns, missing_pars = build_operating_point!(sys, op, u0map, pmap, defs, dvs, ps) floatT = float_type_from_varmap(op) u0_constructor = p_constructor = identity @@ -607,7 +630,7 @@ function SciMLBase.remake_initialization_data( typeof(newp.initials), floatT, Val(1), Val(length(vals)), vals...) end kws = maybe_build_initialization_problem( - sys, SciMLBase.isinplace(odefn), op, u0map, pmap, t0, defs, guesses, + sys, SciMLBase.isinplace(odefn), op, t0, defs, guesses, missing_unknowns; time_dependent_init, use_scc, initialization_eqs, floatT, u0_constructor, p_constructor, allow_incomplete = true) diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl index 9448f1930e..bd77c8c519 100644 --- a/src/systems/parameter_buffer.jl +++ b/src/systems/parameter_buffer.jl @@ -27,7 +27,7 @@ This requires that `complete` has been called on the system (usually via the default behavior). """ function MTKParameters( - sys::AbstractSystem, p, u0 = Dict(); tofloat = false, + sys::AbstractSystem, op; tofloat = false, t0 = nothing, substitution_limit = 1000, floatT = nothing, p_constructor = identity) ic = if has_index_cache(sys) && get_index_cache(sys) !== nothing @@ -40,17 +40,17 @@ function MTKParameters( dvs = unknowns(sys) ps = parameters(sys; initial_parameters = true) - u0 = to_varmap(u0, dvs) - symbols_to_symbolics!(sys, u0) - p = to_varmap(p, ps) - symbols_to_symbolics!(sys, p) + op = to_varmap(op, ps) + symbols_to_symbolics!(sys, op) defs = add_toterms(recursive_unwrap(defaults(sys))) - is_time_dependent(sys) && add_observed!(sys, u0) - add_parameter_dependencies!(sys, p) + is_time_dependent(sys) && add_observed!(sys, op) + add_parameter_dependencies!(sys, op) - op, missing_unknowns, missing_pars = build_operating_point!(sys, - u0, p, defs, dvs, ps) + u0map = anydict() + pmap = anydict() + missing_unknowns, missing_pars = build_operating_point!(sys, op, + u0map, pmap, defs, dvs, ps) if t0 !== nothing op[get_iv(sys)] = t0 diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 8cff86f9d3..42379a30c1 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -335,28 +335,39 @@ Return an array of values where the `i`th element corresponds to the value of `v in `varmap`. Does not perform symbolic substitution in the values of `varmap`. Keyword arguments: -- `tofloat`: Convert values to floating point numbers using `float`. -- `container_type`: The type of container to use for the values. -- `toterm`: The `toterm` method to use for converting symbolics. -- `promotetoconcrete`: whether the promote to a concrete buffer (respecting - `tofloat`). Defaults to `container_type <: AbstractArray`. -- `check`: Error if any variables in `vars` do not have a mapping in `varmap`. Uses - [`missingvars`](@ref) to perform the check. -- `allow_symbolic` allows the returned array to contain symbolic values. If this is `true`, - `promotetoconcrete` is set to `false`. -- `is_initializeprob, guesses`: Used to determine whether the system is missing guesses. +- `container_type`: The type of the returned container. +- `allow_symbolic`: Whether the returned container of values can have symbolic expressions. +- `buffer_eltype`: The `eltype` of the returned container if `!allow_symbolic`. If + `Nothing`, automatically promotes the values in the container to a common `eltype`. +- `tofloat`: Whether to promote values to floating point numbers if + `buffer_eltype == Nothing`. +- `use_union`: Whether to allow using a `Union` as the `eltype` if + `buffer_eltype == Nothing`. +- `toterm`: The `toterm` function for canonicalizing keys of `varmap`. A value of `nothing` + disables this process. +- `check`: Whether to check if all of `vars` are keys of `varmap`. +- `is_initializeprob`: Whether an initialization problem is being constructed. Used for + better error messages. """ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector; - tofloat = true, container_type = Array, floatT = Nothing, - toterm = default_toterm, promotetoconcrete = nothing, check = true, - allow_symbolic = false, is_initializeprob = false) + tofloat = true, use_union = false, container_type = Array, buffer_eltype = Nothing, + toterm = default_toterm, check = true, allow_symbolic = false, + is_initializeprob = false) isempty(vars) && return nothing varmap = recursive_unwrap(varmap) - add_toterms!(varmap; toterm) + if toterm !== nothing + add_toterms!(varmap; toterm) + end if check missing_vars = missingvars(varmap, vars; toterm) - isempty(missing_vars) || throw(MissingVariablesError(missing_vars)) + if !isempty(missing_vars) + if is_initializeprob + throw(MissingGuessError(collect(missing_vars), collect(missing_vars))) + else + throw(MissingVariablesError(missing_vars)) + end + end end vals = map(x -> varmap[x], vars) if !allow_symbolic @@ -372,20 +383,17 @@ function better_varmap_to_vars(varmap::AbstractDict, vars::Vector; is_initializeprob ? throw(MissingGuessError(missingsyms, missingvals)) : throw(UnexpectedSymbolicValueInVarmap(missingsyms[1], missingvals[1])) end - if tofloat && !(floatT == Nothing) - vals = floatT.(vals) + if buffer_eltype == Nothing + vals = promote_to_concrete(vals; tofloat, use_union) + else + vals = buffer_eltype.(vals) end end - if container_type <: Union{AbstractDict, Tuple, Nothing, SciMLBase.NullParameters} + if container_type <: Union{AbstractDict, Nothing, SciMLBase.NullParameters} container_type = Array end - promotetoconcrete === nothing && (promotetoconcrete = container_type <: AbstractArray) - if promotetoconcrete && !allow_symbolic - vals = promote_to_concrete(vals; tofloat = tofloat, use_union = false) - end - if isempty(vals) return nothing elseif container_type <: Tuple @@ -540,8 +548,8 @@ Also updates `u0map` and `pmap` in-place to contain all the initial conditions i by unknowns and parameters respectively. """ function build_operating_point!(sys::AbstractSystem, - u0map::AbstractDict, pmap::AbstractDict, defs::AbstractDict, dvs, ps) - op = add_toterms(u0map) + op::AbstractDict, u0map::AbstractDict, pmap::AbstractDict, defs::AbstractDict, dvs, ps) + add_toterms!(op) missing_unknowns = add_fallbacks!(op, dvs, defs) for (k, v) in defs haskey(op, k) && continue @@ -594,7 +602,7 @@ function build_operating_point!(sys::AbstractSystem, pmap[k] = v end - return op, missing_unknowns, missing_pars + return missing_unknowns, missing_pars end """ @@ -890,13 +898,9 @@ $(TYPEDFIELDS) """ struct InitializationMetadata{R <: ReconstructInitializeprob, GUU, SIU} """ - The `u0map` used to construct the initialization. + The operating point used to construct the initialization. """ - u0map::Dict{Any, Any} - """ - The `pmap` used to construct the initialization. - """ - pmap::Dict{Any, Any} + op::Dict{Any, Any} """ The `guesses` used to construct the initialization. """ @@ -1011,14 +1015,14 @@ end $(TYPEDSIGNATURES) Build and return the initialization problem and associated data as a `NamedTuple` to be passed -to the `SciMLFunction` constructor. Requires the system `sys`, operating point `op`, -user-provided `u0map` and `pmap`, initial time `t`, system defaults `defs`, user-provided -`guesses`, and list of unknowns which don't have a value in `op`. The keyword `implicit_dae` -denotes whether the `SciMLProblem` being constructed is in implicit DAE form (`DAEProblem`). -All other keyword arguments are forwarded to `InitializationProblem`. +to the `SciMLFunction` constructor. Requires the system `sys`, operating point `op`, initial +time `t`, system defaults `defs`, user-provided `guesses`, and list of unknowns which don't +have a value in `op`. The keyword `implicit_dae` denotes whether the `SciMLProblem` being +constructed is in implicit DAE form (`DAEProblem`). All other keyword arguments are forwarded +to `InitializationProblem`. """ function maybe_build_initialization_problem( - sys::AbstractSystem, iip, op::AbstractDict, u0map, pmap, t, defs, + sys::AbstractSystem, iip, op::AbstractDict, t, defs, guesses, missing_unknowns; implicit_dae = false, time_dependent_init = is_time_dependent(sys), u0_constructor = identity, p_constructor = identity, floatT = Float64, initialization_eqs = [], @@ -1030,7 +1034,7 @@ function maybe_build_initialization_problem( end initializeprob = ModelingToolkit.InitializationProblem{iip}( - sys, t, u0map, pmap; guesses, time_dependent_init, initialization_eqs, + sys, t, op; guesses, time_dependent_init, initialization_eqs, use_scc, u0_constructor, p_constructor, kwargs...) if state_values(initializeprob) !== nothing _u0 = state_values(initializeprob) @@ -1064,7 +1068,7 @@ function maybe_build_initialization_problem( nothing end meta = InitializationMetadata( - u0map, pmap, guesses, Vector{Equation}(initialization_eqs), + copy(op), copy(guesses), Vector{Equation}(initialization_eqs), use_scc, time_dependent_init, ReconstructInitializeprob( sys, initializeprob.f.sys; u0_constructor, p_constructor), @@ -1100,7 +1104,7 @@ function maybe_build_initialization_problem( end for p in punknowns - is_parameter_solvable(p, pmap, defs, guesses) || continue + is_parameter_solvable(p, op, defs, guesses) || continue get(op, p, missing) === missing || continue p = unwrap(p) op[p] = getu(initializeprob, p)(initializeprob) @@ -1180,8 +1184,8 @@ Keyword arguments: - `build_initializeprob`: If `false`, avoids building the initialization problem. - `t`: The initial time of the `ODEProblem`. If this is not provided, the initialization problem cannot be built. -- `implicit_dae`: Also build a mapping of derivatives of states to values for implicit DAEs, - using `du0map`. Changes the return value of this function to `(f, du0, u0, p)` instead of +- `implicit_dae`: Also build a mapping of derivatives of states to values for implicit DAEs. + Changes the return value of this function to `(f, du0, u0, p)` instead of `(f, u0, p)`. - `guesses`: The guesses for variables in the system, used as initial values for the initialization problem. @@ -1193,11 +1197,13 @@ Keyword arguments: - `fully_determined`: Override whether the initialization system is fully determined. - `check_initialization_units`: Enable or disable unit checks when constructing the initialization problem. -- `tofloat`, `is_initializeprob`: Passed to [`better_varmap_to_vars`](@ref) for building `u0` (and possibly `p`). +- `tofloat`: Passed to [`better_varmap_to_vars`](@ref) when building the parameter vector of + a non-split system. +- `u0_eltype`: The `eltype` of the `u0` vector. If `nothing`, finds the promoted floating point + type from `op`. - `u0_constructor`: A function to apply to the `u0` value returned from `better_varmap_to_vars` to construct the final `u0` value. - `p_constructor`: A function to apply to each array buffer created when constructing the parameter object. -- `du0map`: A map of derivatives to values. See `implicit_dae`. - `check_length`: Whether to check the number of equations along with number of unknowns and length of `u0` vector for consistency. If `false`, do not check with equations. This is forwarded to `check_eqs_u0` @@ -1220,13 +1226,13 @@ Keyword arguments: All other keyword arguments are passed as-is to `constructor`. """ function process_SciMLProblem( - constructor, sys::AbstractSystem, u0map, pmap; + constructor, sys::AbstractSystem, op; build_initializeprob = supports_initialization(sys), implicit_dae = false, t = nothing, guesses = AnyDict(), warn_initialize_determined = true, initialization_eqs = [], eval_expression = false, eval_module = @__MODULE__, fully_determined = nothing, - check_initialization_units = false, tofloat = true, - u0_constructor = identity, p_constructor = identity, du0map = nothing, + check_initialization_units = false, u0_eltype = nothing, tofloat = true, + u0_constructor = identity, p_constructor = identity, check_length = true, symbolic_u0 = false, warn_cyclic_dependency = false, circular_dependency_max_cycle_length = length(all_symbols(sys)), circular_dependency_max_cycles = 10, @@ -1240,15 +1246,12 @@ function process_SciMLProblem( check_array_equations_unknowns(eqs, dvs) - u0Type = typeof(u0map) - pType = typeof(pmap) + u0Type = pType = typeof(op) - u0map = to_varmap(u0map, dvs) - symbols_to_symbolics!(sys, u0map) - pmap = to_varmap(pmap, parameters(sys)) - symbols_to_symbolics!(sys, pmap) + op = to_varmap(op, dvs) + symbols_to_symbolics!(sys, op) - check_inputmap_keys(sys, u0map, pmap) + check_inputmap_keys(sys, op) defs = add_toterms(recursive_unwrap(defaults(sys))) kwargs = NamedTuple(kwargs) @@ -1259,7 +1262,9 @@ function process_SciMLProblem( obs, _ = unhack_observed(observed(sys), Equation[x for x in eqs if x isa Equation]) end - op, missing_unknowns, missing_pars = build_operating_point!(sys, + u0map = anydict() + pmap = anydict() + missing_unknowns, missing_pars = build_operating_point!(sys, op, u0map, pmap, defs, dvs, ps) floatT = Bool @@ -1269,8 +1274,10 @@ function process_SciMLProblem( floatT = float_type_from_varmap(op, floatT) end + u0_eltype = something(u0_eltype, floatT) + if !is_time_dependent(sys) || is_initializesystem(sys) - add_observed_equations!(u0map, obs) + add_observed_equations!(op, obs) end if u0_constructor === identity && u0Type <: StaticArray u0_constructor = vals -> SymbolicUtils.Code.create_array( @@ -1284,7 +1291,7 @@ function process_SciMLProblem( if build_initializeprob kws = maybe_build_initialization_problem( sys, constructor <: SciMLBase.AbstractSciMLFunction{true}, - op, u0map, pmap, t, defs, guesses, missing_unknowns; + op, t, defs, guesses, missing_unknowns; implicit_dae, warn_initialize_determined, initialization_eqs, eval_expression, eval_module, fully_determined, warn_cyclic_dependency, check_units = check_initialization_units, @@ -1318,7 +1325,7 @@ function process_SciMLProblem( evaluate_varmap!(op, dvs; limit = substitution_limit) u0 = better_varmap_to_vars( - op, dvs; tofloat, floatT, + op, dvs; buffer_eltype = u0_eltype, container_type = u0Type, allow_symbolic = symbolic_u0, is_initializeprob) if u0 !== nothing @@ -1351,11 +1358,9 @@ function process_SciMLProblem( p = p_constructor(better_varmap_to_vars(op, ps; tofloat, container_type = pType)) end - if implicit_dae && du0map !== nothing + if implicit_dae ddvs = map(Differential(iv), dvs) - du0map = to_varmap(du0map, ddvs) - merge!(op, du0map) - du0 = varmap_to_vars(op, ddvs; toterm = identity, + du0 = varmap_to_vars(op, ddvs; toterm = default_toterm, tofloat) kwargs = merge(kwargs, (; ddvs)) else @@ -1388,22 +1393,17 @@ end # Check that the keys of a u0map or pmap are valid # (i.e. are symbolic keys, and are defined for the system.) -function check_inputmap_keys(sys, u0map, pmap) +function check_inputmap_keys(sys, op) badvarkeys = Any[] - for k in keys(u0map) + for k in keys(op) if symbolic_type(k) === NotSymbolic() push!(badvarkeys, k) end end - badparamkeys = Any[] - for k in keys(pmap) - if symbolic_type(k) === NotSymbolic() - push!(badparamkeys, k) - end + if !isempty(badvarkeys) + throw(InvalidKeyError(collect(badvarkeys))) end - (isempty(badvarkeys) && isempty(badparamkeys)) || - throw(InvalidKeyError(collect(badvarkeys), collect(badparamkeys))) end const BAD_KEY_MESSAGE = """ @@ -1413,13 +1413,11 @@ const BAD_KEY_MESSAGE = """ struct InvalidKeyError <: Exception vars::Any - params::Any end function Base.showerror(io::IO, e::InvalidKeyError) println(io, BAD_KEY_MESSAGE) - println(io, "u0map: $(join(e.vars, ", "))") - println(io, "pmap: $(join(e.params, ", "))") + println(io, join(e.vars, ", ")) end function SciMLBase.detect_cycles(sys::AbstractSystem, varmap::Dict{Any, Any}, vars) @@ -1566,7 +1564,7 @@ macro fallback_iip_specialize(ex) fn_sarr = nothing if occursin("Problem", string(fnname_name)) # args should at least contain an argument for the `u0map` - @assert length(args) > 3 + @assert length(args) > 2 u0_arg = args[3] # should not have a type-annotation @assert !Meta.isexpr(u0_arg, :(::)) diff --git a/test/basic_transformations.jl b/test/basic_transformations.jl index 661b508861..6880f77479 100644 --- a/test/basic_transformations.jl +++ b/test/basic_transformations.jl @@ -14,14 +14,14 @@ D = Differential(t) u0 = [x => 1.0, y => 1.0] p = [α => 1.5, β => 1.0, δ => 3.0, γ => 1.0] tspan = (0.0, 10.0) - prob = ODEProblem(sys, u0, tspan, p) + prob = ODEProblem(sys, [u0; p], tspan) sol = solve(prob, Tsit5()) sys2 = liouville_transform(sys) sys2 = complete(sys2) u0 = [x => 1.0, y => 1.0, sys2.trJ => 1.0] - prob = ODEProblem(sys2, u0, tspan, p, jac = true) + prob = ODEProblem(sys2, [u0; p], tspan, jac = true) sol = solve(prob, Tsit5()) @test sol[end, end] ≈ 1.0742818931017244 end @@ -53,8 +53,8 @@ end M1 = mtkcompile(M1; allow_symbolic = true) M2 = mtkcompile(M2; allow_symbolic = true) - prob1 = ODEProblem(M1, [M1.s => 1.0], (1.0, 4.0), []) - prob2 = ODEProblem(M2, [], (1.0, 2.0), []) + prob1 = ODEProblem(M1, [M1.s => 1.0], (1.0, 4.0)) + prob2 = ODEProblem(M2, [], (1.0, 2.0)) sol1 = solve(prob1, Tsit5(); reltol = 1e-10, abstol = 1e-10) sol2 = solve(prob2, Tsit5(); reltol = 1e-10, abstol = 1e-10) ts = range(0.0, 1.0, length = 50) @@ -125,7 +125,7 @@ end Dx = Differential(Mx.x) u0 = [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0] p = [v => 10.0] - prob = ODEProblem(Mx, u0, (0.0, 20.0), p) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities + prob = ODEProblem(Mx, [u0; p], (0.0, 20.0)) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities sol = solve(prob, Tsit5(); reltol = 1e-5) @test all(isapprox.(sol[Mx.y], sol[Mx.x - g * (Mx.t)^2 / 2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2) end @@ -138,7 +138,7 @@ end Mx = mtkcompile(Mx; allow_symbolic = true) Dx = Differential(Mx.x) u0 = [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t => 0.0, Mx.xˍt => 10.0] - prob = ODEProblem(Mx, u0, (0.0, 20.0), []) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities + prob = ODEProblem(Mx, u0, (0.0, 20.0)) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities sol = solve(prob, Tsit5(); reltol = 1e-5) @test all(isapprox.(sol[Mx.y], sol[Mx.x - g * (Mx.t)^2 / 2]; atol = 1e-10)) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2) end @@ -200,7 +200,7 @@ end _f = LinearInterpolation([1.0, 1.0], [-100.0, +100.0]) # constant value 1 M2s = mtkcompile(M2; allow_symbolic = true) - prob = ODEProblem(M2s, [M2s.y => 0.0], (1.0, 4.0), [fc => _f, f => _f]) + prob = ODEProblem(M2s, [M2s.y => 0.0, fc => _f, f => _f], (1.0, 4.0)) sol = solve(prob, Tsit5(); abstol = 1e-5) @test isapprox(sol(4.0, idxs = M2.y), 12.0; atol = 1e-5) # Anal solution is D(y) ~ 12 => y(t) ~ 12*t + C => y(x) ~ 12*√(x) + C. With y(x=1)=0 => 12*(√(x)-1), so y(x=4) ~ 12 end @@ -227,7 +227,7 @@ end Mx = mtkcompile(Mx; allow_symbolic = true) Dx = Differential(Mx.x) u0 = [Mx.y => 0.0, Dx(Mx.y) => 1.0, Mx.t_units => 0.0, Mx.xˍt_units => 10.0] - prob = ODEProblem(Mx, u0, (0.0, 20.0), []) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities + prob = ODEProblem(Mx, u0, (0.0, 20.0)) # 1 = dy/dx = (dy/dt)/(dx/dt) means equal initial horizontal and vertical velocities sol = solve(prob, Tsit5(); reltol = 1e-5) # compare to analytical solution (x(t) = v*t, y(t) = v*t - g*t^2/2) @test all(isapprox.(sol[Mx.y], sol[Mx.x - g * (Mx.t_units)^2 / 2]; atol = 1e-10)) @@ -300,7 +300,7 @@ end new_sys = change_independent_variable(sys, sys.x; add_old_diff = true) ss_new_sys = mtkcompile(new_sys; allow_symbolic = true) u0 = [new_sys.y => 0.5, new_sys.t => 0.0] - prob = ODEProblem(ss_new_sys, u0, (0.0, 0.5), []) + prob = ODEProblem(ss_new_sys, u0, (0.0, 0.5)) sol = solve(prob, Tsit5(); reltol = 1e-5) @test sol[new_sys.y][end] ≈ 0.75 end diff --git a/test/bvproblem.jl b/test/bvproblem.jl index bafb3e9674..abc925ccfa 100644 --- a/test/bvproblem.jl +++ b/test/bvproblem.jl @@ -21,11 +21,11 @@ daesolvers = [Ascher2, Ascher4, Ascher6] tspan = (0.0, 10.0) @mtkcompile lotkavolterra = System(eqs, t) - op = ODEProblem(lotkavolterra, u0map, tspan, parammap) + op = ODEProblem(lotkavolterra, [u0map; parammap], tspan) osol = solve(op, Vern9()) bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( - lotkavolterra, u0map, tspan, parammap) + lotkavolterra, [u0map; parammap], tspan) for solver in solvers sol = solve(bvp, solver(), dt = 0.01) @@ -35,7 +35,7 @@ daesolvers = [Ascher2, Ascher4, Ascher6] # Test out of place bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}( - lotkavolterra, u0map, tspan, parammap) + lotkavolterra, [u0map; parammap], tspan) for solver in solvers sol = solve(bvp2, solver(), dt = 0.01) @@ -58,10 +58,11 @@ end parammap = [:L => 1.0, :g => 9.81] tspan = (0.0, 6.0) - op = ODEProblem(pend, u0map, tspan, parammap) + op = ODEProblem(pend, [u0map; parammap], tspan) osol = solve(op, Vern9()) - bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + pend, [u0map; parammap], tspan) for solver in solvers sol = solve(bvp, solver(), dt = 0.01) @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) @@ -70,7 +71,7 @@ end # Test out-of-place bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}( - pend, u0map, tspan, parammap) + pend, [u0map; parammap], tspan) for solver in solvers sol = solve(bvp2, solver(), dt = 0.01) @@ -289,8 +290,8 @@ end @mtkcompile lksys = System(eqs, t; costs, consolidate) @test_throws ModelingToolkit.SystemCompatibilityError ODEProblem( - lksys, u0map, tspan, parammap) - prob = ODEProblem(lksys, u0map, tspan, parammap; check_compatibility = false) + lksys, [u0map; parammap], tspan) + prob = ODEProblem(lksys, [u0map; parammap], tspan; check_compatibility = false) sol = solve(prob, Tsit5()) costfn = ModelingToolkit.generate_cost( lksys; expression = Val{false}, wrap_gfw = Val{true}) @@ -304,7 +305,7 @@ end @mtkcompile lksys = System(eqs, t; costs, consolidate) @test t_c ∈ Set(parameters(lksys)) push!(parammap, t_c => 0.56) - prob = ODEProblem(lksys, u0map, tspan, parammap; check_compatibility = false) + prob = ODEProblem(lksys, [u0map; parammap], tspan; check_compatibility = false) sol = solve(prob, Tsit5()) costfn = ModelingToolkit.generate_cost( lksys; expression = Val{false}, wrap_gfw = Val{true}) diff --git a/test/clock.jl b/test/clock.jl index 7fc9e7ddd2..a50026b38f 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -120,22 +120,22 @@ eqs = [yd ~ Sample(dt)(y) @test_skip begin Tf = 1.0 - prob = ODEProblem(ss, [x => 0.1], (0.0, Tf), - [kp => 1.0; ud(k - 1) => 2.1; ud(k - 2) => 2.0]) + prob = ODEProblem( + ss, [x => 0.1, kp => 1.0; ud(k - 1) => 2.1; ud(k - 2) => 2.0], (0.0, Tf)) # 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.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(k - 1) => 2.1; ud(k - 2) => 2.0]) # recreate problem to empty saved values + prob = ODEProblem( + ss, [x => 0.1, kp => 1.0; ud(k - 1) => 2.1; ud(k - 2) => 2.0], (0.0, Tf)) # recreate problem to empty saved values sol = solve(prob, Tsit5(), kwargshandle = KeywordArgSilent) ss_nosplit = mtkcompile(sys; split = false) - prob_nosplit = ODEProblem(ss_nosplit, [x => 0.1], (0.0, Tf), - [kp => 1.0; ud(k - 1) => 2.1; ud(k - 2) => 2.0]) + prob_nosplit = ODEProblem( + ss_nosplit, [x => 0.1, kp => 1.0; ud(k - 1) => 2.1; ud(k - 2) => 2.0], (0.0, Tf)) int = init(prob_nosplit, Tsit5(); kwargshandle = KeywordArgSilent) @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(k - 1) => 2.1; ud(k - 2) => 2.0]) # recreate problem to empty saved values + prob_nosplit = ODEProblem( + ss_nosplit, [x => 0.1, kp => 1.0; ud(k - 1) => 2.1; ud(k - 2) => 2.0], (0.0, Tf)) # 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. @@ -299,8 +299,8 @@ eqs = [yd ~ Sample(dt)(y) ss_nosplit = mtkcompile(cl; split = false) if VERSION >= v"1.7" - prob = ODEProblem(ss, [x => 0.0], (0.0, 1.0), [kp => 1.0]) - prob_nosplit = ODEProblem(ss_nosplit, [x => 0.0], (0.0, 1.0), [kp => 1.0]) + prob = ODEProblem(ss, [x => 0.0, kp => 1.0], (0.0, 1.0)) + prob_nosplit = ODEProblem(ss_nosplit, [x => 0.0, kp => 1.0], (0.0, 1.0)) sol = solve(prob, Tsit5(), kwargshandle = KeywordArgSilent) sol_nosplit = solve(prob_nosplit, Tsit5(), kwargshandle = KeywordArgSilent) @@ -535,8 +535,8 @@ eqs = [yd ~ Sample(dt)(y) @test int.ps[x] == 2.0 @test int.ps[x(k - 1)] == 1.0 - @test_throws ErrorException ODEProblem(sys, [], (0.0, 10.0), [x => 2.0]) - prob = ODEProblem(sys, [], (0.0, 10.0), [x(k - 1) => 2.0]) + @test_throws ErrorException ODEProblem(sys, [x => 2.0], (0.0, 10.0)) + prob = ODEProblem(sys, [x(k - 1) => 2.0], (0.0, 10.0)) int = init(prob, Tsit5(); kwargshandle = KeywordArgSilent) @test int.ps[x] == 3.0 @test int.ps[x(k - 1)] == 2.0 diff --git a/test/code_generation.jl b/test/code_generation.jl index 9a3805ce21..2d5925cbca 100644 --- a/test/code_generation.jl +++ b/test/code_generation.jl @@ -59,7 +59,7 @@ end @variables x(t) @parameters p[0:2] (f::Function)(..) @mtkcompile sys = System(D(x) ~ p[0] * x + p[1] * t + p[2] + f(p), t) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => [1.0, 2.0, 3.0], f => sum]) + prob = ODEProblem(sys, [x => 1.0, p => [1.0, 2.0, 3.0], f => sum], (0.0, 1.0)) @test prob.ps[p] == [1.0, 2.0, 3.0] @test prob.ps[p[0]] == 1.0 sol = solve(prob, Tsit5()) @@ -72,8 +72,9 @@ end [D(x[0]) ~ p[1] * x[0] + x[2], D(x[1]) ~ p[2] * f(x) + x[2]], t) sys = mtkcompile(sys, inputs = [x[2]], outputs = []) @test is_parameter(sys, x[2]) - prob = ODEProblem(sys, [x[0] => 1.0, x[1] => 1.0], (0.0, 1.0), - [p => ones(2), f => sum, x[2] => 2.0]) + prob = ODEProblem( + sys, [x[0] => 1.0, x[1] => 1.0, x[2] => 2.0, p => ones(2), f => sum], + (0.0, 1.0)) sol = solve(prob, Tsit5()) @test SciMLBase.successful_retcode(sol) end diff --git a/test/components.jl b/test/components.jl index 689173d61f..e71ef3fa45 100644 --- a/test/components.jl +++ b/test/components.jl @@ -119,15 +119,16 @@ include("common/serial_inductor.jl") u0 = unknowns(sys) .=> 0 @test_nowarn ODEProblem( sys, [], (0, 10.0), guesses = u0, warn_initialize_determined = false) - prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0, [], (0, 0.5), guesses = u0) + prob = DAEProblem(sys, D.(unknowns(sys)) .=> 0, (0, 0.5), guesses = u0) sol = solve(prob, DFBDF()) @test sol.retcode == SciMLBase.ReturnCode.Success sys2 = mtkcompile(ll2_model) @test length(equations(sys2)) == 3 - u0 = unknowns(sys) .=> 0 + u0 = [sys.inductor2.i => 0] prob = ODEProblem(sys, u0, (0, 10.0)) - @test_nowarn sol = solve(prob, FBDF()) + sol = solve(prob, FBDF()) + @test SciMLBase.successful_retcode(sol) end @testset "Compose/extend" begin diff --git a/test/constants.jl b/test/constants.jl index ec98dab9b4..ce5c7e6e8e 100644 --- a/test/constants.jl +++ b/test/constants.jl @@ -12,7 +12,7 @@ UMT = ModelingToolkit.UnitfulUnitCheck D = Differential(t) eqs = [D(x) ~ a] @named sys = System(eqs, t) -prob = ODEProblem(complete(sys), [0], [0.0, 1.0], []) +prob = ODEProblem(complete(sys), [0], [0.0, 1.0]) sol = solve(prob, Tsit5()) # Test mtkcompile substitutions & observed values @@ -43,7 +43,7 @@ simp = mtkcompile(sys) @mtkcompile fol_model = System(eqs, MT.t_nounits) - prob = ODEProblem(fol_model, [], (0.0, 10.0), [h => 1]) + prob = ODEProblem(fol_model, [h => 1], (0.0, 10.0)) @test prob[x] ≈ 1 @test prob.ps[τ] ≈ 0.5 end diff --git a/test/dae_jacobian.jl b/test/dae_jacobian.jl index 94f15cbb7c..17b14b39e4 100644 --- a/test/dae_jacobian.jl +++ b/test/dae_jacobian.jl @@ -43,12 +43,12 @@ u0 = [u1 => 1.0, tspan = (0.0, 10.0) -du0 = [0.5, -2.0] +du0 = [D(u1) => 0.5, D(u2) => -2.0] p = [p1 => 1.5, p2 => 3.0] -prob = DAEProblem(complete(sys), du0, u0, tspan, p, jac = true, sparse = true) +prob = DAEProblem(complete(sys), [du0; u0; p], tspan, jac = true, sparse = true) sol = solve(prob, IDA(linear_solver = :KLU)) @test maximum(sol - sol1) < 1e-12 diff --git a/test/discrete_system.jl b/test/discrete_system.jl index 2db2aab8d3..dc0281c8bf 100644 --- a/test/discrete_system.jl +++ b/test/discrete_system.jl @@ -51,7 +51,7 @@ reorderer = getu(syss, [S, I, R]) u0 = [S(k - 1) => 990.0, I(k - 1) => 10.0, R(k - 1) => 0.0] p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1, nsteps => 400] tspan = (0.0, ModelingToolkit.value(substitute(nsteps, p))) # value function (from Symbolics) is used to convert a Num to Float64 -prob_map = DiscreteProblem(syss, u0, tspan, p) +prob_map = DiscreteProblem(syss, [u0; p], tspan) @test prob_map.f.sys === syss # Solution diff --git a/test/distributed.jl b/test/distributed.jl index e7edc17a66..8c3ca4dcfa 100644 --- a/test/distributed.jl +++ b/test/distributed.jl @@ -16,10 +16,10 @@ addprocs(2) @everywhere @named de = System(eqs, t) @everywhere de = complete(de) -@everywhere u0 = [19.0, 20.0, 50.0] -@everywhere params = [16.0, 45.92, 4] +@everywhere u0 = unknowns(de) .=> [19.0, 20.0, 50.0] +@everywhere params = parameters(de) .=> [16.0, 45.92, 4] -@everywhere ode_prob = ODEProblem(de, u0, (0.0, 10.0), params) +@everywhere ode_prob = ODEProblem(de, [u0; params], (0.0, 10.0)) @everywhere begin using OrdinaryDiffEq diff --git a/test/dq_units.jl b/test/dq_units.jl index 74209c3584..615de1d642 100644 --- a/test/dq_units.jl +++ b/test/dq_units.jl @@ -216,7 +216,7 @@ end p = [pend.g => 1.0, pend.L => 1.0] guess = [pend.λ => 0.0] @test prob = ODEProblem( - pend, u0, (0.0, 1.0), p; guesses = guess, check_units = false) isa Any + pend, [u0; p], (0.0, 1.0); guesses = guess, check_units = false) isa Any end @parameters p [unit = u"L/s"] d [unit = u"s^(-1)"] diff --git a/test/extensions/dynamic_optimization.jl b/test/extensions/dynamic_optimization.jl index b9a512707c..62553121fe 100644 --- a/test/extensions/dynamic_optimization.jl +++ b/test/extensions/dynamic_optimization.jl @@ -29,7 +29,7 @@ const M = ModelingToolkit jprob = JuMPDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01) @test JuMP.num_constraints(jprob.model) == 2 # initials jsol = solve(jprob, Ipopt.Optimizer, constructRK4, silent = true) - oprob = ODEProblem(sys, u0map, tspan, parammap) + oprob = ODEProblem(sys, [u0map; parammap], tspan) osol = solve(oprob, SimpleRK4(), dt = 0.01) cprob = CasADiDynamicOptProblem(sys, u0map, tspan, parammap, dt = 0.01) csol = solve(cprob, "ipopt", constructRK4) @@ -141,7 +141,7 @@ end @parameters (u_interp::ConstantInterpolation)(..) @mtkcompile block_ode = System([D(x(t)) ~ v(t), D(v(t)) ~ u_interp(t)], t) spline = ctrl_to_spline(jsol.input_sol, ConstantInterpolation) - oprob = ODEProblem(block_ode, u0map, tspan, [u_interp => spline]) + oprob = ODEProblem(block_ode, [u0map; [u_interp => spline]], tspan) osol = solve(oprob, Vern8(), dt = 0.01, adaptive = false) @test ≈(jsol.sol.u, osol.u, rtol = 0.05) @test ≈(csol.sol.u, osol.u, rtol = 0.05) @@ -185,10 +185,9 @@ end D(q(t)) ~ -ν * q(t) + c * (1 - α_interp(t)) * s * w(t)] @mtkcompile beesys_ode = System(eqs, t) oprob = ODEProblem(beesys_ode, - u0map, - tspan, - merge(Dict(pmap), - Dict(α_interp => ctrl_to_spline(jsol.input_sol, LinearInterpolation)))) + merge(Dict(u0map), Dict(pmap), + Dict(α_interp => ctrl_to_spline(jsol.input_sol, LinearInterpolation))), + tspan) osol = solve(oprob, Tsit5(); dt = 0.01, adaptive = false) @test ≈(osol.u, jsol.sol.u, rtol = 0.01) @test ≈(osol.u, csol.sol.u, rtol = 0.01) @@ -239,13 +238,13 @@ end D(m(t)) ~ -T_interp(t) / c] @mtkcompile rocket_ode = System(eqs, t) interpmap = Dict(T_interp => ctrl_to_spline(jsol.input_sol, CubicSpline)) - oprob = ODEProblem(rocket_ode, u0map, (ts, te), merge(Dict(pmap), interpmap)) + oprob = ODEProblem(rocket_ode, merge(Dict(u0map), Dict(pmap), interpmap), (ts, te)) osol = solve(oprob, RadauIIA5(); adaptive = false, dt = 0.001) @test ≈(jsol.sol.u, osol.u, rtol = 0.02) @test ≈(csol.sol.u, osol.u, rtol = 0.02) interpmap1 = Dict(T_interp => ctrl_to_spline(isol.input_sol, CubicSpline)) - oprob1 = ODEProblem(rocket_ode, u0map, (ts, te), merge(Dict(pmap), interpmap1)) + oprob1 = ODEProblem(rocket_ode, merge(Dict(u0map), Dict(pmap), interpmap1), (ts, te)) osol1 = solve(oprob1, ImplicitEuler(); adaptive = false, dt = 0.001) @test ≈(isol.sol.u, osol1.u, rtol = 0.01) end diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index f9abeaa449..a4af57c8fe 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -61,8 +61,8 @@ end 0 ~ y * z + 4x^2 + wrapper(r)] @mtkcompile sys = System(eqs) - prob = HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0], - [p => 2.0, q => 4, r => Wrapper([1.0 1.0; 0.0 0.0])]) + prob = HomotopyContinuationProblem(sys, + [x => 1.0, y => 1.0, z => 1.0, p => 2.0, q => 4, r => Wrapper([1.0 1.0; 0.0 0.0])]) @test prob.ps[p] == 2.0 @test prob.ps[q] == 4 @test prob.ps[r].x == [1.0 1.0; 0.0 0.0] @@ -79,7 +79,7 @@ end _x = collect(x) eqs = collect(0 .~ vec(sum(_x * _x'; dims = 2)) + collect(p)) @mtkcompile sys = System(eqs) - prob = HomotopyContinuationProblem(sys, [x => ones(3)], [p => 1:3]) + prob = HomotopyContinuationProblem(sys, [x => ones(3), p => 1:3]) @test prob[x] == ones(3) @test prob[p + x] == [2, 3, 4] prob[x] = 2ones(3) diff --git a/test/fmi/fmi.jl b/test/fmi/fmi.jl index 29b2dc1224..e0eaf373d0 100644 --- a/test/fmi/fmi.jl +++ b/test/fmi/fmi.jl @@ -127,8 +127,9 @@ end systems = [adder]) # c will be solved for by initialization # this tests that initialization also works with FMUs - prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0], - (0.0, 1.0), [sys.adder.value => 2.0]) + prob = ODEProblem( + sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0, sys.adder.value => 2.0], + (0.0, 1.0)) return sys, prob end @@ -182,7 +183,7 @@ end ) prob = ODEProblem( - sys, [sys.sspace.x => 1.0], (0.0, 1.0), [sys.sspace.A => 2.0]; use_scc = false) + sys, [sys.sspace.x => 1.0, sys.sspace.A => 2.0], (0.0, 1.0); use_scc = false) return sys, prob end diff --git a/test/guess_propagation.jl b/test/guess_propagation.jl index 0be9506bb9..62a6ae9871 100644 --- a/test/guess_propagation.jl +++ b/test/guess_propagation.jl @@ -13,7 +13,7 @@ initialization_eqs = [1 ~ exp(1 + x)] @named sys = System(eqs, t; initialization_eqs) sys = complete(mtkcompile(sys)) tspan = (0.0, 0.2) -prob = ODEProblem(sys, [], tspan, []) +prob = ODEProblem(sys, [], tspan) @test prob.f.initializeprob[y] == 2.0 @test prob.f.initializeprob[x] == 2.0 @@ -30,7 +30,7 @@ initialization_eqs = [1 ~ exp(1 + x)] @named sys = System(eqs, t; initialization_eqs) sys = complete(mtkcompile(sys)) tspan = (0.0, 0.2) -prob = ODEProblem(sys, [], tspan, []) +prob = ODEProblem(sys, [], tspan) @test prob.f.initializeprob[x] == 2.0 @test prob.f.initializeprob[y] == 2.0 @@ -49,7 +49,7 @@ initialization_eqs = [1 ~ exp(1 + x)] sys = complete(mtkcompile(sys)) tspan = (0.0, 0.2) -prob = ODEProblem(sys, [], tspan, []) +prob = ODEProblem(sys, [], tspan) @test prob.f.initializeprob[x] == -1.0 sol = solve(prob.f.initializeprob; show_trace = Val(true)) @@ -69,7 +69,7 @@ initialization_eqs = [1 ~ exp(1 + x)] sys = complete(mtkcompile(sys)) tspan = (0.0, 0.2) -prob = ODEProblem(sys, [], tspan, []) +prob = ODEProblem(sys, [], tspan) @test prob.f.initializeprob[x] == -1.0 sol = solve(prob.f.initializeprob; show_trace = Val(true)) @@ -81,7 +81,7 @@ sol = solve(prob.f.initializeprob; show_trace = Val(true)) @variables x(t) @variables y(t) = x @mtkcompile sys = System([x ~ x0, D(y) ~ x], t) -prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) +prob = ODEProblem(sys, [x0 => 1.0], (0.0, 1.0)) @test prob[x] == 1.0 @test prob[y] == 1.0 @@ -89,7 +89,7 @@ prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) @variables x(t) @variables y(t) = x0 @mtkcompile sys = System([x ~ x0, D(y) ~ x], t) -prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) +prob = ODEProblem(sys, [x0 => 1.0], (0.0, 1.0)) @test prob[x] == 1.0 @test prob[y] == 1.0 @@ -97,7 +97,7 @@ prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) @variables x(t) @variables y(t) = x0 @mtkcompile sys = System([x ~ y, D(y) ~ x], t) -prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) +prob = ODEProblem(sys, [x0 => 1.0], (0.0, 1.0)) @test prob[x] == 1.0 @test prob[y] == 1.0 @@ -105,6 +105,6 @@ prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) @variables x(t) = x0 @variables y(t) = x @mtkcompile sys = System([x ~ y, D(y) ~ x], t) -prob = ODEProblem(sys, [], (0.0, 1.0), [x0 => 1.0]) +prob = ODEProblem(sys, [x0 => 1.0], (0.0, 1.0)) @test prob[x] == 1.0 @test prob[y] == 1.0 diff --git a/test/implicit_discrete_system.jl b/test/implicit_discrete_system.jl index 59e4a94977..bbfb045b1f 100644 --- a/test/implicit_discrete_system.jl +++ b/test/implicit_discrete_system.jl @@ -28,7 +28,7 @@ rng = StableRNG(22525) @test prob.u0 == [1.0, 1.0] @variables x(t) @mtkcompile sys = System([x(k) ~ x(k) * x(k - 1) - 3], t) - @test_throws ModelingToolkit.MissingVariablesError prob=ImplicitDiscreteProblem( + @test_throws ModelingToolkit.MissingGuessError prob=ImplicitDiscreteProblem( sys, [], tspan) end diff --git a/test/initial_values.jl b/test/initial_values.jl index 347f794e8d..ef6c404440 100644 --- a/test/initial_values.jl +++ b/test/initial_values.jl @@ -64,7 +64,7 @@ tspan = (0.0, 1.0) ps = [k1 => 1.0, k2 => 5.0] # Broken since we need both X1 and X2 to initialize Γ but this makes the initialization system # overdetermined because parameter initialization isn't in yet -@test_warn "overdetermined" oprob=ODEProblem(osys_m, u0, tspan, ps) +@test_warn "overdetermined" oprob=ODEProblem(osys_m, [u0; ps], tspan) # Make sure it doesn't error on array variables with unspecified size @parameters p::Vector{Real} q[1:3] @@ -93,7 +93,7 @@ eqs = [D(D(z)) ~ ones(2, 2)] B1 => A1, B2 => A2 ]) -prob = ODEProblem(sys, [], (0.0, 1.0), [A1 => 0.3]) +prob = ODEProblem(sys, [A1 => 0.3], (0.0, 1.0)) @test prob.ps[B1] == 0.3 @test prob.ps[B2] == 0.7 @@ -109,7 +109,7 @@ end @variables x(t) @parameters p @mtkcompile sys = System([D(x) ~ p], t; defaults = [x => t, p => 2t]) -prob = ODEProblem(sys, [], (1.0, 2.0), []) +prob = ODEProblem(sys, [], (1.0, 2.0)) @test prob[x] == 1.0 @test prob.ps[p] == 2.0 @@ -140,7 +140,7 @@ end @mtkcompile osys = System(eqs, t, sps, ps) u0map = [x => 1.0] pmap = [c1 => 5.0, c2 => 1.0, c3 => 1.2] - oprob = ODEProblem(osys, u0map, (0.0, 10.0), pmap) + oprob = ODEProblem(osys, [u0map; pmap], (0.0, 10.0)) end @testset "Cyclic dependency checking and substitution limits" begin @@ -160,12 +160,12 @@ end [D(x) ~ x * p, D(y) ~ y * q], t; guesses = [p => 1.0, q => 2.0]) # "unknowns" because they are initialization unknowns @test_warn ["Cycle", "unknowns", "p", "q"] try - ODEProblem(sys, [x => 1, y => 2], (0.0, 1.0), - [p => 2q, q => 3p]; warn_cyclic_dependency = true) + ODEProblem(sys, [x => 1, y => 2, p => 2q, q => 3p], + (0.0, 1.0); warn_cyclic_dependency = true) catch end @test_throws ModelingToolkit.MissingGuessError ODEProblem( - sys, [x => 1, y => 2], (0.0, 1.0), [p => 2q, q => 3p]) + sys, [x => 1, y => 2, p => 2q, q => 3p], (0.0, 1.0)) end @testset "`add_fallbacks!` checks scalarized array parameters correctly" begin @@ -201,8 +201,8 @@ end @mtkcompile pend = System(eqs, t) @test_throws ModelingToolkit.MissingGuessError ODEProblem( - pend, [x => 1], (0, 1), [g => 1], guesses = [y => λ, λ => y + 1]) - ODEProblem(pend, [x => 1], (0, 1), [g => 1], guesses = [y => λ, λ => 0.5]) + pend, [x => 1, g => 1], (0, 1), guesses = [y => λ, λ => y + 1]) + ODEProblem(pend, [x => 1, g => 1], (0, 1), guesses = [y => λ, λ => 0.5]) # Throw multiple if multiple are missing @variables a(t) b(t) c(t) d(t) e(t) @@ -221,7 +221,7 @@ end @mtkcompile sys = System(D(x) ~ interp(t), t) - prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), [interp => spline]) + prob = ODEProblem(sys, [x => 0.0, interp => spline], (0.0, 1.0)) spline2 = LinearInterpolation(ts .^ 2, ts .^ 2) p_new = [interp => spline2] prob2 = remake(prob; p = p_new) @@ -241,7 +241,7 @@ end # Creates the `NonlinearProblem`. u0 = [X => [1.0, 2.0]] ps = [p => [4.0, 5.0]] - @test_nowarn NonlinearProblem(nlsys, u0, ps) + @test_nowarn NonlinearProblem(nlsys, [u0; ps]) end @testset "Issue#3553: Retain `Float32` initial values" begin @@ -251,7 +251,7 @@ end @mtkcompile osys = System(eqs, t) u0 = [X => 1.0f0] ps = [p => 1.0f0, d => 2.0f0] - oprob = ODEProblem(osys, u0, (0.0f0, 1.0f0), ps) + oprob = ODEProblem(osys, [u0; ps], (0.0f0, 1.0f0)) sol = solve(oprob) @test eltype(oprob.u0) == Float32 @test eltype(eltype(sol.u)) == Float32 @@ -268,7 +268,7 @@ end @test Initial(D(x[i])) in ps end @test p in ps - prob = ODEProblem(sys, [x => ones(2)], (0.0, 1.0), [p => 1.0]) + prob = ODEProblem(sys, [x => ones(2), p => 1.0], (0.0, 1.0)) @test prob.p isa Vector{Float64} @test length(prob.p) == 5 end @@ -295,7 +295,7 @@ end ρ => 10.0] tspan = (0.0, 100.0) - prob = ODEProblem(sys, u0, tspan, p, jac = true, guesses = [w2 => -1.0], + prob = ODEProblem(sys, [u0; p], tspan, jac = true, guesses = [w2 => -1.0], warn_initialize_determined = false) @test prob[w2] ≈ -1.0 @test prob.ps[β] ≈ 8 / 3 @@ -321,7 +321,7 @@ end β => 8.0f0 / 3] tspan = (0.0f0, 100.0f0) - prob = ODEProblem(sys, u0, tspan, p) + prob = ODEProblem(sys, [u0; p], tspan) @test prob.p.tunable isa SVector @test prob.p.initials isa SVector end diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 8c1634e4e4..244e4ea4d3 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -13,7 +13,7 @@ eqs = [D(D(x)) ~ λ * x x^2 + y^2 ~ 1] @mtkcompile pend = System(eqs, t) -initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [], [g => 1]; +initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [g => 1]; guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2]) conditions = getfield.(equations(initprob.f.sys), :rhs) @@ -23,11 +23,11 @@ sol = solve(initprob) @test maximum(abs.(sol[conditions])) < 1e-14 @test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem( - pend, 0.0, [], [g => 1]; + pend, 0.0, [g => 1]; guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2], fully_determined = true) -initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0], [g => 1]; +initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0, g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) @test initprob isa NonlinearLeastSquaresProblem sol = solve(initprob) @@ -36,17 +36,17 @@ sol = solve(initprob) @test maximum(abs.(sol[conditions])) < 1e-14 initprob = ModelingToolkit.InitializationProblem( - pend, 0.0, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) + pend, 0.0, [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) @test initprob isa NonlinearLeastSquaresProblem sol = solve(initprob) @test !SciMLBase.successful_retcode(sol) || sol.retcode == SciMLBase.ReturnCode.StalledSuccess @test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem( - pend, 0.0, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend), + pend, 0.0, [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend), fully_determined = true) -prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], +prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5), guesses = ModelingToolkit.missing_variable_defaults(pend)) prob.f.initializeprob isa NonlinearProblem sol = solve(prob.f.initializeprob) @@ -54,7 +54,7 @@ sol = solve(prob.f.initializeprob) sol = solve(prob, Rodas5P()) @test maximum(abs.(sol[conditions][1])) < 1e-14 -prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], +prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5), guesses = ModelingToolkit.missing_variable_defaults(pend)) prob.f.initializeprob isa NonlinearLeastSquaresProblem sol = solve(prob.f.initializeprob) @@ -63,7 +63,7 @@ sol = solve(prob, Rodas5P()) @test maximum(abs.(sol[conditions][1])) < 1e-14 @test_throws ModelingToolkit.ExtraVariablesSystemException ODEProblem( - pend, [x => 1], (0.0, 1.5), [g => 1], + pend, [x => 1, g => 1], (0.0, 1.5), guesses = ModelingToolkit.missing_variable_defaults(pend), fully_determined = true) @@ -360,7 +360,7 @@ p = [σ => 28.0, tspan = (0.0, 100.0) @test_throws ModelingToolkit.IncompleteInitializationError prob=ODEProblem( - sys, u0, tspan, p, jac = true) + sys, [u0; p], tspan, jac = true) # DAE Initialization on ODE with nonlinear system for initial conditions # https://github.com/SciML/ModelingToolkit.jl/issues/2508 @@ -428,7 +428,7 @@ sol = solve(prob, Tsit5()) β => 8 / 3] tspan = (0.0, 0.2) - prob_mtk = ODEProblem(sys, u0, tspan, p) + prob_mtk = ODEProblem(sys, [u0; p], tspan) sol = solve(prob_mtk, Tsit5()) @test sol[x * (ρ - z) - y][1] == 0.0 @@ -474,10 +474,10 @@ eqs = [D(D(x)) ~ λ * x x^2 + y^2 ~ 1] @mtkcompile pend = System(eqs, t) -prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], +prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5), guesses = [λ => 0, y => 1], initialization_eqs = [y ~ 1]) -unsimp = generate_initializesystem(pend; u0map = [x => 1], initialization_eqs = [y ~ 1]) +unsimp = generate_initializesystem(pend; op = [x => 1], initialization_eqs = [y ~ 1]) sys = mtkcompile(unsimp; fully_determined = false) @test length(equations(sys)) in (3, 4) # could be either depending on tearing @@ -495,8 +495,8 @@ sys = extend(sysx, sysy) @variables x(t) y(t) @named sys = System([x^2 + y^2 ~ 25, D(x) ~ 1], t) ssys = mtkcompile(sys) - @test_throws ModelingToolkit.MissingVariablesError ODEProblem( - ssys, [x => 3], (0, 1), []) # y should have a guess + @test_throws ModelingToolkit.MissingGuessError ODEProblem( + ssys, [x => 3], (0, 1)) # y should have a guess end # https://github.com/SciML/ModelingToolkit.jl/issues/3025 @@ -506,7 +506,7 @@ end # system 1 should solve to x = 1 ics1 = [x => 1] sys1 = System([D(x) ~ 0], t; defaults = ics1, name = :sys1) |> mtkcompile - prob1 = ODEProblem(sys1, [], (0.0, 1.0), []) + prob1 = ODEProblem(sys1, [], (0.0, 1.0)) sol1 = solve(prob1, Tsit5()) @test all(sol1[x] .== 1) @@ -516,7 +516,7 @@ end System([D(y) ~ 0], t; initialization_eqs = [y ~ 2], name = :sys2) ) |> mtkcompile ics2 = unknowns(sys1) .=> 2 # should be equivalent to "ics2 = [x => 2]" - prob2 = ODEProblem(sys2, ics2, (0.0, 1.0), []; fully_determined = true) + prob2 = ODEProblem(sys2, ics2, (0.0, 1.0); fully_determined = true) sol2 = solve(prob2, Tsit5()) @test all(sol2[x] .== 2) && all(sol2[y] .== 2) end @@ -527,12 +527,12 @@ end sys = System( [D(D(x)) ~ 0], t; initialization_eqs = [x ~ 0, D(x) ~ 1], name = :sys) |> mtkcompile - @test_nowarn ODEProblem(sys, [], (0.0, 1.0), []) + @test_nowarn ODEProblem(sys, [], (0.0, 1.0)) sys = System( [D(D(x)) ~ 0], t; initialization_eqs = [x ~ 0, D(D(x)) ~ 0], name = :sys) |> mtkcompile - @test_nowarn ODEProblem(sys, [D(x) => 1.0], (0.0, 1.0), []) + @test_nowarn ODEProblem(sys, [D(x) => 1.0], (0.0, 1.0)) end # https://github.com/SciML/ModelingToolkit.jl/issues/3049 @@ -543,7 +543,7 @@ end [D(D(x)) ~ 0], t; initialization_eqs = [D(x)^2 ~ 1, x ~ 0], guesses = [D(x) => sign], name = :sys ) |> mtkcompile - prob = ODEProblem(sys, [], (0.0, 1.0), []) + prob = ODEProblem(sys, [], (0.0, 1.0)) sol = solve(prob, Tsit5()) @test sol(1.0, idxs = sys.x) ≈ sign # system with D(x(0)) = ±1 should solve to x(1) = ±1 end @@ -566,10 +566,10 @@ u0_2nd_order_2 = [Y => 2.0, D(Y) => 0.5] tspan = (0.0, 10.0) ps = [ω => 0.5, k1 => 2.0, k2 => 3.0] -oprob_1st_order_1 = ODEProblem(sys_1st_order, u0_1st_order_1, tspan, ps) -oprob_1st_order_2 = ODEProblem(sys_1st_order, u0_1st_order_2, tspan, ps) -oprob_2nd_order_1 = ODEProblem(sys_2nd_order, u0_2nd_order_1, tspan, ps) # gives sys_2nd_order -oprob_2nd_order_2 = ODEProblem(sys_2nd_order, u0_2nd_order_2, tspan, ps) +oprob_1st_order_1 = ODEProblem(sys_1st_order, [u0_1st_order_1; ps], tspan) +oprob_1st_order_2 = ODEProblem(sys_1st_order, [u0_1st_order_2; ps], tspan) +oprob_2nd_order_1 = ODEProblem(sys_2nd_order, [u0_2nd_order_1; ps], tspan) # gives sys_2nd_order +oprob_2nd_order_2 = ODEProblem(sys_2nd_order, [u0_2nd_order_2; ps], tspan) @test solve(oprob_1st_order_1, Rosenbrock23()).retcode == SciMLBase.ReturnCode.InitialFailure @@ -584,7 +584,7 @@ sol = solve(oprob_2nd_order_2, Rosenbrock23()) # retcode: Success @variables x(t)[1:5] y(t)[1:5] @named sys = System([D(x) ~ x, D(y) ~ y], t; initialization_eqs = [y ~ -x]) sys = mtkcompile(sys) - prob = ODEProblem(sys, [sys.x => ones(5)], (0.0, 1.0), []) + prob = ODEProblem(sys, [sys.x => ones(5)], (0.0, 1.0)) sol = solve(prob, Tsit5(), reltol = 1e-4) @test all(sol(1.0, idxs = sys.x) .≈ +exp(1)) && all(sol(1.0, idxs = sys.y) .≈ -exp(1)) end @@ -638,7 +638,7 @@ end @mtkcompile sys = System( [D(x) ~ x * q + rhss[1], D(y) ~ y * p + rhss[2]], t; defaults = [p => missing], guesses = [p => 1.0]) pmap[p] = 2q - prob = Problem(sys, u0map, (0.0, 1.0), pmap; u0_constructor, p_constructor) + prob = Problem(sys, merge(u0map, pmap), (0.0, 1.0); u0_constructor, p_constructor) test_parameter(prob, p, 2.0) prob2 = remake(prob; u0 = u0map, p = pmap) prob2 = remake(prob2; p = setp_oop(prob2, p)(prob2, 0.0)) @@ -657,7 +657,7 @@ end @mtkcompile sys = System( [D(x) ~ x * q + rhss[1], D(y) ~ y * p + rhss[2]], t; defaults = [p => 2q], guesses = [p => 1.0]) pmap[p] = missing - prob = Problem(sys, u0map, (0.0, 1.0), pmap; u0_constructor, p_constructor) + prob = Problem(sys, merge(u0map, pmap), (0.0, 1.0); u0_constructor, p_constructor) test_parameter(prob, p, 2.0) test_initializesystem(prob, p, p ~ 2q) prob2 = remake(prob; u0 = u0map, p = pmap) @@ -666,7 +666,7 @@ end # `missing` to Problem, provided guess @mtkcompile sys = System( [D(x) ~ x + rhss[1], p ~ x + y + rhss[2]], t; guesses = [p => 0.0]) - prob = Problem(sys, u0map, (0.0, 1.0), pmap; u0_constructor, p_constructor) + prob = Problem(sys, merge(u0map, pmap), (0.0, 1.0); u0_constructor, p_constructor) test_parameter(prob, p, 2.0) test_initializesystem(prob, p, p ~ x + y) prob2 = remake(prob; u0 = u0map, p = pmap) @@ -677,7 +677,7 @@ end @mtkcompile sys = System( [D(x) ~ x * q + rhss[1], D(y) ~ y * p + rhss[2]], t; defaults = [p => 2q], guesses = [p => 0.0]) delete!(pmap, p) - prob = Problem(sys, u0map, (0.0, 1.0), pmap; u0_constructor, p_constructor) + prob = Problem(sys, merge(u0map, pmap), (0.0, 1.0); u0_constructor, p_constructor) test_parameter(prob, p, 2.0) test_initializesystem(prob, p, p ~ 2q) prob2 = remake(prob; u0 = u0map, p = pmap) @@ -688,14 +688,14 @@ end @mtkcompile sys = System( [D(x) ~ q * x + rhss[1], D(y) ~ y * p + rhss[2]], t; defaults = [p => 2q], guesses = [p => 1.0]) _pmap = merge(pmap, Dict(p => q)) - prob = Problem(sys, u0map, (0.0, 1.0), _pmap; u0_constructor, p_constructor) + prob = Problem(sys, merge(u0map, _pmap), (0.0, 1.0); u0_constructor, p_constructor) test_parameter(prob, p, _pmap[q]) test_initializesystem(prob, p, p ~ q) # Problem dependent value with guess, no `missing` @mtkcompile sys = System( [D(x) ~ y * q + p + rhss[1], D(y) ~ x * p + q + rhss[2]], t; guesses = [p => 0.0]) _pmap = merge(pmap, Dict(p => 3q)) - prob = Problem(sys, u0map, (0.0, 1.0), _pmap; u0_constructor, p_constructor) + prob = Problem(sys, merge(u0map, _pmap), (0.0, 1.0); u0_constructor, p_constructor) test_parameter(prob, p, 3pmap[q]) # Should not be solved for: @@ -703,7 +703,7 @@ end @mtkcompile sys = System( [D(x) ~ q * x + rhss[1], D(y) ~ y * p + rhss[2]], t; defaults = [p => 2q], guesses = [p => 1.0]) _pmap = merge(pmap, Dict(p => 1.0)) - prob = Problem(sys, u0map, (0.0, 1.0), _pmap; u0_constructor, p_constructor) + prob = Problem(sys, merge(u0map, _pmap), (0.0, 1.0); u0_constructor, p_constructor) @test prob.ps[p] ≈ 1.0 initsys = prob.f.initialization_data.initializeprob.f.sys @test is_parameter(initsys, p) @@ -712,7 +712,8 @@ end @parameters r::Int s::Int @mtkcompile sys = System( [D(x) ~ s * x + rhss[1], D(y) ~ y * r + rhss[2]], t; defaults = [s => 2r], guesses = [s => 1.0]) - prob = Problem(sys, u0map, (0.0, 1.0), [r => 1]; u0_constructor, p_constructor) + prob = Problem( + sys, merge(u0map, Dict(r => 1)), (0.0, 1.0); u0_constructor, p_constructor) @test prob.ps[r] == 1 @test prob.ps[s] == 2 initsys = prob.f.initialization_data.initializeprob.f.sys @@ -725,8 +726,8 @@ end sys, [x => 1.0, y => 1.0], (0.0, 1.0)) # Unsatisfiable initialization - prob = Problem(sys, [x => 1.0, y => 1.0], (0.0, 1.0), - [p => 2.0]; initialization_eqs = [x^2 + y^2 ~ 3], u0_constructor, p_constructor) + prob = Problem(sys, [x => 1.0, y => 1.0, p => 2.0], (0.0, 1.0); + initialization_eqs = [x^2 + y^2 ~ 3], u0_constructor, p_constructor) @test prob.f.initialization_data !== nothing @test solve(prob, alg).retcode == ReturnCode.InitialFailure cache = init(prob, alg) @@ -737,7 +738,7 @@ end @variables x(t) y(t) s(t) @parameters x0 y0 @mtkcompile sys = System([x ~ x0, y ~ y0, s ~ x + y], t; guesses = [y0 => 0.0]) - prob = ODEProblem(sys, [s => 1.0], (0.0, 1.0), [x0 => 0.3, y0 => missing]) + prob = ODEProblem(sys, [s => 1.0, x0 => 0.3, y0 => missing], (0.0, 1.0)) # trivial initialization run immediately @test prob.ps[y0] ≈ 0.7 @test init(prob, Tsit5()).ps[y0] ≈ 0.7 @@ -763,7 +764,7 @@ end t; systems = [fixed, spring, mass, gravity, constant, damper], guesses = [spring.s_rel0 => 1.0]) - prob = ODEProblem(sys, [], (0.0, 1.0), [spring.s_rel0 => missing]) + prob = ODEProblem(sys, [spring.s_rel0 => missing], (0.0, 1.0)) # trivial initialization run immediately @test prob.ps[spring.s_rel0] ≈ -3.905 @test init(prob, Tsit5()).ps[spring.s_rel0] ≈ -3.905 @@ -891,7 +892,7 @@ end ] @mtkcompile sys = System( [D(x) ~ x + rhss[1], p ~ x + y + rhss[2]], t; guesses = [x => 0.0, p => 0.0]) - prob = Problem(sys, [y => 1.0], (0.0, 1.0), [p => 3.0]) + prob = Problem(sys, [y => 1.0, p => 3.0], (0.0, 1.0)) @test prob.f.initialization_data.initializeprob.ps[p] ≈ 3.0 @test init(prob, alg)[x] ≈ 2.0 prob.ps[p] = 2.0 @@ -921,7 +922,7 @@ end @mtkcompile sys = System( [D(x) ~ 2x + r + rhss], t; parameter_dependencies = [r ~ p + 2q, q ~ p + 3], guesses = [p => 1.0]) - prob = Problem(sys, [x => 1.0], (0.0, 1.0), [p => missing]) + prob = Problem(sys, [x => 1.0, p => missing], (0.0, 1.0)) @test length(equations(ModelingToolkit.get_parent(prob.f.initialization_data.initializeprob.f.sys))) == 4 integ = init(prob, alg) @@ -978,7 +979,7 @@ end @mtkcompile sys = System( [D(x) ~ x * p + y^2 * q + rhss; alge_eqs], t; guesses = [x => 0.0, y => 0.0, z => 0.0, p => 0.0, q => 0.0]) - prob = Problem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0, q => missing]) + prob = Problem(sys, [x => 1.0, p => 1.0, q => missing], (0.0, 1.0)) @test is_variable(prob.f.initialization_data.initializeprob, q) ps = prob.p newps = SciMLStructures.replace(Tunable(), ps, ForwardDiff.Dual.(ps.tunable)) @@ -1034,12 +1035,12 @@ end @mtkcompile sys = System( [D(x) ~ x + p * y^2 + rhss; alge_eqs], t; guesses = [ y => 1.0, p => 1.0]) - prob = Problem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0]) + prob = Problem(sys, [x => 1.0, p => 1.0], (0.0, 1.0)) @test is_variable(prob.f.initialization_data.initializeprob, y) prob2 = @test_nowarn remake(prob; p = [p => 3.0]) # ensure no over/under-determined warning @test is_variable(prob.f.initialization_data.initializeprob, y) - prob = Problem(sys, [y => 1.0, x => 2.0], (0.0, 1.0), [p => missing]) + prob = Problem(sys, [y => 1.0, x => 2.0, p => missing], (0.0, 1.0)) @test is_variable(prob.f.initialization_data.initializeprob, p) prob2 = @test_nowarn remake(prob; u0 = [y => 0.5]) @test is_variable(prob.f.initialization_data.initializeprob, p) @@ -1062,7 +1063,7 @@ end @parameters foo(::Real, ::Real) p @mtkcompile sys = System([D(x) ~ t, 0 ~ foo(x, y)], t; parameter_dependencies = [foo ~ Multiplier(p, 2p)], guesses = [y => -1.0]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0]) + prob = ODEProblem(sys, [x => 1.0, p => 1.0], (0.0, 1.0)) integ = init(prob, Rosenbrock23()) @test integ[y] ≈ -0.5 end @@ -1101,7 +1102,7 @@ end x^2 + y^2 ~ L] @mtkcompile pend = System(eqs, t) - prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], + prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5), guesses = ModelingToolkit.missing_variable_defaults(pend)) sol = solve(prob, Rodas5P()) @test SciMLBase.successful_retcode(sol) @@ -1185,18 +1186,19 @@ end x^2 + y^2 ~ 1] @mtkcompile pend = System(eqs, t) @test_warn ["structurally singular", "initialization", "Guess", "heuristic"] ODEProblem( - pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = [λ => 1]) + pend, [x => 1, y => 0, g => 1], (0.0, 1.5), guesses = [λ => 1]) end @testset "DAEProblem initialization" begin @variables x(t) [guess = 1.0] y(t) [guess = 1.0] @parameters p=missing [guess = 1.0] q=missing [guess = 1.0] @mtkcompile sys = System( - [D(x) ~ p * y + q * t, x^3 + y^3 ~ 5], t; initialization_eqs = [p^2 + q^3 ~ 3]) + [D(x) ~ p * y + q, x^3 + y^3 ~ 5], t; initialization_eqs = [p^2 + q^3 ~ 3]) # FIXME: solve for du0 prob = DAEProblem( - sys, [D(x) => cbrt(4), D(y) => -1 / cbrt(4)], [x => 1.0], (0.0, 1.0), [p => 1.0]) + sys, [D(x) => cbrt(4) + cbrt(2), D(y) => -1 / cbrt(4), x => 1.0, p => 1.0], ( + 0.0, 1.0)) integ = init(prob, DImplicitEuler()) @test integ[x] ≈ 1.0 @@ -1210,7 +1212,7 @@ end @parameters p q=3x @mtkcompile sys = System([D(x) ~ x * p + q, x^3 + y^3 ~ 3], t) prob = ODEProblem( - sys, [], (0.0, 1.0), [p => 1.0]; guesses = [x => 1.0, y => 1.0, q => 1.0]) + sys, [p => 1.0], (0.0, 1.0); guesses = [x => 1.0, y => 1.0, q => 1.0]) @test prob[x] == 1.0 @test prob[y] == 2.0 @test prob.ps[p] == 1.0 @@ -1240,7 +1242,7 @@ end @parameters p [guess = 1.0] q [guess = 1.0] @mtkcompile sys = System( [D(x) ~ p * x + q * y, y ~ 2x], t; parameter_dependencies = [q ~ 2p]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => 1.0]) + prob = ODEProblem(sys, [x => 1.0, p => 1.0], (0.0, 1.0)) test_dummy_initialization_equation(prob, x) prob2 = remake(prob; u0 = [x => 2.0]) @test prob2[x] == 2.0 @@ -1261,7 +1263,7 @@ end @parameters p [guess = 1.0] q [guess = 1.0] @mtkcompile sys = System( [D(x) ~ p * x + q * y, y ~ 2x], t; parameter_dependencies = [q ~ 2p]) - prob = ODEProblem(sys, [:x => 1.0], (0.0, 1.0), [p => 1.0]) + prob = ODEProblem(sys, [:x => 1.0, p => 1.0], (0.0, 1.0)) test_dummy_initialization_equation(prob, x) prob2 = remake(prob; u0 = [:x => 2.0]) test_dummy_initialization_equation(prob2, x) @@ -1276,7 +1278,7 @@ end @named sys = System([D(x) ~ 0, D(y) ~ x + a], t; initialization_eqs = [y ~ a]) ssys = mtkcompile(sys) - prob = ODEProblem(ssys, [], (0, 1), []) + prob = ODEProblem(ssys, [], (0, 1)) @test SciMLBase.successful_retcode(solve(prob)) @@ -1300,7 +1302,7 @@ end u0_vals = [X => 4, Y => 5.0] tspan = (0.0, 10.0) p_vals = [p => 1.0, d => 0.1] - oprob = ODEProblem(osys, u0_vals, tspan, p_vals) + oprob = ODEProblem(osys, [u0_vals; p_vals], tspan) integ = init(oprob) @test integ[X] ≈ 4.0 @test integ[Y] ≈ 5.0 @@ -1326,7 +1328,7 @@ end u0s = [I => 1, R => 0] ps = [S0 => 999, β => 0.01, γ => 0.001] - jprob = JumpProblem(js, u0s, (0.0, 10.0), ps) + jprob = JumpProblem(js, [u0s; ps], (0.0, 10.0)) sol = solve(jprob, SSAStepper()) @test sol[S, 1] ≈ 999 @test SciMLBase.successful_retcode(sol) @@ -1339,7 +1341,7 @@ end D(x) ~ p[1] + p[2] + q, t; defaults = [p[1] => q, p[2] => 2q], guesses = [p[1] => q, p[2] => 2q]) @test ModelingToolkit.is_parameter_solvable(p, Dict(), defaults(sys), guesses(sys)) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [q => 2.0]) + prob = ODEProblem(sys, [x => 1.0, q => 2.0], (0.0, 1.0)) initsys = prob.f.initialization_data.initializeprob.f.sys @test length(ModelingToolkit.observed(initsys)) == 4 sol = solve(prob, Tsit5()) @@ -1351,7 +1353,7 @@ end @parameters p[1:2, 1:2] @mtkcompile sys = System( [D(x) ~ x, D(y) ~ p * y], t; initialization_eqs = [x^2 + y[1]^2 + y[2]^2 ~ 4]) - prob = ODEProblem(sys, [x => 1.0, y[1] => 1], (0.0, 1.0), [p => 2ones(2, 2)]) + prob = ODEProblem(sys, [x => 1.0, y[1] => 1, p => 2ones(2, 2)], (0.0, 1.0)) integ = init(prob, Tsit5()) @test integ[x] ≈ 1.0 @test integ[y] ≈ [1.0, sqrt(2.0)] @@ -1376,7 +1378,7 @@ end [y ~ 0.5] => (; f = stop!) ]) sys = mtkcompile(sys) - prob0 = ODEProblem(sys, [x => NaN], (0.0, 1.0), []) + prob0 = ODEProblem(sys, [x => NaN], (0.0, 1.0)) # final_x(x0) is equivalent to x0 + 0.5 function final_x(x0) @@ -1399,7 +1401,7 @@ end @parameters p @variables x(t) @mtkcompile sys = System(x ~ p * t, t) - prob = @test_nowarn ODEProblem(sys, [], (0.0, 1.0), [p => 1.0]) + prob = @test_nowarn ODEProblem(sys, [p => 1.0], (0.0, 1.0)) @test_nowarn remake(prob, p = [p => 1.0]) @test_nowarn remake(prob, p = [p => ForwardDiff.Dual(1.0)]) end @@ -1416,7 +1418,7 @@ end u0 = [X1 => 1.0, X2 => 2.0] ps = [k1 => 0.1, k2 => 0.2] - oprob1 = ODEProblem(osys, u0, 1.0, ps) + oprob1 = ODEProblem(osys, [u0; ps], 1.0) oprob2 = remake(oprob1, u0 = [X1 => 10.0]) integ1 = init(oprob1) @test integ1[X1] ≈ 1.0 @@ -1443,7 +1445,7 @@ end @mtkcompile sys = System([x^2 + y^2 ~ p1, (x - 1)^2 + (y - 1)^2 ~ p2]; parameter_dependencies = [p2 ~ 2p1], guesses = [p1 => 0.0], defaults = [p1 => missing]) - prob = Problem(sys, [x => 1.0, y => 1.0], [p2 => 6.0]) + prob = Problem(sys, [x => 1.0, y => 1.0, p2 => 6.0]) @test prob.ps[p1] ≈ 3.0 end end @@ -1464,7 +1466,7 @@ end @testset "throws if initialization_eqs contain unknowns" begin u0 = [X1 => 1.0, X2 => 2.0] ps = [k1 => 0.1, k2 => 0.2] - @test_throws ArgumentError Problem(nlsys, u0, ps) + @test_throws ArgumentError Problem(nlsys, [u0; ps]) end eqs = [0 ~ k1 * (Γ[1] - X1) - k2 * X1 @@ -1477,22 +1479,21 @@ end @testset "solves initialization" begin u0 = [X1 => 1.0, X2 => 2.0] ps = [k1 => 0.1, k2 => 0.2] - prob = Problem(nlsys, u0, ps) + prob = Problem(nlsys, [u0; ps]) @test state_values(prob.f.initialization_data.initializeprob) === nothing @test prob.ps[Γ[1]] ≈ 3.0 end @testset "respects explicitly provided value" begin - u0 = [] ps = [k1 => 0.1, k2 => 0.2, Γ => [5.0]] - prob = Problem(nlsys, u0, ps) + prob = Problem(nlsys, ps) @test prob.ps[Γ[1]] ≈ 5.0 end @testset "fails initialization if inconsistent explicit value" begin u0 = [X1 => 1.0, X2 => 2.0] ps = [k1 => 0.1, k2 => 0.2, Γ => [5.0]] - prob = Problem(nlsys, u0, ps) + prob = Problem(nlsys, [u0; ps]) sol = solve(prob) @test sol.retcode == SciMLBase.ReturnCode.InitialFailure end @@ -1500,7 +1501,7 @@ end @testset "Ignores initial equation if given insufficient u0" begin u0 = [X2 => 2.0] ps = [k1 => 0.1, k2 => 0.2, Γ => [5.0]] - prob = Problem(nlsys, u0, ps) + prob = Problem(nlsys, [u0; ps]) sol = solve(prob) @test SciMLBase.successful_retcode(sol) @test sol.ps[Γ[1]] ≈ 5.0 @@ -1511,7 +1512,7 @@ end @variables x(t) y(t) @parameters c1 c2 @mtkcompile sys = System([D(x) ~ -c1 * x + c2 * y, D(y) ~ c1 * x - c2 * y], t) - prob1 = ODEProblem(sys, [1.0, 2.0], (0.0, 1.0), [c1 => 1.0, c2 => 2.0]) + prob1 = ODEProblem(sys, [x => 1.0, y => 2.0, c1 => 1.0, c2 => 2.0], (0.0, 1.0)) prob2 = remake(prob1, u0 = [2.0, 3.0]) prob3 = remake(prob1, u0 = [2.0, 3.0], p = [c1 => 2.0]) integ1 = init(prob1, Tsit5()) @@ -1542,8 +1543,8 @@ end @mtkcompile pend = System(eqs, t) prob = ODEProblem( - pend, [x => (√2 / 2), D(x) => 0.0], (0.0, 1.5), - [g => 1], guesses = [λ => 1, y => √2 / 2]) + pend, [x => (√2 / 2), D(x) => 0.0, g => 1], (0.0, 1.5), + guesses = [λ => 1, y => √2 / 2]) sol = solve(prob) @testset "Guesses of initialization problem copied to algebraic variables" begin @@ -1555,8 +1556,8 @@ end @testset "Initial values for algebraic variables are retained" begin prob2 = ODEProblem( - pend, [x => (√2 / 2), D(y) => 0.0], (0.0, 1.5), - [g => 1], guesses = [λ => 1, y => √2 / 2]) + pend, [x => (√2 / 2), D(y) => 0.0, g => 1], (0.0, 1.5), + guesses = [λ => 1, y => √2 / 2]) sol = solve(prob) @test SciMLBase.successful_retcode(sol) prob3 = DiffEqBase.get_updated_symbolic_problem( @@ -1607,7 +1608,7 @@ end tspan = (0.0, 100.0) getter = getsym(sys, Initial.(unknowns(sys))) - prob = ODEProblem(sys, u0, tspan, p; guesses = [w2 => 3.0]) + prob = ODEProblem(sys, [u0; p], tspan; guesses = [w2 => 3.0]) new_u0, new_p, _ = SciMLBase.get_initial_values( prob, prob, prob.f, SciMLBase.OverrideInit(), Val(true); nlsolve_alg = NewtonRaphson(), abstol = 1e-6, reltol = 1e-3) @@ -1626,8 +1627,8 @@ end D(D(y)) ~ λ * y - g x^2 + y^2 ~ 1] @mtkcompile pend=System(eqs, t) split=false - prob = ODEProblem(pend, [x => 1.0, D(x) => 0.0], (0.0, 1.0), - [g => 1.0]; guesses = [y => 1.0, λ => 1.0]) + prob = ODEProblem( + pend, [x => 1.0, D(x) => 0.0, g => 1.0], (0.0, 1.0); guesses = [y => 1.0, λ => 1.0]) @test !ModelingToolkit.is_split(prob.f.initialization_data.initializeprob.f.sys) end @@ -1638,8 +1639,8 @@ end D(D(y)) ~ λ * y - g x^2 + y^2 ~ 1] @mtkcompile pend = System(eqs, t) - prob = ODEProblem(pend, SA[x => 1.0, D(x) => 0.0], (0.0, 1.0), - SA[g => 1.0]; guesses = [y => 1.0, λ => 1.0]) + prob = ODEProblem(pend, SA[x => 1.0, D(x) => 0.0, g => 1.0], + (0.0, 1.0); guesses = [y => 1.0, λ => 1.0]) @test !SciMLBase.isinplace(prob) @test !SciMLBase.isinplace(prob.f.initialization_data.initializeprob) end diff --git a/test/jacobiansparsity.jl b/test/jacobiansparsity.jl index 56e8edc183..7173ab9e1a 100644 --- a/test/jacobiansparsity.jl +++ b/test/jacobiansparsity.jl @@ -95,7 +95,7 @@ prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = true) u0 = [x => 1, y => 0] prob = ODEProblem( - pend, u0, (0, 11.5), [g => 1], guesses = [λ => 1], sparse = true, jac = true) + pend, [u0; [g => 1]], (0, 11.5), guesses = [λ => 1], sparse = true, jac = true) jac, jac! = generate_jacobian(pend; expression = Val{false}, sparse = true) jac_prototype = ModelingToolkit.jacobian_sparsity(pend) W_prototype = ModelingToolkit.W_sparsity(pend) @@ -126,7 +126,7 @@ end D(D(y)) ~ λ * y - g x^2 + y^2 ~ 1] @mtkcompile pend = System(eqs, t) - prob = ODEProblem(pend, [x => 0.0, D(x) => 1.0], (0.0, 1.0), [g => 1.0]; + prob = ODEProblem(pend, [x => 0.0, D(x) => 1.0, g => 1.0], (0.0, 1.0); guesses = [y => 1.0, λ => 1.0], jac = true, sparse = true) J = deepcopy(prob.f.jac_prototype) prob.f.jac(J, prob.u0, prob.p, 1.0) diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl index 95cc29e7b0..ca3d56359f 100644 --- a/test/jumpsystem.jl +++ b/test/jumpsystem.jl @@ -69,7 +69,7 @@ u₀ = [999, 1, 0]; tspan = (0.0, 250.0); u₀map = [S => 999, I => 1, R => 0] parammap = [β => 0.1 / 1000, γ => 0.01] -jprob = JumpProblem(js2, u₀map, tspan, parammap; aggregator = Direct(), +jprob = JumpProblem(js2, [u₀map; parammap], tspan; aggregator = Direct(), save_positions = (false, false), rng) p = parameter_values(jprob) @test jprob.prob isa DiscreteProblem @@ -85,7 +85,7 @@ end m = getmean(jprob, Nsims) # test auto-alg selection works -jprobb = JumpProblem(js2, u₀map, tspan, parammap; save_positions = (false, false), rng) +jprobb = JumpProblem(js2, [u₀map; parammap], tspan; save_positions = (false, false), rng) mb = getmean(jprobb, Nsims; use_stepper = false) @test abs(m - mb) / m < 0.01 @@ -93,14 +93,14 @@ mb = getmean(jprobb, Nsims; use_stepper = false) obs = [S2 ~ 2 * S] @named js2b = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ, h], observed = obs) js2b = complete(js2b) -jprob = JumpProblem(js2b, u₀map, tspan, parammap; aggregator = Direct(), +jprob = JumpProblem(js2b, [u₀map; parammap], tspan; aggregator = Direct(), save_positions = (false, false), rng) @test jprob.prob isa DiscreteProblem sol = solve(jprob, SSAStepper(); saveat = tspan[2] / 10) @test all(2 .* sol[S] .== sol[S2]) # test save_positions is working -jprob = JumpProblem(js2, u₀map, tspan, parammap; aggregator = Direct(), +jprob = JumpProblem(js2, [u₀map; parammap], tspan; aggregator = Direct(), save_positions = (false, false), rng) sol = solve(jprob, SSAStepper(); saveat = 1.0) @test all((sol.t) .== collect(0.0:tspan[2])) @@ -149,7 +149,7 @@ maj1 = MassActionJump(2 * β / 2, [S => 1, I => 1], [S => -1, I => 1]) maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1]) @named js3 = JumpSystem([maj1, maj2], t, [S, I, R], [β, γ]) js3 = complete(js3) -jprob = JumpProblem(js3, u₀map, tspan, parammap; aggregator = Direct(), rng) +jprob = JumpProblem(js3, [u₀map; parammap], tspan; aggregator = Direct(), rng) @test jprob.prob isa DiscreteProblem m3 = getmean(jprob, Nsims) @test abs(m - m3) / m < 0.01 @@ -157,11 +157,11 @@ m3 = getmean(jprob, Nsims) # maj jump test with various dep graphs @named js3b = JumpSystem([maj1, maj2], t, [S, I, R], [β, γ]) js3b = complete(js3b) -jprobb = JumpProblem(js3b, u₀map, tspan, parammap; aggregator = NRM(), rng) +jprobb = JumpProblem(js3b, [u₀map; parammap], tspan; aggregator = NRM(), rng) @test jprobb.prob isa DiscreteProblem m4 = getmean(jprobb, Nsims) @test abs(m - m4) / m < 0.01 -jprobc = JumpProblem(js3b, u₀map, tspan, parammap; aggregator = RSSA(), rng) +jprobc = JumpProblem(js3b, [u₀map; parammap], tspan; aggregator = RSSA(), rng) @test jprobc.prob isa DiscreteProblem m4 = getmean(jprobc, Nsims) @test abs(m - m4) / m < 0.01 @@ -172,7 +172,7 @@ maj2 = MassActionJump(γ, [S => 1], [S => -1]) @named js4 = JumpSystem([maj1, maj2], t, [S], [β, γ]) js4 = complete(js4) jprob = JumpProblem( - js4, [S => 999], (0, 1000.0), [β => 100.0, γ => 0.01]; aggregator = Direct(), rng) + js4, [S => 999, β => 100.0, γ => 0.01], (0, 1000.0); aggregator = Direct(), rng) @test jprob.prob isa DiscreteProblem m4 = getmean(jprob, Nsims) @test abs(m4 - 2.0 / 0.01) * 0.01 / 2.0 < 0.01 @@ -183,7 +183,7 @@ maj2 = MassActionJump(γ, [S => 2], [S => -1]) @named js4 = JumpSystem([maj1, maj2], t, [S], [β, γ]) js4 = complete(js4) jprob = JumpProblem( - js4, [S => 999], (0, 1000.0), [β => 100.0, γ => 0.01]; aggregator = Direct(), rng) + js4, [S => 999, β => 100.0, γ => 0.01], (0, 1000.0); aggregator = Direct(), rng) @test jprob.prob isa DiscreteProblem sol = solve(jprob, SSAStepper()); @@ -208,7 +208,7 @@ end u₀ = [A => 100, B => 0] tspan = (0.0, 2000.0) jprob = JumpProblem( - js5, u₀, tspan, p; aggregator = Direct(), save_positions = (false, false), rng) + js5, [u₀; p], tspan; aggregator = Direct(), save_positions = (false, false), rng) @test jprob.prob isa DiscreteProblem @test all(jprob.massaction_jump.scaled_rates .== [1.0, 0.0]) @@ -274,9 +274,9 @@ u0 = [X => 10] tspan = (0.0, 1.0) ps = [k => 1.0] -@test_nowarn jp1 = JumpProblem(js1, u0, tspan, ps; aggregator = Direct()) +@test_nowarn jp1 = JumpProblem(js1, [u0; ps], tspan; aggregator = Direct()) @test_nowarn jp2 = JumpProblem(js2, u0, tspan; aggregator = Direct()) -@test_nowarn jp3 = JumpProblem(js3, u0, tspan, ps; aggregator = Direct()) +@test_nowarn jp3 = JumpProblem(js3, [u0; ps], tspan; aggregator = Direct()) @test_nowarn jp4 = JumpProblem(js4, u0, tspan; aggregator = Direct()) # Ensure `mtkcompile` (and `@mtkcompile`) works on JumpSystem (by doing nothing) @@ -302,7 +302,7 @@ j1 = ConstantRateJump(k, [X ~ Pre(X) - 1]) for (N, algtype) in zip(Nv, algtypes) @named jsys = JumpSystem([deepcopy(j1) for _ in 1:N], t, [X], [k]) jsys = complete(jsys) - jprob = JumpProblem(jsys, [X => 10], (0.0, 10.0), [k => 1]) + jprob = JumpProblem(jsys, [X => 10, k => 1], (0.0, 10.0)) @test jprob.aggregator isa algtype end end @@ -316,7 +316,7 @@ end vrj = VariableRateJump(k * (sin(t) + 1), [A ~ Pre(A) + 1, C ~ Pre(C) + 2]) js = complete(JumpSystem([vrj], t, [A, C], [k]; name = :js, observed = [B ~ C * A])) jprob = JumpProblem( - js, [A => 0, C => 0], (0.0, 10.0), [k => 1.0]; aggregator = Direct(), rng) + js, [A => 0, C => 0, k => 1], (0.0, 10.0); aggregator = Direct(), rng) @test jprob.prob isa ODEProblem sol = solve(jprob, Tsit5()) @@ -449,7 +449,7 @@ end k2val = 20.0 p = [k1 => k1val, k2 => k2val] tspan = (0.0, 10.0) - jprob = JumpProblem(jsys, u0, tspan, p; rng, save_positions = (false, false)) + jprob = JumpProblem(jsys, [u0; p], tspan; rng, save_positions = (false, false)) times = range(0.0, tspan[2], length = 100) Nsims = 4000 @@ -488,7 +488,7 @@ end u0map = [X => p.X₀, Y => p.Y₀] pmap = [α => p.α, β => p.β] tspan = (0.0, 20.0) - jprob = JumpProblem(jsys, u0map, tspan, pmap; rng, save_positions = (false, false)) + jprob = JumpProblem(jsys, [u0map; pmap], tspan; rng, save_positions = (false, false)) times = range(0.0, tspan[2], length = 100) Nsims = 4000 Xv = zeros(length(times)) @@ -526,7 +526,7 @@ end continuous_events = cevents) jsys = complete(jsys) tspan = (0.0, 200.0) - jprob = JumpProblem(jsys, u0map, tspan, pmap; rng, save_positions = (false, false)) + jprob = JumpProblem(jsys, [u0map; pmap], tspan; rng, save_positions = (false, false)) Xsamp = 0.0 Nsims = 4000 for n in 1:Nsims @@ -552,7 +552,7 @@ end # Works. @mtkcompile js = JumpSystem([j1, j2], t, [X], [p, d]) jprob = JumpProblem( - js, [X => 15], (0.0, 10.0), [p => 2.0, d => 0.5]; aggregator = Direct()) + js, [X => 15, p => 2.0, d => 0.5], (0.0, 10.0); aggregator = Direct(), u0_eltype = Int) sol = solve(jprob, SSAStepper()) @test eltype(sol[X]) === Int64 end @@ -566,7 +566,7 @@ end crj = ConstantRateJump(rate, affect) @named jsys = JumpSystem([crj, eq], t, [X], [a, b]) jsys = complete(jsys) - jprob = JumpProblem(jsys, [:X => 1.0], (0.0, 10.0), [:a => 1.0, :b => 0.5]) + jprob = JumpProblem(jsys, [:X => 1.0, :a => 1.0, :b => 0.5], (0.0, 10.0)) jprob2 = remake(jprob; u0 = [:X => 10.0]) @test jprob2[X] ≈ 10.0 end diff --git a/test/mass_matrix.jl b/test/mass_matrix.jl index 2e828cda64..82d1cf86a5 100644 --- a/test/mass_matrix.jl +++ b/test/mass_matrix.jl @@ -17,12 +17,10 @@ M = calculate_massmatrix(sys) 0 1 0 0 0 0] -prob_mm = ODEProblem(sys, [y => [1.0, 0.0, 0.0]], (0.0, 1e5), - [k => [0.04, 3e7, 1e4]]) +prob_mm = ODEProblem(sys, [y => [1.0, 0.0, 0.0], k => [0.04, 3e7, 1e4]], (0.0, 1e5)) @test prob_mm.f.mass_matrix isa Diagonal{Float64, Vector{Float64}} sol = solve(prob_mm, Rodas5(), reltol = 1e-8, abstol = 1e-8) -prob_mm = ODEProblem(sys, SA[y => [1.0, 0.0, 0.0]], (0.0, 1e5), - [k => [0.04, 3e7, 1e4]]) +prob_mm = ODEProblem(sys, SA[y => [1.0, 0.0, 0.0], k => [0.04, 3e7, 1e4]], (0.0, 1e5)) @test prob_mm.f.mass_matrix isa Diagonal{Float64, SVector{3, Float64}} function rober(du, u, p, t) @@ -57,8 +55,8 @@ eqs = [D(y[1]) ~ y[1], D(y[2]) ~ y[2], D(y[3]) ~ y[3]] @named sys = System(eqs, t, collect(y), [k]) @named sys = SDESystem(sys, [1, 1, 0]) sys = complete(sys) - prob = SDEProblem(sys, [y => [1.0, 0.0, 0.0]], (0.0, 1e5), [k => [0.04, 3e7, 1e4]]) + prob = SDEProblem(sys, [y => [1.0, 0.0, 0.0], k => [0.04, 3e7, 1e4]], (0.0, 1e5)) @test prob.f.mass_matrix isa Diagonal{Float64, Vector{Float64}} - prob = SDEProblem(sys, SA[y => [1.0, 0.0, 0.0]], (0.0, 1e5), [k => [0.04, 3e7, 1e4]]) + prob = SDEProblem(sys, SA[y => [1.0, 0.0, 0.0], k => [0.04, 3e7, 1e4]], (0.0, 1e5)) @test prob.f.mass_matrix isa Diagonal{Float64, SVector{3, Float64}} end diff --git a/test/modelingtoolkitize.jl b/test/modelingtoolkitize.jl index f6b4dc8a8d..3f2db92344 100644 --- a/test/modelingtoolkitize.jl +++ b/test/modelingtoolkitize.jl @@ -57,7 +57,8 @@ p = [1.0, 100.0] prob = OptimizationProblem(rosenbrock, x0, p) sys = complete(modelingtoolkitize(prob)) # symbolicitize me captain! -prob = OptimizationProblem(sys, x0, p, grad = true, hess = true) +prob = OptimizationProblem( + sys, [unknowns(sys) .=> x0; parameters(sys) .=> p], grad = true, hess = true) sol = solve(prob, NelderMead()) @test sol.objective < 1e-8 @@ -155,7 +156,7 @@ problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫) problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫) sys = complete(modelingtoolkitize(problem)) -fast_problem = ODEProblem(sys, ℬ, 𝒯, parameters(sys) .=> 𝒫) +fast_problem = ODEProblem(sys, [unknowns(sys) .=> ℬ; parameters(sys) .=> 𝒫], 𝒯) @time solution = solve(fast_problem, Tsit5(), saveat = 1:final_time) ## Issue #778 diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl index afd97a5984..2f52dad30b 100644 --- a/test/mtkparameters.jl +++ b/test/mtkparameters.jl @@ -107,7 +107,7 @@ end @test getp(sys, g)(newps) isa Vector{Float32} @testset "Type-stability of `remake_buffer`" begin - prob = ODEProblem(sys, [], (0.0, 1.0), ivs) + prob = ODEProblem(sys, ivs, (0.0, 1.0)) idxs = (a, c, d, e, f, g, h) vals = (1.0, 2.0, 3, ones(3), ones(Int, 3, 3), ones(2), "a") @@ -144,7 +144,7 @@ eq = D(X) ~ p[1] - p[2] * X u0 = [X => 1.0] ps = [p => [2.0, 0.1]] -p = MTKParameters(osys, ps, u0) +p = MTKParameters(osys, [ps; u0]) @test p.tunable == [2.0, 0.1] # Ensure partial update promotes the buffer @@ -166,8 +166,8 @@ u0 = [X => 1.0] tspan = (0.0, 100.0) ps = [p => 1.0] # Value for `d` is missing -@test_throws ModelingToolkit.MissingParametersError ODEProblem(sys, u0, tspan, ps) -@test_nowarn ODEProblem(sys, u0, tspan, [ps..., d => 1.0]) +@test_throws ModelingToolkit.MissingParametersError ODEProblem(sys, [u0; ps], tspan) +@test_nowarn ODEProblem(sys, [u0; ps; [d => 1.0]], tspan) # JET tests @@ -250,7 +250,7 @@ eqs = [D(x) ~ (α - β * y) * x D(y) ~ (δ * x - γ) * y] @mtkcompile odesys = System(eqs, t) odeprob = ODEProblem( - odesys, [x => 1.0, y => 1.0], (0.0, 10.0), [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0]) + odesys, [x => 1.0, y => 1.0, α => 1.5, β => 1.0, γ => 3.0, δ => 1.0], (0.0, 10.0)) tunables, _... = canonicalize(Tunable(), odeprob.p) @test tunables isa AbstractVector{Float64} @@ -328,7 +328,7 @@ end u0 = [V => [10.0, 20.0]] ps_vec = [k => [2.0, 3.0, 4.0, 5.0]] ps_scal = [k[1] => 1.0, k[2] => 2.0, k[3] => 3.0, k[4] => 4.0] - oprob_scal_scal = ODEProblem(osys_scal, u0, 1.0, ps_scal) + oprob_scal_scal = ODEProblem(osys_scal, [u0; ps_scal], 1.0) newoprob = remake(oprob_scal_scal; p = ps_vec, build_initializeprob = false) @test newoprob.ps[k] == [2.0, 3.0, 4.0, 5.0] end diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl index 8575fbc36a..bf05e317ce 100644 --- a/test/nonlinearsystem.jl +++ b/test/nonlinearsystem.jl @@ -80,16 +80,15 @@ sH = calculate_hessian(ns) @test getfield.(ModelingToolkit.hessian_sparsity(ns), :rowval) == getfield.(sparse.(sH), :rowval) -prob = NonlinearProblem(ns, ones(3), [σ => 1.0, ρ => 1.0, β => 1.0]) +prob = NonlinearProblem(ns, [x => 1.0, y => 1.0, z => 1.0, σ => 1.0, ρ => 1.0, β => 1.0]) @test prob.f.sys === ns sol = solve(prob, NewtonRaphson()) @test sol.u[1] ≈ sol.u[2] -prob = NonlinearProblem(ns, ones(3), [σ => 1.0, ρ => 1.0, β => 1.0], jac = true) +prob = NonlinearProblem( + ns, [x => 1.0, y => 1.0, z => 1.0, σ => 1.0, ρ => 1.0, β => 1.0], jac = true) @test_nowarn solve(prob, NewtonRaphson()) -@test_throws ArgumentError NonlinearProblem(ns, ones(4), [σ => 1.0, ρ => 1.0, β => 1.0]) - @variables u F s a eqs1 = [ 0 ~ σ * (y - x) * h + F, @@ -100,7 +99,7 @@ eqs1 = [ lorenz = name -> System(eqs1, [x, y, z, u, F], [σ, ρ, β, h], name = name) lorenz1 = lorenz(:lorenz1) -@test_throws ArgumentError NonlinearProblem(complete(lorenz1), zeros(5), zeros(3)) +@test_throws ArgumentError NonlinearProblem(complete(lorenz1), zeros(4)) lorenz2 = lorenz(:lorenz2) @named connected = System( [s ~ a + lorenz1.x @@ -119,7 +118,7 @@ D = Differential(t) @named sys = System([D(subsys.x) ~ subsys.x + subsys.x], t, systems = [subsys]) sys = mtkcompile(sys) u0 = [subsys.x => 1, subsys.z => 2.0, subsys.y => 1.0] -prob = ODEProblem(sys, u0, (0, 1.0), [subsys.σ => 1, subsys.ρ => 2, subsys.β => 3]) +prob = ODEProblem(sys, [u0; [subsys.σ => 1, subsys.ρ => 2, subsys.β => 3]], (0, 1.0)) sol = solve(prob, FBDF(), reltol = 1e-7, abstol = 1e-7) @test sol[subsys.x] + sol[subsys.y] - sol[subsys.z]≈sol[subsys.u] atol=1e-7 @test_throws ArgumentError convert_system_indepvar(sys, t) @@ -133,7 +132,7 @@ eqs = [0 ~ σ * (y - x), 0 ~ x * y - β * z * h] @named ns = System(eqs, [x, y, z], [σ, ρ, β, h]) np = NonlinearProblem( - complete(ns), [0, 0, 0], [σ => 1, ρ => 2, β => 3], jac = true, sparse = true) + complete(ns), [x => 0, y => 0, z => 0, σ => 1, ρ => 2, β => 3], jac = true, sparse = true) @test calculate_jacobian(ns, sparse = true) isa SparseMatrixCSC # issue #819 @@ -275,12 +274,12 @@ sys = mtkcompile(ns; conservative = true) @test isequal(calculate_jacobian(ns), [(-1 - z + ρ)*σ -x*σ 2x*(-z + ρ) -β-(x^2)]) # solve without analytical jacobian - prob = NonlinearProblem(ns, guesses, ps) + prob = NonlinearProblem(ns, [guesses; ps]) sol = solve(prob, NewtonRaphson()) @test sol.retcode == ReturnCode.Success # solve with analytical jacobian - prob = NonlinearProblem(ns, guesses, ps, jac = true) + prob = NonlinearProblem(ns, [guesses; ps], jac = true) sol = solve(prob, NewtonRaphson()) @test sol.retcode == ReturnCode.Success @@ -327,7 +326,7 @@ end @test length(equations(sys)) == 1 @test length(unknowns(sys)) == 0 T = typeof(ForwardDiff.Dual(1.0)) - prob = NonlinearProblem(sys, [], [p => ForwardDiff.Dual(1.0)]; check_length = false) + prob = NonlinearProblem(sys, [p => ForwardDiff.Dual(1.0)]; check_length = false) @test prob.f(Float64[], prob.p) isa Vector{T} @test prob.f.resid_prototype isa Vector{T} @test_nowarn solve(prob) diff --git a/test/odesystem.jl b/test/odesystem.jl index 463f459df9..c0e11dc718 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -217,15 +217,15 @@ push!(u0, y₂ => 0.0) push!(u0, y₃ => 0.0) p = [k₁ => 0.04, k₃ => 1e4] -p2 = (k₁ => 0.04, +p2 = [k₁ => 0.04, k₂ => 3e7, - k₃ => 1e4) + k₃ => 1e4] tspan = (0.0, 100000.0) -prob1 = ODEProblem(sys, u0, tspan, p) +prob1 = ODEProblem(sys, [u0; p], tspan) @test prob1.f.sys == sys -prob12 = ODEProblem(sys, u0, tspan, [k₁ => 0.04, k₂ => 3e7, k₃ => 1e4]) -prob13 = ODEProblem(sys, u0, tspan, (k₁ => 0.04, k₂ => 3e7, k₃ => 1e4)) -prob14 = ODEProblem(sys, u0, tspan, p2) +prob12 = ODEProblem(sys, [u0; [k₁ => 0.04, k₂ => 3e7, k₃ => 1e4]], tspan) +prob13 = ODEProblem(sys, [u0; [k₁ => 0.04, k₂ => 3e7, k₃ => 1e4]], tspan) +prob14 = ODEProblem(sys, [u0; p2], tspan) for p in [prob1, prob14] @test p.p isa MTKParameters p.ps[k₁] ≈ 0.04 @@ -280,12 +280,12 @@ sol_dpmap = solve(prob_dpmap, Rodas5()) end # test kwargs -prob2 = ODEProblem(sys, u0, tspan, p, jac = true) -prob3 = ODEProblem(sys, u0, tspan, p, jac = true, sparse = true) #SparseMatrixCSC need to handle +prob2 = ODEProblem(sys, [u0; p], tspan, jac = true) +prob3 = ODEProblem(sys, [u0; p], tspan, jac = true, sparse = true) #SparseMatrixCSC need to handle @test prob3.f.jac_prototype isa SparseMatrixCSC -prob3 = ODEProblem(sys, u0, tspan, p, jac = true, sparsity = true) +prob3 = ODEProblem(sys, [u0; p], tspan, jac = true, sparsity = true) @test prob3.f.sparsity isa SparseMatrixCSC -@test_throws ArgumentError ODEProblem(sys, zeros(5), tspan, p) +@test_throws ArgumentError ODEProblem(sys, zeros(5), tspan) for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 1e-12)] local sol sol = solve(prob, Rodas5()) @@ -295,8 +295,8 @@ end du0 = [D(y₁) => -0.04 D(y₂) => 0.04 D(y₃) => 0.0] -prob4 = DAEProblem(sys, du0, u0, tspan, p2) -prob5 = eval(DAEProblem(sys, du0, u0, tspan, p2; expression = Val{true})) +prob4 = DAEProblem(sys, [du0; u0; p2], tspan) +prob5 = eval(DAEProblem(sys, [du0; u0; p2], tspan; expression = Val{true})) for prob in [prob4, prob5] local sol @test prob.differential_vars == [true, true, false] @@ -470,7 +470,7 @@ eqs = [D(x) ~ foo(x, ms); D(ms) ~ bar(ms, p)] @named emptysys = System(Equation[], t) @mtkcompile outersys = compose(emptysys, sys) prob = ODEProblem( - outersys, [sys.x => 1.0, sys.ms => 1:3], (0.0, 1.0), [sys.p => ones(3, 3)]) + outersys, [sys.x => 1.0, sys.ms => 1:3, sys.p => ones(3, 3)], (0.0, 1.0)) @test_nowarn solve(prob, Tsit5()) obsfn = ModelingToolkit.build_explicit_observed_function( outersys, bar(3outersys.sys.ms, 3outersys.sys.p)) @@ -576,18 +576,18 @@ let @named sys = System(eqs, t, vcat(x, [y]), [k], defaults = Dict(x .=> 0)) sys = mtkcompile(sys) - u0 = [0.5, 0] - du0 = 0 .* copy(u0) - prob = DAEProblem(sys, du0, u0, (0, 50)) - @test prob.u0 ≈ u0 - @test prob.du0 ≈ du0 + u0 = unknowns(sys) .=> [0.5, 0] + du0 = D.(unknowns(sys)) .=> 0.0 + prob = DAEProblem(sys, [du0; u0], (0, 50)) + @test prob.u0 ≈ [0.5, 0.0] + @test prob.du0 ≈ [0.0, 0.0] @test prob.p isa MTKParameters @test prob.ps[k] ≈ 1 sol = solve(prob, IDA()) @test sol[y] ≈ 0.9 * sol[x[1]] + sol[x[2]] @test isapprox(sol[x[1]][end], 1, atol = 1e-3) - prob = DAEProblem(sys, [D(y) => 0, D(x[1]) => 0, D(x[2]) => 0], Pair[x[1] => 0.5], + prob = DAEProblem(sys, [D(y) => 0, D(x[1]) => 0, D(x[2]) => 0, x[1] => 0.5], (0, 50)) @test prob.u0 ≈ [0.5, 0] @test prob.du0 ≈ [0, 0] @@ -596,8 +596,8 @@ let sol = solve(prob, IDA()) @test isapprox(sol[x[1]][end], 1, atol = 1e-3) - prob = DAEProblem(sys, [D(y) => 0, D(x[1]) => 0, D(x[2]) => 0], Pair[x[1] => 0.5], - (0, 50), [k => 2]) + prob = DAEProblem(sys, [D(y) => 0, D(x[1]) => 0, D(x[2]) => 0, x[1] => 0.5, k => 2], + (0, 50)) @test prob.u0 ≈ [0.5, 0] @test prob.du0 ≈ [0, 0] @test prob.p isa MTKParameters @@ -607,7 +607,7 @@ let # no initial conditions for D(x[1]) and D(x[2]) provided @test_throws ModelingToolkit.MissingVariablesError prob=DAEProblem( - sys, Pair[], Pair[], (0, 50)) + sys, Pair[], (0, 50)) prob = ODEProblem(sys, Pair[x[1] => 0], (0, 50)) sol = solve(prob, Rosenbrock23()) @@ -621,30 +621,12 @@ let eqs = [D(A) ~ -k1 * k2 * A] @named sys = System(eqs, t) sys = complete(sys) - u0map = [A => 1.0] - pmap = (k1 => 1.0, k2 => 1) + ivmap = [A => 1.0, k1 => 1.0, k2 => 1.0] tspan = (0.0, 1.0) - prob = ODEProblem(sys, u0map, tspan, pmap; tofloat = false) + prob = ODEProblem(sys, ivmap, tspan; tofloat = false) @test prob.p isa MTKParameters @test prob.ps[k1] ≈ 1.0 @test prob.ps[k2] == 1 && prob.ps[k2] isa Int - - prob = ODEProblem(sys, u0map, tspan, pmap) - @test vcat(prob.p...) isa Vector{Float64} - - pmap = [k1 => 1, k2 => 1] - tspan = (0.0, 1.0) - prob = ODEProblem(sys, u0map, tspan, pmap) - @test eltype(vcat(prob.p...)) === Float64 - - prob = ODEProblem(sys, u0map, tspan, pmap) - @test vcat(prob.p...) isa Vector{Float64} - - # No longer supported, Tuple used instead - # pmap = Pair{Any, Union{Int, Float64}}[k1 => 1, k2 => 1.0] - # tspan = (0.0, 1.0) - # prob = ODEProblem(sys, u0map, tspan, pmap) - # @test eltype(prob.p) === Union{Float64, Int} end let @@ -715,7 +697,7 @@ let eqs = [D(A) ~ -k * A] @named osys = System(eqs, t) osys = complete(osys) - oprob = ODEProblem(osys, [A => 1.0], (0.0, 10.0), [k => 1.0]; check_length = false) + oprob = ODEProblem(osys, [A => 1.0, k => 1.0], (0.0, 10.0); check_length = false) @test_nowarn sol = solve(oprob, Tsit5()) end @@ -876,14 +858,14 @@ prob = ODEProblem(sys, [x => 1.0], (0.0, 10.0)) @mtkcompile sys = System( eqs, t; continuous_events = [[norm(x) ~ 3.0] => [x ~ ones(3)]]) # array affect equations used to not work - prob1 = @test_nowarn ODEProblem(sys, [x => ones(3)], (0.0, 10.0), [p => ones(3, 3)]) + prob1 = @test_nowarn ODEProblem(sys, [x => ones(3), p => ones(3, 3)], (0.0, 10.0)) sol1 = @test_nowarn solve(prob1, Tsit5()) # array condition equations also used to not work @mtkcompile sys = System( eqs, t; continuous_events = [[x ~ sqrt(3) * ones(3)] => [x ~ ones(3)]]) # array affect equations used to not work - prob2 = @test_nowarn ODEProblem(sys, [x => ones(3)], (0.0, 10.0), [p => ones(3, 3)]) + prob2 = @test_nowarn ODEProblem(sys, [x => ones(3), p => ones(3, 3)], (0.0, 10.0)) sol2 = @test_nowarn solve(prob2, Tsit5()) @test sol1.u ≈ sol2.u[2:end] @@ -894,7 +876,7 @@ end @variables x(t)[1:3] y(t) @parameters p[1:3, 1:3] @test_nowarn @mtkcompile sys = System([D(x) ~ p * x, D(y) ~ x' * p * x], t) - @test_nowarn ODEProblem(sys, [x => ones(3), y => 2], (0.0, 10.0), [p => ones(3, 3)]) + @test_nowarn ODEProblem(sys, [x => ones(3), y => 2, p => ones(3, 3)], (0.0, 10.0)) end @parameters g L @@ -907,7 +889,7 @@ eqs = [D(D(q₁)) ~ -λ * q₁, @named pend = System(eqs, t) @test_nowarn generate_initializesystem( - pend, u0map = [q₁ => 1.0, q₂ => 0.0], guesses = [λ => 1]) + pend; op = [q₁ => 1.0, q₂ => 0.0], guesses = [λ => 1]) # https://github.com/SciML/ModelingToolkit.jl/issues/2618 @parameters σ ρ β @@ -928,9 +910,9 @@ p = [σ => 28.0, ρ => 10.0, β => 8 / 3] -prob = SteadyStateProblem(sys, u0, p) +prob = SteadyStateProblem(sys, [u0; p]) @test prob isa SteadyStateProblem -prob = SteadyStateProblem(ODEProblem(sys, u0, (0.0, 10.0), p)) +prob = SteadyStateProblem(ODEProblem(sys, [u0; p], (0.0, 10.0))) @test prob isa SteadyStateProblem # Issue#2344 @@ -1048,7 +1030,7 @@ end sys = mtkcompile(System([D(x) ~ P], t, [x], [P]; name = :sys)) function x_at_1(P) - prob = ODEProblem(sys, [x => P], (0.0, 1.0), [sys.P => P]) + prob = ODEProblem(sys, [x => P, sys.P => P], (0.0, 1.0)) return solve(prob, Tsit5())(1.0) end @@ -1095,7 +1077,7 @@ end initialization_eqs = [x ~ T] guesses = [x => 0.0] @named sys2 = System(eqs, T; initialization_eqs, guesses) - prob2 = ODEProblem(mtkcompile(sys2), [], (1.0, 2.0), []) + prob2 = ODEProblem(mtkcompile(sys2), [], (1.0, 2.0)) sol2 = solve(prob2) @test all(sol2[x] .== 1.0) end @@ -1128,7 +1110,7 @@ end defaults = [x => 1, z => y] @named sys = System(eqs, t; defaults) ssys = mtkcompile(sys) - prob = ODEProblem(ssys, [], (0.0, 1.0), []) + prob = ODEProblem(ssys, [], (0.0, 1.0)) @test prob[x] == prob[y] == prob[z] == 1.0 @parameters y0 @@ -1137,7 +1119,7 @@ end defaults = [y0 => 1, x => 1, z => y] @named sys = System(eqs, t; defaults) ssys = mtkcompile(sys) - prob = ODEProblem(ssys, [], (0.0, 1.0), []) + prob = ODEProblem(ssys, [], (0.0, 1.0)) @test prob[x] == prob[y] == prob[z] == 1.0 end @@ -1196,7 +1178,7 @@ end u0 = [sys.y => -1.0, sys.modela.x => -1.0] p = defaults(sys) - prob = ODEProblem(sys, u0, (0.0, 1.0), p) + prob = ODEProblem(sys, merge(p, Dict(u0)), (0.0, 1.0)) # evaluate u0_v = prob.u0 @@ -1336,7 +1318,7 @@ end @variables x(t) @parameters p[1:2] q @mtkcompile sys = System(D(x) ~ sum(p) * x + q * t, t) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p => ones(2), q => 2]) + prob = ODEProblem(sys, [x => 1.0, p => ones(2), q => 2], (0.0, 1.0)) obsfn = ModelingToolkit.build_explicit_observed_function( sys, [p..., q], return_inplace = true)[2] buf = zeros(3) @@ -1391,7 +1373,7 @@ end [c]; discrete_events = [SymbolicDiscreteCallback( 1.0 => [c ~ Pre(c) + 1], discrete_parameters = [c])]) - prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) + prob = ODEProblem(sys, [x => 0.0, c => 1.0], (0.0, 2pi)) sol = solve(prob, Tsit5()) @test sol[obs] ≈ 1:7 end @@ -1400,7 +1382,7 @@ end @variables x(t)=1.0 y(t) [guess = 1.0] @parameters p[1:2] = [1.0, 2.0] @mtkcompile sys = System([D(x) ~ x, y^2 ~ x + sum(p)], t) - prob = DAEProblem(sys, [D(x) => x, D(y) => D(x) / 2y], [], (0.0, 1.0)) + prob = DAEProblem(sys, [D(x) => x, D(y) => D(x) / 2y], (0.0, 1.0)) sol = solve(prob, DFBDF(), abstol = 1e-8, reltol = 1e-8) @test sol[x]≈sol[y^2 - sum(p)] atol=1e-5 end @@ -1423,7 +1405,7 @@ end @mtkcompile sys = System([D(x) ~ p * x + q * t + sum(r), y^3 ~ 2x + 1], t; tstops = [0.5p, [0.1, 0.2], [p + 2q], r]) prob = DAEProblem( - sys, [D(y) => 2D(x) / 3y^2, D(x) => p * x + q * t + sum(r)], [], (0.0, 5.0)) + sys, [D(y) => 2D(x) / 3y^2, D(x) => p * x + q * t + sum(r)], (0.0, 5.0)) sol = solve(prob, DImplicitEuler()) expected_tstops = unique!(sort!(vcat(0.0:0.075:5.0, 0.1, 0.2, 0.65, 0.35, 0.45))) @test all(x -> any(isapprox(x, atol = 1e-6), sol.t), expected_tstops) @@ -1555,7 +1537,7 @@ end β => 8.0f0 / 3.0f0] tspan = (0.0f0, 100.0f0) - prob = ODEProblem{false}(sys, u0, tspan, p) + prob = ODEProblem{false}(sys, [u0; p], tspan) sol = solve(prob, Tsit5()) @test SciMLBase.successful_retcode(sol) end diff --git a/test/optimizationsystem.jl b/test/optimizationsystem.jl index b586d47d4e..0ff85a518d 100644 --- a/test/optimizationsystem.jl +++ b/test/optimizationsystem.jl @@ -28,8 +28,7 @@ using ModelingToolkit, SparseArrays, Test, Optimization, OptimizationOptimJL, generate_cost_hessian(combinedsys) hess_sparsity = ModelingToolkit.cost_hessian_sparsity(sys1) sparse_prob = OptimizationProblem(complete(sys1), - [x => 1, y => 1], - [a => 0.0, b => 0.0], + [x => 1, y => 1, a => 0.0, b => 0.0], grad = true, sparse = true) @test sparse_prob.f.hess_prototype.rowval == hess_sparsity.rowval @@ -46,7 +45,8 @@ using ModelingToolkit, SparseArrays, Test, Optimization, OptimizationOptimJL, sys2.b => 9.0 β => 10.0] - prob = OptimizationProblem(combinedsys, u0, p, grad = true, hess = true, cons_j = true, + prob = OptimizationProblem( + combinedsys, [u0; p], grad = true, hess = true, cons_j = true, cons_h = true) @test prob.f.sys === combinedsys sol = solve(prob, Ipopt.Optimizer(); print_level = 0) @@ -62,7 +62,7 @@ end ] @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons) sys = complete(sys) - prob = OptimizationProblem(sys, [x => 0.0, y => 0.0], [a => 1.0, b => 1.0], + prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, a => 1.0, b => 1.0], grad = true, hess = true, cons_j = true, cons_h = true) @test prob.f.sys === sys sol = solve(prob, IPNewton()) @@ -70,7 +70,7 @@ end sol = solve(prob, Ipopt.Optimizer(); print_level = 0) @test sol.objective < 1.0 - prob = OptimizationProblem(sys, [x => 0.0, y => 0.0], [a => 1.0, b => 1.0], + prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, a => 1.0, b => 1.0], grad = false, hess = false, cons_j = false, cons_h = false) sol = solve(prob, AmplNLWriter.Optimizer(Ipopt_jll.amplexe)) @test_skip sol.objective < 1.0 @@ -85,7 +85,7 @@ end z^2 + y^2 ≲ 1.0] @named sys = OptimizationSystem(loss, [x, y, z], [a, b], constraints = cons) sys = mtkcompile(sys) - prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, z => 0.0], [a => 1.0, b => 1.0], + prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, z => 0.0, a => 1.0, b => 1.0], grad = true, hess = true, cons_j = true, cons_h = true) sol = solve(prob, IPNewton()) @test sol.objective < 1.0 @@ -96,9 +96,9 @@ end @test sol.u≈[0.808, -0.064] atol=1e-3 @test sol[x]^2 + sol[y]^2 ≈ 1.0 - prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, z => 0.0], [a => 1.0, b => 1.0], + prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, z => 0.0, a => 1.0, b => 1.0], grad = false, hess = false, cons_j = false, cons_h = false) - sol = solve(prob, AmplNLWriter.Optimizer(Ipopt_jll.amplexe)) + @test_broken sol = solve(prob, AmplNLWriter.Optimizer(Ipopt_jll.amplexe)) @test_skip sol.objective < 1.0 @test_skip sol.u≈[0.808, -0.064] atol=1e-3 @test_skip sol[x]^2 + sol[y]^2 ≈ 1.0 @@ -109,9 +109,11 @@ end x0 = zeros(2) p = [1.0, 100.0] f = OptimizationFunction(rosenbrock, Optimization.AutoSymbolics()) - prob = OptimizationProblem(f, x0, p) - sol = solve(prob, Newton()) - @test sol.u ≈ [1.0, 1.0] + @test_broken begin + prob = OptimizationProblem(f, x0, p) + sol = solve(prob, Newton()) + @test sol.u ≈ [1.0, 1.0] + end end # issue #819 @@ -208,7 +210,8 @@ end sys2.b => 9.0 β => 10.0] - prob = OptimizationProblem(combinedsys, u0, p, grad = true, hess = true, cons_j = true, + prob = OptimizationProblem( + combinedsys, [u0; p], grad = true, hess = true, cons_j = true, cons_h = true) @test prob.f.sys === combinedsys @test_broken SciMLBase.successful_retcode(solve(prob, @@ -235,7 +238,7 @@ end x[1]^2 + x[2]^2 ≲ 2.0 ]) - prob = OptimizationProblem(complete(sys), [x[1] => 2.0, x[2] => 0.0], [], grad = true, + prob = OptimizationProblem(complete(sys), [x[1] => 2.0, x[2] => 0.0], grad = true, hess = true, cons_j = true, cons_h = true) sol = Optimization.solve(prob, Ipopt.Optimizer(); print_level = 0) @test sol.u ≈ [1, 0] @@ -249,7 +252,7 @@ end @parameters a b loss = (a - x)^2 + b * (y - x^2)^2 @named sys = OptimizationSystem(loss, [x, y], [a, b, c]) - prob = OptimizationProblem(complete(sys), [x => 0.0, y => 0.0], [a => 1.0, b => 100.0]) + prob = OptimizationProblem(complete(sys), [x => 0.0, y => 0.0, a => 1.0, b => 100.0]) @test prob.lb == [-Inf, 0.0] @test prob.ub == [Inf, Inf] end @@ -267,11 +270,11 @@ end name = :sys1, constraints = cons)) - prob1 = OptimizationProblem(sys1, [x₁ => 0.0, x₂ => 0.0], [α₁ => 1.0, α₂ => 100.0], + prob1 = OptimizationProblem(sys1, [x₁ => 0.0, x₂ => 0.0, α₁ => 1.0, α₂ => 100.0], grad = true, hess = true, cons_j = true, cons_h = true) sys2 = complete(modelingtoolkitize(prob1)) - prob2 = OptimizationProblem(sys2, [x₁ => 0.0, x₂ => 0.0], [α₁ => 1.0, α₂ => 100.0], + prob2 = OptimizationProblem(sys2, [x₁ => 0.0, x₂ => 0.0, α₁ => 1.0, α₂ => 100.0], grad = true, hess = true, cons_j = true, cons_h = true) sol1 = Optimization.solve(prob1, Ipopt.Optimizer()) @@ -287,7 +290,7 @@ end @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = [x^2 + y^2 ≲ 0.0]) sys = complete(sys) - prob = OptimizationProblem(sys, [x => 0.0, y => 0.0], [a => 1.0, b => 100.0]) + prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, a => 1.0, b => 100.0]) @test prob.f.expr isa Symbolics.Symbolic @test all(prob.f.cons_expr[i].lhs isa Symbolics.Symbolic for i in 1:length(prob.f.cons_expr)) @@ -300,7 +303,7 @@ end cons2 = [x^2 + y^2 ~ 0, y * sin(x) - x ~ 0] sys = complete(OptimizationSystem( loss, [x, y], [a, b], name = :sys2, constraints = cons2)) - prob = OptimizationProblem(sys, [x => 0.0, y => 0.0], [a => 1.0, b => 100.0], + prob = OptimizationProblem(sys, [x => 0.0, y => 0.0, a => 1.0, b => 100.0], grad = true, hess = true, cons_j = true, cons_h = true) G1 = Array{Float64}(undef, 2) @@ -354,7 +357,7 @@ end @parameters p f(::Real) @mtkcompile sys = OptimizationSystem( x^2 + f(x) * p, [x], [f, p]; constraints = [2.0 ≲ f(x) + p]) - prob = OptimizationProblem(sys, [x => 1.0], [p => 1.0, f => (x -> 2x)]) + prob = OptimizationProblem(sys, [x => 1.0, p => 1.0, f => (x -> 2x)]) @test abs(prob.f.cons(prob.u0, prob.p)[1]) ≈ 1.0 end diff --git a/test/parameter_dependencies.jl b/test/parameter_dependencies.jl index 97937b8811..0abfd2e63e 100644 --- a/test/parameter_dependencies.jl +++ b/test/parameter_dependencies.jl @@ -33,7 +33,7 @@ using NonlinearSolve @test prob.ps[p1] == 1.0 @test prob.ps[p2] == 2.0 @test SciMLBase.successful_retcode(solve(prob, Tsit5())) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.5), [p1 => 1.0], jac = true) + prob = ODEProblem(sys, [x => 1.0, p1 => 1.0], (0.0, 1.5), jac = true) @test prob.ps[p1] == 1.0 @test prob.ps[p2] == 2.0 integ = init(prob, Tsit5()) @@ -220,23 +220,26 @@ end @test_skip begin Tf = 1.0 - prob = ODEProblem(sys, [x => 0.0, y => 0.0], (0.0, Tf), - [kp => 1.0; z(k - 1) => 3.0; yd(k - 1) => 0.0; z(k - 2) => 4.0; - yd(k - 2) => 2.0]) + prob = ODEProblem(sys, + [x => 0.0, y => 0.0, kp => 1.0, z(k - 1) => 3.0, + yd(k - 1) => 0.0, z(k - 2) => 4.0, yd(k - 2) => 2.0], + (0.0, Tf)) @test_nowarn solve(prob, Tsit5()) @mtkcompile sys = System(eqs, t; parameter_dependencies = [kq => 2kp], discrete_events = [SymbolicDiscreteCallback( [0.5] => [kp ~ 2.0], discrete_parameters = [kp])]) - prob = ODEProblem(sys, [x => 0.0, y => 0.0], (0.0, Tf), - [kp => 1.0; z(k - 1) => 3.0; yd(k - 1) => 0.0; z(k - 2) => 4.0; - yd(k - 2) => 2.0]) + prob = ODEProblem(sys, + [x => 0.0, y => 0.0, kp => 1.0, z(k - 1) => 3.0, + yd(k - 1) => 0.0, z(k - 2) => 4.0, yd(k - 2) => 2.0], + (0.0, Tf)) @test prob.ps[kp] == 1.0 @test prob.ps[kq] == 2.0 @test_nowarn solve(prob, Tsit5()) - prob = ODEProblem(sys, [x => 0.0, y => 0.0], (0.0, Tf), - [kp => 1.0; z(k - 1) => 3.0; yd(k - 1) => 0.0; z(k - 2) => 4.0; - yd(k - 2) => 2.0]) + prob = ODEProblem(sys, + [x => 0.0, y => 0.0, kp => 1.0, z(k - 1) => 3.0, + yd(k - 1) => 0.0, z(k - 2) => 4.0, yd(k - 2) => 2.0], + (0.0, Tf)) integ = init(prob, Tsit5()) @test integ.ps[kp] == 1.0 @test integ.ps[kq] == 2.0 @@ -265,7 +268,7 @@ end @test ρ in Set(full_parameters(sdesys)) prob = SDEProblem( - sdesys, [x => 1.0, y => 0.0, z => 0.0], (0.0, 100.0), [σ => 10.0, β => 2.33]) + sdesys, [x => 1.0, y => 0.0, z => 0.0, σ => 10.0, β => 2.33], (0.0, 100.0)) @test prob.ps[ρ] == 2prob.ps[σ] @test_nowarn solve(prob, SRIW1()) @@ -275,7 +278,7 @@ end [10.0] => [σ ~ 15.0], discrete_parameters = [σ])]) sdesys = complete(sdesys) prob = SDEProblem( - sdesys, [x => 1.0, y => 0.0, z => 0.0], (0.0, 100.0), [σ => 10.0, β => 2.33]) + sdesys, [x => 1.0, y => 0.0, z => 0.0, σ => 10.0, β => 2.33], (0.0, 100.0)) integ = init(prob, SRIW1()) @test integ.ps[σ] == 10.0 @test integ.ps[ρ] == 20.0 @@ -303,7 +306,7 @@ end tspan = (0.0, 250.0) u₀map = [S => 999, I => 1, R => 0] parammap = [γ => 0.01] - jprob = JumpProblem(js2, u₀map, tspan, parammap; aggregator = Direct(), + jprob = JumpProblem(js2, [u₀map; parammap], tspan; aggregator = Direct(), save_positions = (false, false), rng = rng) @test jprob.ps[γ] == 0.01 @test jprob.ps[β] == 0.0001 @@ -314,7 +317,7 @@ end discrete_events = [SymbolicDiscreteCallback( [10.0] => [γ ~ 0.02], discrete_parameters = [γ])]) js2 = complete(js2) - jprob = JumpProblem(js2, u₀map, tspan, parammap; aggregator = Direct(), + jprob = JumpProblem(js2, [u₀map; parammap], tspan; aggregator = Direct(), save_positions = (false, false), rng = rng) integ = init(jprob, SSAStepper()) @test integ.ps[γ] == 0.01 @@ -335,7 +338,7 @@ end @test prob.ps[p1] == 1.0 @test prob.ps[p2] == 2.0 @test_nowarn solve(prob, NewtonRaphson()) - prob = NonlinearProblem(sys, [x => 1.0], [p1 => 2.0]) + prob = NonlinearProblem(sys, [x => 1.0, p1 => 2.0]) @test prob.ps[p1] == 2.0 @test prob.ps[p2] == 4.0 end @@ -355,7 +358,7 @@ end t; parameter_dependencies = [p2 => 2p1] ) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.5), [p1 => 1.0], jac = true) + prob = ODEProblem(sys, [x => 1.0, p1 => 1.0], (0.0, 1.5), jac = true) prob.ps[p1] = 3.0 @test prob.ps[p1] == 3.0 @test prob.ps[p2] == 6.0 diff --git a/test/reduction.jl b/test/reduction.jl index af9e510646..43fef0fb63 100644 --- a/test/reduction.jl +++ b/test/reduction.jl @@ -115,10 +115,10 @@ u0 = [lorenz1.x => 1.0 lorenz2.x => 1.0 lorenz2.y => 0.0 lorenz2.z => 0.0] -prob1 = ODEProblem(reduced_system, u0, (0.0, 100.0), pp) +prob1 = ODEProblem(reduced_system, [u0; pp], (0.0, 100.0)) solve(prob1, Rodas5()) -prob2 = SteadyStateProblem(reduced_system, u0, pp) +prob2 = SteadyStateProblem(reduced_system, [u0; pp]) @test prob2.f.observed(lorenz2.u, prob2.u0, prob2.p) === 1.0 # issue #724 and #716 @@ -160,7 +160,7 @@ reducedsys = mtkcompile(sys) u0 = [u2 => 1] pp = [2] -nlprob = NonlinearProblem(reducedsys, u0, [p => pp[1]]) +nlprob = NonlinearProblem(reducedsys, [u0; [p => pp[1]]]) reducedsol = solve(nlprob, NewtonRaphson()) residual = fill(100.0, length(unknowns(reducedsys))) nlprob.f(residual, reducedsol.u, pp) diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl index 278b1a90d1..ded9852670 100644 --- a/test/scc_nonlinear_problem.jl +++ b/test/scc_nonlinear_problem.jl @@ -86,7 +86,7 @@ end @parameters p1[1:6] p2 eqs = 0 .~ collect(nlf(u, (u0, (p1, p2)))) @mtkcompile sys = System(eqs, [u], [p1, p2]) - sccprob = SCCNonlinearProblem(sys, [u => u0], [p1 => p[1], p2 => p[2][]]) + sccprob = SCCNonlinearProblem(sys, [u => u0, p1 => p[1], p2 => p[2][]]) sccsol = solve(sccprob, SimpleNewtonRaphson(); abstol = 1e-9) @test SciMLBase.successful_retcode(sccsol) @test norm(sccsol.resid) < norm(sol.resid) @@ -142,10 +142,10 @@ end subrules = Dict(Symbolics.unwrap(D(y[i])) => ((y[i] - u0[i]) / dt) for i in 1:8) eqs = substitute.(eqs, (subrules,)) @mtkcompile sys = System(eqs) - prob = NonlinearProblem(sys, [y => u0], [t => t0]) + prob = NonlinearProblem(sys, [y => u0, t => t0]) sol = solve(prob, NewtonRaphson(); abstol = 1e-12) - sccprob = SCCNonlinearProblem(sys, [y => u0], [t => t0]) + sccprob = SCCNonlinearProblem(sys, [y => u0, t => t0]) sccsol = solve(sccprob, NewtonRaphson(); abstol = 1e-12) @test sol.u≈sccsol.u atol=1e-10 @@ -266,7 +266,7 @@ end @parameters (f::Function)(..) @mtkcompile sys = System([ 0 ~ x[1]^2 - 9, x[2] ~ 2x[1], 0 ~ x[3]^2 - x[1]^2 + f(x)]) - prob = SCCNonlinearProblem(sys, [x => ones(3)], [f => sum]) + prob = SCCNonlinearProblem(sys, [x => ones(3), f => sum]) sol = solve(prob, NewtonRaphson()) @test SciMLBase.successful_retcode(sol) end @@ -286,7 +286,7 @@ end ρ => 10.0, β => 8 / 3] - sccprob = SCCNonlinearProblem(fullsys, u0, p) + sccprob = SCCNonlinearProblem(fullsys, [u0; p]) @test isequal(parameters(fullsys), parameters(sccprob.f.sys)) end @@ -295,6 +295,6 @@ end @parameters p[1:2] (f::Function)(..) @mtkcompile sys = System([x^2 - p[1]^2 ~ 0, y^2 ~ f(p)]) - prob = SCCNonlinearProblem(sys, [x => 1.0, y => 1.0], [p => ones(2), f => sum]) + prob = SCCNonlinearProblem(sys, [x => 1.0, y => 1.0, p => ones(2), f => sum]) @test_nowarn solve(prob, NewtonRaphson()) end diff --git a/test/sciml_problem_inputs.jl b/test/sciml_problem_inputs.jl index 249991b631..c02bd9dea1 100644 --- a/test/sciml_problem_inputs.jl +++ b/test/sciml_problem_inputs.jl @@ -103,7 +103,7 @@ end # Perform ODE simulations (singular and ensemble). let # Creates normal and ensemble problems. - base_oprob = ODEProblem(osys, u0_alts[1], tspan, p_alts[1]) + base_oprob = ODEProblem(osys, [u0_alts[1]; p_alts[1]], tspan) base_sol = solve(base_oprob, Tsit5(); saveat = 1.0) base_eprob = EnsembleProblem(base_oprob) base_esol = solve(base_eprob, Tsit5(); trajectories = 2, saveat = 1.0) @@ -120,7 +120,7 @@ end # Solves a nonlinear problem (EnsembleProblems are not possible for these). let - base_nlprob = NonlinearProblem(nsys, u0_alts[1], p_alts[1]) + base_nlprob = NonlinearProblem(nsys, [u0_alts[1]; p_alts[1]]) base_sol = solve(base_nlprob, NewtonRaphson()) # Solves problems for all input types, checking that identical solutions are found. for u0 in u0_alts, p in p_alts @@ -132,7 +132,7 @@ end # Perform steady state simulations (singular and ensemble). let # Creates normal and ensemble problems. - base_ssprob = SteadyStateProblem(osys, u0_alts[1], p_alts[1]) + base_ssprob = SteadyStateProblem(osys, [u0_alts[1]; p_alts[1]]) base_sol = solve(base_ssprob, DynamicSS(Tsit5())) base_eprob = EnsembleProblem(base_ssprob) base_esol = solve(base_eprob, DynamicSS(Tsit5()); trajectories = 2) diff --git a/test/sdesystem.jl b/test/sdesystem.jl index bcec551faf..d2a8f77636 100644 --- a/test/sdesystem.jl +++ b/test/sdesystem.jl @@ -27,11 +27,14 @@ f = eval(generate_diffusion_function(de)[1]) @test f(ones(3), rand(3), nothing) == 0.1ones(3) f = SDEFunction(de) -prob = SDEProblem(de, [1.0, 0.0, 0.0], (0.0, 100.0), [10.0, 26.0, 2.33]) +prob = SDEProblem( + de, [unknowns(de) .=> [1.0, 0.0, 0.0]; parameters(de) .=> [10.0, 26.0, 2.33]], + (0.0, 100.0)) sol = solve(prob, SRIW1(), seed = 1) -probexpr = SDEProblem(de, [1.0, 0.0, 0.0], (0.0, 100.0), - [10.0, 26.0, 2.33]) +probexpr = SDEProblem( + de, [unknowns(de) .=> [1.0, 0.0, 0.0]; parameters(de) .=> [10.0, 26.0, 2.33]], + (0.0, 100.0)) solexpr = solve(eval(probexpr), SRIW1(), seed = 1) @test all(x -> x == 0, Array(sol - solexpr)) @@ -54,8 +57,8 @@ f(du, [1, 2, 3.0], p, nothing) 0.1 0.01*2 0.02*1*3 0.2 0.3 0.01*3] -prob = SDEProblem(de, [1.0, 0.0, 0.0], (0.0, 100.0), (σ => 10.0, ρ => 26.0, β => 2.33), - noise_rate_prototype = zeros(3, 3)) +prob = SDEProblem(de, [x => 1.0, y => 0.0, z => 0.0, σ => 10.0, ρ => 26.0, β => 2.33], + (0.0, 100.0), noise_rate_prototype = zeros(3, 3)) sol = solve(prob, EM(), dt = 0.001) u0map = [ @@ -70,13 +73,13 @@ parammap = [ ρ => 2.33 ] -prob = SDEProblem(de, u0map, (0.0, 100.0), parammap) +prob = SDEProblem(de, [u0map; parammap], (0.0, 100.0)) @test prob.f.sys === de @test size(prob.noise_rate_prototype) == (3, 3) @test prob.noise_rate_prototype isa Matrix sol = solve(prob, EM(), dt = 0.001) -prob = SDEProblem(de, u0map, (0.0, 100.0), parammap, sparsenoise = true) +prob = SDEProblem(de, [u0map; parammap], (0.0, 100.0), sparsenoise = true) @test size(prob.noise_rate_prototype) == (3, 3) @test prob.noise_rate_prototype isa SparseMatrixCSC sol = solve(prob, EM(), dt = 0.001) @@ -506,14 +509,14 @@ noiseeqs = [0.1 * x] @named de = SDESystem(eqs, noiseeqs, tt, [x], [α, β], observed = [weight ~ x * 10]) de = complete(de) - prob = SDEProblem(de, u0map, (0.0, 1.0), parammap) + prob = SDEProblem(de, [u0map; parammap], (0.0, 1.0)) sol = solve(prob, EM(), dt = dt) @test observed(de) == [weight ~ x * 10] @test sol[weight] == 10 * sol[x] @named ode = System(eqs, tt, [x], [α, β], observed = [weight ~ x * 10]) ode = complete(ode) - odeprob = ODEProblem(ode, u0map, (0.0, 1.0), parammap) + odeprob = ODEProblem(ode, [u0map; parammap], (0.0, 1.0)) solode = solve(odeprob, Tsit5()) @test observed(ode) == [weight ~ x * 10] @test solode[weight] == 10 * solode[x] @@ -545,7 +548,7 @@ end β => 1.0 ] - prob = SDEProblem(de, u0map, (0.0, 1.0), parammap) + prob = SDEProblem(de, [u0map; parammap], (0.0, 1.0)) function prob_func(prob, i, repeat) remake(prob, seed = seeds[i]) @@ -567,7 +570,7 @@ end u = x demod = complete(ModelingToolkit.Girsanov_transform(de, u; θ0 = 0.1)) - probmod = SDEProblem(demod, u0map, (0.0, 1.0), parammap) + probmod = SDEProblem(demod, [u0map; parammap], (0.0, 1.0)) ensemble_probmod = EnsembleProblem(probmod; output_func = (sol, i) -> (g(sol[x, end]) * @@ -613,8 +616,8 @@ sys2 = complete(sys2) @set! sys2.parent = nothing @test sys1 == sys2 -prob = SDEProblem(sys1, sts .=> [1.0, 0.0, 0.0], - (0.0, 100.0), ps .=> (10.0, 26.0)) +prob = SDEProblem(sys1, [sts .=> [1.0, 0.0, 0.0]; ps .=> [10.0, 26.0]], + (0.0, 100.0)) solve(prob, LambaEulerHeun(), seed = 1) # Test ill-formed due to more equations than states in noise equations @@ -640,7 +643,7 @@ eqs = [D(x) ~ p - d * x + a * sqrt(p)] u0 = @SVector[x => 10.0] tspan = (0.0, 10.0) ps = @SVector[p => 5.0, d => 0.5] -sprob = SDEProblem(sys, u0, tspan, ps) +sprob = SDEProblem(sys, [u0; ps], tspan) @test !isinplace(sprob) @test !isinplace(sprob.f) @test_nowarn solve(sprob, ImplicitEM()) @@ -654,7 +657,7 @@ eqs = [D(x) ~ p - d * x + a * sqrt(p) u0 = @SVector[x => 10.0, y => 20.0] tspan = (0.0, 10.0) ps = @SVector[p => 5.0, d => 0.5] -sprob = SDEProblem(sys, u0, tspan, ps) +sprob = SDEProblem(sys, [u0; ps], tspan) @test sprob.f.g(sprob.u0, sprob.p, sprob.tspan[1]) isa SVector{2, Float64} @test_nowarn solve(sprob, ImplicitEM()) @@ -679,7 +682,7 @@ let β => 26.0, ρ => 2.33 ] - prob = SDEProblem(de, u0map, (0.0, 100.0), parammap) + prob = SDEProblem(de, [u0map; parammap], (0.0, 100.0)) # TODO: re-enable this when we support scalar noise @test solve(prob, SOSRI()).retcode == ReturnCode.Success end @@ -691,7 +694,7 @@ let # test to make sure that scalar noise always receive the same kicks D(y) ~ a] @mtkcompile de = System(eqs, t) - prob = SDEProblem(de, [x => 0, y => 0], (0.0, 10.0), []) + prob = SDEProblem(de, [x => 0, y => 0], (0.0, 10.0)) sol = solve(prob, SOSRI()) @test sol.u[end][1] == sol.u[end][2] end @@ -718,7 +721,7 @@ let # test that diagonal noise is correctly handled ρ => 2.33 ] - prob = SDEProblem(de, u0map, (0.0, 100.0), parammap) + prob = SDEProblem(de, [u0map; parammap], (0.0, 100.0)) # SOSRI only works for diagonal and scalar noise @test solve(prob, SOSRI()).retcode == ReturnCode.Success end @@ -744,7 +747,7 @@ end ρ => 2.33 ] - prob = SDEProblem(de, u0map, (0.0, 100.0), parammap) + prob = SDEProblem(de, [u0map; parammap], (0.0, 100.0)) # SOSRI only works for diagonal and scalar noise @test_throws ErrorException solve(prob, SOSRI()).retcode==ReturnCode.Success # ImplicitEM does work for non-diagonal noise @@ -773,7 +776,7 @@ end ρ => 2.33 ] - prob = SDEProblem(de, u0map, (0.0, 100.0), parammap) + prob = SDEProblem(de, [u0map; parammap], (0.0, 100.0)) @test solve(prob, SOSRI()).retcode == ReturnCode.Success end @@ -800,7 +803,7 @@ end sys = System(eqs, t, sts, ps, browns; name = :name) sys = mtkcompile(sys) @test ModelingToolkit.get_noise_eqs(sys) ≈ [1.0] - prob = SDEProblem(sys, [], (0.0, 1.0), []) + prob = SDEProblem(sys, [], (0.0, 1.0)) @test_nowarn solve(prob, RKMil()) end @@ -964,7 +967,7 @@ end parammap = [σ => 10.0, β => 26.0, ρ => 2.33] @test_throws ["Brownian", "mtkcompile"] SDEProblem( - de, u0map, (0.0, 100.0), parammap) + de, [u0map; parammap], (0.0, 100.0)) de = mtkcompile(de) - @test SDEProblem(de, u0map, (0.0, 100.0), parammap) isa SDEProblem + @test SDEProblem(de, [u0map; parammap], (0.0, 100.0)) isa SDEProblem end diff --git a/test/serialization.jl b/test/serialization.jl index 83e68f5770..41ea9c2c15 100644 --- a/test/serialization.jl +++ b/test/serialization.jl @@ -6,10 +6,8 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @named sys = System([D(x) ~ -0.5 * x], t, defaults = Dict(x => 1.0)) sys = complete(sys) for prob in [ - eval(ModelingToolkit.ODEProblem{false}(sys, nothing, nothing, - SciMLBase.NullParameters())), - eval(ModelingToolkit.ODEProblem{false}(sys, nothing, nothing, - SciMLBase.NullParameters(); expression = Val{true})) + eval(ModelingToolkit.ODEProblem{false}(sys, nothing, nothing)), + eval(ModelingToolkit.ODEProblem{false}(sys, nothing, nothing; expression = Val{true})) ] _fn = tempname() diff --git a/test/split_parameters.jl b/test/split_parameters.jl index ac819351b5..f76ae4f01c 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -84,7 +84,7 @@ eqs = [y ~ src.output.u s = complete(sys) sys = mtkcompile(sys) prob = ODEProblem( - sys, [], (0.0, t_end), [s.src.interpolator => Interpolator(x, dt)]; + sys, [s.src.interpolator => Interpolator(x, dt)], (0.0, t_end); tofloat = false) sol = solve(prob, ImplicitEuler()); @test sol.retcode == ReturnCode.Success @@ -111,7 +111,7 @@ eqs = [D(y) ~ dy * a sys = mtkcompile(model; split = false) tspan = (0.0, t_end) -prob = ODEProblem(sys, [], tspan, []; build_initializeprob = false) +prob = ODEProblem(sys, [], tspan; build_initializeprob = false) @test prob.p isa Vector{Float64} sol = solve(prob, ImplicitEuler()); @@ -120,7 +120,7 @@ sol = solve(prob, ImplicitEuler()); # ------------------------ Mixed Type Conserved prob = ODEProblem( - sys, [], tspan, []; tofloat = false, build_initializeprob = false) + sys, [], tspan; tofloat = false, build_initializeprob = false) sol = solve(prob, ImplicitEuler()); @test sol.retcode == ReturnCode.Success @@ -248,7 +248,7 @@ end sol = solve(prob, Tsit5(); abstol = 1e-10, reltol = 1e-10) @test sol.u[end][] ≈ 2.0 - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [fn => Foo()]) + prob = ODEProblem(sys, [x => 1.0, fn => Foo()], (0.0, 1.0)) @inferred getter(prob) @inferred Vector{<:Real} prob.f(prob.u0, prob.p, prob.tspan[1]) sol = solve(prob; abstol = 1e-10, reltol = 1e-10) @@ -263,7 +263,7 @@ end @mtkcompile sys = System(D(x) ~ fn(x), t) @test is_parameter(sys, fn) getter = getp(sys, fn) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [fn => interp]) + prob = ODEProblem(sys, [x => 1.0, fn => interp], (0.0, 1.0)) @inferred getter(prob) @inferred prob.f(prob.u0, prob.p, prob.tspan[1]) @test_nowarn sol = solve(prob, Tsit5()) diff --git a/test/state_selection.jl b/test/state_selection.jl index f3971bdbe3..97741110e2 100644 --- a/test/state_selection.jl +++ b/test/state_selection.jl @@ -122,9 +122,9 @@ let sys.supply_pipe.v => 0.1, sys.return_pipe.v => 0.1, D(supply_pipe.v) => 0.0, D(return_pipe.fluid_port_a.m) => 0.0, D(supply_pipe.fluid_port_a.m) => 0.0] - prob1 = ODEProblem(sys, [], (0.0, 10.0), [], guesses = u0) - prob2 = ODEProblem(sys, [], (0.0, 10.0), [], guesses = u0) - prob3 = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, [], (0.0, 10.0), guesses = u0) + prob1 = ODEProblem(sys, [], (0.0, 10.0), guesses = u0) + prob2 = ODEProblem(sys, [], (0.0, 10.0), guesses = u0) + prob3 = DAEProblem(sys, D.(unknowns(sys)) .=> 0.0, (0.0, 10.0), guesses = u0) @test solve(prob1, FBDF()).retcode == ReturnCode.Success #@test solve(prob2, FBDF()).retcode == ReturnCode.Success @test solve(prob3, DFBDF()).retcode == ReturnCode.Success @@ -274,6 +274,6 @@ let # solution ------------------------------------------------------------------- @named catapult = System(eqs, t, vars, params, defaults = defs) sys = mtkcompile(catapult) - prob = ODEProblem(sys, [], (0.0, 0.1), [l_2f => 0.55, damp => 1e7]; jac = true) + prob = ODEProblem(sys, [l_2f => 0.55, damp => 1e7], (0.0, 0.1); jac = true) @test solve(prob, Rodas4()).retcode == ReturnCode.Success end diff --git a/test/static_arrays.jl b/test/static_arrays.jl index edb6eeff7d..a23eeddde1 100644 --- a/test/static_arrays.jl +++ b/test/static_arrays.jl @@ -11,17 +11,16 @@ eqs = [D(D(x)) ~ σ * (y - x), @named sys = System(eqs, t) sys = mtkcompile(sys) -u0 = @SVector [D(x) => 2.0, +ivs = @SVector [D(x) => 2.0, x => 1.0, y => 0.0, - z => 0.0] - -p = @SVector [σ => 28.0, + z => 0.0, + σ => 28.0, ρ => 10.0, β => 8 / 3] tspan = (0.0, 100.0) -prob_mtk = ODEProblem(sys, u0, tspan, p) +prob_mtk = ODEProblem(sys, ivs, tspan) @test !SciMLBase.isinplace(prob_mtk) @test prob_mtk.u0 isa SArray diff --git a/test/steadystatesystems.jl b/test/steadystatesystems.jl index a99e84d9d2..505e7da890 100644 --- a/test/steadystatesystems.jl +++ b/test/steadystatesystems.jl @@ -14,10 +14,10 @@ for factor in [1e-1, 1e0, 1e10], u0 = [x => factor * u0_p[1]] p = [r => factor * u0_p[2]] - ss_prob = SteadyStateProblem(de, u0, p) + ss_prob = SteadyStateProblem(de, [u0; p]) sol = solve(ss_prob, SSRootfind()).u[1] @test abs(sol^2 - factor * u0_p[2]) < 1e-8 - ss_prob = SteadyStateProblem(de, u0, p) + ss_prob = SteadyStateProblem(de, [u0; p]) sol_expr = solve(ss_prob, SSRootfind()).u[1] @test all(x -> x == 0, sol - sol_expr) end diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl index 6ab8dde338..94dd1f884b 100644 --- a/test/structural_transformation/index_reduction.jl +++ b/test/structural_transformation/index_reduction.jl @@ -53,16 +53,14 @@ let sys = mtkcompile(pendulum2) @test length(equations(sys)) == 5 @test length(unknowns(sys)) == 5 - u0 = [ + ivs = [ x => sqrt(2) / 2, - y => sqrt(2) / 2 - ] - p = [ + y => sqrt(2) / 2, L => 1.0, g => 9.8 ] - prob_auto = ODEProblem(sys, u0, (0.0, 0.5), p, guesses = [T => 0.0]) + prob_auto = ODEProblem(sys, ivs, (0.0, 0.5), guesses = [T => 0.0]) sol = solve(prob_auto, FBDF()) @test sol.retcode == ReturnCode.Success @test norm(sol[x] .^ 2 + sol[y] .^ 2 .- 1) < 1e-2 @@ -78,7 +76,7 @@ let @named pend = System(eqs, t) sys = complete(mtkcompile(pend; dummy_derivative = false)) prob = ODEProblem( - sys, [x => 1, y => 0, D(x) => 0.0], (0.0, 10.0), [g => 1], guesses = [λ => 0.0]) + sys, [x => 1, y => 0, D(x) => 0.0, g => 1], (0.0, 10.0), guesses = [λ => 0.0]) sol = solve(prob, Rodas5P()) @test SciMLBase.successful_retcode(sol) @test sol[x^2 + y^2][end] < 1.1 diff --git a/test/structural_transformation/tearing.jl b/test/structural_transformation/tearing.jl index e9b9909871..7f18cbb703 100644 --- a/test/structural_transformation/tearing.jl +++ b/test/structural_transformation/tearing.jl @@ -150,7 +150,7 @@ newdaesys = mtkcompile(daesys) @test equations(newdaesys) == [D(x) ~ h * z; 0 ~ y + sin(z) - p * t] @test equations(tearing_substitution(newdaesys)) == [D(x) ~ h * z; 0 ~ x + sin(z) - p * t] @test isequal(unknowns(newdaesys), [x, z]) -prob = ODEProblem(newdaesys, [x => 1.0, z => -0.5π], (0, 1.0), [p => 0.2]) +prob = ODEProblem(newdaesys, [x => 1.0, z => -0.5π, p => 0.2], (0, 1.0)) du = [0.0, 0.0]; u = [1.0, -0.5π]; pr = prob.p; @@ -161,7 +161,7 @@ prob.f(du, u, pr, tt) # test the initial guess is respected @named sys = System(eqs, t, defaults = Dict(z => NaN)) -infprob = ODEProblem(mtkcompile(sys), [x => 1.0], (0, 1.0), [p => 0.2]) +infprob = ODEProblem(mtkcompile(sys), [x => 1.0, p => 0.2], (0, 1.0)) infprob.f(du, infprob.u0, pr, tt) @test any(isnan, du) diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 02298eb480..4d5dead59f 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -55,7 +55,7 @@ end @test length(observed(sys)) == 7 @test any(obs -> isequal(obs, y), observables(sys)) @test any(obs -> isequal(obs, z), observables(sys)) - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn]) + prob = ODEProblem(sys, [x => 1.0, foo => _tmp_fn], (0.0, 1.0)) @test_nowarn prob.f(prob.u0, prob.p, 0.0) isys = ModelingToolkit.generate_initializesystem(sys) @@ -78,7 +78,7 @@ end @mtkcompile sys = System([D(x) ~ y[1] + y[2], y ~ foo(x)], t) @test length(equations(sys)) == 1 @test length(observed(sys)) == 4 - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn2]) + prob = ODEProblem(sys, [x => 1.0, foo => _tmp_fn2], (0.0, 1.0)) val[] = 0 @test_nowarn prob.f(prob.u0, prob.p, 0.0) @test val[] == 1 @@ -99,7 +99,7 @@ end @test length(equations(sys)) == 5 @test length(observed(sys)) == 2 prob = ODEProblem( - sys, [y => ones(2), z => 2ones(2), x => 3.0], (0.0, 1.0), [foo => _tmp_fn2]) + sys, [y => ones(2), z => 2ones(2), x => 3.0, foo => _tmp_fn2], (0.0, 1.0)) val[] = 0 @test_nowarn prob.f(prob.u0, prob.p, 0.0) @test val[] == 2 @@ -397,7 +397,7 @@ end defaults = [p => missing], guesses = [p => 1.0, y => 1.0]) @test length(equations(sys)) == 2 @test length(parameters(sys)) == 2 - prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [q => 2.0]) + prob = ODEProblem(sys, [x => 1.0, q => 2.0], (0.0, 1.0)) integ = init(prob, Rodas5P(); abstol = 1e-10, reltol = 1e-8) @test integ.ps[p]≈1.0 atol=1e-6 @test integ[y]≈0.0 atol=1e-5 @@ -413,7 +413,7 @@ end @test length(observed(sys)) == 1 @test observed(sys)[1].lhs in Set([x, y]) @test length(parameters(sys)) == 2 - prob = NonlinearProblem(sys, [x => 1.0, y => 1.0], [q => 1.0]) + prob = NonlinearProblem(sys, [x => 1.0, y => 1.0, q => 1.0]) integ = init(prob, NewtonRaphson()) @test prob.ps[p] ≈ 2.0 end @@ -431,13 +431,13 @@ end @test isempty(brownians(sys)) neqs = ModelingToolkit.get_noise_eqs(sys) @test issetequal(sum.(eachrow(neqs)), [q, 1 + p]) - prob = SDEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0), [q => 1.0]) + prob = SDEProblem(sys, [x => 1.0, y => 1.0, q => 1.0], (0.0, 1.0)) integ = init(prob, ImplicitEM()) @test integ.ps[p] ≈ 3.0 end end -@testset "Deprecated `structural_simplify` and `@mtkbuild`" begin +@testset "Deprecated `mtkcompile` and `@mtkcompile`" begin @variables x(t) @test_deprecated @mtkbuild sys = System([D(x) ~ x], t) @named sys = System([D(x) ~ x], t) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 94a443fc4d..e61ea9779c 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -454,7 +454,7 @@ end function testsol( sys, probtype, solver, u0, p, tspan; tstops = Float64[], paramtotest = nothing, kwargs...) - prob = probtype(complete(sys), u0, tspan, p; kwargs...) + prob = probtype(complete(sys), [u0; p], tspan; kwargs...) sol = solve(prob, solver(); tstops = tstops, abstol = 1e-10, reltol = 1e-10) @test isapprox(sol(1.0000000001)[1] - sol(0.999999999)[1], 1.0; rtol = 1e-6) paramtotest === nothing || (@test sol.ps[paramtotest] == [0.0, 1.0]) @@ -521,7 +521,7 @@ end @named osys4 = System(eqs, t, [A], [k, t1], discrete_events = [cb1, cb2‵‵]) @named ssys4 = SDESystem(eqs, [0.0], t, [A], [k, t1], discrete_events = [cb1, cb2‵‵]) - oprob4 = ODEProblem(complete(osys4), u0, tspan, p) + oprob4 = ODEProblem(complete(osys4), [u0; p], tspan) testsol(osys4, ODEProblem, Tsit5, u0, p, tspan; tstops = [1.0], paramtotest = k) testsol(ssys4, SDEProblem, RI5, u0, p, tspan; tstops = [1.0], paramtotest = k) @@ -558,7 +558,7 @@ end function testsol(jsys, u0, p, tspan; tstops = Float64[], paramtotest = nothing, N = 40000, kwargs...) jsys = complete(jsys) - jprob = JumpProblem(jsys, u0, tspan, p; aggregator = Direct(), kwargs...) + jprob = JumpProblem(jsys, [u0; p], tspan; aggregator = Direct(), kwargs...) sol = solve(jprob, SSAStepper(); tstops = tstops) @test (sol(1.000000000001)[1] - sol(0.99999999999)[1]) == 1 paramtotest === nothing || (@test sol.ps[paramtotest] == [0.0, 1.0]) @@ -634,7 +634,7 @@ end @named oneosc_ce = compose(eqs_sys, oscce) oneosc_ce_simpl = mtkcompile(oneosc_ce) - prob = ODEProblem(oneosc_ce_simpl, [], (0.0, 2.0), []) + prob = ODEProblem(oneosc_ce_simpl, [], (0.0, 2.0)) sol = solve(prob, Tsit5(), saveat = 0.1) @test typeof(oneosc_ce_simpl) == System @@ -887,7 +887,7 @@ end @mtkcompile sys = System(D(x) ~ cos(t), t, [x], [a, b, c]; continuous_events = [cb1, cb2], discrete_events = [cb3]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 2pi), [a => 1.0, b => 2.0, c => 0.0]) + prob = ODEProblem(sys, [x => 1.0, a => 1.0, b => 2.0, c => 0.0], (0.0, 2pi)) @test sort(canonicalize(Discrete(), prob.p)[1]) == [0.0, 1.0, 2.0] sol = solve(prob, Tsit5()) @@ -1071,7 +1071,7 @@ end cb1 = ModelingToolkit.SymbolicContinuousCallback( [x ~ 0], nothing, initialize = [x ~ 1.5], finalize = f) @mtkcompile sys = System(D(x) ~ -1, t, [x], []; continuous_events = [cb1]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) + prob = ODEProblem(sys, [x => 1.0], (0.0, 2)) sol = solve(prob, Tsit5(); dtmax = 0.01) @test sol[x][1] ≈ 1.0 @test sol[x][2] ≈ 1.5 # the initialize affect has been applied @@ -1089,7 +1089,7 @@ end cb2 = ModelingToolkit.SymbolicContinuousCallback( [x ~ 0.1], nothing, initialize = a, finalize = b) @mtkcompile sys = System(D(x) ~ -1, t, [x], []; continuous_events = [cb1, cb2]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) + prob = ODEProblem(sys, [x => 1.0], (0.0, 2)) sol = solve(prob, Tsit5()) @test sol[x][1] ≈ 1.0 @test sol[x][2] ≈ 1.5 # the initialize affect has been applied @@ -1103,7 +1103,7 @@ end cb3 = ModelingToolkit.SymbolicDiscreteCallback( 1.0, [x ~ 2], initialize = a, finalize = b) @mtkcompile sys = System(D(x) ~ -1, t, [x], []; discrete_events = [cb3]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) + prob = ODEProblem(sys, [x => 1.0], (0.0, 2)) sol = solve(prob, Tsit5()) @test inited == true @test finaled == true @@ -1116,7 +1116,7 @@ end finaled = false cb3 = ModelingToolkit.SymbolicDiscreteCallback(1.0, f, initialize = a, finalize = b) @mtkcompile sys = System(D(x) ~ -1, t, [x], []; discrete_events = [cb3]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) + prob = ODEProblem(sys, [x => 1.0], (0.0, 2)) sol = solve(prob, Tsit5()) @test seen == true @test inited == true @@ -1127,7 +1127,7 @@ end finaled = false cb3 = ModelingToolkit.SymbolicDiscreteCallback([1.0], f, initialize = a, finalize = b) @mtkcompile sys = System(D(x) ~ -1, t, [x], []; discrete_events = [cb3]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) + prob = ODEProblem(sys, [x => 1.0], (0.0, 2)) sol = solve(prob, Tsit5()) @test seen == true @test inited == true @@ -1140,7 +1140,7 @@ end cb3 = ModelingToolkit.SymbolicDiscreteCallback( t == 1.0, f, initialize = a, finalize = b) @mtkcompile sys = System(D(x) ~ -1, t, [x], []; discrete_events = [cb3]) - prob = ODEProblem(sys, [x => 1.0], (0.0, 2), []) + prob = ODEProblem(sys, [x => 1.0], (0.0, 2)) sol = solve(prob, Tsit5(); tstops = 1.0) @test seen == true @test inited == true @@ -1183,7 +1183,7 @@ end end end @mtkcompile decay = DECAY() - prob = ODEProblem(decay, [], (0.0, 10.0), []) + prob = ODEProblem(decay, [], (0.0, 10.0)) @test_nowarn solve(prob, Tsit5(), tstops = [1.0]) end @@ -1253,7 +1253,7 @@ end x^2 + y^2 ~ 1] c_evt = [t ~ 5.0] => [x ~ Pre(x) + 0.1] @mtkcompile pend = System(eqs, t, continuous_events = c_evt) - prob = ODEProblem(pend, [x => -1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + prob = ODEProblem(pend, [x => -1, y => 0, g => 1], (0.0, 10.0), guesses = [λ => 1]) sol = solve(prob, FBDF()) @test ≈(sol(5.000001, idxs = x) - sol(4.999999, idxs = x), 0.1, rtol = 1e-4) @test ≈(sol(5.000001, idxs = x)^2 + sol(5.000001, idxs = y)^2, 1, rtol = 1e-4) @@ -1261,7 +1261,7 @@ end # Implicit affect with Pre c_evt = [t ~ 5.0] => [x ~ Pre(x) + y^2] @mtkcompile pend = System(eqs, t, continuous_events = c_evt) - prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 10.0), guesses = [λ => 1]) sol = solve(prob, FBDF()) @test ≈(sol(5.000001, idxs = y)^2 + sol(4.999999, idxs = x), sol(5.000001, idxs = x), rtol = 1e-4) @@ -1270,7 +1270,7 @@ end # Impossible affect errors c_evt = [t ~ 5.0] => [x ~ Pre(x) + 2] @mtkcompile pend = System(eqs, t, continuous_events = c_evt) - prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 10.0), guesses = [λ => 1]) @test_throws UnsolvableCallbackError sol=solve(prob, FBDF()) # Changing both variables and parameters in the same affect. @@ -1281,7 +1281,7 @@ end c_evt = SymbolicContinuousCallback( [t ~ 5.0], [x ~ Pre(x) + 0.1, g ~ Pre(g) + 1], discrete_parameters = [g], iv = t) @mtkcompile pend = System(eqs, t, continuous_events = c_evt) - prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 10.0), [g => 1], guesses = [λ => 1]) + prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 10.0), guesses = [λ => 1]) sol = solve(prob, FBDF()) @test sol.ps[g] ≈ [1, 2] @test ≈(sol(5.0000001, idxs = x) - sol(4.999999, idxs = x), 0.1, rtol = 1e-4) @@ -1291,7 +1291,7 @@ end c_evt = SymbolicContinuousCallback( [t ~ 5.0], [x ~ Pre(x) + 1, g ~ Pre(g) + 1], discrete_parameters = [g], iv = t) @mtkcompile sys = System(eqs, t, continuous_events = c_evt) - prob = ODEProblem(sys, [x => 1.0], (0.0, 10.0), [g => 2]) + prob = ODEProblem(sys, [x => 1.0, g => 2], (0.0, 10.0)) sol = solve(prob, FBDF()) @test sol.ps[g] ≈ [2.0, 3.0] @test ≈(sol(5.00000001, idxs = x) - sol(4.9999999, idxs = x), 1; rtol = 1e-4) @@ -1300,7 +1300,7 @@ end # Parameters that don't appear in affects should not be mutated. c_evt = [t ~ 5.0] => [x ~ Pre(x) + 1] @mtkcompile sys = System(eqs, t, continuous_events = c_evt) - prob = ODEProblem(sys, [x => 0.5], (0.0, 10.0), [g => 2], guesses = [y => 0]) + prob = ODEProblem(sys, [x => 0.5, g => 2], (0.0, 10.0), guesses = [y => 0]) sol = solve(prob, FBDF()) @test prob.ps[g] == sol.ps[g] end diff --git a/test/symbolic_indexing_interface.jl b/test/symbolic_indexing_interface.jl index f71fe93842..81d444668e 100644 --- a/test/symbolic_indexing_interface.jl +++ b/test/symbolic_indexing_interface.jl @@ -37,7 +37,7 @@ using SciMLStructures: Tunable @test default_values(odesys)[y] == 2.0 @test isequal(default_values(odesys)[xy], x + y) - prob = ODEProblem(odesys, [], (0.0, 1.0), [a => 1.0, b => 2.0]) + prob = ODEProblem(odesys, [a => 1.0, b => 2.0], (0.0, 1.0)) getter = getu(odesys, (x + 1, x + 2)) @test getter(prob) isa Tuple @test_nowarn @inferred getter(prob) @@ -122,7 +122,7 @@ end @test pobs(ps) == [3.0, 5.0] prob = NonlinearProblem( - ns, [x => 1.0, y => 2.0, z => 3.0], [σ => 1.0, ρ => 2.0, β => 3.0]) + ns, [x => 1.0, y => 2.0, z => 3.0, σ => 1.0, ρ => 2.0, β => 3.0]) getter = getu(ns, (x + 1, x + 2)) @test getter(prob) isa Tuple @test_nowarn @inferred getter(prob) @@ -197,7 +197,7 @@ end @parameters p1 p2[1:2, 1:2] @mtkcompile sys = System([D(x) ~ x * t + p1, y ~ 2x, D(z) ~ p2 * z], t) prob = ODEProblem( - sys, [x => 1.0, z => ones(2)], (0.0, 1.0), [p1 => 2.0, p2 => ones(2, 2)]) + sys, [x => 1.0, z => ones(2), p1 => 2.0, p2 => ones(2, 2)], (0.0, 1.0)) @test getu(prob, x)(prob) == getu(prob, :x)(prob) @test getu(prob, [x, y])(prob) == getu(prob, [:x, :y])(prob) @test getu(prob, z)(prob) == getu(prob, :z)(prob) diff --git a/test/symbolic_parameters.jl b/test/symbolic_parameters.jl index 15b19b320a..7a2fef3500 100644 --- a/test/symbolic_parameters.jl +++ b/test/symbolic_parameters.jl @@ -26,7 +26,7 @@ resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns), defaults = ModelingToolkit.defaults(ns)) @test resolved == [1, 0.1 + 1, (0.1 + 1) * 1.1] -prob = NonlinearProblem(complete(ns), [u => 1.0], Pair[]) +prob = NonlinearProblem(complete(ns), [u => 1.0]) @test prob.u0 == [1.0, 1.1, 0.9] sol = solve(prob, NewtonRaphson()) @@ -41,7 +41,7 @@ res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top), @test res == [0.5, 1, 0.1 + 1, (0.1 + 1) * 1.1] top = complete(top) -prob = NonlinearProblem(top, [unknowns(ns, u) => 1.0, a => 1.0], []) +prob = NonlinearProblem(top, [unknowns(ns, u) => 1.0, a => 1.0]) @test prob.u0 == [1.0, 0.5, 1.1, 0.9] sol = solve(prob, NewtonRaphson()) @@ -61,12 +61,8 @@ der = Differential(t) eqs = [der(x) ~ x] @named sys = System(eqs, t, vars, [x0]) sys = complete(sys) -pars = [ - x0 => 10.0 -] -initialValues = [ - x => x0 -] +initialValues = [x => x0 + x0 => 10.0] tspan = (0.0, 1.0) -problem = ODEProblem(sys, initialValues, tspan, pars) +problem = ODEProblem(sys, initialValues, tspan) @test problem.u0 isa Vector{Float64}