Contract removed equations away

Author: YingboMa

src/structural_transformation/partial_state_selection.jl

@@ -248,6 +248,7 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
         end
     end
     if diff_va !== nothing
+        # differentiated alias
         n_dummys = length(dummy_derivatives)
         needed = count(x -> x isa Int, diff_to_eq) - n_dummys
         n = 0

src/structural_transformation/symbolics_tearing.jl

@@ -217,7 +217,7 @@ function check_diff_graph(var_to_diff, fullvars)
 end
 =#
 
-function tearing_reassemble(state::TearingState, var_eq_matching; simplify = false)
+function tearing_reassemble(state::TearingState, var_eq_matching, ag = nothing; simplify = false)
     @unpack fullvars, sys, structure = state
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure
 
@@ -557,7 +557,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
     end
     invvarsperm = [diff_vars;
                    setdiff!(setdiff(1:ndsts(graph), diff_vars_set),
-                            BitSet(solved_variables))]
+                            union!(BitSet(solved_variables), keys(ag)))]
     varsperm = zeros(Int, ndsts(graph))
     for (i, v) in enumerate(invvarsperm)
         varsperm[v] = i
@@ -632,6 +632,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify = fal
         isdiffeq(eq) || continue
         obs_sub[eq.lhs] = eq.rhs
     end
+    # TODO: compute the dependency correctly so that we don't have to do this
     obs = substitute.([oldobs; subeqs], (obs_sub,))
     @set! sys.observed = obs
     @set! state.sys = sys
@@ -688,7 +689,7 @@ end
 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, kwargs...)
+function dummy_derivative(sys, state = TearingState(sys), ag = nothing; simplify = false, kwargs...)
     jac = let state = state
         (eqs, vars) -> begin
             symeqs = EquationsView(state)[eqs]
@@ -710,6 +711,6 @@ function dummy_derivative(sys, state = TearingState(sys); simplify = false, kwar
             p
         end
     end
-    var_eq_matching = dummy_derivative_graph!(state, jac; state_priority, kwargs...)
-    tearing_reassemble(state, var_eq_matching; simplify = simplify)
+    var_eq_matching = dummy_derivative_graph!(state, jac, (ag, nothing); state_priority, kwargs...)
+    tearing_reassemble(state, var_eq_matching, ag; simplify = simplify)
 end

src/structural_transformation/utils.jl

@@ -46,12 +46,13 @@ end
 function check_consistency(state::TearingState, ag = nothing)
     fullvars = state.fullvars
     @unpack graph, var_to_diff = state.structure
-    #n_highest_vars = count(v -> length(outneighbors(var_to_diff, v)) == 0 && !isempty(𝑑neighbors(graph, v)),
-    #                       vertices(var_to_diff))
-    #n_highest_vars = count(v -> length(outneighbors(var_to_diff, v)) == 0 && !haskey(ag, v),
-    #                       vertices(var_to_diff))
-    n_highest_vars = count(v -> length(outneighbors(var_to_diff, v)) == 0,
-                           vertices(var_to_diff))
+    if ag === nothing
+        n_highest_vars = count(v -> length(outneighbors(var_to_diff, v)) == 0,
+                            vertices(var_to_diff))
+    else
+        n_highest_vars = count(v -> length(outneighbors(var_to_diff, v)) == 0 && !haskey(ag, v),
+                            vertices(var_to_diff))
+    end
     neqs = nsrcs(graph)
     is_balanced = n_highest_vars == neqs
 
@@ -73,7 +74,7 @@ function check_consistency(state::TearingState, ag = nothing)
     # details, check the equation (15) of the original paper.
     extended_graph = (@set graph.fadjlist = Vector{Int}[graph.fadjlist;
                                                         map(collect, edges(var_to_diff))])
-    extended_var_eq_matching = maximal_matching(extended_graph, eq->true, v->!haskey(ag, v))
+    extended_var_eq_matching = maximal_matching(extended_graph, eq->true, v->ag !== nothing || !haskey(ag, v))
 
     unassigned_var = []
     for (vj, eq) in enumerate(extended_var_eq_matching)

src/systems/abstractsystem.jl

@@ -1041,7 +1041,7 @@ function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false
     #has_io && markio!(state, io..., check = false)
     check_consistency(state, ag)
     #find_solvables!(state; kwargs...)
-    sys = dummy_derivative(sys, state; simplify)
+    sys = dummy_derivative(sys, state, ag; simplify)
     fullstates = [map(eq -> eq.lhs, observed(sys)); states(sys)]
     @set! sys.observed = topsort_equations(observed(sys), fullstates)
     invalidate_cache!(sys)

src/systems/alias_elimination.jl

@@ -167,7 +167,6 @@ function alias_elimination!(state::TearingState)
             diff_to_var[j] === nothing && push!(newstates, fullvars[j])
         end
     end
-    #=
     new_graph = BipartiteGraph(n_new_eqs, ndsts(graph))
     new_solvable_graph = BipartiteGraph(n_new_eqs, ndsts(graph))
     new_eq_to_diff = DiffGraph(n_new_eqs)
@@ -182,8 +181,8 @@ function alias_elimination!(state::TearingState)
     state.structure.solvable_graph = new_solvable_graph
     state.structure.eq_to_diff = new_eq_to_diff
     @show length(new_eq_to_diff), nsrcs(new_graph), nsrcs(new_solvable_graph), length(eqs)
-    =#
 
+    #=
     new_graph = BipartiteGraph(n_new_eqs, n_new_vars)
     new_solvable_graph = BipartiteGraph(n_new_eqs, n_new_vars)
     new_eq_to_diff = DiffGraph(n_new_eqs)
@@ -207,6 +206,7 @@ function alias_elimination!(state::TearingState)
     state.structure.eq_to_diff = complete(new_eq_to_diff)
     state.structure.var_to_diff = complete(new_var_to_diff)
     state.fullvars = new_fullvars
+    =#
 
     sys = state.sys
     @set! sys.eqs = eqs