From e7bb024990fe27dff8771db9c1ee05c5b44049f1 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Fri, 21 Mar 2025 02:55:57 -0100
Subject: [PATCH 01/27] Change allow_symbolic to default to true

From the longer conversation. The reason for not defaulting to it before was a scare that the eliminated expression may be zero. This was found before, leading to NaNs in the resulting system evaluations. However, such zeros would also be problematic to the solver as well, since leaving the uneliminated term in leads to the structural jacobian not matching the true jacobian, which can make the structural DAE index different from the true index. Because of this phenomena, it is no less safe to eliminate the extra variables. But it can lead to some numerical improvements. For example, `p * y ~ p * z` is trivial at `p = 0`, but `y ~ z` is not, and therefore eliminating the `p` is more numerically robust.
---
 src/systems/systems.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 0f8633f31f..52153b3059 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -22,12 +22,12 @@ topological sort of the observed equations in `sys`.
 
 ### Optional Keyword Arguments:
 + When `simplify=true`, the `simplify` function will be applied during the tearing process.
-+ `allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``.
++ `allow_symbolic=true`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``.
 + `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations.
 """
 function structural_simplify(
         sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true,
-        allow_symbolic = false, allow_parameter = true, conservative = false, fully_determined = true,
+        allow_symbolic = true, allow_parameter = true, conservative = false, fully_determined = true,
         kwargs...)
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
     newsys′ = __structural_simplify(sys, io; simplify,

From 83ec981e45df9f86f6d2836bbe97b2c57d455715 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 24 Mar 2025 13:01:35 +0530
Subject: [PATCH 02/27] test: update some tests to account for `use_symbolic =
 true`

---
 test/components.jl                               |  2 +-
 test/implicit_discrete_system.jl                 |  2 +-
 .../structural_transformation/index_reduction.jl |  4 ++--
 test/structural_transformation/utils.jl          | 16 ++++++++--------
 4 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/test/components.jl b/test/components.jl
index 8ac40f6fbb..28e952b015 100644
--- a/test/components.jl
+++ b/test/components.jl
@@ -42,7 +42,7 @@ end
 completed_rc_model = complete(rc_model)
 @test isequal(completed_rc_model.resistor.n.i, resistor.n.i)
 @test ModelingToolkit.n_expanded_connection_equations(capacitor) == 2
-@test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 2
+@test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 1
 sys = structural_simplify(rc_model)
 @test_throws ModelingToolkit.RepeatedStructuralSimplificationError structural_simplify(sys)
 @test length(equations(sys)) == 1
diff --git a/test/implicit_discrete_system.jl b/test/implicit_discrete_system.jl
index 932b6c6981..5d137eba94 100644
--- a/test/implicit_discrete_system.jl
+++ b/test/implicit_discrete_system.jl
@@ -7,7 +7,7 @@ rng = StableRNG(22525)
 
 @testset "Correct ImplicitDiscreteFunction" begin
     @variables x(t) = 1
-    @mtkbuild sys = ImplicitDiscreteSystem([x(k) ~ x(k) * x(k - 1) - 3], t)
+    @mtkbuild sys = ImplicitDiscreteSystem([x(k) ~ abs(x(k)) * x(k - 1) - 3], t)
     tspan = (0, 10)
 
     # u[2] - u_next[1]
diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl
index d7f19e1fa2..b0687de2cc 100644
--- a/test/structural_transformation/index_reduction.jl
+++ b/test/structural_transformation/index_reduction.jl
@@ -134,8 +134,8 @@ let pss_pendulum = partial_state_selection(pendulum)
 end
 
 let sys = structural_simplify(pendulum2)
-    @test length(equations(sys)) == 5
-    @test length(unknowns(sys)) == 5
+    @test length(equations(sys)) == 4
+    @test length(unknowns(sys)) == 4
 
     u0 = [
         x => sqrt(2) / 2,
diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl
index 6dfc107cc9..09cc514ea4 100644
--- a/test/structural_transformation/utils.jl
+++ b/test/structural_transformation/utils.jl
@@ -224,20 +224,20 @@ end
             xtt = default_toterm(unwrap(D(D(x))))
             @test mapping[x] == (0 ~ 1 - x^2 - y^2)
             @test mapping[y] == (D(y) ~ yt)
-            @test mapping[D(y)] == (D(yt) ~ -g + y * λ)
-            @test mapping[D(x)] == (0 ~ -2xt * x - 2yt * y)
-            @test mapping[D(D(x))] == (xtt ~ x * λ)
-            @test length(mapping) == 5
+            @test mapping[D(y)] == (D(yt) ~ (2xtt * x + 2xt^2 + 2yt^2) / (-2y))
+            @test mapping[D(x)] == (xt ~ -yt * y / x)
+            @test mapping[D(D(x))] == (0 ~ -xtt + x * λ)
+            @test length(mapping) == 6
 
             @testset "`rename_dummy_derivatives = false`" begin
                 mapping = map_variables_to_equations(sys; rename_dummy_derivatives = false)
 
                 @test mapping[x] == (0 ~ 1 - x^2 - y^2)
                 @test mapping[y] == (D(y) ~ yt)
-                @test mapping[yt] == (D(yt) ~ -g + y * λ)
-                @test mapping[xt] == (0 ~ -2xt * x - 2yt * y)
-                @test mapping[xtt] == (xtt ~ x * λ)
-                @test length(mapping) == 5
+                @test mapping[yt] == (D(yt) ~ (2xtt * x + 2xt^2 + 2yt^2) / (-2y))
+                @test mapping[xt] == (xt ~ -yt * y / x)
+                @test mapping[xtt] == (0 ~ -xtt + x * λ)
+                @test length(mapping) == 6
             end
         end
         @testset "DDEs" begin

From 1049bad34ee62331cddfbf6e0126f4f4c0b2e3ef Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 27 Mar 2025 15:49:19 +0530
Subject: [PATCH 03/27] feat: add `denominators` field to `SystemStructure`

---
 src/systems/systemstructure.jl | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index d27e5c93a1..a114653ed5 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -149,8 +149,16 @@ Base.@kwdef mutable struct SystemStructure
     # or as `torn` to assert that tearing has run.
     """Graph that connects equations to variables that appear in them."""
     graph::BipartiteGraph{Int, Nothing}
-    """Graph that connects equations to the variable they will be solved for during simplification."""
+    """
+    Graph that connects equations to the variables that they can be analytically solved
+    for.
+    """
     solvable_graph::Union{BipartiteGraph{Int, Nothing}, Nothing}
+    """
+    Dict mapping `eq => var` edges in `solvable_graph` to the variables that occur in the
+    denominator when `eq` is analytically solved for `var`.
+    """
+    denominators::Dict{Pair{Int, Int}, Vector{Int}}
     """Variable types (brownian, variable, parameter) in the system."""
     var_types::Union{Vector{VariableType}, Nothing}
     """Whether the system is discrete."""
@@ -160,7 +168,7 @@ end
 function Base.copy(structure::SystemStructure)
     var_types = structure.var_types === nothing ? nothing : copy(structure.var_types)
     SystemStructure(copy(structure.var_to_diff), copy(structure.eq_to_diff),
-        copy(structure.graph), copy(structure.solvable_graph),
+        copy(structure.graph), copy(structure.solvable_graph), copy(structure.denominators),
         var_types, structure.only_discrete)
 end
 
@@ -437,7 +445,7 @@ function TearingState(sys; quick_cancel = false, check = true)
 
     ts = TearingState(sys, fullvars,
         SystemStructure(complete(var_to_diff), complete(eq_to_diff),
-            complete(graph), nothing, var_types, sys isa AbstractDiscreteSystem),
+            complete(graph), nothing, Dict(), var_types, sys isa AbstractDiscreteSystem),
         Any[])
     if sys isa DiscreteSystem
         ts = shift_discrete_system(ts)

From cc8fa2bdbab1282ee2528b69ca068a13fb74a829 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 27 Mar 2025 15:49:47 +0530
Subject: [PATCH 04/27] feat: populate `state.structure.denominators` in
 `find_eq_solvables!`

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

diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
index 14628f2958..c60187986f 100644
--- a/src/structural_transformation/utils.jl
+++ b/src/structural_transformation/utils.jl
@@ -206,9 +206,9 @@ end
 ### Structural and symbolic utilities
 ###
 
-function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = nothing;
+function find_eq_solvables!(state::TearingState, ieq::Int, to_rm = Int[], coeffs = nothing;
         may_be_zero = false,
-        allow_symbolic = false, allow_parameter = true,
+        allow_symbolic = false, allow_parameter = true, allow_algebraic = true,
         conservative = false,
         kwargs...)
     fullvars = state.fullvars
@@ -218,6 +218,7 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
     all_int_vars = true
     coeffs === nothing || empty!(coeffs)
     empty!(to_rm)
+    varsbuf = Set()
     for j in 𝑠neighbors(graph, ieq)
         var = fullvars[j]
         isirreducible(var) && (all_int_vars = false; continue)
@@ -229,13 +230,18 @@ function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = no
         if a isa Symbolic
             all_int_vars = false
             if !allow_symbolic
-                if allow_parameter
-                    all(
-                        x -> ModelingToolkit.isparameter(x) || ModelingToolkit.isconstant(x),
-                        vars(a)) || continue
-                else
+                allow_parameter || allow_algebraic || continue
+                empty!(varsbuf)
+                vars!(varsbuf, a)
+                denomvars = Int[]
+                for v in varsbuf
+                    idx = findfirst(isequal(v), fullvars)
+                    idx === nothing || push!(denomvars, idx)
+                end
+                if !allow_algebraic && !isempty(denomvars)
                     continue
                 end
+                state.structure.denominators[ieq => j] = denomvars
             end
             add_edge!(solvable_graph, ieq, j)
             continue

From a0fc653a46f2915122b48fa07c86794a457f6564 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 27 Mar 2025 15:50:16 +0530
Subject: [PATCH 05/27] feat: add `allow_algebraic`, default `allow_symbolic`
 to false

---
 .../partial_state_selection.jl                | 19 ++++++++-----
 .../symbolics_tearing.jl                      |  8 ++++--
 src/structural_transformation/utils.jl        | 16 +++++++++++
 src/systems/systems.jl                        | 27 ++++++++++++-------
 src/systems/systemstructure.jl                | 17 +++++++-----
 5 files changed, 63 insertions(+), 24 deletions(-)

diff --git a/src/structural_transformation/partial_state_selection.jl b/src/structural_transformation/partial_state_selection.jl
index 8a0ae5276e..bc949b74c8 100644
--- a/src/structural_transformation/partial_state_selection.jl
+++ b/src/structural_transformation/partial_state_selection.jl
@@ -170,11 +170,13 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
 end
 
 function dummy_derivative_graph!(state::TransformationState, jac = nothing;
-        state_priority = nothing, log = Val(false), kwargs...)
-    state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
+        state_priority = nothing, log = Val(false), allow_symbolic = false, kwargs...)
+    state.structure.solvable_graph === nothing &&
+        find_solvables!(state; allow_symbolic, kwargs...)
     complete!(state.structure)
-    var_eq_matching = complete(pantelides!(state; kwargs...))
-    dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority, log)
+    var_eq_matching = complete(pantelides!(state; allow_symbolic, kwargs...))
+    dummy_derivative_graph!(
+        state.structure, var_eq_matching, jac, state_priority, log; allow_symbolic)
 end
 
 struct DummyDerivativeSummary
@@ -184,7 +186,7 @@ end
 
 function dummy_derivative_graph!(
         structure::SystemStructure, var_eq_matching, jac = nothing,
-        state_priority = nothing, ::Val{log} = Val(false)) where {log}
+        state_priority = nothing, ::Val{log} = Val(false); allow_symbolic = false) where {log}
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -342,7 +344,12 @@ function dummy_derivative_graph!(
         @warn "The number of dummy derivatives ($n_dummys) does not match the number of differentiated equations ($n_diff_eqs)."
     end
 
-    ret = tearing_with_dummy_derivatives(structure, BitSet(dummy_derivatives))
+    dummy_derivatives_set = BitSet(dummy_derivatives)
+    if !allow_symbolic
+        make_differential_denominators_unsolvable!(structure, dummy_derivatives_set)
+    end
+
+    ret = tearing_with_dummy_derivatives(structure, dummy_derivatives_set)
     if log
         (ret..., DummyDerivativeSummary(var_dummy_scc, var_state_priority))
     else
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 552c6d13c3..75bc30631f 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -992,8 +992,12 @@ end
     ndims = ndims(arr)
 end
 
-function tearing(state::TearingState; kwargs...)
-    state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
+function tearing(state::TearingState; allow_symbolic = false, kwargs...)
+    state.structure.solvable_graph === nothing &&
+        find_solvables!(state; allow_symbolic, kwargs...)
+    if !allow_symbolic
+        make_differential_denominators_unsolvable!(state.structure)
+    end
     complete!(state.structure)
     tearing_with_dummy_derivatives(state.structure, ())
 end
diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
index c60187986f..3a9f144803 100644
--- a/src/structural_transformation/utils.jl
+++ b/src/structural_transformation/utils.jl
@@ -275,6 +275,22 @@ function find_eq_solvables!(state::TearingState, ieq::Int, to_rm = Int[], coeffs
     all_int_vars, term
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Remove edges in `structure.solvable_graph` that require differential variables in the
+denominator to solve. `additional_algevars` is a collection of integers corresponding to
+differential variables that should be considered as algebraic for the purpose of this
+transformation.
+"""
+function make_differential_denominators_unsolvable!(
+        structure::SystemStructure, additional_algevars = ())
+    for ((eqi, vari), denoms) in structure.denominators
+        all(i -> isalgvar(structure, i) || i in additional_algevars, denoms) && continue
+        rem_edge!(structure.solvable_graph, eqi, vari)
+    end
+end
+
 function find_solvables!(state::TearingState; kwargs...)
     @assert state.structure.solvable_graph === nothing
     eqs = equations(state)
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 52153b3059..9f26b527a4 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -18,21 +18,30 @@ Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations in `sys`.
 
 ### Optional Arguments:
