From eb976a85b0fe65f59382132a4c6bd176b434b7f7 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 14:16:49 +0530
Subject: [PATCH 01/18] docs: fix FunctionAffect usage and documentation

---
 docs/src/basics/Events.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md
index 91df993f70..61d59863e2 100644
--- a/docs/src/basics/Events.md
+++ b/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)

From 0f6931a425f06f133a8473c5b91ab4fe560b01f7 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 14:18:58 +0530
Subject: [PATCH 02/18] docs: remove usages of ModelingToolkitDesigner

---
 docs/Project.toml                        |  2 --
 docs/src/tutorials/domain_connections.md | 40 +-----------------------
 2 files changed, 1 insertion(+), 41 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 49fba20527..4c11aca431 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,7 +6,6 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 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"
@@ -26,7 +25,6 @@ DifferentialEquations = "7.6"
 Distributions = "0.25"
 Documenter = "1"
 ModelingToolkit = "8.33, 9"
-ModelingToolkitDesigner = "1"
 NonlinearSolve = "0.3, 1, 2, 3"
 Optim = "1.7"
 Optimization = "3.9"
diff --git a/docs/src/tutorials/domain_connections.md b/docs/src/tutorials/domain_connections.md
index b381668544..d6dc2d8781 100644
--- a/docs/src/tutorials/domain_connections.md
+++ b/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

From 82a570f8c77cadb7ee18b9959e41b6f3c3f9f1e5 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 14:20:52 +0530
Subject: [PATCH 03/18] docs: remove StructuralIndentifiability dependency

Doc examples also commented out, to be added when
StructuralIndentifiability works with v9
---
 docs/Project.toml                               |  2 --
 docs/src/tutorials/parameter_identifiability.md | 10 +++++-----
 2 files changed, 5 insertions(+), 7 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 4c11aca431..cfb4dc6e22 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -13,7 +13,6 @@ 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"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
@@ -32,7 +31,6 @@ OptimizationOptimJL = "0.1"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
 StochasticDiffEq = "6"
-StructuralIdentifiability = "0.4, 0.5"
 SymbolicUtils = "1"
 Symbolics = "5"
 Unitful = "1.12"
diff --git a/docs/src/tutorials/parameter_identifiability.md b/docs/src/tutorials/parameter_identifiability.md
index bcd027dcf0..54c3d0f42d 100644
--- a/docs/src/tutorials/parameter_identifiability.md
+++ b/docs/src/tutorials/parameter_identifiability.md
@@ -28,7 +28,7 @@ To define the ode system in Julia, we use `ModelingToolkit.jl`.
 
 We first define the parameters, variables, differential equations and the output equations.
 
-```@example SI
+```julia
 using StructuralIdentifiability, ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
@@ -66,7 +66,7 @@ end
 
 After that, we are ready to check the system for local identifiability:
 
-```@example SI
+```julia
 # query local identifiability
 # we pass the ode-system
 local_id_all = assess_local_identifiability(de, p = 0.99)
@@ -76,7 +76,7 @@ We can see that all unknowns (except $x_7$) and all parameters are locally ident
 
 Let's try to check specific parameters and their combinations
 
-```@example SI
+```julia
 to_check = [de.k5, de.k7, de.k10 / de.k9, de.k5 + de.k6]
 local_id_some = assess_local_identifiability(de, funcs_to_check = to_check, p = 0.99)
 ```
@@ -103,7 +103,7 @@ We will run a global identifiability check on this enzyme dynamics[^3] model. We
 
 Global identifiability needs information about local identifiability first, but the function we chose here will take care of that extra step for us.
 
-```@example SI2
+```julia
 using StructuralIdentifiability, ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
@@ -144,7 +144,7 @@ We can see that only parameters `a, g` are unidentifiable, and everything else c
 
 Let us consider the same system but with two inputs, and we will find out identifiability with probability `0.9` for parameters `c` and `b`:
 
-```@example SI3
+```julia
 using StructuralIdentifiability, ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 

From 6e010f948e95a5535531048b17b761d1b3beeed4 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 14:17:02 +0530
Subject: [PATCH 04/18] docs: use DynamicQuantities instead of Unitful

---
 docs/Project.toml             |  2 ++
 docs/src/basics/Validation.md | 23 +++++++++++------------
 2 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index cfb4dc6e22..4c90d35b77 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,6 +4,7 @@ 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"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
@@ -24,6 +25,7 @@ DifferentialEquations = "7.6"
 Distributions = "0.25"
 Documenter = "1"
 ModelingToolkit = "8.33, 9"
+DynamicQuantities = "^0.11.2"
 NonlinearSolve = "0.3, 1, 2, 3"
 Optim = "1.7"
 Optimization = "3.9"
diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md
index 346110f6d5..ada0e57810 100644
--- a/docs/src/basics/Validation.md
+++ b/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)
@@ -65,7 +65,7 @@ ModelingToolkit.get_unit(eqs[1].rhs)
 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 +81,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
@@ -105,12 +104,12 @@ sys = ODESystem(eqs, t, [sts...;], [ps...;], name = :sys)
 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 +123,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 +137,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 +148,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
 

From 9f9d8bf27fec018e00577cdc504dc431fd00b1df Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 15:41:43 +0530
Subject: [PATCH 05/18] docs: use `@brownian` in the stochastic_diffeq tutorial

---
 docs/src/tutorials/stochastic_diffeq.md | 18 +++++++-----------
 1 file changed, 7 insertions(+), 11 deletions(-)

diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md
index 91ec05ef01..85c2061af0 100644
--- a/docs/src/tutorials/stochastic_diffeq.md
+++ b/docs/src/tutorials/stochastic_diffeq.md
@@ -5,7 +5,7 @@ to the model: randomness. This is a
 [stochastic differential equation](https://en.wikipedia.org/wiki/Stochastic_differential_equation)
 which has a deterministic (drift) component and a stochastic (diffusion)
 component. Let's take the Lorenz equation from the first tutorial and extend
-it to have multiplicative noise.
+it to have multiplicative noise by creating `@brownian` variables in the equations.
 
 ```@example SDE
 using ModelingToolkit, StochasticDiffEq
@@ -14,16 +14,12 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 # Define some variables
 @parameters σ ρ β
 @variables x(t) y(t) z(t)
+@brownian a
+eqs = [D(x) ~ σ * (y - x) + 0.1a * x,
+    D(y) ~ x * (ρ - z) - y + 0.1a * y,
+    D(z) ~ x * y - β * z + 0.1a * z]
 
