From e56de20aa260ffbd29bbacddcb98389c595996f7 Mon Sep 17 00:00:00 2001
From: Brad Carman <brad_carman@instron.com>
Date: Mon, 25 Sep 2023 10:18:54 -0400
Subject: [PATCH 01/50] split parameters linearize bug

---
 test/split_parameters.jl | 73 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 73 insertions(+)

diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index ef8f434dca..11608e1a7f 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -77,3 +77,76 @@ prob = ODEProblem(sys, [], tspan, []; tofloat = false)
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int64}}
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
+
+
+
+
+
+
+# ------------------------- Bug
+using ModelingToolkit, LinearAlgebra
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkitStandardLibrary.Blocks: t
+using ModelingToolkit: connect
+
+"A wrapper function to make symbolic indexing easier"
+function wr(sys)
+    ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems=[sys], name=:a_wrapper)
+end
+indexof(sym,syms) = findfirst(isequal(sym),syms)
+
+# Parameters
+m1 = 1.
+m2 = 1.
+k = 10. # Spring stiffness
+c = 3.  # Damping coefficient
+
+@named inertia1 = Inertia(; J = m1)
+@named inertia2 = Inertia(; J = m2)
+@named spring = Spring(; c = k)
+@named damper = Damper(; d = c)
+@named torque = Torque(use_support=false)
+
+function SystemModel(u=nothing; name=:model)
+    eqs = [
+        connect(torque.flange, inertia1.flange_a)
+        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
+        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
+    ]
+    if u !== nothing
+        push!(eqs, connect(torque.tau, u.output))
+        return @named model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper, u])
+    end
+    ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
+end
+
+
+model = SystemModel() # Model with load disturbance
+@named d = Step(start_time=1., duration=10., offset=0., height=1.) # Disturbance
+model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi] # This is the state realization we want to control
+inputs = [model.torque.tau.u]
+matrices, ssys = ModelingToolkit.linearize(wr(model), inputs, model_outputs)
+
+# Design state-feedback gain using LQR
+# Define cost matrices
+x_costs = [
+    model.inertia1.w =>   1.
+    model.inertia2.w =>   1.
+    model.inertia1.phi => 1.
+    model.inertia2.phi => 1.
+]
+L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller
+@named state_feedback = MatrixGain(K=-L) # Build negative feedback into the feedback matrix
+@named add = Add(;k1=1., k2=1.) # To add the control signal and the disturbance
+
+connections = [
+    [state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4]
+    connect(d.output, :d, add.input1)
+    connect(add.input2, state_feedback.output)
+    connect(add.output, :u, model.torque.tau)
+]
+closed_loop = ODESystem(connections, t, systems=[model, state_feedback, add, d], name=:closed_loop)
+S = get_sensitivity(closed_loop, :u)
+
+

From cee69ed07ad7b04dbb5977b9053b59d59bb86fed Mon Sep 17 00:00:00 2001
From: Brad Carman <brad_carman@instron.com>
Date: Mon, 25 Sep 2023 13:31:30 -0400
Subject: [PATCH 02/50] added old MatrixGain definition for comparison

---
 test/split_parameters.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index 11608e1a7f..3ce5f7e2f0 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -137,6 +137,16 @@ x_costs = [
     model.inertia2.phi => 1.
 ]
 L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller
+
+# This old definition of MatrixGain will work because the parameter space does not include K (an Array term)
+# @component function MatrixGainAlt(K::AbstractArray; name)
+#     nout, nin = size(K, 1), size(K, 2)
+#     @named input = RealInput(; nin = nin)
+#     @named output = RealOutput(; nout = nout)
+#     eqs = [output.u[i] ~ sum(K[i, j] * input.u[j] for j in 1:nin) for i in 1:nout]
+#     compose(ODESystem(eqs, t, [], []; name = name), [input, output])
+# end
+
 @named state_feedback = MatrixGain(K=-L) # Build negative feedback into the feedback matrix
 @named add = Add(;k1=1., k2=1.) # To add the control signal and the disturbance
 

From b4f821eeb21b4e5ac6f92a0c522e0abd3362544d Mon Sep 17 00:00:00 2001
From: Brad Carman <brad_carman@instron.com>
Date: Mon, 25 Sep 2023 19:14:07 -0400
Subject: [PATCH 03/50] new promote_to_concrete

---
 src/utils.jl             | 56 ++++++++++++++++++----------------------
 test/split_parameters.jl | 28 ++++++++++++++++++++
 2 files changed, 53 insertions(+), 31 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index abea65ab21..452323c6ff 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -657,46 +657,40 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
     end
     T = eltype(vs)
     if Base.isconcretetype(T) && (!tofloat || T === float(T)) # nothing to do
-        vs
+        return vs
     else
         sym_vs = filter(x -> SymbolicUtils.issym(x) || SymbolicUtils.istree(x), vs)
         isempty(sym_vs) || throw_missingvars_in_sys(sym_vs)
-        C = typeof(first(vs))
-        I = Int8
-        has_int = false
-        has_array = false
-        has_bool = false
-        array_T = nothing
+        
+        C = nothing
         for v in vs
-            if v isa AbstractArray
-                has_array = true
-                array_T = typeof(v)
-            end
-            E = eltype(v)
-            C = promote_type(C, E)
-            if E <: Integer
-                has_int = true
-                I = promote_type(I, E)
+            E = typeof(v)
+            if E <: Number
+                if tofloat
+                    E = float(E)
+                end
             end
-            if E <: Bool
-                has_bool = true
+            if use_union
+                if C === nothing
+                    C = E
+                else
+                    C = Union{C, E}
+                end
+            else
+                C = E
             end
         end
-        if tofloat && !has_array
-            C = float(C)
-        elseif has_array || (use_union && has_int && C !== I)
-            if has_array
-                C = Union{C, array_T}
-            end
-            if has_int
-                C = Union{C, I}
-            end
-            if has_bool
-                C = Union{C, Bool}
+
+        y = similar(vs, C)
+        for i in eachindex(vs)
+            if (vs[i] isa Number) & tofloat
+                y[i] = float(vs[i])
+            else
+                y[i] = vs[i]
             end
-            return copyto!(similar(vs, C), vs)
         end
-        convert.(C, vs)
+
+        return y
     end
 end
 
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index 3ce5f7e2f0..1ea7161ee0 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -2,6 +2,34 @@ using ModelingToolkit, Test
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq
 
+
+x = [1, 2.0, false, [1,2,3], Parameter(1.0)]
+
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Union{Float64, Parameter{Float64}, Vector{Int64}}
+
+y = ModelingToolkit.promote_to_concrete(x; tofloat=false)
+@test eltype(y) == Union{Bool, Float64, Int64, Parameter{Float64}, Vector{Int64}}
+
+
+x = [1, 2.0, false, [1,2,3]]
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Union{Float64, Vector{Int64}}
+
+x = Any[1, 2.0, false]
+y = ModelingToolkit.promote_to_concrete(x; tofloat=false)
+@test eltype(y) == Union{Bool, Float64, Int64}
+
+y = ModelingToolkit.promote_to_concrete(x; use_union=false)
+@test eltype(y) == Float64
+
+x = Float16[1., 2., 3.]
+y = ModelingToolkit.promote_to_concrete(x)
+@test eltype(y) == Float16
+
+
+
+
 # ------------------------ Mixed Single Values and Vector
 
 dt = 4e-4

From e62e3ee4562414a662cb59a65d0bf90be2ad4c26 Mon Sep 17 00:00:00 2001
From: Brad Carman <brad_carman@instron.com>
Date: Tue, 26 Sep 2023 05:58:01 -0400
Subject: [PATCH 04/50] improvements to protmote_to_concrete

---
 src/utils.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index 452323c6ff..0a6c5089a3 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -670,13 +670,13 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
                     E = float(E)
                 end
             end
+            if C === nothing
+                C = E
+            end
             if use_union
-                if C === nothing
-                    C = E
-                else
-                    C = Union{C, E}
-                end
+                C = Union{C, E}
             else
+                @assert C == E "`promote_to_concrete` can't make type $E uniform with $C"
                 C = E
             end
         end
@@ -684,9 +684,9 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
         y = similar(vs, C)
         for i in eachindex(vs)
             if (vs[i] isa Number) & tofloat
-                y[i] = float(vs[i])
+                y[i] = float(vs[i]) #needed because copyto! can't convert Int to Float automatically
             else
-                y[i] = vs[i]
+                y[i] = vs[i] 
             end
         end
 

From ee4deb98dcc7cf83d3ccb6d4cf0edb544614af92 Mon Sep 17 00:00:00 2001
From: Brad Carman <brad_carman@instron.com>
Date: Tue, 26 Sep 2023 05:58:14 -0400
Subject: [PATCH 05/50] remake issue

---
 test/split_parameters.jl | 18 +++++++++++++++++-
 1 file changed, 17 insertions(+), 1 deletion(-)

diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index 1ea7161ee0..3f231e2dd7 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -74,12 +74,28 @@ eqs = [y ~ src.output.u
 @named sys = ODESystem(eqs, t, vars, []; systems = [int, src])
 s = complete(sys)
 sys = structural_simplify(sys)
-prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x])
+prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat=false)
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int}, Vector{Vector{Float64}}}
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
 @test sol[y][end] == x[end]
 
+#TODO: remake becomes more complicated now, how to improve?
+defs = ModelingToolkit.defaults(sys)
+defs[s.src.data] = 2x
+p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)
+p′, = ModelingToolkit.split_parameters_by_type(p′) #NOTE: we need to ensure this is called now before calling remake()
+prob′ = remake(prob; p=p′)
+sol = solve(prob′, ImplicitEuler());
+@test sol.retcode == ReturnCode.Success
+@test sol[y][end] == 2x[end]
+
+prob′′ = remake(prob; p=[s.src.data => x])
+@test prob′′.p isa Tuple
+
+
+
+
 # ------------------------ Mixed Type Converted to float (default behavior)
 
 vars = @variables y(t)=1 dy(t)=0 ddy(t)=0

From a9c56b616b0a1bb57e1b8ba96fa81045eee03e4e Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Tue, 26 Sep 2023 10:10:19 -0400
Subject: [PATCH 06/50] Suppose heterogeneous parameters for linearize and
 remake

---
 Project.toml                             |  2 +-
 src/systems/abstractsystem.jl            | 18 ++++-
 src/systems/diffeqs/abstractodesystem.jl | 10 ++-
 src/systems/diffeqs/odesystem.jl         | 10 ++-
 src/utils.jl                             |  6 +-
 src/variables.jl                         | 26 ++++---
 test/split_parameters.jl                 | 90 +++++++++---------------
 7 files changed, 85 insertions(+), 77 deletions(-)

diff --git a/Project.toml b/Project.toml
index 58d73f2921..f5035f44cb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -85,7 +85,7 @@ PrecompileTools = "1"
 RecursiveArrayTools = "2.3"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SciMLBase = "2.0.1"
+SciMLBase = "1, 2.0.1"
 Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "0.1.0"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 48d647f621..df0ff4773f 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -230,7 +230,8 @@ for prop in [:eqs
     :metadata
     :gui_metadata
     :discrete_subsystems
-    :unknown_states]
+    :unknown_states
+    :split_idxs]
     fname1 = Symbol(:get_, prop)
     fname2 = Symbol(:has_, prop)
     @eval begin
@@ -1274,14 +1275,25 @@ See also [`linearize`](@ref) which provides a higher-level interface.
 function linearization_function(sys::AbstractSystem, inputs,
     outputs; simplify = false,
     initialize = true,
+    op = Dict(),
+    p = DiffEqBase.NullParameters(),
     kwargs...)
     sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify,
         kwargs...)
+    x0 = merge(defaults(sys), op)
+    u0, p, _ = get_u0_p(sys, x0, p; use_union = false, tofloat = true)
+    p, split_idxs = split_parameters_by_type(p)
+    ps = parameters(sys)
+    if p isa Tuple
+        ps = Base.Fix1(getindex, ps).(split_idxs)
+        ps = (ps...,) #if p is Tuple, ps should be Tuple
+    end
+
     lin_fun = let diff_idxs = diff_idxs,
         alge_idxs = alge_idxs,
         input_idxs = input_idxs,
         sts = states(sys),
-        fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys),
+        fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys, states(sys), ps; p = p),
         h = build_explicit_observed_function(sys, outputs),
         chunk = ForwardDiff.Chunk(input_idxs)
 
@@ -1600,11 +1612,11 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
     allow_input_derivatives = false,
     zero_dummy_der = false,
     kwargs...)
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; kwargs...)
     if zero_dummy_der
         dummyder = setdiff(states(ssys), states(sys))
         op = merge(op, Dict(x => 0.0 for x in dummyder))
     end
+    lin_fun, ssys = linearization_function(sys, inputs, outputs; op, kwargs...)
     linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys
 end
 
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 8528dd3c12..f11d1283c6 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -354,6 +354,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s
     checkbounds = false,
     sparsity = false,
     analytic = nothing,
+    split_idxs = nothing,
     kwargs...) where {iip, specialize}
     f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression},
         expression_module = eval_module, checkbounds = checkbounds,
@@ -508,6 +509,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s
         nothing
     end
 
+    @set! sys.split_idxs = split_idxs
     ODEFunction{iip, specialize}(f;
         sys = sys,
         jac = _jac === nothing ? nothing : _jac,
@@ -765,7 +767,7 @@ Take dictionaries with initial conditions and parameters and convert them to num
 """
 function get_u0_p(sys,
     u0map,
-    parammap;
+    parammap = nothing;
     use_union = true,
     tofloat = true,
     symbolic_u0 = false)
@@ -773,7 +775,9 @@ function get_u0_p(sys,
     ps = parameters(sys)
 
     defs = defaults(sys)
-    defs = mergedefaults(defs, parammap, ps)
+    if parammap !== nothing
+        defs = mergedefaults(defs, parammap, ps)
+    end
     defs = mergedefaults(defs, u0map, dvs)
 
     if symbolic_u0
@@ -835,7 +839,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap;
     f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac,
         checkbounds = checkbounds, p = p,
         linenumbers = linenumbers, parallel = parallel, simplify = simplify,
-        sparse = sparse, eval_expression = eval_expression,
+        sparse = sparse, eval_expression = eval_expression, split_idxs,
         kwargs...)
     implicit_dae ? (f, du0, u0, p) : (f, u0, p)
 end
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index dcbecfcfb0..3b828702c2 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -139,6 +139,10 @@ struct ODESystem <: AbstractODESystem
     used for ODAEProblem.
     """
     unknown_states::Union{Nothing, Vector{Any}}
+    """
+    split_idxs: a vector of vectors of indices for the split parameters.
+    """
+    split_idxs::Union{Nothing, Vector{Vector{Int}}}
 
     function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad,
         jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults,
@@ -146,8 +150,8 @@ struct ODESystem <: AbstractODESystem
         devents, metadata = nothing, gui_metadata = nothing,
         tearing_state = nothing,
         substitutions = nothing, complete = false,
-        discrete_subsystems = nothing, unknown_states = nothing;
-        checks::Union{Bool, Int} = true)
+        discrete_subsystems = nothing, unknown_states = nothing,
+        split_idxs = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
             check_parameters(ps, iv)
@@ -161,7 +165,7 @@ struct ODESystem <: AbstractODESystem
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,
             connector_type, preface, cevents, devents, metadata, gui_metadata,
             tearing_state, substitutions, complete, discrete_subsystems,
-            unknown_states)
+            unknown_states, split_idxs)
     end
 end
 
diff --git a/src/utils.jl b/src/utils.jl
index 0a6c5089a3..ab17bd3d8d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -661,7 +661,7 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
     else
         sym_vs = filter(x -> SymbolicUtils.issym(x) || SymbolicUtils.istree(x), vs)
         isempty(sym_vs) || throw_missingvars_in_sys(sym_vs)
-        
+
         C = nothing
         for v in vs
             E = typeof(v)
@@ -676,7 +676,7 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
             if use_union
                 C = Union{C, E}
             else
-                @assert C == E "`promote_to_concrete` can't make type $E uniform with $C"
+                @assert C==E "`promote_to_concrete` can't make type $E uniform with $C"
                 C = E
             end
         end
@@ -686,7 +686,7 @@ function promote_to_concrete(vs; tofloat = true, use_union = true)
             if (vs[i] isa Number) & tofloat
                 y[i] = float(vs[i]) #needed because copyto! can't convert Int to Float automatically
             else
-                y[i] = vs[i] 
+                y[i] = vs[i]
             end
         end
 
diff --git a/src/variables.jl b/src/variables.jl
index 4cffd3fdeb..854680998e 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -145,15 +145,23 @@ function SciMLBase.process_p_u0_symbolic(prob::Union{SciMLBase.AbstractDEProblem
                                 " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order."))
     end
 
-    # assemble defaults
-    defs = defaults(prob.f.sys)
-    defs = mergedefaults(defs, prob.p, parameters(prob.f.sys))
-    defs = mergedefaults(defs, p, parameters(prob.f.sys))
-    defs = mergedefaults(defs, prob.u0, states(prob.f.sys))
-    defs = mergedefaults(defs, u0, states(prob.f.sys))
-
-    u0 = varmap_to_vars(u0, states(prob.f.sys); defaults = defs, tofloat = true)
-    p = varmap_to_vars(p, parameters(prob.f.sys); defaults = defs)
+    sys = prob.f.sys
+    defs = defaults(sys)
+    ps = parameters(sys)
+    if has_split_idxs(sys) && (split_idxs = get_split_idxs(sys)) !== nothing
+        for (i, idxs) in enumerate(split_idxs)
+            defs = mergedefaults(defs, prob.p[i], ps[idxs])
+        end
+    else
+        # assemble defaults
+        defs = defaults(sys)
+        defs = mergedefaults(defs, prob.p, ps)
+    end
+    defs = mergedefaults(defs, p, ps)
+    sts = states(sys)
+    defs = mergedefaults(defs, prob.u0, sts)
+    defs = mergedefaults(defs, u0, sts)
+    u0, p, defs = get_u0_p(sys, defs)
 
     return p, u0
 end
diff --git a/test/split_parameters.jl b/test/split_parameters.jl
index 3f231e2dd7..d37b9e2c13 100644
--- a/test/split_parameters.jl
+++ b/test/split_parameters.jl
@@ -2,34 +2,29 @@ using ModelingToolkit, Test
 using ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq
 
-
-x = [1, 2.0, false, [1,2,3], Parameter(1.0)]
+x = [1, 2.0, false, [1, 2, 3], Parameter(1.0)]
 
 y = ModelingToolkit.promote_to_concrete(x)
 @test eltype(y) == Union{Float64, Parameter{Float64}, Vector{Int64}}
 
-y = ModelingToolkit.promote_to_concrete(x; tofloat=false)
+y = ModelingToolkit.promote_to_concrete(x; tofloat = false)
 @test eltype(y) == Union{Bool, Float64, Int64, Parameter{Float64}, Vector{Int64}}
 
-
-x = [1, 2.0, false, [1,2,3]]
+x = [1, 2.0, false, [1, 2, 3]]
 y = ModelingToolkit.promote_to_concrete(x)
 @test eltype(y) == Union{Float64, Vector{Int64}}
 
 x = Any[1, 2.0, false]
-y = ModelingToolkit.promote_to_concrete(x; tofloat=false)
+y = ModelingToolkit.promote_to_concrete(x; tofloat = false)
 @test eltype(y) == Union{Bool, Float64, Int64}
 
-y = ModelingToolkit.promote_to_concrete(x; use_union=false)
+y = ModelingToolkit.promote_to_concrete(x; use_union = false)
 @test eltype(y) == Float64
 
-x = Float16[1., 2., 3.]
+x = Float16[1.0, 2.0, 3.0]
 y = ModelingToolkit.promote_to_concrete(x)
 @test eltype(y) == Float16
 
-
-
-
 # ------------------------ Mixed Single Values and Vector
 
 dt = 4e-4
@@ -74,7 +69,7 @@ eqs = [y ~ src.output.u
 @named sys = ODESystem(eqs, t, vars, []; systems = [int, src])
 s = complete(sys)
 sys = structural_simplify(sys)
-prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat=false)
+prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat = false)
 @test prob.p isa Tuple{Vector{Float64}, Vector{Int}, Vector{Vector{Float64}}}
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
@@ -83,18 +78,15 @@ sol = solve(prob, ImplicitEuler());
 #TODO: remake becomes more complicated now, how to improve?
 defs = ModelingToolkit.defaults(sys)
 defs[s.src.data] = 2x
-p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false)
+p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false)
 p′, = ModelingToolkit.split_parameters_by_type(p′) #NOTE: we need to ensure this is called now before calling remake()
-prob′ = remake(prob; p=p′)
+prob′ = remake(prob; p = p′)
 sol = solve(prob′, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
 @test sol[y][end] == 2x[end]
 
-prob′′ = remake(prob; p=[s.src.data => x])
-@test prob′′.p isa Tuple
-
-
-
+prob′′ = remake(prob; p = [s.src.data => x])
+@test_broken prob′′.p isa Tuple
 
 # ------------------------ Mixed Type Converted to float (default behavior)
 
@@ -122,11 +114,6 @@ prob = ODEProblem(sys, [], tspan, []; tofloat = false)
 sol = solve(prob, ImplicitEuler());
 @test sol.retcode == ReturnCode.Success
 
-
-
-
-
-
 # ------------------------- Bug
 using ModelingToolkit, LinearAlgebra
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
@@ -136,51 +123,48 @@ using ModelingToolkit: connect
 
 "A wrapper function to make symbolic indexing easier"
 function wr(sys)
-    ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems=[sys], name=:a_wrapper)
+    ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems = [sys], name = :a_wrapper)
 end
-indexof(sym,syms) = findfirst(isequal(sym),syms)
+indexof(sym, syms) = findfirst(isequal(sym), syms)
 
 # Parameters
-m1 = 1.
-m2 = 1.
-k = 10. # Spring stiffness
-c = 3.  # Damping coefficient
+m1 = 1.0
+m2 = 1.0
+k = 10.0 # Spring stiffness
+c = 3.0  # Damping coefficient
 
 @named inertia1 = Inertia(; J = m1)
 @named inertia2 = Inertia(; J = m2)
 @named spring = Spring(; c = k)
 @named damper = Damper(; d = c)
-@named torque = Torque(use_support=false)
+@named torque = Torque(use_support = false)
 
-function SystemModel(u=nothing; name=:model)
-    eqs = [
-        connect(torque.flange, inertia1.flange_a)
+function SystemModel(u = nothing; name = :model)
+    eqs = [connect(torque.flange, inertia1.flange_a)
         connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
-        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
-    ]
+        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
     if u !== nothing
         push!(eqs, connect(torque.tau, u.output))
-        return @named model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper, u])
+        return @named model = ODESystem(eqs,
+            t;
+            systems = [torque, inertia1, inertia2, spring, damper, u])
     end
     ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
 end
 
-
 model = SystemModel() # Model with load disturbance
-@named d = Step(start_time=1., duration=10., offset=0., height=1.) # Disturbance
+@named d = Step(start_time = 1.0, duration = 10.0, offset = 0.0, height = 1.0) # Disturbance
 model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi] # This is the state realization we want to control
 inputs = [model.torque.tau.u]
 matrices, ssys = ModelingToolkit.linearize(wr(model), inputs, model_outputs)
 
 # Design state-feedback gain using LQR
 # Define cost matrices