-+ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all `inputs` to parameters and allow them to be unconnected, i.e., simplification will allow models where `n_unknowns = n_equations - n_inputs`.
++ optional argument `io` may take a tuple `(inputs, outputs)`. This will convert all
+  `inputs` to parameters and allow them to be unconnected, i.e., simplification will
+  allow models where `n_unknowns = n_equations - n_inputs`.
 
 ### Optional Keyword Arguments:
-+ When `simplify=true`, the `simplify` function will be applied during the tearing process.
-+ `allow_symbolic=true`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``.
-+ `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations.
++ When `simplify=true`, the `simplify` function will be applied during the tearing
+  process.
++ `allow_symbolic=false`, `allow_algebraic=true`, `allow_parameter=true`, and
+  `conservative=false` limit the coefficient types during tearing. In particular,
+  `conservative=true` limits tearing to only solve for trivial linear systems where
+  the coefficient has the absolute value of ``1``. `allow_symbolic` allows arbitrary
+  symbolic coefficients. If it is false, `allow_algebraic` allows symbolic coefficients
+  involving only algebraic variables and parameters. Otherwise, `allow_parameter` only
+  allows coefficients containing parameters.
++ `fully_determined=true` controls whether or not an error will be thrown if the number
+  of equations don't match the number of inputs, outputs, and equations.
 """
 function structural_simplify(
         sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true,
-        allow_symbolic = true, allow_parameter = true, conservative = false, fully_determined = true,
-        kwargs...)
+        allow_symbolic = false, allow_algebraic = true, allow_parameter = true, conservative = false,
+        fully_determined = true, kwargs...)
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
-    newsys′ = __structural_simplify(sys, io; simplify,
-        allow_symbolic, allow_parameter, conservative, fully_determined,
-        kwargs...)
+    newsys′ = __structural_simplify(sys, io; simplify, allow_symbolic, allow_algebraic,
+        allow_parameter, conservative, fully_determined, kwargs...)
     if newsys′ isa Tuple
         @assert length(newsys′) == 2
         newsys = newsys′[1]
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index a114653ed5..e475c2ff98 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -693,7 +693,7 @@ end
 
 function _structural_simplify!(state::TearingState, io; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = false,
-        dummy_derivative = true,
+        dummy_derivative = true, allow_symbolic = false,
         kwargs...)
     if fully_determined isa Bool
         check_consistency &= fully_determined
@@ -710,24 +710,27 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     else
         input_idxs = 0:-1 # Empty range
     end
-    sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
+    sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)
     if check_consistency
         fully_determined = ModelingToolkit.check_consistency(
             state, orig_inputs; nothrow = fully_determined === nothing)
     end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
     elseif fully_determined
-        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
+        var_eq_matching = pantelides!(state; finalize = false, allow_symbolic, kwargs...)
+        if !allow_symbolic
+            StructuralTransformations.make_differential_denominators_unsolvable!(state.structure)
+        end
         sys = pantelides_reassemble(state, var_eq_matching)
         state = TearingState(sys)
-        sys, mm = ModelingToolkit.alias_elimination!(state; kwargs...)
+        sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
     else
         sys = ModelingToolkit.tearing(
-            sys, state; simplify, mm, check_consistency, kwargs...)
+            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
     end
     fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)]
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)

From ae0c750ba2308c0965c4945ed35c1a5e52ab3b82 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 27 Mar 2025 15:51:00 +0530
Subject: [PATCH 06/27] Revert "test: update some tests to account for
 `use_symbolic = true`"

This reverts commit cff3906f3b52cfcc0e5bd4ea3abfc103b14ea6b2.
---
 test/components.jl                               |  2 +-
 test/implicit_discrete_system.jl                 |  2 +-
 .../structural_transformation/index_reduction.jl |  4 ++--
 test/structural_transformation/utils.jl          | 16 ++++++++--------
 4 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/test/components.jl b/test/components.jl
index 28e952b015..8ac40f6fbb 100644
--- a/test/components.jl
+++ b/test/components.jl
@@ -42,7 +42,7 @@ end
 completed_rc_model = complete(rc_model)
 @test isequal(completed_rc_model.resistor.n.i, resistor.n.i)
 @test ModelingToolkit.n_expanded_connection_equations(capacitor) == 2
-@test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 1
+@test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 2
 sys = structural_simplify(rc_model)
 @test_throws ModelingToolkit.RepeatedStructuralSimplificationError structural_simplify(sys)
 @test length(equations(sys)) == 1
diff --git a/test/implicit_discrete_system.jl b/test/implicit_discrete_system.jl
index 5d137eba94..932b6c6981 100644
--- a/test/implicit_discrete_system.jl
+++ b/test/implicit_discrete_system.jl
@@ -7,7 +7,7 @@ rng = StableRNG(22525)
 
 @testset "Correct ImplicitDiscreteFunction" begin
     @variables x(t) = 1
-    @mtkbuild sys = ImplicitDiscreteSystem([x(k) ~ abs(x(k)) * x(k - 1) - 3], t)
+    @mtkbuild sys = ImplicitDiscreteSystem([x(k) ~ x(k) * x(k - 1) - 3], t)
     tspan = (0, 10)
 
     # u[2] - u_next[1]
diff --git a/test/structural_transformation/index_reduction.jl b/test/structural_transformation/index_reduction.jl
index b0687de2cc..d7f19e1fa2 100644
--- a/test/structural_transformation/index_reduction.jl
+++ b/test/structural_transformation/index_reduction.jl
@@ -134,8 +134,8 @@ let pss_pendulum = partial_state_selection(pendulum)
 end
 
 let sys = structural_simplify(pendulum2)
-    @test length(equations(sys)) == 4
-    @test length(unknowns(sys)) == 4
+    @test length(equations(sys)) == 5
+    @test length(unknowns(sys)) == 5
 
     u0 = [
         x => sqrt(2) / 2,
diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl
index 09cc514ea4..6dfc107cc9 100644
--- a/test/structural_transformation/utils.jl
+++ b/test/structural_transformation/utils.jl
@@ -224,20 +224,20 @@ end
             xtt = default_toterm(unwrap(D(D(x))))
             @test mapping[x] == (0 ~ 1 - x^2 - y^2)
             @test mapping[y] == (D(y) ~ yt)
-            @test mapping[D(y)] == (D(yt) ~ (2xtt * x + 2xt^2 + 2yt^2) / (-2y))
-            @test mapping[D(x)] == (xt ~ -yt * y / x)
-            @test mapping[D(D(x))] == (0 ~ -xtt + x * λ)
-            @test length(mapping) == 6
+            @test mapping[D(y)] == (D(yt) ~ -g + y * λ)
+            @test mapping[D(x)] == (0 ~ -2xt * x - 2yt * y)
+            @test mapping[D(D(x))] == (xtt ~ x * λ)
+            @test length(mapping) == 5
 
             @testset "`rename_dummy_derivatives = false`" begin
                 mapping = map_variables_to_equations(sys; rename_dummy_derivatives = false)
 
                 @test mapping[x] == (0 ~ 1 - x^2 - y^2)
                 @test mapping[y] == (D(y) ~ yt)
-                @test mapping[yt] == (D(yt) ~ (2xtt * x + 2xt^2 + 2yt^2) / (-2y))
-                @test mapping[xt] == (xt ~ -yt * y / x)
-                @test mapping[xtt] == (0 ~ -xtt + x * λ)
-                @test length(mapping) == 6
+                @test mapping[yt] == (D(yt) ~ -g + y * λ)
+                @test mapping[xt] == (0 ~ -2xt * x - 2yt * y)
+                @test mapping[xtt] == (xtt ~ x * λ)
+                @test length(mapping) == 5
             end
         end
         @testset "DDEs" begin

From 828bcdb00750de649c242ff17cc22935cf98eaa2 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 27 Mar 2025 18:18:14 +0530
Subject: [PATCH 07/27] feat: propagate and implement `allow_algebraic`

---
 src/structural_transformation/partial_state_selection.jl | 4 +---
 src/structural_transformation/symbolics_tearing.jl       | 9 +++------
 src/systems/systemstructure.jl                           | 6 ++----
 3 files changed, 6 insertions(+), 13 deletions(-)

diff --git a/src/structural_transformation/partial_state_selection.jl b/src/structural_transformation/partial_state_selection.jl
index bc949b74c8..ed71e540c8 100644
--- a/src/structural_transformation/partial_state_selection.jl
+++ b/src/structural_transformation/partial_state_selection.jl
@@ -345,9 +345,7 @@ function dummy_derivative_graph!(
     end
 
     dummy_derivatives_set = BitSet(dummy_derivatives)
-    if !allow_symbolic
-        make_differential_denominators_unsolvable!(structure, dummy_derivatives_set)
-    end
+    make_differential_denominators_unsolvable!(structure, dummy_derivatives_set)
 
     ret = tearing_with_dummy_derivatives(structure, dummy_derivatives_set)
     if log
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 75bc30631f..62886e2d95 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -992,12 +992,9 @@ end
     ndims = ndims(arr)
 end
 
-function tearing(state::TearingState; allow_symbolic = false, kwargs...)
-    state.structure.solvable_graph === nothing &&
-        find_solvables!(state; allow_symbolic, kwargs...)
-    if !allow_symbolic
-        make_differential_denominators_unsolvable!(state.structure)
-    end
+function tearing(state::TearingState; kwargs...)
+    state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
+    make_differential_denominators_unsolvable!(state.structure)
     complete!(state.structure)
     tearing_with_dummy_derivatives(state.structure, ())
 end
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index e475c2ff98..187db79ac4 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -719,10 +719,8 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
         sys = ModelingToolkit.dummy_derivative(
             sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
     elseif fully_determined
-        var_eq_matching = pantelides!(state; finalize = false, allow_symbolic, kwargs...)
-        if !allow_symbolic
-            StructuralTransformations.make_differential_denominators_unsolvable!(state.structure)
-        end
+        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
+        StructuralTransformations.make_differential_denominators_unsolvable!(state.structure)
         sys = pantelides_reassemble(state, var_eq_matching)
         state = TearingState(sys)
         sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)

From 8a304b34248e252169a34d91a1118489954941cd Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 27 Mar 2025 18:07:37 +0530
Subject: [PATCH 08/27] test: update tests

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

diff --git a/test/components.jl b/test/components.jl
index 8ac40f6fbb..428f226f91 100644
--- a/test/components.jl
+++ b/test/components.jl
@@ -42,7 +42,8 @@ end
 completed_rc_model = complete(rc_model)
 @test isequal(completed_rc_model.resistor.n.i, resistor.n.i)
 @test ModelingToolkit.n_expanded_connection_equations(capacitor) == 2
-@test length(equations(structural_simplify(rc_model, allow_parameter = false))) == 2
+@test length(equations(structural_simplify(
+    rc_model, allow_algebraic = false, allow_parameter = false))) == 2
 sys = structural_simplify(rc_model)
 @test_throws ModelingToolkit.RepeatedStructuralSimplificationError structural_simplify(sys)
 @test length(equations(sys)) == 1

From b20fc7f32e976ae59bc0da4c2466ce731593d5b2 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 28 Mar 2025 10:37:06 +0530
Subject: [PATCH 09/27] fix: account for `allow_algebraic` in
 `partial_state_selection_graph!`

---
 src/structural_transformation/partial_state_selection.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/structural_transformation/partial_state_selection.jl b/src/structural_transformation/partial_state_selection.jl
index ed71e540c8..713c48456d 100644
--- a/src/structural_transformation/partial_state_selection.jl
+++ b/src/structural_transformation/partial_state_selection.jl
@@ -1,6 +1,6 @@
 function partial_state_selection_graph!(state::TransformationState)
     find_solvables!(state; allow_symbolic = true)
-    var_eq_matching = complete(pantelides!(state))
+    var_eq_matching = complete(pantelides!(state; allow_algebraic = false))
     complete!(state.structure)
     partial_state_selection_graph!(state.structure, var_eq_matching)
 end

From c7e923d3e586cbe84917680cd8f7ae5a62277faf Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 12:58:29 +0530
Subject: [PATCH 10/27] fix: remove CSE hack

---
 .../symbolics_tearing.jl                      | 89 +++----------------
 test/structural_transformation/utils.jl       |  4 +-
 2 files changed, 12 insertions(+), 81 deletions(-)

diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 62886e2d95..6ffeb28373 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -708,7 +708,7 @@ Update the system equations, unknowns, and observables after simplification.
 """
 function update_simplified_system!(
         state::TearingState, neweqs, solved_eqs, dummy_sub, var_eq_matching, extra_unknowns;
-        cse_hack = true, array_hack = true)
+        array_hack = true)
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure
     diff_to_var = invview(var_to_diff)
 
@@ -732,8 +732,8 @@ function update_simplified_system!(
     unknowns = [unknowns; extra_unknowns]
     @set! sys.unknowns = unknowns
 
-    obs, subeqs, deps = cse_and_array_hacks(
-        sys, obs, solved_eqs, unknowns, neweqs; cse = cse_hack, array = array_hack)
+    obs, subeqs, deps = array_var_hack(
+        sys, obs, solved_eqs, unknowns, neweqs; array = array_hack)
 
     @set! sys.eqs = neweqs
     @set! sys.observed = obs
@@ -775,7 +775,7 @@ appear in the system. Algebraic variables are variables that are not
 differential variables.
 """
 function tearing_reassemble(state::TearingState, var_eq_matching,
-        full_var_eq_matching = nothing; simplify = false, mm = nothing, cse_hack = true, array_hack = true)
+        full_var_eq_matching = nothing; simplify = false, mm = nothing, array_hack = true)
     extra_vars = Int[]
     if full_var_eq_matching !== nothing
         for v in 𝑑vertices(state.structure.graph)
@@ -811,7 +811,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
         state, var_eq_matching, eq_ordering, var_ordering, nelim_eq, nelim_var)
 
     sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_eq_matching,
-        extra_unknowns; cse_hack, array_hack)
+        extra_unknowns; array_hack)
 
     @set! state.sys = sys
     @set! sys.tearing_state = state
@@ -819,13 +819,6 @@ function tearing_reassemble(state::TearingState, var_eq_matching,
 end
 
 """
-# HACK 1
-
-Since we don't support array equations, any equation of the sort `x[1:n] ~ f(...)[1:n]`
-gets turned into `x[1] ~ f(...)[1], x[2] ~ f(...)[2]`. Repeatedly calling `f` gets
-_very_ expensive. this hack performs a limited form of CSE specifically for this case to
-avoid the unnecessary cost. This and the below hack are implemented simultaneously
-
 # HACK 2
 
 Add equations for array observed variables. If `p[i] ~ (...)` are equations, add an
@@ -834,12 +827,7 @@ if all `p[i]` are present and the unscalarized form is used in any equation (obs
 not) we first count the number of times the scalarized form of each observed variable
 occurs in observed equations (and unknowns if it's split).
 """
-function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, array = true)
-    # HACK 1
-    # mapping of rhs to temporary CSE variable
-    # `f(...) => tmpvar` in above example
-    rhs_to_tempvar = Dict()
-
+function array_var_hack(sys, obs, subeqs, unknowns, neweqs; array = true)
     # HACK 2
     # map of array observed variable (unscalarized) to number of its
     # scalarized terms that appear in observed equations
@@ -851,36 +839,6 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr
         rhs = eq.rhs
         vars!(all_vars, rhs)
 
-        # HACK 1
-        if cse && is_getindexed_array(rhs)
-            rhs_arr = arguments(rhs)[1]
-            iscall(rhs_arr) && operation(rhs_arr) isa Symbolics.Operator && continue
-            if !haskey(rhs_to_tempvar, rhs_arr)
-                tempvar = gensym(Symbol(lhs))
-                N = length(rhs_arr)
-                tempvar = unwrap(Symbolics.variable(
-                    tempvar; T = Symbolics.symtype(rhs_arr)))
-                tempvar = setmetadata(
-                    tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr))
-                tempeq = tempvar ~ rhs_arr
-                rhs_to_tempvar[rhs_arr] = tempvar
-                push!(obs, tempeq)
-                push!(subeqs, tempeq)
-            end
-
-            # getindex_wrapper is used because `observed2graph` treats `x` and `x[i]` as different,
-            # so it doesn't find a dependency between this equation and `tempvar ~ rhs_arr`
-            # which fails the topological sort
-            neweq = lhs ~ getindex_wrapper(
-                rhs_to_tempvar[rhs_arr], Tuple(arguments(rhs)[2:end]))
-            obs[i] = neweq
-            subeqi = findfirst(isequal(eq), subeqs)
-            if subeqi !== nothing
-                subeqs[subeqi] = neweq
-            end
-        end
-        # end HACK 1
-
         array || continue
         iscall(lhs) || continue
         operation(lhs) === getindex || continue
@@ -891,33 +849,6 @@ function cse_and_array_hacks(sys, obs, subeqs, unknowns, neweqs; cse = true, arr
         continue
     end
 
-    # Also do CSE for `equations(sys)`
-    if cse
-        for (i, eq) in enumerate(neweqs)
-            (; lhs, rhs) = eq
-            is_getindexed_array(rhs) || continue
-            rhs_arr = arguments(rhs)[1]
-            if !haskey(rhs_to_tempvar, rhs_arr)
-                tempvar = gensym(Symbol(lhs))
-                N = length(rhs_arr)
-                tempvar = unwrap(Symbolics.variable(
-                    tempvar; T = Symbolics.symtype(rhs_arr)))
-                tempvar = setmetadata(
-                    tempvar, Symbolics.ArrayShapeCtx, Symbolics.shape(rhs_arr))
-                vars!(all_vars, rhs_arr)
-                tempeq = tempvar ~ rhs_arr
-                rhs_to_tempvar[rhs_arr] = tempvar
-                push!(obs, tempeq)
-                push!(subeqs, tempeq)
-            end
-            # don't need getindex_wrapper, but do it anyway to know that this
-            # hack took place
-            neweq = lhs ~ getindex_wrapper(
-                rhs_to_tempvar[rhs_arr], Tuple(arguments(rhs)[2:end]))
-            neweqs[i] = neweq
-        end
-    end
-
     # count variables in unknowns if they are scalarized forms of variables
     # also present as observed. e.g. if `x[1]` is an unknown and `x[2] ~ (..)`
     # is an observed equation.
@@ -1007,10 +938,10 @@ new residual equations after tearing. End users are encouraged to call [`structu
 instead, which calls this function internally.
 """
 function tearing(sys::AbstractSystem, state = TearingState(sys); mm = nothing,
-        simplify = false, cse_hack = true, array_hack = true, kwargs...)
+        simplify = false, array_hack = true, kwargs...)
     var_eq_matching, full_var_eq_matching = tearing(state)
     invalidate_cache!(tearing_reassemble(
-        state, var_eq_matching, full_var_eq_matching; mm, simplify, cse_hack, array_hack))
+        state, var_eq_matching, full_var_eq_matching; mm, simplify, array_hack))
 end
 
 """
