Construct linear matrix and solvable graph at the same time

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
+                       AliasGraph, filter_kwargs, lower_varname, setio, SparseMatrixCLIL
 
 using ModelingToolkit.BipartiteGraphs
 import .BipartiteGraphs: invview, complete
@@ -43,7 +43,8 @@ using NonlinearSolve
 export tearing, partial_state_selection, dae_index_lowering, check_consistency
 export dummy_derivative
 export build_torn_function, build_observed_function, ODAEProblem
-export sorted_incidence_matrix, pantelides!, tearing_reassemble, find_solvables!
+export sorted_incidence_matrix, pantelides!, tearing_reassemble, find_solvables!,
+       linear_subsys_adjmat!
 export tearing_assignments, tearing_substitution
 export torn_system_jacobian_sparsity
 export full_equations

src/structural_transformation/symbolics_tearing.jl

@@ -62,9 +62,8 @@ function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int)
     eq_diff = eq_derivative_graph!(s, ieq)
 
     sys = ts.sys
-    D = Differential(get_iv(sys))
     eq = equations(ts)[ieq]
-    eq = ModelingToolkit.expand_derivatives(0 ~ D(eq.rhs - eq.lhs))
+    eq = 0 ~ ModelingToolkit.derivative(eq.rhs - eq.lhs, get_iv(sys))
     push!(equations(ts), eq)
     # Analyze the new equation and update the graph/solvable_graph
     # First, copy the previous incidence and add the derivative terms.

src/structural_transformation/utils.jl

@@ -160,20 +160,24 @@ end
 ### Structural and symbolic utilities
 ###
 