-x_costs = [
-    model.inertia1.w =>   1.
-    model.inertia2.w =>   1.
-    model.inertia1.phi => 1.
-    model.inertia2.phi => 1.
-]
-L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller
+x_costs = [model.inertia1.w => 1.0
+    model.inertia2.w => 1.0
+    model.inertia1.phi => 1.0
+    model.inertia2.phi => 1.0]
+L = randn(1, 4) # Post-multiply by `C` to get the correct input to the controller
 
 # This old definition of MatrixGain will work because the parameter space does not include K (an Array term)
 # @component function MatrixGainAlt(K::AbstractArray; name)
@@ -191,16 +175,12 @@ L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller
 #     compose(ODESystem(eqs, t, [], []; name = name), [input, output])
 # end
 
-@named state_feedback = MatrixGain(K=-L) # Build negative feedback into the feedback matrix
-@named add = Add(;k1=1., k2=1.) # To add the control signal and the disturbance
+@named state_feedback = MatrixGain(K = -L) # Build negative feedback into the feedback matrix
+@named add = Add(; k1 = 1.0, k2 = 1.0) # To add the control signal and the disturbance
 
-connections = [
-    [state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4]
+connections = [[state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4]
     connect(d.output, :d, add.input1)
     connect(add.input2, state_feedback.output)
-    connect(add.output, :u, model.torque.tau)
-]
-closed_loop = ODESystem(connections, t, systems=[model, state_feedback, add, d], name=:closed_loop)
+    connect(add.output, :u, model.torque.tau)]
+@named closed_loop = ODESystem(connections, t, systems = [model, state_feedback, add, d])
 S = get_sensitivity(closed_loop, :u)
-
-

From 1a96915f65de389430566c990e8aa308a4ae523e Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Tue, 26 Sep 2023 10:42:49 -0400
Subject: [PATCH 07/50] Fix zero_dummy_der

Co-authored-by: Fredrik Bagge Carlson <baggepinnen@gmail.com>
---
 src/systems/abstractsystem.jl | 22 ++++++++++++++++------
 test/input_output_handling.jl |  4 ++--
 2 files changed, 18 insertions(+), 8 deletions(-)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index df0ff4773f..c452417409 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1277,9 +1277,18 @@ function linearization_function(sys::AbstractSystem, inputs,
     initialize = true,
     op = Dict(),
     p = DiffEqBase.NullParameters(),
+    zero_dummy_der = false,
     kwargs...)
-    sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify,
+    ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs;
+        simplify,
         kwargs...)
+    if zero_dummy_der
+        dummyder = setdiff(states(ssys), states(sys))
+        defs = Dict(x => 0.0 for x in dummyder)
+        @set! ssys.defaults = merge(defs, defaults(ssys))
+        op = merge(defs, op)
+    end
+    sys = ssys
     x0 = merge(defaults(sys), op)
     u0, p, _ = get_u0_p(sys, x0, p; use_union = false, tofloat = true)
     p, split_idxs = split_parameters_by_type(p)
@@ -1612,11 +1621,12 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
     allow_input_derivatives = false,
     zero_dummy_der = false,
     kwargs...)
-    if zero_dummy_der
-        dummyder = setdiff(states(ssys), states(sys))
-        op = merge(op, Dict(x => 0.0 for x in dummyder))
-    end
-    lin_fun, ssys = linearization_function(sys, inputs, outputs; op, kwargs...)
+    lin_fun, ssys = linearization_function(sys,
+        inputs,
+        outputs;
+        zero_dummy_der,
+        op,
+        kwargs...)
     linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys
 end
 
diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl
index 12f68f9f04..249a94e9d0 100644
--- a/test/input_output_handling.jl
+++ b/test/input_output_handling.jl
@@ -124,8 +124,8 @@ using ModelingToolkitStandardLibrary.Mechanical.Rotational
 t = ModelingToolkitStandardLibrary.Mechanical.Rotational.t
 @named inertia1 = Inertia(; J = 1)
 @named inertia2 = Inertia(; J = 1)
-@named spring = Spring(; c = 10)
-@named damper = Damper(; d = 3)
+@named spring = Rotational.Spring(; c = 10)
+@named damper = Rotational.Damper(; d = 3)
 @named torque = Torque(; use_support = false)
 @variables y(t) = 0
 eqs = [connect(torque.flange, inertia1.flange_a)

From 023e024f3ebd82991f34ea3900b9689572ac80ae Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Tue, 26 Sep 2023 14:25:27 -0400
Subject: [PATCH 08/50] Relax tests

---
 test/odesystem.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/odesystem.jl b/test/odesystem.jl
index 56d436d76c..336eca11e7 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -241,7 +241,7 @@ p2 = (k₁ => 0.04,
     k₃ => 1e4)
 tspan = (0.0, 100000.0)
 prob1 = ODEProblem(sys, u0, tspan, p)
-@test prob1.f.sys === sys
+@test prob1.f.sys == sys
 prob12 = ODEProblem(sys, u0, tspan, [0.04, 3e7, 1e4])
 prob13 = ODEProblem(sys, u0, tspan, (0.04, 3e7, 1e4))
 prob14 = ODEProblem(sys, u0, tspan, p2)

From 119cb27546c31ed22309240c72d14a952926f7b6 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Tue, 26 Sep 2023 16:50:06 -0400
Subject: [PATCH 09/50] Don't use SciMLBase 2

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

diff --git a/Project.toml b/Project.toml
index f5035f44cb..bafa1835a3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -85,7 +85,7 @@ PrecompileTools = "1"
 RecursiveArrayTools = "2.3"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SciMLBase = "1, 2.0.1"
+SciMLBase = "1"
 Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "0.1.0"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"

From 29aec7b3d19f6f333c863956dba07d02b070f5bf Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Mon, 2 Oct 2023 09:38:17 -0400
Subject: [PATCH 10/50] Fix eq ordering

---
 test/stream_connectors.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl
index 9502b55be9..45d58c8552 100644
--- a/test/stream_connectors.jl
+++ b/test/stream_connectors.jl
@@ -157,8 +157,8 @@ eqns = [domain_connect(fluid, n1m1.port_a)
     pipe.port_b.P ~ sink.port.P]
 @test sort(equations(expand_connections(n1m1Test)), by = string) ==
       [0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
-    0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ n1m1.source.port1.m_flow - n1m1.port_a.m_flow
+    0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ pipe.port_b.m_flow + sink.port.m_flow
     fluid.m_flow ~ 0
     n1m1.port_a.P ~ pipe.port_a.P

From addf3357af1f173a876fa735c473dd2b3ce5cba1 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Mon, 2 Oct 2023 09:53:16 -0400
Subject: [PATCH 11/50] Change bounds back

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

diff --git a/Project.toml b/Project.toml
index a52124a694..470adc92b6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -85,7 +85,7 @@ PrecompileTools = "1"
 RecursiveArrayTools = "2.3"
 Reexport = "0.2, 1"
 RuntimeGeneratedFunctions = "0.5.9"
-SciMLBase = "1"
+SciMLBase = "2.0.1"
 Setfield = "0.7, 0.8, 1"
 SimpleNonlinearSolve = "0.1.0"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"

From 2e4c1bbd6f31249a8f7045c68cde346885b04025 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Mon, 2 Oct 2023 09:55:29 -0400
Subject: [PATCH 12/50] Format

---
 docs/src/basics/MTKModel_Connector.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/basics/MTKModel_Connector.md b/docs/src/basics/MTKModel_Connector.md
index 8cc106b05d..ac285d0253 100644
--- a/docs/src/basics/MTKModel_Connector.md
+++ b/docs/src/basics/MTKModel_Connector.md
@@ -130,6 +130,7 @@ end
 `@connector`s accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`.
 
 !!! note
+    
     For more examples of usage, checkout [ModelingToolkitStandardLibrary.jl](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/)
 
 ### What's a `structure` dictionary?

From 2a553a7479da55a41d44a89b98633a2433351710 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Mon, 2 Oct 2023 14:08:07 -0400
Subject: [PATCH 13/50] Be careful with string sorting

---
 test/odesystem.jl         |  7 +++---
 test/stream_connectors.jl | 45 ++++++++++++++++++++-------------------
 2 files changed, 27 insertions(+), 25 deletions(-)

diff --git a/test/odesystem.jl b/test/odesystem.jl
index 336eca11e7..39bcda58de 100644
--- a/test/odesystem.jl
+++ b/test/odesystem.jl
@@ -23,7 +23,8 @@ eqs = [D(x) ~ σ * (y - x),
 ModelingToolkit.toexpr.(eqs)[1]
 @named de = ODESystem(eqs; defaults = Dict(x => 1))
 subed = substitute(de, [σ => k])
-@test isequal(sort(parameters(subed), by = string), [k, β, ρ])
+ssort(eqs) = sort(eqs, by = string)
+@test isequal(ssort(parameters(subed)), [k, β, ρ])
 @test isequal(equations(subed),
     [D(x) ~ k * (y - x)
         D(y) ~ (ρ - z) * x - y
@@ -348,14 +349,14 @@ eqs = [0 ~ x + z
     D(accumulation_y) ~ y
     D(accumulation_z) ~ z
     D(x) ~ y]
-@test sort(equations(asys), by = string) == eqs
+@test ssort(equations(asys)) == ssort(eqs)
 @variables ac(t)
 asys = add_accumulations(sys, [ac => (x + y)^2])
 eqs = [0 ~ x + z
     0 ~ x - y
     D(ac) ~ (x + y)^2
     D(x) ~ y]
-@test sort(equations(asys), by = string) == eqs
+@test ssort(equations(asys)) == ssort(eqs)
 
 sys2 = ode_order_lowering(sys)
 M = ModelingToolkit.calculate_massmatrix(sys2)
diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl
index 45d58c8552..e4365876cf 100644
--- a/test/stream_connectors.jl
+++ b/test/stream_connectors.jl
@@ -136,27 +136,28 @@ eqns = [domain_connect(fluid, n1m1.port_a)
 
 @test_nowarn structural_simplify(n1m1Test)
 @unpack source, port_a = n1m1
-@test sort(equations(expand_connections(n1m1)), by = string) == [0 ~ port_a.m_flow
+ssort(eqs) = sort(eqs, by = string)
+@test ssort(equations(expand_connections(n1m1))) == ssort([0 ~ port_a.m_flow
     0 ~ source.port1.m_flow - port_a.m_flow
     source.port1.P ~ port_a.P
     source.port1.P ~ source.P
     source.port1.h_outflow ~ port_a.h_outflow
-    source.port1.h_outflow ~ source.h]
+    source.port1.h_outflow ~ source.h])
 @unpack port_a, port_b = pipe
-@test sort(equations(expand_connections(pipe)), by = string) ==
-      [0 ~ -port_a.m_flow - port_b.m_flow
+@test ssort(equations(expand_connections(pipe))) ==
+      ssort([0 ~ -port_a.m_flow - port_b.m_flow
     0 ~ port_a.m_flow
     0 ~ port_b.m_flow
     port_a.P ~ port_b.P
     port_a.h_outflow ~ instream(port_b.h_outflow)
-    port_b.h_outflow ~ instream(port_a.h_outflow)]
-@test sort(equations(expand_connections(sys)), by = string) ==
-      [0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
+    port_b.h_outflow ~ instream(port_a.h_outflow)])
+@test ssort(equations(expand_connections(sys))) ==
+      ssort([0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ pipe.port_b.m_flow + sink.port.m_flow
     n1m1.port_a.P ~ pipe.port_a.P
-    pipe.port_b.P ~ sink.port.P]
-@test sort(equations(expand_connections(n1m1Test)), by = string) ==
-      [0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
+    pipe.port_b.P ~ sink.port.P])
+@test ssort(equations(expand_connections(n1m1Test))) ==
+      ssort([0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow
     0 ~ n1m1.source.port1.m_flow - n1m1.port_a.m_flow
     0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow
     0 ~ pipe.port_b.m_flow + sink.port.m_flow
@@ -172,7 +173,7 @@ eqns = [domain_connect(fluid, n1m1.port_a)
     pipe.port_b.h_outflow ~ n1m1.port_a.h_outflow
     sink.port.P ~ sink.P
     sink.port.h_outflow ~ sink.h_in
-    sink.port.m_flow ~ -sink.m_flow_in]
+    sink.port.m_flow ~ -sink.m_flow_in])
 
 # N1M2 model and test code.
 function N1M2(; name,
@@ -255,12 +256,12 @@ eqns = [connect(source.port, n2m2.port_a)
 @named sp2 = TwoPhaseFluidPort()
 @named sys = ODESystem([connect(sp1, sp2)], t)
 sys_exp = expand_connections(compose(sys, [sp1, sp2]))
-@test sort(equations(sys_exp), by = string) == [0 ~ -sp1.m_flow - sp2.m_flow
+@test ssort(equations(sys_exp)) == ssort([0 ~ -sp1.m_flow - sp2.m_flow
     0 ~ sp1.m_flow
     0 ~ sp2.m_flow
     sp1.P ~ sp2.P
     sp1.h_outflow ~ ModelingToolkit.instream(sp2.h_outflow)
-    sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)]
+    sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)])
 
 # array var
 @connector function VecPin(; name)
@@ -274,15 +275,15 @@ end
 
 @named simple = ODESystem([connect(vp1, vp2, vp3)], t)
 sys = expand_connections(compose(simple, [vp1, vp2, vp3]))
-@test sort(equations(sys), by = string) == sort([0 .~ collect(vp1.i)
-        0 .~ collect(vp2.i)
-        0 .~ collect(vp3.i)
-        vp1.v[1] ~ vp2.v[1]
-        vp1.v[2] ~ vp2.v[2]
-        vp1.v[1] ~ vp3.v[1]
-        vp1.v[2] ~ vp3.v[2]
-        0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1]
-        0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]], by = string)
+@test ssort(equations(sys)) == ssort([0 .~ collect(vp1.i)
+    0 .~ collect(vp2.i)
+    0 .~ collect(vp3.i)
+    vp1.v[1] ~ vp2.v[1]
+    vp1.v[2] ~ vp2.v[2]
+    vp1.v[1] ~ vp3.v[1]
+    vp1.v[2] ~ vp3.v[2]
+    0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1]
+    0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]])
 
 @connector function VectorHeatPort(; name, N = 100, T0 = 0.0, Q0 = 0.0)
     @variables (T(t))[1:N]=T0 (Q(t))[1:N]=Q0 [connect = Flow]

From f9ea3a6a294cd7ddcd2ef3fda871f7109a6d9626 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Mon, 2 Oct 2023 15:21:45 -0400
Subject: [PATCH 14/50] Update LaTeXify ref tests

---
 test/latexify/10.tex | 6 +++---
 test/latexify/20.tex | 4 ++--
 test/latexify/30.tex | 4 ++--
 3 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/test/latexify/10.tex b/test/latexify/10.tex
index 9a8335bd34..09ec74ab72 100644
--- a/test/latexify/10.tex
+++ b/test/latexify/10.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\sigma \left(  - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left(  - y\left( t \right) + x\left( t \right) \right)}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\
-0 =&  - y\left( t \right) + \frac{1}{10} \sigma \left( \rho - z\left( t \right) \right) x\left( t \right) \\
-\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - \beta z\left( t \right)
+\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\left(  - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left( x\left( t \right) - y\left( t \right) \right) \sigma}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\
+0 =&  - y\left( t \right) + \frac{1}{10} x\left( t \right) \left(  - z\left( t \right) + \rho \right) \sigma \\
+\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - z\left( t \right) \beta
 \end{align}
diff --git a/test/latexify/20.tex b/test/latexify/20.tex
index e32ea81eb7..9de213aafd 100644
--- a/test/latexify/20.tex
+++ b/test/latexify/20.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
-0 =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left(  - u(t)_1 + u(t)_2 \right) \\
+0 =&  - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\
 \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}
diff --git a/test/latexify/30.tex b/test/latexify/30.tex
index bc46f41e9b..b83feeba72 100644
--- a/test/latexify/30.tex
+++ b/test/latexify/30.tex
@@ -1,5 +1,5 @@
 \begin{align}
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
-\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left(  - u(t)_1 + u(t)_2 \right) \\
+\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =&  - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\
 \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}

From c2f92f127c94c5369fc61178106cb1d3bded838b Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 4 Oct 2023 11:57:52 -0400
Subject: [PATCH 15/50] init

---
 Project.toml                |  2 ++
 ext/MTKBifurcationKitExt.jl | 44 +++++++++++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+)
 create mode 100644 ext/MTKBifurcationKitExt.jl

diff --git a/Project.toml b/Project.toml
index 470adc92b6..ae899f1220 100644
--- a/Project.toml
+++ b/Project.toml
@@ -51,9 +51,11 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 
 [extensions]
+MTKBifurcationKitExt = "BifurcationKit"
 MTKDeepDiffsExt = "DeepDiffs"
 
 [compat]
diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
new file mode 100644
index 0000000000..5aa37e7117
--- /dev/null
+++ b/ext/MTKBifurcationKitExt.jl
@@ -0,0 +1,44 @@
+module MTKBifurcationKitExt
+
+println("BifurcationKit extension loaded")
+
+### Preparations ###
+
+# Imports
+using ModelingToolkit, Setfield
+import BifurcationKit: BifurcationProblem
+
+### Creates BifurcationProblem Overloads ###
+
+# When input is a NonlinearProblem.
+function BifurcationKit.BifurcationProblem(nprob::NonlinearProblem, u0, p, bif_par, args...; plot_variable=nothing, record_from_solution=record_sol_default, kwargs...)
+    # check if jac does nto exist and modify.
+    F = nprob.f.f
+    J = nprob.f.jac.f
+
+    bif_idx = findfirst(isequal(bif_par), parameters(nsys))
+    if isnothing(plot_variable) 
+        plot_idx = findfirst(isequal(plot_variable), variables(nsys))
+        record_from_solution = (x, p) -> x[plot_idx]
+    end
+    return BifurcationProblem(F, u0, p, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
+end
+
+# When input is a NonlinearSystem.
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0, p, args...; kwargs...)
+    return BifurcationProblem(NonlinearProblem(nsys, u0, p; jac=true), args...; kwargs...)
+end
+
+# When input is a ODEProblem.
+function BifurcationKit.BifurcationProblem(oprob::ODEProblem, u0, p, args...; kwargs...)
+    return BifurcationProblem(convert(NonlinearProblem, oprob), args...; kwargs...)
+end
+
+# When input is a ODESystem.
+function BifurcationKit.BifurcationProblem(osys::ODESystem, u0, p, args...; kwargs...)
+    return BifurcationProblem(convert(NonlinearProblem, osys), args...; kwargs...)
+end
+
+
+
+end # module

From b05bc5f694a1cb20f77a80fc00631185d0dc222f Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 4 Oct 2023 12:44:11 -0400
Subject: [PATCH 16/50] makefunction

---
 ext/MTKBifurcationKitExt.jl | 30 +++++++++++++-----------------
 1 file changed, 13 insertions(+), 17 deletions(-)

diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index 5aa37e7117..60ee91b9df 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -6,32 +6,28 @@ println("BifurcationKit extension loaded")
 
 # Imports
 using ModelingToolkit, Setfield
-import BifurcationKit: BifurcationProblem
+import BifurcationKit
 
 ### Creates BifurcationProblem Overloads ###
 
-# When input is a NonlinearProblem.
-function BifurcationKit.BifurcationProblem(nprob::NonlinearProblem, u0, p, bif_par, args...; plot_variable=nothing, record_from_solution=record_sol_default, kwargs...)
-    # check if jac does nto exist and modify.
-    F = nprob.f.f
-    J = nprob.f.jac.f
+# When input is a NonlinearSystem.
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_start, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, kwargs...)
+    # Creates F and J functions.
+    ofun = NonlinearFunction(nsys; jac=true)
+    F = ofun.f
+    J = ofun.jac
 
+    # Computes bifurcation parameter and plot var indexes.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))
-    if isnothing(plot_variable) 
-        plot_idx = findfirst(isequal(plot_variable), variables(nsys))
+    if !isnothing(plot_var) 
+        plot_idx = findfirst(isequal(plot_var), states(nsys))
         record_from_solution = (x, p) -> x[plot_idx]
     end
-    return BifurcationProblem(F, u0, p, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
-end
 
-# When input is a NonlinearSystem.
-function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0, p, args...; kwargs...)
-    return BifurcationProblem(NonlinearProblem(nsys, u0, p; jac=true), args...; kwargs...)
-end
+    # Converts the input state guess.
+    u0_guess = ModelingToolkit.varmap_to_vars(u0_guess, states(nsys))
 
-# When input is a ODEProblem.
-function BifurcationKit.BifurcationProblem(oprob::ODEProblem, u0, p, args...; kwargs...)
-    return BifurcationProblem(convert(NonlinearProblem, oprob), args...; kwargs...)
+    return BifurcationKit.BifurcationProblem(F, u0_guess, [p_start], (@lens _[1]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
 end
 
 # When input is a ODESystem.

From f2d80b025965983d7d7f5fd65b6ead68ba245ee0 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Wed, 4 Oct 2023 15:27:49 -0400
Subject: [PATCH 17/50] 1.6 has a different term ordering

---
 test/runtests.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index e6387e2cd9..b60d82a6ce 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -56,4 +56,6 @@ end
 @safetestset "FuncAffect Test" include("funcaffect.jl")
 @safetestset "Constants Test" include("constants.jl")
 # Reference tests go Last
-@safetestset "Latexify recipes Test" include("latexify.jl")
+if VERSION >= v"1.9"
+    @safetestset "Latexify recipes Test" include("latexify.jl")
+end

From ae050413d35737edc0e7a404708c4159e978aac4 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Wed, 4 Oct 2023 17:00:00 -0400
Subject: [PATCH 18/50] Fix DiscreteProblem construction

---
 src/systems/jumps/jumpsystem.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl
index 8e87d145c2..b37aad5d19 100644
--- a/src/systems/jumps/jumpsystem.jl
+++ b/src/systems/jumps/jumpsystem.jl
@@ -293,7 +293,7 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothing},
     parammap = DiffEqBase.NullParameters();
     checkbounds = false,
-    use_union = false,
+    use_union = true,
     kwargs...)
     dvs = states(sys)
     ps = parameters(sys)

From fb7c3af551c64d48e338d398892086c6dc90b383 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Wed, 4 Oct 2023 18:02:32 -0400
Subject: [PATCH 19/50] Update Project.toml

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

diff --git a/Project.toml b/Project.toml
index 470adc92b6..3834488ae3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "8.71.1"
+version = "8.71.2"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

From 6e795b8574c4de02681b7a771c1ab3e3f1185926 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Fri, 6 Oct 2023 14:41:01 +0200
Subject: [PATCH 20/50] Fix recursive structure and setup precompilation

The recursive structure of SystemStructure was giving precompilation an issue so I removed that module. Then I setup precompilation on the basic interface.  Test case:

```julia
@time using ModelingToolkit
@time using DifferentialEquations
@time using Plots
@time begin
    @variables t x(t)
    D = Differential(t)
    @named sys = ODESystem([D(x) ~ -x])
end;
@time prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], jac = true);
@time sol = solve(prob);
@time plot(sol, idxs=[x]);
```

Before:

```
11.082586 seconds (19.32 M allocations: 1.098 GiB, 4.09% gc time, 0.71% compilation time: 87% of which was recompilation)
0.639738 seconds (661.39 k allocations: 101.321 MiB, 4.33% gc time, 6.46% compilation time)
3.703724 seconds (5.71 M allocations: 322.840 MiB, 5.22% gc time, 9.92% compilation time: 86% of which was recompilation)
7.795297 seconds (8.25 M allocations: 483.041 MiB, 2.50% gc time, 99.88% compilation time)
21.719376 seconds (44.11 M allocations: 2.485 GiB, 5.68% gc time, 99.48% compilation time)
2.602250 seconds (4.04 M allocations: 253.058 MiB, 4.60% gc time, 99.90% compilation time)
2.450509 seconds (5.17 M allocations: 332.101 MiB, 5.89% gc time, 99.41% compilation time: 30% of which was recompilation)
```

After:

```
9.129141 seconds (22.77 M allocations: 1.291 GiB, 4.65% gc time, 0.62% compilation time: 87% of which was recompilation)
0.784464 seconds (667.59 k allocations: 101.524 MiB, 3.95% gc time, 4.16% compilation time)
3.111142 seconds (5.42 M allocations: 305.594 MiB, 3.82% gc time, 6.39% compilation time: 82% of which was recompilation)
0.105567 seconds (157.39 k allocations: 10.522 MiB, 8.81% gc time, 95.49% compilation time: 74% of which was recompilation)
1.993642 seconds (4.03 M allocations: 218.310 MiB, 2.69% gc time, 96.95% compilation time: 82% of which was recompilation)
1.806758 seconds (4.06 M allocations: 254.371 MiB, 4.44% gc time, 99.91% compilation time)
1.694666 seconds (5.27 M allocations: 339.088 MiB, 6.18% gc time, 99.39% compilation time: 31% of which was recompilation)
```

And that's on v1.9, so v1.10 should be even better. 20 seconds off hehe.
---
 ext/MTKDeepDiffsExt.jl                           |  2 +-
 src/ModelingToolkit.jl                           | 16 ++++++++++++----
 .../StructuralTransformations.jl                 | 10 +++++++---
 src/systems/systemstructure.jl                   | 11 +++--------
 test/clock.jl                                    |  5 ++---
 5 files changed, 25 insertions(+), 19 deletions(-)

diff --git a/ext/MTKDeepDiffsExt.jl b/ext/MTKDeepDiffsExt.jl
index 92bc9ba2b9..1d361f96a3 100644
--- a/ext/MTKDeepDiffsExt.jl
+++ b/ext/MTKDeepDiffsExt.jl
@@ -4,7 +4,7 @@ using DeepDiffs, ModelingToolkit
 using ModelingToolkit.BipartiteGraphs: Label,
     BipartiteAdjacencyList, unassigned,
     HighlightInt
-using ModelingToolkit.SystemStructures: SystemStructure,
+using ModelingToolkit: SystemStructure,
     MatchedSystemStructure,
     SystemStructurePrintMatrix
 
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index d5c4172340..1a6e0356a5 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -51,8 +51,8 @@ using PrecompileTools, Reexport
     using MLStyle
 
     using Reexport
+    using Symbolics
     using Symbolics: degree
-    @reexport using Symbolics
     using Symbolics: _parse_vars, value, @derivatives, get_variables,
         exprs_occur_in, solve_for, build_expr, unwrap, wrap,
         VariableSource, getname, variable, Connection, connect,
@@ -70,9 +70,10 @@ using PrecompileTools, Reexport
     import OrdinaryDiffEq
 
     import Graphs: SimpleDiGraph, add_edge!, incidence_matrix
-
-    @reexport using UnPack
 end
+
+@reexport using Symbolics
+@reexport using UnPack
 RuntimeGeneratedFunctions.init(@__MODULE__)
 
 export @derivatives
@@ -156,7 +157,6 @@ include("systems/dependency_graphs.jl")
 include("clock.jl")
 include("discretedomain.jl")
 include("systems/systemstructure.jl")
-using .SystemStructures
 include("systems/clock_inference.jl")
 include("systems/systems.jl")
 
@@ -172,6 +172,14 @@ for S in subtypes(ModelingToolkit.AbstractSystem)
     @eval convert_system(::Type{<:$S}, sys::$S) = sys
 end
 
+PrecompileTools.@compile_workload begin
+    using ModelingToolkit
+    @variables t x(t)
+    D = Differential(t)
+    @named sys = ODESystem([D(x) ~ -x])
+    prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], jac = true)
+end
+
 export AbstractTimeDependentSystem,
     AbstractTimeIndependentSystem,
     AbstractMultivariateSystem
diff --git a/src/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl
index 47b140b7d1..7289df4232 100644
--- a/src/structural_transformation/StructuralTransformations.jl
+++ b/src/structural_transformation/StructuralTransformations.jl
@@ -23,14 +23,18 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
     IncrementalCycleTracker, add_edge_checked!, topological_sort,
     invalidate_cache!, Substitutions, get_or_construct_tearing_state,
     filter_kwargs, lower_varname, setio, SparseMatrixCLIL,
-    fast_substitute, get_fullvars, has_equations
+    fast_substitute, get_fullvars, has_equations, observed
 
 using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview, complete
 import ModelingToolkit: var_derivative!, var_derivative_graph!
 using Graphs
-using ModelingToolkit.SystemStructures
-using ModelingToolkit.SystemStructures: algeqs, EquationsView
+using ModelingToolkit: algeqs, EquationsView,
+    SystemStructure, TransformationState, TearingState, structural_simplify!,
+    isdiffvar, isdervar, isalgvar, isdiffeq, algeqs, is_only_discrete,
+    dervars_range, diffvars_range, algvars_range,
+    DiffGraph, complete!,
+    get_fullvars, system_subset
 
 using ModelingToolkit.DiffEqBase
 using ModelingToolkit.StaticArrays
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 4e330619d5..18e50afc9b 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -1,7 +1,5 @@
-module SystemStructures
-
 using DataStructures
-using Symbolics: linear_expansion, unwrap
+using Symbolics: linear_expansion, unwrap, Connection
 using SymbolicUtils: istree, operation, arguments, Symbolic
 using SymbolicUtils: quick_cancel, similarterm
 using ..ModelingToolkit
@@ -10,7 +8,7 @@ import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten,
     isparameter, isconstant,
     independent_variables, SparseMatrixCLIL, AbstractSystem,
     equations, isirreducible, input_timedomain, TimeDomain,
-    VariableType, getvariabletype, has_equations
+    VariableType, getvariabletype, has_equations, ODESystem
 using ..BipartiteGraphs
 import ..BipartiteGraphs: invview, complete
 using Graphs
@@ -27,8 +25,7 @@ function quick_cancel_expr(expr)
 end
 
 export SystemStructure, TransformationState, TearingState, structural_simplify!
-export initialize_system_structure, find_linear_equations
-export isdiffvar, isdervar, isalgvar, isdiffeq, isalgeq, algeqs, is_only_discrete
+export isdiffvar, isdervar, isalgvar, isdiffeq, algeqs, is_only_discrete
 export dervars_range, diffvars_range, algvars_range
 export DiffGraph, complete!
 export get_fullvars, system_subset
@@ -620,5 +617,3 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullstates)
     ModelingToolkit.invalidate_cache!(sys), input_idxs
 end
-
-end # module
diff --git a/test/clock.jl b/test/clock.jl
index 42074b65ff..74946d21d7 100644
--- a/test/clock.jl
+++ b/test/clock.jl
@@ -62,13 +62,12 @@ By inference:
 =>   Shift(x, 0, dt) := (Shift(x, -1, dt) + dt) / (1 - dt) # Discrete system
 =#
 
-using ModelingToolkit.SystemStructures
 ci, varmap = infer_clocks(sys)
 eqmap = ci.eq_domain
 tss, inputs = ModelingToolkit.split_system(deepcopy(ci))
-sss, = SystemStructures._structural_simplify!(deepcopy(tss[1]), (inputs[1], ()))
+sss, = ModelingToolkit._structural_simplify!(deepcopy(tss[1]), (inputs[1], ()))
 @test equations(sss) == [D(x) ~ u - x]
-sss, = SystemStructures._structural_simplify!(deepcopy(tss[2]), (inputs[2], ()))
+sss, = ModelingToolkit._structural_simplify!(deepcopy(tss[2]), (inputs[2], ()))
 @test isempty(equations(sss))
 @test observed(sss) == [yd ~ Sample(t, dt)(y); r ~ 1.0; ud ~ kp * (r - yd)]
 

From c5bd8210b20eada9994f851f83ea702dc6d85ac2 Mon Sep 17 00:00:00 2001
From: Venkateshprasad <32921645+ven-k@users.noreply.github.com>
Date: Wed, 27 Sep 2023 14:23:03 +0530
Subject: [PATCH 21/50] docs(refactor): update acausal components tutorial with
 `@mtkmodel`

---
 docs/src/tutorials/acausal_components.md | 301 ++++++++++++-----------
 1 file changed, 155 insertions(+), 146 deletions(-)

diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index cb2031d91a..b920c07df3 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -10,7 +10,7 @@ component variables. We then simplify this to an ODE by eliminating the
 equalities before solving. Let's see this in action.
 
 !!! note
-    
+
     This tutorial teaches how to build the entire RC circuit from scratch.
     However, to simulate electric components with more ease, check out the
     [ModelingToolkitStandardLibrary.jl](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/)
@@ -23,77 +23,88 @@ equalities before solving. Let's see this in action.
 using ModelingToolkit, Plots, DifferentialEquations
 
 @variables t
-@connector function Pin(; name)
-    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
-    ODESystem(Equation[], t, sts, []; name = name)
+@connector Pin begin
+    v(t) = 1.0
+    i(t) = 1.0, [connect = Flow]
 end
 
-@component function Ground(; name)
-    @named g = Pin()
-    eqs = [g.v ~ 0]
-    compose(ODESystem(eqs, t, [], []; name = name), g)
+@mtkmodel Ground begin
+    @components begin
+        g = Pin()
+    end
+    @equations begin
+        g.v ~ 0
+    end
 end
 
-@component function OnePort(; name)
-    @named p = Pin()
-    @named n = Pin()
-    sts = @variables v(t)=1.0 i(t)=1.0
-    eqs = [v ~ p.v - n.v
+@mtkmodel OnePort begin
+    @components begin
+        p = Pin()
+        n = Pin()
+    end
+    @variables begin
+        v(t) = 1.0
+        i(t) = 1.0
+    end
+    @equations begin
+        v ~ p.v - n.v
         0 ~ p.i + n.i
-        i ~ p.i]
-    compose(ODESystem(eqs, t, sts, []; name = name), p, n)
+        i ~ p.i
+    end
+end
+
+@mtkmodel Resistor begin
+    @extend v, i = oneport = OnePort()
+    @parameters begin
+        R = 1.0
+    end
+    @equations begin
+        v ~ i * R
+    end
 end
 
-@component function Resistor(; name, R = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters R = R
-    eqs = [
-        v ~ i * R,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+D = Differential(t)
+
+@mtkmodel Capacitor begin
+    @extend v, i = oneport = OnePort()
+    @parameters begin
+        C = 1.0
+    end
+    @equations begin
+        D(v) ~ i / C
+    end
 end
 
-@component function Capacitor(; name, C = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters C = C
-    D = Differential(t)
-    eqs = [
-        D(v) ~ i / C,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel ConstantVoltage begin
+    @extend (v,) = oneport = OnePort()
+    @parameters begin
+        V = 1.0
+    end
+    @equations begin
+        V ~ v
+    end
 end
 
-@component function ConstantVoltage(; name, V = 1.0)
-    @named oneport = OnePort()
-    @unpack v = oneport
-    ps = @parameters V = V
-    eqs = [
-        V ~ v,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel RCModel begin
+    @components begin
+        resistor = Resistor(R = 1.0)
+        capacitor = Capacitor(C = 1.0)
+        source = ConstantVoltage(V = 1.0)
+        ground = Ground()
+    end
+    @equations begin
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n)
+        connect(capacitor.n, ground.g)
+    end
 end
 
-R = 1.0
-C = 1.0
-V = 1.0
-@named resistor = Resistor(R = R)
-@named capacitor = Capacitor(C = C)
-@named source = ConstantVoltage(V = V)
-@named ground = Ground()
-
-rc_eqs = [connect(source.p, resistor.p)
-    connect(resistor.n, capacitor.p)
-    connect(capacitor.n, source.n)
-    connect(capacitor.n, ground.g)]
-
-@named _rc_model = ODESystem(rc_eqs, t)
-@named rc_model = compose(_rc_model,
-    [resistor, capacitor, source, ground])
+@named rc_model = RCModel(resistor.R = 2.0)
+rc_model = complete(rc_model)
 sys = structural_simplify(rc_model)
 u0 = [
-    capacitor.v => 0.0,
+    rc_model.capacitor.v => 0.0,
 ]
 prob = ODAEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Tsit5())
@@ -108,7 +119,7 @@ We wish to build the following RC circuit by building individual components and
 
 ### Building the Component Library
 
-For each of our components, we use a Julia function which emits an `ODESystem`.
+For each of our components, we use ModelingToolkit `Model` that emits an `ODESystem`.
 At the top, we start with defining the fundamental qualities of an electric
 circuit component. At every input and output pin, a circuit component has
 two values: the current at the pin and the voltage. Thus we define the `Pin`
@@ -119,40 +130,35 @@ i.e. that currents sum to zero and voltages across the pins are equal.
 default, variables are equal in a connection.
 
 ```@example acausal