@@ -1032,7 +963,7 @@ Perform index reduction and use the dummy derivative technique to ensure that
 the system is balanced.
 """
 function dummy_derivative(sys, state = TearingState(sys); simplify = false,
-        mm = nothing, cse_hack = true, array_hack = true, kwargs...)
+        mm = nothing, array_hack = true, kwargs...)
     jac = let state = state
         (eqs, vars) -> begin
             symeqs = EquationsView(state)[eqs]
@@ -1056,5 +987,5 @@ function dummy_derivative(sys, state = TearingState(sys); simplify = false,
     end
     var_eq_matching = dummy_derivative_graph!(state, jac; state_priority,
         kwargs...)
-    tearing_reassemble(state, var_eq_matching; simplify, mm, cse_hack, array_hack)
+    tearing_reassemble(state, var_eq_matching; simplify, mm, array_hack)
 end
diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl
index 6dfc107cc9..47dcdafdce 100644
--- a/test/structural_transformation/utils.jl
+++ b/test/structural_transformation/utils.jl
@@ -121,7 +121,7 @@ end
         @named sys = ODESystem(
             [D(x) ~ z[1] + z[2] + foo(z)[1], y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t)
 
-        sys1 = structural_simplify(sys; cse_hack = false)
+        sys1 = structural_simplify(sys)
         @test length(observed(sys1)) == 6
         @test !any(observed(sys1)) do eq
             iscall(eq.rhs) &&
@@ -142,7 +142,7 @@ end
         @named sys = ODESystem(
             [D(x) ~ z[1] + z[2] + foo(z)[1] + w, y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t)
 
-        sys1 = structural_simplify(sys; cse_hack = false, fully_determined = false)
+        sys1 = structural_simplify(sys; fully_determined = false)
         @test length(observed(sys1)) == 6
         @test !any(observed(sys1)) do eq
             iscall(eq.rhs) &&

From 77bea8a783105c1d3d1c7dd10e2861c6f9bd9657 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 13:48:56 +0530
Subject: [PATCH 11/27] feat: add `remove_denominators`

---
 src/utils.jl | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/src/utils.jl b/src/utils.jl
index 1884a91c19..7a29301b4d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1299,3 +1299,16 @@ function var_in_varlist(var, varlist::AbstractSet, iv)
            # delayed variables
            (isdelay(var, iv) && var_in_varlist(operation(var)(iv), varlist, iv))
 end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Transform `expr` to have a common denominator and remove it.
+"""
+function remove_denominators(expr)
+    expr = simplify_fractions(expr)
+    if iscall(expr) && operation(expr) == (/)
+        expr = first(arguments(expr))
+    end
+    return expr
+end