-function find_eq_solvables!(state::TearingState, ieq; may_be_zero = false,
+function find_eq_solvables!(state::TearingState, ieq, to_rm = Int[], coeffs = nothing;
+                            may_be_zero = false,
                             allow_symbolic = false, allow_parameter = true, kwargs...)
     fullvars = state.fullvars
     @unpack graph, solvable_graph = state.structure
     eq = equations(state)[ieq]
     term = value(eq.rhs - eq.lhs)
-    to_rm = Int[]
+    all_int_vars = true
+    coeffs === nothing || empty!(coeffs)
+    empty!(to_rm)
     for j in 𝑠neighbors(graph, ieq)
         var = fullvars[j]
-        isirreducible(var) && continue
+        isirreducible(var) && (all_int_vars = false; continue)
         a, b, islinear = linear_expansion(term, var)
-        a = unwrap(a)
-        islinear || continue
+        a, b = unwrap(a), unwrap(b)
+        islinear || (all_int_vars = false; continue)
         if a isa Symbolic
+            all_int_vars = false
             if !allow_symbolic
                 if allow_parameter
                     all(ModelingToolkit.isparameter, vars(a)) || continue
@@ -184,7 +188,18 @@ function find_eq_solvables!(state::TearingState, ieq; may_be_zero = false,
             add_edge!(solvable_graph, ieq, j)
             continue
         end
-        (a isa Number) || continue
+        if !(a isa Number)
+            all_int_vars = false
+            continue
+        end
+        # When the expression is linear with numeric `a`, then we can safely
+        # only consider `b` for the following iterations.
+        term = b
+        if isone(abs(a))
+            coeffs === nothing || push!(coeffs, convert(Int, a))
+        else
+            all_int_vars = false
+        end
         if a != 0
             add_edge!(solvable_graph, ieq, j)
         else
@@ -198,19 +213,54 @@ function find_eq_solvables!(state::TearingState, ieq; may_be_zero = false,
     for j in to_rm
         rem_edge!(graph, ieq, j)
     end
+    all_int_vars, term
 end
 
 function find_solvables!(state::TearingState; kwargs...)
     @assert state.structure.solvable_graph === nothing
     eqs = equations(state)
     graph = state.structure.graph
     state.structure.solvable_graph = BipartiteGraph(nsrcs(graph), ndsts(graph))
+    to_rm = Int[]
     for ieq in 1:length(eqs)
-        find_eq_solvables!(state, ieq; kwargs...)
+        find_eq_solvables!(state, ieq, to_rm; kwargs...)
     end
     return nothing
 end
 
+function linear_subsys_adjmat!(state::TransformationState)
+    graph = state.structure.graph
+    if state.structure.solvable_graph === nothing
+        state.structure.solvable_graph = BipartiteGraph(nsrcs(graph), ndsts(graph))
+    end
+    linear_equations = Int[]
+    eqs = equations(state.sys)
+    eadj = Vector{Int}[]
+    cadj = Vector{Int}[]
+    coeffs = Int[]
+    to_rm = Int[]
+    for i in eachindex(eqs)
+        all_int_vars, rhs = find_eq_solvables!(state, i, to_rm, coeffs)
+
+        # Check if all states in the equation is both linear and homogeneous,
+        # i.e. it is in the form of
+        #
+        #       ``∑ c_i * v_i = 0``,
+        #
+        # where ``c_i`` ∈ ℤ and ``v_i`` denotes states.
+        if all_int_vars && Symbolics._iszero(rhs)
+            push!(linear_equations, i)
+            push!(eadj, copy(𝑠neighbors(graph, i)))
+            push!(cadj, copy(coeffs))
+        end
+    end
+
+    mm = SparseMatrixCLIL(nsrcs(graph),
+                          ndsts(graph),
+                          linear_equations, eadj, cadj)
+    return mm
+end
+
 highest_order_variable_mask(ts) =
     let v2d = ts.structure.var_to_diff
         v -> isempty(outneighbors(v2d, v))

src/systems/alias_elimination.jl

@@ -5,7 +5,7 @@ using Graphs.Experimental.Traversals
 const KEEP = typemin(Int)
 
 function alias_eliminate_graph!(state::TransformationState)
-    mm = linear_subsys_adjmat(state)
+    mm = linear_subsys_adjmat!(state)
     if size(mm, 1) == 0
         ag = AliasGraph(ndsts(state.structure.graph))
         return ag, mm, ag, mm, BitSet() # No linear subsystems

src/systems/systemstructure.jl

@@ -351,56 +351,6 @@ function TearingState(sys; quick_cancel = false, check = true)
                                         complete(graph), nothing), Any[])
 end
 
-function linear_subsys_adjmat(state::TransformationState)
-    fullvars = state.fullvars
-    graph = state.structure.graph
-    is_linear_equations = falses(nsrcs(graph))
-    eqs = equations(state.sys)
-    eadj = Vector{Int}[]
-    cadj = Vector{Int}[]
-    coeffs = Int[]
-    for (i, eq) in enumerate(eqs)
-        empty!(coeffs)
-        linear_term = 0
-        all_int_vars = true
-
-        term = value(eq.rhs - eq.lhs)
-        for j in 𝑠neighbors(graph, i)
-            var = fullvars[j]
-            a, b, islinear = linear_expansion(term, var)
-            a = unwrap(a)
-            if islinear && !(a isa Symbolic) && a isa Number && !isirreducible(var)
-                if a == 1 || a == -1
-                    a = convert(Integer, a)
-                    linear_term += a * var
-                    push!(coeffs, a)
-                else
-                    all_int_vars = false
-                end
-            end
-        end
-
-        # Check if all states in the equation is both linear and homogeneous,
-        # i.e. it is in the form of
-        #
-        #       ``∑ c_i * v_i = 0``,
-        #
-        # where ``c_i`` ∈ ℤ and ``v_i`` denotes states.
-        if all_int_vars && isequal(linear_term, term)
-            is_linear_equations[i] = true
-            push!(eadj, copy(𝑠neighbors(graph, i)))
-            push!(cadj, copy(coeffs))
-        else
-            is_linear_equations[i] = false
-        end
-    end
-
-    linear_equations = findall(is_linear_equations)
-    return SparseMatrixCLIL(nsrcs(graph),
-                            ndsts(graph),
-                            linear_equations, eadj, cadj)
-end
-
 using .BipartiteGraphs: Label, BipartiteAdjacencyList
 struct SystemStructurePrintMatrix <:
        AbstractMatrix{Union{Label, Int, BipartiteAdjacencyList}}