-@connector function Pin(; name)
-    sts = @variables v(t)=1.0 i(t)=1.0 [connect = Flow]
-    ODESystem(Equation[], t, sts, []; name = name)
+@connector Pin begin
+    v(t) = 1.0
+    i(t) = 1.0, [connect = Flow]
 end
 ```
 
 Note that this is an incompletely specified ODESystem: it cannot be simulated
 on its own because the equations for `v(t)` and `i(t)` are unknown. Instead,
 this just gives a common syntax for receiving this pair with some default
-values. Notice that in a component, we define the `name` as a keyword argument:
-this is because later we will generate different `Pin` objects with different
-names to correspond to duplicates of this topology with unique variables.
-One can then construct a `Pin` like:
-
-```@example acausal
-Pin(name = :mypin1)
-```
-
-or equivalently, using the `@named` helper macro:
+values.
+One can then construct a `Pin` using the `@named` helper macro:
 
 ```@example acausal
 @named mypin1 = Pin()
 ```
 
 Next, we build our ground node. A ground node is just a pin that is connected
-to a constant voltage reservoir, typically taken to be `V=0`. Thus to define
+to a constant voltage reservoir, typically taken to be `V = 0`. Thus to define
 this component, we generate an `ODESystem` with a `Pin` subcomponent and specify
 that the voltage in such a `Pin` is equal to zero. This gives:
 
 ```@example acausal
-@component function Ground(; name)
-    @named g = Pin()
-    eqs = [g.v ~ 0]
-    compose(ODESystem(eqs, t, [], []; name = name), g)
+@mtkmodel Ground begin
+    @components begin
+        g = Pin()
+    end
+    @equations begin
+        g.v ~ 0
+    end
 end
 ```
 
@@ -163,14 +169,20 @@ zero, and the current of the component equals to the current of the positive
 pin.
 
 ```@example acausal
-@component function OnePort(; name)
-    @named p = Pin()
-    @named n = Pin()
-    sts = @variables v(t)=1.0 i(t)=1.0
-    eqs = [v ~ p.v - n.v
+@mtkmodel OnePort begin
+    @components begin
+        p = Pin()
+        n = Pin()
+    end
+    @variables begin
+        v(t) = 1.0
+        i(t) = 1.0
+    end
+    @equations begin
+        v ~ p.v - n.v
         0 ~ p.i + n.i
-        i ~ p.i]
-    compose(ODESystem(eqs, t, sts, []; name = name), p, n)
+        i ~ p.i
+    end
 end
 ```
 
@@ -182,38 +194,37 @@ of charge we know that the current in must equal the current out, which means
 zero. This leads to our resistor equations:
 
 ```@example acausal