-eqs = [D(x) ~ σ * (y - x),
-    D(y) ~ x * (ρ - z) - y,
-    D(z) ~ x * y - β * z]
-
-noiseeqs = [0.1 * x,
-    0.1 * y,
-    0.1 * z]
-
-@mtkbuild de = SDESystem(eqs, noiseeqs, t, [x, y, z], [σ, ρ, β])
+@mtkbuild de = System(eqs, t)
 
 u0map = [
     x => 1.0,
@@ -38,5 +34,5 @@ parammap = [
 ]
 
 prob = SDEProblem(de, u0map, (0.0, 100.0), parammap)
-sol = solve(prob, SOSRI())
+sol = solve(prob, LambdaEulerHeun())
 ```

From 940d637f0e7765af54e5b8c607cd9f0812a3c7fe Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 15:55:05 +0530
Subject: [PATCH 06/18] refactor: add warning that `substitute` does not update
 events

---
 src/systems/abstractsystem.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 7bb3c84818..c9d0a8fb0b 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2113,6 +2113,10 @@ end
 
 keytype(::Type{<:Pair{T, V}}) where {T, V} = T
 function Symbolics.substitute(sys::AbstractSystem, rules::Union{Vector{<:Pair}, Dict})
+    if has_continuous_domain(sys) && get_continuous_events(sys) !== nothing ||
+       has_discrete_events(sys) && get_discrete_events(sys) !== nothing
+        @warn "`substitute` only supports performing substitutions in equations. This system has events, which will not be updated."
+    end
     if keytype(eltype(rules)) <: Symbol
         dict = todict(rules)
         systems = get_systems(sys)

From 0272392f2d606a6e65df87c73815064e7cc4cfd3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 16:04:22 +0530
Subject: [PATCH 07/18] refactor: warn that initial values for observed
 variables will not be respected

---
 src/systems/diffeqs/abstractodesystem.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 34899e3361..a734845118 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -775,6 +775,16 @@ function get_u0_p(sys,
     if parammap !== nothing
         defs = mergedefaults(defs, parammap, ps)
     end
+    if u0map isa Vector && eltype(u0map) <: Pair
+        u0map = Dict(u0map)
+    end
+    if u0map isa Dict
+        allobs = Set(getproperty.(observed(sys), :lhs))
+        if any(in(allobs), keys(u0map))
+            u0s_in_obs = filter(in(allobs), keys(u0map))
+            @warn "Observed variables cannot assigned initial values. Initial values for $u0s_in_obs will be ignored."
+        end
+    end
     defs = mergedefaults(defs, u0map, dvs)
     for (k, v) in defs
         if Symbolics.isarraysymbolic(k)

From 5cbd523383128744778daa2fbc786a99b009f8a4 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 18:02:11 +0530
Subject: [PATCH 08/18] feat: add dump_variable_metadata, dump_parameters,
 dump_variables

Also make `istunable` default to `true` for parameters without tunable
metadata. `tunable_parameters` is also updated accordingly
---
 src/systems/abstractsystem.jl  | 20 ++++++++++++
 src/variables.jl               | 59 +++++++++++++++++++++++++++++++---
 test/test_variable_metadata.jl | 36 +++++++++++++++++----
 3 files changed, 104 insertions(+), 11 deletions(-)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index c9d0a8fb0b..0d1684bfb6 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2146,3 +2146,23 @@ function process_parameter_dependencies(pdeps, ps)
     end
     return pdeps, ps
 end
+
+function dump_parameters(sys::AbstractSystem)
+    defs = defaults(sys)
+    map(dump_variable_metadata.(parameters(sys))) do meta
+        if haskey(defs, meta.var)
+            meta = merge(meta, (; default = defs[meta.var]))
+        end
+        meta
+    end
+end
+
+function dump_unknowns(sys::AbstractSystem)
+    defs = defaults(sys)
+    map(dump_variable_metadata.(unknowns(sys))) do meta
+        if haskey(defs, meta.var)
+            meta = merge(meta, (; default = defs[meta.var]))
+        end
+        meta
+    end
+end
diff --git a/src/variables.jl b/src/variables.jl
index ca6d1b6954..5b2205cbea 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -15,6 +15,57 @@ Symbolics.option_to_metadata_type(::Val{:irreducible}) = VariableIrreducible
 Symbolics.option_to_metadata_type(::Val{:state_priority}) = VariableStatePriority
 Symbolics.option_to_metadata_type(::Val{:misc}) = VariableMisc
 
+function dump_variable_metadata(var)
+    uvar = unwrap(var)
+    vartype, name = get(uvar.metadata, VariableSource, (:unknown, :unknown))
+    shape = Symbolics.shape(var)
+    if shape == ()
+        shape = nothing
+    end
+    unit = get(uvar.metadata, VariableUnit, nothing)
+    connect = get(uvar.metadata, VariableConnectType, nothing)
+    noise = get(uvar.metadata, VariableNoiseType, nothing)
+    input = isinput(uvar) || nothing
+    output = isoutput(uvar) || nothing
+    irreducible = get(uvar.metadata, VariableIrreducible, nothing)
+    state_priority = get(uvar.metadata, VariableStatePriority, nothing)
+    misc = get(uvar.metadata, VariableMisc, nothing)
+    bounds = hasbounds(uvar) ? getbounds(uvar) : nothing
+    desc = getdescription(var)
+    if desc == ""
+        desc = nothing
+    end
+    guess = getguess(uvar)
+    disturbance = isdisturbance(uvar) || nothing
+    tunable = istunable(uvar, isparameter(uvar))
+    dist = getdist(uvar)
+    type = symtype(uvar)
+
+    meta = (
+        var = var,
+        vartype,
+        name,
+        shape,
+        unit,
+        connect,
+        noise,
+        input,
+        output,
+        irreducible,
+        state_priority,
+        misc,
+        bounds,
+        desc,
+        guess,
+        disturbance,
+        tunable,
+        dist,
+        type
+    )
+
+    return NamedTuple(k => v for (k, v) in pairs(meta) if v !== nothing)
+end
+
 abstract type AbstractConnectType end
 struct Equality <: AbstractConnectType end # Equality connection
 struct Flow <: AbstractConnectType end     # sum to 0
@@ -243,7 +294,7 @@ Symbolics.option_to_metadata_type(::Val{:tunable}) = VariableTunable
 istunable(x::Num, args...) = istunable(Symbolics.unwrap(x), args...)
 
 """
-    istunable(x, default = false)
+    istunable(x, default = true)
 
 Determine whether symbolic variable `x` is marked as a tunable for an automatic tuning algorithm.
 
@@ -257,7 +308,7 @@ Create a tunable parameter by
 
 See also [`tunable_parameters`](@ref), [`getbounds`](@ref)
 """
-function istunable(x, default = false)
+function istunable(x, default = true)
     p = Symbolics.getparent(x, nothing)
     p === nothing || (x = p)
     Symbolics.getmetadata(x, VariableTunable, default)
@@ -302,7 +353,7 @@ end
 ## System interface
 
 """
-    tunable_parameters(sys, p = parameters(sys); default=false)
+    tunable_parameters(sys, p = parameters(sys); default=true)
 
 Get all parameters of `sys` that are marked as `tunable`.
 
@@ -316,7 +367,7 @@ Create a tunable parameter by
 
 See also [`getbounds`](@ref), [`istunable`](@ref)
 """
-function tunable_parameters(sys, p = parameters(sys); default = false)
+function tunable_parameters(sys, p = parameters(sys); default = true)
     filter(x -> istunable(x, default), p)
 end
 
diff --git a/test/test_variable_metadata.jl b/test/test_variable_metadata.jl
index 886458302e..b79df5e72d 100644
--- a/test/test_variable_metadata.jl
+++ b/test/test_variable_metadata.jl
@@ -4,31 +4,43 @@ using ModelingToolkit
 @variables u [bounds = (-1, 1)]
 @test getbounds(u) == (-1, 1)
 @test hasbounds(u)
+@test ModelingToolkit.dump_variable_metadata(u).bounds == (-1, 1)
 
 @variables y
 @test !hasbounds(y)
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :bounds)
 
 # Guess
 @variables y [guess = 0]
 @test getguess(y) === 0
 @test hasguess(y) === true
+@test ModelingToolkit.dump_variable_metadata(y).guess == 0
 
 @variables y
 @test hasguess(y) === false
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :guess)
 
 # Disturbance
 @variables u [disturbance = true]
 @test isdisturbance(u)
+@test ModelingToolkit.dump_variable_metadata(u).disturbance
 
 @variables y
 @test !isdisturbance(y)
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :disturbance)
 
 # Tunable
 @parameters u [tunable = true]
 @test istunable(u)
+@test ModelingToolkit.dump_variable_metadata(u).tunable
+
+@parameters u2 [tunable = false]
+@test !istunable(u2)
+@test !ModelingToolkit.dump_variable_metadata(u2).tunable
 
 @parameters y
-@test !istunable(y)
+@test istunable(y)
+@test ModelingToolkit.dump_variable_metadata(y).tunable
 
 # Distributions
 struct FakeNormal end
@@ -36,20 +48,28 @@ d = FakeNormal()
 @parameters u [dist = d]
 @test hasdist(u)
 @test getdist(u) == d
+@test ModelingToolkit.dump_variable_metadata(u).dist == d
 
 @parameters y
 @test !hasdist(y)
+@test !haskey(ModelingToolkit.dump_variable_metadata(y), :dist)
 
 ## System interface
 @parameters t
 Dₜ = Differential(t)
 @variables x(t)=0 [bounds = (-10, 10)] u(t)=0 [input = true] y(t)=0 [output = true]
-@parameters T [tunable = true, bounds = (0, Inf)]
+@parameters T [bounds = (0, Inf)]
 @parameters k [tunable = true, bounds = (0, Inf)]
-@parameters k2
+@parameters k2 [tunable = false]
 eqs = [Dₜ(x) ~ (-k2 * x + k * u) / T
        y ~ x]
 sys = ODESystem(eqs, t, name = :tunable_first_order)
+unk_meta = ModelingToolkit.dump_unknowns(sys)
+@test length(unk_meta) == 3
+@test all(iszero, meta.default for meta in unk_meta)
+param_meta = ModelingToolkit.dump_parameters(sys)
+@test length(param_meta) == 3
+@test all(!haskey(meta, :default) for meta in param_meta)
 
 p = tunable_parameters(sys)
 sp = Set(p)
@@ -68,21 +88,23 @@ b = getbounds(sys)
 b = getbounds(sys, unknowns(sys))
 @test b[x] == (-10, 10)
 
-p = tunable_parameters(sys, default = true)
+p = tunable_parameters(sys, default = false)
 sp = Set(p)
 @test k ∈ sp
-@test T ∈ sp
-@test k2 ∈ sp
-@test length(p) == 3
+@test T ∉ sp
+@test k2 ∉ sp
+@test length(p) == 1
 
 ## Descriptions
 @variables u [description = "This is my input"]
 @test getdescription(u) == "This is my input"
 @test hasdescription(u)
+@test ModelingToolkit.dump_variable_metadata(u).desc == "This is my input"
 
 @variables u
 @test getdescription(u) == ""
 @test !hasdescription(u)
+@test !haskey(ModelingToolkit.dump_variable_metadata(u), :desc)
 
 @parameters t
 @variables u(t) [description = "A short description of u"]

From 47aed59e37157f965a39c5f64e51c76ee099f827 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 18:48:18 +0530
Subject: [PATCH 09/18] fix: fix repack when new value does not alias
 canonicalized array

---
 src/systems/parameter_buffer.jl | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl
index 897225c47f..62e3c39f72 100644
--- a/src/systems/parameter_buffer.jl
+++ b/src/systems/parameter_buffer.jl
@@ -128,7 +128,7 @@ end
 function split_into_buffers(raw::AbstractArray, buf; recurse = true)
     idx = 1
     function _helper(buf_v; recurse = true)
-        if eltype(buf_v) isa AbstractArray && recurse
+        if eltype(buf_v) <: AbstractArray && recurse
             return _helper.(buf_v; recurse = false)
         else
             res = reshape(raw[idx:(idx + length(buf_v) - 1)], size(buf_v))
@@ -139,6 +139,19 @@ function split_into_buffers(raw::AbstractArray, buf; recurse = true)
     return Tuple(_helper(buf_v; recurse) for buf_v in buf)
 end
 
+function update_tuple_of_buffers(raw::AbstractArray, buf)
+    idx = 1
+    function _helper(buf_v)
+        if eltype(buf_v) <: AbstractArray
+            _helper.(buf_v)
+        else
+            copyto!(buf_v, view(raw, idx:(idx + length(buf_v) - 1)))
+            idx += length(buf_v)
+        end
+    end
+    _helper.(buf)
+end
+
 SciMLStructures.isscimlstructure(::MTKParameters) = true
 
 SciMLStructures.ismutablescimlstructure(::MTKParameters) = true
@@ -151,7 +164,7 @@ for (Portion, field) in [(SciMLStructures.Tunable, :tunable)
         repack = let as_vector = as_vector, p = p
             function (new_val)
                 if new_val !== as_vector
-                    p.$field = split_into_buffers(new_val, p.$field)
+                    update_tuple_of_buffers(new_val, p.$field)
                 end
                 if p.dependent_update_iip !== nothing
                     p.dependent_update_iip(ArrayPartition(p.dependent), p...)

From 16d63110e49deaef7670b9d4120cd66757d35071 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 19:07:22 +0530
Subject: [PATCH 10/18] docs: document new `dump_*` functions

---
 docs/src/basics/Variable_metadata.md |  9 +++++++
 src/systems/abstractsystem.jl        | 38 ++++++++++++++++++++++++++++
 src/variables.jl                     | 12 +++++++++
 3 files changed, 59 insertions(+)

diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/basics/Variable_metadata.md
index a973832538..c2c1f6d4ce 100644
--- a/docs/src/basics/Variable_metadata.md
+++ b/docs/src/basics/Variable_metadata.md
@@ -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
+```
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 0d1684bfb6..afc23b7cd3 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2147,6 +2147,25 @@ function process_parameter_dependencies(pdeps, ps)
     return pdeps, ps
 end
 
+"""
+    dump_parameters(sys::AbstractSystem)
+
+Return an array of `NamedTuple`s containing the metadata associated with each parameter in
+`sys`. Also includes the default value of the parameter, if provided.
+
+```@example
+using ModelingToolkit
+using DynamicQuantities
+using ModelingToolkit: t, D
+
+@parameters p = 1.0, [description = "My parameter", tunable = false] q = 2.0, [description = "Other parameter"]
+@variables x(t) = 3.0 [unit = u"m"]
+@named sys = ODESystem(Equation[], t, [x], [p, q])
+ModelingToolkit.dump_parameters(sys)
+```
+
+See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_unknowns`](@ref)
+"""
 function dump_parameters(sys::AbstractSystem)
     defs = defaults(sys)
     map(dump_variable_metadata.(parameters(sys))) do meta
@@ -2157,6 +2176,25 @@ function dump_parameters(sys::AbstractSystem)
     end
 end
 
+"""
+    dump_unknowns(sys::AbstractSystem)
+
+Return an array of `NamedTuple`s containing the metadata associated with each unknown in
+`sys`. Also includes the default value of the unknown, if provided.
+
+```@example
+using ModelingToolkit
+using DynamicQuantities
+using ModelingToolkit: t, D
+
+@parameters p = 1.0, [description = "My parameter", tunable = false] q = 2.0, [description = "Other parameter"]
+@variables x(t) = 3.0 [unit = u"m"]
+@named sys = ODESystem(Equation[], t, [x], [p, q])
+ModelingToolkit.dump_unknowns(sys)
+```
+
+See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref)
+"""
 function dump_unknowns(sys::AbstractSystem)
     defs = defaults(sys)
     map(dump_variable_metadata.(unknowns(sys))) do meta
diff --git a/src/variables.jl b/src/variables.jl
index 5b2205cbea..52db3339ef 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -15,6 +15,18 @@ Symbolics.option_to_metadata_type(::Val{:irreducible}) = VariableIrreducible
 Symbolics.option_to_metadata_type(::Val{:state_priority}) = VariableStatePriority
 Symbolics.option_to_metadata_type(::Val{:misc}) = VariableMisc
 
+"""
+    dump_variable_metadata(var)
+
+Return all the metadata associated with symbolic variable `var` as a `NamedTuple`.
+
+```@example
+using ModelingToolkit
+
+@parameters p::Int [description = "My description", bounds = (0.5, 1.5)]
+ModelingToolkit.dump_variable_metadata(p)
+```
+"""
 function dump_variable_metadata(var)
     uvar = unwrap(var)
     vartype, name = get(uvar.metadata, VariableSource, (:unknown, :unknown))

From 48ec10b53139bd6761037eec8848012d034b674e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 20:18:40 +0530
Subject: [PATCH 11/18] refactor: remove `integer` and `binary` variable
 metadata

---
 src/ModelingToolkit.jl                        |  3 +-
 src/systems/model_parsing.jl                  |  4 +--
 .../optimization/optimizationsystem.jl        | 14 ++++----
 src/variables.jl                              | 34 -------------------
 test/model_parsing.jl                         |  8 ++---
 test/test_variable_metadata.jl                | 16 ---------
 6 files changed, 14 insertions(+), 65 deletions(-)

diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 224ed22724..7e0964b541 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -222,8 +222,7 @@ export connect, domain_connect, @connector, Connection, Flow, Stream, instream
 export @component, @mtkmodel, @mtkbuild
 export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbance,
        istunable, getdist, hasdist,
-       tunable_parameters, isirreducible, getdescription, hasdescription, isbinaryvar,
-       isintegervar
+       tunable_parameters, isirreducible, getdescription, hasdescription
 export ode_order_lowering, dae_order_lowering, liouville_transform
 export PDESystem
 export Differential, expand_derivatives, @derivatives
diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index a6bd0ce8c9..0f5e82aad6 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -120,9 +120,7 @@ function parse_variable_def!(dict, mod, arg, varclass, kwargs;
         (:misc, VariableMisc),
         (:disturbance, VariableDisturbance),
         (:tunable, VariableTunable),
-        (:dist, VariableDistribution),
-        (:binary, VariableBinary),
-        (:integer, VariableInteger)]
+        (:dist, VariableDistribution)]
 
     arg isa LineNumberNode && return
     MLStyle.@match arg begin
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 12b44074e6..71b37ea978 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -258,8 +258,9 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
     if isnothing(lb) && isnothing(ub) # use the symbolically specified bounds
         lb = first.(getbounds.(dvs))
         ub = last.(getbounds.(dvs))