From 148d3dd968ca37e69bc49e7391aef662caf69d6d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 13:48:38 +0530
Subject: [PATCH 12/27] refactor: remove denominators in `generate_function`
 and `calculate_jacobian`

---
 src/systems/diffeqs/abstractodesystem.jl | 17 ++++++++++++++---
 1 file changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index b89dc79722..83b2377f1f 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -82,8 +82,14 @@ function calculate_jacobian(sys::AbstractODESystem;
                 jac[i, j] = 0
             end
         end
+
+        (Is, Js, Vs) = findnz(jac)
+        for (i, j) in zip(Is, Js)
+            jac[i, j] = remove_denominators(jac[i, j])
+        end
     else
         jac = jacobian(rhs, dvs, simplify = simplify)
+        jac = remove_denominators.(jac)
     end
 
     if isequal(dvs, unknowns(sys))
@@ -198,12 +204,17 @@ function generate_function(sys::AbstractODESystem, dvs = unknowns(sys),
                nothing,
         isdde = false,
         kwargs...)
-    eqs = [eq for eq in equations(sys)]
-    if !implicit_dae
+    eqs = full_equations(sys)
+    if implicit_dae
+        rhss = [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs]
+        rhss = remove_denominators.(rhss)
+    else
         check_operator_variables(eqs, Differential)
         check_lhs(eqs, Differential, Set(dvs))
+        alge_idxs = is_alg_equation.(eqs)
+        rhss = [eq.rhs for eq in eqs]
+        rhss[alge_idxs] .= remove_denominators.(rhss[alge_idxs])
     end
-
     rhss = implicit_dae ? [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] :
            [eq.rhs for eq in eqs]
 

From f5d037ebaaf15c1f59f9597209f0d4e41421ca15 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 13:49:07 +0530
Subject: [PATCH 13/27] refactor: remove denominators in `NonlinearSystem`

---
 src/systems/nonlinear/nonlinearsystem.jl | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index 856822492b..58819ce38f 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -277,8 +277,13 @@ function calculate_jacobian(sys::NonlinearSystem; sparse = false, simplify = fal
 
     if sparse
         jac = sparsejacobian(rhs, vals, simplify = simplify)
+        (Is, Js, Vs) = findnz(jac)
+        for (i, j) in zip(Is, Js)
+            jac[i, j] = remove_denominators(jac[i, j])
+        end
     else
         jac = jacobian(rhs, vals, simplify = simplify)
+        jac = remove_denominators.(jac)
     end
     get_jac(sys)[] = jac, (sparse, simplify)
     return jac
@@ -318,7 +323,8 @@ function generate_function(
         sys::NonlinearSystem, dvs = unknowns(sys), ps = parameters(
             sys; initial_parameters = true);
         scalar = false, kwargs...)
-    rhss = [deq.rhs for deq in equations(sys)]
+    rhss = [deq.rhs for deq in full_equations(sys; allow_singular = true)]
+    rhss = remove_denominators.(rhss)
     dvs′ = value.(dvs)
     if scalar
         rhss = only(rhss)

From e7a6245a7dcf4b3712544cc22b0debfd683359c7 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 16:48:19 +0530
Subject: [PATCH 14/27] test: update tests to account for removal of CSE hack

---
 test/structural_transformation/utils.jl | 33 ++++---------------------
 1 file changed, 5 insertions(+), 28 deletions(-)

diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl
index 47dcdafdce..bbcde0c4fa 100644
--- a/test/structural_transformation/utils.jl
+++ b/test/structural_transformation/utils.jl
@@ -50,7 +50,7 @@ end
     @mtkbuild sys = ODESystem(
         [D(x) ~ z[1] + z[2] + foo(z)[1], y[1] ~ 2t, y[2] ~ 3t, z ~ foo(y)], t)
     @test length(equations(sys)) == 1
-    @test length(observed(sys)) == 7
+    @test length(observed(sys)) == 6
     @test any(eq -> isequal(eq.lhs, y), observed(sys))
     @test any(eq -> isequal(eq.lhs, z), observed(sys))
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn])
@@ -75,7 +75,7 @@ end
     end
     @mtkbuild sys = ODESystem([D(x) ~ y[1] + y[2], y ~ foo(x)], t)
     @test length(equations(sys)) == 1
