diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index dceee5f1d2..14bea532b3 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -31,8 +31,8 @@ jobs: JULIA_DEBUG: "Documenter" run: DISPLAY=:0 xvfb-run -s '-screen 0 1024x768x24' julia --project=docs/ --code-coverage=user docs/make.jl - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v5 with: - file: lcov.info + files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} fail_ci_if_error: true diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 7a7556efa3..f0aefa87f6 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -68,8 +68,8 @@ jobs: exit(0) # Exit immediately, as a success end - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v5 with: - file: lcov.info + files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} - fail_ci_if_error: true + fail_ci_if_error: false diff --git a/.github/workflows/Tests.yml b/.github/workflows/Tests.yml index 609ab2fb3d..a95f19a085 100644 --- a/.github/workflows/Tests.yml +++ b/.github/workflows/Tests.yml @@ -32,6 +32,7 @@ jobs: group: - InterfaceI - InterfaceII + - Initialization - SymbolicIndexingInterface - Extended - Extensions diff --git a/Project.toml b/Project.toml index 0e804e1f63..1621669c2c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.50.0" +version = "9.54.0" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" @@ -64,6 +64,7 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" +InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" [extensions] @@ -72,6 +73,7 @@ MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKHomotopyContinuationExt = "HomotopyContinuation" MTKLabelledArraysExt = "LabelledArrays" +MTKInfiniteOptExt = "InfiniteOpt" [compat] AbstractTrees = "0.3, 0.4" @@ -104,6 +106,7 @@ FunctionWrappers = "1.1" FunctionWrappersWrappers = "0.1" Graphs = "1.5.2" HomotopyContinuation = "2.11" +InfiniteOpt = "0.5" InteractiveUtils = "1" JuliaFormatter = "1.0.47" JumpProcesses = "9.13.1" @@ -123,7 +126,7 @@ REPL = "1" RecursiveArrayTools = "3.26" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "2.57.1" +SciMLBase = "2.64" SciMLStructures = "1.0" Serialization = "1" Setfield = "0.7, 0.8, 1" @@ -131,9 +134,9 @@ SimpleNonlinearSolve = "0.1.0, 1, 2" SparseArrays = "1" SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" StaticArrays = "0.10, 0.11, 0.12, 1.0" -SymbolicIndexingInterface = "0.3.31" +SymbolicIndexingInterface = "0.3.35" SymbolicUtils = "3.7" -Symbolics = "6.15.4" +Symbolics = "6.19" URIs = "1" UnPack = "0.1, 1.0" Unitful = "1.1" diff --git a/docs/Project.toml b/docs/Project.toml index 781508f9d2..24f00b4db2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -33,8 +33,8 @@ DynamicQuantities = "^0.11.2, 0.12, 1" ModelingToolkit = "8.33, 9" NonlinearSolve = "3, 4" Optim = "1.7" -Optimization = "3.9" -OptimizationOptimJL = "0.1" +Optimization = "3.9, 4" +OptimizationOptimJL = "0.1, 0.4" OrdinaryDiffEq = "6.31" Plots = "1.36" SciMLStructures = "1.1" diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index bb3d01fc30..10671299c6 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -82,6 +82,16 @@ parameter_index(sys, sym) Note that while the variable index will be an integer, the parameter index is a struct of type `ParameterIndex` whose internals should not be relied upon. +## Can I index with strings? + +Strings are not considered symbolic variables, and thus cannot directly be used for symbolic +indexing. However, ModelingToolkit does provide a method to parse the string representation of +a variable, given the system in which that variable exists. + +```@docs +ModelingToolkit.parse_variable +``` + ## Transforming value maps to arrays ModelingToolkit.jl allows (and recommends) input maps like `[x => 2.0, y => 3.0]` diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/basics/Variable_metadata.md index 6462c3a469..e9bdabfe38 100644 --- a/docs/src/basics/Variable_metadata.md +++ b/docs/src/basics/Variable_metadata.md @@ -86,6 +86,44 @@ hasbounds(u) getbounds(u) ``` +Bounds can also be specified for array variables. A scalar array bound is applied to each +element of the array. A bound may also be specified as an array, in which case the size of +the array must match the size of the symbolic variable. + +```@example metadata +@variables x[1:2, 1:2] [bounds = (-1, 1)] +hasbounds(x) +``` + +```@example metadata +getbounds(x) +``` + +```@example metadata +getbounds(x[1, 1]) +``` + +```@example metadata +getbounds(x[1:2, 1]) +``` + +```@example metadata +@variables x[1:2] [bounds = (-Inf, [1.0, Inf])] +hasbounds(x) +``` + +```@example metadata +getbounds(x) +``` + +```@example metadata +getbounds(x[2]) +``` + +```@example metadata +hasbounds(x[2]) +``` + ## Guess Specify an initial guess for custom initial conditions of an `ODESystem`. diff --git a/docs/src/examples/higher_order.md b/docs/src/examples/higher_order.md index 7dafe758dc..fac707525f 100644 --- a/docs/src/examples/higher_order.md +++ b/docs/src/examples/higher_order.md @@ -3,7 +3,7 @@ ModelingToolkit has a system for transformations of mathematical systems. These transformations allow for symbolically changing the representation of the model to problems that are easier to -numerically solve. One simple to demonstrate transformation is the +numerically solve. One simple to demonstrate transformation, is `structural_simplify`, which does a lot of tricks, one being the transformation that turns an Nth order ODE into N coupled 1st order ODEs. @@ -15,16 +15,28 @@ We utilize the derivative operator twice here to define the second order: using ModelingToolkit, OrdinaryDiffEq using ModelingToolkit: t_nounits as t, D_nounits as D -@parameters σ ρ β -@variables x(t) y(t) z(t) - -eqs = [D(D(x)) ~ σ * (y - x), - D(y) ~ x * (ρ - z) - y, - D(z) ~ x * y - β * z] - -@named sys = ODESystem(eqs, t) +@mtkmodel SECOND_ORDER begin + @parameters begin + σ = 28.0 + ρ = 10.0 + β = 8 / 3 + end + @variables begin + x(t) = 1.0 + y(t) = 0.0 + z(t) = 0.0 + end + @equations begin + D(D(x)) ~ σ * (y - x) + D(y) ~ x * (ρ - z) - y + D(z) ~ x * y - β * z + end +end +@mtkbuild sys = SECOND_ORDER() ``` +The second order ODE has been automatically transformed to two first order ODEs. + Note that we could've used an alternative syntax for 2nd order, i.e. `D = Differential(t)^2` and then `D(x)` would be the second derivative, and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compose @@ -33,28 +45,17 @@ and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compos Now let's transform this into the `ODESystem` of first order components. We do this by calling `structural_simplify`: -```@example orderlowering -sys = structural_simplify(sys) -``` - Now we can directly numerically solve the lowered system. Note that, following the original problem, the solution requires knowing the -initial condition for `x'`, and thus we include that in our input -specification: +initial condition for both `x` and `D(x)`. +The former already got assigned a default value in the `@mtkmodel`, +but we still have to provide a value for the latter. ```@example orderlowering -u0 = [D(x) => 2.0, - x => 1.0, - y => 0.0, - z => 0.0] - -p = [σ => 28.0, - ρ => 10.0, - β => 8 / 3] - +u0 = [D(sys.x) => 2.0] tspan = (0.0, 100.0) -prob = ODEProblem(sys, u0, tspan, p, jac = true) +prob = ODEProblem(sys, u0, tspan, [], jac = true) sol = solve(prob, Tsit5()) -using Plots; -plot(sol, idxs = (x, y)); +using Plots +plot(sol, idxs = (sys.x, sys.y)) ``` diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md index cd1e843f55..407248709d 100644 --- a/docs/src/examples/perturbation.md +++ b/docs/src/examples/perturbation.md @@ -39,7 +39,8 @@ eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2) !!! note The 0-th order equation can be solved analytically, but ModelingToolkit does currently not feature automatic analytical solution of ODEs, so we proceed with solving it numerically. - These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities: + +These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities: ```@example perturbation @mtkbuild sys = ODESystem(eqs_pert, t) diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index c3dbf8bbb0..c01310eb06 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -113,7 +113,7 @@ the `initialization_eqs` keyword argument, for example: ```@example init prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1], - initialization_eqs = [y ~ 1]) + initialization_eqs = [y ~ 0]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) ``` diff --git a/ext/MTKHomotopyContinuationExt.jl b/ext/MTKHomotopyContinuationExt.jl index 19c21483c3..fa5d1bfcd4 100644 --- a/ext/MTKHomotopyContinuationExt.jl +++ b/ext/MTKHomotopyContinuationExt.jl @@ -54,16 +54,73 @@ end PolynomialData() = PolynomialData(BasicSymbolic[], NonPolynomialReason.T[], false) +abstract type PolynomialTransformationError <: Exception end + +struct MultivarTerm <: PolynomialTransformationError + term::Any + vars::Any +end + +function Base.showerror(io::IO, err::MultivarTerm) + println(io, + "Cannot convert system to polynomial: Found term $(err.term) which is a function of multiple unknowns $(err.vars).") +end + +struct MultipleTermsOfSameVar <: PolynomialTransformationError + terms::Any + var::Any +end + +function Base.showerror(io::IO, err::MultipleTermsOfSameVar) + println(io, + "Cannot convert system to polynomial: Found multiple non-polynomial terms $(err.terms) involving the same unknown $(err.var).") +end + +struct SymbolicSolveFailure <: PolynomialTransformationError + term::Any + var::Any +end + +function Base.showerror(io::IO, err::SymbolicSolveFailure) + println(io, + "Cannot convert system to polynomial: Unable to symbolically solve $(err.term) for $(err.var).") +end + +struct NemoNotLoaded <: PolynomialTransformationError end + +function Base.showerror(io::IO, err::NemoNotLoaded) + println(io, + "ModelingToolkit may be able to solve this system as a polynomial system if `Nemo` is loaded. Run `import Nemo` and try again.") +end + +struct VariablesAsPolyAndNonPoly <: PolynomialTransformationError + vars::Any +end + +function Base.showerror(io::IO, err::VariablesAsPolyAndNonPoly) + println(io, + "Cannot convert convert system to polynomial: Variables $(err.vars) occur in both polynomial and non-polynomial terms in the system.") +end + struct NotPolynomialError <: Exception - eq::Equation - data::PolynomialData + transformation_err::Union{PolynomialTransformationError, Nothing} + eq::Vector{Equation} + data::Vector{PolynomialData} end function Base.showerror(io::IO, err::NotPolynomialError) - println(io, - "Equation $(err.eq) is not a polynomial in the unknowns for the following reasons:") - for (term, reason) in zip(err.data.non_polynomial_terms, err.data.reasons) - println(io, display_reason(reason, term)) + if err.transformation_err !== nothing + Base.showerror(io, err.transformation_err) + end + for (eq, data) in zip(err.eq, err.data) + if isempty(data.non_polynomial_terms) + continue + end + println(io, + "Equation $(eq) is not a polynomial in the unknowns for the following reasons:") + for (term, reason) in zip(data.non_polynomial_terms, data.reasons) + println(io, display_reason(reason, term)) + end end end @@ -86,7 +143,9 @@ function process_polynomial!(data::PolynomialData, x, wrt) any(isequal(x), wrt) && return true if operation(x) in (*, +, -, /) - return all(y -> is_polynomial!(data, y, wrt), arguments(x)) + # `map` because `all` will early exit, but we want to search + # through everything to get all the non-polynomial terms + return all(map(y -> is_polynomial!(data, y, wrt), arguments(x))) end if operation(x) == (^) b, p = arguments(x) @@ -105,10 +164,6 @@ function process_polynomial!(data::PolynomialData, x, wrt) push!(data.reasons, NonPolynomialReason.ExponentContainsUnknowns) end base_polynomial = is_polynomial!(data, b, wrt) - if !base_polynomial - push!(data.non_polynomial_terms, x) - push!(data.reasons, NonPolynomialReason.BaseNotPolynomial) - end return base_polynomial && !exponent_has_unknowns && is_pow_integer end push!(data.non_polynomial_terms, x) @@ -234,6 +289,12 @@ end SymbolicIndexingInterface.parameter_values(s::MTKHomotopySystem) = s.p +struct PolynomialTransformationData + new_var::BasicSymbolic + term::BasicSymbolic + inv_term::Vector +end + """ $(TYPEDSIGNATURES) @@ -265,18 +326,95 @@ function MTK.HomotopyContinuationProblem( # CSE/hashconsing. eqs = full_equations(sys) - denoms = [] - has_parametric_exponents = false - eqs2 = map(eqs) do eq + polydata = map(eqs) do eq data = PolynomialData() process_polynomial!(data, eq.lhs, dvs) process_polynomial!(data, eq.rhs, dvs) - has_parametric_exponents |= data.has_parametric_exponent - if !isempty(data.non_polynomial_terms) - throw(NotPolynomialError(eq, data)) + data + end + + has_parametric_exponents = any(d -> d.has_parametric_exponent, polydata) + + all_non_poly_terms = mapreduce(d -> d.non_polynomial_terms, vcat, polydata) + unique!(all_non_poly_terms) + + var_to_nonpoly = Dict{BasicSymbolic, PolynomialTransformationData}() + + is_poly = true + transformation_err = nothing + for t in all_non_poly_terms + # if the term involves multiple unknowns, we can't invert it + dvs_in_term = map(x -> occursin(x, t), dvs) + if count(dvs_in_term) > 1 + transformation_err = MultivarTerm(t, dvs[dvs_in_term]) + is_poly = false + break + end + # we already have a substitution solving for `var` + var = dvs[findfirst(dvs_in_term)] + if haskey(var_to_nonpoly, var) && !isequal(var_to_nonpoly[var].term, t) + transformation_err = MultipleTermsOfSameVar([t, var_to_nonpoly[var].term], var) + is_poly = false + break + end + # we want to solve `term - new_var` for `var` + new_var = gensym(Symbol(var)) + new_var = unwrap(only(@variables $new_var)) + invterm = Symbolics.ia_solve( + t - new_var, var; complex_roots = false, periodic_roots = false, warns = false) + # if we can't invert it, quit + if invterm === nothing || isempty(invterm) + transformation_err = SymbolicSolveFailure(t, var) + is_poly = false + break + end + # `ia_solve` returns lazy terms i.e. `asin(1.0)` instead of `pi/2` + # this just evaluates the constant expressions + invterm = Symbolics.substitute.(invterm, (Dict(),)) + # RootsOf implies Symbolics couldn't solve the inner polynomial because + # `Nemo` wasn't loaded. + if any(x -> MTK.iscall(x) && MTK.operation(x) == Symbolics.RootsOf, invterm) + transformation_err = NemoNotLoaded() + is_poly = false + break end - num, den = handle_rational_polynomials(eq.rhs - eq.lhs, dvs) + var_to_nonpoly[var] = PolynomialTransformationData(new_var, t, invterm) + end + + if !is_poly + throw(NotPolynomialError(transformation_err, eqs, polydata)) + end + + subrules = Dict() + combinations = Vector[] + new_dvs = [] + for x in dvs + if haskey(var_to_nonpoly, x) + _data = var_to_nonpoly[x] + subrules[_data.term] = _data.new_var + push!(combinations, _data.inv_term) + push!(new_dvs, _data.new_var) + else + push!(combinations, [x]) + push!(new_dvs, x) + end + end + all_solutions = collect.(collect(Iterators.product(combinations...))) + denoms = [] + eqs2 = map(eqs) do eq + t = eq.rhs - eq.lhs + t = Symbolics.fixpoint_sub(t, subrules; maxiters = length(dvs)) + # the substituted variable occurs outside the substituted term + poly_and_nonpoly = map(dvs) do x + haskey(var_to_nonpoly, x) && occursin(x, t) + end + if any(poly_and_nonpoly) + throw(NotPolynomialError( + VariablesAsPolyAndNonPoly(dvs[poly_and_nonpoly]), eqs, polydata)) + end + + num, den = handle_rational_polynomials(t, new_dvs) # make factors different elements, otherwise the nonzero factors artificially # inflate the error of the zero factor. if iscall(den) && operation(den) == * @@ -292,16 +430,19 @@ function MTK.HomotopyContinuationProblem( end sys2 = MTK.@set sys.eqs = eqs2 + MTK.@set! sys2.unknowns = new_dvs # remove observed equations to avoid adding them in codegen MTK.@set! sys2.observed = Equation[] MTK.@set! sys2.substitutions = nothing - nlfn, u0, p = MTK.process_SciMLProblem(NonlinearFunction{true}, sys2, u0map, parammap; - jac = true, eval_expression, eval_module) + _, u0, p = MTK.process_SciMLProblem( + MTK.EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module) + nlfn = NonlinearFunction{true}(sys2; jac = true, eval_expression, eval_module) - denominator = MTK.build_explicit_observed_function(sys, denoms) + denominator = MTK.build_explicit_observed_function(sys2, denoms) + unpack_solution = MTK.build_explicit_observed_function(sys2, all_solutions) - hvars = symbolics_to_hc.(dvs) + hvars = symbolics_to_hc.(new_dvs) mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(eqs)) obsfn = MTK.ObservedFunctionCache(sys; eval_expression, eval_module) @@ -319,7 +460,7 @@ function MTK.HomotopyContinuationProblem( solver_and_starts = HomotopyContinuation.solver_startsolutions(mtkhsys; kwargs...) end return MTK.HomotopyContinuationProblem( - u0, mtkhsys, denominator, sys, obsfn, solver_and_starts) + u0, mtkhsys, denominator, sys, obsfn, solver_and_starts, unpack_solution) end """ @@ -353,25 +494,35 @@ function CommonSolve.solve(prob::MTK.HomotopyContinuationProblem, if isempty(realsols) u = state_values(prob) retcode = SciMLBase.ReturnCode.ConvergenceFailure + resid = prob.homotopy_continuation_system(u) else T = eltype(state_values(prob)) - distance, idx = findmin(realsols) do result + distance = T(Inf) + u = state_values(prob) + resid = nothing + for result in realsols if any(<=(denominator_abstol), prob.denominator(real.(result.solution), parameter_values(prob))) - return T(Inf) + continue + end + for truesol in prob.unpack_solution(result.solution, parameter_values(prob)) + dist = norm(truesol - state_values(prob)) + if dist < distance + distance = dist + u = T.(real.(truesol)) + resid = T.(real.(prob.homotopy_continuation_system(result.solution))) + end end - norm(result.solution - state_values(prob)) end # all roots cause denominator to be zero if isinf(distance) u = state_values(prob) + resid = prob.homotopy_continuation_system(u) retcode = SciMLBase.ReturnCode.Infeasible else - u = real.(realsols[idx].solution) retcode = SciMLBase.ReturnCode.Success end end - resid = prob.homotopy_continuation_system(u) return SciMLBase.build_solution( prob, :HomotopyContinuation, u, resid; retcode, original = sol) diff --git a/ext/MTKInfiniteOptExt.jl b/ext/MTKInfiniteOptExt.jl new file mode 100644 index 0000000000..3b482a1f03 --- /dev/null +++ b/ext/MTKInfiniteOptExt.jl @@ -0,0 +1,26 @@ +module MTKInfiniteOptExt +import ModelingToolkit +import SymbolicUtils +import NaNMath +import InfiniteOpt +import InfiniteOpt: JuMP, GeneralVariableRef + +# This file contains method definitions to make it possible to trace through functions generated by MTK using JuMP variables + +for ff in [acos, log1p, acosh, log2, asin, tan, atanh, cos, log, sin, log10, sqrt] + f = nameof(ff) + # These need to be defined so that JuMP can trace through functions built by Symbolics + @eval NaNMath.$f(x::GeneralVariableRef) = Base.$f(x) +end + +# JuMP variables and Symbolics variables never compare equal. When tracing through dynamics, a function argument can be either a JuMP variable or A Symbolics variable, it can never be both. +function Base.isequal(::SymbolicUtils.Symbolic, + ::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr, InfiniteOpt.AbstractInfOptExpr}) + false +end +function Base.isequal( + ::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr, InfiniteOpt.AbstractInfOptExpr}, + ::SymbolicUtils.Symbolic) + false +end +end diff --git a/src/inputoutput.jl b/src/inputoutput.jl index 9ed0984050..4a99ec11f5 100644 --- a/src/inputoutput.jl +++ b/src/inputoutput.jl @@ -235,18 +235,20 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu # TODO: add an optional check on the ordering of observed equations u = map(x -> time_varying_as_func(value(x), sys), dvs) p = map(x -> time_varying_as_func(value(x), sys), ps) + p = reorder_parameters(sys, p) t = get_iv(sys) # pre = has_difference ? (ex -> ex) : get_postprocess_fbody(sys) - args = (u, inputs, p, t) + args = (u, inputs, p..., t) if implicit_dae ddvs = map(Differential(get_iv(sys)), dvs) args = (ddvs, args...) end process = get_postprocess_fbody(sys) f = build_function(rhss, args...; postprocess_fbody = process, - expression = Val{true}, wrap_code = wrap_array_vars(sys, rhss; dvs, ps) .∘ + expression = Val{true}, wrap_code = wrap_mtkparameters(sys, false, 3) .∘ + wrap_array_vars(sys, rhss; dvs, ps, inputs) .∘ wrap_parameter_dependencies(sys, false), kwargs...) f = eval_or_rgf.(f; eval_expression, eval_module) @@ -426,7 +428,7 @@ function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kw augmented_sys = ODESystem(eqs, t, systems = [dsys], name = gensym(:outer)) augmented_sys = extend(augmented_sys, sys) - (f_oop, f_ip), dvs, p = generate_control_function(augmented_sys, all_inputs, + (f_oop, f_ip), dvs, p, io_sys = generate_control_function(augmented_sys, all_inputs, [d]; kwargs...) - (f_oop, f_ip), augmented_sys, dvs, p + (f_oop, f_ip), augmented_sys, dvs, p, io_sys end diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index a854acb9b1..8161a9572b 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -88,6 +88,14 @@ function tearing_sub(expr, dict, s) s ? simplify(expr) : expr end +function tearing_substitute_expr(sys::AbstractSystem, expr; simplify = false) + empty_substitutions(sys) && return expr + substitutions = get_substitutions(sys) + @unpack subs = substitutions + solved = Dict(eq.lhs => eq.rhs for eq in subs) + return tearing_sub(expr, solved, simplify) +end + """ $(TYPEDSIGNATURES) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl index 8b2ac1f69c..bd24a1d017 100644 --- a/src/structural_transformation/utils.jl +++ b/src/structural_transformation/utils.jl @@ -58,7 +58,41 @@ end ### ### Structural check ### -function check_consistency(state::TransformationState, orig_inputs) + +""" + $(TYPEDSIGNATURES) + +Check if the `state` represents a singular system, and return the unmatched variables. +""" +function singular_check(state::TransformationState) + @unpack graph, var_to_diff = state.structure + fullvars = get_fullvars(state) + # This is defined to check if Pantelides algorithm terminates. For more + # details, check the equation (15) of the original paper. + extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist; + map(collect, edges(var_to_diff))]) + extended_var_eq_matching = maximal_matching(extended_graph) + + nvars = ndsts(graph) + unassigned_var = [] + for (vj, eq) in enumerate(extended_var_eq_matching) + vj > nvars && break + if eq === unassigned && !isempty(𝑑neighbors(graph, vj)) + push!(unassigned_var, fullvars[vj]) + end + end + return unassigned_var +end + +""" + $(TYPEDSIGNATURES) + +Check the consistency of `state`, given the inputs `orig_inputs`. If `nothrow == false`, +throws an error if the system is under-/over-determined or singular. In this case, if the +function returns it will return `true`. If `nothrow == true`, it will return `false` +instead of throwing an error. The singular case will print a warning. +""" +function check_consistency(state::TransformationState, orig_inputs; nothrow = false) fullvars = get_fullvars(state) neqs = n_concrete_eqs(state) @unpack graph, var_to_diff = state.structure @@ -72,6 +106,7 @@ function check_consistency(state::TransformationState, orig_inputs) is_balanced = n_highest_vars == neqs if neqs > 0 && !is_balanced + nothrow && return false varwhitelist = var_to_diff .== nothing var_eq_matching = maximal_matching(graph, eq -> true, v -> varwhitelist[v]) # not assigned # Just use `error_reporting` to do conditional @@ -85,22 +120,12 @@ function check_consistency(state::TransformationState, orig_inputs) error_reporting(state, bad_idxs, n_highest_vars, iseqs, orig_inputs) end - # This is defined to check if Pantelides algorithm terminates. For more - # details, check the equation (15) of the original paper. - extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist; - map(collect, edges(var_to_diff))]) - extended_var_eq_matching = maximal_matching(extended_graph) - - nvars = ndsts(graph) - unassigned_var = [] - for (vj, eq) in enumerate(extended_var_eq_matching) - vj > nvars && break - if eq === unassigned && !isempty(𝑑neighbors(graph, vj)) - push!(unassigned_var, fullvars[vj]) - end - end + unassigned_var = singular_check(state) if !isempty(unassigned_var) || !is_balanced + if nothrow + return false + end io = IOBuffer() Base.print_array(io, unassigned_var) unassigned_var_str = String(take!(io)) @@ -110,7 +135,7 @@ function check_consistency(state::TransformationState, orig_inputs) throw(InvalidSystemException(errmsg)) end - return nothing + return true end ### diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index e715dc6eea..3a9bbd33e1 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -731,7 +731,7 @@ end function has_observed_with_lhs(sys, sym) has_observed(sys) || return false if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing - return any(isequal(sym), ic.observed_syms) + return haskey(ic.observed_syms_to_timeseries, sym) else return any(isequal(sym), [eq.lhs for eq in observed(sys)]) end @@ -740,7 +740,7 @@ end function has_parameter_dependency_with_lhs(sys, sym) has_parameter_dependencies(sys) || return false if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing - return any(isequal(sym), ic.dependent_pars) + return haskey(ic.dependent_pars_to_timeseries, unwrap(sym)) else return any(isequal(sym), [eq.lhs for eq in parameter_dependencies(sys)]) end @@ -762,11 +762,28 @@ for traitT in [ allsyms = vars(sym; op = Symbolics.Operator) for s in allsyms s = unwrap(s) - if is_variable(sys, s) || is_independent_variable(sys, s) || - has_observed_with_lhs(sys, s) + if is_variable(sys, s) || is_independent_variable(sys, s) push!(ts_idxs, ContinuousTimeseries()) elseif is_timeseries_parameter(sys, s) push!(ts_idxs, timeseries_parameter_index(sys, s).timeseries_idx) + elseif is_time_dependent(sys) && iscall(s) && issym(operation(s)) && + is_variable(sys, operation(s)(get_iv(sys))) + # DDEs case, to detect x(t - k) + push!(ts_idxs, ContinuousTimeseries()) + else + if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing + if (ts = get(ic.observed_syms_to_timeseries, s, nothing)) !== nothing + union!(ts_idxs, ts) + elseif (ts = get(ic.dependent_pars_to_timeseries, s, nothing)) !== + nothing + union!(ts_idxs, ts) + end + else + # for split=false systems + if has_observed_with_lhs(sys, sym) + push!(ts_idxs, ContinuousTimeseries()) + end + end end end end @@ -818,6 +835,8 @@ function SymbolicIndexingInterface.is_observed(sys::AbstractSystem, sym) !is_independent_variable(sys, sym) && symbolic_type(sym) != NotSymbolic() end +SymbolicIndexingInterface.supports_tuple_observed(::AbstractSystem) = true + function SymbolicIndexingInterface.observed( sys::AbstractSystem, sym; eval_expression = false, eval_module = @__MODULE__) if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing @@ -827,7 +846,8 @@ function SymbolicIndexingInterface.observed( throw(ArgumentError("Symbol $sym does not exist in the system")) end sym = _sym - elseif sym isa AbstractArray && symbolic_type(sym) isa NotSymbolic && + elseif (sym isa Tuple || + (sym isa AbstractArray && symbolic_type(sym) isa NotSymbolic)) && any(x -> x isa Symbol, sym) sym = map(sym) do s if s isa Symbol @@ -1039,6 +1059,7 @@ for prop in [:eqs :split_idxs :parent :is_dde + :tstops :index_cache :is_scalar_noise :isscheduled] @@ -1119,7 +1140,7 @@ function Base.propertynames(sys::AbstractSystem; private = false) return fieldnames(typeof(sys)) else if has_parent(sys) && (parent = get_parent(sys); parent !== nothing) - sys = parent + return propertynames(parent; private) end names = Symbol[] for s in get_systems(sys) @@ -1141,7 +1162,7 @@ end function Base.getproperty(sys::AbstractSystem, name::Symbol; namespace = !iscomplete(sys)) if has_parent(sys) && (parent = get_parent(sys); parent !== nothing) - sys = parent + return getproperty(parent, name; namespace) end wrap(getvar(sys, name; namespace = namespace)) end @@ -1365,6 +1386,14 @@ function namespace_initialization_equations( map(eq -> namespace_equation(eq, sys; ivs), eqs) end +function namespace_tstops(sys::AbstractSystem) + tstops = symbolic_tstops(sys) + isempty(tstops) && return tstops + map(tstops) do val + namespace_expr(val, sys) + end +end + function namespace_equation(eq::Equation, sys, n = nameof(sys); @@ -1620,6 +1649,14 @@ function initialization_equations(sys::AbstractSystem) end end +function symbolic_tstops(sys::AbstractSystem) + tstops = get_tstops(sys) + systems = get_systems(sys) + isempty(systems) && return tstops + tstops = [tstops; reduce(vcat, namespace_tstops.(get_systems(sys)); init = [])] + return tstops +end + function preface(sys::AbstractSystem) has_preface(sys) || return nothing pre = get_preface(sys) @@ -3273,6 +3310,97 @@ function dump_unknowns(sys::AbstractSystem) end end +""" + $(TYPEDSIGNATURES) + +Return the variable in `sys` referred to by its string representation `str`. +Roughly supports the following CFG: + +``` +varname = "D(" varname ")" | "Differential(" iv ")(" varname ")" | arrvar | maybe_dummy_var +arrvar = maybe_dummy_var "[idxs...]" +idxs = int | int "," idxs +maybe_dummy_var = namespacedvar | namespacedvar "(" iv ")" | + namespacedvar "(" iv ")" "ˍ" ts | namespacedvar "ˍ" ts | + namespacedvar "ˍ" ts "(" iv ")" +ts = iv | iv ts +namespacedvar = ident "₊" namespacedvar | ident "." namespacedvar | ident +``` + +Where `iv` is the independent variable, `int` is an integer and `ident` is an identifier. +""" +function parse_variable(sys::AbstractSystem, str::AbstractString) + iv = has_iv(sys) ? string(getname(get_iv(sys))) : nothing + + # I'd write a regex to validate `str`, but https://xkcd.com/1171/ + str = strip(str) + derivative_level = 0 + while ((cond1 = startswith(str, "D(")) || startswith(str, "Differential(")) && + endswith(str, ")") + if cond1 + derivative_level += 1 + str = _string_view_inner(str, 2, 1) + continue + end + _tmpstr = _string_view_inner(str, 13, 1) + if !startswith(_tmpstr, "$iv)(") + throw(ArgumentError("Expected differential with respect to independent variable $iv in $str")) + end + derivative_level += 1 + str = _string_view_inner(_tmpstr, length(iv) + 2, 0) + end + + arr_idxs = nothing + if endswith(str, ']') + open_idx = only(findfirst('[', str)) + idxs_range = nextind(str, open_idx):prevind(str, lastindex(str)) + idxs_str = view(str, idxs_range) + str = view(str, firstindex(str):prevind(str, open_idx)) + arr_idxs = map(Base.Fix1(parse, Int), eachsplit(idxs_str, ",")) + end + + if iv !== nothing && endswith(str, "($iv)") + str = _string_view_inner(str, 0, 2 + length(iv)) + end + + dummyderivative_level = 0 + if iv !== nothing && (dd_idx = findfirst('ˍ', str)) !== nothing + t_idx = findnext(iv, str, dd_idx) + while t_idx !== nothing + dummyderivative_level += 1 + t_idx = findnext(iv, str, nextind(str, last(t_idx))) + end + str = view(str, firstindex(str):prevind(str, dd_idx)) + end + + if iv !== nothing && endswith(str, "($iv)") + str = _string_view_inner(str, 0, 2 + length(iv)) + end + + cur = sys + for ident in eachsplit(str, ('.', NAMESPACE_SEPARATOR)) + ident = Symbol(ident) + hasproperty(cur, ident) || + throw(ArgumentError("System $(nameof(cur)) does not have a subsystem/variable named $(ident)")) + cur = getproperty(cur, ident) + end + + if arr_idxs !== nothing + cur = cur[arr_idxs...] + end + + for i in 1:(derivative_level + dummyderivative_level) + cur = Differential(get_iv(sys))(cur) + end + + return cur +end + +function _string_view_inner(str, startoffset, endoffset) + view(str, + nextind(str, firstindex(str), startoffset):prevind(str, lastindex(str), endoffset)) +end + ### Functions for accessing algebraic/differential equations in systems ### """ diff --git a/src/systems/callbacks.jl b/src/systems/callbacks.jl index b9423b8f07..0cb7f16c9f 100644 --- a/src/systems/callbacks.jl +++ b/src/systems/callbacks.jl @@ -1,7 +1,9 @@ #################################### system operations ##################################### -get_continuous_events(sys::AbstractSystem) = SymbolicContinuousCallback[] -get_continuous_events(sys::AbstractODESystem) = getfield(sys, :continuous_events) has_continuous_events(sys::AbstractSystem) = isdefined(sys, :continuous_events) +function get_continuous_events(sys::AbstractSystem) + has_continuous_events(sys) || return SymbolicContinuousCallback[] + getfield(sys, :continuous_events) +end has_discrete_events(sys::AbstractSystem) = isdefined(sys, :discrete_events) function get_discrete_events(sys::AbstractSystem) @@ -795,8 +797,8 @@ function compile_affect(eqs::Vector{Equation}, cb, sys, dvs, ps; outputidxs = no end end -function generate_rootfinding_callback(sys::AbstractODESystem, dvs = unknowns(sys), - ps = parameters(sys); kwargs...) +function generate_rootfinding_callback(sys::AbstractTimeDependentSystem, + dvs = unknowns(sys), ps = parameters(sys); kwargs...) cbs = continuous_events(sys) isempty(cbs) && return nothing generate_rootfinding_callback(cbs, sys, dvs, ps; kwargs...) @@ -806,7 +808,7 @@ Generate a single rootfinding callback; this happens if there is only one equati generate_rootfinding_callback and thus we can produce a ContinuousCallback instead of a VectorContinuousCallback. """ function generate_single_rootfinding_callback( - eq, cb, sys::AbstractODESystem, dvs = unknowns(sys), + eq, cb, sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = parameters(sys); kwargs...) if !isequal(eq.lhs, 0) eq = 0 ~ eq.lhs - eq.rhs @@ -849,7 +851,7 @@ function generate_single_rootfinding_callback( end function generate_vector_rootfinding_callback( - cbs, sys::AbstractODESystem, dvs = unknowns(sys), + cbs, sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = parameters(sys); rootfind = SciMLBase.RightRootFind, reinitialization = SciMLBase.CheckInit(), kwargs...) eqs = map(cb -> flatten_equations(cb.eqs), cbs) @@ -950,7 +952,7 @@ end """ Compile a single continuous callback affect function(s). """ -function compile_affect_fn(cb, sys::AbstractODESystem, dvs, ps, kwargs) +function compile_affect_fn(cb, sys::AbstractTimeDependentSystem, dvs, ps, kwargs) eq_aff = affects(cb) eq_neg_aff = affect_negs(cb) affect = compile_affect(eq_aff, cb, sys, dvs, ps; expression = Val{false}, kwargs...) @@ -974,8 +976,8 @@ function compile_affect_fn(cb, sys::AbstractODESystem, dvs, ps, kwargs) (affect = affect, affect_neg = affect_neg, initialize = initialize, finalize = finalize) end -function generate_rootfinding_callback(cbs, sys::AbstractODESystem, dvs = unknowns(sys), - ps = parameters(sys); kwargs...) +function generate_rootfinding_callback(cbs, sys::AbstractTimeDependentSystem, + dvs = unknowns(sys), ps = parameters(sys); kwargs...) eqs = map(cb -> flatten_equations(cb.eqs), cbs) num_eqs = length.(eqs) total_eqs = sum(num_eqs) @@ -1315,12 +1317,12 @@ merge_cb(x, ::Nothing) = x merge_cb(x, y) = CallbackSet(x, y) function process_events(sys; callback = nothing, kwargs...) - if has_continuous_events(sys) + if has_continuous_events(sys) && !isempty(continuous_events(sys)) contin_cb = generate_rootfinding_callback(sys; kwargs...) else contin_cb = nothing end - if has_discrete_events(sys) + if has_discrete_events(sys) && !isempty(discrete_events(sys)) discrete_cb = generate_discrete_callbacks(sys; kwargs...) else discrete_cb = nothing diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 518b48e929..af96c4fcfe 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -71,8 +71,6 @@ function calculate_jacobian(sys::AbstractODESystem; rhs = [eq.rhs - eq.lhs for eq in full_equations(sys)] #need du terms on rhs for differentiating wrt du - iv = get_iv(sys) - if sparse jac = sparsejacobian(rhs, dvs, simplify = simplify) else @@ -94,8 +92,6 @@ function calculate_control_jacobian(sys::AbstractODESystem; end rhs = [eq.rhs for eq in full_equations(sys)] - - iv = get_iv(sys) ctrls = controls(sys) if sparse @@ -221,9 +217,11 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys), pre, sol_states = get_substitutions_and_solved_unknowns(sys) if implicit_dae + # inputs = [] makes `wrap_array_vars` offset by 1 since there is an extra + # argument build_function(rhss, ddvs, u, p..., t; postprocess_fbody = pre, states = sol_states, - wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps) .∘ + wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs, ps, inputs = []) .∘ wrap_parameter_dependencies(sys, false), kwargs...) else @@ -500,6 +498,8 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) checkbounds = false, initializeprob = nothing, initializeprobmap = nothing, + initializeprobpmap = nothing, + update_initializeprob! = nothing, kwargs...) where {iip} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEFunction`") @@ -553,7 +553,9 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = unknowns(sys) jac_prototype = jac_prototype, observed = observedfun, initializeprob = initializeprob, - initializeprobmap = initializeprobmap) + initializeprobmap = initializeprobmap, + initializeprobpmap = initializeprobpmap, + update_initializeprob! = update_initializeprob!) end function DiffEqBase.DDEFunction(sys::AbstractODESystem, args...; kwargs...) @@ -756,6 +758,39 @@ function DAEFunctionExpr(sys::AbstractODESystem, args...; kwargs...) DAEFunctionExpr{true}(sys, args...; kwargs...) end +struct SymbolicTstops{F} + fn::F +end + +function (st::SymbolicTstops)(p, tspan) + unique!(sort!(reduce(vcat, st.fn(p..., tspan...)))) +end + +function SymbolicTstops( + sys::AbstractSystem; eval_expression = false, eval_module = @__MODULE__) + tstops = symbolic_tstops(sys) + isempty(tstops) && return nothing + t0 = gensym(:t0) + t1 = gensym(:t1) + tstops = map(tstops) do val + if is_array_of_symbolics(val) || val isa AbstractArray + collect(val) + else + term(:, t0, unwrap(val), t1; type = AbstractArray{Real}) + end + end + rps = reorder_parameters(sys, parameters(sys)) + tstops, _ = build_function(tstops, + rps..., + t0, + t1; + expression = Val{true}, + wrap_code = wrap_array_vars(sys, tstops; dvs = nothing) .∘ + wrap_parameter_dependencies(sys, false)) + tstops = eval_or_rgf(tstops; eval_expression, eval_module) + return SymbolicTstops(tstops) +end + """ ```julia DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem, u0map, tspan, @@ -790,12 +825,6 @@ function DiffEqBase.ODEProblem{false}(sys::AbstractODESystem, args...; kwargs... ODEProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) end -struct DiscreteSaveAffect{F, S} <: Function - f::F - s::S -end -(d::DiscreteSaveAffect)(args...) = d.f(args..., d.s) - function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], tspan = get_tspan(sys), parammap = DiffEqBase.NullParameters(); @@ -821,6 +850,11 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = kwargs1 = merge(kwargs1, (callback = cbs,)) end + tstops = SymbolicTstops(sys; eval_expression, eval_module) + if tstops !== nothing + kwargs1 = merge(kwargs1, (; tstops)) + end + return ODEProblem{iip}(f, u0, tspan, p, pt; kwargs1..., kwargs...) end get_callback(prob::ODEProblem) = prob.kwargs[:callback] @@ -847,7 +881,7 @@ end function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, parammap = DiffEqBase.NullParameters(); warn_initialize_determined = true, - check_length = true, kwargs...) where {iip} + check_length = true, eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating a `DAEProblem`") end @@ -860,8 +894,15 @@ function DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan differential_vars = map(Base.Fix2(in, diffvars), sts) kwargs = filter_kwargs(kwargs) + kwargs1 = (;) + + tstops = SymbolicTstops(sys; eval_expression, eval_module) + if tstops !== nothing + kwargs1 = merge(kwargs1, (; tstops)) + end + DAEProblem{iip}(f, du0, u0, tspan, p; differential_vars = differential_vars, - kwargs...) + kwargs..., kwargs1...) end function generate_history(sys::AbstractODESystem, u0; expression = Val{false}, kwargs...) @@ -1258,7 +1299,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, check_length = true, warn_initialize_determined = true, initialization_eqs = [], - fully_determined = false, + fully_determined = nothing, check_units = true, kwargs...) where {iip, specialize} if !iscomplete(sys) @@ -1269,11 +1310,24 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, elseif isempty(u0map) && get_initializesystem(sys) === nothing isys = structural_simplify( generate_initializesystem( - sys; initialization_eqs, check_units, pmap = parammap); fully_determined) + sys; initialization_eqs, check_units, pmap = parammap, guesses); fully_determined) else isys = structural_simplify( generate_initializesystem( - sys; u0map, initialization_eqs, check_units, pmap = parammap); fully_determined) + sys; u0map, initialization_eqs, check_units, pmap = parammap, guesses); fully_determined) + end + + ts = get_tearing_state(isys) + if warn_initialize_determined && + (unassigned_vars = StructuralTransformations.singular_check(ts); !isempty(unassigned_vars)) + errmsg = """ + The initialization system is structurally singular. Guess values may \ + significantly affect the initial values of the ODE. The problematic variables \ + are $unassigned_vars. + + Note that the identification of problematic variables is a best-effort heuristic. + """ + @warn errmsg end uninit = setdiff(unknowns(sys), [unknowns(isys); getfield.(observed(isys), :lhs)]) @@ -1319,6 +1373,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, u0T = promote_type(u0T, typeof(fullmap[eq.lhs])) end if u0T != Union{} + u0T = eltype(u0T) u0map = Dict(k => if symbolic_type(v) == NotSymbolic() && !is_array_of_symbolics(v) v isa AbstractArray ? u0T.(v) : u0T(v) else diff --git a/src/systems/diffeqs/modelingtoolkitize.jl b/src/systems/diffeqs/modelingtoolkitize.jl index 68a970d8b5..2e0623bea1 100644 --- a/src/systems/diffeqs/modelingtoolkitize.jl +++ b/src/systems/diffeqs/modelingtoolkitize.jl @@ -277,10 +277,24 @@ function modelingtoolkitize(prob::DiffEqBase.SDEProblem; kwargs...) else Vector(vec(params)) end + sts = Vector(vec(vars)) + default_u0 = Dict(sts .=> vec(collect(prob.u0))) + default_p = if has_p + if prob.p isa AbstractDict + Dict(v => prob.p[k] for (k, v) in pairs(_params)) + elseif prob.p isa MTKParameters + Dict(params .=> reduce(vcat, prob.p)) + else + Dict(params .=> vec(collect(prob.p))) + end + else + Dict() + end - de = SDESystem(deqs, neqs, t, Vector(vec(vars)), params; + de = SDESystem(deqs, neqs, t, sts, params; name = gensym(:MTKizedSDE), tspan = prob.tspan, + defaults = merge(default_u0, default_p), kwargs...) de diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index e211d6abd2..4040b7a646 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -149,6 +149,11 @@ struct ODESystem <: AbstractODESystem """ is_dde::Bool """ + A list of points to provide to the solver as tstops. Uses the same syntax as discrete + events. + """ + tstops::Vector{Any} + """ Cache for intermediate tearing state. """ tearing_state::Any @@ -187,7 +192,7 @@ struct ODESystem <: AbstractODESystem connector_type, preface, cevents, devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing, is_dde = false, - tearing_state = nothing, + tstops = [], tearing_state = nothing, substitutions = nothing, complete = false, index_cache = nothing, discrete_subsystems = nothing, solved_unknowns = nothing, split_idxs = nothing, parent = nothing; checks::Union{Bool, Int} = true) @@ -206,7 +211,7 @@ struct ODESystem <: AbstractODESystem ctrl_jac, Wfact, Wfact_t, name, description, systems, defaults, guesses, torn_matching, initializesystem, initialization_eqs, schedule, connector_type, preface, cevents, devents, parameter_dependencies, metadata, - gui_metadata, is_dde, tearing_state, substitutions, complete, index_cache, + gui_metadata, is_dde, tstops, tearing_state, substitutions, complete, index_cache, discrete_subsystems, solved_unknowns, split_idxs, parent) end end @@ -233,7 +238,8 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; checks = true, metadata = nothing, gui_metadata = nothing, - is_dde = nothing) + is_dde = nothing, + tstops = []) name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro")) @assert all(control -> any(isequal.(control, ps)), controls) "All controls must also be parameters." @@ -299,7 +305,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; defaults, guesses, nothing, initializesystem, initialization_eqs, schedule, connector_type, preface, cont_callbacks, disc_callbacks, parameter_dependencies, - metadata, gui_metadata, is_dde, checks = checks) + metadata, gui_metadata, is_dde, tstops, checks = checks) end function ODESystem(eqs, iv; kwargs...) @@ -402,6 +408,7 @@ function flatten(sys::ODESystem, noeqs = false) description = description(sys), initialization_eqs = initialization_equations(sys), is_dde = is_dde(sys), + tstops = symbolic_tstops(sys), checks = false) end end @@ -409,34 +416,48 @@ end ODESystem(eq::Equation, args...; kwargs...) = ODESystem([eq], args...; kwargs...) """ -$(SIGNATURES) + build_explicit_observed_function(sys, ts; kwargs...) -> Function(s) + +Generates a function that computes the observed value(s) `ts` in the system `sys`, while making the assumption that there are no cycles in the equations. + +## Arguments +- `sys`: The system for which to generate the function +- `ts`: The symbolic observed values whose value should be computed + +## Keywords +- `return_inplace = false`: If true and the observed value is a vector, then return both the in place and out of place methods. +- `expression = false`: Generates a Julia `Expr`` computing the observed value if `expression` is true +- `eval_expression = false`: If true and `expression = false`, evaluates the returned function in the module `eval_module` +- `output_type = Array` the type of the array generated by a out-of-place vector-valued function +- `param_only = false` if true, only allow the generated function to access system parameters +- `inputs = nothing` additinoal symbolic variables that should be provided to the generated function +- `checkbounds = true` checks bounds if true when destructuring parameters +- `op = Operator` sets the recursion terminator for the walk done by `vars` to identify the variables that appear in `ts`. See the documentation for `vars` for more detail. +- `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist. +- `mkarray`; only used if the output is an array (that is, `!isscalar(ts)` and `ts` is not a tuple, in which case the result will always be a tuple). Called as `mkarray(ts, output_type)` where `ts` are the expressions to put in +the array and `output_type` is the argument of the same name passed to build_explicit_observed_function. + +## Returns -Generates a function that computes the observed value(s) `ts` in the system `sys` assuming that there are no cycles in the equations. - The return value will be either: -* a single function if the input is a scalar or if the input is a Vector but `return_inplace` is false -* the out of place and in-place functions `(ip, oop)` if `return_inplace` is true and the input is a `Vector` +* a single function `f_oop` if the input is a scalar or if the input is a Vector but `return_inplace` is false +* the out of place and in-place functions `(f_ip, f_oop)` if `return_inplace` is true and the input is a `Vector` -The function(s) will be: +The function(s) `f_oop` (and potentially `f_ip`) will be: * `RuntimeGeneratedFunction`s by default, * A Julia `Expr` if `expression` is true, -* A directly evaluated Julia function in the module `eval_module` if `eval_expression` is true +* A directly evaluated Julia function in the module `eval_module` if `eval_expression` is true and `expression` is false. The signatures will be of the form `g(...)` with arguments: -* `output` for in-place functions -* `unknowns` if `params_only` is `false` -* `inputs` if `inputs` is an array of symbolic inputs that should be available in `ts` -* `p...` unconditionally; note that in the case of `MTKParameters` more than one parameters argument may be present, so it must be splatted -* `t` if the system is time-dependent; for example `NonlinearSystem` will not have `t` -For example, a function `g(op, unknowns, p, inputs, t)` will be the in-place function generated if `return_inplace` is true, `ts` is a vector, an array of inputs `inputs` is given, and `params_only` is false for a time-dependent system. - -Options not otherwise specified are: -* `output_type = Array` the type of the array generated by the out-of-place vector-valued function -* `checkbounds = true` checks bounds if true when destructuring parameters -* `op = Operator` sets the recursion terminator for the walk done by `vars` to identify the variables that appear in `ts`. See the documentation for `vars` for more detail. -* `throw = true` if true, throw an error when generating a function for `ts` that reference variables that do not exist -* `drop_expr` is deprecated. -* `array_type`; only used if the output is an array (that is, `!isscalar(ts)`). If it is `Vector`, then it will generate an array, if `Tuple` then it will generate a tuple. + +- `output` for in-place functions +- `unknowns` if `param_only` is `false` +- `inputs` if `inputs` is an array of symbolic inputs that should be available in `ts` +- `p...` unconditionally; note that in the case of `MTKParameters` more than one parameters argument may be present, so it must be splatted +- `t` if the system is time-dependent; for example `NonlinearSystem` will not have `t` + +For example, a function `g(op, unknowns, p..., inputs, t)` will be the in-place function generated if `return_inplace` is true, `ts` is a vector, +an array of inputs `inputs` is given, and `param_only` is false for a time-dependent system. """ function build_explicit_observed_function(sys, ts; inputs = nothing, @@ -445,13 +466,16 @@ function build_explicit_observed_function(sys, ts; eval_module = @__MODULE__, output_type = Array, checkbounds = true, - drop_expr = drop_expr, ps = parameters(sys), return_inplace = false, param_only = false, op = Operator, throw = true, - array_type = Vector) + mkarray = MakeArray) + is_tuple = ts isa Tuple + if is_tuple + ts = collect(ts) + end if (isscalar = symbolic_type(ts) !== NotSymbolic()) ts = [ts] end @@ -473,7 +497,16 @@ function build_explicit_observed_function(sys, ts; ivs = independent_variables(sys) dep_vars = scalarize(setdiff(vars, ivs)) - obs = param_only ? Equation[] : observed(sys) + obs = observed(sys) + if param_only + if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing + obs = filter(obs) do eq + !(ContinuousTimeseries() in ic.observed_syms_to_timeseries[eq.lhs]) + end + else + obs = Equation[] + end + end cs = collect_constants(obs) if !isempty(cs) > 0 @@ -600,8 +633,17 @@ function build_explicit_observed_function(sys, ts; (array_type <: Vector ? MakeArray(ts, output_type) : MakeTuple(ts)) # Need to keep old method of building the function since it uses `output_type`, # which can't be provided to `build_function` - oop_fn = Func(args, [], pre(Let(obsexprs, output_expr, false))) |> array_wrapper[1] |> - oop_mtkp_wrapper |> toexpr + return_value = if isscalar + ts[1] + elseif is_tuple + MakeTuple(Tuple(ts)) + else + mkarray(ts, output_type) + end + oop_fn = Func(args, [], + pre(Let(obsexprs, + return_value, + false))) |> array_wrapper[1] |> oop_mtkp_wrapper |> toexpr oop_fn = expression ? oop_fn : eval_or_rgf(oop_fn; eval_expression, eval_module) if !isscalar diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index bd295055fa..4b0a342078 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -38,29 +38,33 @@ const UnknownIndexMap = Dict{ BasicSymbolic, Union{Int, UnitRange{Int}, AbstractArray{Int}}} const TunableIndexMap = Dict{BasicSymbolic, Union{Int, UnitRange{Int}, Base.ReshapedArray{Int, N, UnitRange{Int}} where {N}}} +const TimeseriesSetType = Set{Union{ContinuousTimeseries, Int}} + +const SymbolicParam = Union{BasicSymbolic, CallWithMetadata} struct IndexCache unknown_idx::UnknownIndexMap # sym => (bufferidx, idx_in_buffer) - discrete_idx::Dict{BasicSymbolic, DiscreteIndex} + discrete_idx::Dict{SymbolicParam, DiscreteIndex} # sym => (clockidx, idx_in_clockbuffer) callback_to_clocks::Dict{Any, Vector{Int}} tunable_idx::TunableIndexMap constant_idx::ParamIndexMap nonnumeric_idx::NonnumericMap - observed_syms::Set{BasicSymbolic} - dependent_pars::Set{Union{BasicSymbolic, CallWithMetadata}} + observed_syms_to_timeseries::Dict{BasicSymbolic, TimeseriesSetType} + dependent_pars_to_timeseries::Dict{ + Union{BasicSymbolic, CallWithMetadata}, TimeseriesSetType} discrete_buffer_sizes::Vector{Vector{BufferTemplate}} tunable_buffer_size::BufferTemplate constant_buffer_sizes::Vector{BufferTemplate} nonnumeric_buffer_sizes::Vector{BufferTemplate} - symbol_to_variable::Dict{Symbol, Union{BasicSymbolic, CallWithMetadata}} + symbol_to_variable::Dict{Symbol, SymbolicParam} end function IndexCache(sys::AbstractSystem) unks = solved_unknowns(sys) unk_idxs = UnknownIndexMap() - symbol_to_variable = Dict{Symbol, Union{BasicSymbolic, CallWithMetadata}}() + symbol_to_variable = Dict{Symbol, SymbolicParam}() let idx = 1 for sym in unks @@ -91,23 +95,9 @@ function IndexCache(sys::AbstractSystem) end end - observed_syms = Set{BasicSymbolic}() - for eq in observed(sys) - if symbolic_type(eq.lhs) != NotSymbolic() - sym = eq.lhs - ttsym = default_toterm(sym) - rsym = renamespace(sys, sym) - rttsym = renamespace(sys, ttsym) - push!(observed_syms, sym) - push!(observed_syms, ttsym) - push!(observed_syms, rsym) - push!(observed_syms, rttsym) - end - end - tunable_buffers = Dict{Any, Set{BasicSymbolic}}() constant_buffers = Dict{Any, Set{BasicSymbolic}}() - nonnumeric_buffers = Dict{Any, Set{Union{BasicSymbolic, CallWithMetadata}}}() + nonnumeric_buffers = Dict{Any, Set{SymbolicParam}}() function insert_by_type!(buffers::Dict{Any, S}, sym, ctype) where {S} sym = unwrap(sym) @@ -115,10 +105,10 @@ function IndexCache(sys::AbstractSystem) push!(buf, sym) end - disc_param_callbacks = Dict{BasicSymbolic, Set{Int}}() + disc_param_callbacks = Dict{SymbolicParam, Set{Int}}() events = vcat(continuous_events(sys), discrete_events(sys)) for (i, event) in enumerate(events) - discs = Set{BasicSymbolic}() + discs = Set{SymbolicParam}() affs = affects(event) if !(affs isa AbstractArray) affs = [affs] @@ -142,26 +132,32 @@ function IndexCache(sys::AbstractSystem) isequal(only(arguments(sym)), get_iv(sys)) clocks = get!(() -> Set{Int}(), disc_param_callbacks, sym) push!(clocks, i) - else + elseif is_variable_floatingpoint(sym) insert_by_type!(constant_buffers, sym, symtype(sym)) + else + stype = symtype(sym) + if stype <: FnType + stype = fntype_to_function_type(stype) + end + insert_by_type!(nonnumeric_buffers, sym, stype) end end end clock_partitions = unique(collect(values(disc_param_callbacks))) disc_symtypes = unique(symtype.(keys(disc_param_callbacks))) disc_symtype_idx = Dict(disc_symtypes .=> eachindex(disc_symtypes)) - disc_syms_by_symtype = [BasicSymbolic[] for _ in disc_symtypes] + disc_syms_by_symtype = [SymbolicParam[] for _ in disc_symtypes] for sym in keys(disc_param_callbacks) push!(disc_syms_by_symtype[disc_symtype_idx[symtype(sym)]], sym) end - disc_syms_by_symtype_by_partition = [Vector{BasicSymbolic}[] for _ in disc_symtypes] + disc_syms_by_symtype_by_partition = [Vector{SymbolicParam}[] for _ in disc_symtypes] for (i, buffer) in enumerate(disc_syms_by_symtype) for partition in clock_partitions push!(disc_syms_by_symtype_by_partition[i], [sym for sym in buffer if disc_param_callbacks[sym] == partition]) end end - disc_idxs = Dict{BasicSymbolic, DiscreteIndex}() + disc_idxs = Dict{SymbolicParam, DiscreteIndex}() callback_to_clocks = Dict{ Union{SymbolicContinuousCallback, SymbolicDiscreteCallback}, Set{Int}}() for (typei, disc_syms_by_partition) in enumerate(disc_syms_by_symtype_by_partition) @@ -203,6 +199,7 @@ function IndexCache(sys::AbstractSystem) end haskey(disc_idxs, p) && continue haskey(constant_buffers, ctype) && p in constant_buffers[ctype] && continue + haskey(nonnumeric_buffers, ctype) && p in nonnumeric_buffers[ctype] && continue insert_by_type!( if ctype <: Real || ctype <: AbstractArray{<:Real} if istunable(p, true) && Symbolics.shape(p) != Symbolics.Unknown() && @@ -267,29 +264,68 @@ function IndexCache(sys::AbstractSystem) end end - for sym in Iterators.flatten((keys(unk_idxs), keys(disc_idxs), keys(tunable_idxs), - keys(const_idxs), keys(nonnumeric_idxs), - observed_syms, independent_variable_symbols(sys))) - if hasname(sym) && (!iscall(sym) || operation(sym) !== getindex) - symbol_to_variable[getname(sym)] = sym - end - end - - dependent_pars = Set{Union{BasicSymbolic, CallWithMetadata}}() + dependent_pars_to_timeseries = Dict{ + Union{BasicSymbolic, CallWithMetadata}, TimeseriesSetType}() for eq in parameter_dependencies(sys) sym = eq.lhs + vs = vars(eq.rhs) + timeseries = TimeseriesSetType() + if is_time_dependent(sys) + for v in vs + if (idx = get(disc_idxs, v, nothing)) !== nothing + push!(timeseries, idx.clock_idx) + end + end + end ttsym = default_toterm(sym) rsym = renamespace(sys, sym) rttsym = renamespace(sys, ttsym) - for s in [sym, ttsym, rsym, rttsym] - push!(dependent_pars, s) + for s in (sym, ttsym, rsym, rttsym) + dependent_pars_to_timeseries[s] = timeseries if hasname(s) && (!iscall(s) || operation(s) != getindex) symbol_to_variable[getname(s)] = sym end end end + observed_syms_to_timeseries = Dict{BasicSymbolic, TimeseriesSetType}() + for eq in observed(sys) + if symbolic_type(eq.lhs) != NotSymbolic() + sym = eq.lhs + vs = vars(eq.rhs; op = Nothing) + timeseries = TimeseriesSetType() + if is_time_dependent(sys) + for v in vs + if (idx = get(disc_idxs, v, nothing)) !== nothing + push!(timeseries, idx.clock_idx) + elseif haskey(observed_syms_to_timeseries, v) + union!(timeseries, observed_syms_to_timeseries[v]) + elseif haskey(dependent_pars_to_timeseries, v) + union!(timeseries, dependent_pars_to_timeseries[v]) + end + end + if isempty(timeseries) + push!(timeseries, ContinuousTimeseries()) + end + end + ttsym = default_toterm(sym) + rsym = renamespace(sys, sym) + rttsym = renamespace(sys, ttsym) + for s in (sym, ttsym, rsym, rttsym) + observed_syms_to_timeseries[s] = timeseries + end + end + end + + for sym in Iterators.flatten((keys(unk_idxs), keys(disc_idxs), keys(tunable_idxs), + keys(const_idxs), keys(nonnumeric_idxs), + keys(observed_syms_to_timeseries), independent_variable_symbols(sys))) + if hasname(sym) && (!iscall(sym) || operation(sym) !== getindex) + symbol_to_variable[getname(sym)] = sym + end + end + return IndexCache( unk_idxs, disc_idxs, @@ -297,8 +333,8 @@ function IndexCache(sys::AbstractSystem) tunable_idxs, const_idxs, nonnumeric_idxs, - observed_syms, - dependent_pars, + observed_syms_to_timeseries, + dependent_pars_to_timeseries, disc_buffer_templates, BufferTemplate(Real, tunable_buffer_size), const_buffer_sizes, diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl index 9fa073b4b0..e5e17fb5f9 100644 --- a/src/systems/jumps/jumpsystem.jl +++ b/src/systems/jumps/jumpsystem.jl @@ -88,6 +88,11 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem """ connector_type::Any """ + A `Vector{SymbolicContinuousCallback}` that model events. + The integrator will use root finding to guarantee that it steps at each zero crossing. + """ + continuous_events::Vector{SymbolicContinuousCallback} + """ A `Vector{SymbolicDiscreteCallback}` that models events. Symbolic analog to `SciMLBase.DiscreteCallback` that executes an affect when a given condition is true at the end of an integration step. Note, one must make sure to call @@ -120,8 +125,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem function JumpSystem{U}( tag, ap::U, iv, unknowns, ps, var_to_name, observed, name, description, - systems, - defaults, connector_type, devents, parameter_dependencies, + systems, defaults, connector_type, cevents, devents, parameter_dependencies, metadata = nothing, gui_metadata = nothing, complete = false, index_cache = nothing, isscheduled = false; checks::Union{Bool, Int} = true) where {U <: ArrayPartition} @@ -136,8 +140,8 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem end new{U}(tag, ap, iv, unknowns, ps, var_to_name, observed, name, description, systems, defaults, - connector_type, devents, parameter_dependencies, metadata, gui_metadata, - complete, index_cache, isscheduled) + connector_type, cevents, devents, parameter_dependencies, metadata, + gui_metadata, complete, index_cache, isscheduled) end end function JumpSystem(tag, ap, iv, states, ps, var_to_name, args...; kwargs...) @@ -194,7 +198,8 @@ function JumpSystem(eqs, iv, unknowns, ps; # this and the treatment of continuous events are the only part # unique to JumpSystems eqs = scalarize.(eqs) - ap = ArrayPartition(MassActionJump[], ConstantRateJump[], VariableRateJump[]) + ap = ArrayPartition( + MassActionJump[], ConstantRateJump[], VariableRateJump[], Equation[]) for eq in eqs if eq isa MassActionJump push!(ap.x[1], eq) @@ -202,18 +207,19 @@ function JumpSystem(eqs, iv, unknowns, ps; push!(ap.x[2], eq) elseif eq isa VariableRateJump push!(ap.x[3], eq) + elseif eq isa Equation + push!(ap.x[4], eq) else - error("JumpSystem equations must contain MassActionJumps, ConstantRateJumps, or VariableRateJumps.") + error("JumpSystem equations must contain MassActionJumps, ConstantRateJumps, VariableRateJumps, or Equations.") end end - (continuous_events === nothing) || - error("JumpSystems currently only support discrete events.") + cont_callbacks = SymbolicContinuousCallbacks(continuous_events) disc_callbacks = SymbolicDiscreteCallbacks(discrete_events) JumpSystem{typeof(ap)}(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), ap, iv′, us′, ps′, var_to_name, observed, name, description, systems, - defaults, connector_type, disc_callbacks, parameter_dependencies, + defaults, connector_type, cont_callbacks, disc_callbacks, parameter_dependencies, metadata, gui_metadata, checks = checks) end @@ -245,6 +251,7 @@ end has_massactionjumps(js::JumpSystem) = !isempty(equations(js).x[1]) has_constantratejumps(js::JumpSystem) = !isempty(equations(js).x[2]) has_variableratejumps(js::JumpSystem) = !isempty(equations(js).x[3]) +has_equations(js::JumpSystem) = !isempty(equations(js).x[4]) function generate_rate_function(js::JumpSystem, rate) consts = collect_constants(rate) @@ -281,7 +288,7 @@ function assemble_vrj( outputidxs = [unknowntoid[var] for var in outputvars] affect = eval_or_rgf(generate_affect_function(js, vrj.affect!, outputidxs); eval_expression, eval_module) - VariableRateJump(rate, affect) + VariableRateJump(rate, affect; save_positions = vrj.save_positions) end function assemble_vrj_expr(js, vrj, unknowntoid) @@ -390,6 +397,11 @@ function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, if !iscomplete(sys) error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`") end + + if has_equations(sys) || (!isempty(continuous_events(sys))) + error("The passed in JumpSystem contains `Equation`s or continuous events, please use a problem type that supports these features, such as ODEProblem.") + end + _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap; t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, check_length = false) f = DiffEqBase.DISCRETE_INPLACE_DEFAULT @@ -478,14 +490,24 @@ function DiffEqBase.ODEProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothi error("A completed `JumpSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `DiscreteProblem`") end - _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap; - t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, check_length = false) - - observedfun = ObservedFunctionCache(sys; eval_expression, eval_module) - - f = (du, u, p, t) -> (du .= 0; nothing) - df = ODEFunction(f; sys, observed = observedfun) - ODEProblem(df, u0, tspan, p; kwargs...) + # forward everything to be an ODESystem but the jumps and discrete events + if has_equations(sys) + osys = ODESystem(equations(sys).x[4], get_iv(sys), unknowns(sys), parameters(sys); + observed = observed(sys), name = nameof(sys), description = description(sys), + systems = get_systems(sys), defaults = defaults(sys), + parameter_dependencies = parameter_dependencies(sys), + metadata = get_metadata(sys), gui_metadata = get_gui_metadata(sys)) + osys = complete(osys) + return ODEProblem(osys, u0map, tspan, parammap; check_length = false, kwargs...) + else + _, u0, p = process_SciMLProblem(EmptySciMLFunction, sys, u0map, parammap; + t = tspan === nothing ? nothing : tspan[1], use_union, tofloat = false, + check_length = false) + f = (du, u, p, t) -> (du .= 0; nothing) + observedfun = ObservedFunctionCache(sys; eval_expression, eval_module) + df = ODEFunction(f; sys, observed = observedfun) + return ODEProblem(df, u0, tspan, p; kwargs...) + end end """ @@ -521,8 +543,11 @@ function JumpProcesses.JumpProblem(js::JumpSystem, prob, for j in eqs.x[2]] vrjs = VariableRateJump[assemble_vrj(js, j, unknowntoid; eval_expression, eval_module) for j in eqs.x[3]] - ((prob isa DiscreteProblem) && !isempty(vrjs)) && - error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps") + if prob isa DiscreteProblem + if (!isempty(vrjs) || has_equations(js) || !isempty(continuous_events(js))) + error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps, coupled differential equations, or continuous events.") + end + end jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, majs) # dep graphs are only for constant rate jumps diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl index 5b79eaa91b..5402659b1c 100644 --- a/src/systems/model_parsing.jl +++ b/src/systems/model_parsing.jl @@ -923,7 +923,7 @@ function parse_variable_arg(dict, mod, arg, varclass, kwargs, where_types) end push!(varexpr.args, metadata_expr) - return vv isa Num ? name : :($name...), varexpr + return symbolic_type(vv) == ScalarSymbolic() ? name : :($name...), varexpr else return vv end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index eefe393acc..726d171bd0 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -30,7 +30,8 @@ function generate_initializesystem(sys::ODESystem; # 1) process dummy derivatives and u0map into initialization system eqs_ics = eqs[idxs_alge] # start equation list with algebraic equations defs = copy(defaults(sys)) # copy so we don't modify sys.defaults - guesses = merge(get_guesses(sys), todict(guesses)) + additional_guesses = anydict(guesses) + guesses = merge(get_guesses(sys), additional_guesses) schedule = getfield(sys, :schedule) if !isnothing(schedule) for x in filter(x -> !isnothing(x[1]), schedule.dummy_sub) @@ -178,7 +179,7 @@ function generate_initializesystem(sys::ODESystem; for k in keys(defs) defs[k] = substitute(defs[k], paramsubs) end - meta = InitializationSystemMetadata(Dict{Any, Any}(u0map), Dict{Any, Any}(pmap)) + meta = InitializationSystemMetadata(anydict(u0map), anydict(pmap), additional_guesses) return NonlinearSystem(eqs_ics, vars, pars; @@ -193,6 +194,7 @@ end struct InitializationSystemMetadata u0map::Dict{Any, Any} pmap::Dict{Any, Any} + additional_guesses::Dict{Any, Any} end function is_parameter_solvable(p, pmap, defs, guesses) @@ -208,24 +210,23 @@ function is_parameter_solvable(p, pmap, defs, guesses) _val1 === nothing && _val2 !== nothing)) && _val3 !== nothing end -function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p) +function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p, newu0, newp) if u0 === missing && p === missing - return odefn.initializeprob, odefn.update_initializeprob!, odefn.initializeprobmap, - odefn.initializeprobpmap + return odefn.initialization_data end if !(eltype(u0) <: Pair) && !(eltype(p) <: Pair) oldinitprob = odefn.initializeprob - if oldinitprob === nothing || !SciMLBase.has_sys(oldinitprob.f) || - !(oldinitprob.f.sys isa NonlinearSystem) - return oldinitprob, odefn.update_initializeprob!, odefn.initializeprobmap, - odefn.initializeprobpmap + oldinitprob === nothing && return nothing + if !SciMLBase.has_sys(oldinitprob.f) || !(oldinitprob.f.sys isa NonlinearSystem) + return SciMLBase.OverrideInitData(oldinitprob, odefn.update_initializeprob!, + odefn.initializeprobmap, odefn.initializeprobpmap) end pidxs = ParameterIndex[] pvals = [] u0idxs = Int[] u0vals = [] for sym in variable_symbols(oldinitprob) - if is_variable(sys, sym) + if is_variable(sys, sym) || has_observed_with_lhs(sys, sym) u0 !== missing || continue idx = variable_index(oldinitprob, sym) push!(u0idxs, idx) @@ -247,81 +248,91 @@ function SciMLBase.remake_initializeprob(sys::ODESystem, odefn, u0, t0, p) end end end - newu0 = remake_buffer(oldinitprob.f.sys, state_values(oldinitprob), u0idxs, u0vals) - newp = remake_buffer(oldinitprob.f.sys, parameter_values(oldinitprob), pidxs, pvals) - initprob = remake(oldinitprob; u0 = newu0, p = newp) - return initprob, odefn.update_initializeprob!, odefn.initializeprobmap, - odefn.initializeprobpmap - end - if u0 === missing || isempty(u0) - u0 = Dict() - elseif !(eltype(u0) <: Pair) - u0 = Dict(unknowns(sys) .=> u0) - end - if p === missing - p = Dict() - end - if t0 === nothing - t0 = 0.0 + if isempty(u0idxs) + newu0 = state_values(oldinitprob) + else + newu0 = remake_buffer( + oldinitprob.f.sys, state_values(oldinitprob), u0idxs, u0vals) + end + if isempty(pidxs) + newp = parameter_values(oldinitprob) + else + newp = remake_buffer( + oldinitprob.f.sys, parameter_values(oldinitprob), pidxs, pvals) + end + if oldinitprob.f.resid_prototype === nothing + newf = oldinitprob.f + else + newf = NonlinearFunction{ + SciMLBase.isinplace(oldinitprob.f), SciMLBase.specialization(oldinitprob.f)}( + oldinitprob.f; + resid_prototype = calculate_resid_prototype( + length(oldinitprob.f.resid_prototype), newu0, newp)) + end + initprob = remake(oldinitprob; f = newf, u0 = newu0, p = newp) + return SciMLBase.OverrideInitData(initprob, odefn.update_initializeprob!, + odefn.initializeprobmap, odefn.initializeprobpmap) end - u0 = todict(u0) + dvs = unknowns(sys) + ps = parameters(sys) + u0map = to_varmap(u0, dvs) + symbols_to_symbolics!(sys, u0map) + pmap = to_varmap(p, ps) + symbols_to_symbolics!(sys, pmap) + guesses = Dict() defs = defaults(sys) - varmap = merge(defs, u0) - for k in collect(keys(varmap)) - if varmap[k] === nothing - delete!(varmap, k) + if SciMLBase.has_initializeprob(odefn) + oldsys = odefn.initializeprob.f.sys + meta = get_metadata(oldsys) + if meta isa InitializationSystemMetadata + u0map = merge(meta.u0map, u0map) + pmap = merge(meta.pmap, pmap) + merge!(guesses, meta.additional_guesses) end - end - varmap = canonicalize_varmap(varmap) - missingvars = setdiff(unknowns(sys), collect(keys(varmap))) - setobserved = filter(keys(varmap)) do var - has_observed_with_lhs(sys, var) || has_observed_with_lhs(sys, default_toterm(var)) - end - p = todict(p) - guesses = ModelingToolkit.guesses(sys) - solvablepars = [par - for par in parameters(sys) - if is_parameter_solvable(par, p, defs, guesses)] - pvarmap = merge(defs, p) - setparobserved = filter(keys(pvarmap)) do var - has_parameter_dependency_with_lhs(sys, var) - end - if (((!isempty(missingvars) || !isempty(solvablepars) || - !isempty(setobserved) || !isempty(setparobserved)) && - ModelingToolkit.get_tearing_state(sys) !== nothing) || - !isempty(initialization_equations(sys))) - if SciMLBase.has_initializeprob(odefn) - oldsys = odefn.initializeprob.f.sys - meta = get_metadata(oldsys) - if meta isa InitializationSystemMetadata - u0 = merge(meta.u0map, u0) - p = merge(meta.pmap, p) + else + # there is no initializeprob, so the original problem construction + # had no solvable parameters and had the differential variables + # specified in `u0map`. + if u0 === missing + # the user didn't pass `u0` to `remake`, so they want to retain + # existing values. Fill the differential variables in `u0map`, + # initialization will either be elided or solve for the algebraic + # variables + diff_idxs = isdiffeq.(equations(sys)) + for i in eachindex(dvs) + diff_idxs[i] || continue + u0map[dvs[i]] = newu0[i] end end - for k in collect(keys(u0)) - if u0[k] === nothing - delete!(u0, k) + if p === missing + # the user didn't pass `p` to `remake`, so they want to retain + # existing values. Fill all parameters in `pmap` so that none of + # them are solvable. + for p in ps + pmap[p] = getp(sys, p)(newp) end end - for k in collect(keys(p)) - if p[k] === nothing - delete!(p, k) - end + # all non-solvable parameters need values regardless + for p in ps + haskey(pmap, p) && continue + is_parameter_solvable(p, pmap, defs, guesses) && continue + pmap[p] = getp(sys, p)(newp) end - - initprob = InitializationProblem(sys, t0, u0, p) - initprobmap = getu(initprob, unknowns(sys)) - punknowns = [p for p in all_variable_symbols(initprob) if is_parameter(sys, p)] - getpunknowns = getu(initprob, punknowns) - setpunknowns = setp(sys, punknowns) - initprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns) - reqd_syms = parameter_symbols(initprob) - update_initializeprob! = UpdateInitializeprob( - getu(sys, reqd_syms), setu(initprob, reqd_syms)) - return initprob, update_initializeprob!, initprobmap, initprobpmap - else - return nothing, nothing, nothing, nothing end + if t0 === nothing + t0 = 0.0 + end + filter_missing_values!(u0map) + filter_missing_values!(pmap) + f, _ = process_SciMLProblem(EmptySciMLFunction, sys, u0map, pmap; guesses, t = t0) + kws = f.kwargs + initprob = get(kws, :initializeprob, nothing) + if initprob === nothing + return nothing + end + return SciMLBase.OverrideInitData(initprob, get(kws, :update_initializeprob!, nothing), + get(kws, :initializeprobmap, nothing), + get(kws, :initializeprobpmap, nothing)) end """ diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index c6b2610ae7..1289388197 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -283,6 +283,16 @@ function hessian_sparsity(sys::NonlinearSystem) unknowns(sys)) for eq in equations(sys)] end +function calculate_resid_prototype(N, u0, p) + u0ElType = u0 === nothing ? Float64 : eltype(u0) + if SciMLStructures.isscimlstructure(p) + u0ElType = promote_type( + eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]), + u0ElType) + end + return zeros(u0ElType, N) +end + """ ```julia SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(sys), @@ -337,13 +347,7 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s if length(dvs) == length(equations(sys)) resid_prototype = nothing else - u0ElType = u0 === nothing ? Float64 : eltype(u0) - if SciMLStructures.isscimlstructure(p) - u0ElType = promote_type( - eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]), - u0ElType) - end - resid_prototype = zeros(u0ElType, length(equations(sys))) + resid_prototype = calculate_resid_prototype(length(equations(sys)), u0, p) end NonlinearFunction{iip}(f, @@ -690,7 +694,7 @@ A type of Nonlinear problem which specializes on polynomial systems and uses HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to create and solve. """ -struct HomotopyContinuationProblem{uType, H, D, O, SS} <: +struct HomotopyContinuationProblem{uType, H, D, O, SS, U} <: SciMLBase.AbstractNonlinearProblem{uType, true} """ The initial values of states in the system. If there are multiple real roots of @@ -721,6 +725,12 @@ struct HomotopyContinuationProblem{uType, H, D, O, SS} <: `HomotopyContinuation.solver_startsystems`. """ solver_and_starts::SS + """ + A function which takes a solution of the transformed system, and returns a vector + of solutions for the original system. This is utilized when converting systems + to polynomials. + """ + unpack_solution::U end function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...) diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl index a0fa07be47..43e9294dd3 100644 --- a/src/systems/optimization/optimizationsystem.jl +++ b/src/systems/optimization/optimizationsystem.jl @@ -101,8 +101,8 @@ function OptimizationSystem(op, unknowns, ps; gui_metadata = nothing) name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro")) - constraints = value.(scalarize(constraints)) - unknowns′ = value.(scalarize(unknowns)) + constraints = value.(reduce(vcat, scalarize(constraints); init = [])) + unknowns′ = value.(reduce(vcat, scalarize(unknowns); init = [])) ps′ = value.(ps) op′ = value(scalarize(op)) @@ -144,6 +144,34 @@ function OptimizationSystem(op, unknowns, ps; checks = checks) end +function OptimizationSystem(objective; constraints = [], kwargs...) + allunknowns = OrderedSet() + ps = OrderedSet() + collect_vars!(allunknowns, ps, objective, nothing) + for cons in constraints + collect_vars!(allunknowns, ps, cons, nothing) + end + for ssys in get(kwargs, :systems, OptimizationSystem[]) + collect_scoped_vars!(allunknowns, ps, ssys, nothing) + end + new_ps = OrderedSet() + for p in ps + if iscall(p) && operation(p) === getindex + par = arguments(p)[begin] + if Symbolics.shape(Symbolics.unwrap(par)) !== Symbolics.Unknown() && + all(par[i] in ps for i in eachindex(par)) + push!(new_ps, par) + else + push!(new_ps, p) + end + else + push!(new_ps, p) + end + end + return OptimizationSystem( + objective, collect(allunknowns), collect(new_ps); constraints, kwargs...) +end + function flatten(sys::OptimizationSystem) systems = get_systems(sys) isempty(systems) && return sys @@ -284,6 +312,7 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map, linenumbers = true, parallel = SerialForm(), eval_expression = false, eval_module = @__MODULE__, use_union = false, + checks = true, kwargs...) where {iip} if !iscomplete(sys) error("A completed `OptimizationSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `OptimizationProblem`") @@ -393,12 +422,17 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map, observedfun = ObservedFunctionCache(sys; eval_expression, eval_module) if length(cstr) > 0 - @named cons_sys = ConstraintsSystem(cstr, dvs, ps) + @named cons_sys = ConstraintsSystem(cstr, dvs, ps; checks) cons_sys = complete(cons_sys) cons, lcons_, ucons_ = generate_function(cons_sys, checkbounds = checkbounds, linenumbers = linenumbers, expression = Val{true}) - cons = eval_or_rgf.(cons; eval_expression, eval_module) + cons = let (cons_oop, cons_iip) = eval_or_rgf.(cons; eval_expression, eval_module) + _cons(u, p) = cons_oop(u, p) + _cons(resid, u, p) = cons_iip(resid, u, p) + _cons(u, p::MTKParameters) = cons_oop(u, p...) + _cons(resid, u, p::MTKParameters) = cons_iip(resid, u, p...) + end if cons_j _cons_j = let (cons_jac_oop, cons_jac_iip) = eval_or_rgf.( generate_jacobian(cons_sys; @@ -464,7 +498,7 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map, grad = _grad, hess = _hess, hess_prototype = hess_prototype, - cons = cons[2], + cons = cons, cons_j = _cons_j, cons_h = _cons_h, cons_jac_prototype = cons_jac_prototype, diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 142d2f6d71..b837eef98e 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -4,12 +4,13 @@ const AnyDict = Dict{Any, Any} $(TYPEDSIGNATURES) If called without arguments, return `Dict{Any, Any}`. Otherwise, interpret the input -as a symbolic map and turn it into a `Dict{Any, Any}`. Handles `SciMLBase.NullParameters` -and `nothing`. +as a symbolic map and turn it into a `Dict{Any, Any}`. Handles `SciMLBase.NullParameters`, +`missing` and `nothing`. """ anydict() = AnyDict() anydict(::SciMLBase.NullParameters) = AnyDict() anydict(::Nothing) = AnyDict() +anydict(::Missing) = AnyDict() anydict(x::AnyDict) = x anydict(x) = AnyDict(x) @@ -51,6 +52,42 @@ function add_toterms(varmap::AbstractDict; toterm = default_toterm) return cp end +""" + $(TYPEDSIGNATURES) + +Turn any `Symbol` keys in `varmap` to the appropriate symbolic variables in `sys`. Any +symbols that cannot be converted are ignored. +""" +function symbols_to_symbolics!(sys::AbstractSystem, varmap::AbstractDict) + if is_split(sys) + ic = get_index_cache(sys) + for k in collect(keys(varmap)) + k isa Symbol || continue + newk = get(ic.symbol_to_variable, k, nothing) + newk === nothing && continue + varmap[newk] = varmap[k] + delete!(varmap, k) + end + else + syms = all_symbols(sys) + for k in collect(keys(varmap)) + k isa Symbol || continue + idx = findfirst(syms) do sym + hasname(sym) || return false + name = getname(sym) + return name == k + end + idx === nothing && continue + newk = syms[idx] + if iscall(newk) && operation(newk) === getindex + newk = arguments(newk)[1] + end + varmap[newk] = varmap[k] + delete!(varmap, k) + end + end +end + """ $(TYPEDSIGNATURES) @@ -388,6 +425,15 @@ function evaluate_varmap!(varmap::AbstractDict, vars; limit = 100) end end +""" + $(TYPEDSIGNATURES) + +Remove keys in `varmap` whose values are `nothing`. +""" +function filter_missing_values!(varmap::AbstractDict) + filter!(kvp -> kvp[2] !== nothing, varmap) +end + struct GetUpdatedMTKParameters{G, S} # `getu` functor which gets parameters that are unknowns during initialization getpunknowns::G @@ -431,12 +477,16 @@ end $(TYPEDEF) A simple utility meant to be used as the `constructor` passed to `process_SciMLProblem` in -case constructing a SciMLFunction is not required. +case constructing a SciMLFunction is not required. The arguments passed to it are available +in the `args` field, and the keyword arguments in the `kwargs` field. """ -struct EmptySciMLFunction end +struct EmptySciMLFunction{A, K} + args::A + kwargs::K +end function EmptySciMLFunction(args...; kwargs...) - return nothing + return EmptySciMLFunction{typeof(args), typeof(kwargs)}(args, kwargs) end """ @@ -498,7 +548,7 @@ function process_SciMLProblem( constructor, sys::AbstractSystem, u0map, pmap; build_initializeprob = true, implicit_dae = false, t = nothing, guesses = AnyDict(), warn_initialize_determined = true, initialization_eqs = [], - eval_expression = false, eval_module = @__MODULE__, fully_determined = false, + eval_expression = false, eval_module = @__MODULE__, fully_determined = nothing, check_initialization_units = false, tofloat = true, use_union = false, u0_constructor = identity, du0map = nothing, check_length = true, symbolic_u0 = false, warn_cyclic_dependency = false, @@ -516,8 +566,10 @@ function process_SciMLProblem( pType = typeof(pmap) _u0map = u0map u0map = to_varmap(u0map, dvs) + symbols_to_symbolics!(sys, u0map) _pmap = pmap pmap = to_varmap(pmap, ps) + symbols_to_symbolics!(sys, pmap) defs = add_toterms(recursive_unwrap(defaults(sys))) cmap, cs = get_cmap(sys) kwargs = NamedTuple(kwargs) @@ -641,8 +693,7 @@ function process_SciMLProblem( ddvs = map(Differential(iv), dvs) du0map = to_varmap(du0map, ddvs) merge!(op, du0map) - - du0 = varmap_to_vars(du0map, ddvs; toterm = identity, + du0 = varmap_to_vars(op, ddvs; toterm = identity, tofloat = true) kwargs = merge(kwargs, (; ddvs)) else diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 848980605f..a54206d1dd 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -152,8 +152,10 @@ function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = fal noise_eqs = sorted_g_rows is_scalar_noise = false end + + noise_eqs = StructuralTransformations.tearing_substitute_expr(ode_sys, noise_eqs) return SDESystem(full_equations(ode_sys), noise_eqs, get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys); - name = nameof(ode_sys), is_scalar_noise) + name = nameof(ode_sys), is_scalar_noise, observed = observed(ode_sys)) end end diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 71bdf6fe30..1c0ad8b8ae 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -677,7 +677,11 @@ function _structural_simplify!(state::TearingState, io; simplify = false, check_consistency = true, fully_determined = true, warn_initialize_determined = false, dummy_derivative = true, kwargs...) - check_consistency &= fully_determined + if fully_determined isa Bool + check_consistency &= fully_determined + else + check_consistency = true + end has_io = io !== nothing orig_inputs = Set() if has_io @@ -690,7 +694,8 @@ function _structural_simplify!(state::TearingState, io; simplify = false, end sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...) if check_consistency - ModelingToolkit.check_consistency(state, orig_inputs) + fully_determined = ModelingToolkit.check_consistency( + state, orig_inputs; nothrow = fully_determined === nothing) end if fully_determined && dummy_derivative sys = ModelingToolkit.dummy_derivative( diff --git a/src/utils.jl b/src/utils.jl index 830ec98e44..416efd8f2c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -516,6 +516,15 @@ function collect_scoped_vars!(unknowns, parameters, sys, iv; depth = 1, op = Dif end end end + if has_constraints(sys) + for eq in get_constraints(sys) + eqtype_supports_collect_vars(eq) || continue + collect_vars!(unknowns, parameters, eq, iv; depth, op) + end + end + if has_op(sys) + collect_vars!(unknowns, parameters, get_op(sys), iv; depth, op) + end newdepth = depth == -1 ? depth : depth + 1 for ssys in get_systems(sys) collect_scoped_vars!(unknowns, parameters, ssys, iv; depth = newdepth, op) @@ -544,9 +553,10 @@ Can be dispatched by higher-level libraries to indicate support. """ eqtype_supports_collect_vars(eq) = false eqtype_supports_collect_vars(eq::Equation) = true +eqtype_supports_collect_vars(eq::Inequality) = true eqtype_supports_collect_vars(eq::Pair) = true -function collect_vars!(unknowns, parameters, eq::Equation, iv; +function collect_vars!(unknowns, parameters, eq::Union{Equation, Inequality}, iv; depth = 0, op = Differential) collect_vars!(unknowns, parameters, eq.lhs, iv; depth, op) collect_vars!(unknowns, parameters, eq.rhs, iv; depth, op) diff --git a/src/variables.jl b/src/variables.jl index fcfc7f9b1e..d11c6c1834 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -254,7 +254,6 @@ end ## Bounds ====================================================================== struct VariableBounds end Symbolics.option_to_metadata_type(::Val{:bounds}) = VariableBounds -getbounds(x::Num) = getbounds(Symbolics.unwrap(x)) """ getbounds(x) @@ -266,10 +265,35 @@ Create parameters with bounds like this @parameters p [bounds=(-1, 1)] ``` """ -function getbounds(x) +function getbounds(x::Union{Num, Symbolics.Arr, SymbolicUtils.Symbolic}) + x = unwrap(x) p = Symbolics.getparent(x, nothing) - p === nothing || (x = p) - Symbolics.getmetadata(x, VariableBounds, (-Inf, Inf)) + if p === nothing + bounds = Symbolics.getmetadata(x, VariableBounds, (-Inf, Inf)) + if symbolic_type(x) == ArraySymbolic() && Symbolics.shape(x) != Symbolics.Unknown() + bounds = map(bounds) do b + b isa AbstractArray && return b + return fill(b, size(x)) + end + end + else + # if we reached here, `x` is the result of calling `getindex` + bounds = @something Symbolics.getmetadata(x, VariableBounds, nothing) getbounds(p) + idxs = arguments(x)[2:end] + bounds = map(bounds) do b + if b isa AbstractArray + if Symbolics.shape(p) != Symbolics.Unknown() && size(p) != size(b) + throw(DimensionMismatch("Expected array variable $p with shape $(size(p)) to have bounds of identical size. Found $bounds of size $(size(bounds)).")) + end + return b[idxs...] + elseif symbolic_type(x) == ArraySymbolic() + return fill(b, size(x)) + else + return b + end + end + end + return bounds end """ @@ -280,7 +304,7 @@ See also [`getbounds`](@ref). """ function hasbounds(x) b = getbounds(x) - isfinite(b[1]) || isfinite(b[2]) + any(isfinite.(b[1]) .|| isfinite.(b[2])) end ## Disturbance ================================================================= diff --git a/test/components.jl b/test/components.jl index a284c9b0c8..298f9fceb9 100644 --- a/test/components.jl +++ b/test/components.jl @@ -350,3 +350,23 @@ end @test_throws ArgumentError simp.inner₊p @test_throws ArgumentError outer.inner₊p end + +@testset "`getproperty` on `structural_simplify(complete(sys))`" begin + @mtkmodel Foo begin + @variables begin + x(t) + end + end + @mtkmodel Bar begin + @components begin + foo = Foo() + end + @equations begin + D(foo.x) ~ foo.x + end + end + @named bar = Bar() + cbar = complete(bar) + ss = structural_simplify(cbar) + @test isequal(cbar.foo.x, ss.foo.x) +end diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index d11b88a6e9..5f7afe222a 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -4,8 +4,12 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ChainRulesTestUtils = "cdddcdb0-9152-4a09-a978-84456f9df70a" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" +InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 172ee58890..81c252c84a 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -75,19 +75,59 @@ end @testset "Polynomial check and warnings" begin @variables x = 1.0 @mtkbuild sys = NonlinearSystem([x^1.5 + x^2 - 1 ~ 0]) - @test_throws ["Exponent", "not an integer", "not a polynomial"] HomotopyContinuationProblem( + @test_throws ["Cannot convert", "Unable", "symbolically solve", + "Exponent", "not an integer", "not a polynomial"] HomotopyContinuationProblem( sys, []) @mtkbuild sys = NonlinearSystem([x^x - x ~ 0]) - @test_throws ["Exponent", "unknowns", "not a polynomial"] HomotopyContinuationProblem( + @test_throws ["Cannot convert", "Unable", "symbolically solve", + "Exponent", "unknowns", "not a polynomial"] HomotopyContinuationProblem( sys, []) @mtkbuild sys = NonlinearSystem([((x^2) / sin(x))^2 + x ~ 0]) - @test_throws ["recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( + @test_throws ["Cannot convert", "both polynomial", "non-polynomial", + "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( sys, []) @variables y = 2.0 @mtkbuild sys = NonlinearSystem([x^2 + y^2 + 2 ~ 0, y ~ sin(x)]) - @test_throws ["recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( + @test_throws ["Cannot convert", "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( sys, []) + + @mtkbuild sys = NonlinearSystem([x^2 + y^2 - 2 ~ 0, sin(x + y) ~ 0]) + @test_throws ["Cannot convert", "function of multiple unknowns"] HomotopyContinuationProblem( + sys, []) + + @mtkbuild sys = NonlinearSystem([sin(x)^2 + 1 ~ 0, cos(y) - cos(x) - 1 ~ 0]) + @test_throws ["Cannot convert", "multiple non-polynomial terms", "same unknown"] HomotopyContinuationProblem( + sys, []) + + @mtkbuild sys = NonlinearSystem([sin(x^2)^2 + sin(x^2) - 1 ~ 0]) + @test_throws ["import Nemo"] HomotopyContinuationProblem(sys, []) +end + +import Nemo + +@testset "With Nemo" begin + @variables x = 2.0 + @mtkbuild sys = NonlinearSystem([sin(x^2)^2 + sin(x^2) - 1 ~ 0]) + prob = HomotopyContinuationProblem(sys, []) + @test prob[1] ≈ 2.0 + sol = solve(prob; threading = false) + _x = sol[1] + @test sin(_x^2)^2 + sin(_x^2) - 1≈0.0 atol=1e-12 +end + +@testset "Function of polynomial" begin + @variables x=0.25 y=0.125 + a = sin(x^2 - 4x + 1) + b = cos(3log(y) + 4) + @mtkbuild sys = NonlinearSystem([(a^2 - 4a * b + 4b^2) / (a - 0.25) ~ 0 + (a^2 - 0.75a + 0.125) ~ 0]) + prob = HomotopyContinuationProblem(sys, []) + @test prob[x] ≈ 0.25 + @test prob[y] ≈ 0.125 + sol = solve(prob; threading = false) + @test sol[a]≈0.5 atol=1e-6 + @test sol[b]≈0.25 atol=1e-6 end @testset "Rational functions" begin diff --git a/test/extensions/test_infiniteopt.jl b/test/extensions/test_infiniteopt.jl new file mode 100644 index 0000000000..e45aa0f2fd --- /dev/null +++ b/test/extensions/test_infiniteopt.jl @@ -0,0 +1,102 @@ +using ModelingToolkit, InfiniteOpt, JuMP, Ipopt +using ModelingToolkit: D_nounits as D, t_nounits as t, varmap_to_vars + +@mtkmodel Pendulum begin + @parameters begin + g = 9.8 + L = 0.4 + K = 1.2 + m = 0.3 + end + @variables begin + θ(t) # state + ω(t) # state + τ(t) = 0 # input + y(t) # output + end + @equations begin + D(θ) ~ ω + D(ω) ~ -g / L * sin(θ) - K / m * ω + τ / m / L^2 + y ~ θ * 180 / π + end +end +@named model = Pendulum() +model = complete(model) + +inputs = [model.τ] +(f_oop, f_ip), dvs, psym, io_sys = ModelingToolkit.generate_control_function( + model, inputs, split = false) + +outputs = [model.y] +f_obs = ModelingToolkit.build_explicit_observed_function(io_sys, outputs; inputs = inputs) + +expected_state_order = [model.θ, model.ω] +permutation = [findfirst(isequal(x), expected_state_order) for x in dvs] # This maps our expected state order to the actual state order + +## + +ub = varmap_to_vars([model.θ => 2pi, model.ω => 10], dvs) +lb = varmap_to_vars([model.θ => -2pi, model.ω => -10], dvs) +xf = varmap_to_vars([model.θ => pi, model.ω => 0], dvs) +nx = length(dvs) +nu = length(inputs) +ny = length(outputs) + +## +m = InfiniteModel(optimizer_with_attributes(Ipopt.Optimizer, + "print_level" => 0, "acceptable_tol" => 1e-3, "constr_viol_tol" => 1e-5, "max_iter" => 1000, + "tol" => 1e-5, "mu_strategy" => "monotone", "nlp_scaling_method" => "gradient-based", + "alpha_for_y" => "safer-min-dual-infeas", "bound_mult_init_method" => "mu-based", "print_user_options" => "yes")); + +@infinite_parameter(m, τ in [0, 1], num_supports=51, + derivative_method=OrthogonalCollocation(4)) # Time variable +guess_xs = [t -> pi, t -> 0.1][permutation] +guess_us = [t -> 0.1] +InfiniteOpt.@variables(m, + begin + # state variables + (lb[i] <= x[i = 1:nx] <= ub[i], Infinite(τ), start = guess_xs[i]) # state variables + -10 <= u[i = 1:nu] <= 10, Infinite(τ), (start = guess_us[i]) # control variables + 0 <= tf <= 10, (start = 5) # Final time + 0.2 <= L <= 0.6, (start = 0.4) # Length parameter + end) + +# Trace the dynamics +x0, p = ModelingToolkit.get_u0_p(io_sys, [model.θ => 0, model.ω => 0], [model.L => L]) + +xp = f_oop(x, u, p, τ) +cp = f_obs(x, u, p, τ) # Test that it's possible to trace through an observed function + +@objective(m, Min, tf) +@constraint(m, [i = 1:nx], x[i](0)==x0[i]) # Initial condition +@constraint(m, [i = 1:nx], x[i](1)==xf[i]) # Terminal state + +x_scale = varmap_to_vars([model.θ => 1 + model.ω => 1], dvs) + +# Add dynamics constraints +@constraint(m, [i = 1:nx], (∂(x[i], τ) - tf * xp[i]) / x_scale[i]==0) + +optimize!(m) + +# Extract the optimal solution +opt_tf = value(tf) +opt_time = opt_tf * value(τ) +opt_x = [value(x[i]) for i in permutation] +opt_u = [value(u[i]) for i in 1:nu] +opt_L = value(L) + +# Plot the results +# using Plots +# plot(opt_time, opt_x[1], label = "θ", xlabel = "Time [s]", layout=3) +# plot!(opt_time, opt_x[2], label = "ω", sp=2) +# plot!(opt_time, opt_u[1], label = "τ", sp=3) + +using Test +@test opt_x[1][end]≈pi atol=1e-3 +@test opt_x[2][end]≈0 atol=1e-3 + +@test opt_x[1][1]≈0 atol=1e-3 +@test opt_x[2][1]≈0 atol=1e-3 + +@test opt_L≈0.2 atol=1e-3 # Smallest permissible length is optimal diff --git a/test/index_cache.jl b/test/index_cache.jl index c479075797..455203d759 100644 --- a/test/index_cache.jl +++ b/test/index_cache.jl @@ -92,3 +92,29 @@ end reorder_dimension_by_tunables!(dst, sys, src, [r, q, p]; dim = 2) @test dst ≈ stack([vcat(4ones(4), 3ones(3), 1.0) for i in 1:5]; dims = 1) end + +mutable struct ParamTest + y::Any +end +(pt::ParamTest)(x) = pt.y - x +@testset "Issue#3215: Callable discrete parameter" begin + function update_affect!(integ, u, p, ctx) + integ.p[p.p_1].y = integ.t + end + + tp1 = typeof(ParamTest(1)) + @parameters (p_1::tp1)(..) = ParamTest(1) + @variables x(ModelingToolkit.t_nounits) = 0 + + event1 = [1.0, 2, 3] => (update_affect!, [], [p_1], [p_1], nothing) + + @named sys = ODESystem([ + ModelingToolkit.D_nounits(x) ~ p_1(x) + ], + ModelingToolkit.t_nounits; + discrete_events = [event1] + ) + ss = @test_nowarn complete(sys) + @test length(parameters(ss)) == 1 + @test !is_timeseries_parameter(ss, p_1) +end diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index c9783c9659..f3015f7db0 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -877,3 +877,175 @@ end sol = solve(prob, Rodas5P()) @test SciMLBase.successful_retcode(sol) end + +@testset "Issue#3205" begin + using ModelingToolkitStandardLibrary.Electrical + import ModelingToolkitStandardLibrary.Mechanical.Rotational as MR + using ModelingToolkitStandardLibrary.Blocks + using SciMLBase + + function dc_motor(R1 = 0.5) + R = R1 # [Ohm] armature resistance + L = 4.5e-3 # [H] armature inductance + k = 0.5 # [N.m/A] motor constant + J = 0.02 # [kg.m²] inertia + f = 0.01 # [N.m.s/rad] friction factor + tau_L_step = -0.3 # [N.m] amplitude of the load torque step + + @named ground = Ground() + @named source = Voltage() + @named ref = Blocks.Step(height = 0.2, start_time = 0) + @named pi_controller = Blocks.LimPI(k = 1.1, T = 0.035, u_max = 10, Ta = 0.035) + @named feedback = Blocks.Feedback() + @named R1 = Resistor(R = R) + @named L1 = Inductor(L = L) + @named emf = EMF(k = k) + @named fixed = MR.Fixed() + @named load = MR.Torque() + @named load_step = Blocks.Step(height = tau_L_step, start_time = 3) + @named inertia = MR.Inertia(J = J) + @named friction = MR.Damper(d = f) + @named speed_sensor = MR.SpeedSensor() + + connections = [connect(fixed.flange, emf.support, friction.flange_b) + connect(emf.flange, friction.flange_a, inertia.flange_a) + connect(inertia.flange_b, load.flange) + connect(inertia.flange_b, speed_sensor.flange) + connect(load_step.output, load.tau) + connect(ref.output, feedback.input1) + connect(speed_sensor.w, :y, feedback.input2) + connect(feedback.output, pi_controller.err_input) + connect(pi_controller.ctr_output, :u, source.V) + connect(source.p, R1.p) + connect(R1.n, L1.p) + connect(L1.n, emf.p) + connect(emf.n, source.n, ground.g)] + + @named model = ODESystem(connections, t, + systems = [ + ground, + ref, + pi_controller, + feedback, + source, + R1, + L1, + emf, + fixed, + load, + load_step, + inertia, + friction, + speed_sensor + ]) + end + + model = dc_motor() + sys = structural_simplify(model) + + prob = ODEProblem(sys, [sys.L1.i => 0.0], (0, 6.0)) + + @test_nowarn remake(prob, p = prob.p) +end + +@testset "Singular initialization prints a warning" begin + @parameters g + @variables x(t) y(t) [state_priority = 10] λ(t) + eqs = [D(D(x)) ~ λ * x + D(D(y)) ~ λ * y - g + x^2 + y^2 ~ 1] + @mtkbuild pend = ODESystem(eqs, t) + @test_warn ["structurally singular", "initialization", "Guess", "heuristic"] ODEProblem( + pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], 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] + @mtkbuild sys = ODESystem( + [D(x) ~ p * y + q * t, 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]) + + integ = init(prob, DImplicitEuler()) + @test integ[x] ≈ 1.0 + @test integ[y]≈cbrt(4) rtol=1e-6 + @test integ.ps[p] ≈ 1.0 + @test integ.ps[q]≈cbrt(2) rtol=1e-6 +end + +@testset "Guesses provided to `ODEProblem` are used in `remake`" begin + @variables x(t) y(t)=2x + @parameters p q=3x + @mtkbuild sys = ODESystem([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]) + @test prob[x] == 0.0 + @test prob[y] == 0.0 + @test prob.ps[p] == 1.0 + @test prob.ps[q] == 0.0 + integ = init(prob) + @test integ[x] ≈ 1 / cbrt(3) + @test integ[y] ≈ 2 / cbrt(3) + @test integ.ps[p] == 1.0 + @test integ.ps[q] ≈ 3 / cbrt(3) + prob2 = remake(prob; u0 = [y => 3x], p = [q => 2x]) + integ2 = init(prob2) + @test integ2[x] ≈ cbrt(3 / 28) + @test integ2[y] ≈ 3cbrt(3 / 28) + @test integ2.ps[p] == 1.0 + @test integ2.ps[q] ≈ 2cbrt(3 / 28) +end + +@testset "Remake problem with no initializeprob" begin + @variables x(t) [guess = 1.0] y(t) [guess = 1.0] + @parameters p [guess = 1.0] q [guess = 1.0] + @mtkbuild sys = ODESystem( + [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]) + @test prob.f.initialization_data === nothing + prob2 = remake(prob; u0 = [x => 2.0]) + @test prob2[x] == 2.0 + @test prob2.f.initialization_data === nothing + prob3 = remake(prob; u0 = [y => 2.0]) + @test prob3.f.initialization_data !== nothing + @test init(prob3)[x] ≈ 1.0 + prob4 = remake(prob; p = [p => 1.0]) + @test prob4.f.initialization_data === nothing + prob5 = remake(prob; p = [p => missing, q => 2.0]) + @test prob5.f.initialization_data !== nothing + @test init(prob5).ps[p] ≈ 1.0 +end + +@testset "Variables provided as symbols" begin + @variables x(t) [guess = 1.0] y(t) [guess = 1.0] + @parameters p [guess = 1.0] q [guess = 1.0] + @mtkbuild sys = ODESystem( + [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]) + @test prob.f.initialization_data === nothing + prob2 = remake(prob; u0 = [:x => 2.0]) + @test prob2.f.initialization_data === nothing + prob3 = remake(prob; u0 = [:y => 1.0]) + @test prob3.f.initialization_data !== nothing + @test init(prob3)[x] ≈ 0.5 +end + +@testset "Issue#3246: type promotion with parameter dependent initialization_eqs" begin + @variables x(t)=1 y(t)=1 + @parameters a = 1 + @named sys = ODESystem([D(x) ~ 0, D(y) ~ x + a], t; initialization_eqs = [y ~ a]) + + ssys = structural_simplify(sys) + prob = ODEProblem(ssys, [], (0, 1), []) + + @test SciMLBase.successful_retcode(solve(prob)) + + seta = setsym_oop(prob, [a]) + (newu0, newp) = seta(prob, ForwardDiff.Dual{ForwardDiff.Tag{:tag, Float64}}.([1.0], 1)) + newprob = remake(prob, u0 = newu0, p = newp) + + @test SciMLBase.successful_retcode(solve(newprob)) +end diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl index 07a26d2d1a..9550a87f31 100644 --- a/test/input_output_handling.jl +++ b/test/input_output_handling.jl @@ -140,8 +140,7 @@ if VERSION >= v"1.8" # :opaque_closure not supported before A, B, C, D = matrices obsf = ModelingToolkit.build_explicit_observed_function(ssys, [y], - inputs = [torque.tau.u], - drop_expr = identity) + inputs = [torque.tau.u]) x = randn(size(A, 1)) u = randn(size(B, 2)) p = (getindex.( @@ -161,12 +160,12 @@ eqs = [ ] @named sys = ODESystem(eqs, t) -f, dvs, ps = ModelingToolkit.generate_control_function(sys, simplify = true) +f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys, simplify = true) @test isequal(dvs[], x) @test isempty(ps) -p = [] +p = nothing x = [rand()] u = [rand()] @test f[1](x, u, p, 1) == -x + u @@ -222,10 +221,10 @@ eqs = [connect_sd(sd, mass1, mass2) @named _model = ODESystem(eqs, t) @named model = compose(_model, mass1, mass2, sd); -f, dvs, ps = ModelingToolkit.generate_control_function(model, simplify = true) +f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(model, simplify = true) @test length(dvs) == 4 @test length(ps) == length(parameters(model)) -p = ModelingToolkit.varmap_to_vars(ModelingToolkit.defaults(model), ps) +p = MTKParameters(io_sys, [io_sys.u => NaN]) x = ModelingToolkit.varmap_to_vars( merge(ModelingToolkit.defaults(model), Dict(D.(unknowns(model)) .=> 0.0)), dvs) @@ -289,7 +288,7 @@ model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.i @named dmodel = Blocks.StateSpace([0.0], [1.0], [1.0], [0.0]) # An integrating disturbance @named dist = ModelingToolkit.DisturbanceModel(model.torque.tau.u, dmodel) -(f_oop, f_ip), outersys, dvs, p = ModelingToolkit.add_input_disturbance(model, dist) +(f_oop, f_ip), outersys, dvs, p, io_sys = ModelingToolkit.add_input_disturbance(model, dist) @unpack u, d = outersys matrices, ssys = linearize(outersys, [u, d], model_outputs) @@ -303,7 +302,7 @@ x_add = ModelingToolkit.varmap_to_vars(merge(Dict(dvs .=> 0), Dict(dstate => 1)) x0 = randn(5) x1 = copy(x0) + x_add # add disturbance state perturbation u = randn(1) -pn = ModelingToolkit.varmap_to_vars(def, p) +pn = MTKParameters(io_sys, []) xp0 = f_oop(x0, u, pn, 0) xp1 = f_oop(x1, u, pn, 0) @@ -402,3 +401,15 @@ end f, dvs, ps = ModelingToolkit.generate_control_function(sys, simplify = true) @test f[1]([0.5], nothing, nothing, 0.0) == [1.0] end + +@testset "With callable symbolic" begin + @variables x(t)=0 u(t)=0 [input = true] + @parameters p(::Real) = (x -> 2x) + eqs = [D(x) ~ -x + p(u)] + @named sys = ODESystem(eqs, t) + f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys, simplify = true) + p = MTKParameters(io_sys, []) + u = [1.0] + x = [1.0] + @test_nowarn f[1](x, u, p, 0.0) +end diff --git a/test/jumpsystem.jl b/test/jumpsystem.jl index 30cb636226..0fd6dc4af0 100644 --- a/test/jumpsystem.jl +++ b/test/jumpsystem.jl @@ -1,5 +1,5 @@ using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra -using Random, StableRNGs +using Random, StableRNGs, NonlinearSolve using OrdinaryDiffEq using ModelingToolkit: t_nounits as t, D_nounits as D MT = ModelingToolkit @@ -422,3 +422,114 @@ let @test issetequal(us, [x5]) @test issetequal(ps, [p5]) end + +# PDMP test +let + seed = 1111 + Random.seed!(rng, seed) + @variables X(t) Y(t) + @parameters k1 k2 + vrj1 = VariableRateJump(k1 * X, [X ~ X - 1]; save_positions = (false, false)) + vrj2 = VariableRateJump(k1, [Y ~ Y + 1]; save_positions = (false, false)) + eqs = [D(X) ~ k2, D(Y) ~ -k2 / 10 * Y] + @named jsys = JumpSystem([vrj1, vrj2, eqs[1], eqs[2]], t, [X, Y], [k1, k2]) + jsys = complete(jsys) + X0 = 0.0 + Y0 = 3.0 + u0 = [X => X0, Y => Y0] + k1val = 1.0 + k2val = 20.0 + p = [k1 => k1val, k2 => k2val] + tspan = (0.0, 10.0) + oprob = ODEProblem(jsys, u0, tspan, p) + jprob = JumpProblem(jsys, oprob; rng, save_positions = (false, false)) + + times = range(0.0, tspan[2], length = 100) + Nsims = 4000 + Xv = zeros(length(times)) + Yv = zeros(length(times)) + for n in 1:Nsims + sol = solve(jprob, Tsit5(); saveat = times, seed) + Xv .+= sol[1, :] + Yv .+= sol[2, :] + seed += 1 + end + Xv ./= Nsims + Yv ./= Nsims + + Xact(t) = X0 * exp(-k1val * t) + (k2val / k1val) * (1 - exp(-k1val * t)) + function Yact(t) + Y0 * exp(-k2val / 10 * t) + (k1val / (k2val / 10)) * (1 - exp(-k2val / 10 * t)) + end + @test all(abs.(Xv .- Xact.(times)) .<= 0.05 .* Xv) + @test all(abs.(Yv .- Yact.(times)) .<= 0.1 .* Yv) +end + +# that mixes ODEs and jump types, and then contin events +let + seed = 1111 + Random.seed!(rng, seed) + @variables X(t) Y(t) + @parameters α β + vrj = VariableRateJump(β * X, [X ~ X - 1]; save_positions = (false, false)) + crj = ConstantRateJump(β * Y, [Y ~ Y - 1]) + maj = MassActionJump(α, [0 => 1], [Y => 1]) + eqs = [D(X) ~ α * (1 + Y)] + @named jsys = JumpSystem([maj, crj, vrj, eqs[1]], t, [X, Y], [α, β]) + jsys = complete(jsys) + p = (α = 6.0, β = 2.0, X₀ = 2.0, Y₀ = 1.0) + u0map = [X => p.X₀, Y => p.Y₀] + pmap = [α => p.α, β => p.β] + tspan = (0.0, 20.0) + oprob = ODEProblem(jsys, u0map, tspan, pmap) + jprob = JumpProblem(jsys, oprob; rng, save_positions = (false, false)) + times = range(0.0, tspan[2], length = 100) + Nsims = 4000 + Xv = zeros(length(times)) + Yv = zeros(length(times)) + for n in 1:Nsims + sol = solve(jprob, Tsit5(); saveat = times, seed) + Xv .+= sol[1, :] + Yv .+= sol[2, :] + seed += 1 + end + Xv ./= Nsims + Yv ./= Nsims + + function Yf(t, p) + local α, β, X₀, Y₀ = p + return (α / β) + (Y₀ - α / β) * exp(-β * t) + end + function Xf(t, p) + local α, β, X₀, Y₀ = p + return (α / β) + (α^2 / β^2) + α * (Y₀ - α / β) * t * exp(-β * t) + + (X₀ - α / β - α^2 / β^2) * exp(-β * t) + end + Xact = [Xf(t, p) for t in times] + Yact = [Yf(t, p) for t in times] + @test all(abs.(Xv .- Xact) .<= 0.05 .* Xv) + @test all(abs.(Yv .- Yact) .<= 0.05 .* Yv) + + function affect!(integ, u, p, ctx) + savevalues!(integ, true) + terminate!(integ) + nothing + end + cevents = [t ~ 0.2] => (affect!, [], [], [], nothing) + @named jsys = JumpSystem([maj, crj, vrj, eqs[1]], t, [X, Y], [α, β]; + continuous_events = cevents) + jsys = complete(jsys) + tspan = (0.0, 200.0) + oprob = ODEProblem(jsys, u0map, tspan, pmap) + jprob = JumpProblem(jsys, oprob; rng, save_positions = (false, false)) + Xsamp = 0.0 + Nsims = 4000 + for n in 1:Nsims + sol = solve(jprob, Tsit5(); saveat = tspan[2], seed) + @test sol.retcode == ReturnCode.Terminated + Xsamp += sol[1, end] + seed += 1 + end + Xsamp /= Nsims + @test abs(Xsamp - Xf(0.2, p) < 0.05 * Xf(0.2, p)) +end diff --git a/test/model_parsing.jl b/test/model_parsing.jl index fd223800e3..6ab7636b32 100644 --- a/test/model_parsing.jl +++ b/test/model_parsing.jl @@ -979,3 +979,14 @@ end @test MultipleExtend.structure[:extend][1] == [:inmodel, :b, :inmodel_b] @test tosymbol.(parameters(multiple_extend)) == [:b, :inmodel_b₊p, :inmodel₊p] end + +struct CustomStruct end +@testset "Nonnumeric parameters" begin + @mtkmodel MyModel begin + @parameters begin + p::CustomStruct + end + end + @named sys = MyModel(p = CustomStruct()) + @test ModelingToolkit.defaults(sys)[@nonamespace sys.p] == CustomStruct() +end diff --git a/test/modelingtoolkitize.jl b/test/modelingtoolkitize.jl index 621bf530b7..3bcfdecc8b 100644 --- a/test/modelingtoolkitize.jl +++ b/test/modelingtoolkitize.jl @@ -441,3 +441,31 @@ prob = NonlinearLeastSquaresProblem( sys = modelingtoolkitize(prob) @test length(equations(sys)) == 3 @test length(equations(structural_simplify(sys; fully_determined = false))) == 0 + +@testset "`modelingtoolkitize(::SDEProblem)` sets defaults" begin + function sdeg!(du, u, p, t) + du[1] = 0.3 * u[1] + du[2] = 0.3 * u[2] + du[3] = 0.3 * u[3] + end + function sdef!(du, u, p, t) + x, y, z = u + sigma, rho, beta = p + du[1] = sigma * (y - x) + du[2] = x * (rho - z) - y + du[3] = x * y - beta * z + end + u0 = [1.0, 0.0, 0.0] + tspan = (0.0, 100.0) + p = [10.0, 28.0, 2.66] + sprob = SDEProblem(sdef!, sdeg!, u0, tspan, p) + sys = complete(modelingtoolkitize(sprob)) + @test length(ModelingToolkit.defaults(sys)) == 6 + sprob2 = SDEProblem(sys, [], tspan) + + truevals = similar(u0) + sprob.f(truevals, u0, p, tspan[1]) + mtkvals = similar(u0) + sprob2.f(mtkvals, sprob2.u0, sprob2.p, tspan[1]) + @test mtkvals ≈ truevals +end diff --git a/test/odesystem.jl b/test/odesystem.jl index d7a2658f63..85d135b338 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1504,3 +1504,51 @@ end sys2 = complete(sys; split = false) @test ModelingToolkit.get_index_cache(sys2) === nothing end + +# https://github.com/SciML/SciMLBase.jl/issues/786 +@testset "Observed variables dependent on discrete parameters" begin + @variables x(t) obs(t) + @parameters c(t) + @mtkbuild sys = ODESystem( + [D(x) ~ c * cos(x), obs ~ c], t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]]) + prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) + sol = solve(prob, Tsit5()) + @test sol[obs] ≈ 1:7 +end + +@testset "DAEProblem with array parameters" begin + @variables x(t)=1.0 y(t) [guess = 1.0] + @parameters p[1:2] = [1.0, 2.0] + @mtkbuild sys = ODESystem([D(x) ~ x, y^2 ~ x + sum(p)], t) + 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 + +@testset "Symbolic tstops" begin + @variables x(t) = 1.0 + @parameters p=0.15 q=0.25 r[1:2]=[0.35, 0.45] + @mtkbuild sys = ODESystem( + [D(x) ~ p * x + q * t + sum(r)], t; tstops = [0.5p, [0.1, 0.2], [p + 2q], r]) + prob = ODEProblem(sys, [], (0.0, 5.0)) + sol = solve(prob) + 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) + prob2 = remake(prob; tspan = (0.0, 10.0)) + sol2 = solve(prob2) + expected_tstops = unique!(sort!(vcat(0.0:0.075:10.0, 0.1, 0.2, 0.65, 0.35, 0.45))) + @test all(x -> any(isapprox(x, atol = 1e-6), sol2.t), expected_tstops) + + @variables y(t) [guess = 1.0] + @mtkbuild sys = ODESystem([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)) + 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) + prob2 = remake(prob; tspan = (0.0, 10.0)) + sol2 = solve(prob2, DImplicitEuler()) + expected_tstops = unique!(sort!(vcat(0.0:0.075:10.0, 0.1, 0.2, 0.65, 0.35, 0.45))) + @test all(x -> any(isapprox(x, atol = 1e-6), sol2.t), expected_tstops) +end diff --git a/test/optimizationsystem.jl b/test/optimizationsystem.jl index dfa11fca37..c20613441a 100644 --- a/test/optimizationsystem.jl +++ b/test/optimizationsystem.jl @@ -358,4 +358,33 @@ end @mtkbuild sys = OptimizationSystem(obj, [x, y, z], []; constraints = cons) @test is_variable(sys, z) @test !is_variable(sys, y) + + @variables x[1:3] [bounds = ([-Inf, -1.0, -2.0], [Inf, 1.0, 2.0])] + obj = x[1]^2 + x[2]^2 + x[3]^2 + cons = [x[2] ~ 2x[1] + 3, x[3] ~ x[1] + x[2]] + @mtkbuild sys = OptimizationSystem(obj, [x], []; constraints = cons) + @test length(unknowns(sys)) == 2 + @test !is_variable(sys, x[1]) + @test is_variable(sys, x[2]) + @test is_variable(sys, x[3]) +end + +@testset "Constraints work with nonnumeric parameters" begin + @variables x + @parameters p f(::Real) + @mtkbuild 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)]) + @test abs(prob.f.cons(prob.u0, prob.p)[1]) ≈ 1.0 +end + +@testset "Variable discovery" begin + @variables x1 x2 + @parameters p1 p2 + @named sys1 = OptimizationSystem(x1^2; constraints = [p1 * x1 ≲ 2.0]) + @named sys2 = OptimizationSystem(x2^2; constraints = [p2 * x2 ≲ 2.0], systems = [sys1]) + @test isequal(only(unknowns(sys1)), x1) + @test isequal(only(parameters(sys1)), p1) + @test all(y -> any(x -> isequal(x, y), unknowns(sys2)), [x2, sys1.x1]) + @test all(y -> any(x -> isequal(x, y), parameters(sys2)), [p2, sys1.p1]) end diff --git a/test/runtests.jl b/test/runtests.jl index 192ead935c..44846eed57 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,9 +33,6 @@ end @safetestset "Dynamic Quantities Test" include("dq_units.jl") @safetestset "Unitful Quantities Test" include("units.jl") @safetestset "Mass Matrix Test" include("mass_matrix.jl") - @safetestset "InitializationSystem Test" include("initializationsystem.jl") - @safetestset "Guess Propagation" include("guess_propagation.jl") - @safetestset "Hierarchical Initialization Equations" include("hierarchical_initialization_eqs.jl") @safetestset "Reduction Test" include("reduction.jl") @safetestset "Split Parameters Test" include("split_parameters.jl") @safetestset "StaticArrays Test" include("static_arrays.jl") @@ -58,12 +55,18 @@ end @safetestset "Constants Test" include("constants.jl") @safetestset "Parameter Dependency Test" include("parameter_dependencies.jl") @safetestset "Generate Custom Function Test" include("generate_custom_function.jl") - @safetestset "Initial Values Test" include("initial_values.jl") @safetestset "Equation Type Accessors Test" include("equation_type_accessors.jl") @safetestset "Equations with complex values" include("complex.jl") end end + if GROUP == "All" || GROUP == "Initialization" + @safetestset "Guess Propagation" include("guess_propagation.jl") + @safetestset "Hierarchical Initialization Equations" include("hierarchical_initialization_eqs.jl") + @safetestset "InitializationSystem Test" include("initializationsystem.jl") + @safetestset "Initial Values Test" include("initial_values.jl") + end + if GROUP == "All" || GROUP == "InterfaceII" @testset "InterfaceII" begin @safetestset "IndexCache Test" include("index_cache.jl") @@ -113,5 +116,6 @@ end @safetestset "Auto Differentiation Test" include("extensions/ad.jl") @safetestset "LabelledArrays Test" include("labelledarrays.jl") @safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl") + @safetestset "InfiniteOpt Extension Test" include("extensions/test_infiniteopt.jl") end end diff --git a/test/sdesystem.jl b/test/sdesystem.jl index c258a4142b..749aca86a7 100644 --- a/test/sdesystem.jl +++ b/test/sdesystem.jl @@ -780,3 +780,32 @@ end prob = @test_nowarn SDEProblem(sys, nothing, (0.0, 1.0)) @test_nowarn solve(prob, ImplicitEM()) end + +@testset "Issue#3212: Noise dependent on observed" begin + sts = @variables begin + x(t) = 1.0 + input(t) + [input = true] + end + ps = @parameters a = 2 + @brownian η + + eqs = [D(x) ~ -a * x + (input + 1) * η + input ~ 0.0] + + sys = System(eqs, t, sts, ps; name = :name) + sys = structural_simplify(sys) + @test ModelingToolkit.get_noiseeqs(sys) ≈ [1.0] + prob = SDEProblem(sys, [], (0.0, 1.0), []) + @test_nowarn solve(prob, RKMil()) +end + +@testset "Observed variables retained after `structural_simplify`" begin + @variables x(t) y(t) z(t) + @brownian a + @mtkbuild sys = System([D(x) ~ x + a, D(y) ~ y + a, z ~ x + y], t) + @test sys isa SDESystem + @test length(observed(sys)) == 1 + prob = SDEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0)) + @test prob[z] ≈ 2.0 +end diff --git a/test/symbolic_indexing_interface.jl b/test/symbolic_indexing_interface.jl index 0f42ed7d34..e46c197dde 100644 --- a/test/symbolic_indexing_interface.jl +++ b/test/symbolic_indexing_interface.jl @@ -8,6 +8,7 @@ using SciMLStructures: Tunable eqs = [D(x) ~ a * y + t, D(y) ~ b * t] @named odesys = ODESystem(eqs, t, [x, y], [a, b]; observed = [xy ~ x + y]) odesys = complete(odesys) + @test SymbolicIndexingInterface.supports_tuple_observed(odesys) @test all(is_variable.((odesys,), [x, y, 1, 2, :x, :y])) @test all(.!is_variable.((odesys,), [a, b, t, 3, 0, :a, :b])) @test variable_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) == @@ -33,6 +34,14 @@ 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]) + getter = getu(odesys, (x + 1, x + 2)) + @test getter(prob) isa Tuple + @test_nowarn @inferred getter(prob) + getter = getp(odesys, (a + 1, a + 2)) + @test getter(prob) isa Tuple + @test_nowarn @inferred getter(prob) + @named odesys = ODESystem( eqs, t, [x, y], [a, b]; defaults = [xy => 3.0], observed = [xy ~ x + y]) odesys = complete(odesys) @@ -99,6 +108,7 @@ end 0 ~ x * y - β * z] @named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β]) ns = complete(ns) + @test SymbolicIndexingInterface.supports_tuple_observed(ns) @test !is_time_dependent(ns) ps = ModelingToolkit.MTKParameters(ns, [σ => 1.0, ρ => 2.0, β => 3.0]) pobs = parameter_observed(ns, σ + ρ) @@ -107,6 +117,15 @@ end pobs = parameter_observed(ns, [σ + ρ, ρ + β]) @test isempty(get_all_timeseries_indexes(ns, [σ + ρ, ρ + β])) @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]) + getter = getu(ns, (x + 1, x + 2)) + @test getter(prob) isa Tuple + @test_nowarn @inferred getter(prob) + getter = getp(ns, (σ + 1, σ + 2)) + @test getter(prob) isa Tuple + @test_nowarn @inferred getter(prob) end @testset "PDESystem" begin @@ -194,3 +213,14 @@ end prob = ODEProblem(sys, [], (0.0, 1.0)) @test prob.ps[b] == prob.ps[:b] end + +@testset "`get_all_timeseries_indexes` with non-split systems" begin + @variables x(t) y(t) z(t) + @parameters a + @named sys = ODESystem([D(x) ~ a * x, y ~ 2x, z ~ 0.0], t) + sys = structural_simplify(sys, split = false) + for sym in [x, y, z, x + y, x + a, y / x] + @test only(get_all_timeseries_indexes(sys, sym)) == ContinuousTimeseries() + end + @test isempty(get_all_timeseries_indexes(sys, a)) +end diff --git a/test/test_variable_metadata.jl b/test/test_variable_metadata.jl index 7f9799edb3..7b093be366 100644 --- a/test/test_variable_metadata.jl +++ b/test/test_variable_metadata.jl @@ -10,6 +10,47 @@ using ModelingToolkit @test !hasbounds(y) @test !haskey(ModelingToolkit.dump_variable_metadata(y), :bounds) +@variables y[1:3] +@test !hasbounds(y) +@test getbounds(y)[1] == [-Inf, -Inf, -Inf] +@test getbounds(y)[2] == [Inf, Inf, Inf] +for i in eachindex(y) + @test !hasbounds(y[i]) + b = getbounds(y[i]) + @test b[1] == -Inf && b[2] == Inf +end + +@variables y[1:3] [bounds = (-1, 1)] +@test hasbounds(y) +@test getbounds(y)[1] == -ones(3) +@test getbounds(y)[2] == ones(3) +for i in eachindex(y) + @test hasbounds(y[i]) + b = getbounds(y[i]) + @test b[1] == -1.0 && b[2] == 1.0 +end +@test getbounds(y[1:2])[1] == -ones(2) +@test getbounds(y[1:2])[2] == ones(2) + +@variables y[1:2, 1:2] [bounds = (-1, [1.0 Inf; 2.0 3.0])] +@test hasbounds(y) +@test getbounds(y)[1] == [-1 -1; -1 -1] +@test getbounds(y)[2] == [1.0 Inf; 2.0 3.0] +for i in eachindex(y) + @test hasbounds(y[i]) + b = getbounds(y[i]) + @test b[1] == -1 && b[2] == [1.0 Inf; 2.0 3.0][i] +end + +@variables y[1:2] [bounds = (-Inf, [1.0, Inf])] +@test hasbounds(y) +@test getbounds(y)[1] == [-Inf, -Inf] +@test getbounds(y)[2] == [1.0, Inf] +@test hasbounds(y[1]) +@test getbounds(y[1]) == (-Inf, 1.0) +@test !hasbounds(y[2]) +@test getbounds(y[2]) == (-Inf, Inf) + # Guess @variables y [guess = 0] @test getguess(y) === 0 diff --git a/test/variable_utils.jl b/test/variable_utils.jl index d76f2f1209..ecc2421955 100644 --- a/test/variable_utils.jl +++ b/test/variable_utils.jl @@ -1,6 +1,7 @@ using ModelingToolkit, Test -using ModelingToolkit: value, vars +using ModelingToolkit: value, vars, parse_variable using SymbolicUtils: <ₑ + @parameters α β δ expr = (((1 / β - 1) + δ) / α)^(1 / (α - 1)) ref = sort([β, δ, α], lt = <ₑ) @@ -41,3 +42,106 @@ ts = collect_ivs([eq]) res = vars(fn([x, y], z)) @test length(res) == 3 end + +@testset "parse_variable with iv: $iv" for iv in [t, only(@independent_variables tt)] + D = Differential(iv) + function Lorenz(; name) + @variables begin + x(iv) + y(iv) + z(iv) + end + @parameters begin + σ + ρ + β + end + sys = ODESystem( + [D(D(x)) ~ σ * (y - x) + D(y) ~ x * (ρ - z) - y + D(z) ~ x * y - β * z], iv; name) + end + function ArrSys(; name) + @variables begin + x(iv)[1:2] + end + @parameters begin + p[1:2, 1:2] + end + sys = ODESystem([D(D(x)) ~ p * x], iv; name) + end + function Outer(; name) + @named 😄 = Lorenz() + @named arr = ArrSys() + sys = ODESystem(Equation[], iv; name, systems = [😄, arr]) + end + + @mtkbuild sys = Outer() + for (str, var) in [ + # unicode system, scalar variable + ("😄.x", sys.😄.x), + ("😄.x($iv)", sys.😄.x), + ("😄₊x", sys.😄.x), + ("😄₊x($iv)", sys.😄.x), + # derivative + ("D(😄.x)", D(sys.😄.x)), + ("D(😄.x($iv))", D(sys.😄.x)), + ("D(😄₊x)", D(sys.😄.x)), + ("D(😄₊x($iv))", D(sys.😄.x)), + ("Differential($iv)(😄.x)", D(sys.😄.x)), + ("Differential($iv)(😄.x($iv))", D(sys.😄.x)), + ("Differential($iv)(😄₊x)", D(sys.😄.x)), + ("Differential($iv)(😄₊x($iv))", D(sys.😄.x)), + # other derivative + ("😄.xˍ$iv", D(sys.😄.x)), + ("😄.x($iv)ˍ$iv", D(sys.😄.x)), + ("😄₊xˍ$iv", D(sys.😄.x)), + ("😄₊x($iv)ˍ$iv", D(sys.😄.x)), + # scalar parameter + ("😄.σ", sys.😄.σ), + ("😄₊σ", sys.😄.σ), + # array variable + ("arr.x", sys.arr.x), + ("arr₊x", sys.arr.x), + ("arr.x($iv)", sys.arr.x), + ("arr₊x($iv)", sys.arr.x), + # getindex + ("arr.x[1]", sys.arr.x[1]), + ("arr₊x[1]", sys.arr.x[1]), + ("arr.x($iv)[1]", sys.arr.x[1]), + ("arr₊x($iv)[1]", sys.arr.x[1]), + # derivative + ("D(arr.x($iv))", D(sys.arr.x)), + ("D(arr₊x($iv))", D(sys.arr.x)), + ("D(arr.x[1])", D(sys.arr.x[1])), + ("D(arr₊x[1])", D(sys.arr.x[1])), + ("D(arr.x($iv)[1])", D(sys.arr.x[1])), + ("D(arr₊x($iv)[1])", D(sys.arr.x[1])), + ("Differential($iv)(arr.x($iv))", D(sys.arr.x)), + ("Differential($iv)(arr₊x($iv))", D(sys.arr.x)), + ("Differential($iv)(arr.x[1])", D(sys.arr.x[1])), + ("Differential($iv)(arr₊x[1])", D(sys.arr.x[1])), + ("Differential($iv)(arr.x($iv)[1])", D(sys.arr.x[1])), + ("Differential($iv)(arr₊x($iv)[1])", D(sys.arr.x[1])), + # other derivative + ("arr.xˍ$iv", D(sys.arr.x)), + ("arr₊xˍ$iv", D(sys.arr.x)), + ("arr.xˍ$iv($iv)", D(sys.arr.x)), + ("arr₊xˍ$iv($iv)", D(sys.arr.x)), + ("arr.xˍ$iv[1]", D(sys.arr.x[1])), + ("arr₊xˍ$iv[1]", D(sys.arr.x[1])), + ("arr.xˍ$iv($iv)[1]", D(sys.arr.x[1])), + ("arr₊xˍ$iv($iv)[1]", D(sys.arr.x[1])), + ("arr.x($iv)ˍ$iv", D(sys.arr.x)), + ("arr₊x($iv)ˍ$iv", D(sys.arr.x)), + ("arr.x($iv)ˍ$iv[1]", D(sys.arr.x[1])), + ("arr₊x($iv)ˍ$iv[1]", D(sys.arr.x[1])), + # array parameter + ("arr.p", sys.arr.p), + ("arr₊p", sys.arr.p), + ("arr.p[1, 2]", sys.arr.p[1, 2]), + ("arr₊p[1, 2]", sys.arr.p[1, 2]) + ] + @test isequal(parse_variable(sys, str), var) + end +end