-        lb[isbinaryvar.(dvs)] .= 0
-        ub[isbinaryvar.(dvs)] .= 1
+        isboolean = symtype.(unwrap.(dvs)) .<: Bool
+        lb[isboolean] .= 0
+        ub[isboolean] .= 1
     else # use the user supplied variable bounds
         xor(isnothing(lb), isnothing(ub)) &&
             throw(ArgumentError("Expected both `lb` and `ub` to be supplied"))
@@ -269,7 +270,7 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
             throw(ArgumentError("Expected both `ub` to be of the same length as the vector of optimization variables"))
     end
 
-    int = isintegervar.(dvs) .| isbinaryvar.(dvs)
+    int = symtype.(unwrap.(dvs)) .<: Integer
 
     defs = defaults(sys)
     defs = mergedefaults(defs, parammap, ps)
@@ -476,8 +477,9 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0map,
     if isnothing(lb) && isnothing(ub) # use the symbolically specified bounds
         lb = first.(getbounds.(dvs))
         ub = last.(getbounds.(dvs))
-        lb[isbinaryvar.(dvs)] .= 0
-        ub[isbinaryvar.(dvs)] .= 1
+        isboolean = symtype.(unwrap.(dvs)) .<: Bool
+        lb[isboolean] .= 0
+        ub[isboolean] .= 1
     else # use the user supplied variable bounds
         xor(isnothing(lb), isnothing(ub)) &&
             throw(ArgumentError("Expected both `lb` and `ub` to be supplied"))