-    @test length(observed(sys)) == 3
+    @test length(observed(sys)) == 2
     prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [foo => _tmp_fn2])
     val[] = 0
     @test_nowarn prob.f(prob.u0, prob.p, 0.0)
@@ -88,32 +88,9 @@ end
         iscall(eq.rhs) && operation(eq.rhs) in [StructuralTransformations.getindex_wrapper,
             StructuralTransformations.change_origin]
     end
-
-    @testset "CSE hack in equations(sys)" begin
-        val[] = 0
-        @variables z(t)[1:2]
-        @mtkbuild sys = ODESystem(
-            [D(y) ~ foo(x), D(x) ~ sum(y), zeros(2) ~ foo(prod(z))], t)
-        @test length(equations(sys)) == 5
-        @test length(observed(sys)) == 2
-        prob = ODEProblem(
-            sys, [y => ones(2), z => 2ones(2), x => 3.0], (0.0, 1.0), [foo => _tmp_fn2])
-        val[] = 0
-        @test_nowarn prob.f(prob.u0, prob.p, 0.0)
-        @test val[] == 2
-
-        isys = ModelingToolkit.generate_initializesystem(sys)
-        @test length(unknowns(isys)) == 5
-        @test length(equations(isys)) == 2
-        @test !any(equations(isys)) do eq
-            iscall(eq.rhs) &&
-                operation(eq.rhs) in [StructuralTransformations.getindex_wrapper,
-                    StructuralTransformations.change_origin]
-        end
-    end
 end
 
-@testset "array and cse hacks can be disabled" begin
+@testset "array hacks can be disabled" begin
     @testset "fully_determined = true" begin
         @variables x(t) y(t)[1:2] z(t)[1:2]
         @parameters foo(::AbstractVector)[1:2]
@@ -129,7 +106,7 @@ end
         end
 
         sys2 = structural_simplify(sys; array_hack = false)
-        @test length(observed(sys2)) == 5
+        @test length(observed(sys2)) == 4
         @test !any(observed(sys2)) do eq
             iscall(eq.rhs) && operation(eq.rhs) == StructuralTransformations.change_origin
         end
@@ -150,7 +127,7 @@ end
         end
 
         sys2 = structural_simplify(sys; array_hack = false, fully_determined = false)
-        @test length(observed(sys2)) == 5
+        @test length(observed(sys2)) == 4
         @test !any(observed(sys2)) do eq
             iscall(eq.rhs) && operation(eq.rhs) == StructuralTransformations.change_origin
         end

From 8ce87ea43d92d4c9ae731c28f2a27614b76d6ecb Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 16:48:33 +0530
Subject: [PATCH 15/27] test: update test to account for `allow_algebraic`

---
 test/nonlinearsystem.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index a315371141..187d6b162a 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -300,8 +300,7 @@ sys = structural_simplify(ns; conservative = true)
     ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
     @mtkbuild ns = NonlinearSystem(eqs)
 
-    @test isequal(calculate_jacobian(ns), [(-1 - z + ρ)*σ -x*σ
-                                           2x*(-z + ρ) -β-(x^2)])
+    @test isequal(calculate_jacobian(ns), [2x])
     # solve without analytical jacobian
     prob = NonlinearProblem(ns, guesses, ps)
     sol = solve(prob, NewtonRaphson())

From 4a7b12d11890f8d813c9f72e2e36fb0578d6b2e9 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 23:13:28 +0530
Subject: [PATCH 16/27] feat: allow `full_equations` on a singular system