-@component function Resistor(; name, R = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters R = R
-    eqs = [
-        v ~ i * R,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel Resistor begin
+    @extend v, i = oneport = OnePort()
+    @parameters begin
+        R = 1.0
+    end
+    @equations begin
+        v ~ i * R
+    end
 end
 ```
 
 Notice that we have created this system with a default parameter `R` for the
 resistor's resistance. By doing so, if the resistance of this resistor is not
 overridden by a higher level default or overridden at `ODEProblem` construction
-time, this will be the value of the resistance. Also, note the use of `@unpack`
-and `extend`. For the `Resistor`, we want to simply inherit `OnePort`'s
-equations and states and extend them with a new equation. ModelingToolkit makes
-a new namespaced variable `oneport₊v(t)` when using the syntax `oneport.v`, and
-we can use `@unpack` to avoid the namespacing.
+time, this will be the value of the resistance. Also, note the use of `@extend`.
+For the `Resistor`, we want to simply inherit `OnePort`'s
+equations and states and extend them with a new equation. Note that `v`, `i` are not namespaced as `oneport.v` or `oneport.i`.
 
 Using our knowledge of circuits, we similarly construct the `Capacitor`:
 
 ```@example acausal
-@component function Capacitor(; name, C = 1.0)
-    @named oneport = OnePort()
-    @unpack v, i = oneport
-    ps = @parameters C = C
-    D = Differential(t)
-    eqs = [
-        D(v) ~ i / C,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+D = Differential(t)
+
+@mtkmodel Capacitor begin
+    @extend v, i = oneport = OnePort()
+    @parameters begin
+        C = 1.0
+    end
+    @equations begin
+        D(v) ~ i / C
+    end
 end
 ```
 
@@ -223,55 +234,52 @@ constant voltage, essentially generating the electric current. We would then
 model this as:
 
 ```@example acausal
-@component function ConstantVoltage(; name, V = 1.0)
-    @named oneport = OnePort()
-    @unpack v = oneport
-    ps = @parameters V = V
-    eqs = [
-        V ~ v,
-    ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+@mtkmodel ConstantVoltage begin
+    @extend (v,) = oneport = OnePort()
+    @parameters begin
+        V = 1.0
+    end
+    @equations begin
+        V ~ v
+    end
 end
 ```
 
+Note that as we are extending only `v` from `OnePort`, it is explicitly specified as a tuple.
+
 ### Connecting and Simulating Our Electric Circuit
 
 Now we are ready to simulate our circuit. Let's build our four components:
 a `resistor`, `capacitor`, `source`, and `ground` term. For simplicity, we will
-make all of our parameter values 1. This is done by:
+make all of our parameter values 1.0. As `resistor`, `capacitor`, `source` lists
+`R`, `C`, `V` in their argument list, they are promoted as arguments of RCModel as
+`resistor.R`, `capacitor.C`, `source.V`
 
 ```@example acausal
-R = 1.0
-C = 1.0
-V = 1.0
-@named resistor = Resistor(R = R)
-@named capacitor = Capacitor(C = C)
-@named source = ConstantVoltage(V = V)
-@named ground = Ground()
-```
-
-Subsequently, we will connect the pieces of our circuit together. Let's connect the
-positive pin of the resistor to the source, the negative pin of the resistor
-to the capacitor, and the negative pin of the capacitor to a junction between
-the source and the ground. This would mean our connection equations are:
-
-```@example acausal
-rc_eqs = [connect(source.p, resistor.p)
-    connect(resistor.n, capacitor.p)
-    connect(capacitor.n, source.n)
-    connect(capacitor.n, ground.g)]
+@mtkmodel RCModel begin
+    @components begin
+        resistor = Resistor(R = 1.0)
+        capacitor = Capacitor(C = 1.0)
+        source = ConstantVoltage(V = 1.0)
+        ground = Ground()
+    end
+    @equations begin
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n)
+        connect(capacitor.n, ground.g)
+    end
+end
 ```
 
-Finally, we build our four-component model with these connection rules:
+We can create a RCModel component with `@named`. And using `subcomponent_name.parameter` we can set
+the parameters or defaults values of variables of subcomponents.
 
 ```@example acausal
-@named _rc_model = ODESystem(rc_eqs, t)
-@named rc_model = compose(_rc_model,
-    [resistor, capacitor, source, ground])
+@named rc_model = RCModel(resistor.R = 2.0)
 ```
 
-Note that we can also specify the subsystems in a vector. This model is acausal
-because we have not specified anything about the causality of the model. We have
+This model is acausal because we have not specified anything about the causality of the model. We have
 simply specified what is true about each of the variables. This forms a system
 of differential-algebraic equations (DAEs) which define the evolution of each
 state of the system. The equations are:
@@ -320,8 +328,9 @@ DAE solver](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/#OrdinaryD
 This is done as follows:
 
 ```@example acausal
-u0 = [capacitor.v => 0.0
-    capacitor.p.i => 0.0]
+u0 = [rc_model.capacitor.v => 0.0
+    rc_model.capacitor.p.i => 0.0]
+
 prob = ODEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
 plot(sol)
@@ -333,7 +342,7 @@ letter `A`):
 
 ```@example acausal
 u0 = [
-    capacitor.v => 0.0,
+    rc_model.capacitor.v => 0.0,
 ]
 prob = ODAEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
@@ -344,8 +353,8 @@ Notice that this solves the whole system by only solving for one variable!
 
 However, what if we wanted to plot the timeseries of a different variable? Do
 not worry, that information was not thrown away! Instead, transformations
-like `structural_simplify` simply change state variables into `observed`
-variables. Let's see what our observed variables are:
+like `structural_simplify` simply change state variables into observables which are
+defined by `observed` equations.
 
 ```@example acausal
 observed(sys)
@@ -360,11 +369,11 @@ The solution object can be accessed via its symbols. For example, let's retrieve
 the voltage of the resistor over time:
 
 ```@example acausal
-sol[resistor.v]
+sol[rc_model.resistor.v]
 ```
 
 or we can plot the timeseries of the resistor's voltage:
 
 ```@example acausal
-plot(sol, idxs = [resistor.v])
+plot(sol, idxs = [rc_model.resistor.v])
 ```

From 6ac545d10b476881695dfd9a6926090925de163c Mon Sep 17 00:00:00 2001
From: Venkateshprasad <32921645+ven-k@users.noreply.github.com>
Date: Wed, 27 Sep 2023 14:24:46 +0530
Subject: [PATCH 22/50] docs(refactor): update structural identifiability
 tutorial with `@mtkmodel`

---
 .../tutorials/parameter_identifiability.md    | 151 ++++++++++++------
 1 file changed, 102 insertions(+), 49 deletions(-)

diff --git a/docs/src/tutorials/parameter_identifiability.md b/docs/src/tutorials/parameter_identifiability.md
index 670c3397cc..eb0837360b 100644
--- a/docs/src/tutorials/parameter_identifiability.md
+++ b/docs/src/tutorials/parameter_identifiability.md
@@ -31,24 +31,40 @@ We first define the parameters, variables, differential equations and the output
 ```@example SI
 using StructuralIdentifiability, ModelingToolkit
 
-# define parameters and variables
-@variables t x4(t) x5(t) x6(t) x7(t) y1(t) y2(t)
-@parameters k5 k6 k7 k8 k9 k10
+@variables t
 D = Differential(t)
 
-# define equations
-eqs = [
-    D(x4) ~ -k5 * x4 / (k6 + x4),
-    D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6),
-    D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10,
-    D(x7) ~ k9 * x6 * (k10 - x6) / k10,
-]
-
-# define the output functions (quantities that can be measured)
-measured_quantities = [y1 ~ x4, y2 ~ x5]
+@mtkmodel Biohydrogenation begin
+    @variables begin
+        x4(t)
+        x5(t)
+        x6(t)
+        x7(t)
+        y1(t), [output = true]
+        y2(t), [output = true]
+    end
+    @parameters begin
+        k5
+        k6
+        k7
+        k8
+        k9
+        k10
+    end
+    # define equations
+    @equations begin
+        D(x4) ~ -k5 * x4 / (k6 + x4)
+        D(x5) ~ k5 * x4 / (k6 + x4) - k7 * x5 / (k8 + x5 + x6)
+        D(x6) ~ k7 * x5 / (k8 + x5 + x6) - k9 * x6 * (k10 - x6) / k10
+        D(x7) ~ k9 * x6 * (k10 - x6) / k10
+        y1 ~ x4
+        y2 ~ x5
+    end
+end
 
 # define the system
-de = ODESystem(eqs, t, name = :Biohydrogenation)
+@named de = Biohydrogenation()
+de = complete(de)
 ```
 
 After that, we are ready to check the system for local identifiability:
@@ -56,8 +72,7 @@ After that, we are ready to check the system for local identifiability:
 ```@example SI
 # query local identifiability
 # we pass the ode-system
-local_id_all = assess_local_identifiability(de, measured_quantities = measured_quantities,
-    p = 0.99)
+local_id_all = assess_local_identifiability(de, p = 0.99)
 ```
 
 We can see that all states (except $x_7$) and all parameters are locally identifiable with probability 0.99.
@@ -65,9 +80,8 @@ We can see that all states (except $x_7$) and all parameters are locally identif
 Let's try to check specific parameters and their combinations
 
 ```@example SI
-to_check = [k5, k7, k10 / k9, k5 + k6]
-local_id_some = assess_local_identifiability(de, measured_quantities = measured_quantities,
-    funcs_to_check = to_check, p = 0.99)
+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)
 ```
 
 Notice that in this case, everything (except the state variable $x_7$) is locally identifiable, including combinations such as $k_{10}/k_9, k_5+k_6$
@@ -92,26 +106,43 @@ 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.
 
-__Note__: as of writing this tutorial, UTF-symbols such as Greek characters are not supported by one of the project's dependencies, see [this issue](https://github.com/SciML/StructuralIdentifiability.jl/issues/43).
-
 ```@example SI2
 using StructuralIdentifiability, ModelingToolkit
-@parameters b c a beta g delta sigma
-@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t)
-D = Differential(t)
-
-eqs = [
-    D(x1) ~ -b * x1 + 1 / (c + x4),
-    D(x2) ~ a * x1 - beta * x2,
-    D(x3) ~ g * x2 - delta * x3,
-    D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3,
-]
 
-measured_quantities = [y ~ x1 + x2, y2 ~ x2]
-
-ode = ODESystem(eqs, t, name = :GoodwinOsc)
+@variables t
+D = Differential(t)
 
-global_id = assess_identifiability(ode, measured_quantities = measured_quantities)
+@mtkmodel GoodwinOsc begin
+    @parameters begin
+        b
+        c
+        α
+        β
+        γ
+        δ
+        σ
+    end
+    @variables begin
+        x1(t)
+        x2(t)
+        x3(t)
+        x4(t)
+        y(t), [output = true]
+        y2(t), [output = true]
+    end
+    @equations begin
+        D(x1) ~ -b * x1 + 1 / (c + x4)
+        D(x2) ~ α * x1 - β * x2
+        D(x3) ~ γ * x2 - δ * x3
+        D(x4) ~ σ * x4 * (γ * x2 - δ * x3) / x3
+        y ~ x1 + x2
+        y2 ~ x2
+    end
+end
+
+@named ode = GoodwinOsc()
+
+global_id = assess_identifiability(ode)
 ```
 
 We can see that only parameters `a, g` are unidentifiable, and everything else can be uniquely recovered.
@@ -120,25 +151,47 @@ Let us consider the same system but with two inputs, and we will find out identi
 
 ```@example SI3
 using StructuralIdentifiability, ModelingToolkit
-@parameters b c a beta g delta sigma
-@variables t x1(t) x2(t) x3(t) x4(t) y(t) y2(t) u1(t) [input = true] u2(t) [input = true]
+
+@variables t
 D = Differential(t)
 
-eqs = [
-    D(x1) ~ -b * x1 + 1 / (c + x4),
-    D(x2) ~ a * x1 - beta * x2 - u1,
-    D(x3) ~ g * x2 - delta * x3 + u2,
-    D(x4) ~ sigma * x4 * (g * x2 - delta * x3) / x3,
-]
-measured_quantities = [y ~ x1 + x2, y2 ~ x2]
+@mtkmodel GoodwinOscillator begin
+    @parameters begin
+        b
+        c
+        α
+        β
+        γ
+        δ
+        σ
+    end
+    @variables begin
+        x1(t)
+        x2(t)
+        x3(t)
+        x4(t)
+        y(t), [output = true]
+        y2(t), [output = true]
+        u1(t), [input = true]
+        u2(t), [input = true]
+    end
+    @equations begin
+        D(x1) ~ -b * x1 + 1 / (c + x4)
+        D(x2) ~ α * x1 - β * x2 - u1
+        D(x3) ~ γ * x2 - δ * x3 + u2
+        D(x4) ~ σ * x4 * (γ * x2 - δ * x3) / x3
+        y ~ x1 + x2
+        y2 ~ x2
+    end
+end
+
+@named ode = GoodwinOscillator()
+ode = complete(ode)
 
 # check only 2 parameters
-to_check = [b, c]
-
-ode = ODESystem(eqs, t, name = :GoodwinOsc)
+to_check = [ode.b, ode.c]
 
-global_id = assess_identifiability(ode, measured_quantities = measured_quantities,
-    funcs_to_check = to_check, p = 0.9)
+global_id = assess_identifiability(ode, funcs_to_check = to_check, p = 0.9)
 ```
 
 Both parameters `b, c` are globally identifiable with probability `0.9` in this case.

From 278698bad460781cdb62837ea1381f54bb0de5f4 Mon Sep 17 00:00:00 2001
From: Venkateshprasad <32921645+ven-k@users.noreply.github.com>
Date: Thu, 28 Sep 2023 12:21:43 +0530
Subject: [PATCH 23/50] docs(fix): set a `u0` within the constraint for the
 optimization problem

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

diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md
index 24b74c737b..b4ee03141b 100644
--- a/docs/src/tutorials/optimization.md
+++ b/docs/src/tutorials/optimization.md
@@ -69,8 +69,8 @@ cons = [
     x^2 + y^2 ≲ 1,
 ]
 @named sys = OptimizationSystem(loss, [x, y], [a, b], constraints = cons)
-u0 = [x => 1.0
-    y => 2.0]
+u0 = [x => 0.14
+    y => 0.14]
 prob = OptimizationProblem(sys, u0, grad = true, hess = true, cons_j = true, cons_h = true)
 solve(prob, IPNewton())
 ```

From 94fb749315d35f89ebfe4de6c0b06bd6d27c6a80 Mon Sep 17 00:00:00 2001
From: Venkateshprasad <32921645+ven-k@users.noreply.github.com>
Date: Tue, 3 Oct 2023 15:23:22 +0530
Subject: [PATCH 24/50] docs(refactor): refactor "Getting Started" section with
 `@mtkmodel`

- All sections, except building hierarchal models, use `@mtkmodel`.
- "Non-DSL way of defining an ODESystem" section is added.
---
 docs/src/tutorials/ode_modeling.md | 221 ++++++++++++++++++++++++-----
 1 file changed, 184 insertions(+), 37 deletions(-)

diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index 5763eeabf6..dd5875da7c 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -9,21 +9,34 @@ illustrate the basic user-facing functionality.
 A much deeper tutorial with forcing functions and sparse Jacobians is below.
 But if you want to just see some code and run, here's an example:
 
-```@example ode
+```@example first-mtkmodel
 using ModelingToolkit
 
-@variables t x(t)   # independent and dependent variables
-@parameters τ       # parameters
-@constants h = 1    # constants have an assigned value
-D = Differential(t) # define an operator for the differentiation w.r.t. time
-
-# your first ODE, consisting of a single equation, the equality indicated by ~
-@named fol = ODESystem([D(x) ~ (h - x) / τ])
+@variables t
+D = Differential(t)
+
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+    end
+    @structural_parameters begin
+        h = 1
+    end
+    @equations begin
+        D(x) ~ (h - x) / τ
+    end
+end
 
 using DifferentialEquations: solve
 
-prob = ODEProblem(fol, [x => 0.0], (0.0, 10.0), [τ => 3.0])
-# parameter `τ` can be assigned a value, but constant `h` cannot
+@named fol = FOL()
+fol = complete(fol)
+
+prob = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
+# parameter `τ` can be assigned a value, but structural parameter `h` cannot'.
 sol = solve(prob)
 
 using Plots
@@ -43,24 +56,50 @@ first-order lag element:
 
 Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) state
 variable, ``f(t)`` is an external forcing function, and ``\tau`` is a
-parameter. In MTK, this system can be modelled as follows. For simplicity, we
-first set the forcing function to a time-independent value.
+parameter.
+In MTK, this system can be modelled as follows. For simplicity, we
+first set the forcing function to a time-independent value ``h``. And the
+independent variable ``t`` is automatically added by ``@mtkmodel``.
 
 ```@example ode2
 using ModelingToolkit
 
-@variables t x(t)  # independent and dependent variables
-@parameters τ       # parameters
-@constants h = 1    # constants
-D = Differential(t) # define an operator for the differentiation w.r.t. time
+@variables t
+D = Differential(t)
+
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+    end
+    @structural_parameters begin
+        h = 1
+    end
+    @equations begin
+        D(x) ~ (h - x) / τ
+    end
+end
 
-# your first ODE, consisting of a single equation, indicated by ~
-@named fol_model = ODESystem(D(x) ~ (h - x) / τ)
+@named fol_incomplete = FOL()
+fol = complete(fol_incomplete)
 ```
 
 Note that equations in MTK use the tilde character (`~`) as equality sign.
-Also note that the `@named` macro simply ensures that the symbolic name
-matches the name in the REPL. If omitted, you can directly set the `name` keyword.
+
+`@named` creates an instance of `FOL` named as `fol`. Before creating an
+ODEProblem with `fol` run `complete`. Once the system is complete, it will no
+longer namespace its subsystems or variables. This is necessary to correctly pass
+the intial values of states and parameters to the ODEProblem.
+
+```julia
+julia> fol_incomplete.x
+fol_incomplete₊x(t)
+
+julia> fol.x
+x(t)
+```
 
 After construction of the ODE, you can solve it using [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/):
 
@@ -68,7 +107,7 @@ After construction of the ODE, you can solve it using [DifferentialEquations.jl]
 using DifferentialEquations
 using Plots
 
-prob = ODEProblem(fol_model, [x => 0.0], (0.0, 10.0), [τ => 3.0])
+prob = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
 plot(solve(prob))
 ```
 
@@ -76,15 +115,51 @@ The initial state and the parameter values are specified using a mapping
 from the actual symbolic elements to their values, represented as an array
 of `Pair`s, which are constructed using the `=>` operator.
 
+## Non-DSL way of defining an ODESystem
+
+Using `@mtkmodel` is the preferred way of defining ODEs with MTK. However, let us
+look at how we can define the same system without `@mtkmodel`. This is useful for
+defining PDESystem etc.
+
+```@example first-mtkmodel
+@variables t x(t)   # independent and dependent variables
+@parameters τ       # parameters
+@constants h = 1    # constants
+D = Differential(t) # define an operator for the differentiation w.r.t. time
+
+# your first ODE, consisting of a single equation, indicated by ~
+@named fol_model = ODESystem(D(x) ~ (h - x) / τ)
+```
+
 ## Algebraic relations and structural simplification
 
 You could separate the calculation of the right-hand side, by introducing an
 intermediate variable `RHS`:
 
 ```@example ode2
-@variables RHS(t)
-@named fol_separate = ODESystem([RHS ~ (h - x) / τ,
-    D(x) ~ RHS])
+using ModelingToolkit
+
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+        RHS(t)
+    end
+    @structural_parameters begin
+        h = 1
+    end
+    begin
+        D = Differential(t)
+    end
+    @equations begin
+        RHS ~ (h - x) / τ
+        D(x) ~ RHS
+    end
+end
+
+@named fol_separate = FOL()
 ```
 
 To directly solve this system, you would have to create a Differential-Algebraic
@@ -94,12 +169,12 @@ transformed into the single ODE we used in the first example above. MTK achieves
 this by structural simplification:
 
 ```@example ode2
-fol_simplified = structural_simplify(fol_separate)
+fol_simplified = structural_simplify(complete(fol_separate))
 equations(fol_simplified)
 ```
 
 ```@example ode2
-equations(fol_simplified) == equations(fol_model)
+equations(fol_simplified) == equations(fol)
 ```
 
 You can extract the equations from a system using `equations` (and, in the same
@@ -114,9 +189,12 @@ along with the state variable. Note that this has to be requested explicitly,
 through:
 
 ```@example ode2
-prob = ODEProblem(fol_simplified, [x => 0.0], (0.0, 10.0), [τ => 3.0])
+prob = ODEProblem(fol_simplified,
+    [fol_simplified.x => 0.0],
+    (0.0, 10.0),
+    [fol_simplified.τ => 3.0])
 sol = solve(prob)
-plot(sol, vars = [x, RHS])
+plot(sol, vars = [fol_simplified.x, fol_simplified.RHS])
 ```
 
 By default, `structural_simplify` also replaces symbolic `constants` with
@@ -136,8 +214,27 @@ What if the forcing function (the “external input”) ``f(t)`` is not constant
 Obviously, one could use an explicit, symbolic function of time:
 
 ```@example ode2
-@variables f(t)
-@named fol_variable_f = ODESystem([f ~ sin(t), D(x) ~ (f - x) / τ])
+@mtkmodel FOL begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+        f(t)
+    end
+    @structural_parameters begin
+        h = 1
+    end
+    begin
+        D = Differential(t)
+    end
+    @equations begin
+        f ~ sin(t)
+        D(x) ~ (f - x) / τ
+    end
+end
+
+@named fol_variable_f = FOL()
 ```
 
 But often this function might not be available in an explicit form.
@@ -154,11 +251,35 @@ value_vector = randn(10)
 f_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t)) + 1]
 @register_symbolic f_fun(t)
 
-@named fol_external_f = ODESystem([f ~ f_fun(t), D(x) ~ (f - x) / τ])
-prob = ODEProblem(structural_simplify(fol_external_f), [x => 0.0], (0.0, 10.0), [τ => 0.75])
+@mtkmodel FOLExternalFunction begin
+    @parameters begin
+        τ # parameters
+    end
+    @variables begin
+        x(t) # dependent variables
+        f(t)
+    end
+    @structural_parameters begin
+        h = 1
+    end
+    begin
+        D = Differential(t)
+    end
+    @equations begin
+        f ~ f_fun(t)
+        D(x) ~ (f - x) / τ
+    end
+end
+
+@named fol_external_f = FOLExternalFunction()
+fol_external_f = complete(fol_external_f)
+prob = ODEProblem(structural_simplify(fol_external_f),
+    [fol_external_f.x => 0.0],
+    (0.0, 10.0),
+    [fol_external_f.τ => 0.75])
 
 sol = solve(prob)
-plot(sol, vars = [x, f])
+plot(sol, vars = [fol_external_f.x, fol_external_f.f])
 ```
 
 ## Building component-based, hierarchical models
@@ -235,19 +356,43 @@ plot(solve(prob))
 
 More on this topic may be found in [Composing Models and Building Reusable Components](@ref acausal).
 
-## Defaults
+## Inital Guess
 
 It is often a good idea to specify reasonable values for the initial state and the
 parameters of a model component. Then, these do not have to be explicitly specified when constructing the `ODEProblem`.
 
 ```@example ode2
-function unitstep_fol_factory(; name)
+@mtkmodel UnitstepFOLFactory begin
+    @parameters begin
+        τ = 1.0
+    end
+    @variables begin
+        x(t) = 0.0
+    end
+    @equations begin
+        D(x) ~ (1 - x) / τ
+    end
+end
+```
+
+While defining the model `UnitstepFOLFactory`, an initial guess of 0.0 is assigned to `x(t)` and 1.0 to `τ`.
+Additionaly, these initial guesses can be modified while creating instances of `UnitstepFOLFactory` by passing arguements.
+
+```@example ode2
+@named fol = UnitstepFOLFactory(; x = 0.1)
+sol = ODEProblem(fol, [], (0.0, 5.0), []) |> solve
+```
+
+In non-DSL definitions, one can pass `defaults` dictionary to set the initial guess of the symbolic variables.
+
+```@example ode3
+using ModelingToolkit
+
+function UnitstepFOLFactory(; name)
     @parameters τ
     @variables t x(t)
     ODESystem(D(x) ~ (1 - x) / τ; name, defaults = Dict(x => 0.0, τ => 1.0))
 end
-
-ODEProblem(unitstep_fol_factory(name = :fol), [], (0.0, 5.0), []) |> solve
 ```
 
 Note that the defaults can be functions of the other variables, which is then
@@ -316,6 +461,8 @@ Where to go next?
 
   - Not sure how MTK relates to similar tools and packages? Read
     [Comparison of ModelingToolkit vs Equation-Based and Block Modeling Languages](@ref).
+  - For a more detailed explanation of `@mtkmodel` checkout
+    [Defining components with `@mtkmodel` and connectors with `@connectors`](@ref mtkmodel_connector)
   - Depending on what you want to do with MTK, have a look at some of the other
     **Symbolic Modeling Tutorials**.
   - If you want to automatically convert an existing function to a symbolic

From 438a22547d889165185ec573dbba8e6fc4872af1 Mon Sep 17 00:00:00 2001
From: Venkateshprasad <32921645+ven-k@users.noreply.github.com>
Date: Tue, 3 Oct 2023 15:24:52 +0530
Subject: [PATCH 25/50] docs(refactor): rename the `@mtkmodel` docs page as
 "Defining components with `@mtkmodel` and connectors with `@connectors`" to
 indicate that it describes `@connector`s too.

---
 docs/src/basics/MTKModel_Connector.md | 21 +++++++++++++++++++--
 1 file changed, 19 insertions(+), 2 deletions(-)

diff --git a/docs/src/basics/MTKModel_Connector.md b/docs/src/basics/MTKModel_Connector.md
index ac285d0253..902780366d 100644
--- a/docs/src/basics/MTKModel_Connector.md
+++ b/docs/src/basics/MTKModel_Connector.md
@@ -1,6 +1,23 @@
-# Defining components with `@mtkmodel`
+# [Defining components with `@mtkmodel` and connectors with `@connectors`](@id mtkmodel-connectors)
 
-`@mtkmodel` is a convenience macro to define ModelingToolkit components. It returns `ModelingToolkit.Model`, which includes a constructor that returns an ODESystem, a `structure` dictionary with metadata and flag `isconnector` which is set to `false`.
+## MTK Model
+
+MTK represents components and connectors with `Model`.
+
+```@docs
+ModelingToolkit.Model
+```
+
+## Components
+
+Components are models from various domains. These models contain states and their
+equations.
+
+### [Defining components with `@mtkmodel`](@id mtkmodel)
+
+`@mtkmodel` is a convenience macro to define components. It returns
+`ModelingToolkit.Model`, which includes a constructor that returns the ODESystem, a
+`structure` dictionary with metadata, and flag `isconnector` which is set to `false`.
 
 ## What can an MTK-Model definition have?
 

From ff2a778798db46f0691d94f357e4c7a0e11ab25c Mon Sep 17 00:00:00 2001
From: Venkateshprasad <32921645+ven-k@users.noreply.github.com>
Date: Tue, 3 Oct 2023 20:07:22 +0530
Subject: [PATCH 26/50] docs: reword Components and Connectors and elaborate
 Connectors and Model.structure.

---
 docs/src/basics/MTKModel_Connector.md | 68 ++++++++++++++++++++++-----
 1 file changed, 56 insertions(+), 12 deletions(-)

diff --git a/docs/src/basics/MTKModel_Connector.md b/docs/src/basics/MTKModel_Connector.md
index 902780366d..cee466eb37 100644
--- a/docs/src/basics/MTKModel_Connector.md
+++ b/docs/src/basics/MTKModel_Connector.md
@@ -1,4 +1,4 @@
-# [Defining components with `@mtkmodel` and connectors with `@connectors`](@id mtkmodel-connectors)
+# [Components and Connectors](@id mtkmodel_connector)
 
 ## MTK Model
 
@@ -73,9 +73,9 @@ end
 
   - Parameters and variables are declared with respective begin blocks.
   - Variables must be functions of an independent variable.
-  - Optionally, default values and metadata can be specified for these parameters and variables. See `ModelB` in the above example.
+  - Optionally, initial guess and metadata can be specified for these parameters and variables. See `ModelB` in the above example.
   - Along with creating parameters and variables, keyword arguments of same name with default value `nothing` are created.
-  - Whenever a parameter or variable has default value, for example `v(t) = 0.0`, a symbolic variable named `v` with default value 0.0 and a keyword argument `v`, with default value `nothing` are created. <br> This way, users can optionally pass new value of `v` while creating a component.
+  - Whenever a parameter or variable has initial value, for example `v(t) = 0.0`, a symbolic variable named `v` with initial value 0.0 and a keyword argument `v`, with default value `nothing` are created. <br> This way, users can optionally pass new value of `v` while creating a component.
 
 ```julia
 julia > @named model_c = ModelC(; v = 2.0);
@@ -89,7 +89,7 @@ julia > ModelingToolkit.getdefault(model_c.v)
   - This block is for non symbolic input arguements. These are for inputs that usually are not meant to be part of components; but influence how they are defined. One can list inputs like boolean flags, functions etc... here.
   - Whenever default values are specified, unlike parameters/variables, they are reflected in the keyword argument list.
 
-### `@extend` block
+#### `@extend` begin block
 
 To extend a partial system,
 
@@ -103,7 +103,8 @@ julia> @named model_c = ModelC(; p1 = 2.0)
 
 ```
 
-However, as `p2` isn't listed in the model definition, its default can't be modified by users.
+However, as `p2` isn't listed in the model definition, its initial guess can't
+specified while creating an instance of `ModelC`.
 
 ### `@components` begin block
 
@@ -127,13 +128,26 @@ And as `k2` isn't listed in the sub-component definition of `ModelC`, its defaul
 
   - Any other Julia operations can be included with dedicated begin blocks.
 
-## Defining connectors with `@connector`
+## Connectors
 
-`@connector` returns `ModelingToolkit.Model`. It includes a constructor that returns a connector ODESystem, a `structure` dictionary with metadata and flag `isconnector` which is set to `true`.
+Connectors are special models that can be used to connect different components together.
+MTK provides 3 distinct connectors:
+
+  - `DomainConnector`: A connector which has only one state which is of `Flow` type,
+    specified by `[connect = Flow]`.
+  - `StreamConnector`: A connector which has atleast one stream variable, specified by
+    `[connect = Stream]`. A `StreamConnector` must have exactly one flow variable.
+  - `RegularConnector`: Connectors that don't fall under above categories.
+
+### [Defining connectors with `@connector`](@id connector)
+
+`@connector` returns `ModelingToolkit.Model`. It includes a constructor that returns
+a connector ODESystem, a `structure` dictionary with metadata, and flag `isconnector`
+which is set to `true`.
 
 A simple connector can be defined with syntax similar to following example:
 
-```julia
+```@example connector
 using ModelingToolkit
 
 @connector Pin begin
@@ -142,17 +156,47 @@ using ModelingToolkit
 end
 ```
 
-  - Variables (as function of independent variable) are listed out in the definition. These variables can optionally have default values and metadata like `descrption`, `connect` and so on.
+Variables (as functions of independent variable) are listed out in the definition. These variables can optionally have initial values and metadata like `description`, `connect` and so on. For more details on setting metadata, check out [Symbolic Metadata](@id symbolic_metadata).
 
-`@connector`s accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`.
+Similar to `@mtkmodel`, `@connector` accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`.
+For example, the following `HydraulicFluid` connector is defined with parameters, variables and equations.
+
+```@example connector
+@connector HydraulicFluid begin
+    @parameters begin
+        ρ = 997
+        β = 2.09e9
+        μ = 0.0010016
+        n = 1
+        let_gas = 1
+        ρ_gas = 0.0073955
+        p_gas = -1000
+    end
+    @variables begin
+        dm(t) = 0.0, [connect = Flow]
+    end
+    @equations begin
+        dm ~ 0
+    end
+end
+```
 
 !!! note
     
     For more examples of usage, checkout [ModelingToolkitStandardLibrary.jl](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/)
 