@@ -487,7 +489,7 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0map,
             throw(ArgumentError("Expected `ub` to be of the same length as the vector of optimization variables"))
     end
 
-    int = isintegervar.(dvs) .| isbinaryvar.(dvs)
+    int = symtype.(unwrap.(dvs)) .<: Integer
 
     defs = defaults(sys)
     defs = mergedefaults(defs, parammap, ps)
diff --git a/src/variables.jl b/src/variables.jl
index 52db3339ef..9a9dbe4364 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -439,40 +439,6 @@ function hasdescription(x)
     getdescription(x) != ""
 end
 
-## binary variables =================================================================
-struct VariableBinary end
-Symbolics.option_to_metadata_type(::Val{:binary}) = VariableBinary
-
-isbinaryvar(x::Num) = isbinaryvar(Symbolics.unwrap(x))
-
-"""
-    isbinaryvar(x)
-
-Determine if a variable is binary.
-"""
-function isbinaryvar(x)
-    p = Symbolics.getparent(x, nothing)
-    p === nothing || (x = p)
-    return Symbolics.getmetadata(x, VariableBinary, false)
-end
-
-## integer variables =================================================================
-struct VariableInteger end
-Symbolics.option_to_metadata_type(::Val{:integer}) = VariableInteger
-
-isintegervar(x::Num) = isintegervar(Symbolics.unwrap(x))
-
-"""
-    isintegervar(x)
-
-Determine if a variable is an integer.
-"""
-function isintegervar(x)
-    p = Symbolics.getparent(x, nothing)
-    p === nothing || (x = p)
-    return Symbolics.getmetadata(x, VariableInteger, false)
-end
-
 ## Brownian
 """
     tobrownian(s::Sym)
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 46f455497e..70b90b669f 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -273,14 +273,14 @@ end
 
 @testset "Metadata in variables" begin
     metadata = Dict(:description => "Variable to test metadata in the Model.structure",
-        :input => true, :bounds => (-1, 1), :connection_type => :Flow, :integer => true,
-        :binary => false, :tunable => false, :disturbance => true, :dist => Normal(1, 1))
+        :input => true, :bounds => (-1, 1), :connection_type => :Flow,
+        :tunable => false, :disturbance => true, :dist => Normal(1, 1))
 
     @connector MockMeta begin
         m(t),
         [description = "Variable to test metadata in the Model.structure",
-            input = true, bounds = (-1, 1), connect = Flow, integer = true,
-            binary = false, tunable = false, disturbance = true, dist = Normal(1, 1)]
+            input = true, bounds = (-1, 1), connect = Flow,
+            tunable = false, disturbance = true, dist = Normal(1, 1)]
     end
 
     for (k, v) in metadata
diff --git a/test/test_variable_metadata.jl b/test/test_variable_metadata.jl
index b79df5e72d..cdaeae8b93 100644
--- a/test/test_variable_metadata.jl
+++ b/test/test_variable_metadata.jl
@@ -112,19 +112,3 @@ sp = Set(p)
 @named sys = ODESystem([u ~ p], t)
 
 @test_nowarn show(stdout, "text/plain", sys)
-
-@testset "binary" begin
-    @parameters t
-    @variables u(t) [binary = true]
-    @parameters p [binary = true]
-    @test isbinaryvar(u)
-    @test isbinaryvar(p)
-end
-
-@testset "integer" begin
-    @parameters t
-    @variables u(t) [integer = true]
-    @parameters p [integer = true]
-    @test isintegervar(u)
-    @test isintegervar(p)
-end

From 3953073a5d9e8f913ae0bf2185016c6f7b136db4 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Fri, 23 Feb 2024 20:22:29 +0530
Subject: [PATCH 12/18] docs: disable OptimizationSystem doc examples

---
 docs/src/tutorials/optimization.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md
index 3aa4a25596..9f71b592d8 100644
--- a/docs/src/tutorials/optimization.md
+++ b/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

From 7aece73e94d1eb95cbbe82aa8e9c1a6d2ccf85ba Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 26 Feb 2024 11:13:38 +0530
Subject: [PATCH 13/18] docs: fix various example block errors

---
 docs/Project.toml                                 | 2 ++
 docs/src/basics/Validation.md                     | 9 +++++++--
 docs/src/basics/Variable_metadata.md              | 2 +-
 docs/src/examples/parsing.md                      | 5 +++--
 docs/src/tutorials/ode_modeling.md                | 2 +-
 docs/src/tutorials/programmatically_generating.md | 2 +-
 docs/src/tutorials/stochastic_diffeq.md           | 2 +-
 7 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 4c90d35b77..8aaabc6054 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -14,6 +14,7 @@ OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
@@ -33,6 +34,7 @@ OptimizationOptimJL = "0.1"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
 StochasticDiffEq = "6"
+SymbolicIndexingInterface = "0.3.1"
 SymbolicUtils = "1"
 Symbolics = "5"
 Unitful = "1.12"
diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md
index ada0e57810..c9c662f13b 100644
--- a/docs/src/basics/Validation.md
+++ b/docs/src/basics/Validation.md
@@ -59,7 +59,11 @@ 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.
@@ -100,7 +104,8 @@ 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)
 ```
 
diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/basics/Variable_metadata.md
index c2c1f6d4ce..5e25f33373 100644
--- a/docs/src/basics/Variable_metadata.md
+++ b/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)
 ```
diff --git a/docs/src/examples/parsing.md b/docs/src/examples/parsing.md
index 78ad824d6e..66a4e4d82f 100644
--- a/docs/src/examples/parsing.md
+++ b/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())
 ```
diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index e749f42362..ef354071d9 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/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
 ```
 
diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md
index 079f50632d..8b8b03027a 100644
--- a/docs/src/tutorials/programmatically_generating.md
+++ b/docs/src/tutorials/programmatically_generating.md
@@ -41,7 +41,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 eqs = [D(x) ~ (h - x) / τ] # create an array of equations
 
 # your first ODE, consisting of a single equation, indicated by ~
-@named fol_model = ODESystem(eqs, t)
+@named model = ODESystem(eqs, t)
 
 # Perform the standard transformations and mark the model complete
 # Note: Complete models cannot be subsystems of other models!
diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md
index 85c2061af0..fd03d51810 100644
--- a/docs/src/tutorials/stochastic_diffeq.md
+++ b/docs/src/tutorials/stochastic_diffeq.md
@@ -34,5 +34,5 @@ parammap = [
 ]
 
 prob = SDEProblem(de, u0map, (0.0, 100.0), parammap)
-sol = solve(prob, LambdaEulerHeun())
+sol = solve(prob, LambaEulerHeun())
 ```

From 6c5d0809b35adc20cb15d629af3338d1d7161aac Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 26 Feb 2024 16:07:51 +0530
Subject: [PATCH 14/18] feat: add generate_custom_function

---
 src/ModelingToolkit.jl                   |   2 +-
 src/systems/abstractsystem.jl            | 156 ++++++++++++++++++++++-
 src/systems/diffeqs/abstractodesystem.jl |  22 +---
 test/generate_custom_function.jl         |  53 ++++++++
 test/runtests.jl                         |   1 +
 5 files changed, 211 insertions(+), 23 deletions(-)
 create mode 100644 test/generate_custom_function.jl

diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 7e0964b541..2099c06c2e 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -233,7 +233,7 @@ export independent_variable, equations, controls,
        observed, full_equations
 export structural_simplify, expand_connections, linearize, linearization_function
 
-export calculate_jacobian, generate_jacobian, generate_function
+export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function
 export calculate_control_jacobian, generate_control_jacobian
 export calculate_tgrad, generate_tgrad
 export calculate_gradient, generate_gradient
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index afc23b7cd3..a1a98317e2 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -82,7 +82,7 @@ function calculate_hessian end
 
 """
 ```julia
-generate_tgrad(sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_tgrad(sys::AbstractTimeDependentSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                expression = Val{true}; kwargs...)
 ```
 