---
 src/structural_transformation/symbolics_tearing.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 6ffeb28373..9023cd42c6 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -104,7 +104,7 @@ These equations matches generated numerical code.
 
 See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref).
 """
-function full_equations(sys::AbstractSystem; simplify = false)
+function full_equations(sys::AbstractSystem; simplify = false, allow_singular = false)
     empty_substitutions(sys) && return equations(sys)
     substitutions = get_substitutions(sys)
     substitutions.subed_eqs === nothing || return substitutions.subed_eqs
@@ -119,7 +119,7 @@ function full_equations(sys::AbstractSystem; simplify = false)
                 eq = 0 ~ eq.rhs - eq.lhs
             end
             rhs = tearing_sub(eq.rhs, solved, simplify)
-            if rhs isa Symbolic
+            if rhs isa Symbolic || allow_singular
                 return 0 ~ rhs
             else # a number
                 error("tearing failed because the system is singular")

From c988e2bd738d4e0c3e4a73fe5651b13dbeefdc29 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 23:08:46 +0530
Subject: [PATCH 17/27] refactor: default `allow_algebraic` to
 `fully_determined`

---
 src/systems/diffeqs/abstractodesystem.jl |  3 ++-
 src/systems/systems.jl                   |  7 ++++---
 src/systems/systemstructure.jl           | 11 +++++++++--
 3 files changed, 15 insertions(+), 6 deletions(-)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index 83b2377f1f..a56e524c82 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -1455,6 +1455,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
         allow_incomplete = false,
         force_time_independent = false,
         algebraic_only = false,
+        allow_algebraic = nothing,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
@@ -1482,7 +1483,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
     end
 
     if simplify_system
-        isys = structural_simplify(isys; fully_determined)
+        isys = structural_simplify(isys; fully_determined, allow_algebraic)
     end
 
     meta = get_metadata(isys)
diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index 9f26b527a4..bb2eaab58f 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -25,19 +25,20 @@ topological sort of the observed equations in `sys`.
 ### Optional Keyword Arguments:
 + When `simplify=true`, the `simplify` function will be applied during the tearing
   process.
-+ `allow_symbolic=false`, `allow_algebraic=true`, `allow_parameter=true`, and
++ `allow_symbolic=false`, `allow_algebraic=nothing`, `allow_parameter=true`, and
   `conservative=false` limit the coefficient types during tearing. In particular,
   `conservative=true` limits tearing to only solve for trivial linear systems where
   the coefficient has the absolute value of ``1``. `allow_symbolic` allows arbitrary
   symbolic coefficients. If it is false, `allow_algebraic` allows symbolic coefficients
   involving only algebraic variables and parameters. Otherwise, `allow_parameter` only
-  allows coefficients containing parameters.
+  allows coefficients containing parameters. `allow_algebraic` defaults to
+  `fully_determined` if `nothing`.
 + `fully_determined=true` controls whether or not an error will be thrown if the number
   of equations don't match the number of inputs, outputs, and equations.
 """
 function structural_simplify(
         sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true,
-        allow_symbolic = false, allow_algebraic = true, allow_parameter = true, conservative = false,
+        allow_symbolic = false, allow_algebraic = nothing, allow_parameter = true, conservative = false,
         fully_determined = true, kwargs...)
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
     newsys′ = __structural_simplify(sys, io; simplify, allow_symbolic, allow_algebraic,
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 187db79ac4..2223551eb5 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -693,7 +693,7 @@ end
 
 function _structural_simplify!(state::TearingState, io; simplify = false,
         check_consistency = true, fully_determined = true, warn_initialize_determined = false,
-        dummy_derivative = true, allow_symbolic = false,
+        dummy_derivative = true, allow_symbolic = false, allow_algebraic = nothing,
         kwargs...)
     if fully_determined isa Bool
         check_consistency &= fully_determined
@@ -710,11 +710,18 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     else
         input_idxs = 0:-1 # Empty range
     end
-    sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)
+    # use `false` for alias elimination since it doesn't really care about `allow_alebraic`
+    # anyway
+    _allow_algebraic = something(allow_algebraic, false)
+    sys, mm = ModelingToolkit.alias_elimination!(
+        state; allow_symbolic, allow_algebraic = _allow_algebraic, kwargs...)
     if check_consistency
         fully_determined = ModelingToolkit.check_consistency(
             state, orig_inputs; nothrow = fully_determined === nothing)
     end
+    if allow_algebraic === nothing
+        allow_algebraic = fully_determined
+    end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(
             sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)

From 8f8bf396fd6017b9baa4c3482d5bb33f665e9035 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 23:12:26 +0530
Subject: [PATCH 18/27] refactor: format

---
 src/ModelingToolkit.jl                   |  3 ++-
 src/systems/diffeqs/abstractodesystem.jl | 23 +++++++++++++----------
 test/jacobiansparsity.jl                 | 10 ++++++----
 3 files changed, 21 insertions(+), 15 deletions(-)

diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 76783d6be0..9f69458528 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -303,7 +303,8 @@ export structural_simplify, expand_connections, linearize, linearization_functio
        LinearizationProblem
 export solve
 
-export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function, generate_W
+export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function,
+       generate_W
 export calculate_control_jacobian, generate_control_jacobian
 export calculate_tgrad, generate_tgrad
 export calculate_gradient, generate_gradient
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index a56e524c82..c991e0eb7e 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -147,23 +147,25 @@ end
 
 function assert_jac_length_header(sys)
     W = W_sparsity(sys)
-    identity, expr -> Func([expr.args...], [], LiteralExpr(quote
-        @assert $(findnz)($(expr.args[1]))[1:2] == $(findnz)($W)[1:2]
-        $(expr.body)
-    end))
+    identity,
+    expr -> Func([expr.args...], [],
+        LiteralExpr(quote
+            @assert $(findnz)($(expr.args[1]))[1:2] == $(findnz)($W)[1:2]
+            $(expr.body)
+        end))
 end
 
-function generate_W(sys::AbstractODESystem, γ = 1., dvs = unknowns(sys),
-        ps = parameters(sys; initial_parameters = true); 
+function generate_W(sys::AbstractODESystem, γ = 1.0, dvs = unknowns(sys),
+        ps = parameters(sys; initial_parameters = true);
         simplify = false, sparse = false, kwargs...)
     @variables ˍ₋gamma
     M = calculate_massmatrix(sys; simplify)
     sparse && (M = SparseArrays.sparse(M))
     J = calculate_jacobian(sys; simplify, sparse, dvs)
-    W = ˍ₋gamma*M + J
+    W = ˍ₋gamma * M + J
 
     p = reorder_parameters(sys, ps)
-    return build_function_wrapper(sys, W, 
+    return build_function_wrapper(sys, W,
         dvs,
         p...,
         ˍ₋gamma,
@@ -313,11 +315,12 @@ function jacobian_dae_sparsity(sys::AbstractODESystem)
     J1 + J2
 end
 
-function W_sparsity(sys::AbstractODESystem) 
+function W_sparsity(sys::AbstractODESystem)
     jac_sparsity = jacobian_sparsity(sys)
     (n, n) = size(jac_sparsity)
     M = calculate_massmatrix(sys)
-    M_sparsity = M isa UniformScaling ? sparse(I(n)) : SparseMatrixCSC{Bool, Int64}((!iszero).(M))
+    M_sparsity = M isa UniformScaling ? sparse(I(n)) :
+                 SparseMatrixCSC{Bool, Int64}((!iszero).(M))
     jac_sparsity .| M_sparsity
 end
 
diff --git a/test/jacobiansparsity.jl b/test/jacobiansparsity.jl
index f6f9853137..efa1534979 100644
--- a/test/jacobiansparsity.jl
+++ b/test/jacobiansparsity.jl
@@ -74,7 +74,7 @@ f = DiffEqBase.ODEFunction(sys, u0 = nothing, sparse = true, jac = false)
 # test when u0 is not Float64
 u0 = similar(init_brusselator_2d(xyd_brusselator), Float32)
 prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
-                                     u0, (0.0, 11.5), p)
+    u0, (0.0, 11.5), p)
 sys = complete(modelingtoolkitize(prob_ode_brusselator_2d))
 
 prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = false)
@@ -94,7 +94,8 @@ prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = true)
     @mtkbuild pend = ODESystem(eqs, t)
 
     u0 = [x => 1, y => 0]
-    prob = ODEProblem(pend, u0, (0, 11.5), [g => 1], guesses = [λ => 1], sparse = true, jac = true)
+    prob = ODEProblem(
+        pend, u0, (0, 11.5), [g => 1], guesses = [λ => 1], sparse = true, jac = true)
     jac, jac! = generate_jacobian(pend; expression = Val{false}, sparse = true)
     jac_prototype = ModelingToolkit.jacobian_sparsity(pend)
     W_prototype = ModelingToolkit.W_sparsity(pend)