-### What's a `structure` dictionary?
+## More on `Model.structure`
+
+`structure` stores metadata that describes composition of a model. It includes:
 
-For components defined with `@mtkmodel` or `@connector`, a dictionary with metadata is created. It lists `:components` (sub-component list), `:extend` (the extended states and base system), `:parameters`, `:variables`, ``:kwargs`` (list of keyword arguments), `:independent_variable`, `:equations`.
+  - `:components`: List of sub-components in the form of [[name, sub_component_name],...].
+  - `:extend`: The list of extended states, name given to the base system, and name of the base system.
+  - `:structural_parameters`: Dictionary of structural parameters mapped to their default values.
+  - `:parameters`: Dictionary of symbolic parameters mapped to their metadata.
+  - `:variables`: Dictionary of symbolic variables mapped to their metadata.
+  - `:kwargs`: Dictionary of keyword arguments mapped to their default values.
+  - `:independent_variable`: Independent variable, which is added while generating the Model.
+  - `:equations`: List of equations (represented as strings).
 
 For example, the structure of `ModelC` is:
 

From 2ef4752004276c1f04b8f61302a2bf11706769c7 Mon Sep 17 00:00:00 2001
From: Venkateshprasad <32921645+ven-k@users.noreply.github.com>
Date: Tue, 3 Oct 2023 20:10:10 +0530
Subject: [PATCH 27/50] docs: add `connect` section to Variable metadata

- while here, update the description on setting unit metadata in   Contextual-Variables
---
 docs/src/basics/ContextualVariables.md   | 12 +++++++-----
 docs/src/basics/MTKModel_Connector.md    |  2 +-
 docs/src/basics/Variable_metadata.md     | 18 +++++++++++++++++-
 docs/src/tutorials/acausal_components.md |  2 +-
 4 files changed, 26 insertions(+), 8 deletions(-)

diff --git a/docs/src/basics/ContextualVariables.md b/docs/src/basics/ContextualVariables.md
index dcaebe95d7..de2802ce40 100644
--- a/docs/src/basics/ContextualVariables.md
+++ b/docs/src/basics/ContextualVariables.md
@@ -48,7 +48,7 @@ for which its arguments must be specified each time it is used. This is useful w
 PDEs for example, where one may need to use `u(t, x)` in the equations, but will
 need to be able to write `u(t, 0.0)` to define a boundary condition at `x = 0`.
 
-## Variable metadata [Experimental/TODO]
+## Variable metadata
 
 In many engineering systems, some variables act like “flows” while others do not.
 For example, in circuit models you have current which flows, and the related
@@ -69,15 +69,17 @@ the metadata. One can get and set metadata by
 ```julia
 julia> @variables x [unit = u"m^3/s"];
 
-julia> hasmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+julia> hasmetadata(x, VariableUnit)
 true
 
-julia> getmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+julia> ModelingToolkit.get_unit(x)
 m³ s⁻¹
 
-julia> x = setmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)), u"m/s")
+julia> x = setmetadata(x, VariableUnit, u"m/s")
 x
 
-julia> getmetadata(x, Symbolics.option_to_metadata_type(Val(:unit)))
+julia> ModelingToolkit.get_unit(x)
 m s⁻¹
 ```
+
+See [Symbolic Metadata](@ref symbolic_metadata) for more details on variable metadata.
diff --git a/docs/src/basics/MTKModel_Connector.md b/docs/src/basics/MTKModel_Connector.md
index cee466eb37..9645260ed6 100644
--- a/docs/src/basics/MTKModel_Connector.md
+++ b/docs/src/basics/MTKModel_Connector.md
@@ -156,7 +156,7 @@ using ModelingToolkit
 end
 ```
 
-Variables (as functions of independent variable) are listed out in the definition. These variables can optionally have initial values and metadata like `description`, `connect` and so on. For more details on setting metadata, check out [Symbolic Metadata](@id symbolic_metadata).
+Variables (as functions of independent variable) are listed out in the definition. These variables can optionally have initial values and metadata like `description`, `connect` and so on. For more details on setting metadata, check out [Symbolic Metadata](@ref symbolic_metadata).
 
 Similar to `@mtkmodel`, `@connector` accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`.
 For example, the following `HydraulicFluid` connector is defined with parameters, variables and equations.
diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/basics/Variable_metadata.md
index 53b8eb1fcf..d186cdf959 100644
--- a/docs/src/basics/Variable_metadata.md
+++ b/docs/src/basics/Variable_metadata.md
@@ -1,4 +1,4 @@
-# Symbolic metadata
+# [Symbolic Metadata](@id symbolic_metadata)
 
 It is possible to add metadata to symbolic variables, the metadata will be displayed when calling help on a variable.
 
@@ -39,6 +39,22 @@ help?> u
   Symbolics.VariableSource: (:variables, :u)
 ```
 
+## Connect
+
+Variables in connectors can have `connect` metadata which describes the type of connections.
+
+`Flow` is used for variables that represent physical quantities that "flow" ex:
+current in a resistor. These variables sum up to zero in connections.
+
+`Stream` can be specified for variables that flow bi-directionally.
+
+```@example connect
+using ModelingToolkit
+
+@variables t, i(t) [connect = Flow]
+@variables k(t) [connect = Stream]
+```
+
 ## Input or output
 
 Designate a variable as either an input or an output using the following
diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index b920c07df3..c2a469eea3 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -10,7 +10,7 @@ component variables. We then simplify this to an ODE by eliminating the
 equalities before solving. Let's see this in action.
 
 !!! note
-
+    
     This tutorial teaches how to build the entire RC circuit from scratch.
     However, to simulate electric components with more ease, check out the
     [ModelingToolkitStandardLibrary.jl](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/)

From c4868e3093d22912148b18e8d74ea18d3bf6bf0d Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 8 Oct 2023 19:26:35 +0200
Subject: [PATCH 28/50] Add `@mtkbuild` greatly simplify the tutorials
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Good thing is that the tutorials are much simpler.

Bad thing is that there's an issue with completing a structurally simplified model which only allows for using the flattened variables, and thus requires using:

```julia
u0 = [
    rc_model.capacitor₊v => 0.0,
]
```

by the user, which is not nice as it should be:

```julia
u0 = [
    rc_model.capacitor.v => 0.0,
]
```

which is something we should probably fix ASAP
---
 docs/pages.jl                                 |   1 +
 docs/src/tutorials/acausal_components.md      |  59 +++----
 docs/src/tutorials/ode_modeling.md            | 145 +++++++-----------
 .../tutorials/programmatically_generating.md  | 106 +++++++++++++
 src/ModelingToolkit.jl                        |   2 +-
 src/systems/abstractsystem.jl                 |   9 ++
 6 files changed, 196 insertions(+), 126 deletions(-)
 create mode 100644 docs/src/tutorials/programmatically_generating.md

diff --git a/docs/pages.jl b/docs/pages.jl
index 0b6d0555cd..f6f5d7af50 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -5,6 +5,7 @@ pages = [
         "tutorials/nonlinear.md",
         "tutorials/optimization.md",
         "tutorials/modelingtoolkitize.md",
+        "tutorials/programmatically_generating.md",
         "tutorials/stochastic_diffeq.md",
         "tutorials/parameter_identifiability.md",
         "tutorials/domain_connections.md"],
diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index c2a469eea3..4a7e4ddd1c 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -24,8 +24,8 @@ using ModelingToolkit, Plots, DifferentialEquations
 
 @variables t
 @connector Pin begin
-    v(t) = 1.0
-    i(t) = 1.0, [connect = Flow]
+    v(t)
+    i(t), [connect = Flow]
 end
 
 @mtkmodel Ground begin
@@ -43,8 +43,8 @@ end
         n = Pin()
     end
     @variables begin
-        v(t) = 1.0
-        i(t) = 1.0
+        v(t)
+        i(t)
     end
     @equations begin
         v ~ p.v - n.v
@@ -56,7 +56,7 @@ end
 @mtkmodel Resistor begin
     @extend v, i = oneport = OnePort()
     @parameters begin
-        R = 1.0
+        R = 1.0 # Sets the default resistance
     end
     @equations begin
         v ~ i * R
@@ -100,13 +100,11 @@ end
     end
 end
 
-@named rc_model = RCModel(resistor.R = 2.0)
-rc_model = complete(rc_model)
-sys = structural_simplify(rc_model)
+@mtkbuild rc_model = RCModel(resistor.R = 2.0)
 u0 = [
-    rc_model.capacitor.v => 0.0,
+    rc_model.capacitor₊v => 0.0,
 ]
-prob = ODAEProblem(sys, u0, (0, 10.0))
+prob = ODAEProblem(rc_model, u0, (0, 10.0))
 sol = solve(prob, Tsit5())
 plot(sol)
 ```
@@ -131,8 +129,8 @@ default, variables are equal in a connection.
 
 ```@example acausal
 @connector Pin begin
-    v(t) = 1.0
-    i(t) = 1.0, [connect = Flow]
+    v(t)
+    i(t), [connect = Flow]
 end
 ```
 
@@ -175,8 +173,8 @@ pin.
         n = Pin()
     end
     @variables begin
-        v(t) = 1.0
-        i(t) = 1.0
+        v(t)
+        i(t)
     end
     @equations begin
         v ~ p.v - n.v
@@ -276,7 +274,7 @@ We can create a RCModel component with `@named`. And using `subcomponent_name.pa
 the parameters or defaults values of variables of subcomponents.
 
 ```@example acausal
-@named rc_model = RCModel(resistor.R = 2.0)
+@mtkbuild rc_model = RCModel(resistor.R = 2.0)
 ```
 
 This model is acausal because we have not specified anything about the causality of the model. We have
@@ -300,26 +298,15 @@ and the parameters are:
 parameters(rc_model)
 ```
 
-## Simplifying and Solving this System
-
-This system could be solved directly as a DAE using [one of the DAE solvers
-from DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/).
-However, let's take a second to symbolically simplify the system before doing the
-solve. Although we can use ODE solvers that handle mass matrices to solve the
-above system directly, we want to run the `structural_simplify` function first,
-as it eliminates many unnecessary variables to build the leanest numerical
-representation of the system. Let's see what it does here:
+The observed equations are:
 
 ```@example acausal
-sys = structural_simplify(rc_model)
-equations(sys)
+observed(rc_model)
 ```
 
-```@example acausal
-states(sys)
-```
+## Solving this System
 
-After structural simplification, we are left with a system of only two equations
+We are left with a system of only two equations
 with two state variables. One of the equations is a differential equation,
 while the other is an algebraic equation. We can then give the values for the
 initial conditions of our states, and solve the system by converting it to
@@ -328,21 +315,21 @@ DAE solver](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/#OrdinaryD
 This is done as follows:
 
 ```@example acausal
-u0 = [rc_model.capacitor.v => 0.0
-    rc_model.capacitor.p.i => 0.0]
+u0 = [rc_model.capacitor₊v => 0.0
+    rc_model.capacitor₊p₊i => 0.0]
 
 prob = ODEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
 plot(sol)
 ```
 
-Since we have run `structural_simplify`, MTK can numerically solve all the
+MTK can numerically solve all the
 unreduced algebraic equations using the `ODAEProblem` (note the
 letter `A`):
 
 ```@example acausal
 u0 = [
-    rc_model.capacitor.v => 0.0,
+    rc_model.capacitor₊v => 0.0,
 ]
 prob = ODAEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
@@ -369,11 +356,11 @@ The solution object can be accessed via its symbols. For example, let's retrieve
 the voltage of the resistor over time:
 
 ```@example acausal
-sol[rc_model.resistor.v]
+sol[rc_model.resistor₊v]
 ```
 
 or we can plot the timeseries of the resistor's voltage:
 
 ```@example acausal
-plot(sol, idxs = [rc_model.resistor.v])
+plot(sol, idxs = [rc_model.resistor₊v])
 ```
diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index dd5875da7c..b318a6c955 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -1,8 +1,17 @@
 # Getting Started with ModelingToolkit.jl
 
-This is an introductory tutorial for ModelingToolkit (MTK).
-Some examples of Ordinary Differential Equations (ODE) are used to
-illustrate the basic user-facing functionality.
+This is an introductory tutorial for ModelingToolkit (MTK). We will demonstrate
+the basics of the package by demonstrating how to define and simulate simple
+Ordinary Differential Equation (ODE) systems.
+
+## Installing ModelingToolkit
+
+To install ModelingToolkit, use the Julia package manager. This can be done as follows:
+
+```julia
+using Pkg
+Pkg.add("ModelingToolkit")
+```
 
 ## Copy-Pastable Simplified Example
 
@@ -22,21 +31,14 @@ D = Differential(t)
     @variables begin
         x(t) # dependent variables
     end
-    @structural_parameters begin
-        h = 1
-    end
     @equations begin
-        D(x) ~ (h - x) / τ
+        D(x) ~ (1 - x) / τ
     end
 end
 
 using DifferentialEquations: solve
-
-@named fol = FOL()
-fol = complete(fol)
-
+@mtkbuild fol = FOL()
 prob = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
-# parameter `τ` can be assigned a value, but structural parameter `h` cannot'.
 sol = solve(prob)
 
 using Plots
@@ -58,7 +60,7 @@ Here, ``t`` is the independent variable (time), ``x(t)`` is the (scalar) state
 variable, ``f(t)`` is an external forcing function, and ``\tau`` is a
 parameter.
 In MTK, this system can be modelled as follows. For simplicity, we
-first set the forcing function to a time-independent value ``h``. And the
+first set the forcing function to a time-independent value ``1``. And the
 independent variable ``t`` is automatically added by ``@mtkmodel``.
 
 ```@example ode2
@@ -74,32 +76,17 @@ D = Differential(t)
     @variables begin
         x(t) # dependent variables
     end
-    @structural_parameters begin
-        h = 1
-    end
     @equations begin
-        D(x) ~ (h - x) / τ
+        D(x) ~ (1 - x) / τ
     end
 end
 
-@named fol_incomplete = FOL()
-fol = complete(fol_incomplete)
+@mtkbuild fol = FOL()
 ```
 
 Note that equations in MTK use the tilde character (`~`) as equality sign.
 
-`@named` creates an instance of `FOL` named as `fol`. Before creating an
-ODEProblem with `fol` run `complete`. Once the system is complete, it will no
-longer namespace its subsystems or variables. This is necessary to correctly pass
-the intial values of states and parameters to the ODEProblem.
-
-```julia
-julia> fol_incomplete.x
-fol_incomplete₊x(t)
-
-julia> fol.x
-x(t)
-```
+`@mtkbuild` creates an instance of `FOL` named as `fol`. 
 
 After construction of the ODE, you can solve it using [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/):
 
@@ -115,22 +102,6 @@ The initial state and the parameter values are specified using a mapping
 from the actual symbolic elements to their values, represented as an array
 of `Pair`s, which are constructed using the `=>` operator.
 
-## Non-DSL way of defining an ODESystem
-
-Using `@mtkmodel` is the preferred way of defining ODEs with MTK. However, let us
-look at how we can define the same system without `@mtkmodel`. This is useful for
-defining PDESystem etc.
-
-```@example first-mtkmodel
-@variables t x(t)   # independent and dependent variables
-@parameters τ       # parameters
-@constants h = 1    # constants
-D = Differential(t) # define an operator for the differentiation w.r.t. time
-
-# your first ODE, consisting of a single equation, indicated by ~
-@named fol_model = ODESystem(D(x) ~ (h - x) / τ)
-```
-
 ## Algebraic relations and structural simplification
 
 You could separate the calculation of the right-hand side, by introducing an
@@ -147,46 +118,40 @@ using ModelingToolkit
         x(t) # dependent variables
         RHS(t)
     end
-    @structural_parameters begin
-        h = 1
-    end
     begin
         D = Differential(t)
     end
     @equations begin
-        RHS ~ (h - x) / τ
+        RHS ~ (1 - x) / τ
         D(x) ~ RHS
     end
 end
 
-@named fol_separate = FOL()
+@mtkbuild fol = FOL()
 ```
 
-To directly solve this system, you would have to create a Differential-Algebraic
-Equation (DAE) problem, since besides the differential equation, there is an
-additional algebraic equation now. However, this DAE system can obviously be
-transformed into the single ODE we used in the first example above. MTK achieves
-this by structural simplification:
+You can look at the equations by using the command `equations`:
 
 ```@example ode2
-fol_simplified = structural_simplify(complete(fol_separate))
-equations(fol_simplified)
+equations(fol)
 ```
 
+Notice that there is only one equation in this system, `Differential(t)(x(t)) ~ RHS(t)`.
+The other equation was removed from the system and was transformed into an `observed`
+variable. Observed equations are variables which can be computed on-demand but are not
+necessary for the solution of the system, and thus MTK tracks it separately. One can
+check the observed equations via the `observed` function:
+
 ```@example ode2
-equations(fol_simplified) == equations(fol)
+observed(fol)
 ```
 
-You can extract the equations from a system using `equations` (and, in the same
-way, `states` and `parameters`). The simplified equation is exactly the same
-as the original one, so the simulation performance will also be the same.
-However, there is one difference. MTK does keep track of the eliminated
-algebraic variables as "observables" (see
-[Observables and Variable Elimination](@ref)).
-That means, MTK still knows how to calculate them out of the information available
+For more information on this process, see [Observables and Variable Elimination](@ref).
+
+MTK still knows how to calculate them out of the information available
 in a simulation result. The intermediate variable `RHS` therefore can be plotted
-along with the state variable. Note that this has to be requested explicitly,
-through:
+along with the state variable. Note that this has to be requested explicitly
+like as follows:
 
 ```@example ode2
 prob = ODEProblem(fol_simplified,
@@ -197,16 +162,26 @@ sol = solve(prob)
 plot(sol, vars = [fol_simplified.x, fol_simplified.RHS])
 ```
 
-By default, `structural_simplify` also replaces symbolic `constants` with
-their default values. This allows additional simplifications not possible
-when using `parameters` (e.g., solution of linear equations by dividing out
-the constant's value, which cannot be done for parameters, since they may
-be zero).
+## Named Indexing of Solutions
+
+Note that the indexing of the solution similarly works via the names, and so to get
+the time series for `x`, one would do:
+
+```@example ode2
+sol[fol.x]
+```
+
+or to get the second value in the time series for `x`:
+
+```@example ode2
+sol[fol.x,2]
+```
 
-Note that the indexing of the solution similarly works via the names, and so
-`sol[x]` gives the time-series for `x`, `sol[x,2:10]` gives the 2nd through 10th
-values of `x` matching `sol.t`, etc. Note that this works even for variables
-which have been eliminated, and thus `sol[RHS]` retrieves the values of `RHS`.
+Similarly, the time series for `RHS` can be retrieved using the same indexing:
+
+```@example ode2
+sol[fol.RHS]
+```
 
 ## Specifying a time-variable forcing function
 
@@ -222,9 +197,6 @@ Obviously, one could use an explicit, symbolic function of time:
         x(t) # dependent variables
         f(t)
     end
-    @structural_parameters begin
-        h = 1
-    end
     begin
         D = Differential(t)
     end
@@ -445,14 +417,9 @@ using the structural information. For more information, see the
 
 Here are some notes that may be helpful during your initial steps with MTK:
 
-  - Sometimes, the symbolic engine within MTK cannot correctly identify the
-    independent variable (e.g. time) out of all variables. In such a case, you
-    usually get an error that some variable(s) is "missing from variable map". In
-    most cases, it is then sufficient to specify the independent variable as second
-    argument to `ODESystem`, e.g. `ODESystem(eqs, t)`.
-  - A completely macro-free usage of MTK is possible and is discussed in a
-    separate tutorial. This is for package developers, since the macros are only
-    essential for automatic symbolic naming for modelers.
+  - The `@mtkmodel` macro is for high-level usage of MTK. However, in many cases you
+    may need to programmatically generate `ODESystem`s. If that's the case, check out
+    the [Programmatically Generating and Scripting ODESystems Tutorial](@ref programmatically).
   - Vector-valued parameters and variables are possible. A cleaner, more
     consistent treatment of these is still a work in progress, however. Once finished,
     this introductory tutorial will also cover this feature.
diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md
new file mode 100644
index 0000000000..5ec7425d39
--- /dev/null
+++ b/docs/src/tutorials/programmatically_generating.md
@@ -0,0 +1,106 @@
+# [Programmatically Generating and Scripting ODESystems](@id programmatically)
+
+In the following tutorial we will discuss how to programmatically generate `ODESystem`s.
+This is for cases where one is writing functions that generating `ODESystem`s, for example
+if implementing a reader which parses some file format to generate an `ODESystem` (for example,
+SBML), or for writing functions that transform an `ODESystem` (for example, if you write a
+function that log-transforms a variable in an `ODESystem`).
+
+## The Representation of a ModelingToolkit System
+
+ModelingToolkit is built on [Symbolics.jl](https://symbolics.juliasymbolics.org/dev/),
+a symbolic Computer Algebra System (CAS) developed in Julia. As such, all CAS functionality
+is available on ModelingToolkit systems, such as symbolic differentiation, Groebner basis
+calculations, and whatever else you can think of. Under the hood, all ModelingToolkit
+variables and expressions are Symbolics.jl variables and expressions. Thus when scripting
+a ModelingToolkit system, one simply needs to generate Symbolics.jl variables and equations
+as demonstrated in the Symbolics.jl documentation. This looks like:
+
+```@example scripting
+using Symbolics
+@variables t x(t) y(t) # Define variables
+D = Differential(t) # Define a differential operator
+eqs = [D(x) ~ y
+       D(y) ~ x] # Define an array of equations
+```
+
+## The Non-DSL (non-`@mtkmodel`) Way of Defining an ODESystem
+
+Using `@mtkmodel` is the preferred way of defining ODEs with MTK. However, let us
+look at how we can define the same system without `@mtkmodel`. This is useful for
+defining PDESystem etc.
+
+```@example scripting
+using ModelingToolkit
+
+@variables t x(t)   # independent and dependent variables
+@parameters τ       # parameters
+@constants h = 1    # constants
+D = Differential(t) # define an operator for the differentiation w.r.t. time
+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)
+
+# Perform the standard transformations and mark the model complete
+# Note: Complete models cannot be subsystems of other models!
+fol_model = complete(structural_simplify(fol_model))
+```
+
+As you can see, generating an ODESystem is as simple as creating the array of equations
+and passing it to the `ODESystem` constructor.
+
+## Understanding the Difference Between the Julia Variable and the Symbolic Variable
+
+In the most basic usage of ModelingToolkit and Symbolics, the name of the Julia variable
+and the symbolic variable are the same. For example, when we do:
+
+```@example scripting
+@variables a
+```
+
+the name of the symbolic variable is `a` and same with the Julia variable. However, we can
+de-couple these by setting `a` to a new symbolic variable, for example:
+
+```@example scripting
+b = only(@variables(a))
+```
+
+Now the Julia variable `b` refers to the variable named `a`. However, the downside of this current
+approach is that it requires that the user writing the script knows the name `a` that they want to
+place to the variable. But what if for example we needed to get the variable's name from a file?
+
+To do this, one can interpolate a symbol into the `@variables` macro using `$`. For example:
+
+```@example scripting
+a = :c
+b = only(@variables($a))
+```
+
+In this example, `@variables($a)` created a variable named `c`, and set this variable to `b`. 
+
+Variables are not the only thing with names. For example, when you build a system, it knows its name
+that name is used in the namespacing. In the standard usage, again the Julia variable and the
+symbolic name are made the same via:
+
+```@example scripting
+@named fol_model = ODESystem(eqs, t)
+```
+
+However, one can decouple these two properties by noting that `@named` is simply shorthand for the
+following:
+
+```@example scripting
+fol_model = ODESystem(eqs, t; name = :fol_model)
+```
+
+Thus if we had read a name from a file and wish to populate an `ODESystem` with said name, we could do:
+
+```@example scripting
+namesym = :name_from_file
+fol_model = ODESystem(eqs, t; name = namesym)
+```
+
+## Warning About Mutation
+
+Be advsied that it's never a good idea to mutate an `ODESystem`, or any `AbstractSystem`.
\ No newline at end of file
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 1a6e0356a5..de63a67837 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -201,7 +201,7 @@ export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, Flow, Stream, instream
-export @component, @mtkmodel
+export @component, @mtkmodel, @mtkbuild
 export isinput, isoutput, getbounds, hasbounds, isdisturbance, istunable, getdist, hasdist,
     tunable_parameters, isirreducible, getdescription, hasdescription, isbinaryvar,
     isintegervar
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index ab68cc9e7e..9c278c9804 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1180,6 +1180,15 @@ macro component(expr)
     esc(component_post_processing(expr, false))
 end
 
+macro mtkbuild(expr)
+    named_expr = ModelingToolkit.named_expr(expr)
+    name = named_expr.args[1]
+    esc(quote
+        $named_expr
+        $name = complete(structural_simplify($name))
+    end)
+end
+
 """
 $(SIGNATURES)
 

From ac75e770e8478fedd09de17c1290f782bca31bbe Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Sun, 8 Oct 2023 14:42:10 -0400
Subject: [PATCH 29/50] Add `parent` to hierarchical systems and set it in
 `structural_simplify`

---
 src/systems/abstractsystem.jl                  | 8 ++++++--
 src/systems/diffeqs/odesystem.jl               | 8 ++++++--
 src/systems/diffeqs/sdesystem.jl               | 8 ++++++--
 src/systems/discrete_system/discrete_system.jl | 8 ++++++--
 src/systems/nonlinear/nonlinearsystem.jl       | 9 +++++++--
 src/systems/optimization/optimizationsystem.jl | 9 +++++++--
 src/systems/systems.jl                         | 6 ++++++
 7 files changed, 44 insertions(+), 12 deletions(-)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 9c278c9804..37f4a49e23 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -230,7 +230,8 @@ for prop in [:eqs
     :gui_metadata
     :discrete_subsystems
     :unknown_states
-    :split_idxs]
+    :split_idxs
+    :parent]
     fname1 = Symbol(:get_, prop)
     fname2 = Symbol(:has_, prop)
     @eval begin
@@ -307,6 +308,9 @@ function Base.propertynames(sys::AbstractSystem; private = false)
 end
 
 function Base.getproperty(sys::AbstractSystem, name::Symbol; namespace = !iscomplete(sys))
+    if has_parent(sys) && (parent = get_parent(sys); parent !== nothing)
+        sys = parent
+    end
     wrap(getvar(sys, name; namespace = namespace))
 end
 function getvar(sys::AbstractSystem, name::Symbol; namespace = !iscomplete(sys))
@@ -1185,7 +1189,7 @@ macro mtkbuild(expr)
     name = named_expr.args[1]
     esc(quote
         $named_expr
-        $name = complete(structural_simplify($name))
+        $name = $structural_simplify($name)
     end)
 end
 
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
index 3b828702c2..9259953dfc 100644
--- a/src/systems/diffeqs/odesystem.jl
+++ b/src/systems/diffeqs/odesystem.jl
@@ -143,6 +143,10 @@ struct ODESystem <: AbstractODESystem
     split_idxs: a vector of vectors of indices for the split parameters.
     """
     split_idxs::Union{Nothing, Vector{Vector{Int}}}
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad,
         jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults,
@@ -151,7 +155,7 @@ struct ODESystem <: AbstractODESystem
         tearing_state = nothing,
         substitutions = nothing, complete = false,
         discrete_subsystems = nothing, unknown_states = nothing,
-        split_idxs = nothing; checks::Union{Bool, Int} = true)
+        split_idxs = nothing, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
             check_parameters(ps, iv)
@@ -165,7 +169,7 @@ struct ODESystem <: AbstractODESystem
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,
             connector_type, preface, cevents, devents, metadata, gui_metadata,
             tearing_state, substitutions, complete, discrete_subsystems,
-            unknown_states, split_idxs)
+            unknown_states, split_idxs, parent)
     end
 end
 
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
index 30d0e10666..8c42f7ea0d 100644
--- a/src/systems/diffeqs/sdesystem.jl
+++ b/src/systems/diffeqs/sdesystem.jl
@@ -115,13 +115,17 @@ struct SDESystem <: AbstractODESystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function SDESystem(tag, deqs, neqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed,
         tgrad,
         jac,
         ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connector_type,
         cevents, devents, metadata = nothing, gui_metadata = nothing,
-        complete = false;
+        complete = false, parent = nothing;
         checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
@@ -135,7 +139,7 @@ struct SDESystem <: AbstractODESystem
         new(tag, deqs, neqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac,
             Wfact, Wfact_t, name, systems, defaults, connector_type, cevents, devents,
-            metadata, gui_metadata, complete)
+            metadata, gui_metadata, complete, parent)
     end
 end
 
diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl
index c0307c586b..ecd0cae8c6 100644
--- a/src/systems/discrete_system/discrete_system.jl
+++ b/src/systems/discrete_system/discrete_system.jl
@@ -86,6 +86,10 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function DiscreteSystem(tag, discreteEqs, iv, dvs, ps, tspan, var_to_name, ctrls,
         observed,
@@ -93,7 +97,7 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
         systems, defaults, preface, connector_type,
         metadata = nothing, gui_metadata = nothing,
         tearing_state = nothing, substitutions = nothing,
-        complete = false; checks::Union{Bool, Int} = true)
+        complete = false, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckComponents) > 0
             check_variables(dvs, iv)
             check_parameters(ps, iv)
@@ -105,7 +109,7 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
             systems,
             defaults,
             preface, connector_type, metadata, gui_metadata,
-            tearing_state, substitutions, complete)
+            tearing_state, substitutions, complete, parent)
     end
 end
 
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index d33a198710..192089b7da 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -75,18 +75,23 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function NonlinearSystem(tag, eqs, states, ps, var_to_name, observed, jac, name,
         systems,
         defaults, connector_type, metadata = nothing,
         gui_metadata = nothing,
         tearing_state = nothing, substitutions = nothing,
-        complete = false; checks::Union{Bool, Int} = true)
+        complete = false, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
             all_dimensionless([states; ps]) || check_units(eqs)
         end
         new(tag, eqs, states, ps, var_to_name, observed, jac, name, systems, defaults,
-            connector_type, metadata, gui_metadata, tearing_state, substitutions, complete)
+            connector_type, metadata, gui_metadata, tearing_state, substitutions, complete,
+            parent)
     end
 end
 
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
index 2d935aa030..7a9b10cd49 100644
--- a/src/systems/optimization/optimizationsystem.jl
+++ b/src/systems/optimization/optimizationsystem.jl
@@ -56,10 +56,14 @@ struct OptimizationSystem <: AbstractOptimizationSystem
     complete: if a model `sys` is complete, then `sys.x` no longer performs namespacing.
     """
     complete::Bool
+    """
+    parent: the hierarchical parent system before simplification.
+    """
+    parent::Any
 
     function OptimizationSystem(tag, op, states, ps, var_to_name, observed,
         constraints, name, systems, defaults, metadata = nothing,
-        gui_metadata = nothing, complete = false;
+        gui_metadata = nothing, complete = false, parent = nothing;
         checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
             unwrap(op) isa Symbolic && check_units(op)
@@ -67,7 +71,8 @@ struct OptimizationSystem <: AbstractOptimizationSystem
             all_dimensionless([states; ps]) || check_units(constraints)
         end
         new(tag, op, states, ps, var_to_name, observed,
-            constraints, name, systems, defaults, metadata, gui_metadata, complete)
+            constraints, name, systems, defaults, metadata, gui_metadata, complete,
+            parent)
     end
 end
 
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 7aaaabb917..aabdb45ce0 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -17,6 +17,12 @@ This will convert all `inputs` to parameters and allow them to be unconnected, i
 simplification will allow models where `n_states = n_equations - n_inputs`.
 """
 function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
+    kwargs...)
+    newsys = __structural_simplify(sys, io; simplify, kwargs...)
+    @set! newsys.parent = complete(sys)
+    return complete(newsys)
+end
+function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
     kwargs...)
     sys = expand_connections(sys)
     sys isa DiscreteSystem && return sys

From b0c2938bbead2c43112c4d4bd58c5cf9a1e6afca Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Sun, 8 Oct 2023 14:43:29 -0400
Subject: [PATCH 30/50] Format

---
 docs/src/tutorials/ode_modeling.md                | 4 ++--
 docs/src/tutorials/programmatically_generating.md | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index b318a6c955..2782e6357f 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -86,7 +86,7 @@ end
 
 Note that equations in MTK use the tilde character (`~`) as equality sign.
 
-`@mtkbuild` creates an instance of `FOL` named as `fol`. 
+`@mtkbuild` creates an instance of `FOL` named as `fol`.
 
 After construction of the ODE, you can solve it using [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/):
 
@@ -174,7 +174,7 @@ sol[fol.x]
 or to get the second value in the time series for `x`:
 
 ```@example ode2
-sol[fol.x,2]
+sol[fol.x, 2]
 ```
 
 Similarly, the time series for `RHS` can be retrieved using the same indexing:
diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md
index 5ec7425d39..ca11243c84 100644
--- a/docs/src/tutorials/programmatically_generating.md
+++ b/docs/src/tutorials/programmatically_generating.md
@@ -21,7 +21,7 @@ using Symbolics
 @variables t x(t) y(t) # Define variables
 D = Differential(t) # Define a differential operator
 eqs = [D(x) ~ y
-       D(y) ~ x] # Define an array of equations
+    D(y) ~ x] # Define an array of equations
 ```
 
 ## The Non-DSL (non-`@mtkmodel`) Way of Defining an ODESystem
@@ -77,7 +77,7 @@ a = :c
 b = only(@variables($a))
 ```
 
-In this example, `@variables($a)` created a variable named `c`, and set this variable to `b`. 
+In this example, `@variables($a)` created a variable named `c`, and set this variable to `b`.
 
 Variables are not the only thing with names. For example, when you build a system, it knows its name
 that name is used in the namespacing. In the standard usage, again the Julia variable and the
@@ -103,4 +103,4 @@ fol_model = ODESystem(eqs, t; name = namesym)
 
 ## Warning About Mutation
 
-Be advsied that it's never a good idea to mutate an `ODESystem`, or any `AbstractSystem`.
\ No newline at end of file
+Be advsied that it's never a good idea to mutate an `ODESystem`, or any `AbstractSystem`.

From c3c030e1fb39924516a265c97e302bfdbe4eb97e Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Sun, 8 Oct 2023 14:50:04 -0400
Subject: [PATCH 31/50] Add "docs/src/assets/*.toml" to gitignore

---
 .gitignore | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitignore b/.gitignore