@@ -93,7 +93,7 @@ function generate_tgrad end
 
 """
 ```julia
-generate_gradient(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_gradient(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                   expression = Val{true}; kwargs...)
 ```
 
@@ -104,7 +104,7 @@ function generate_gradient end
 
 """
 ```julia
-generate_jacobian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_jacobian(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                   expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -115,7 +115,7 @@ function generate_jacobian end
 
 """
 ```julia
-generate_factorized_W(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_factorized_W(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                       expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -126,7 +126,7 @@ function generate_factorized_W end
 
 """
 ```julia
-generate_hessian(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_hessian(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                  expression = Val{true}; sparse = false, kwargs...)
 ```
 
@@ -137,7 +137,7 @@ function generate_hessian end
 
 """
 ```julia
-generate_function(sys::AbstractSystem, dvs = unknowns(sys), ps = parameters(sys),
+generate_function(sys::AbstractSystem, dvs = unknowns(sys), ps = full_parameters(sys),
                   expression = Val{true}; kwargs...)
 ```
 
@@ -145,6 +145,150 @@ Generate a function to evaluate the system's equations.
 """
 function generate_function end
 
+function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
+        ps = parameters(sys); wrap_code = nothing, kwargs...)
+    p = reorder_parameters(sys, ps)
+    isscalar = !(exprs isa AbstractArray)
+    if wrap_code === nothing
+        wrap_code = isscalar ? identity : (identity, identity)
+    end
+    pre, sol_states = get_substitutions_and_solved_unknowns(sys)
+
+    if is_time_dependent(sys)
+        return build_function(exprs,
+            dvs,
+            p...,
+            get_iv(sys);
+            kwargs...,
+            postprocess_fbody = pre,
+            states = sol_states,
+            wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
+                        wrap_array_vars(sys, exprs; dvs)
+        )
+    else
+        return build_function(exprs,
+            dvs,
+            p...;
+            kwargs...,
+            postprocess_fbody = pre,
+            states = sol_states,
+            wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘
+                        wrap_array_vars(sys, exprs; dvs)
+        )
+    end
+end
+
+function wrap_array_vars(sys::AbstractSystem, exprs; dvs = unknowns(sys))
+    isscalar = !(exprs isa AbstractArray)
+    allvars = if isscalar
+        Set(get_variables(exprs))
+    else
+        union(get_variables.(exprs)...)
+    end
+    array_vars = Dict{Any, AbstractArray{Int}}()
+    for (j, x) in enumerate(dvs)
+        if istree(x) && operation(x) == getindex
+            arg = arguments(x)[1]
+            arg in allvars || continue
+            inds = get!(() -> Int[], array_vars, arg)
+            push!(inds, j)
+        end
+    end
+    for (k, inds) in array_vars
+        if inds == (inds′ = inds[1]:inds[end])
+            array_vars[k] = inds′
+        end
+    end
+    if isscalar
+        function (expr)
+            Func(
+                expr.args,
+                [],
+                Let(
+                    [k ← :(view($(expr.args[1].name), $v)) for (k, v) in array_vars],
+                    expr.body,
+                    false
+                )
+            )
+        end
+    else
+        function (expr)
+            Func(
+                expr.args,
+                [],
+                Let(
+                    [k ← :(view($(expr.args[1].name), $v)) for (k, v) in array_vars],
+                    expr.body,
+                    false
+                )
+            )
+        end,
+        function (expr)
+            Func(
+                expr.args,
+                [],
+                Let(
+                    [k ← :(view($(expr.args[2].name), $v)) for (k, v) in array_vars],
+                    expr.body,
+                    false
+                )
+            )
+        end
+    end
+end
+
+function wrap_mtkparameters(sys::AbstractSystem, isscalar::Bool)
+    if has_index_cache(sys) && get_index_cache(sys) !== nothing
+        offset = Int(is_time_dependent(sys))
+
+        if isscalar
+            function (expr)
+                p = gensym(:p)
+                Func(
+                    [
+                        expr.args[1],
+                        DestructuredArgs(
+                            [arg.name for arg in expr.args[2:(end - offset)]], p),
+                        (isone(offset) ? (expr.args[end],) : ())...
+                    ],
+                    [],
+                    Let(expr.args[2:(end - offset)], expr.body, false)
+                )
+            end
+        else
+            function (expr)
+                p = gensym(:p)
+                Func(
+                    [
+                        expr.args[1],
+                        DestructuredArgs(
+                            [arg.name for arg in expr.args[2:(end - offset)]], p),
+                        (isone(offset) ? (expr.args[end],) : ())...
+                    ],
+                    [],
+                    Let(expr.args[2:(end - offset)], expr.body, false)
+                )
+            end,
+            function (expr)
+                p = gensym(:p)
+                Func(
+                    [
+                        expr.args[1],
+                        expr.args[2],
+                        DestructuredArgs(
+                            [arg.name for arg in expr.args[3:(end - offset)]], p),
+                        (isone(offset) ? (expr.args[end],) : ())...
+                    ],
+                    [],
+                    Let(expr.args[3:(end - offset)], expr.body, false)
+                )
+            end
+        end
+    else
+        identity
+    end
+end
+
 mutable struct Substitutions
     subs::Vector{Equation}
     deps::Vector{Vector{Int}}
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index a734845118..b6e8c50563 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -152,22 +152,8 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
         ddvs = implicit_dae ? map(Differential(get_iv(sys)), dvs) :
                nothing,
         isdde = false,
+        wrap_code = nothing,
         kwargs...)
-    array_vars = Dict{Any, Vector{Int}}()
-    for (j, x) in enumerate(dvs)
-        if istree(x) && operation(x) == getindex
-            arg = arguments(x)[1]
-            inds = get!(() -> Int[], array_vars, arg)
-            push!(inds, j)
-        end
-    end
-    subs = Dict()
-    for (k, inds) in array_vars
-        if inds == (inds′ = inds[1]:inds[end])
-            inds = inds′
-        end
-        subs[k] = term(view, Sym{Any}(Symbol("ˍ₋arg1")), inds)
-    end
     if isdde
         eqs = delay_to_function(sys)
     else
@@ -180,13 +166,15 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
     # substitute x(t) by just x
     rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
            [eq.rhs for eq in eqs]
-    rhss = fast_substitute(rhss, subs)
 
     # 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), reorder_parameters(sys, ps))
     t = get_iv(sys)
 
+    if wrap_code === nothing
+        wrap_code = (identity, identity)
+    end
     if isdde
         build_function(rhss, u, DDE_HISTORY_FUN, p..., t; kwargs...)
     else
@@ -195,10 +183,12 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
         if implicit_dae
             build_function(rhss, ddvs, u, p..., t; postprocess_fbody = pre,
                 states = sol_states,
+                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs),
                 kwargs...)
         else
             build_function(rhss, u, p..., t; postprocess_fbody = pre,
                 states = sol_states,
+                wrap_code = wrap_code .∘ wrap_array_vars(sys, rhss; dvs),
                 kwargs...)
         end
     end
diff --git a/test/generate_custom_function.jl b/test/generate_custom_function.jl
new file mode 100644
index 0000000000..c4ea7134f5
--- /dev/null
+++ b/test/generate_custom_function.jl
@@ -0,0 +1,53 @@
+using ModelingToolkit
+using ModelingToolkit: t_nounits as t, D_nounits as D
+using IfElse
+
+@variables x(t) y(t)[1:3]
+@parameters p1=1.0 p2[1:3]=[1.0, 2.0, 3.0] p3::Int=1 p4::Bool=false
+
+sys = complete(ODESystem(Equation[], t, [x; y], [p1, p2, p3, p4]; name = :sys))
+u0 = [1.0, 2.0, 3.0, 4.0]
+p = ModelingToolkit.MTKParameters(sys, [])
+
+fn1 = generate_custom_function(sys, x + y[1] + p1 + p2[1] + p3 * t; expression = Val(false))
+@test fn1(u0, p, 0.0) == 5.0
+
+fn2 = generate_custom_function(
+    sys, x + y[1] + p1 + p2[1] + p3 * t, [x], [p1, p2, p3]; expression = Val(false))
+@test fn1(u0, p, 0.0) == 5.0
+
+fn3_oop, fn3_iip = generate_custom_function(
+    sys, [x + y[2], y[3] + p2[2], p1 + p3, 3t]; expression = Val(false))
+
+buffer = zeros(4)
+fn3_iip(buffer, u0, p, 1.0)
+@test buffer == [4.0, 6.0, 2.0, 3.0]
+@test fn3_oop(u0, p, 1.0) == [4.0, 6.0, 2.0, 3.0]
+
+fn4 = generate_custom_function(sys, IfElse.ifelse(p4, p1, p2[2]); expression = Val(false))
+@test fn4(u0, p, 1.0) == 2.0
+fn5 = generate_custom_function(sys, IfElse.ifelse(!p4, p1, p2[2]); expression = Val(false))
+@test fn5(u0, p, 1.0) == 1.0
+
+@variables x y[1:3]
+sys = complete(NonlinearSystem(Equation[], [x; y], [p1, p2, p3, p4]; name = :sys))
+
+fn1 = generate_custom_function(sys, x + y[1] + p1 + p2[1] + p3; expression = Val(false))
+@test fn1(u0, p) == 6.0
+
+fn2 = generate_custom_function(
+    sys, x + y[1] + p1 + p2[1] + p3, [x], [p1, p2, p3]; expression = Val(false))
+@test fn1(u0, p) == 6.0
+
+fn3_oop, fn3_iip = generate_custom_function(
+    sys, [x + y[2], y[3] + p2[2], p1 + p3]; expression = Val(false))
+
+buffer = zeros(3)
+fn3_iip(buffer, u0, p)
+@test buffer == [4.0, 6.0, 2.0]
+@test fn3_oop(u0, p, 1.0) == [4.0, 6.0, 2.0]
+
+fn4 = generate_custom_function(sys, IfElse.ifelse(p4, p1, p2[2]); expression = Val(false))
+@test fn4(u0, p, 1.0) == 2.0
+fn5 = generate_custom_function(sys, IfElse.ifelse(!p4, p1, p2[2]); expression = Val(false))
+@test fn5(u0, p, 1.0) == 1.0
diff --git a/test/runtests.jl b/test/runtests.jl
index 610ed55354..96bdbbf06e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -63,6 +63,7 @@ end
         @safetestset "FuncAffect Test" include("funcaffect.jl")
         @safetestset "Constants Test" include("constants.jl")
         @safetestset "Parameter Dependency Test" include("parameter_dependencies.jl")
+        @safetestset "Generate Custom Function Test" include("generate_custom_function.jl")
     end
 
     if GROUP == "All" || GROUP == "InterfaceII"

From 9b6f8618a8dbbad3acfe3ffeca6843d9b2e31e22 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 26 Feb 2024 16:08:01 +0530
Subject: [PATCH 15/18] test: remove ControlSystemsMTK from test dependencies

---
 Project.toml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index b1f9d68908..4085f91af9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -116,7 +116,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 +138,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"]

From f5848156f74a7058f838d5cef8dfd61050f0151d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 26 Feb 2024 18:32:37 +0530
Subject: [PATCH 16/18] feat: support indexing with `Symbol`s in IndexCache

---
 src/systems/abstractsystem.jl | 19 ++++++++++++++++---
 src/systems/index_cache.jl    | 18 ++++++++++++++----
 2 files changed, 30 insertions(+), 7 deletions(-)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index a1a98317e2..3adaf7d528 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -346,6 +346,9 @@ function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return haskey(ic.unknown_idx, hash(sym))
+    end
     return any(isequal(sym), getname.(variable_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
            count(isequal(sym), Symbol.(nameof(sys), :₊, getname.(variable_symbols(sys)))) ==
@@ -377,6 +380,9 @@ function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return get(ic.unknown_idx, h, nothing)
+    end
     idx = findfirst(isequal(sym), getname.(variable_symbols(sys)))
     if idx !== nothing
         return idx
@@ -401,12 +407,14 @@ function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym)
         ic = get_index_cache(sys)
         h = getsymbolhash(sym)
         return if haskey(ic.param_idx, h) || haskey(ic.discrete_idx, h) ||
-                  haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h)
+                  haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h) ||
+                  haskey(ic.nonnumeric_idx, h)
             true
         else
             h = getsymbolhash(default_toterm(sym))
             haskey(ic.param_idx, h) || haskey(ic.discrete_idx, h) ||
-                haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h)
+                haskey(ic.constant_idx, h) || haskey(ic.dependent_idx, h) ||
+                haskey(ic.nonnumeric_idx, h)
         end
     end
     return any(isequal(sym), parameter_symbols(sys)) ||
@@ -414,6 +422,9 @@ function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return ParameterIndex(ic, sym) !== nothing
+    end
     return any(isequal(sym), getname.(parameter_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
            count(isequal(sym),
@@ -426,7 +437,6 @@ function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym)
     end
     if has_index_cache(sys) && get_index_cache(sys) !== nothing
         ic = get_index_cache(sys)
-        h = getsymbolhash(sym)
         return if (idx = ParameterIndex(ic, sym)) !== nothing
             idx
         elseif (idx = ParameterIndex(ic, default_toterm(sym))) !== nothing
@@ -444,6 +454,9 @@ function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym)
 end
 
 function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym::Symbol)
+    if has_index_cache(sys) && (ic = get_index_cache(sys)) !== nothing
+        return ParameterIndex(ic, sym)
+    end
     idx = findfirst(isequal(sym), getname.(parameter_symbols(sys)))
     if idx !== nothing
         return idx
diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl
index 9228320c36..a7e6c8b157 100644
--- a/src/systems/index_cache.jl
+++ b/src/systems/index_cache.jl
@@ -40,10 +40,16 @@ function IndexCache(sys::AbstractSystem)
     let idx = 1
         for sym in unks
             h = getsymbolhash(sym)
-            if Symbolics.isarraysymbolic(sym)
-                unk_idxs[h] = idx:(idx + length(sym) - 1)
+            sym_idx = if Symbolics.isarraysymbolic(sym)
+                idx:(idx + length(sym) - 1)
             else
-                unk_idxs[h] = idx
+                idx
+            end
+            unk_idxs[h] = sym_idx
+
+            if hasname(sym)
+                h = hash(getname(sym))
+                unk_idxs[h] = sym_idx
             end
             idx += length(sym)
         end
@@ -122,6 +128,10 @@ function IndexCache(sys::AbstractSystem)
                 idxs[h] = (i, j)
                 h = getsymbolhash(default_toterm(p))
                 idxs[h] = (i, j)
+                if hasname(p)
+                    h = hash(getname(p))
+                    idxs[h] = (i, j)
+                end
             end
             push!(buffer_sizes, BufferTemplate(T, length(buf)))
         end
@@ -151,7 +161,7 @@ end
 
 function ParameterIndex(ic::IndexCache, p, sub_idx = ())
     p = unwrap(p)
-    h = getsymbolhash(p)
+    h = p isa Symbol ? hash(p) : getsymbolhash(p)
     return if haskey(ic.param_idx, h)
         ParameterIndex(SciMLStructures.Tunable(), (ic.param_idx[h]..., sub_idx...))
     elseif haskey(ic.discrete_idx, h)

From 87df269cbcb7a51dc170669e565ba1f30f2918af Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 26 Feb 2024 20:08:18 +0530
Subject: [PATCH 17/18] docs: add doc for `generate_custom_function`

---
 docs/src/internals.md         |  2 +-
 src/systems/abstractsystem.jl | 18 ++++++++++++++++++
 2 files changed, 19 insertions(+), 1 deletion(-)

diff --git a/docs/src/internals.md b/docs/src/internals.md
index 79ded5c57f..a0d9ebfdf6 100644
--- a/docs/src/internals.md
+++ b/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`.
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 3adaf7d528..f95cb1fe6f 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -145,8 +145,26 @@ Generate a function to evaluate the system's equations.
 """
 function generate_function end
 
+"""
+```julia
+generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
+                         ps = full_parameters(sys); kwargs...)
+```
+
+Generate a function to evaluate `exprs`. `exprs` is a symbolic expression or
+array of symbolic expression involving symbolic variables in `sys`. The symbolic variables
+may be subsetted using `dvs` and `ps`. All `kwargs` except `postprocess_fbody` and `states`
+are passed to the internal [`build_function`](@ref) call. The returned function can be called
+as `f(u, p, t)` or `f(du, u, p, t)` for time-dependent systems and `f(u, p)` or `f(du, u, p)`
+for time-independent systems. If `split=true` (the default) was passed to [`complete`](@ref),
+[`structural_simplify`](@ref) or [`@mtkbuild`](@ref), `p` is expected to be an `MTKParameters`
+object.
+"""
 function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys),
         ps = parameters(sys); wrap_code = nothing, kwargs...)