@@ -109,8 +110,9 @@ prob = ODEProblem(sys, u0, (0, 11.5), sparse = true, jac = true)
     @test_throws AssertionError jac!(similar(jac_prototype, Float64), u, p, t)
 
     W, W! = generate_W(pend; expression = Val{false}, sparse = true)
-    γ = .1
+    γ = 0.1
     M = sparse(calculate_massmatrix(pend))
     @test_throws AssertionError W!(similar(jac_prototype, Float64), u, p, γ, t)
-    @test W!(similar(W_prototype, Float64), u, p, γ, t) == 0.1 * M + jac!(similar(W_prototype, Float64), u, p, t)
+    @test W!(similar(W_prototype, Float64), u, p, γ, t) ==
+          0.1 * M + jac!(similar(W_prototype, Float64), u, p, t)
 end

From 8a177ff5741bff0c3154e1837f7a87d5d7dc127f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Wed, 9 Apr 2025 23:47:14 +0530
Subject: [PATCH 19/27] fixup! refactor: default `allow_algebraic` to
 `fully_determined`

---
 src/systems/systemstructure.jl | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 2223551eb5..617a4c1abc 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -710,9 +710,7 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     else
         input_idxs = 0:-1 # Empty range
     end
-    # use `false` for alias elimination since it doesn't really care about `allow_alebraic`
-    # anyway
-    _allow_algebraic = something(allow_algebraic, false)
+    _allow_algebraic = something(allow_algebraic, true)
     sys, mm = ModelingToolkit.alias_elimination!(
         state; allow_symbolic, allow_algebraic = _allow_algebraic, kwargs...)
     if check_consistency

From 3423435d52df79898d6d7dcee6d47c18a8c80dae Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 12:36:36 +0530
Subject: [PATCH 20/27] fix: fix usage of
 `make_differential_denominators_unsolvable!`

