docs: update docs for MTKv9

NEWS.md

@@ -49,3 +49,4 @@
     `(single_parameter => expression_involving_other_parameters)` and a `Vector` of these can be passed to
     the `parameter_dependencies` keyword argument of `ODESystem`, `SDESystem` and `JumpSystem`. The dependent
     parameters are updated whenever other parameters are modified, e.g. in callbacks.
+  - Support for `IfElse.jl` has been dropped. `Base.ifelse` can be used instead.

Project.toml

@@ -23,7 +23,6 @@ FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
@@ -81,7 +80,6 @@ FindFirstFunctions = "1"
 ForwardDiff = "0.10.3"
 FunctionWrappersWrappers = "0.1"
 Graphs = "1.5.2"
-IfElse = "0.1"
 InteractiveUtils = "1"
 JuliaFormatter = "1.0.47"
 JumpProcesses = "9.1"
@@ -116,7 +114,6 @@ julia = "1.9"
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
-ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
@@ -139,4 +136,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg"]
+test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg"]

docs/Project.toml

@@ -4,17 +4,17 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-ModelingToolkitDesigner = "23d639d0-9462-4d1e-84fe-d700424865b8"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
@@ -26,15 +26,15 @@ DifferentialEquations = "7.6"
 Distributions = "0.25"
 Documenter = "1"
 ModelingToolkit = "8.33, 9"
-ModelingToolkitDesigner = "1"
+DynamicQuantities = "^0.11.2"
 NonlinearSolve = "0.3, 1, 2, 3"
 Optim = "1.7"
 Optimization = "3.9"
 OptimizationOptimJL = "0.1"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
 StochasticDiffEq = "6"
-StructuralIdentifiability = "0.4, 0.5"
+SymbolicIndexingInterface = "0.3.1"
 SymbolicUtils = "1"
 Symbolics = "5"
 Unitful = "1.12"

docs/src/basics/Events.md

@@ -144,12 +144,12 @@ ModelingToolkit therefore supports regular Julia functions as affects: instead
 of one or more equations, an affect is defined as a `tuple`:
 
 ```julia
-[x ~ 0] => (affect!, [v, x], [p, q], ctx)
+[x ~ 0] => (affect!, [v, x], [p, q], [discretes...], ctx)
 ```
 
 where, `affect!` is a Julia function with the signature: `affect!(integ, u, p, ctx)`; `[u,v]` and `[p,q]` are the symbolic unknowns (variables) and parameters
-that are accessed by `affect!`, respectively; and `ctx` is any context that is
-passed to `affect!` as the `ctx` argument.
+that are accessed by `affect!`, respectively; `discretes` are the parameters modified by `affect!`, if any;
+and `ctx` is any context that is passed to `affect!` as the `ctx` argument.
 
 `affect!` receives a [DifferentialEquations.jl
 integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)
@@ -172,7 +172,7 @@ When accessing variables of a sub-system, it can be useful to rename them
 (alternatively, an affect function may be reused in different contexts):
 
 ```julia
-[x ~ 0] => (affect!, [resistor₊v => :v, x], [p, q => :p2], ctx)
+[x ~ 0] => (affect!, [resistor₊v => :v, x], [p, q => :p2], [], ctx)
 ```
 
 Here, the symbolic variable `resistor₊v` is passed as `v` while the symbolic
@@ -191,7 +191,7 @@ function bb_affect!(integ, u, p, ctx)
     integ.u[u.v] = -integ.u[u.v]
 end
 
-reflect = [x ~ 0] => (bb_affect!, [v], [], nothing)
+reflect = [x ~ 0] => (bb_affect!, [v], [], [], nothing)
 
 @mtkbuild bb_sys = ODESystem(bb_eqs, t, sts, par,
     continuous_events = reflect)

docs/src/basics/FAQ.md

@@ -24,12 +24,11 @@ pnew = varmap_to_vars([β => 3.0, c => 10.0, γ => 2.0], parameters(sys))
 
 ## How do I handle `if` statements in my symbolic forms?
 