+    if !iscomplete(sys)
+        error("A completed system is required. Call `complete` or `structural_simplify` on the system.")
+    end
     p = reorder_parameters(sys, ps)
     isscalar = !(exprs isa AbstractArray)
     if wrap_code === nothing

From 9535ba7b607cbc886fba10dcf7cad7cdcaf434b9 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@juliahub.com>
Date: Mon, 26 Feb 2024 21:25:59 +0530
Subject: [PATCH 18/18] feat: remove `IfElse` support

---
 NEWS.md                          |  1 +
 Project.toml                     |  2 --
 docs/src/basics/FAQ.md           | 11 +++++------
 src/ModelingToolkit.jl           |  1 -
 src/systems/unit_check.jl        |  2 +-
 src/systems/validation.jl        |  2 +-
 test/domain_connectors.jl        |  1 -
 test/downstream/linearize.jl     |  2 +-
 test/dq_units.jl                 |  2 +-
 test/generate_custom_function.jl |  9 ++++-----
 test/state_selection.jl          |  6 +++---
 test/units.jl                    |  4 ++--
 12 files changed, 19 insertions(+), 24 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 94b8c507ca..1a1c6c4126 100644
--- a/NEWS.md
+++ b/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.
diff --git a/Project.toml b/Project.toml
index 4085f91af9..aa0466e91d 100644
--- a/Project.toml
+++ b/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"
diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md
index 19fcf68080..9d12b619fb 100644
--- a/docs/src/basics/FAQ.md
+++ b/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?
 
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 2099c06c2e..b70074f6fa 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -27,7 +27,6 @@ using PrecompileTools, Reexport
     using DocStringExtensions
     using Base: RefValue
     using Combinatorics