---
 src/structural_transformation/symbolics_tearing.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index 9023cd42c6..c5b11938fc 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -925,8 +925,8 @@ end
 
 function tearing(state::TearingState; kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
-    make_differential_denominators_unsolvable!(state.structure)
     complete!(state.structure)
+    make_differential_denominators_unsolvable!(state.structure)
     tearing_with_dummy_derivatives(state.structure, ())
 end
 

From f019275498d23c9215665ad6fa49ec9bbcbaac7d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 12:36:58 +0530
Subject: [PATCH 21/27] fix: don't accept `allow_algebraic` as a kwarg to
 `InitializationProblem`

---
 src/systems/diffeqs/abstractodesystem.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
index c991e0eb7e..f2dd6f3e24 100644
--- a/src/systems/diffeqs/abstractodesystem.jl
+++ b/src/systems/diffeqs/abstractodesystem.jl
@@ -1458,7 +1458,6 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
         allow_incomplete = false,
         force_time_independent = false,
         algebraic_only = false,
-        allow_algebraic = nothing,
         kwargs...) where {iip, specialize}
     if !iscomplete(sys)
         error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`")
@@ -1486,7 +1485,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractSystem,
     end
 
     if simplify_system
-        isys = structural_simplify(isys; fully_determined, allow_algebraic)
+        isys = structural_simplify(isys; fully_determined)
     end
 
     meta = get_metadata(isys)

From 38f2b0af57c4ea36e6ad5d6d9d9459e52f11063e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 12:39:45 +0530
Subject: [PATCH 22/27] fix: fix defaulting `allow_algebraic =
 fully_determined`

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

diff --git a/src/systems/systems.jl b/src/systems/systems.jl
index bb2eaab58f..a59b507c5e 100644
--- a/src/systems/systems.jl
+++ b/src/systems/systems.jl
@@ -25,21 +25,21 @@ topological sort of the observed equations in `sys`.
 ### Optional Keyword Arguments:
 + When `simplify=true`, the `simplify` function will be applied during the tearing
   process.
-+ `allow_symbolic=false`, `allow_algebraic=nothing`, `allow_parameter=true`, and
++ `allow_symbolic=false`, `allow_algebraic=fully_determined`, `allow_parameter=true`, and
   `conservative=false` limit the coefficient types during tearing. In particular,
   `conservative=true` limits tearing to only solve for trivial linear systems where
   the coefficient has the absolute value of ``1``. `allow_symbolic` allows arbitrary
   symbolic coefficients. If it is false, `allow_algebraic` allows symbolic coefficients
   involving only algebraic variables and parameters. Otherwise, `allow_parameter` only
   allows coefficients containing parameters. `allow_algebraic` defaults to
-  `fully_determined` if `nothing`.
+  `fully_determined`.
 + `fully_determined=true` controls whether or not an error will be thrown if the number
   of equations don't match the number of inputs, outputs, and equations.
 """
 function structural_simplify(
         sys::AbstractSystem, io = nothing; additional_passes = [], simplify = false, split = true,
-        allow_symbolic = false, allow_algebraic = nothing, allow_parameter = true, conservative = false,
-        fully_determined = true, kwargs...)
+        fully_determined = true, allow_symbolic = false, allow_algebraic = fully_determined,
+        allow_parameter = true, conservative = false, kwargs...)
     isscheduled(sys) && throw(RepeatedStructuralSimplificationError())
     newsys′ = __structural_simplify(sys, io; simplify, allow_symbolic, allow_algebraic,
         allow_parameter, conservative, fully_determined, kwargs...)

From 0022f198c5f755962ef684ccce4f697d56133e3d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 12:40:02 +0530
Subject: [PATCH 23/27] test: account for `allow_algebraic` in tests

---
 test/nonlinearsystem.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index 187d6b162a..43d8c16583 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -300,7 +300,7 @@ sys = structural_simplify(ns; conservative = true)
     ps = [σ => 10.0, ρ => 26.0, β => 8 / 3]
     @mtkbuild ns = NonlinearSystem(eqs)
 
-    @test isequal(calculate_jacobian(ns), [2x])
+    @test isequal(calculate_jacobian(ns), [2x;;])
     # solve without analytical jacobian
     prob = NonlinearProblem(ns, guesses, ps)
     sol = solve(prob, NewtonRaphson())
@@ -374,7 +374,7 @@ end
     end
 
     @variables y
-    @mtkbuild sys = NonlinearSystem([0 ~ x * x - p * x + p, 0 ~ x * y + p])
+    @mtkbuild sys=NonlinearSystem([0 ~ x * x - p * x + p, 0 ~ x * y + p]) allow_algebraic=false
     @test_throws ["single equation", "unknown"] IntervalNonlinearProblem(sys, (0.0, 1.0))
     @test_throws ["single equation", "unknown"] IntervalNonlinearFunction(sys, (0.0, 1.0))
     @test_throws ["single equation", "unknown"] IntervalNonlinearProblemExpr(

From 76d9045ffb3d18e6b3924f1ea03fe23ccaec48b0 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 13:43:50 +0530
Subject: [PATCH 24/27] test: account for `allow_algebraic` in HC tests

---
 test/extensions/homotopy_continuation.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl
index 65f7d765a4..bfc2414360 100644
--- a/test/extensions/homotopy_continuation.jl
+++ b/test/extensions/homotopy_continuation.jl
@@ -33,7 +33,7 @@ end
     eqs = [0 ~ x^2 + y^2 + 2x * y
            0 ~ x^2 + 4x + 4
            0 ~ y * z + 4x^2]
-    @mtkbuild sys = NonlinearSystem(eqs)
+    @mtkbuild sys=NonlinearSystem(eqs) allow_algebraic=false
     u0 = [x => 1.0, y => 1.0, z => 1.0]
     prob = HomotopyContinuationProblem(sys, u0)
     @test prob isa NonlinearProblem
@@ -60,7 +60,7 @@ end
            0 ~ x^2 + 4x + q
            0 ~ y * z + 4x^2 + wrapper(r)]
 
-    @mtkbuild sys = NonlinearSystem(eqs)
+    @mtkbuild sys=NonlinearSystem(eqs) allow_algebraic=false
     prob = HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0],
         [p => 2.0, q => 4, r => Wrapper([1.0 1.0; 0.0 0.0])])
     @test prob.ps[p] == 2.0
@@ -78,7 +78,7 @@ end
     @parameters p[1:3]
     _x = collect(x)
     eqs = collect(0 .~ vec(sum(_x * _x'; dims = 2)) + collect(p))
-    @mtkbuild sys = NonlinearSystem(eqs)
+    @mtkbuild sys=NonlinearSystem(eqs) allow_algebraic=false
     prob = HomotopyContinuationProblem(sys, [x => ones(3)], [p => 1:3])
     @test prob[x] == ones(3)
     @test prob[p + x] == [2, 3, 4]

From ad8730c37a11a76615d6af6056d556a8085f0c7f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 14:41:02 +0530
Subject: [PATCH 25/27] fix: better handle `allow_algebraic=false`

---
 .../partial_state_selection.jl                  | 15 +++++++++------
 .../symbolics_tearing.jl                        |  5 +++--
 src/structural_transformation/utils.jl          |  6 ++++--
 src/systems/systemstructure.jl                  | 17 +++++++++++------
 4 files changed, 27 insertions(+), 16 deletions(-)

diff --git a/src/structural_transformation/partial_state_selection.jl b/src/structural_transformation/partial_state_selection.jl
index 713c48456d..1fbee42770 100644
--- a/src/structural_transformation/partial_state_selection.jl
+++ b/src/structural_transformation/partial_state_selection.jl
@@ -170,13 +170,14 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
 end
 
 function dummy_derivative_graph!(state::TransformationState, jac = nothing;
-        state_priority = nothing, log = Val(false), allow_symbolic = false, kwargs...)
+        state_priority = nothing, log = Val(false), allow_symbolic = false, allow_algebraic = true, kwargs...)
     state.structure.solvable_graph === nothing &&
-        find_solvables!(state; allow_symbolic, kwargs...)
+        find_solvables!(state; allow_symbolic, allow_algebraic, kwargs...)
     complete!(state.structure)
-    var_eq_matching = complete(pantelides!(state; allow_symbolic, kwargs...))
+    var_eq_matching = complete(pantelides!(
+        state; allow_symbolic, allow_algebraic, kwargs...))
     dummy_derivative_graph!(
-        state.structure, var_eq_matching, jac, state_priority, log; allow_symbolic)
+        state.structure, var_eq_matching, jac, state_priority, log; allow_symbolic, allow_algebraic)
 end
 
 struct DummyDerivativeSummary
@@ -186,7 +187,8 @@ end
 
 function dummy_derivative_graph!(
         structure::SystemStructure, var_eq_matching, jac = nothing,
-        state_priority = nothing, ::Val{log} = Val(false); allow_symbolic = false) where {log}
+        state_priority = nothing, ::Val{log} = Val(false); allow_symbolic = false,
+        allow_algebraic = true) where {log}
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -345,7 +347,8 @@ function dummy_derivative_graph!(
     end
 
     dummy_derivatives_set = BitSet(dummy_derivatives)
-    make_differential_denominators_unsolvable!(structure, dummy_derivatives_set)
+    make_differential_denominators_unsolvable!(
+        structure, dummy_derivatives_set; allow_algebraic)
 
     ret = tearing_with_dummy_derivatives(structure, dummy_derivatives_set)
     if log
diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl
index c5b11938fc..936098a0ed 100644
--- a/src/structural_transformation/symbolics_tearing.jl
+++ b/src/structural_transformation/symbolics_tearing.jl
@@ -926,7 +926,8 @@ end
 function tearing(state::TearingState; kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
-    make_differential_denominators_unsolvable!(state.structure)
+    make_differential_denominators_unsolvable!(
+        state.structure; allow_algebraic = get(kwargs, :allow_algebraic, true))
     tearing_with_dummy_derivatives(state.structure, ())
 end
 
@@ -939,7 +940,7 @@ instead, which calls this function internally.
 """
 function tearing(sys::AbstractSystem, state = TearingState(sys); mm = nothing,
         simplify = false, array_hack = true, kwargs...)
-    var_eq_matching, full_var_eq_matching = tearing(state)
+    var_eq_matching, full_var_eq_matching = tearing(state; kwargs...)
     invalidate_cache!(tearing_reassemble(
         state, var_eq_matching, full_var_eq_matching; mm, simplify, array_hack))
 end
diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
index 3a9f144803..388824557a 100644
--- a/src/structural_transformation/utils.jl
+++ b/src/structural_transformation/utils.jl
@@ -284,9 +284,11 @@ differential variables that should be considered as algebraic for the purpose of
 transformation.
 """
 function make_differential_denominators_unsolvable!(
-        structure::SystemStructure, additional_algevars = ())
+        structure::SystemStructure, additional_algevars = (); allow_algebraic)
     for ((eqi, vari), denoms) in structure.denominators
-        all(i -> isalgvar(structure, i) || i in additional_algevars, denoms) && continue
+        if allow_algebraic && all(i -> isalgvar(structure, i) || i in additional_algevars, denoms)
+            continue
+        end
         rem_edge!(structure.solvable_graph, eqi, vari)
     end
 end
diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl
index 617a4c1abc..299e572a65 100644
--- a/src/systems/systemstructure.jl
+++ b/src/systems/systemstructure.jl
@@ -722,18 +722,23 @@ function _structural_simplify!(state::TearingState, io; simplify = false,
     end
     if fully_determined && dummy_derivative
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     elseif fully_determined
-        var_eq_matching = pantelides!(state; finalize = false, kwargs...)
-        StructuralTransformations.make_differential_denominators_unsolvable!(state.structure)
+        var_eq_matching = pantelides!(state; finalize = false, allow_algebraic, kwargs...)
+        StructuralTransformations.make_differential_denominators_unsolvable!(
+            state.structure; allow_algebraic)
         sys = pantelides_reassemble(state, var_eq_matching)
         state = TearingState(sys)
-        sys, mm = ModelingToolkit.alias_elimination!(state; allow_symbolic, kwargs...)
+        sys, mm = ModelingToolkit.alias_elimination!(
+            state; allow_symbolic, allow_algebraic, kwargs...)
         sys = ModelingToolkit.dummy_derivative(
-            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     else
         sys = ModelingToolkit.tearing(
-            sys, state; simplify, mm, check_consistency, allow_symbolic, kwargs...)
+            sys, state; simplify, mm, check_consistency,
+            allow_symbolic, allow_algebraic, kwargs...)
     end
     fullunknowns = [map(eq -> eq.lhs, observed(sys)); unknowns(sys)]
     @set! sys.observed = ModelingToolkit.topsort_equations(observed(sys), fullunknowns)

From 66034d6d330b055f8aa6a9fbffe0b1a7d1835e3d Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 14:49:55 +0530
Subject: [PATCH 26/27] fix: check edge existence before removing

---
 src/structural_transformation/utils.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl
index 388824557a..51802fc88f 100644
--- a/src/structural_transformation/utils.jl
+++ b/src/structural_transformation/utils.jl
@@ -286,7 +286,9 @@ transformation.
 function make_differential_denominators_unsolvable!(
         structure::SystemStructure, additional_algevars = (); allow_algebraic)
     for ((eqi, vari), denoms) in structure.denominators
-        if allow_algebraic && all(i -> isalgvar(structure, i) || i in additional_algevars, denoms)
+        if allow_algebraic &&
+           all(i -> isalgvar(structure, i) || i in additional_algevars, denoms) ||
+           !has_edge(structure.solvable_graph, BipartiteEdge(eqi, vari))
             continue
         end
         rem_edge!(structure.solvable_graph, eqi, vari)

From 25917bdcd9d92ae9c56e8d47487977220fdc292f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 11 Apr 2025 14:50:33 +0530
Subject: [PATCH 27/27] test: account for `allow_algebraic` in tests

---
 test/initializationsystem.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl
index 0f40d4eaf1..dfefd4dcbe 100644
--- a/test/initializationsystem.jl
+++ b/test/initializationsystem.jl
@@ -769,7 +769,7 @@ end
         eqs = [0 ~ σ * (y - x),
             0 ~ x * (ρ - z) - y,
             0 ~ x * y - β * z]
-        @mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+        @mtkbuild ns=NonlinearSystem(eqs, [x, y, z], [σ, ρ, β]) allow_algebraic=false
 
         prob = NonlinearProblem(ns, [])
         @test prob.f.initialization_data.update_initializeprob! === nothing
@@ -806,7 +806,7 @@ end
         eqs = [0 ~ p * (y - x),
             0 ~ x * (q - z) - y,
             0 ~ x * y - c * z]
-        @mtkbuild sys = NonlinearSystem(eqs; initialization_eqs = [p^2 + q^2 + 2p * q ~ 0])
+        @mtkbuild sys=NonlinearSystem(eqs; initialization_eqs = [p^2 + q^2 + 2p * q ~ 0]) allow_algebraic=false
         # @mtkbuild sys = NonlinearSystem(
         #     [p * x^2 + q * y^3 ~ 0, x - q ~ 0]; defaults = [q => missing],
         #     guesses = [q => 1.0], initialization_eqs = [p^2 + q^2 + 2p * q ~ 0])