Remove dead code

Author: YingboMa

src/structural_transformation/StructuralTransformations.jl

@@ -22,7 +22,7 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di
                        get_postprocess_fbody, vars!,
                        IncrementalCycleTracker, add_edge_checked!, topological_sort,
                        invalidate_cache!, Substitutions, get_or_construct_tearing_state,
-                       AliasGraph, filter_kwargs, lower_varname, setio, SparseMatrixCLIL,
+                       filter_kwargs, lower_varname, setio, SparseMatrixCLIL,
                        fast_substitute, get_fullvars, has_equations
 
 using ModelingToolkit.BipartiteGraphs

src/structural_transformation/pantelides.jl

@@ -70,34 +70,20 @@ function pantelides_reassemble(state::TearingState, var_eq_matching)
 end
 
 """
-    computed_highest_diff_variables(var_to_diff, ag)
+    computed_highest_diff_variables(var_to_diff)
 
 Computes which variables are the "highest-differentiated" for purposes of
 pantelides. Ordinarily this is relatively straightforward. However, in our
-case, there are two complicating  conditions:
+case, there is one complicating condition:
 
- 1. We allow variables in the structure graph that don't appear in the
+    We allow variables in the structure graph that don't appear in the
     system at all. What we are interested in is the highest-differentiated
     variable that actually appears in the system.
 
- 2. We have an alias graph. The alias graph implicitly contributes an
-    alias equation, so it doesn't actually whitelist any additional variables,
-    but it may change which variable is considered the highest differentiated one.
-    Consider the following situation:
-
-    Vars: x, y
-    Eqs: 0 = f(x)
-    Alias: ẋ = ẏ
-
-    In the absence of the alias, we would consider `x` to be the highest
-    differentiated variable. However, because of the alias (and because there
-    is no alias for `x=y`), we actually need to take `ẋ` as the highest
-    differentiated variable.
-
 This function takes care of these complications are returns a boolean array
 for every variable, indicating whether it is considered "highest-differentiated".
 """
-function computed_highest_diff_variables(structure, ag::Union{AliasGraph, Nothing})
+function computed_highest_diff_variables(structure)
     @unpack graph, var_to_diff = structure
 
     nvars = length(var_to_diff)