-    import IfElse
     import Distributions
     import FunctionWrappersWrappers
     using URIs: URI
diff --git a/src/systems/unit_check.jl b/src/systems/unit_check.jl
index 6ff59f08fa..677d29895f 100644
--- a/src/systems/unit_check.jl
+++ b/src/systems/unit_check.jl
@@ -1,5 +1,5 @@
 #For dispatching get_unit
-const Conditional = Union{typeof(ifelse), typeof(IfElse.ifelse)}
+const Conditional = Union{typeof(ifelse)}
 const Comparison = Union{typeof.([==, !=, ≠, <, <=, ≤, >, >=, ≥])...}
 
 struct ValidationError <: Exception
diff --git a/src/systems/validation.jl b/src/systems/validation.jl
index f0c5fa95fd..ea458d8b7e 100644
--- a/src/systems/validation.jl
+++ b/src/systems/validation.jl
@@ -1,6 +1,6 @@
 module UnitfulUnitCheck
 
-using ..ModelingToolkit, Symbolics, SciMLBase, Unitful, IfElse, RecursiveArrayTools
+using ..ModelingToolkit, Symbolics, SciMLBase, Unitful, RecursiveArrayTools
 using ..ModelingToolkit: ValidationError,
                          ModelingToolkit, Connection, instream, JumpType, VariableUnit,
                          get_systems,
