Merge pull request #1701 from SciML/myb/alias

Simplify alias elimination logic

Author: YingboMa

src/systems/alias_elimination.jl

@@ -450,20 +450,28 @@ count_nonzeros(a::AbstractArray) = count(!iszero, a)
 # Here we have a guarantee that they won't, so we can make this identification
 count_nonzeros(a::SparseVector) = nnz(a)
 
-function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL, only_algebraic = true,
-                      irreducibles = ())
+function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL, irreducibles = ())
     mm = copy(mm_orig)
     is_linear_equations = falses(size(AsSubMatrix(mm_orig), 1))
     for e in mm_orig.nzrows
         is_linear_equations[e] = true
     end
 
-    is_not_potential_state = isnothing.(var_to_diff)
+    # If linear highest differentiated variables cannot be assigned to a pivot,
+    # then we can set it to zero. We use `rank1` to track this.
+    #
+    # We only use alias graph to record reducible variables. We use `rank2` to
+    # track this.
+    #
+    # For all the other variables, we can update the original system with
+    # Bareiss'ed coefficients as Gaussian elimination is nullspace perserving
+    # and we are only working on linear homogeneous subsystem.
+    is_linear_variables = isnothing.(var_to_diff)
+    is_reducible = trues(length(var_to_diff))
     for v in irreducibles
-        is_not_potential_state[v] = false
+        is_linear_variables[v] = false
+        is_reducible[v] = false
     end
-    is_linear_variables = only_algebraic ? copy(is_not_potential_state) :
-                          is_not_potential_state
     for i in 𝑠vertices(graph)
         is_linear_equations[i] && continue
         for j in 𝑠neighbors(graph, i)
@@ -481,12 +489,10 @@ function aag_bareiss!(graph, var_to_diff, mm_orig::SparseMatrixCLIL, only_algebr
                 r !== nothing && return r
                 rank1 = k - 1
             end
-            if only_algebraic
-                if rank2 === nothing
-                    r = find_masked_pivot(is_not_potential_state, M, k)
-                    r !== nothing && return r
-                    rank2 = k - 1
-                end
+            if rank2 === nothing
+                r = find_masked_pivot(is_reducible, M, k)
+                r !== nothing && return r
+                rank2 = k - 1
             end
             # TODO: It would be better to sort the variables by
             # derivative order here to enable more elimination
@@ -524,25 +530,9 @@ function lss(mm, pivots, ag)
     end
 end
 
-function simple_aliases!(ag, graph, var_to_diff, mm_orig, only_algebraic, irreducibles = ())
-    # Let `m = the number of linear equations` and `n = the number of
-    # variables`.
-    #
-    # `do_bareiss` conceptually gives us this system:
-    # rank1 | [ M₁₁  M₁₂ | M₁₃ ]   [v₁] = [0]
-    # rank2 | [ 0    M₂₂ | M₂₃ ] P [v₂] = [0]
-    # -------------------|-------------------
-    #         [ 0    0   | 0   ]   [v₃] = [0]
-
-    # Where `v₁` are the purely linear algebraic variables (i.e. those that only
-    # appear in linear algebraic equations), `v₂` are the variables that may be
-    # potentially solved by the linear system, and `v₃` are the variables that
-    # contribute to the equations, but are not solved by the linear system. Note
-    # that the complete system may be larger than the linear subsystem and
-    # include variables that do not appear here.
+function simple_aliases!(ag, graph, var_to_diff, mm_orig, irreducibles = ())
     mm, solvable_variables, (rank1, rank2, rank3, pivots) = aag_bareiss!(graph, var_to_diff,
                                                                          mm_orig,
-                                                                         only_algebraic,
                                                                          irreducibles)
 
     # Step 2: Simplify the system using the Bareiss factorization
@@ -585,7 +575,7 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     #
     nvars = ndsts(graph)
     ag = AliasGraph(nvars)
-    mm = simple_aliases!(ag, graph, var_to_diff, mm_orig, false)
+    mm = simple_aliases!(ag, graph, var_to_diff, mm_orig)
 
     # Step 3: Handle differentiated variables
     # At this point, `var_to_diff` and `ag` form a tree structure like the
@@ -694,7 +684,7 @@ function alias_eliminate_graph!(graph, var_to_diff, mm_orig::SparseMatrixCLIL)
     if !isempty(irreducibles)
         ag = newag
         mm_orig2 = isempty(ag) ? mm_orig : reduce!(copy(mm_orig), ag)
-        mm = simple_aliases!(ag, graph, var_to_diff, mm_orig2, true, irreducibles)
+        mm = simple_aliases!(ag, graph, var_to_diff, mm_orig2, irreducibles)
     end
 
     # for (v, (c, a)) in ag

test/components.jl

@@ -19,6 +19,8 @@ function check_contract(sys)
 end
 
 function check_rc_sol(sol)
+    rpi = sol[rc_model.resistor.p.i]
+    @test any(!isequal(rpi[1]), rpi) # test that we don't have a constant system
     @test sol[rc_model.resistor.p.i] == sol[resistor.p.i] == sol[capacitor.p.i]
     @test sol[rc_model.resistor.n.i] == sol[resistor.n.i] == -sol[capacitor.p.i]
     @test sol[rc_model.capacitor.n.i] == sol[capacitor.n.i] == -sol[capacitor.p.i]
@@ -31,8 +33,9 @@ end
 include("../examples/rc_model.jl")
 
 @test ModelingToolkit.n_extra_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 length(equations(sys)) == 1
 check_contract(sys)
 @test !isempty(ModelingToolkit.defaults(sys))
 u0 = [capacitor.v => 0.0
@@ -141,6 +144,7 @@ sol = solve(prob, Tsit5())
 
 include("../examples/serial_inductor.jl")
 sys = structural_simplify(ll_model)
+@test length(equations(sys)) == 2
 check_contract(sys)
 u0 = states(sys) .=> 0
 @test_nowarn ODEProblem(sys, u0, (0, 10.0))