@@ -107,56 +93,12 @@ function computed_highest_diff_variables(structure, ag::Union{AliasGraph, Nothin
             # This variable is structurally highest-differentiated, but may not actually appear in the
             # system (complication 1 above). Ascend the differentiation graph to find the highest
             # differentiated variable that does appear in the system or the alias graph).
-            while isempty(𝑑neighbors(graph, var)) && (ag === nothing || !haskey(ag, var))
+            while isempty(𝑑neighbors(graph, var))
                 var′ = invview(var_to_diff)[var]
                 var′ === nothing && break
                 var = var′
             end
-            # If we don't have an alias graph, we are done. If we do have an alias graph, we may
-            # have to keep going along the stem, for as long as our differentiation path
-            # matches that of the stem (see complication 2 above). Note that we may end up
-            # whitelisting multiple differentiation levels of the stem here from different
-            # starting points that all map to the same stem. We clean that up in a post-processing
-            # pass below.
-            if ag !== nothing && haskey(ag, var)
-                (_, stem) = ag[var]
-                stem == 0 && continue
-                # If we have a self-loop in the stem, we could have the
-                # var′ also alias to the original stem. In that case, the
-                # derivative of the stem is highest differentiated, because of the loop
-                loop_found = false
-                var′ = invview(var_to_diff)[var]
-                while var′ !== nothing
-                    if var′ == stem || (haskey(ag, var′) && ag[var′][2] == stem)
-                        dstem = var_to_diff[stem]
-                        @assert dstem !== nothing
-                        varwhitelist[dstem] = true
-                        loop_found = true
-                        break
-                    end
-                    var′ = invview(var_to_diff)[var′]
-                end
-                loop_found && continue
-                # Ascend the stem
-                while isempty(𝑑neighbors(graph, var))
-                    var′ = invview(var_to_diff)[var]
-                    var′ === nothing && break
-                    loop_found = false
-                    cvar = var′
-                    # Invariant from alias elimination: Stem is chosen to have
-                    # the highest differentiation order.
-                    stem′ = invview(var_to_diff)[stem]
-                    @assert stem′ !== nothing
-                    if !haskey(ag, var′) || (ag[var′][2] != stem′)
-                        varwhitelist[stem] = true
-                        break
-                    end
-                    stem = stem′
-                    var = var′
-                end
-            else
-                varwhitelist[var] = true
-            end
+            varwhitelist[var] = true
         end
     end
 
@@ -181,8 +123,7 @@ end
 
 Perform Pantelides algorithm.
 """
-function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing} = nothing;
-                     finalize = true, maxiters = 8000)
+function pantelides!(state::TransformationState; finalize = true, maxiters = 8000)
     @unpack graph, solvable_graph, var_to_diff, eq_to_diff = state.structure
     neqs = nsrcs(graph)
     nvars = nv(var_to_diff)
@@ -193,7 +134,7 @@ function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing}
     nnonemptyeqs = count(eq -> !isempty(𝑠neighbors(graph, eq)) && eq_to_diff[eq] === nothing,
                          1:neqs′)
 
-    varwhitelist = computed_highest_diff_variables(state.structure, ag)
+    varwhitelist = computed_highest_diff_variables(state.structure)
 
     if nnonemptyeqs > count(varwhitelist)
         throw(InvalidSystemException("System is structurally singular"))

src/structural_transformation/partial_state_selection.jl

@@ -1,7 +1,6 @@
-function partial_state_selection_graph!(state::TransformationState;
-                                        ag::Union{AliasGraph, Nothing} = nothing)
+function partial_state_selection_graph!(state::TransformationState)
     find_solvables!(state; allow_symbolic = true)
-    var_eq_matching = complete(pantelides!(state, ag))
+    var_eq_matching = complete(pantelides!(state))
     complete!(state.structure)
     partial_state_selection_graph!(state.structure, var_eq_matching)
 end
@@ -150,14 +149,12 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
     var_eq_matching
 end
 
-function dummy_derivative_graph!(state::TransformationState, jac = nothing,
-                                 (ag, diff_va) = (nothing, nothing);
+function dummy_derivative_graph!(state::TransformationState, jac = nothing;
                                  state_priority = nothing, kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
-    var_eq_matching = complete(pantelides!(state, ag))
-    dummy_derivative_graph!(state.structure, var_eq_matching, jac, (ag, diff_va),
-                            state_priority)
+    var_eq_matching = complete(pantelides!(state))
+    dummy_derivative_graph!(state.structure, var_eq_matching, jac, state_priority)
 end
 
 function compute_diff_level(diff_to_x)
@@ -178,7 +175,7 @@ function compute_diff_level(diff_to_x)
 end
 
 function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac,
-                                 (ag, diff_va), state_priority)
+                                 state_priority)
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -249,61 +246,14 @@ 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
-        for v in diff_va
-            c, a = ag[v]
-            n += 1
-            push!(dummy_derivatives, iszero(c) ? v : a)
-            needed == n && break
-            continue
-        end
-    end
-
     if (n_diff_eqs = count(!isnothing, diff_to_eq)) !=
        (n_dummys = length(dummy_derivatives))
         @warn "The number of dummy derivatives ($n_dummys) does not match the number of differentiated equations ($n_diff_eqs)."
     end
     dummy_derivatives_set = BitSet(dummy_derivatives)
 
-    irreducible_set = BitSet()
-    if ag !== nothing
-        function isreducible(x)
-            # `k` is reducible if all lower differentiated variables are.
-            isred = true
-            while isred
-                if x in dummy_derivatives_set
-                    break
-                end
-                x = diff_to_var[x]
-                x === nothing && break
-                # We deliberately do not check `isempty(𝑑neighbors(graph, x))`
-                # because when `D(x)` appears in the alias graph, and `x`
-                # doesn't appear in any equations nor in the alias graph, `D(x)`
-                # is not reducible. Consider the system `D(x) ~ 0`.
-                if !haskey(ag, x)
-                    isred = false
-                end
-            end
-            isred
-        end
-        for (k, (c, v)) in ag
-            isreducible(k) || push!(irreducible_set, k)
-            iszero(c) && continue
-            isempty(𝑑neighbors(graph, v)) || push!(irreducible_set, v)
-        end
-    end
-
-    is_not_present_non_rec = let graph = graph, irreducible_set = irreducible_set
-        v -> begin
-            not_in_eqs = isempty(𝑑neighbors(graph, v))
-            ag === nothing && return not_in_eqs
-            isirreducible = v in irreducible_set
-            return not_in_eqs && !isirreducible
-        end
+    is_not_present_non_rec = let graph = graph
+        v -> isempty(𝑑neighbors(graph, v))
     end
 
     is_not_present = let var_to_diff = var_to_diff
@@ -330,11 +280,8 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
     # We can eliminate variables that are not a selected state (differential
     # variables). Selected states are differentiated variables that are not
     # dummy derivatives.
-    can_eliminate = let var_to_diff = var_to_diff, ag = ag
+    can_eliminate = let var_to_diff = var_to_diff
         v -> begin
-            if ag !== nothing
-                haskey(ag, v) && return false
-            end
             dv = var_to_diff[v]
             dv === nothing && return true
             is_some_diff(dv) || return true

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, ag = nothing;
+function tearing_reassemble(state::TearingState, var_eq_matching;
                             simplify = false, mm = nothing)
     @unpack fullvars, sys, structure = state
     @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure
@@ -494,7 +494,6 @@ function tearing_reassemble(state::TearingState, var_eq_matching, ag = nothing;
         error("Tearing internal error: lowering DAE into semi-implicit ODE failed!")
     end
     solved_variables_set = BitSet(solved_variables)
-    ag === nothing || union!(solved_variables_set, keys(ag))
     invvarsperm = [diff_vars;
                    setdiff!(setdiff(1:ndsts(graph), diff_vars_set),
                             solved_variables_set)]
@@ -608,7 +607,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), ag = nothing; simplify = false,
+function dummy_derivative(sys, state = TearingState(sys); simplify = false,
                           mm = nothing, kwargs...)
     jac = let state = state
         (eqs, vars) -> begin
@@ -631,7 +630,7 @@ function dummy_derivative(sys, state = TearingState(sys), ag = nothing; simplify
             p
         end
     end
-    var_eq_matching = dummy_derivative_graph!(state, jac, (ag, nothing); state_priority,
+    var_eq_matching = dummy_derivative_graph!(state, jac; state_priority,
                                               kwargs...)
-    tearing_reassemble(state, var_eq_matching, ag; simplify, mm)
+    tearing_reassemble(state, var_eq_matching; simplify, mm)
 end

src/structural_transformation/utils.jl

@@ -58,11 +58,11 @@ end
 ###
 ### Structural check
 ###
-function check_consistency(state::TransformationState, ag, orig_inputs)
+function check_consistency(state::TransformationState, orig_inputs)
     fullvars = get_fullvars(state)
     neqs = n_concrete_eqs(state)
     @unpack graph, var_to_diff = state.structure
-    highest_vars = computed_highest_diff_variables(complete!(state.structure), ag)
+    highest_vars = computed_highest_diff_variables(complete!(state.structure))
     n_highest_vars = 0
     for (v, h) in enumerate(highest_vars)
         h || continue
@@ -89,13 +89,11 @@ function check_consistency(state::TransformationState, ag, orig_inputs)
     # 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 -> ag === nothing || !haskey(ag, v))
+    extended_var_eq_matching = maximal_matching(extended_graph)
 
     unassigned_var = []
     for (vj, eq) in enumerate(extended_var_eq_matching)
-        if eq === unassigned && (ag === nothing || !haskey(ag, vj)) &&
-           !isempty(𝑑neighbors(graph, vj))
+        if eq === unassigned && !isempty(𝑑neighbors(graph, vj))
             push!(unassigned_var, fullvars[vj])
         end
     end