-For statements that are in the `if then else` form, use `IfElse.ifelse` from the
-[IfElse.jl](https://github.com/SciML/IfElse.jl) package to represent the code in a
-functional form. For handling direct `if` statements, you can use equivalent boolean
-mathematical expressions. For example, `if x > 0 ...` can be implemented as just
-`(x > 0) * `, where if `x <= 0` then the boolean will evaluate to `0` and thus the
-term will be excluded from the model.
+For statements that are in the `if then else` form, use `Base.ifelse` from the
+to represent the code in a functional form. For handling direct `if` statements,
+you can use equivalent boolean mathematical expressions. For example, `if x > 0 ...`
+can be implemented as just `(x > 0) * `, where if `x <= 0` then the boolean will
+evaluate to `0` and thus the term will be excluded from the model.
 
 ## ERROR: TypeError: non-boolean (Num) used in boolean context?
 

docs/src/basics/Validation.md

@@ -7,7 +7,7 @@ ModelingToolkit.jl provides extensive functionality for model validation and uni
 Units may be assigned with the following syntax.
 
 ```@example validation
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @variables t [unit = u"s"] x(t) [unit = u"m"] g(t) w(t) [unit = "Hz"]
 
 @variables(t, [unit = u"s"], x(t), [unit = u"m"], g(t), w(t), [unit = "Hz"])
@@ -45,7 +45,7 @@ ModelingToolkit.get_unit
 Example usage below. Note that `ModelingToolkit` does not force unit conversions to preferred units in the event of nonstandard combinations -- it merely checks that the equations are consistent.
 
 ```@example validation
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @parameters τ [unit = u"ms"]
 @variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
@@ -59,13 +59,17 @@ ModelingToolkit.validate(eqs[1])
 ```
 
 ```@example validation
-ModelingToolkit.get_unit(eqs[1].rhs)
+try
+    ModelingToolkit.get_unit(eqs[1].rhs)
+catch e
+    showerror(stdout, e)
+end
 ```
 
 An example of an inconsistent system: at present, `ModelingToolkit` requires that the units of all terms in an equation or sum to be equal-valued (`ModelingToolkit.equivalent(u1,u2)`), rather than simply dimensionally consistent. In the future, the validation stage may be upgraded to support the insertion of conversion factors into the equations.
 
 ```@example validation
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @parameters τ [unit = u"ms"]
 @variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
 D = Differential(t)
@@ -81,10 +85,9 @@ expect an object type, while two-parameter calls expect a function type as the f
 second argument.
 
 ```@example validation2
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
+using ModelingToolkit: t_nounits as t, D_nounits as D
 # Composite type parameter in registered function
-@variables t
-D = Differential(t)
 struct NewType
     f::Any
 end
@@ -101,16 +104,17 @@ end
 sts = @variables a(t)=0 [unit = u"cm"]
 ps = @parameters s=-1 [unit = u"cm"] c=c [unit = u"cm"]
 eqs = [D(a) ~ dummycomplex(c, s);]
-sys = ODESystem(eqs, t, [sts...;], [ps...;], name = :sys)
+sys = ODESystem(
+    eqs, t, [sts...;], [ps...;], name = :sys, checks = ~ModelingToolkit.CheckUnits)
 sys_simple = structural_simplify(sys)
 ```
 
-## `Unitful` Literals
+## `DynamicQuantities` Literals
 
 In order for a function to work correctly during both validation & execution, the function must be unit-agnostic. That is, no unitful literals may be used. Any unitful quantity must either be a `parameter` or `variable`. For example, these equations will not validate successfully.
 
 ```julia
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
 D = Differential(t)
 eqs = [D(E) ~ P - E / 1u"ms"]
@@ -124,7 +128,7 @@ ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
 Instead, they should be parameterized:
 
 ```@example validation3
-using ModelingToolkit, Unitful
+using ModelingToolkit, DynamicQuantities
 @parameters τ [unit = u"ms"]
 @variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
@@ -138,8 +142,8 @@ eqs = [D(E) ~ P - myfunc(E, τ)]
 ModelingToolkit.validate(eqs) #Returns true
 ```
 
-It is recommended *not* to circumvent unit validation by specializing user-defined functions on `Unitful` arguments vs. `Numbers`. This both fails to take advantage of `validate` for ensuring correctness, and may cause in errors in the
-future when `ModelingToolkit` is extended to support eliminating `Unitful` literals from functions.
+It is recommended *not* to circumvent unit validation by specializing user-defined functions on `DynamicQuantities` arguments vs. `Numbers`. This both fails to take advantage of `validate` for ensuring correctness, and may cause in errors in the
+future when `ModelingToolkit` is extended to support eliminating `DynamicQuantities` literals from functions.
 
 ## Other Restrictions
 
@@ -149,7 +153,7 @@ future when `ModelingToolkit` is extended to support eliminating `Unitful` liter
 
 If a system fails to validate due to unit issues, at least one warning message will appear, including a line number as well as the unit types and expressions that were in conflict. Some system constructors re-order equations before the unit checking can be done, in which case the equation numbers may be inaccurate. The printed expression that the problem resides in is always correctly shown.
 
-Symbolic exponents for unitful variables *are* supported (ex: `P^γ` in thermodynamics). However, this means that `ModelingToolkit` cannot reduce such expressions to `Unitful.Unitlike` subtypes at validation time because the exponent value is not available. In this case `ModelingToolkit.get_unit` is type-unstable, yielding a symbolic result, which can still be checked for symbolic equality with `ModelingToolkit.equivalent`.
+Symbolic exponents for unitful variables *are* supported (ex: `P^γ` in thermodynamics). However, this means that `ModelingToolkit` cannot reduce such expressions to `DynamicQuantities.Quantity` subtypes at validation time because the exponent value is not available. In this case `ModelingToolkit.get_unit` is type-unstable, yielding a symbolic result, which can still be checked for symbolic equality with `ModelingToolkit.equivalent`.
 
 ## Parameter & Initial Condition Values
 

docs/src/basics/Variable_metadata.md

@@ -142,7 +142,7 @@ In the example below, we define a system with tunable parameters and extract bou
 @variables x(t)=0 u(t)=0 [input = true] y(t)=0 [output = true]
 @parameters T [tunable = true, bounds = (0, Inf)]
 @parameters k [tunable = true, bounds = (0, Inf)]
-eqs = [Dₜ(x) ~ (-x + k * u) / T # A first-order system with time constant T and gain k
+eqs = [D(x) ~ (-x + k * u) / T # A first-order system with time constant T and gain k
        y ~ x]
 sys = ODESystem(eqs, t, name = :tunable_first_order)
 ```
@@ -159,6 +159,9 @@ lb, ub = getbounds(p) # operating on a vector, we get lower and upper bound vect
 b = getbounds(sys) # Operating on the system, we get a dict
 ```
 
+See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref),
+[`ModelingToolkit.dump_unknowns`](@ref).
+
 ## Index
 
 ```@index
@@ -172,3 +175,9 @@ Modules = [ModelingToolkit]
 Pages = ["variables.jl"]
 Private = false
 ```
+
+```@docs
+ModelingToolkit.dump_variable_metadata
+ModelingToolkit.dump_parameters
+ModelingToolkit.dump_unknowns
+```

docs/src/examples/parsing.md

@@ -23,10 +23,11 @@ From there, we can use ModelingToolkit to transform the symbolic equations into
 nonlinear solve:
 
 ```@example parsing
-using ModelingToolkit, NonlinearSolve
+using ModelingToolkit, SymbolicIndexingInterface, NonlinearSolve
 vars = union(ModelingToolkit.vars.(eqs)...)
 @mtkbuild ns = NonlinearSystem(eqs, vars, [])
 
-prob = NonlinearProblem(ns, [1.0, 1.0, 1.0])
+varmap = Dict(SymbolicIndexingInterface.getname.(vars) .=> vars)
+prob = NonlinearProblem(ns, [varmap[:x] => 1.0, varmap[:y] => 1.0, varmap[:z] => 1.0])
 sol = solve(prob, NewtonRaphson())
 ```

docs/src/internals.md

@@ -32,6 +32,6 @@ The call chain typically looks like this, with the function names in the case of
 
  1. Problem constructor ([`ODEProblem`](@ref))
  2. Build an `DEFunction` ([`process_DEProblem`](@ref) -> [`ODEFunction`](@ref)
- 3. Write actual executable code ([`generate_function`](@ref))
+ 3. Write actual executable code ([`generate_function`](@ref) or [`generate_custom_function`](@ref))
 
 Apart from [`generate_function`](@ref), which generates the dynamics function, `ODEFunction` also builds functions for observed equations (`build_explicit_observed_function`) and Jacobians (`generate_jacobian`) etc. These are all stored in the `ODEFunction`.

docs/src/tutorials/domain_connections.md

@@ -94,7 +94,7 @@ end
 nothing #hide
 ```
 
-When the system is defined we can generate a fluid component and connect it to the system.  Here `fluid` is connected to the `src.port`, but it could also be connected to `vol.port`, any connection in the network is fine.  Note: we can visualize the system using `ModelingToolkitDesigner.jl`, where a dashed line is used to show the `fluid` connection to represent a domain connection that is only transporting parameters and not unknown variables.
+When the system is defined we can generate a fluid component and connect it to the system.  Here `fluid` is connected to the `src.port`, but it could also be connected to `vol.port`, any connection in the network is fine.
 
 ```@example domain
 @component function System(; name)
@@ -115,20 +115,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# code to generate diagrams...
-# using ModelingToolkitDesigner 
-# path = raw"C:\Work\Assets\ModelingToolkit.jl\domain_connections" 
-# design = ODESystemDesign(odesys, path); 
-
-# using CairoMakie
-# CairoMakie.set_theme!(Theme(;fontsize=12))
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "odesys.svg"), fig; resolution=(300,300))
-```
-
-![odesys](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/d19fbcf4-781c-4743-87b7-30bed348ff98)
-
 To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`.  First we assemble the system using `structural_simplify()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid₊ρ`.  Likewise, `src.port.ρ`, will also point to the setter value `fluid₊ρ`.  Therefore, there is now only 1 defined density value `fluid₊ρ` which sets the density for the connected network.
 
 ```@repl domain
@@ -195,14 +181,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# design = ODESystemDesign(actsys2, path);
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "actsys2.svg"), fig; resolution=(500,300))
-```
-
-![actsys2](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/8ed50035-f6ac-48cb-a585-1ef415154a02)
-
 After running `structural_simplify()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a₊ρ` and `act.port_b.ρ` points to `fluid_b₊ρ`.  This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks.  Therefore, we can connect a single fluid to both networks.  This does not interfere with the mathematical equations of the system, since no unknown variables are connected.
 
 ```@example domain
