From 6ea9875b7f7bb850a1e0dd2942cf17281b018dd0 Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Fri, 15 May 2026 11:32:34 -0400
Subject: [PATCH 01/15] feat(alias_elimination): detect `x ~ -y` perfect
 aliases
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Extend `find_perfect_aliases!` to recognize the negated alias pattern
`c1*v1 + c2*v2 ~ 0` with `c1 == c2` (i.e. `x ~ -y`) in addition to the
existing `c1 == -c2` case (`x ~ y`). Each candidate alias edge carries
an `Int8` sign, and a BFS from each group's chosen target propagates
signs to every member -- so transitive chains including double
negation (`x ~ -y` and `y ~ -z` ⇒ `x ~ z`) are resolved correctly.

When a group's signed edges are inconsistent (e.g. `x ~ y` and `x ~ -y`
both appear), every non-sticky variable in that group is substituted
with `0` and a `v ~ 0` equation is added to `additional_observed`;
irreducibles remain unknowns and the surviving alias equations carry
the resulting `irr ~ 0` constraint into the system. Substitution rhs
and derivative-cascade substitutions both honor the propagated sign.

Tests cover: basic `y ~ -x`, mixed chain `x ~ y, y ~ -z`,
double-negation transitivity, negated alias with irreducible target,
and the sign-conflict force-to-zero case.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
---
 src/systems/alias_elimination.jl | 191 +++++++++++++++++++++++--------
 test/reduction.jl                |  65 +++++++++++
 2 files changed, 209 insertions(+), 47 deletions(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index f2b94fb63a..efdf5b2e78 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -92,10 +92,12 @@ function find_perfect_aliases!(
 
     # Not `IntDisjointSet` because we don't want singleton sets for every single variable
     alias_groups = DisjointSet{Int}()
-    # Candidate alias equations `(ieq, v1_idx, v2_idx)`. Removal is decided below once
-    # each group's target is known: equations with a non-target irreducible endpoint
-    # must stay so the remaining irreducibles are still constrained to the target.
-    candidate_eqs = Tuple{Int, Int, Int}[]
+    # Candidate alias equations `(ieq, v1_idx, v2_idx, edge_sign)`. `edge_sign == +1`
+    # encodes `v1 ~ v2` (coefficients sum to zero); `edge_sign == -1` encodes
+    # `v1 ~ -v2` (coefficients are equal). Removal is decided below once each group's
+    # target is known: equations with a non-target irreducible endpoint must stay so
+    # the remaining irreducibles are still constrained to the target.
+    candidate_eqs = Tuple{Int, Int, Int, Int8}[]
 
     for ieq in 1:nsrcs(graph)
         snbors = 𝑠neighbors(graph, ieq)
@@ -106,15 +108,23 @@ function find_perfect_aliases!(
         eq = eqs[ieq]
         v1 = fullvars[snbors[1]]
         v2 = fullvars[snbors[2]]
+        local edge_sign::Int8
         Moshi.Match.@match eq.rhs begin
             BSImpl.AddMul(; variant, coeff, dict) && if variant == SU.AddMulVariant.ADD end => begin
                 SU._iszero(coeff) || continue
                 length(dict) == 2 || continue
-                haskey(dict, v1) && haskey(dict, v2) && SU._iszero(dict[v1] + dict[v2]) || continue
+                haskey(dict, v1) && haskey(dict, v2) || continue
+                if SU._iszero(dict[v1] + dict[v2])
+                    edge_sign = Int8(1)   # v1 ~ v2
+                elseif SU._iszero(dict[v1] - dict[v2])
+                    edge_sign = Int8(-1)  # v1 ~ -v2
+                else
+                    continue
+                end
             end
             _ => continue
         end
-        push!(candidate_eqs, (ieq, snbors[1], snbors[2]))
+        push!(candidate_eqs, (ieq, snbors[1], snbors[2], edge_sign))
         push!(alias_groups, snbors[1])
         push!(alias_groups, snbors[2])
         union!(alias_groups, snbors[1], snbors[2])
@@ -133,65 +143,152 @@ function find_perfect_aliases!(
         group_target[root] = pick_alias_target(fullvars, group_vars, state_priorities, irreducibles, graph)
     end
 
-    # Queue an alias equation for removal only if both endpoints collapse onto the
-    # target after non-irreducibles are substituted -- i.e. the equation becomes
-    # `T ~ T`. Any equation with a non-target irreducible endpoint is kept; when the
-    # other endpoint is a non-irreducible, the existing substitution machinery below
-    # rewrites the kept equation into `I ~ T` form automatically.
-    for (ieq, v1, v2) in candidate_eqs
-        target = group_target[DataStructures.find_root!(alias_groups, v1)]
-        c1 = contains_possibly_indexed_element(irreducibles, fullvars[v1]) || state_priorities[v1] > 0 ? v1 : target
-        c2 = contains_possibly_indexed_element(irreducibles, fullvars[v2]) || state_priorities[v2] > 0 ? v2 : target
-        c1 == c2 && push!(eqs_to_rm, ieq)
+    # Build per-group adjacency from the (signed) candidate edges and run a BFS from
+    # each group's target to assign `var_sign[v] = ±1`, indicating whether `v` is
+    # aliased to `+target` or `-target`. Sign propagation correctly handles
+    # transitive chains, including double negation: `x ~ -y` and `y ~ -z` yield
+    # `var_sign[z] = +1` (i.e. `z ~ x`). If a variable is reached via two paths
+    # implying opposite signs (e.g. both `x ~ y` and `x ~ -y` are present), the
+    # group is inconsistent: every variable in it is implied to be zero. We record
+    # such groups in `conflict_groups` and substitute their non-irreducibles with
+    # the symbolic literal `0`; irreducibles stay and the surviving alias equations
+    # encode the resulting `irr ~ 0` constraint.
+    group_adj = Dict{Int, Vector{Tuple{Int, Int8}}}()
+    for (_, v1, v2, s) in candidate_eqs
+        push!(get!(() -> Tuple{Int, Int8}[], group_adj, v1), (v2, s))
+        push!(get!(() -> Tuple{Int, Int8}[], group_adj, v2), (v1, s))
+    end
+    var_sign = Dict{Int, Int8}()
+    conflict_groups = Set{Int}()
+    bfs_queue = Int[]
+    for (root, _) in alias_sets
+        target = group_target[root]
+        var_sign[target] = Int8(1)
+        empty!(bfs_queue)
+        push!(bfs_queue, target)
+        while !isempty(bfs_queue)
+            u = popfirst!(bfs_queue)
+            us = var_sign[u]
+            nbrs = get(group_adj, u, nothing)
+            nbrs === nothing && continue
+            for (n, s) in nbrs
+                implied = Int8(s * us)
+                if haskey(var_sign, n)
+                    if var_sign[n] != implied
+                        push!(conflict_groups, root)
+                    end
+                else
+                    var_sign[n] = implied
+                    push!(bfs_queue, n)
+                end
+            end
+        end
+    end
+
+    # Queue an alias equation for removal only if both endpoints collapse out after
+    # substitution -- i.e. the equation becomes `0 ~ 0` (or `T ~ T`, which is the
+    # same once written). Any equation with a non-target irreducible endpoint is
+    # kept; when the other endpoint is a non-irreducible the existing substitution
+    # machinery below rewrites the kept equation automatically (into `I ~ T` for
+    # consistent groups, or `I ~ 0` for conflict groups).
+    is_sticky = (v) -> contains_possibly_indexed_element(irreducibles, fullvars[v]) || state_priorities[v] > 0
+    for (ieq, v1, v2, _) in candidate_eqs
+        root = DataStructures.find_root!(alias_groups, v1)
+        if root in conflict_groups
+            # Both endpoints non-sticky -> both become 0 -> equation is `0 ~ 0`.
+            !is_sticky(v1) && !is_sticky(v2) && push!(eqs_to_rm, ieq)
+        else
+            target = group_target[root]
+            c1 = is_sticky(v1) ? v1 : target
+            c2 = is_sticky(v2) ? v2 : target
+            c1 == c2 && push!(eqs_to_rm, ieq)
+        end
     end
 
+    # Symbolic literal zero used as the substitution target for variables in conflict
+    # (sign-inconsistent) groups.
+    zero_sym = Symbolics.COMMON_ZERO
+    signed_sym = (s::Int8, x::SymbolicT) -> s == 1 ? x : (-1) * x
+
     eqs_to_substitute = Int[]
     for (root, group_vars) in alias_sets
         target = group_target[root]
-        state.always_present[target] = true
+        is_conflict = root in conflict_groups
+        # Only meaningful for non-conflict groups: the chosen target survives as an
+        # unknown. For conflict groups every non-sticky variable (target included) is
+        # replaced by `0`, so we don't pin the target with `always_present`.
+        is_conflict || (state.always_present[target] = true)
         for v in group_vars
-            v == target && continue
+            !is_conflict && v == target && continue
             # Irreducibles other than the target stay as unknowns; only non-irreducibles
-            # are eliminated in favor of the target.
-            if contains_possibly_indexed_element(irreducibles, fullvars[v]) || state_priorities[v] > 0
+            # are eliminated. In conflict groups they are replaced with `0`; in
+            # consistent groups they are replaced with `±target` according to the sign
+            # propagated by the BFS above.
+            if is_sticky(v)
                 state.always_present[v] = true
                 continue
             end
 
-            push!(vars_to_rm, v)
-            subs[fullvars[v]] = fullvars[target]
-            push!(state.additional_observed, fullvars[v] ~ fullvars[target])
-            aliases[v] = target
-
-            dnbors = copy(𝑑neighbors(graph, v))
-            for e in dnbors
-                snbors = 𝑠neighbors(graph, e)
-                push!(eqs_to_substitute, e)
-                Graphs.rem_edge!(graph, e, v)
-                Graphs.add_edge!(graph, e, target)
-            end
+            if is_conflict
+                push!(vars_to_rm, v)
+                subs[fullvars[v]] = zero_sym
+                push!(state.additional_observed, fullvars[v] ~ zero_sym)
 
-            dv = var_to_diff[v]
-            dtarget = var_to_diff[target]
-            while dv isa Int
-                if dtarget === nothing
-                    dtarget = StateSelection.var_derivative!(state, target)
+                dnbors = copy(𝑑neighbors(graph, v))
+                for e in dnbors
+                    push!(eqs_to_substitute, e)
+                    Graphs.rem_edge!(graph, e, v)
                 end
-                push!(vars_to_rm, dv)
 
-                subs[fullvars[dv]] = fullvars[dtarget]
-                aliases[dv] = dtarget
+                dv = var_to_diff[v]
+                while dv isa Int
+                    push!(vars_to_rm, dv)
+                    subs[fullvars[dv]] = zero_sym
 
-                dnbors = copy(𝑑neighbors(graph, dv))
+                    dnbors = copy(𝑑neighbors(graph, dv))
+                    for e in dnbors
+                        push!(eqs_to_substitute, e)
+                        Graphs.rem_edge!(graph, e, dv)
+                    end
+
+                    dv = var_to_diff[dv]
+                end
+            else
+                s = var_sign[v]
+                push!(vars_to_rm, v)
+                rhs_sym = signed_sym(s, fullvars[target])
+                subs[fullvars[v]] = rhs_sym
+                push!(state.additional_observed, fullvars[v] ~ rhs_sym)
+                aliases[v] = target
+
+                dnbors = copy(𝑑neighbors(graph, v))
                 for e in dnbors
-                    snbors = 𝑠neighbors(graph, e)
                     push!(eqs_to_substitute, e)
-                    Graphs.rem_edge!(graph, e, dv)
-                    Graphs.add_edge!(graph, e, dtarget)
+                    Graphs.rem_edge!(graph, e, v)
+                    Graphs.add_edge!(graph, e, target)
                 end
 
-                dv = var_to_diff[dv]
-                dtarget = var_to_diff[dtarget]
+                dv = var_to_diff[v]
+                dtarget = var_to_diff[target]
+                while dv isa Int
+                    if dtarget === nothing
+                        dtarget = StateSelection.var_derivative!(state, target)
+                    end
+                    push!(vars_to_rm, dv)
+
+                    subs[fullvars[dv]] = signed_sym(s, fullvars[dtarget])
+                    aliases[dv] = dtarget
+
+                    dnbors = copy(𝑑neighbors(graph, dv))
+                    for e in dnbors
+                        push!(eqs_to_substitute, e)
+                        Graphs.rem_edge!(graph, e, dv)
+                        Graphs.add_edge!(graph, e, dtarget)
+                    end
+
+                    dv = var_to_diff[dv]
+                    dtarget = var_to_diff[dtarget]
+                end
             end
         end
     end
@@ -221,7 +318,7 @@ function find_perfect_aliases!(
     let seen = Set{Tuple{Int,Int}}()
         eqs_rm_set = Set(eqs_to_rm)
         removed_additional_eqs = false
-        for (ieq, _, _) in candidate_eqs
+        for (ieq, _, _, _) in candidate_eqs
             ieq in eqs_rm_set && continue
             snbors = 𝑠neighbors(graph, ieq)
             length(snbors) == 2 || continue
diff --git a/test/reduction.jl b/test/reduction.jl
index c03d2a1bee..9d5e4d129a 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -384,3 +384,68 @@ end
         isequal(eq, 0 ~ 1 - dot(x, Symbolics.SConst([y[1], y[1]])))
     end
 end
+
+@testset "Perfect aliases detect negated form `x ~ -y`" begin
+    # Helper: numerically evaluate the RHS of an observed equation at a chosen
+    # value of the surviving unknown, to verify the implied relationship without
+    # depending on the exact symbolic representation produced by `-`.
+    eval_rhs(rhs, var, val) = Symbolics.value(Symbolics.substitute(rhs, Dict(var => val)))
+
+    @testset "basic negated alias `y ~ -x`" begin
+        @variables x(t) y(t)
+        @mtkcompile sys = System([D(x) ~ -x, y ~ -x], t; state_priorities = [x => 10])
+        @test isequal(only(unknowns(sys)), x)
+        obs_eq = only(filter(eq -> isequal(eq.lhs, y), observed(sys)))
+        @test eval_rhs(obs_eq.rhs, x, 3.0) == -3.0
+    end
+
+    @testset "mixed chain `x ~ y, y ~ -z`" begin
+        @variables x(t) y(t) z(t)
+        @mtkcompile sys = System(
+            [D(x) ~ -x, x ~ y, y ~ -z], t; state_priorities = [x => 10])
+        @test isequal(only(unknowns(sys)), x)
+        obs = observed(sys)
+        y_eq = only(filter(eq -> isequal(eq.lhs, y), obs))
+        z_eq = only(filter(eq -> isequal(eq.lhs, z), obs))
+        # y ~ x  (sign +1),  z ~ -x  (sign -1)
+        @test eval_rhs(y_eq.rhs, x, 3.0) == 3.0
+        @test eval_rhs(z_eq.rhs, x, 3.0) == -3.0
+    end
+
+    @testset "double-negation transitivity `x ~ -y, y ~ -z` ⇒ `x ~ z`" begin
+        @variables x(t) y(t) z(t)
+        @mtkcompile sys = System(
+            [D(x) ~ -x, x ~ -y, y ~ -z], t; state_priorities = [x => 10])
+        @test isequal(only(unknowns(sys)), x)
+        obs = observed(sys)
+        y_eq = only(filter(eq -> isequal(eq.lhs, y), obs))
+        z_eq = only(filter(eq -> isequal(eq.lhs, z), obs))
+        # y ~ -x  (sign -1),  z ~ x  (sign +1: two negations cancel)
+        @test eval_rhs(y_eq.rhs, x, 3.0) == -3.0
+        @test eval_rhs(z_eq.rhs, x, 3.0) == 3.0
+    end
+
+    @testset "negated alias with irreducible target" begin
+        @variables x(t) [irreducible = true]
+        @variables y(t)
+        @mtkcompile sys = System([D(x) ~ -x, y ~ -x], t)
+        # x is irreducible so it must survive; y is eliminated as `y ~ -x`.
+        @test any(isequal(x), unknowns(sys))
+        @test !any(isequal(y), unknowns(sys))
+        obs_eq = only(filter(eq -> isequal(eq.lhs, y), observed(sys)))
+        @test eval_rhs(obs_eq.rhs, x, 2.5) == -2.5
+    end
+
+    @testset "conflicting signs force both to zero" begin
+        @variables x(t) y(t) z(t)
+        # `x ~ y` and `x ~ -y` together imply `x = y = 0`. `z` carries the dynamics
+        # so the compiled system has a well-defined surviving state.
+        @mtkcompile sys = System([D(z) ~ z, x ~ y, x ~ -y], t)
+        @test isequal(only(unknowns(sys)), z)
+        obs = observed(sys)
+        x_eq = only(filter(eq -> isequal(eq.lhs, x), obs))
+        y_eq = only(filter(eq -> isequal(eq.lhs, y), obs))
+        @test iszero(Symbolics.value(x_eq.rhs))
+        @test iszero(Symbolics.value(y_eq.rhs))
+    end
+end

From f029ac5adb440eee5c4d2745cffc2b9a9cf3552a Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Fri, 15 May 2026 14:09:22 -0400
Subject: [PATCH 02/15] refactor(alias_elimination): clean up conflict-group
 handling
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Split the loop over `candidate_eqs` into separate conflict-group and
consistent-group passes. The conflict-group pass now eagerly pins one
irreducible per claimed alias equation by rewriting it to `v ~ 0` (and
stripping the bipartite edge to the other endpoint), which both:

  - correctly handles the previously-incorrect case where both
    endpoints of a conflict-group alias equation are irreducible (the
    surviving `c1*v1 + c2*v2 = 0` did not individually pin either to
    zero), and
  - eliminates the need for a separate `pinned_eqs` data structure
    living across the substitution pass.

The main alias_sets loop now only needs to distinguish "irreducible in
a conflict group" (just mark `always_present`, since the pin is already
written) from non-irreducibles (substitute to 0 / to ±target as
before). Priority>0 non-irreducibles in conflict groups are eliminated
normally instead of being kept as unconstrained unknowns.

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
---
 src/systems/alias_elimination.jl | 91 ++++++++++++++++++++++----------
 1 file changed, 63 insertions(+), 28 deletions(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index efdf5b2e78..71f3f359fe 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -150,9 +150,9 @@ function find_perfect_aliases!(
     # `var_sign[z] = +1` (i.e. `z ~ x`). If a variable is reached via two paths
     # implying opposite signs (e.g. both `x ~ y` and `x ~ -y` are present), the
     # group is inconsistent: every variable in it is implied to be zero. We record
-    # such groups in `conflict_groups` and substitute their non-irreducibles with
-    # the symbolic literal `0`; irreducibles stay and the surviving alias equations
-    # encode the resulting `irr ~ 0` constraint.
+    # such groups in `conflict_groups`; below we substitute their non-irreducibles
+    # with the symbolic literal `0` and rewrite one alias equation per irreducible
+    # to `irr ~ 0`.
     group_adj = Dict{Int, Vector{Tuple{Int, Int8}}}()
     for (_, v1, v2, s) in candidate_eqs
         push!(get!(() -> Tuple{Int, Int8}[], group_adj, v1), (v2, s))
@@ -185,46 +185,81 @@ function find_perfect_aliases!(
         end
     end
 
-    # Queue an alias equation for removal only if both endpoints collapse out after
-    # substitution -- i.e. the equation becomes `0 ~ 0` (or `T ~ T`, which is the
-    # same once written). Any equation with a non-target irreducible endpoint is
-    # kept; when the other endpoint is a non-irreducible the existing substitution
-    # machinery below rewrites the kept equation automatically (into `I ~ T` for
-    # consistent groups, or `I ~ 0` for conflict groups).
     is_sticky = (v) -> contains_possibly_indexed_element(irreducibles, fullvars[v]) || state_priorities[v] > 0
-    for (ieq, v1, v2, _) in candidate_eqs
-        root = DataStructures.find_root!(alias_groups, v1)
-        if root in conflict_groups
-            # Both endpoints non-sticky -> both become 0 -> equation is `0 ~ 0`.
-            !is_sticky(v1) && !is_sticky(v2) && push!(eqs_to_rm, ieq)
-        else
-            target = group_target[root]
-            c1 = is_sticky(v1) ? v1 : target
-            c2 = is_sticky(v2) ? v2 : target
-            c1 == c2 && push!(eqs_to_rm, ieq)
-        end
-    end
+    is_irreducible_v = (v) -> contains_possibly_indexed_element(irreducibles, fullvars[v])
 
     # Symbolic literal zero used as the substitution target for variables in conflict
     # (sign-inconsistent) groups.
     zero_sym = Symbolics.COMMON_ZERO
     signed_sym = (s::Int8, x::SymbolicT) -> s == 1 ? x : (-1) * x
 
+    # Conflict-group alias equations: in conflict groups every variable is
+    # forced to `0`. Each irreducible must stay as an unknown, so claim one
+    # alias equation per irreducible and rewrite it eagerly to `v ~ 0` (and
+    # strip the bipartite edge to the other endpoint). A conflict group has
+    # at least as many alias equations as vertices (it contains a cycle), so
+    # there is always one equation per irreducible. Remaining alias equations
+    # become `0 ~ 0` after the substitution pass and are queued for removal.
+    # Doing the rewrite here -- before the main loop walks
+    # `𝑑neighbors(graph, v)` for non-irreducibles -- keeps the pinned eq out
+    # of `eqs_to_substitute`, so it survives untouched.
+    let claimed = Set{Int}()
+        for (ieq, v1, v2, _) in candidate_eqs
+            root = DataStructures.find_root!(alias_groups, v1)
+            root in conflict_groups || continue
+            v_pin = if is_irreducible_v(v1) && !(v1 in claimed)
+                v1
+            elseif is_irreducible_v(v2) && !(v2 in claimed)
+                v2
+            else
+                nothing
+            end
+            if v_pin === nothing
+                push!(eqs_to_rm, ieq)
+            else
+                push!(claimed, v_pin)
+                for u in copy(𝑠neighbors(graph, ieq))
+                    u == v_pin && continue
+                    Graphs.rem_edge!(graph, ieq, u)
+                end
+                eqs[ieq] = fullvars[v_pin] ~ zero_sym
+                original_eqs[ieq] = fullvars[v_pin] ~ zero_sym
+            end
+        end
+    end
+
+    # Consistent-group alias equations: queue for removal only if both
+    # endpoints collapse out after substitution (the equation becomes `0 ~ 0`).
+    # An equation with a non-target irreducible endpoint must stay so the
+    # remaining irreducible is still pinned to the target after substitution.
+    for (ieq, v1, v2, _) in candidate_eqs
+        root = DataStructures.find_root!(alias_groups, v1)
+        root in conflict_groups && continue
+        target = group_target[root]
+        c1 = is_sticky(v1) ? v1 : target
+        c2 = is_sticky(v2) ? v2 : target
+        c1 == c2 && push!(eqs_to_rm, ieq)
+    end
+
     eqs_to_substitute = Int[]
     for (root, group_vars) in alias_sets
         target = group_target[root]
         is_conflict = root in conflict_groups
         # Only meaningful for non-conflict groups: the chosen target survives as an
-        # unknown. For conflict groups every non-sticky variable (target included) is
-        # replaced by `0`, so we don't pin the target with `always_present`.
+        # unknown. For conflict groups every non-irreducible variable (target
+        # included) is replaced by `0`, so we don't pin the target with
+        # `always_present` here -- irreducibles get marked individually below.
         is_conflict || (state.always_present[target] = true)
         for v in group_vars
             !is_conflict && v == target && continue
-            # Irreducibles other than the target stay as unknowns; only non-irreducibles
-            # are eliminated. In conflict groups they are replaced with `0`; in
-            # consistent groups they are replaced with `±target` according to the sign
-            # propagated by the BFS above.
-            if is_sticky(v)
+            # In consistent groups, sticky vars (irreducible OR priority>0) other
+            # than the target stay as unknowns. In conflict groups, ONLY truly
+            # irreducible vars survive as unknowns; priority>0 non-irreducibles
+            # are eliminated along with everything else. Irreducibles in conflict
+            # groups already had one alias equation rewritten to `v ~ 0` above,
+            # so we just mark them present here -- adding them to `subs` would
+            # clobber that pin.
+            if is_conflict ? is_irreducible_v(v) : is_sticky(v)
                 state.always_present[v] = true
                 continue
             end

From 33f1ba94bc097bb44dab20b13078f3738350deb5 Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Fri, 15 May 2026 14:38:35 -0400
Subject: [PATCH 03/15] test: exercise alias-elim pass directly via
 TearingState

Switch the negated-alias testset from `@mtkcompile` to `TearingState` +
`eliminate_perfect_aliases!` so we verify the behavior of this pass in
isolation rather than the cumulative output of `mtkcompile`. Assertions
now target `state.additional_observed`, `equations(state.sys)`,
`state.fullvars`, and `state.always_present`.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 test/reduction.jl | 117 ++++++++++++++++++++++++++++++----------------
 1 file changed, 78 insertions(+), 39 deletions(-)

diff --git a/test/reduction.jl b/test/reduction.jl
index 9d5e4d129a..e6a14db6e5 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -386,66 +386,105 @@ end
 end
 
 @testset "Perfect aliases detect negated form `x ~ -y`" begin
-    # Helper: numerically evaluate the RHS of an observed equation at a chosen
-    # value of the surviving unknown, to verify the implied relationship without
-    # depending on the exact symbolic representation produced by `-`.
+    # Run only `eliminate_perfect_aliases!` on a fresh `TearingState`, so we
+    # observe the output of *this* pass rather than the cumulative effect of
+    # the rest of `mtkcompile`. The pass mutates `state` in place:
+    #   - `state.additional_observed` collects `v ~ rhs` for each eliminated v
+    #   - surviving equations of the system are returned by `equations(state.sys)`
+    #   - irreducibles forced to zero by sign conflicts have one alias equation
+    #     rewritten in place to `irr ~ 0`.
+
+    # Numerically evaluate the RHS of an observed equation at a chosen value
+    # of the surviving unknown, without depending on whether the symbolic
+    # representation of `-x` is `-x` or `-1*x` etc.
     eval_rhs(rhs, var, val) = Symbolics.value(Symbolics.substitute(rhs, Dict(var => val)))
 
+    # Find the entry of `state.additional_observed` whose lhs matches `v`.
+    obs_for(state, v) = only(filter(eq -> isequal(eq.lhs, v), state.additional_observed))
+
     @testset "basic negated alias `y ~ -x`" begin
         @variables x(t) y(t)
-        @mtkcompile sys = System([D(x) ~ -x, y ~ -x], t; state_priorities = [x => 10])
-        @test isequal(only(unknowns(sys)), x)
-        obs_eq = only(filter(eq -> isequal(eq.lhs, y), observed(sys)))
-        @test eval_rhs(obs_eq.rhs, x, 3.0) == -3.0
+        @named sys = System([D(x) ~ -x, y ~ -x], t; state_priorities = [x => 10])
+        state = TearingState(sys)
+        ModelingToolkit.eliminate_perfect_aliases!(state)
+        # y ~ -x recorded as observed
+        @test eval_rhs(obs_for(state, y).rhs, x, 3.0) == -3.0
+        # The alias equation is gone; the dynamics `D(x) ~ -x` is the only one left.
+        eqs = equations(state.sys)
+        @test length(eqs) == 1
+        @test isequal(only(eqs).lhs, D(x))
+        # `y` is no longer in `fullvars` after `rm_eqs_vars!`.
+        @test !any(isequal(unwrap(y)), state.fullvars)
     end
 
     @testset "mixed chain `x ~ y, y ~ -z`" begin
         @variables x(t) y(t) z(t)
-        @mtkcompile sys = System(
+        @named sys = System(
             [D(x) ~ -x, x ~ y, y ~ -z], t; state_priorities = [x => 10])
-        @test isequal(only(unknowns(sys)), x)
-        obs = observed(sys)
-        y_eq = only(filter(eq -> isequal(eq.lhs, y), obs))
-        z_eq = only(filter(eq -> isequal(eq.lhs, z), obs))
-        # y ~ x  (sign +1),  z ~ -x  (sign -1)
-        @test eval_rhs(y_eq.rhs, x, 3.0) == 3.0
-        @test eval_rhs(z_eq.rhs, x, 3.0) == -3.0
+        state = TearingState(sys)
+        ModelingToolkit.eliminate_perfect_aliases!(state)
+        # y ~ x (sign +1) and z ~ -x (sign -1: y has +1, z has -1 via `y ~ -z`)
+        @test eval_rhs(obs_for(state, y).rhs, x, 3.0) == 3.0
+        @test eval_rhs(obs_for(state, z).rhs, x, 3.0) == -3.0
     end
 
     @testset "double-negation transitivity `x ~ -y, y ~ -z` ⇒ `x ~ z`" begin
         @variables x(t) y(t) z(t)
-        @mtkcompile sys = System(
+        @named sys = System(
             [D(x) ~ -x, x ~ -y, y ~ -z], t; state_priorities = [x => 10])
-        @test isequal(only(unknowns(sys)), x)
-        obs = observed(sys)
-        y_eq = only(filter(eq -> isequal(eq.lhs, y), obs))
-        z_eq = only(filter(eq -> isequal(eq.lhs, z), obs))
-        # y ~ -x  (sign -1),  z ~ x  (sign +1: two negations cancel)
-        @test eval_rhs(y_eq.rhs, x, 3.0) == -3.0
-        @test eval_rhs(z_eq.rhs, x, 3.0) == 3.0
+        state = TearingState(sys)
+        ModelingToolkit.eliminate_perfect_aliases!(state)
+        # y ~ -x (sign -1), z ~ x (sign +1: two negations cancel)
+        @test eval_rhs(obs_for(state, y).rhs, x, 3.0) == -3.0
+        @test eval_rhs(obs_for(state, z).rhs, x, 3.0) == 3.0
     end
 
     @testset "negated alias with irreducible target" begin
         @variables x(t) [irreducible = true]
         @variables y(t)
-        @mtkcompile sys = System([D(x) ~ -x, y ~ -x], t)
-        # x is irreducible so it must survive; y is eliminated as `y ~ -x`.
-        @test any(isequal(x), unknowns(sys))
-        @test !any(isequal(y), unknowns(sys))
-        obs_eq = only(filter(eq -> isequal(eq.lhs, y), observed(sys)))
-        @test eval_rhs(obs_eq.rhs, x, 2.5) == -2.5
+        @named sys = System([D(x) ~ -x, y ~ -x], t)
+        state = TearingState(sys)
+        ModelingToolkit.eliminate_perfect_aliases!(state)
+        # x is irreducible, so y is the eliminated one: y ~ -x is observed.
+        @test eval_rhs(obs_for(state, y).rhs, x, 2.5) == -2.5
+        # x must be marked `always_present` by the pass (so downstream keeps it).
+        x_idx = findfirst(isequal(unwrap(x)), state.fullvars)
+        @test x_idx !== nothing && state.always_present[x_idx]
     end
 
-    @testset "conflicting signs force both to zero" begin
+    @testset "conflicting signs force both to zero (no irreducibles)" begin
         @variables x(t) y(t) z(t)
-        # `x ~ y` and `x ~ -y` together imply `x = y = 0`. `z` carries the dynamics
-        # so the compiled system has a well-defined surviving state.
-        @mtkcompile sys = System([D(z) ~ z, x ~ y, x ~ -y], t)
-        @test isequal(only(unknowns(sys)), z)
-        obs = observed(sys)
-        x_eq = only(filter(eq -> isequal(eq.lhs, x), obs))
-        y_eq = only(filter(eq -> isequal(eq.lhs, y), obs))
-        @test iszero(Symbolics.value(x_eq.rhs))
-        @test iszero(Symbolics.value(y_eq.rhs))
+        # `x ~ y` and `x ~ -y` together imply `x = y = 0`. Neither is irreducible,
+        # so both are substituted to `0` (no alias-equation pin needed) and both
+        # alias equations become `0 ~ 0` and are removed.
+        @named sys = System([D(z) ~ z, x ~ y, x ~ -y], t)
+        state = TearingState(sys)
+        ModelingToolkit.eliminate_perfect_aliases!(state)
+        @test iszero(Symbolics.value(obs_for(state, x).rhs))
+        @test iszero(Symbolics.value(obs_for(state, y).rhs))
+        # Both original alias equations are removed; only `D(z) ~ z` survives.
+        eqs = equations(state.sys)
+        @test length(eqs) == 1
+        @test isequal(only(eqs).lhs, D(z))
+    end
+
+    @testset "conflicting signs with irreducible: alias eq rewritten to `v ~ 0`" begin
+        @variables x(t) [irreducible = true]
+        @variables y(t) z(t)
+        # `x ~ y` and `x ~ -y` with `x` irreducible: `x` cannot be removed, so one
+        # alias equation is rewritten to `x ~ 0` and the other is removed.
+        @named sys = System([D(z) ~ z, x ~ y, x ~ -y], t)
+        state = TearingState(sys)
+        ModelingToolkit.eliminate_perfect_aliases!(state)
+        # y is non-irreducible: substituted to 0 in observed.
+        @test iszero(Symbolics.value(obs_for(state, y).rhs))
+        # x is NOT in additional_observed (it survives as an unknown).
+        @test !any(eq -> isequal(eq.lhs, unwrap(x)), state.additional_observed)
+        # The system carries `x ~ 0` as a real equation now, alongside the
+        # untouched `D(z) ~ z`.
+        eqs = equations(state.sys)
+        @test length(eqs) == 2
+        @test any(eq -> isequal(eq.lhs, unwrap(x)) && iszero(Symbolics.value(eq.rhs)), eqs)
+        @test any(eq -> isequal(eq.lhs, D(z)), eqs)
     end
 end

From db47fc319a4745ae54100a808d732ce61e5f73f8 Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Fri, 15 May 2026 15:36:47 -0400
Subject: [PATCH 04/15] optimizations

---
 .../src/ModelingToolkitBase.jl                |   6 -
 src/ModelingToolkit.jl                        |   6 -
 src/systems/alias_elimination.jl              | 142 ++++++++++--------
 3 files changed, 83 insertions(+), 71 deletions(-)

diff --git a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
index a3b2bed3da..68fb203fc2 100644
--- a/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
+++ b/lib/ModelingToolkitBase/src/ModelingToolkitBase.jl
@@ -38,12 +38,6 @@ using SymbolicIndexingInterface
 using LinearAlgebra, SparseArrays
 using InteractiveUtils
 using DataStructures
-@static if pkgversion(DataStructures) >= v"0.19"
-    import DataStructures: IntDisjointSet
-else
-    import DataStructures: IntDisjointSets
-    const IntDisjointSet = IntDisjointSets
-end
 using Base.Threads
 using ArrayInterface
 using Setfield, ConstructionBase
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index bdba846f91..72f573f95e 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -33,12 +33,6 @@ using SymbolicIndexingInterface
 using LinearAlgebra, SparseArrays
 using InteractiveUtils
 using DataStructures
-@static if pkgversion(DataStructures) >= v"0.19"
-    import DataStructures: IntDisjointSet
-else
-    import DataStructures: IntDisjointSets
-    const IntDisjointSet = IntDisjointSets
-end
 using Base.Threads
 using Setfield, ConstructionBase
 import Libdl
diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 71f3f359fe..93e3e689f0 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -40,6 +40,57 @@ function eliminate_perfect_aliases!(state::TearingState)
     return nothing
 end
 
+"""
+    $TYPEDSIGNATURES
+
+Weighted union-find augmented with sign tracking. Merges the components of `v1`
+and `v2` under the relation `v1 ≡ edge_sign · v2` (`edge_sign ∈ ±1`), maintaining
+the invariant that `parent[v]` is always the current root and `parity[v]` is the
+sign of `v` relative to that root -- so finds are O(1).
+
+Merges are smaller-into-larger using `members` (root → member list), giving
+O((V+E) log V) amortized work across all unions. If an edge creates a contradiction
+within an already-merged component the root is recorded in `conflict_roots`; the
+flag is migrated when a conflict component is absorbed by a larger one.
+"""
+function union_with_sign!(
+        parent::Dict{Int, Int}, parity::Dict{Int, Int8},
+        members::Dict{Int, Vector{Int}}, conflict_roots::Set{Int},
+        v1::Int, v2::Int, edge_sign::Int8
+    )
+    for v in (v1, v2)
+        if !haskey(parent, v)
+            parent[v] = v
+            parity[v] = Int8(1)
+            members[v] = [v]
+        end
+    end
+    r1 = parent[v1]; s1 = parity[v1]
+    r2 = parent[v2]; s2 = parity[v2]
+    if r1 == r2
+        s1 == Int8(edge_sign * s2) || push!(conflict_roots, r1)
+        return
+    end
+    if length(members[r1]) < length(members[r2])
+        r1, r2 = r2, r1
+        s1, s2 = s2, s1
+    end
+    # Edge `v1 = edge_sign · v2` ⇒ sign of r2 relative to r1 = s1 · edge_sign · s2.
+    r2_to_r1 = Int8(s1 * edge_sign * s2)
+    r1_members = members[r1]
+    for m in members[r2]
+        parent[m] = r1
+        parity[m] = Int8(parity[m] * r2_to_r1)
+        push!(r1_members, m)
+    end
+    delete!(members, r2)
+    if r2 in conflict_roots
+        delete!(conflict_roots, r2)
+        push!(conflict_roots, r1)
+    end
+    return
+end
+
 """
     $TYPEDSIGNATURES
 
@@ -90,8 +141,16 @@ function find_perfect_aliases!(
     original_eqs = state.original_eqs
     irreducibles = get_irreducibles(sys)
 
-    # Not `IntDisjointSet` because we don't want singleton sets for every single variable
-    alias_groups = DisjointSet{Int}()
+    # Weighted union-find tracking each variable's sign relative to its
+    # component's current root. `parent[v]` is always the root (no path
+    # walking needed), `parity[v] ∈ ±1` is the sign of `v` relative to that
+    # root, `members` maps root → member list (used for smaller-to-larger
+    # merges in `union_with_sign!`), and `conflict_roots` collects roots
+    # whose edges are sign-inconsistent.
+    parent = Dict{Int, Int}()
+    parity = Dict{Int, Int8}()
+    members = Dict{Int, Vector{Int}}()
+    conflict_roots = Set{Int}()
     # Candidate alias equations `(ieq, v1_idx, v2_idx, edge_sign)`. `edge_sign == +1`
     # encodes `v1 ~ v2` (coefficients sum to zero); `edge_sign == -1` encodes
     # `v1 ~ -v2` (coefficients are equal). Removal is decided below once each group's
@@ -125,65 +184,31 @@ function find_perfect_aliases!(
             _ => continue
         end
         push!(candidate_eqs, (ieq, snbors[1], snbors[2], edge_sign))
-        push!(alias_groups, snbors[1])
-        push!(alias_groups, snbors[2])
-        union!(alias_groups, snbors[1], snbors[2])
-    end
-
-    alias_sets = Dict{Int, Vector{Int}}()
-    sizehint!(alias_sets, DataStructures.num_groups(alias_groups))
-    sizehint!(aliases, length(alias_groups) - DataStructures.num_groups(alias_groups))
-    for var in alias_groups
-        set = get!(() -> Int[], alias_sets, DataStructures.find_root!(alias_groups, var))
-        push!(set, var)
+        union_with_sign!(parent, parity, members, conflict_roots,
+            snbors[1], snbors[2], edge_sign)
     end
 
+    # After the scan above, every variable in any alias group has `parent[v] ==
+    # root` and `parity[v] = ±1` (its sign relative to the root). Pick a target
+    # per group and rebase `parity` so that `var_sign[v] = +1` means `v ~ target`
+    # and `-1` means `v ~ -target`. Sign-inconsistent groups were collected in
+    # `conflict_roots`; their non-irreducibles get substituted with the symbolic
+    # literal `0` below, and one alias equation per irreducible is rewritten to
+    # `irr ~ 0`.
     group_target = Dict{Int, Int}()
-    for (root, group_vars) in alias_sets
+    sizehint!(group_target, length(members))
+    for (root, group_vars) in members
         group_target[root] = pick_alias_target(fullvars, group_vars, state_priorities, irreducibles, graph)
     end
-
-    # Build per-group adjacency from the (signed) candidate edges and run a BFS from
-    # each group's target to assign `var_sign[v] = ±1`, indicating whether `v` is
-    # aliased to `+target` or `-target`. Sign propagation correctly handles
-    # transitive chains, including double negation: `x ~ -y` and `y ~ -z` yield
-    # `var_sign[z] = +1` (i.e. `z ~ x`). If a variable is reached via two paths
-    # implying opposite signs (e.g. both `x ~ y` and `x ~ -y` are present), the
-    # group is inconsistent: every variable in it is implied to be zero. We record
-    # such groups in `conflict_groups`; below we substitute their non-irreducibles
-    # with the symbolic literal `0` and rewrite one alias equation per irreducible
-    # to `irr ~ 0`.
-    group_adj = Dict{Int, Vector{Tuple{Int, Int8}}}()
-    for (_, v1, v2, s) in candidate_eqs
-        push!(get!(() -> Tuple{Int, Int8}[], group_adj, v1), (v2, s))
-        push!(get!(() -> Tuple{Int, Int8}[], group_adj, v2), (v1, s))
-    end
     var_sign = Dict{Int, Int8}()
-    conflict_groups = Set{Int}()
-    bfs_queue = Int[]
-    for (root, _) in alias_sets
-        target = group_target[root]
-        var_sign[target] = Int8(1)
-        empty!(bfs_queue)
-        push!(bfs_queue, target)
-        while !isempty(bfs_queue)
-            u = popfirst!(bfs_queue)
-            us = var_sign[u]
-            nbrs = get(group_adj, u, nothing)
-            nbrs === nothing && continue
-            for (n, s) in nbrs
-                implied = Int8(s * us)
-                if haskey(var_sign, n)
-                    if var_sign[n] != implied
-                        push!(conflict_groups, root)
-                    end
-                else
-                    var_sign[n] = implied
-                    push!(bfs_queue, n)
-                end
-            end
+    sizehint!(var_sign, length(parent))
+    for (root, group_vars) in members
+        target_p = parity[group_target[root]]
+        for v in group_vars
+            var_sign[v] = Int8(parity[v] * target_p)
         end
     end
+    sizehint!(aliases, length(parent) - length(members))
 
     is_sticky = (v) -> contains_possibly_indexed_element(irreducibles, fullvars[v]) || state_priorities[v] > 0
     is_irreducible_v = (v) -> contains_possibly_indexed_element(irreducibles, fullvars[v])
@@ -205,8 +230,7 @@ function find_perfect_aliases!(
     # of `eqs_to_substitute`, so it survives untouched.
     let claimed = Set{Int}()
         for (ieq, v1, v2, _) in candidate_eqs
-            root = DataStructures.find_root!(alias_groups, v1)
-            root in conflict_groups || continue
+            parent[v1] in conflict_roots || continue
             v_pin = if is_irreducible_v(v1) && !(v1 in claimed)
                 v1
             elseif is_irreducible_v(v2) && !(v2 in claimed)
@@ -233,8 +257,8 @@ function find_perfect_aliases!(
     # An equation with a non-target irreducible endpoint must stay so the
     # remaining irreducible is still pinned to the target after substitution.
     for (ieq, v1, v2, _) in candidate_eqs
-        root = DataStructures.find_root!(alias_groups, v1)
-        root in conflict_groups && continue
+        root = parent[v1]
+        root in conflict_roots && continue
         target = group_target[root]
         c1 = is_sticky(v1) ? v1 : target
         c2 = is_sticky(v2) ? v2 : target
@@ -242,9 +266,9 @@ function find_perfect_aliases!(
     end
 
     eqs_to_substitute = Int[]
-    for (root, group_vars) in alias_sets
+    for (root, group_vars) in members
         target = group_target[root]
-        is_conflict = root in conflict_groups
+        is_conflict = root in conflict_roots
         # Only meaningful for non-conflict groups: the chosen target survives as an
         # unknown. For conflict groups every non-irreducible variable (target
         # included) is replaced by `0`, so we don't pin the target with

From 9f831aa132695a17784211be1fb374e832f88556 Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Fri, 15 May 2026 15:43:37 -0400
Subject: [PATCH 05/15] cleanup

---
 src/systems/alias_elimination.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 93e3e689f0..6c121601a2 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -77,12 +77,12 @@ function union_with_sign!(
     end
     # Edge `v1 = edge_sign · v2` ⇒ sign of r2 relative to r1 = s1 · edge_sign · s2.
     r2_to_r1 = Int8(s1 * edge_sign * s2)
-    r1_members = members[r1]
-    for m in members[r2]
+    r2_members = members[r2]
+    for m in r2_members
         parent[m] = r1
         parity[m] = Int8(parity[m] * r2_to_r1)
-        push!(r1_members, m)
     end
+    append!(members[r1], r2_members)
     delete!(members, r2)
     if r2 in conflict_roots
         delete!(conflict_roots, r2)

From e66a0d874edcbd20c104cddc6018c5402ea6875d Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Fri, 15 May 2026 16:55:42 -0400
Subject: [PATCH 06/15] simplifiations

---
 src/systems/alias_elimination.jl | 242 +++++++++++++++----------------
 1 file changed, 114 insertions(+), 128 deletions(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 6c121601a2..b8a134dfe8 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -44,18 +44,18 @@ end
     $TYPEDSIGNATURES
 
 Weighted union-find augmented with sign tracking. Merges the components of `v1`
-and `v2` under the relation `v1 ≡ edge_sign · v2` (`edge_sign ∈ ±1`), maintaining
-the invariant that `parent[v]` is always the current root and `parity[v]` is the
-sign of `v` relative to that root -- so finds are O(1).
-
-Merges are smaller-into-larger using `members` (root → member list), giving
-O((V+E) log V) amortized work across all unions. If an edge creates a contradiction
-within an already-merged component the root is recorded in `conflict_roots`; the
-flag is migrated when a conflict component is absorbed by a larger one.
+and `v2` under the relation `v1 ~ edge_sign · v2` (`edge_sign ∈ ±1`)
+maintains the invariant that `parent[v]` is always the current root and `parity[v]`
+is the sign of `v` relative to that root so finds are O(1).
+
+Merges are smaller-into-larger using `members` (root → member list)
+giving O((V+E) log V) across all unions.
+A contradiction (e.g. `x ~ y, x ~ -y`) zeros out every member's `parity` in the
+affected component; downstream code treats `parity[v] == 0` as "v is forced to 0".
 """
 function union_with_sign!(
         parent::Dict{Int, Int}, parity::Dict{Int, Int8},
-        members::Dict{Int, Vector{Int}}, conflict_roots::Set{Int},
+        members::Dict{Int, Vector{Int}},
         v1::Int, v2::Int, edge_sign::Int8
     )
     for v in (v1, v2)
@@ -68,7 +68,11 @@ function union_with_sign!(
     r1 = parent[v1]; s1 = parity[v1]
     r2 = parent[v2]; s2 = parity[v2]
     if r1 == r2
-        s1 == Int8(edge_sign * s2) || push!(conflict_roots, r1)
+        if s1 != Int8(edge_sign * s2)
+            for m in members[r1]
+                parity[m] = Int8(0)
+            end
+        end
         return
     end
     if length(members[r1]) < length(members[r2])
@@ -76,18 +80,21 @@ function union_with_sign!(
         s1, s2 = s2, s1
     end
     # Edge `v1 = edge_sign · v2` ⇒ sign of r2 relative to r1 = s1 · edge_sign · s2.
+    # If either side is already zero (conflict-tainted), `r2_to_r1` is 0 and we
+    # propagate the zero through both components.
     r2_to_r1 = Int8(s1 * edge_sign * s2)
     r2_members = members[r2]
     for m in r2_members
         parent[m] = r1
         parity[m] = Int8(parity[m] * r2_to_r1)
     end
+    if r2_to_r1 == 0
+        for m in members[r1]
+            parity[m] = Int8(0)
+        end
+    end
     append!(members[r1], r2_members)
     delete!(members, r2)
-    if r2 in conflict_roots
-        delete!(conflict_roots, r2)
-        push!(conflict_roots, r1)
-    end
     return
 end
 
@@ -142,19 +149,17 @@ function find_perfect_aliases!(
     irreducibles = get_irreducibles(sys)
 
     # Weighted union-find tracking each variable's sign relative to its
-    # component's current root. `parent[v]` is always the root (no path
-    # walking needed), `parity[v] ∈ ±1` is the sign of `v` relative to that
-    # root, `members` maps root → member list (used for smaller-to-larger
-    # merges in `union_with_sign!`), and `conflict_roots` collects roots
-    # whose edges are sign-inconsistent.
+    # component's current root. `parent[v]` is always the root.
+    # `parity[v] ∈ ±1` is the sign of `v` relative to that root, or `0` if the
+    # component contains a sign contradiction (forcing every member to 0).
+    # `members` maps root → member list (used for smaller-to-larger merges).
     parent = Dict{Int, Int}()
     parity = Dict{Int, Int8}()
     members = Dict{Int, Vector{Int}}()
-    conflict_roots = Set{Int}()
-    # Candidate alias equations `(ieq, v1_idx, v2_idx, edge_sign)`. `edge_sign == +1`
-    # encodes `v1 ~ v2` (coefficients sum to zero); `edge_sign == -1` encodes
-    # `v1 ~ -v2` (coefficients are equal). Removal is decided below once each group's
-    # target is known: equations with a non-target irreducible endpoint must stay so
+    # Candidate alias equations `(ieq, v1_idx, v2_idx, edge_sign)`. 
+    # `edge_sign == +1` encodes `v1 ~ edge_sign*v2`.
+    # Removal is decided below once each group's target is known:
+    # equations with a non-target irreducible endpoint must stay so
     # the remaining irreducibles are still constrained to the target.
     candidate_eqs = Tuple{Int, Int, Int, Int8}[]
 
@@ -184,53 +189,36 @@ function find_perfect_aliases!(
             _ => continue
         end
         push!(candidate_eqs, (ieq, snbors[1], snbors[2], edge_sign))
-        union_with_sign!(parent, parity, members, conflict_roots,
-            snbors[1], snbors[2], edge_sign)
+        union_with_sign!(parent, parity, members, snbors[1], snbors[2], edge_sign)
     end
+    sizehint!(aliases, length(parent) - length(members))
 
-    # After the scan above, every variable in any alias group has `parent[v] ==
-    # root` and `parity[v] = ±1` (its sign relative to the root). Pick a target
-    # per group and rebase `parity` so that `var_sign[v] = +1` means `v ~ target`
-    # and `-1` means `v ~ -target`. Sign-inconsistent groups were collected in
-    # `conflict_roots`; their non-irreducibles get substituted with the symbolic
-    # literal `0` below, and one alias equation per irreducible is rewritten to
-    # `irr ~ 0`.
+    # After the scan above, every variable in any alias group has
+    # `parent[v] == root` and `parity[v] ∈ {-1, 0, +1}`
+    # (sign relative to root, or 0 if the component is conflict-tainted). 
+    # The main rewriting loop below picks a target per non-conflict group 
+    # and uses `s = parity[v] * parity[target]` (±1) as `v`'s sign relative to the target.
+    # For Conflict groups: every non-irreducible gets substituted with the
+    # symbolic `0`, and one alias eq per irreducible is rewritten to `irr ~ 0`.
     group_target = Dict{Int, Int}()
     sizehint!(group_target, length(members))
-    for (root, group_vars) in members
-        group_target[root] = pick_alias_target(fullvars, group_vars, state_priorities, irreducibles, graph)
-    end
-    var_sign = Dict{Int, Int8}()
-    sizehint!(var_sign, length(parent))
-    for (root, group_vars) in members
-        target_p = parity[group_target[root]]
-        for v in group_vars
-            var_sign[v] = Int8(parity[v] * target_p)
-        end
-    end
-    sizehint!(aliases, length(parent) - length(members))
 
     is_sticky = (v) -> contains_possibly_indexed_element(irreducibles, fullvars[v]) || state_priorities[v] > 0
     is_irreducible_v = (v) -> contains_possibly_indexed_element(irreducibles, fullvars[v])
 
-    # Symbolic literal zero used as the substitution target for variables in conflict
-    # (sign-inconsistent) groups.
+    # Symbolic zero used as substitution target for variables in conflict groups.
     zero_sym = Symbolics.COMMON_ZERO
     signed_sym = (s::Int8, x::SymbolicT) -> s == 1 ? x : (-1) * x
 
-    # Conflict-group alias equations: in conflict groups every variable is
-    # forced to `0`. Each irreducible must stay as an unknown, so claim one
-    # alias equation per irreducible and rewrite it eagerly to `v ~ 0` (and
-    # strip the bipartite edge to the other endpoint). A conflict group has
-    # at least as many alias equations as vertices (it contains a cycle), so
-    # there is always one equation per irreducible. Remaining alias equations
-    # become `0 ~ 0` after the substitution pass and are queued for removal.
-    # Doing the rewrite here -- before the main loop walks
-    # `𝑑neighbors(graph, v)` for non-irreducibles -- keeps the pinned eq out
-    # of `eqs_to_substitute`, so it survives untouched.
+    # Conflict-group alias equations: every variable is forced to `0`.
+    # Each irreducible must stay in unknowns, so claim one alias eq per irreducible
+    # and rewrite it eagerly to `v ~ 0` (stripping the bipartite edge to the other endpoint).
+    # A conflict group contains a cycle so has at least as many alias equations as vertices
+    # so there is always one eq per irreducible.
+    # Remaining alias equations become `0 ~ 0` after the substitution and are queued for removal.
     let claimed = Set{Int}()
         for (ieq, v1, v2, _) in candidate_eqs
-            parent[v1] in conflict_roots || continue
+            parity[v1] == 0 || continue
             v_pin = if is_irreducible_v(v1) && !(v1 in claimed)
                 v1
             elseif is_irreducible_v(v2) && !(v2 in claimed)
@@ -252,49 +240,22 @@ function find_perfect_aliases!(
         end
     end
 
-    # Consistent-group alias equations: queue for removal only if both
-    # endpoints collapse out after substitution (the equation becomes `0 ~ 0`).
-    # An equation with a non-target irreducible endpoint must stay so the
-    # remaining irreducible is still pinned to the target after substitution.
-    for (ieq, v1, v2, _) in candidate_eqs
-        root = parent[v1]
-        root in conflict_roots && continue
-        target = group_target[root]
-        c1 = is_sticky(v1) ? v1 : target
-        c2 = is_sticky(v2) ? v2 : target
-        c1 == c2 && push!(eqs_to_rm, ieq)
-    end
-
     eqs_to_substitute = Int[]
     for (root, group_vars) in members
-        target = group_target[root]
-        is_conflict = root in conflict_roots
-        # Only meaningful for non-conflict groups: the chosen target survives as an
-        # unknown. For conflict groups every non-irreducible variable (target
-        # included) is replaced by `0`, so we don't pin the target with
-        # `always_present` here -- irreducibles get marked individually below.
-        is_conflict || (state.always_present[target] = true)
-        for v in group_vars
-            !is_conflict && v == target && continue
-            # In consistent groups, sticky vars (irreducible OR priority>0) other
-            # than the target stay as unknowns. In conflict groups, ONLY truly
-            # irreducible vars survive as unknowns; priority>0 non-irreducibles
-            # are eliminated along with everything else. Irreducibles in conflict
-            # groups already had one alias equation rewritten to `v ~ 0` above,
-            # so we just mark them present here -- adding them to `subs` would
-            # clobber that pin.
-            if is_conflict ? is_irreducible_v(v) : is_sticky(v)
-                state.always_present[v] = true
-                continue
-            end
-
-            if is_conflict
+        if parity[root] == 0 # is conflict group
+            for v in group_vars
+                # In conflict groups, ONLY irreducible vars survive as unknowns
+                # Irreducibles already had equation rewritten to `v ~ 0` above
+                # so just mark them present here
+                if is_irreducible_v(v)
+                    state.always_present[v] = true
+                    continue
+                end
                 push!(vars_to_rm, v)
                 subs[fullvars[v]] = zero_sym
                 push!(state.additional_observed, fullvars[v] ~ zero_sym)
 
-                dnbors = copy(𝑑neighbors(graph, v))
-                for e in dnbors
+                for e in copy(𝑑neighbors(graph, v))
                     push!(eqs_to_substitute, e)
                     Graphs.rem_edge!(graph, e, v)
                 end
@@ -304,54 +265,79 @@ function find_perfect_aliases!(
                     push!(vars_to_rm, dv)
                     subs[fullvars[dv]] = zero_sym
 
-                    dnbors = copy(𝑑neighbors(graph, dv))
-                    for e in dnbors
+                    for e in copy(𝑑neighbors(graph, dv))
                         push!(eqs_to_substitute, e)
                         Graphs.rem_edge!(graph, e, dv)
                     end
 
                     dv = var_to_diff[dv]
                 end
-            else
-                s = var_sign[v]
-                push!(vars_to_rm, v)
-                rhs_sym = signed_sym(s, fullvars[target])
-                subs[fullvars[v]] = rhs_sym
-                push!(state.additional_observed, fullvars[v] ~ rhs_sym)
-                aliases[v] = target
-
-                dnbors = copy(𝑑neighbors(graph, v))
-                for e in dnbors
-                    push!(eqs_to_substitute, e)
-                    Graphs.rem_edge!(graph, e, v)
-                    Graphs.add_edge!(graph, e, target)
-                end
+            end
+            continue
+        end
+        # For consistent groups pick a target (survives as unknown) and rebase
+        # parities relative to it via `target_p`.
+        target = pick_alias_target(fullvars, group_vars, state_priorities, irreducibles, graph)
+        group_target[root] = target
+        target_p = parity[target]
+        state.always_present[target] = true
+        for v in group_vars
+            v == target && continue
+            # In consistent groups sticky vars (irreducible OR priority>0) other
+            # than the target stay as unknowns.
+            if is_sticky(v)
+                state.always_present[v] = true
+                continue
+            end
 
-                dv = var_to_diff[v]
-                dtarget = var_to_diff[target]
-                while dv isa Int
-                    if dtarget === nothing
-                        dtarget = StateSelection.var_derivative!(state, target)
-                    end
-                    push!(vars_to_rm, dv)
+            s = Int8(parity[v] * target_p)
+            push!(vars_to_rm, v)
+            rhs_sym = signed_sym(s, fullvars[target])
+            subs[fullvars[v]] = rhs_sym
+            push!(state.additional_observed, fullvars[v] ~ rhs_sym)
+            aliases[v] = target
+
+            for e in copy(𝑑neighbors(graph, v))
+                push!(eqs_to_substitute, e)
+                Graphs.rem_edge!(graph, e, v)
+                Graphs.add_edge!(graph, e, target)
+            end
 
-                    subs[fullvars[dv]] = signed_sym(s, fullvars[dtarget])
-                    aliases[dv] = dtarget
+            dv = var_to_diff[v]
+            dtarget = var_to_diff[target]
+            while dv isa Int
+                if dtarget === nothing
+                    dtarget = StateSelection.var_derivative!(state, target)
+                end
+                push!(vars_to_rm, dv)
 
-                    dnbors = copy(𝑑neighbors(graph, dv))
-                    for e in dnbors
-                        push!(eqs_to_substitute, e)
-                        Graphs.rem_edge!(graph, e, dv)
-                        Graphs.add_edge!(graph, e, dtarget)
-                    end
+                subs[fullvars[dv]] = signed_sym(s, fullvars[dtarget])
+                aliases[dv] = dtarget
 
-                    dv = var_to_diff[dv]
-                    dtarget = var_to_diff[dtarget]
+                for e in copy(𝑑neighbors(graph, dv))
+                    push!(eqs_to_substitute, e)
+                    Graphs.rem_edge!(graph, e, dv)
+                    Graphs.add_edge!(graph, e, dtarget)
                 end
+
+                dv = var_to_diff[dv]
+                dtarget = var_to_diff[dtarget]
             end
         end
     end
 
+    # Consistent-group alias equations: queue for removal iff both endpoints
+    # collapse (equation becomes `0 ~ 0`). An equation with a non-target
+    # irreducible endpoint must stay so the remaining irreducible is still
+    # pinned to the target after substitution.
+    for (ieq, v1, v2, _) in candidate_eqs
+        parity[v1] == 0 && continue
+        target = group_target[parent[v1]]
+        c1 = is_sticky(v1) ? v1 : target
+        c2 = is_sticky(v2) ? v2 : target
+        c1 == c2 && push!(eqs_to_rm, ieq)
+    end
+
     # We need to handle unscalarized array variables
     for k in keys(subs)
         k, isarr = split_indexed_var(k)

From fe894ba8fe49474ad67c2039f70290c002dd1669 Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Fri, 15 May 2026 17:55:27 -0400
Subject: [PATCH 07/15] fixes and simplifications

---
 src/systems/alias_elimination.jl | 98 ++++++++++++++------------------
 test/reduction.jl                | 24 ++++++++
 2 files changed, 66 insertions(+), 56 deletions(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index b8a134dfe8..8ba1535587 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -49,7 +49,7 @@ maintains the invariant that `parent[v]` is always the current root and `parity[
 is the sign of `v` relative to that root so finds are O(1).
 
 Merges are smaller-into-larger using `members` (root → member list)
-giving O((V+E) log V) across all unions.
+giving O(E log V) across all unions.
 A contradiction (e.g. `x ~ y, x ~ -y`) zeros out every member's `parity` in the
 affected component; downstream code treats `parity[v] == 0` as "v is forced to 0".
 """
@@ -195,8 +195,8 @@ function find_perfect_aliases!(
 
     # After the scan above, every variable in any alias group has
     # `parent[v] == root` and `parity[v] ∈ {-1, 0, +1}`
-    # (sign relative to root, or 0 if the component is conflict-tainted). 
-    # The main rewriting loop below picks a target per non-conflict group 
+    # (sign relative to root, or 0 if the component is conflict-tainted).
+    # The main rewriting loop below picks a target per non-conflict group
     # and uses `s = parity[v] * parity[target]` (±1) as `v`'s sign relative to the target.
     # For Conflict groups: every non-irreducible gets substituted with the
     # symbolic `0`, and one alias eq per irreducible is rewritten to `irr ~ 0`.
@@ -210,43 +210,16 @@ function find_perfect_aliases!(
     zero_sym = Symbolics.COMMON_ZERO
     signed_sym = (s::Int8, x::SymbolicT) -> s == 1 ? x : (-1) * x
 
-    # Conflict-group alias equations: every variable is forced to `0`.
-    # Each irreducible must stay in unknowns, so claim one alias eq per irreducible
-    # and rewrite it eagerly to `v ~ 0` (stripping the bipartite edge to the other endpoint).
-    # A conflict group contains a cycle so has at least as many alias equations as vertices
-    # so there is always one eq per irreducible.
-    # Remaining alias equations become `0 ~ 0` after the substitution and are queued for removal.
-    let claimed = Set{Int}()
-        for (ieq, v1, v2, _) in candidate_eqs
-            parity[v1] == 0 || continue
-            v_pin = if is_irreducible_v(v1) && !(v1 in claimed)
-                v1
-            elseif is_irreducible_v(v2) && !(v2 in claimed)
-                v2
-            else
-                nothing
-            end
-            if v_pin === nothing
-                push!(eqs_to_rm, ieq)
-            else
-                push!(claimed, v_pin)
-                for u in copy(𝑠neighbors(graph, ieq))
-                    u == v_pin && continue
-                    Graphs.rem_edge!(graph, ieq, u)
-                end
-                eqs[ieq] = fullvars[v_pin] ~ zero_sym
-                original_eqs[ieq] = fullvars[v_pin] ~ zero_sym
-            end
-        end
-    end
-
     eqs_to_substitute = Int[]
+    irrs_by_root = Dict{Int, Vector{Int}}()
     for (root, group_vars) in members
         if parity[root] == 0 # is conflict group
+            # collect the irreducibles that need an equation forcing them to `0`.
+            irrs_by_root[root] = filter(is_irreducible_v, group_vars)
             for v in group_vars
-                # In conflict groups, ONLY irreducible vars survive as unknowns
-                # Irreducibles already had equation rewritten to `v ~ 0` above
-                # so just mark them present here
+                # In conflict groups, ONLY irreducible vars survive as unknowns.
+                # The per-equation cleanup below rewrites one alias eq per
+                # irreducible to `v ~ 0`; here we just mark them present.
                 if is_irreducible_v(v)
                     state.always_present[v] = true
                     continue
@@ -255,20 +228,16 @@ function find_perfect_aliases!(
                 subs[fullvars[v]] = zero_sym
                 push!(state.additional_observed, fullvars[v] ~ zero_sym)
 
-                for e in copy(𝑑neighbors(graph, v))
-                    push!(eqs_to_substitute, e)
-                    Graphs.rem_edge!(graph, e, v)
-                end
+                append!(eqs_to_substitute, 𝑑neighbors(graph, v))
+                BipartiteGraphs.set_neighbors!(BipartiteGraphs.invview(graph), v, ())
 
                 dv = var_to_diff[v]
                 while dv isa Int
                     push!(vars_to_rm, dv)
                     subs[fullvars[dv]] = zero_sym
 
-                    for e in copy(𝑑neighbors(graph, dv))
-                        push!(eqs_to_substitute, e)
-                        Graphs.rem_edge!(graph, e, dv)
-                    end
+                    append!(eqs_to_substitute, 𝑑neighbors(graph, dv))
+                    BipartiteGraphs.set_neighbors!(BipartiteGraphs.invview(graph), dv, ())
 
                     dv = var_to_diff[dv]
                 end
@@ -280,12 +249,10 @@ function find_perfect_aliases!(
         target = pick_alias_target(fullvars, group_vars, state_priorities, irreducibles, graph)
         group_target[root] = target
         target_p = parity[target]
-        state.always_present[target] = true
         for v in group_vars
-            v == target && continue
             # In consistent groups sticky vars (irreducible OR priority>0) other
             # than the target stay as unknowns.
-            if is_sticky(v)
+            if is_sticky(v) || v == target
                 state.always_present[v] = true
                 continue
             end
@@ -326,16 +293,35 @@ function find_perfect_aliases!(
         end
     end
 
-    # Consistent-group alias equations: queue for removal iff both endpoints
-    # collapse (equation becomes `0 ~ 0`). An equation with a non-target
-    # irreducible endpoint must stay so the remaining irreducible is still
-    # pinned to the target after substitution.
+    # Per-equation cleanup. Two cases:
+    #
+    # Conflict groups: every variable is forced to `0`.
+    # Take the first `length(irrs)` eq in the group to force the irrs to 0
+    # (overwriting the eq and replacing its graph row with the edge to `irr`). 
+    # Remaining conflict eqs become `0 ~ 0` and are queued for removal.
+    # Pinned eqs end up in `eqs_to_substitute` from the var-elim pass above
+    # but `subber` leaves `irr ~ 0` untouched (irreducibles aren't in `subs`).
+    #
+    # Consistent groups: queue for removal iff both endpoints collapse. An eq
+    # with a non-target irreducible endpoint must stay so the remaining
+    # irreducible is still pinned to the target after substitution.
     for (ieq, v1, v2, _) in candidate_eqs
-        parity[v1] == 0 && continue
-        target = group_target[parent[v1]]
-        c1 = is_sticky(v1) ? v1 : target
-        c2 = is_sticky(v2) ? v2 : target
-        c1 == c2 && push!(eqs_to_rm, ieq)
+        if parity[v1] == 0 # if conflict group
+            irrs = irrs_by_root[parent[v1]]
+            if isempty(irrs)
+                push!(eqs_to_rm, ieq)
+            else
+                v_pin = pop!(irrs) # irriducable var to pin to 0
+                BipartiteGraphs.set_neighbors!(graph, ieq, [v_pin])
+                eqs[ieq] = fullvars[v_pin] ~ zero_sym
+                original_eqs[ieq] = fullvars[v_pin] ~ zero_sym
+            end
+        else
+            target = group_target[parent[v1]]
+            c1 = is_sticky(v1) ? v1 : target
+            c2 = is_sticky(v2) ? v2 : target
+            c1 == c2 && push!(eqs_to_rm, ieq)
+        end
     end
 
     # We need to handle unscalarized array variables
diff --git a/test/reduction.jl b/test/reduction.jl
index e6a14db6e5..69210235ed 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -468,6 +468,30 @@ end
         @test isequal(only(eqs).lhs, D(z))
     end
 
+    @testset "conflict with 3 irreducibles, one only on a single eq" begin
+        @variables a(t) [irreducible = true]
+        @variables b(t) [irreducible = true]
+        @variables c(t) [irreducible = true]
+        @variables w(t)
+        # All three are irreducible, all in one conflict group via `a ~ b`,
+        # `a ~ c`, `a ~ -c`. Each irreducible needs its own `irr ~ 0` equation.
+        # `b` only appears on `e1 = a ~ b`. A previous greedy claiming
+        # algorithm (prefer v1, fall back to v2) would have stranded `b`:
+        # e1 → claim a, e2=(a,c) → claim c, e3=(a,c) → both claimed, drop. So
+        # b would survive as an unconstrained unknown. With the matching-free
+        # rewrite we now do (any irr ↔ any group eq), all three irreducibles
+        # get their own pinned equation.
+        @named sys = System([D(w) ~ w, a ~ b, a ~ c, a ~ -c], t)
+        state = TearingState(sys)
+        ModelingToolkit.eliminate_perfect_aliases!(state)
+        eqs = equations(state.sys)
+        for v in (a, b, c)
+            @test any(
+                eq -> isequal(eq.lhs, unwrap(v)) && iszero(Symbolics.value(eq.rhs)), eqs)
+        end
+        @test any(eq -> isequal(eq.lhs, D(w)), eqs)
+    end
+
     @testset "conflicting signs with irreducible: alias eq rewritten to `v ~ 0`" begin
         @variables x(t) [irreducible = true]
         @variables y(t) z(t)

From 40c67ee4039428abc449cda8b0cf2a168938305c Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Sat, 16 May 2026 10:27:19 -0400
Subject: [PATCH 08/15] Update src/systems/alias_elimination.jl

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 src/systems/alias_elimination.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 8ba1535587..bb6cf23001 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -311,7 +311,7 @@ function find_perfect_aliases!(
             if isempty(irrs)
                 push!(eqs_to_rm, ieq)
             else
-                v_pin = pop!(irrs) # irriducable var to pin to 0
+                v_pin = pop!(irrs) # irreducible var to pin to 0
                 BipartiteGraphs.set_neighbors!(graph, ieq, [v_pin])
                 eqs[ieq] = fullvars[v_pin] ~ zero_sym
                 original_eqs[ieq] = fullvars[v_pin] ~ zero_sym

From 171c3e1d70c9f92de072afbf2930cee340d73fcf Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Sat, 16 May 2026 10:27:27 -0400
Subject: [PATCH 09/15] Update src/systems/alias_elimination.jl

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 src/systems/alias_elimination.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index bb6cf23001..ba9f17de5d 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -157,7 +157,7 @@ function find_perfect_aliases!(
     parity = Dict{Int, Int8}()
     members = Dict{Int, Vector{Int}}()
     # Candidate alias equations `(ieq, v1_idx, v2_idx, edge_sign)`. 
-    # `edge_sign == +1` encodes `v1 ~ edge_sign*v2`.
+    # `edge_sign == ±1` encodes `v1 ~ edge_sign*v2`.
     # Removal is decided below once each group's target is known:
     # equations with a non-target irreducible endpoint must stay so
     # the remaining irreducibles are still constrained to the target.

From 4280b5baf4fee54bff978d10040be8bb78ad1f9f Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Mon, 18 May 2026 11:52:29 -0400
Subject: [PATCH 10/15] update eq incedence and test

---
 src/systems/alias_elimination.jl |  9 +++++++
 test/reduction.jl                | 46 ++++++++++++++++++++++++++++++++
 2 files changed, 55 insertions(+)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index ba9f17de5d..b36d9a7e55 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -339,6 +339,15 @@ function find_perfect_aliases!(
         # substitutes `x` to `[x[1], x[2]]`. The second substitutes `x[1]` and `x[2]`.
         eqs[e] = subber(subber(eqs[e]))
         original_eqs[e] = subber(subber(original_eqs[e]))
+        # Substitution can drop vars beyond the one substituted: zero
+        # substitution annihilates multiplicative cofactors (`v*w` with `v→0`
+        # also removes `w`); alias substitution can cancel the target
+        # (`v - target` with `v→target` leaves neither). Substitution never
+        # introduces new vars, so the new incidence is a subset of the old —
+        # prune the row to entries actually present in the simplified RHS.
+        new_vars = Symbolics.get_variables(eqs[e], fullvars)
+        new_row = filter(v_idx -> fullvars[v_idx] in new_vars, 𝑠neighbors(graph, e))
+        BipartiteGraphs.set_neighbors!(graph, e, new_row)
     end
 
     # After substitution, alias equations that connected a sticky non-target
diff --git a/test/reduction.jl b/test/reduction.jl
index 69210235ed..49d3f9b174 100644
--- a/test/reduction.jl
+++ b/test/reduction.jl
@@ -512,3 +512,49 @@ end
         @test any(eq -> isequal(eq.lhs, D(z)), eqs)
     end
 end
+
+@testset "Substitution rebuilds equation incidence" begin
+    # Calls `find_perfect_aliases!` directly so we can inspect the bipartite
+    # graph *before* `rm_eqs_vars!` renumbers it.
+    sneighbors = ModelingToolkit.BipartiteGraphs.𝑠neighbors
+
+    @testset "zero substitution drops multiplicatively annihilated var" begin
+        @variables x(t) y(t) a(t) b(t)
+        # `x ~ y` and `x ~ -y` is a sign conflict ⇒ x = y = 0. The non-alias
+        # eq `D(a) ~ x*b + a` simplifies to `D(a) ~ a` after `x → 0`, so `b`
+        # is annihilated and its edge to that eq must be dropped even though
+        # `b` itself wasn't substituted.
+        @named sys = System(
+            [D(a) ~ x * b + a, D(b) ~ b, x ~ y, x ~ -y], t)
+        state = TearingState(sys)
+
+        eqs_before = collect(ModelingToolkit.equations(state))
+        da_ieq = findfirst(eq -> isequal(eq.lhs, D(a)), eqs_before)
+        b_idx = findfirst(isequal(unwrap(b)), state.fullvars)
+        # Precondition: b is incident on the D(a) eq before the pass runs.
+        @test b_idx in sneighbors(state.structure.graph, da_ieq)
+
+        ModelingToolkit.find_perfect_aliases!(state, Int[], Int[])
+
+        @test !(b_idx in sneighbors(state.structure.graph, da_ieq))
+    end
+
+    @testset "alias substitution drops cancelled target var" begin
+        @variables x(t) y(t) w(t)
+        # `x ~ y` is a consistent alias. `x` gets eliminated in favor of `y`.
+        # The non-alias eq `D(w) ~ x - y + w` becomes `D(w) ~ w`, so `y` is no
+        # longer in the equation even though it was the alias *target*.
+        @named sys = System([D(w) ~ x - y + w, x ~ y], t;
+            state_priorities = [y => 10])
+        state = TearingState(sys)
+
+        eqs_before = collect(ModelingToolkit.equations(state))
+        dw_ieq = findfirst(eq -> isequal(eq.lhs, D(w)), eqs_before)
+        y_idx = findfirst(isequal(unwrap(y)), state.fullvars)
+        @test y_idx in sneighbors(state.structure.graph, dw_ieq)
+
+        ModelingToolkit.find_perfect_aliases!(state, Int[], Int[])
+
+        @test !(y_idx in sneighbors(state.structure.graph, dw_ieq))
+    end
+end

From ac7517fc92a4a2f0b02bf5b428a88aa7d3cc3661 Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Tue, 19 May 2026 09:00:41 -0400
Subject: [PATCH 11/15] Update src/systems/alias_elimination.jl

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 src/systems/alias_elimination.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index b36d9a7e55..29e3e7bc19 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -333,7 +333,7 @@ function find_perfect_aliases!(
     end
 
     ir = get_irstructure(sys)
-    subber = SU.IRSubstituter{false}(ir, subs)
+    subber = SU.IRSubstituter{true}(ir, subs)
     for e in eqs_to_substitute
         # Double substitute to handle unscalarized array variables. First one
         # substitutes `x` to `[x[1], x[2]]`. The second substitutes `x[1]` and `x[2]`.

From 04bac966a862f1ab099d7006dd3d5b9ff8a7335b Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Tue, 19 May 2026 09:01:02 -0400
Subject: [PATCH 12/15] Update src/systems/alias_elimination.jl

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 src/systems/alias_elimination.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index 29e3e7bc19..ce5bc62b6b 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -334,6 +334,7 @@ function find_perfect_aliases!(
 
     ir = get_irstructure(sys)
     subber = SU.IRSubstituter{true}(ir, subs)
+    new_vars = SU.IRStructureSearchBuffer(ir, Set{SymbolicT}())
     for e in eqs_to_substitute
         # Double substitute to handle unscalarized array variables. First one
         # substitutes `x` to `[x[1], x[2]]`. The second substitutes `x[1]` and `x[2]`.

From 0d6a00dc6ad383045b6aecb7fc0f67ec806493b1 Mon Sep 17 00:00:00 2001
From: Oscar Smith <oscardssmith@gmail.com>
Date: Wed, 20 May 2026 21:05:09 -0400
Subject: [PATCH 13/15] Update src/systems/alias_elimination.jl

Co-authored-by: Aayush Sabharwal <aayush.sabharwal@gmail.com>
---
 src/systems/alias_elimination.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/systems/alias_elimination.jl b/src/systems/alias_elimination.jl
index ce5bc62b6b..621bfe8379 100644
--- a/src/systems/alias_elimination.jl
+++ b/src/systems/alias_elimination.jl
@@ -346,7 +346,8 @@ function find_perfect_aliases!(
         # (`v - target` with `v→target` leaves neither). Substitution never
         # introduces new vars, so the new incidence is a subset of the old —
         # prune the row to entries actually present in the simplified RHS.
-        new_vars = Symbolics.get_variables(eqs[e], fullvars)
+        empty!(new_vars)
+        SU.search_variables!(new_vars, eqs[e])
         new_row = filter(v_idx -> fullvars[v_idx] in new_vars, 𝑠neighbors(graph, e))
         BipartiteGraphs.set_neighbors!(graph, e, new_row)
     end

From ae13991f27e90312eba39d358c8456e5e7b11f9c Mon Sep 17 00:00:00 2001
From: oscarddssmith <oscar.smith@juliahub.com>
Date: Thu, 21 May 2026 16:06:15 -0400
Subject: [PATCH 14/15] update gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 9193389ab8..bdef6a4e9d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ Manifest.toml
 docs/src/assets/Project.toml
 docs/src/assets/Manifest.toml
 .CondaPkg
+.claude/settings.local.json

From 60ad3fe6769bd98888b5efd37d27a14c0457672c Mon Sep 17 00:00:00 2001
From: oscarddssmith <oscar.smith@juliahub.com>
Date: Thu, 21 May 2026 16:10:31 -0400
Subject: [PATCH 15/15] bump SymbolicUtils compat

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

diff --git a/Project.toml b/Project.toml
index b0d1665bc0..cd1bf6503c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -118,7 +118,7 @@ StateSelection = "1.9.1"
 StaticArrays = "1.9.14"
 StochasticDiffEq = "6.82.0, 7"
 SymbolicIndexingInterface = "0.3.39"
-SymbolicUtils = "4.28"
+SymbolicUtils = "4.31"
 Symbolics = "7"
 UnPack = "0.1, 1.0"
 julia = "1.9"