diff --git a/test/domain_connectors.jl b/test/domain_connectors.jl
index 572f440681..9a43d2938f 100644
--- a/test/domain_connectors.jl
+++ b/test/domain_connectors.jl
@@ -1,7 +1,6 @@
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using Test
-using IfElse: ifelse
 
 @connector function HydraulicPort(; p_int, name)
     pars = @parameters begin
diff --git a/test/downstream/linearize.jl b/test/downstream/linearize.jl
index 302fabed94..c1aa6f6724 100644
--- a/test/downstream/linearize.jl
+++ b/test/downstream/linearize.jl
@@ -163,7 +163,7 @@ end
 function saturation(; y_max, y_min = y_max > 0 ? -y_max : -Inf, name)
     @variables u(t)=0 y(t)=0
     @parameters y_max=y_max y_min=y_min
-    ie = ModelingToolkit.IfElse.ifelse
+    ie = ifelse
     eqs = [
     # The equation below is equivalent to y ~ clamp(u, y_min, y_max)
         y ~ ie(u > y_max, y_max, ie((y_min < u) & (u < y_max), u, y_min))
diff --git a/test/dq_units.jl b/test/dq_units.jl
index a4bf6cba9b..75bbd4c4a9 100644
--- a/test/dq_units.jl
+++ b/test/dq_units.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, IfElse, DynamicQuantities
+using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, DynamicQuantities
 using Test
 MT = ModelingToolkit
 using ModelingToolkit: t, D
diff --git a/test/generate_custom_function.jl b/test/generate_custom_function.jl
index c4ea7134f5..9b64b20c12 100644
--- a/test/generate_custom_function.jl
+++ b/test/generate_custom_function.jl
@@ -1,6 +1,5 @@
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
-using IfElse
 
 @variables x(t) y(t)[1:3]
 @parameters p1=1.0 p2[1:3]=[1.0, 2.0, 3.0] p3::Int=1 p4::Bool=false
@@ -24,9 +23,9 @@ fn3_iip(buffer, u0, p, 1.0)
 @test buffer == [4.0, 6.0, 2.0, 3.0]
 @test fn3_oop(u0, p, 1.0) == [4.0, 6.0, 2.0, 3.0]
 
-fn4 = generate_custom_function(sys, IfElse.ifelse(p4, p1, p2[2]); expression = Val(false))
+fn4 = generate_custom_function(sys, ifelse(p4, p1, p2[2]); expression = Val(false))
 @test fn4(u0, p, 1.0) == 2.0
-fn5 = generate_custom_function(sys, IfElse.ifelse(!p4, p1, p2[2]); expression = Val(false))
+fn5 = generate_custom_function(sys, ifelse(!p4, p1, p2[2]); expression = Val(false))
 @test fn5(u0, p, 1.0) == 1.0
 
 @variables x y[1:3]
@@ -47,7 +46,7 @@ fn3_iip(buffer, u0, p)
 @test buffer == [4.0, 6.0, 2.0]
 @test fn3_oop(u0, p, 1.0) == [4.0, 6.0, 2.0]
 
-fn4 = generate_custom_function(sys, IfElse.ifelse(p4, p1, p2[2]); expression = Val(false))
+fn4 = generate_custom_function(sys, ifelse(p4, p1, p2[2]); expression = Val(false))
 @test fn4(u0, p, 1.0) == 2.0
-fn5 = generate_custom_function(sys, IfElse.ifelse(!p4, p1, p2[2]); expression = Val(false))
+fn5 = generate_custom_function(sys, ifelse(!p4, p1, p2[2]); expression = Val(false))
 @test fn5(u0, p, 1.0) == 1.0
diff --git a/test/state_selection.jl b/test/state_selection.jl
index 34e4767460..c1e4328715 100644
--- a/test/state_selection.jl
+++ b/test/state_selection.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, IfElse, Test
+using ModelingToolkit, OrdinaryDiffEq, Test
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 sts = @variables x1(t) x2(t) x3(t) x4(t)
@@ -249,7 +249,7 @@ let
     Δp2 = p2
 
     eqs = [+flow(xm, Δp1) ~ rho1 * dV1 + drho1 * V1
-           0 ~ IfElse.ifelse(w > 0.5,
+           0 ~ ifelse(w > 0.5,
                (0) - (rho2 * dV2 + drho2 * V2),
                (-flow(xm, Δp2)) - (rho2 * dV2 + drho2 * V2))
            V1 ~ (l_1f + w) * A_1f
@@ -265,7 +265,7 @@ let
            D(w) ~ dw
            D(dw) ~ ddw
            xf ~ 20e-3 * (1 - cos(2 * π * 5 * t))
-           0 ~ IfElse.ifelse(w > 0.5,
+           0 ~ ifelse(w > 0.5,
                (m_total * ddw) - (p1 * A_1f - p2 * A_2f - damp * dw),
                (m_total * ddw) - (p1 * A_1f - p2 * A_2f))]
     # ----------------------------------------------------------------------------
diff --git a/test/units.jl b/test/units.jl
index fd737cf816..e170540270 100644
--- a/test/units.jl
+++ b/test/units.jl
@@ -1,4 +1,4 @@
-using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, IfElse, Unitful
+using ModelingToolkit, OrdinaryDiffEq, JumpProcesses, Unitful
 using Test
 MT = ModelingToolkit
 UMT = ModelingToolkit.UnitfulUnitCheck
@@ -33,7 +33,7 @@ D = Differential(t)
 @test UMT.get_unit(t / τ) == UMT.unitless
 @test UMT.equivalent(UMT.get_unit(P - E / τ), u"MW")
 @test UMT.equivalent(UMT.get_unit(D(D(E))), u"MW/ms")
-@test UMT.get_unit(IfElse.ifelse(t > t, P, E / τ)) == u"MW"
+@test UMT.get_unit(ifelse(t > t, P, E / τ)) == u"MW"
 @test UMT.get_unit(1.0^(t / τ)) == UMT.unitless
 @test UMT.get_unit(exp(t / τ)) == UMT.unitless
 @test UMT.get_unit(sin(t / τ)) == UMT.unitless