@@ -227,14 +205,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# design = ODESystemDesign(actsys1, path);
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "actsys1.svg"), fig; resolution=(500,300))
-```
-
-![actsys1](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/054404eb-dbb7-4b85-95c0-c9503d0c4d00)
-
 ## Special Connection Cases (`domain_connect()`)
 
 In some cases a component will be defined with 2 connectors of the same domain, but they are not connected.  For example the `Restrictor` defined here gives equations to define the behavior of how the 2 connectors `port_a` and `port_b` are physically connected.
@@ -282,14 +252,6 @@ end
 nothing #hide
 ```
 
-```@setup domain
-# design = ODESystemDesign(ressys, path);
-# fig = ModelingToolkitDesigner.view(design, false)
-# save(joinpath(path, "ressys.svg"), fig; resolution=(500,300))
-```
-
-![ressys](https://github.com/SciML/ModelingToolkit.jl/assets/40798837/3740f0e2-7324-4c1f-af8b-eba02cfece81)
-
 When `structural_simplify()` is applied to this system it can be seen that the defaults are missing for `res.port_b` and `vol.port`.
 
 ```@repl domain

docs/src/tutorials/ode_modeling.md

@@ -338,7 +338,7 @@ While defining the model `UnitstepFOLFactory`, an initial guess of 0.0 is assign
 Additionally, these initial guesses can be modified while creating instances of `UnitstepFOLFactory` by passing arguments.
 
 ```@example ode2
-@named fol = UnitstepFOLFactory(; x = 0.1)
+@mtkbuild fol = UnitstepFOLFactory(; x = 0.1)
 sol = ODEProblem(fol, [], (0.0, 5.0), []) |> solve
 ```
 

docs/src/tutorials/optimization.md

@@ -6,13 +6,13 @@ Let's solve the classical _Rosenbrock function_ in two dimensions.
 
 First, we need to make some imports.
 
-```@example rosenbrock_2d
+```julia
 using ModelingToolkit, Optimization, OptimizationOptimJL
 ```
 
 Now we can define our optimization problem.
 
-```@example rosenbrock_2d
+```julia
 @variables begin
     x, [bounds = (-2.0, 2.0)]
     y, [bounds = (-1.0, 3.0)]
@@ -44,7 +44,7 @@ Every optimization problem consists of a set of _optimization variables_. In thi
 
 Next, the actual `OptimizationProblem` can be created. At this stage, an initial guess `u0` for the optimization variables needs to be provided via map, using the symbols from before. Concrete values for the parameters of the system can also be provided or changed. However, if the parameters have default values assigned, they are used automatically.
 
-```@example rosenbrock_2d
+```julia
 u0 = [x => 1.0
       y => 2.0]
 p = [a => 1.0
@@ -56,7 +56,7 @@ solve(prob, GradientDescent())
 
 ## Rosenbrock Function with Constraints
 
-```@example rosenbrock_2d_cstr
+```julia
 using ModelingToolkit, Optimization, OptimizationOptimJL
 
 @variables begin