index 3401a5a4b3..6d305ad348 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,5 @@
 Manifest.toml
 .vscode
 .vscode/*
+docs/src/assets/Project.toml
+docs/src/assets/Manifest.toml

From cd0f04afe895026d398d0d6350f1a1647c1cf0b0 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Sun, 8 Oct 2023 14:55:56 -0400
Subject: [PATCH 32/50] Add tests

---
 test/model_parsing.jl | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 3f380df100..c64bfa7b76 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -1,6 +1,6 @@
 using ModelingToolkit, Test
 using ModelingToolkit: get_gui_metadata,
-    VariableDescription, getdefault, RegularConnector, get_ps
+    VariableDescription, getdefault, RegularConnector, get_ps, getname
 using URIs: URI
 using Distributions
 using Unitful
@@ -146,7 +146,11 @@ l15 0" stroke="black" stroke-width="1" stroke-linejoin="bevel" fill="none"></pat
     C_val = 20
     R_val = 20
     res__R = 100
-    @named rc = RC(; C_val, R_val, resistor.R = res__R)
+    @mtkbuild rc = RC(; C_val, R_val, resistor.R = res__R)
+    resistor = getproperty(rc, :resistor; namespace = false)
+    @test getname(rc.resistor) === getname(resistor)
+    @test getname(rc.resistor.R) === getname(resistor.R)
+    @test getname(rc.resistor.v) === getname(resistor.v)
     # Test that `resistor.R` overrides `R_val` in the argument.
     @test getdefault(rc.resistor.R) == res__R != R_val
     # Test that `C_val` passed via argument is set as default of C.
@@ -166,7 +170,7 @@ l15 0" stroke="black" stroke-width="1" stroke-linejoin="bevel" fill="none"></pat
     @test ModelingToolkit.get_gui_metadata(rc.resistor.p).layout == Pin.structure[:icon] ==
           URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "pin.png"))
 
-    @test length(equations(structural_simplify(rc))) == 1
+    @test length(equations(rc)) == 1
 end
 
 @testset "Parameters and Structural parameters in various modes" begin

From 0d5ab59421ee573a2ff713a58371ecf9e7442b46 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Sun, 8 Oct 2023 15:35:34 -0400
Subject: [PATCH 33/50] Fix CI

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

diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index aabdb45ce0..bdf8f31808 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -18,9 +18,20 @@ simplification will allow models where `n_states = n_equations - n_inputs`.
 """
 function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
     kwargs...)
-    newsys = __structural_simplify(sys, io; simplify, kwargs...)
+    newsys′ = __structural_simplify(sys, io; simplify, kwargs...)
+    if newsys′ isa Tuple
+        @assert length(newsys′) == 2
+        newsys = newsys′[1]
+    else
+        newsys = newsys′
+    end
     @set! newsys.parent = complete(sys)
-    return complete(newsys)
+    newsys = complete(newsys)
+    if newsys′ isa Tuple
+        return newsys, newsys′[2]
+    else
+        return newsys
+    end
 end
 function __structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
     kwargs...)

From 7370a45e5076f7db463c6f89ba9d4e12909e37a6 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 8 Oct 2023 21:49:08 +0200
Subject: [PATCH 34/50] fix tutorial syntax

---
 docs/src/tutorials/acausal_components.md | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index 4a7e4ddd1c..e6d90ac5e5 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -102,7 +102,7 @@ end
 
 @mtkbuild rc_model = RCModel(resistor.R = 2.0)
 u0 = [
-    rc_model.capacitor₊v => 0.0,
+    rc_model.capacitor.v => 0.0,
 ]
 prob = ODAEProblem(rc_model, u0, (0, 10.0))
 sol = solve(prob, Tsit5())
@@ -315,8 +315,8 @@ DAE solver](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/#OrdinaryD
 This is done as follows:
 
 ```@example acausal
-u0 = [rc_model.capacitor₊v => 0.0
-    rc_model.capacitor₊p₊i => 0.0]
+u0 = [rc_model.capacitor.v => 0.0
+    rc_model.capacitor.p.i => 0.0]
 
 prob = ODEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
@@ -329,7 +329,7 @@ letter `A`):
 
 ```@example acausal
 u0 = [
-    rc_model.capacitor₊v => 0.0,
+    rc_model.capacitor.v => 0.0,
 ]
 prob = ODAEProblem(sys, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
@@ -356,11 +356,11 @@ The solution object can be accessed via its symbols. For example, let's retrieve
 the voltage of the resistor over time:
 
 ```@example acausal
-sol[rc_model.resistor₊v]
+sol[rc_model.resistor.v]
 ```
 
 or we can plot the timeseries of the resistor's voltage:
 
 ```@example acausal
-plot(sol, idxs = [rc_model.resistor₊v])
+plot(sol, idxs = [rc_model.resistor.v])
 ```

From 137d3d116add2d9d5304a930d564d195fcea4f2f Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Sun, 8 Oct 2023 16:48:58 -0400
Subject: [PATCH 35/50] Support implicit name unpack in `at extend`

---
 src/systems/model_parsing.jl | 21 ++++++++++++++++-----
 test/model_parsing.jl        |  2 +-
 2 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index 6a95a6e9e8..16d36ca918 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -237,7 +237,7 @@ function parse_model!(exprs, comps, ext, eqs, icon, vs, ps, sps,
     if mname == Symbol("@components")
         parse_components!(exprs, comps, dict, body, kwargs)
     elseif mname == Symbol("@extend")
-        parse_extend!(exprs, ext, dict, body, kwargs)
+        parse_extend!(exprs, ext, dict, mod, body, kwargs)
     elseif mname == Symbol("@variables")
         parse_variables!(exprs, vs, dict, mod, body, :variables, kwargs)
     elseif mname == Symbol("@parameters")
@@ -372,7 +372,7 @@ function extend_args!(a, b, dict, expr, kwargs, varexpr, has_param = false)
     end
 end
 
-function parse_extend!(exprs, ext, dict, body, kwargs)
+function parse_extend!(exprs, ext, dict, mod, body, kwargs)
     expr = Expr(:block)
     varexpr = Expr(:block)
     push!(exprs, varexpr)
@@ -380,16 +380,27 @@ function parse_extend!(exprs, ext, dict, body, kwargs)
     body = deepcopy(body)
     MLStyle.@match body begin
         Expr(:(=), a, b) => begin
-            vars = nothing
             if Meta.isexpr(b, :(=))
                 vars = a
                 if !Meta.isexpr(vars, :tuple)
                     error("`@extend` destructuring only takes an tuple as LHS. Got $body")
                 end
                 a, b = b.args
-                extend_args!(a, b, dict, expr, kwargs, varexpr)
-                vars, a, b
+            else
+                if (model = getproperty(mod, b.args[1])) isa Model
+                    _vars = keys(get(model.structure, :variables, Dict()))
+                    _vars = union(_vars, keys(get(model.structure, :parameters, Dict())))
+                    _vars = union(_vars,
+                        map(first, get(model.structure, :components, Vector{Symbol}[])))
+                    vars = Expr(:tuple)
+                    append!(vars.args, collect(_vars))
+                else
+                    error("Cannot infer the exact `Model` that `@extend $(body)` refers." *
+                          " Please specify the names that it brings into scope by:" *
+                          " `@extend a, b = oneport = OnePort()`.")
+                end
             end
+            extend_args!(a, b, dict, expr, kwargs, varexpr)
             ext[] = a
             push!(b.args, Expr(:kw, :name, Meta.quot(a)))
             push!(expr.args, :($a = $b))
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 3f380df100..f678e30f16 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -101,7 +101,7 @@ l15 0" stroke="black" stroke-width="1" stroke-linejoin="bevel" fill="none"></pat
         @parameters begin
             C, [unit = u"F"]
         end
-        @extend v, i = oneport = OnePort(; v = 0.0)
+        @extend oneport = OnePort(; v = 0.0)
         @icon "https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg"
         @equations begin
             D(v) ~ i / C

From c679dd99374bd681918dafe7a4e75607b03a6c3b Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 8 Oct 2023 22:53:20 +0200
Subject: [PATCH 36/50] fix acausal tutorial

---
 Project.toml                             | 2 +-
 docs/src/tutorials/acausal_components.md | 6 +++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3834488ae3..ca73a90799 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "8.71.2"
+version = "8.72.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index e6d90ac5e5..cbee878ec6 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -318,7 +318,7 @@ This is done as follows:
 u0 = [rc_model.capacitor.v => 0.0
     rc_model.capacitor.p.i => 0.0]
 
-prob = ODEProblem(sys, u0, (0, 10.0))
+prob = ODEProblem(rc_model, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
 plot(sol)
 ```
@@ -331,7 +331,7 @@ letter `A`):
 u0 = [
     rc_model.capacitor.v => 0.0,
 ]
-prob = ODAEProblem(sys, u0, (0, 10.0))
+prob = ODAEProblem(rc_model, u0, (0, 10.0))
 sol = solve(prob, Rodas4())
 plot(sol)
 ```
@@ -344,7 +344,7 @@ like `structural_simplify` simply change state variables into observables which
 defined by `observed` equations.
 
 ```@example acausal
-observed(sys)
+observed(rc_model)
 ```
 
 These are explicit algebraic equations which can then be used to reconstruct

From 3d83c4ce35015febb7e3cb0574366f8da6f06948 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 8 Oct 2023 23:11:30 +0200
Subject: [PATCH 37/50] one more tutorial simplification

---
 docs/src/tutorials/ode_modeling.md | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index 2782e6357f..a68dcefd21 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -243,9 +243,8 @@ f_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t)) + 1]
     end
 end
 
-@named fol_external_f = FOLExternalFunction()
-fol_external_f = complete(fol_external_f)
-prob = ODEProblem(structural_simplify(fol_external_f),
+@mtkbuild fol_external_f = FOLExternalFunction()
+prob = ODEProblem(fol_external_f,
     [fol_external_f.x => 0.0],
     (0.0, 10.0),
     [fol_external_f.τ => 0.75])

From 73f9e0a2824b791b56b5b69bd2c3e9bb91016645 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Sun, 8 Oct 2023 18:16:45 -0400
Subject: [PATCH 38/50] Stop using `at testset` macro

---
 src/systems/model_parsing.jl |   2 +-
 test/model_parsing.jl        | 258 +++++++++++++++++------------------
 2 files changed, 129 insertions(+), 131 deletions(-)

diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index 16d36ca918..1d2bf4da83 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -386,7 +386,7 @@ function parse_extend!(exprs, ext, dict, mod, body, kwargs)
                     error("`@extend` destructuring only takes an tuple as LHS. Got $body")
                 end
                 a, b = b.args
-            else
+            elseif Meta.isexpr(b, :call)
                 if (model = getproperty(mod, b.args[1])) isa Model
                     _vars = keys(get(model.structure, :variables, Dict()))
                     _vars = union(_vars, keys(get(model.structure, :parameters, Dict())))
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 9608b19656..5ae797c5e4 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -7,77 +7,76 @@ using Unitful
 
 ENV["MTK_ICONS_DIR"] = "$(@__DIR__)/icons"
 
-@testset "Comprehensive Test of Parsing Models (with an RC Circuit)" begin
-    @connector RealInput begin
-        u(t), [input = true, unit = u"V"]
+@connector RealInput begin
+    u(t), [input = true, unit = u"V"]
+end
+@connector RealOutput begin
+    u(t), [output = true, unit = u"V"]
+end
+@mtkmodel Constant begin
+    @components begin
+        output = RealOutput()
     end
-    @connector RealOutput begin
-        u(t), [output = true, unit = u"V"]
+    @parameters begin
+        k, [description = "Constant output value of block"]
     end
-    @mtkmodel Constant begin
-        @components begin
-            output = RealOutput()
-        end
-        @parameters begin
-            k, [description = "Constant output value of block"]
-        end
-        @equations begin
-            output.u ~ k
-        end
+    @equations begin
+        output.u ~ k
     end
+end
 
-    @variables t [unit = u"s"]
-    D = Differential(t)
+@variables t [unit = u"s"]
+D = Differential(t)
 
-    @connector Pin begin
-        v(t), [unit = u"V"]                    # Potential at the pin [V]
-        i(t), [connect = Flow, unit = u"A"]    # Current flowing into the pin [A]
-        @icon "pin.png"
-    end
+@connector Pin begin
+    v(t), [unit = u"V"]                    # Potential at the pin [V]
+    i(t), [connect = Flow, unit = u"A"]    # Current flowing into the pin [A]
+    @icon "pin.png"
+end
 
-    @named p = Pin(; v = π)
-    @test getdefault(p.v) == π
-    @test Pin.isconnector == true
+@named p = Pin(; v = π)
+@test getdefault(p.v) == π
+@test Pin.isconnector == true
 
-    @mtkmodel OnePort begin
-        @components begin
-            p = Pin()
-            n = Pin()
-        end
-        @variables begin
-            v(t), [unit = u"V"]
-            i(t), [unit = u"A"]
-        end
-        @icon "oneport.png"
-        @equations begin
-            v ~ p.v - n.v
-            0 ~ p.i + n.i
-            i ~ p.i
-        end
+@mtkmodel OnePort begin
+    @components begin
+        p = Pin()
+        n = Pin()
+    end
+    @variables begin
+        v(t), [unit = u"V"]
+        i(t), [unit = u"A"]
     end
+    @icon "oneport.png"
+    @equations begin
+        v ~ p.v - n.v
+        0 ~ p.i + n.i
+        i ~ p.i
+    end
+end
 
-    @test OnePort.isconnector == false
+@test OnePort.isconnector == false
 
-    @mtkmodel Ground begin
-        @components begin
-            g = Pin()
-        end
-        @icon begin
-            read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
-        end
-        @equations begin
-            g.v ~ 0
-        end
+@mtkmodel Ground begin
+    @components begin
+        g = Pin()
+    end
+    @icon begin
+        read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
     end
+    @equations begin
+        g.v ~ 0
+    end
+end
 
-    resistor_log = "$(@__DIR__)/logo/resistor.svg"
-    @mtkmodel Resistor begin
-        @extend v, i = oneport = OnePort()
-        @parameters begin
-            R, [unit = u"Ω"]
-        end
-        @icon begin
-            """<?xml version="1.0" encoding="UTF-8"?>
+resistor_log = "$(@__DIR__)/logo/resistor.svg"
+@mtkmodel Resistor begin
+    @extend v, i = oneport = OnePort()
+    @parameters begin
+        R, [unit = u"Ω"]
+    end
+    @icon begin
+        """<?xml version="1.0" encoding="UTF-8"?>
 <svg xmlns="http://www.w3.org/2000/svg" width="80" height="30">
   <path d="M10 15
 l15 0
@@ -91,88 +90,87 @@ l2.5 -5
 l15 0" stroke="black" stroke-width="1" stroke-linejoin="bevel" fill="none"></path>
 </svg>
 """
-        end
-        @equations begin
-            v ~ i * R
-        end
     end
+    @equations begin
+        v ~ i * R
+    end
+end
 
-    @mtkmodel Capacitor begin
-        @parameters begin
-            C, [unit = u"F"]
-        end
-        @extend oneport = OnePort(; v = 0.0)
-        @icon "https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg"
-        @equations begin
-            D(v) ~ i / C
-        end
+@mtkmodel Capacitor begin
+    @parameters begin
+        C, [unit = u"F"]
+    end
+    @extend oneport = OnePort(; v = 0.0)
+    @icon "https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg"
+    @equations begin
+        D(v) ~ i / C
     end
+end
 
-    @named capacitor = Capacitor(C = 10, v = 10.0)
-    @test getdefault(capacitor.v) == 10.0
+@named capacitor = Capacitor(C = 10, v = 10.0)
+@test getdefault(capacitor.v) == 10.0
 
-    @mtkmodel Voltage begin
-        @extend v, i = oneport = OnePort()
-        @components begin
-            V = RealInput()
-        end
-        @equations begin
-            v ~ V.u
-        end
+@mtkmodel Voltage begin
+    @extend v, i = oneport = OnePort()
+    @components begin
+        V = RealInput()
     end
+    @equations begin
+        v ~ V.u
+    end
+end
 
-    @mtkmodel RC begin
-        @structural_parameters begin
-            R_val = 10
-            C_val = 10
-            k_val = 10
-        end
-        @components begin
-            resistor = Resistor(; R = R_val)
-            capacitor = Capacitor(; C = C_val)
-            source = Voltage()
-            constant = Constant(; k = k_val)
-            ground = Ground()
-        end
-
-        @equations begin
-            connect(constant.output, source.V)
-            connect(source.p, resistor.p)
-            connect(resistor.n, capacitor.p)
-            connect(capacitor.n, source.n, ground.g)
-        end
+@mtkmodel RC begin
+    @structural_parameters begin
+        R_val = 10
+        C_val = 10
+        k_val = 10
+    end
+    @components begin
+        resistor = Resistor(; R = R_val)
+        capacitor = Capacitor(; C = C_val)
+        source = Voltage()
+        constant = Constant(; k = k_val)
+        ground = Ground()
     end
 
-    C_val = 20
-    R_val = 20
-    res__R = 100
-    @mtkbuild rc = RC(; C_val, R_val, resistor.R = res__R)
-    resistor = getproperty(rc, :resistor; namespace = false)
-    @test getname(rc.resistor) === getname(resistor)
-    @test getname(rc.resistor.R) === getname(resistor.R)
-    @test getname(rc.resistor.v) === getname(resistor.v)
-    # Test that `resistor.R` overrides `R_val` in the argument.
-    @test getdefault(rc.resistor.R) == res__R != R_val
-    # Test that `C_val` passed via argument is set as default of C.
-    @test getdefault(rc.capacitor.C) == C_val
-    # Test that `k`'s default value is unchanged.
-    @test getdefault(rc.constant.k) == RC.structure[:kwargs][:k_val]
-    @test getdefault(rc.capacitor.v) == 0.0
-
-    @test get_gui_metadata(rc.resistor).layout == Resistor.structure[:icon] ==
-          read(joinpath(ENV["MTK_ICONS_DIR"], "resistor.svg"), String)
-    @test get_gui_metadata(rc.ground).layout ==
-          read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
-    @test get_gui_metadata(rc.capacitor).layout ==
-          URI("https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg")
-    @test OnePort.structure[:icon] ==
-          URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "oneport.png"))
-    @test ModelingToolkit.get_gui_metadata(rc.resistor.p).layout == Pin.structure[:icon] ==
-          URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "pin.png"))
-
-    @test length(equations(rc)) == 1
+    @equations begin
+        connect(constant.output, source.V)
+        connect(source.p, resistor.p)
+        connect(resistor.n, capacitor.p)
+        connect(capacitor.n, source.n, ground.g)
+    end
 end
 
+C_val = 20
+R_val = 20
+res__R = 100
+@mtkbuild rc = RC(; C_val, R_val, resistor.R = res__R)
+resistor = getproperty(rc, :resistor; namespace = false)
+@test getname(rc.resistor) === getname(resistor)
+@test getname(rc.resistor.R) === getname(resistor.R)
+@test getname(rc.resistor.v) === getname(resistor.v)
+# Test that `resistor.R` overrides `R_val` in the argument.
+@test getdefault(rc.resistor.R) == res__R != R_val
+# Test that `C_val` passed via argument is set as default of C.
+@test getdefault(rc.capacitor.C) == C_val
+# Test that `k`'s default value is unchanged.
+@test getdefault(rc.constant.k) == RC.structure[:kwargs][:k_val]
+@test getdefault(rc.capacitor.v) == 0.0
+
+@test get_gui_metadata(rc.resistor).layout == Resistor.structure[:icon] ==
+      read(joinpath(ENV["MTK_ICONS_DIR"], "resistor.svg"), String)
+@test get_gui_metadata(rc.ground).layout ==
+      read(abspath(ENV["MTK_ICONS_DIR"], "ground.svg"), String)
+@test get_gui_metadata(rc.capacitor).layout ==
+      URI("https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg")
+@test OnePort.structure[:icon] ==
+      URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "oneport.png"))
+@test ModelingToolkit.get_gui_metadata(rc.resistor.p).layout == Pin.structure[:icon] ==
+      URI("file:///" * abspath(ENV["MTK_ICONS_DIR"], "pin.png"))
+
+@test length(equations(rc)) == 1
+
 @testset "Parameters and Structural parameters in various modes" begin
     @mtkmodel MockModel begin
         @parameters begin

From 718ae818315fd3ffd1126aaf65961d7a77536434 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Mon, 9 Oct 2023 00:21:31 +0200
Subject: [PATCH 39/50] Remove ODAEProblem from acausal tutorial

---
 docs/src/tutorials/acausal_components.md | 21 +++------------------
 1 file changed, 3 insertions(+), 18 deletions(-)

diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index cbee878ec6..7639442162 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -104,8 +104,8 @@ end
 u0 = [
     rc_model.capacitor.v => 0.0,
 ]
-prob = ODAEProblem(rc_model, u0, (0, 10.0))
-sol = solve(prob, Tsit5())
+prob = ODEProblem(rc_model, u0, (0, 10.0))
+sol = solve(prob)
 plot(sol)
 ```
 
@@ -319,25 +319,10 @@ u0 = [rc_model.capacitor.v => 0.0
     rc_model.capacitor.p.i => 0.0]
 
 prob = ODEProblem(rc_model, u0, (0, 10.0))
-sol = solve(prob, Rodas4())
-plot(sol)
-```
-
-MTK can numerically solve all the
-unreduced algebraic equations using the `ODAEProblem` (note the
-letter `A`):
-
-```@example acausal
-u0 = [
-    rc_model.capacitor.v => 0.0,
-]
-prob = ODAEProblem(rc_model, u0, (0, 10.0))
-sol = solve(prob, Rodas4())
+sol = solve(prob)
 plot(sol)
 ```
 
-Notice that this solves the whole system by only solving for one variable!
-
 However, what if we wanted to plot the timeseries of a different variable? Do
 not worry, that information was not thrown away! Instead, transformations
 like `structural_simplify` simply change state variables into observables which are

From 298a9a5fde793e9a29828dc1f69546c9cfd99dd1 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Mon, 9 Oct 2023 05:32:37 +0200
Subject: [PATCH 40/50] patch system names

---
 Project.toml                       | 2 +-
 docs/src/tutorials/ode_modeling.md | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/Project.toml b/Project.toml
index ca73a90799..4dde81dde7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "8.72.0"
+version = "8.72.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index a68dcefd21..ae56d4f44d 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -154,12 +154,12 @@ along with the state variable. Note that this has to be requested explicitly
 like as follows:
 
 ```@example ode2
-prob = ODEProblem(fol_simplified,
-    [fol_simplified.x => 0.0],
+prob = ODEProblem(fol,
+    [fol.x => 0.0],
     (0.0, 10.0),
-    [fol_simplified.τ => 3.0])
+    [fol.τ => 3.0])
 sol = solve(prob)
-plot(sol, vars = [fol_simplified.x, fol_simplified.RHS])
+plot(sol, vars = [fol.x, fol.RHS])
 ```
 
 ## Named Indexing of Solutions

From bdc45d13a9b9176f2c5963e116fb270a8adecbb6 Mon Sep 17 00:00:00 2001
From: Yingbo Ma <mayingbo5@gmail.com>
Date: Mon, 9 Oct 2023 02:26:27 -0400
Subject: [PATCH 41/50] Simplify `at extend`

---
 src/systems/model_parsing.jl | 60 +++++++++++++++++++++++-------------
 test/model_parsing.jl        |  2 +-
 2 files changed, 39 insertions(+), 23 deletions(-)

diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
index 1d2bf4da83..299793c236 100644
--- a/src/systems/model_parsing.jl
+++ b/src/systems/model_parsing.jl
@@ -372,6 +372,28 @@ function extend_args!(a, b, dict, expr, kwargs, varexpr, has_param = false)
     end
 end
 
+const EMPTY_DICT = Dict()
+const EMPTY_VoVoSYMBOL = Vector{Symbol}[]
+
+function Base.names(model::Model)
+    vars = keys(get(model.structure, :variables, EMPTY_DICT))
+    vars = union(vars, keys(get(model.structure, :parameters, EMPTY_DICT)))
+    vars = union(vars,
+        map(first, get(model.structure, :components, EMPTY_VoVoSYMBOL)))
+    collect(vars)
+end
+
+function _parse_extend!(ext, a, b, dict, expr, kwargs, varexpr, vars)
+    extend_args!(a, b, dict, expr, kwargs, varexpr)
+    ext[] = a
+    push!(b.args, Expr(:kw, :name, Meta.quot(a)))
+    push!(expr.args, :($a = $b))
+
+    dict[:extend] = [Symbol.(vars.args), a, b.args[1]]
+
+    push!(expr.args, :(@unpack $vars = $a))
+end
+
 function parse_extend!(exprs, ext, dict, mod, body, kwargs)
     expr = Expr(:block)
     varexpr = Expr(:block)
@@ -386,33 +408,27 @@ function parse_extend!(exprs, ext, dict, mod, body, kwargs)
                     error("`@extend` destructuring only takes an tuple as LHS. Got $body")
                 end
                 a, b = b.args
-            elseif Meta.isexpr(b, :call)
-                if (model = getproperty(mod, b.args[1])) isa Model
-                    _vars = keys(get(model.structure, :variables, Dict()))
-                    _vars = union(_vars, keys(get(model.structure, :parameters, Dict())))
-                    _vars = union(_vars,
-                        map(first, get(model.structure, :components, Vector{Symbol}[])))
-                    vars = Expr(:tuple)
-                    append!(vars.args, collect(_vars))
-                else
-                    error("Cannot infer the exact `Model` that `@extend $(body)` refers." *
-                          " Please specify the names that it brings into scope by:" *
-                          " `@extend a, b = oneport = OnePort()`.")
-                end
+                _parse_extend!(ext, a, b, dict, expr, kwargs, varexpr, vars)
+            else
+                error("When explicitly destructing in `@extend` please use the syntax: `@extend a, b = oneport = OnePort()`.")
             end
-            extend_args!(a, b, dict, expr, kwargs, varexpr)
-            ext[] = a
-            push!(b.args, Expr(:kw, :name, Meta.quot(a)))
-            push!(expr.args, :($a = $b))
-
-            dict[:extend] = [Symbol.(vars.args), a, b.args[1]]
-
-            if vars !== nothing
-                push!(expr.args, :(@unpack $vars = $a))
+        end
+        Expr(:call, a′, _...) => begin
+            a = Symbol(Symbol("#mtkmodel"), :__anonymous__, a′)
+            b = body
+            if (model = getproperty(mod, b.args[1])) isa Model
+                vars = Expr(:tuple)
+                append!(vars.args, names(model))
+                _parse_extend!(ext, a, b, dict, expr, kwargs, varexpr, vars)
+            else
+                error("Cannot infer the exact `Model` that `@extend $(body)` refers." *
+                      " Please specify the names that it brings into scope by:" *
+                      " `@extend a, b = oneport = OnePort()`.")
             end
         end
         _ => error("`@extend` only takes an assignment expression. Got $body")
     end
+    return nothing
 end
 
 function parse_variable_arg!(expr, vs, dict, mod, arg, varclass, kwargs)
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
index 5ae797c5e4..701e785928 100644
--- a/test/model_parsing.jl
+++ b/test/model_parsing.jl
@@ -100,7 +100,7 @@ end
     @parameters begin
         C, [unit = u"F"]
     end
-    @extend oneport = OnePort(; v = 0.0)
+    @extend OnePort(; v = 0.0)
     @icon "https://upload.wikimedia.org/wikipedia/commons/7/78/Capacitor_symbol.svg"
     @equations begin
         D(v) ~ i / C

From 1152f30615a5493139f6c37556d8fd676d700611 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Mon, 9 Oct 2023 07:45:36 +0200
Subject: [PATCH 42/50] update docs for new extension syntax

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

diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index 7639442162..c3f5eed879 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -54,7 +54,7 @@ end
 end
 
 @mtkmodel Resistor begin
-    @extend v, i = oneport = OnePort()
+    @extend oneport = OnePort()
     @parameters begin
         R = 1.0 # Sets the default resistance
     end
@@ -76,7 +76,7 @@ D = Differential(t)
 end
 
 @mtkmodel ConstantVoltage begin
-    @extend (v,) = oneport = OnePort()
+    @extend oneport = OnePort()
     @parameters begin
         V = 1.0
     end
@@ -193,7 +193,7 @@ zero. This leads to our resistor equations:
 
 ```@example acausal
 @mtkmodel Resistor begin
-    @extend v, i = oneport = OnePort()
+    @extend oneport = OnePort()
     @parameters begin
         R = 1.0
     end
@@ -216,7 +216,7 @@ Using our knowledge of circuits, we similarly construct the `Capacitor`:
 D = Differential(t)
 
 @mtkmodel Capacitor begin
-    @extend v, i = oneport = OnePort()
+    @extend oneport = OnePort()
     @parameters begin
         C = 1.0
     end
@@ -233,7 +233,7 @@ model this as:
 
 ```@example acausal
 @mtkmodel ConstantVoltage begin
-    @extend (v,) = oneport = OnePort()
+    @extend oneport = OnePort()
     @parameters begin
         V = 1.0
     end

From 9e45fee9a2052bcf881e2099386206cbc42d852e Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Mon, 9 Oct 2023 08:40:07 +0200
Subject: [PATCH 43/50] fix for simplified extend

---
 docs/src/tutorials/acausal_components.md | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index c3f5eed879..508912c860 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -54,7 +54,7 @@ end
 end
 
 @mtkmodel Resistor begin
-    @extend oneport = OnePort()
+    @extend OnePort()
     @parameters begin
         R = 1.0 # Sets the default resistance
     end
@@ -66,7 +66,7 @@ end
 D = Differential(t)
 
 @mtkmodel Capacitor begin
-    @extend v, i = oneport = OnePort()
+    @extend OnePort()
     @parameters begin
         C = 1.0
     end
@@ -76,7 +76,7 @@ D = Differential(t)
 end
 
 @mtkmodel ConstantVoltage begin
-    @extend oneport = OnePort()
+    @extend OnePort()
     @parameters begin
         V = 1.0
     end
@@ -193,7 +193,7 @@ zero. This leads to our resistor equations:
 
 ```@example acausal
 @mtkmodel Resistor begin
-    @extend oneport = OnePort()
+    @extend OnePort()
     @parameters begin
         R = 1.0
     end
@@ -216,7 +216,7 @@ Using our knowledge of circuits, we similarly construct the `Capacitor`:
 D = Differential(t)
 
 @mtkmodel Capacitor begin
-    @extend oneport = OnePort()
+    @extend OnePort()
     @parameters begin
         C = 1.0
     end
@@ -233,7 +233,7 @@ model this as:
 
 ```@example acausal
 @mtkmodel ConstantVoltage begin
-    @extend oneport = OnePort()
+    @extend OnePort()
     @parameters begin
         V = 1.0
     end

From d4239a77eb6ed91ddf16aa381d870ab3c43d1802 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Mon, 9 Oct 2023 09:48:18 +0200
Subject: [PATCH 44/50] Update Project.toml

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

diff --git a/Project.toml b/Project.toml
index 4dde81dde7..4474aa35ad 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <mayingbo5@gmail.com>", "Chris Rackauckas <accounts@chrisrackauckas.com> and contributors"]
-version = "8.72.1"
+version = "8.72.2"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

From bc93796722531dab1ed1b0978f81c19bf7bb2f97 Mon Sep 17 00:00:00 2001
From: Arno Strouwen <arno.strouwen@telenet.be>
Date: Sat, 14 Oct 2023 11:00:30 +0200
Subject: [PATCH 45/50] non canonical docstrings

---
 docs/src/systems/DiscreteSystem.md  | 2 +-
 docs/src/systems/JumpSystem.md      | 8 +++++---
 docs/src/systems/NonlinearSystem.md | 6 +++---
 docs/src/systems/ODESystem.md       | 4 ++--
 docs/src/systems/SDESystem.md       | 8 +++++---
 5 files changed, 16 insertions(+), 12 deletions(-)

diff --git a/docs/src/systems/DiscreteSystem.md b/docs/src/systems/DiscreteSystem.md
index c9d4b81c48..d9934d5fb7 100644
--- a/docs/src/systems/DiscreteSystem.md
+++ b/docs/src/systems/DiscreteSystem.md
@@ -23,7 +23,7 @@ DiscreteSystem
 
 ```@docs
 DiscreteFunction(sys::ModelingToolkit.DiscreteSystem, args...)
-DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, args...)
+DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, u0map, tspan)
 ```
 
 ## Expression Constructors
diff --git a/docs/src/systems/JumpSystem.md b/docs/src/systems/JumpSystem.md
index 3dbeb8500d..68f0de2a5c 100644
--- a/docs/src/systems/JumpSystem.md
+++ b/docs/src/systems/JumpSystem.md
@@ -15,7 +15,7 @@ JumpSystem
 
 ## Transformations
 
-```@docs
+```@docs; canonical=false
 structural_simplify
 ```
 
@@ -23,7 +23,9 @@ structural_simplify
 
 ## Problem Constructors
 
+```@docs; canonical=false
+DiscreteProblem(sys::ModelingToolkit.DiscreteSystem, u0map, tspan)
+```
 ```@docs
-SciMLBase.DiscreteProblem(sys::JumpSystem,args...)
-JumpProcesses.JumpProblem(sys::JumpSystem,args...)
+JumpProblem(sys::JumpSystem, prob, aggregator)
 ```
diff --git a/docs/src/systems/NonlinearSystem.md b/docs/src/systems/NonlinearSystem.md
index d14e6a035c..dfdb8db508 100644
--- a/docs/src/systems/NonlinearSystem.md
+++ b/docs/src/systems/NonlinearSystem.md
@@ -15,7 +15,7 @@ NonlinearSystem
 
 ## Transformations
 
-```@docs
+```@docs; canonical=false
 structural_simplify
 alias_elimination
 tearing
@@ -23,14 +23,14 @@ tearing
 
 ## Analyses
 
-```@docs
+```@docs; canonical=false
 ModelingToolkit.isaffine
 ModelingToolkit.islinear
 ```
 
 ## Applicable Calculation and Generation Functions
 
-```julia
+```@docs; canonical=false
 calculate_jacobian
 generate_jacobian
 jacobian_sparsity
diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md
index 207ab098e7..f79263436f 100644
--- a/docs/src/systems/ODESystem.md
+++ b/docs/src/systems/ODESystem.md
@@ -35,7 +35,7 @@ ModelingToolkit.isaffine
 
 ## Applicable Calculation and Generation Functions
 
-```julia
+```@docs; canonical=false
 calculate_jacobian
 calculate_tgrad
 calculate_factorized_W
@@ -56,7 +56,7 @@ SteadyStateProblem(sys::ModelingToolkit.AbstractODESystem, args...)
 ## Torn Problem Constructors
 
 ```@docs
-ODAEProblem(sys::ModelingToolkit.AbstractODESystem, args...)
+ODAEProblem
 ```
 
 ## Expression Constructors
diff --git a/docs/src/systems/SDESystem.md b/docs/src/systems/SDESystem.md
index 57b89b304a..17f062c5ae 100644
--- a/docs/src/systems/SDESystem.md
+++ b/docs/src/systems/SDESystem.md
@@ -22,17 +22,19 @@ sde = SDESystem(ode, noiseeqs)
 
 ## Transformations
 
-```@docs
+```@docs; canonical=false
 structural_simplify
 alias_elimination
-Girsanov_transform
+```
+```@docs
+ModelingToolkit.Girsanov_transform
 ```
 
 ## Analyses
 
 ## Applicable Calculation and Generation Functions
 
-```julia
+```@docs; canonical=false
 calculate_jacobian
 calculate_tgrad
 calculate_factorized_W

From b54b05ba58b6eda239d7dd8c1d8161f296f89914 Mon Sep 17 00:00:00 2001
From: Arno Strouwen <arno.strouwen@telenet.be>
Date: Sat, 14 Oct 2023 11:00:49 +0200
Subject: [PATCH 46/50] update contributing docs

---
 docs/src/index.md | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index fc30b36714..1820bb40d6 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -173,13 +173,16 @@ system:
 
   - Please refer to the
     [SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://github.com/SciML/ColPrac/blob/master/README.md)
-    for guidance on PRs, issues, and other matters relating to contributing to ModelingToolkit.
+    for guidance on PRs, issues, and other matters relating to contributing to SciML.
 
+  - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions.
   - There are a few community forums:
     
-      + The #diffeq-bridged channel in the [Julia Slack](https://julialang.org/slack/)
-      + [JuliaDiffEq](https://gitter.im/JuliaDiffEq/Lobby) on Gitter
-      + On the Julia Discourse forums (look for the [modelingtoolkit tag](https://discourse.julialang.org/tag/modelingtoolkit))
+      + The #diffeq-bridged and #sciml-bridged channels in the
+        [Julia Slack](https://julialang.org/slack/)
+      + The #diffeq-bridged and #sciml-bridged channels in the
+        [Julia Zulip](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged)
+      + On the [Julia Discourse forums](https://discourse.julialang.org)
       + See also [SciML Community page](https://sciml.ai/community/)
 
 ## Reproducibility

From 2c86306eca02251e1fe6b43d45833762264d9080 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 14 Oct 2023 11:18:38 -0400
Subject: [PATCH 47/50] updates

---
 Project.toml                                  |  3 +-
 docs/Project.toml                             |  1 +
 docs/pages.jl                                 |  1 +
 .../bifurcation_diagram_computation.md        | 93 +++++++++++++++++++
 ext/MTKBifurcationKitExt.jl                   | 18 ++--
 test/extensions/bifurcationkit.jl             | 29 ++++++
 test/runtests.jl                              |  1 +
 7 files changed, 136 insertions(+), 10 deletions(-)
 create mode 100644 docs/src/tutorials/bifurcation_diagram_computation.md
 create mode 100644 test/extensions/bifurcationkit.jl

diff --git a/Project.toml b/Project.toml
index ae899f1220..86da45a0a7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -102,6 +102,7 @@ julia = "1.6"
 
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
@@ -125,4 +126,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
+test = ["AmplNLWriter", "BifurcationKit", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 81a99e004e..68d33cc4c5 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
diff --git a/docs/pages.jl b/docs/pages.jl
index 0b6d0555cd..0104c5eb48 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -7,6 +7,7 @@ pages = [
         "tutorials/modelingtoolkitize.md",
         "tutorials/stochastic_diffeq.md",
         "tutorials/parameter_identifiability.md",
+        "tutorials/bifurcation_diagram_computation.md",
         "tutorials/domain_connections.md"],
     "Examples" => Any["Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",
diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md
new file mode 100644
index 0000000000..1bb6e33935
--- /dev/null
+++ b/docs/src/tutorials/bifurcation_diagram_computation.md
@@ -0,0 +1,93 @@
+# [Bifurcation Diagrams](@id bifurcation_diagrams)
+Bifurcation diagrams describes how, for a dynamic system, the quantity and quality of its steady states changes with a parameter's value. These can be computed through the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. ModelingToolkit provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `NonlinearSystem`s and `ODESystem`s. All teh features provided by BifurcationKit can then be applied to these systems. This tutorial provides a brief introduction for these features, with BifurcationKit.jl providing [a more extensive documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/).
+
+### Creating a `BifurcationProblem`
+Let us first consider a simple `NonlinearSystem`:
+```@example Bif1
+using ModelingToolkit
+@variables t x(t) y(t)
+@parameters μ α
+eqs = [0 ~ μ*x - x^3 + α*y,
+        0 ~ -y]
+@named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+```
+we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information:
+1. The system for which we wish to compute the bifurcation diagram (`nsys`).
+2. The parameter which we wish to vary (`μ`).
+3. The parameter set for which we want to compute the bifurcation diagram.
+4. An initial guess of the state of the system for which there is a steady state at our provided parameter value.
+5. The variable which value we wish to plot in the bifurcation diagram (this argument is optional, if not provided, BifurcationKit default plot functions are used).
+
+We declare this additional information:
+```@example Bif1
+bif_par = μ
+p_start = [μ => -1.0, α => 1.0]
+u0_guess = [x => 1.0, y => 1.0]
+plot_var = x;
+```
+For the initial state guess (`u0_guess`), typically any value can be provided, however, read BifurcatioKit's documentation for more details.
+
+We can now create our `BifurcationProblem`, which can be provided as input to BifurcationKit's various functions.
+```@example Bif1
+using BifurcationKit
+bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+```
+Here, the `jac` argument (by default set to `true`) sets whenever to provide BifurcationKit with a Jacobian or not.
+
+
+### Computing a bifurcation diagram
+
+Let us consider the `BifurcationProblem` from the last section. If we wish to compute the corresponding bifurcation diagram we must first declare various settings used by BifurcationKit to compute the diagram. These are stored in a `ContinuationPar` structure (which also contain a `NewtonPar` structure). 
+```@example Bif1
+p_span = (-4.0, 6.0)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_min = p_span[1], p_max = p_span[2],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+```
+Here, `p_span` sets the interval over which we wish to compute the diagram.
+
+Next, we can use this as input to our bifurcation diagram, and then plot it.
+```@example Bif1
+bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+```
+Here, the value `2` sets how sub-branches of the diagram that BifurcationKit should compute. Generally, for bifurcation diagrams, it is recommended to use the `bothside=true` argument.
+```@example Bif1
+using Plots
+plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, xguide="μ", yguide = "x")
+```
+Here, the system exhibits a pitchfork bifurcation at *μ=0.0*.
+
+### Using `ODESystem` inputs
+It is also possible to use `ODESystem`s (rather than `NonlinearSystem`s) as input to `BifurcationProblem`. Here follows a brief such example.
+
+```@example Bif2
+using BifurcationKit, ModelingToolkit, Plots
+
+@variables t x(t) y(t)
+@parameters μ
+D = Differential(t)
+eqs = [D(x) ~ μ*x - y - x*(x^2+y^2), 
+        D(y) ~ x + μ*y - y*(x^2+y^2)]
+@named osys = ODESystem(eqs, t)
+
+bif_par = μ
+plot_var = x
+p_start = [μ => 1.0]
+u0_guess = [x => 0.0, y=> 0.0]
+
+bprob = BifurcationProblem(osys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+
+p_span = (-3.0, 3.0)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_max = p_span[2], p_min = p_span[1],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9)
+
+bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+using Plots
+plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, xguide="μ", yguide = "x")
+```
+Here, the value of `x` in the steady state does not change, however, at `μ=0` a Hopf bifurcation occur and the steady state turn unstable.
\ No newline at end of file
diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index 60ee91b9df..9a966ba55c 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -11,11 +11,11 @@ import BifurcationKit
 ### Creates BifurcationProblem Overloads ###
 
 # When input is a NonlinearSystem.
-function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_start, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, kwargs...)
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_bif, ps, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, jac=true, kwargs...)
     # Creates F and J functions.
-    ofun = NonlinearFunction(nsys; jac=true)
+    ofun = NonlinearFunction(nsys; jac=jac)
     F = ofun.f
-    J = ofun.jac
+    J = jac ? ofun.jac : nothing
 
     # Computes bifurcation parameter and plot var indexes.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))
@@ -25,16 +25,16 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_st
     end
 
     # Converts the input state guess.
-    u0_guess = ModelingToolkit.varmap_to_vars(u0_guess, states(nsys))
+    u0_bif = ModelingToolkit.varmap_to_vars(u0_bif, states(nsys))
+    ps = ModelingToolkit.varmap_to_vars(ps, parameters(nsys))
 
-    return BifurcationKit.BifurcationProblem(F, u0_guess, [p_start], (@lens _[1]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
+    return BifurcationKit.BifurcationProblem(F, u0_bif, ps, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
 end
 
 # When input is a ODESystem.
-function BifurcationKit.BifurcationProblem(osys::ODESystem, u0, p, args...; kwargs...)
-    return BifurcationProblem(convert(NonlinearProblem, osys), args...; kwargs...)
+function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
+    nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)], states(osys), parameters(osys); name=osys.name)
+    return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
 end
 
-
-
 end # module
diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl
new file mode 100644
index 0000000000..c9ae8890ad
--- /dev/null
+++ b/test/extensions/bifurcationkit.jl
@@ -0,0 +1,29 @@
+using BifurcationKit, ModelingToolkit, Test
+
+# Checks pitchfork diagram and that there are the correct number of branches (a main one and two children)
+let 
+    @variables t x(t) y(t)
+    @parameters μ α
+    eqs = [0 ~ μ*x - x^3 + α*y,
+            0 ~ -y]
+    @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+
+    bif_par = μ
+    p_start = [μ => -1.0, α => 1.0]
+    u0_guess = [x => 1.0, y => 1.0]
+    plot_var = x;
+
+    using BifurcationKit
+    bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+
+    p_span = (-4.0, 6.0)
+    opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+    opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+        max_steps = 100, nev = 2, newton_options = opt_newton,
+        p_min = p_span[1], p_max = p_span[2],
+        detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+
+    bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+
+    @test length(bf.child) == 2
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index e6387e2cd9..0322bba523 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -55,5 +55,6 @@ end
 @safetestset "OptimizationSystem Test" include("optimizationsystem.jl")
 @safetestset "FuncAffect Test" include("funcaffect.jl")
 @safetestset "Constants Test" include("constants.jl")
+@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl")
 # Reference tests go Last
 @safetestset "Latexify recipes Test" include("latexify.jl")

From 3249aebd18c3006699efb0e35d458b44260f9e9d Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 4 Oct 2023 11:57:52 -0400
Subject: [PATCH 48/50] init

---
 Project.toml                |  2 ++
 ext/MTKBifurcationKitExt.jl | 44 +++++++++++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+)
 create mode 100644 ext/MTKBifurcationKitExt.jl

diff --git a/Project.toml b/Project.toml
index 4474aa35ad..a578cc3974 100644
--- a/Project.toml
+++ b/Project.toml
@@ -51,9 +51,11 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [weakdeps]
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
 
 [extensions]
+MTKBifurcationKitExt = "BifurcationKit"
 MTKDeepDiffsExt = "DeepDiffs"
 
 [compat]
diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
new file mode 100644
index 0000000000..5aa37e7117
--- /dev/null
+++ b/ext/MTKBifurcationKitExt.jl
@@ -0,0 +1,44 @@
+module MTKBifurcationKitExt
+
+println("BifurcationKit extension loaded")
+
+### Preparations ###
+
+# Imports
+using ModelingToolkit, Setfield
+import BifurcationKit: BifurcationProblem
+
+### Creates BifurcationProblem Overloads ###
+
+# When input is a NonlinearProblem.
+function BifurcationKit.BifurcationProblem(nprob::NonlinearProblem, u0, p, bif_par, args...; plot_variable=nothing, record_from_solution=record_sol_default, kwargs...)
+    # check if jac does nto exist and modify.
+    F = nprob.f.f
+    J = nprob.f.jac.f
+
+    bif_idx = findfirst(isequal(bif_par), parameters(nsys))
+    if isnothing(plot_variable) 
+        plot_idx = findfirst(isequal(plot_variable), variables(nsys))
+        record_from_solution = (x, p) -> x[plot_idx]
+    end
+    return BifurcationProblem(F, u0, p, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
+end
+
+# When input is a NonlinearSystem.
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0, p, args...; kwargs...)
+    return BifurcationProblem(NonlinearProblem(nsys, u0, p; jac=true), args...; kwargs...)
+end
+
+# When input is a ODEProblem.
+function BifurcationKit.BifurcationProblem(oprob::ODEProblem, u0, p, args...; kwargs...)
+    return BifurcationProblem(convert(NonlinearProblem, oprob), args...; kwargs...)
+end
+
+# When input is a ODESystem.
+function BifurcationKit.BifurcationProblem(osys::ODESystem, u0, p, args...; kwargs...)
+    return BifurcationProblem(convert(NonlinearProblem, osys), args...; kwargs...)
+end
+
+
+
+end # module

From 7f5d8beb8dcb7a0d461f9dcb9b32163702695865 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 4 Oct 2023 12:44:11 -0400
Subject: [PATCH 49/50] makefunction

---
 ext/MTKBifurcationKitExt.jl | 30 +++++++++++++-----------------
 1 file changed, 13 insertions(+), 17 deletions(-)

diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index 5aa37e7117..60ee91b9df 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -6,32 +6,28 @@ println("BifurcationKit extension loaded")
 
 # Imports
 using ModelingToolkit, Setfield
-import BifurcationKit: BifurcationProblem
+import BifurcationKit
 
 ### Creates BifurcationProblem Overloads ###
 
-# When input is a NonlinearProblem.
-function BifurcationKit.BifurcationProblem(nprob::NonlinearProblem, u0, p, bif_par, args...; plot_variable=nothing, record_from_solution=record_sol_default, kwargs...)
-    # check if jac does nto exist and modify.
-    F = nprob.f.f
-    J = nprob.f.jac.f
+# When input is a NonlinearSystem.
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_start, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, kwargs...)
+    # Creates F and J functions.
+    ofun = NonlinearFunction(nsys; jac=true)
+    F = ofun.f
+    J = ofun.jac
 
+    # Computes bifurcation parameter and plot var indexes.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))
-    if isnothing(plot_variable) 
-        plot_idx = findfirst(isequal(plot_variable), variables(nsys))
+    if !isnothing(plot_var) 
+        plot_idx = findfirst(isequal(plot_var), states(nsys))
         record_from_solution = (x, p) -> x[plot_idx]
     end
-    return BifurcationProblem(F, u0, p, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
-end
 
-# When input is a NonlinearSystem.
-function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0, p, args...; kwargs...)
-    return BifurcationProblem(NonlinearProblem(nsys, u0, p; jac=true), args...; kwargs...)
-end
+    # Converts the input state guess.
+    u0_guess = ModelingToolkit.varmap_to_vars(u0_guess, states(nsys))
 
-# When input is a ODEProblem.
-function BifurcationKit.BifurcationProblem(oprob::ODEProblem, u0, p, args...; kwargs...)
-    return BifurcationProblem(convert(NonlinearProblem, oprob), args...; kwargs...)
+    return BifurcationKit.BifurcationProblem(F, u0_guess, [p_start], (@lens _[1]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
 end
 
 # When input is a ODESystem.

From c4a1a20422ca8d1f25cf75ecdba5ab491dbec621 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 14 Oct 2023 11:18:38 -0400
Subject: [PATCH 50/50] updates

---
 Project.toml                                  |  3 +-
 docs/Project.toml                             |  1 +
 docs/pages.jl                                 |  1 +
 .../bifurcation_diagram_computation.md        | 93 +++++++++++++++++++
 ext/MTKBifurcationKitExt.jl                   | 18 ++--
 test/extensions/bifurcationkit.jl             | 29 ++++++
 test/runtests.jl                              |  1 +
 7 files changed, 136 insertions(+), 10 deletions(-)
 create mode 100644 docs/src/tutorials/bifurcation_diagram_computation.md
 create mode 100644 test/extensions/bifurcationkit.jl

diff --git a/Project.toml b/Project.toml
index a578cc3974..8cceafe504 100644
--- a/Project.toml
+++ b/Project.toml
@@ -102,6 +102,7 @@ julia = "1.6"
 
 [extras]
 AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ControlSystemsMTK = "687d7614-c7e5-45fc-bfc3-9ee385575c88"
 DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6"
@@ -125,4 +126,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
+test = ["AmplNLWriter", "BifurcationKit", "BenchmarkTools", "ControlSystemsMTK", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 81a99e004e..68d33cc4c5 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
diff --git a/docs/pages.jl b/docs/pages.jl
index f6f5d7af50..9d49c477ec 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -8,6 +8,7 @@ pages = [
         "tutorials/programmatically_generating.md",
         "tutorials/stochastic_diffeq.md",
         "tutorials/parameter_identifiability.md",
+        "tutorials/bifurcation_diagram_computation.md",
         "tutorials/domain_connections.md"],
     "Examples" => Any["Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",
diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md
new file mode 100644
index 0000000000..1bb6e33935
--- /dev/null
+++ b/docs/src/tutorials/bifurcation_diagram_computation.md
@@ -0,0 +1,93 @@
+# [Bifurcation Diagrams](@id bifurcation_diagrams)
+Bifurcation diagrams describes how, for a dynamic system, the quantity and quality of its steady states changes with a parameter's value. These can be computed through the [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) package. ModelingToolkit provides a simple interface for creating BifurcationKit compatible `BifurcationProblem`s from `NonlinearSystem`s and `ODESystem`s. All teh features provided by BifurcationKit can then be applied to these systems. This tutorial provides a brief introduction for these features, with BifurcationKit.jl providing [a more extensive documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/).
+
+### Creating a `BifurcationProblem`
+Let us first consider a simple `NonlinearSystem`:
+```@example Bif1
+using ModelingToolkit
+@variables t x(t) y(t)
+@parameters μ α
+eqs = [0 ~ μ*x - x^3 + α*y,
+        0 ~ -y]
+@named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+```
+we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information:
+1. The system for which we wish to compute the bifurcation diagram (`nsys`).
+2. The parameter which we wish to vary (`μ`).
+3. The parameter set for which we want to compute the bifurcation diagram.
+4. An initial guess of the state of the system for which there is a steady state at our provided parameter value.
+5. The variable which value we wish to plot in the bifurcation diagram (this argument is optional, if not provided, BifurcationKit default plot functions are used).
+
+We declare this additional information:
+```@example Bif1
+bif_par = μ
+p_start = [μ => -1.0, α => 1.0]
+u0_guess = [x => 1.0, y => 1.0]
+plot_var = x;
+```
+For the initial state guess (`u0_guess`), typically any value can be provided, however, read BifurcatioKit's documentation for more details.
+
+We can now create our `BifurcationProblem`, which can be provided as input to BifurcationKit's various functions.
+```@example Bif1
+using BifurcationKit
+bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+```
+Here, the `jac` argument (by default set to `true`) sets whenever to provide BifurcationKit with a Jacobian or not.
+
+
+### Computing a bifurcation diagram
+
+Let us consider the `BifurcationProblem` from the last section. If we wish to compute the corresponding bifurcation diagram we must first declare various settings used by BifurcationKit to compute the diagram. These are stored in a `ContinuationPar` structure (which also contain a `NewtonPar` structure). 
+```@example Bif1
+p_span = (-4.0, 6.0)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_min = p_span[1], p_max = p_span[2],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+```
+Here, `p_span` sets the interval over which we wish to compute the diagram.
+
+Next, we can use this as input to our bifurcation diagram, and then plot it.
+```@example Bif1
+bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+```
+Here, the value `2` sets how sub-branches of the diagram that BifurcationKit should compute. Generally, for bifurcation diagrams, it is recommended to use the `bothside=true` argument.
+```@example Bif1
+using Plots
+plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, xguide="μ", yguide = "x")
+```
+Here, the system exhibits a pitchfork bifurcation at *μ=0.0*.
+
+### Using `ODESystem` inputs
+It is also possible to use `ODESystem`s (rather than `NonlinearSystem`s) as input to `BifurcationProblem`. Here follows a brief such example.
+
+```@example Bif2
+using BifurcationKit, ModelingToolkit, Plots
+
+@variables t x(t) y(t)
+@parameters μ
+D = Differential(t)
+eqs = [D(x) ~ μ*x - y - x*(x^2+y^2), 
+        D(y) ~ x + μ*y - y*(x^2+y^2)]
+@named osys = ODESystem(eqs, t)
+
+bif_par = μ
+plot_var = x
+p_start = [μ => 1.0]
+u0_guess = [x => 0.0, y=> 0.0]
+
+bprob = BifurcationProblem(osys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+
+p_span = (-3.0, 3.0)
+opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+	max_steps = 100, nev = 2, newton_options = opt_newton,
+	p_max = p_span[2], p_min = p_span[1],
+	detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9)
+
+bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+using Plots
+plot(bf; putspecialptlegend=false, markersize=2, plotfold=false, xguide="μ", yguide = "x")
+```
+Here, the value of `x` in the steady state does not change, however, at `μ=0` a Hopf bifurcation occur and the steady state turn unstable.
\ No newline at end of file
diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index 60ee91b9df..9a966ba55c 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -11,11 +11,11 @@ import BifurcationKit
 ### Creates BifurcationProblem Overloads ###
 
 # When input is a NonlinearSystem.
-function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_start, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, kwargs...)
+function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_bif, ps, bif_par, args...; plot_var=nothing, record_from_solution=BifurcationKit.record_sol_default, jac=true, kwargs...)
     # Creates F and J functions.
-    ofun = NonlinearFunction(nsys; jac=true)
+    ofun = NonlinearFunction(nsys; jac=jac)
     F = ofun.f
-    J = ofun.jac
+    J = jac ? ofun.jac : nothing
 
     # Computes bifurcation parameter and plot var indexes.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))
@@ -25,16 +25,16 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem, u0_guess, p_st
     end
 
     # Converts the input state guess.
-    u0_guess = ModelingToolkit.varmap_to_vars(u0_guess, states(nsys))
+    u0_bif = ModelingToolkit.varmap_to_vars(u0_bif, states(nsys))
+    ps = ModelingToolkit.varmap_to_vars(ps, parameters(nsys))
 
-    return BifurcationKit.BifurcationProblem(F, u0_guess, [p_start], (@lens _[1]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
+    return BifurcationKit.BifurcationProblem(F, u0_bif, ps, (@lens _[bif_idx]), args...; record_from_solution = record_from_solution, J = J, kwargs...)
 end
 
 # When input is a ODESystem.
-function BifurcationKit.BifurcationProblem(osys::ODESystem, u0, p, args...; kwargs...)
-    return BifurcationProblem(convert(NonlinearProblem, osys), args...; kwargs...)
+function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
+    nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)], states(osys), parameters(osys); name=osys.name)
+    return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
 end
 
-
-
 end # module
diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl
new file mode 100644
index 0000000000..c9ae8890ad
--- /dev/null
+++ b/test/extensions/bifurcationkit.jl
@@ -0,0 +1,29 @@
+using BifurcationKit, ModelingToolkit, Test
+
+# Checks pitchfork diagram and that there are the correct number of branches (a main one and two children)
+let 
+    @variables t x(t) y(t)
+    @parameters μ α
+    eqs = [0 ~ μ*x - x^3 + α*y,
+            0 ~ -y]
+    @named nsys = NonlinearSystem(eqs, [x, y], [μ, α])
+
+    bif_par = μ
+    p_start = [μ => -1.0, α => 1.0]
+    u0_guess = [x => 1.0, y => 1.0]
+    plot_var = x;
+
+    using BifurcationKit
+    bprob = BifurcationProblem(nsys, u0_guess, p_start, bif_par; plot_var=plot_var, jac=false)
+
+    p_span = (-4.0, 6.0)
+    opt_newton = NewtonPar(tol = 1e-9, max_iterations = 20)
+    opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds = 0.01,
+        max_steps = 100, nev = 2, newton_options = opt_newton,
+        p_min = p_span[1], p_max = p_span[2],
+        detect_bifurcation = 3, n_inversion = 4, tol_bisection_eigenvalue = 1e-8, dsmin_bisection = 1e-9);
+
+    bf = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+
+    @test length(bf.child) == 2
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index b60d82a6ce..8e1c124afa 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -55,6 +55,7 @@ end
 @safetestset "OptimizationSystem Test" include("optimizationsystem.jl")
 @safetestset "FuncAffect Test" include("funcaffect.jl")
 @safetestset "Constants Test" include("constants.jl")
+@safetestset "BifurcationKit Extension Test" include("extensions/bifurcationkit.jl")
 # Reference tests go Last
 if VERSION >= v"1.9"
     @safetestset "Latexify recipes Test" include("latexify.jl")