Merge pull request #3041 from SciML/myb/uneven

Tolerate empty equations and uneven matching

Author: YingboMa

src/bipartite_graph.jl

@@ -69,21 +69,29 @@ function Base.setindex!(m::Matching{U}, v::Union{Integer, U}, i::Integer) where
 
         # To maintain the invariant that `m.inv_match[m.match[i]] == i`, we need
         # to unassign the matching at `m.inv_match[v]` if it exists.
-        if v isa Int && (iv = m.inv_match[v]) isa Int
+        if v isa Int && 1 <= v <= length(m.inv_match) && (iv = m.inv_match[v]) isa Int
             m.match[iv] = unassigned
         end
         if isa(oldv, Int)
             @assert m.inv_match[oldv] == i
             m.inv_match[oldv] = unassigned
         end
-        isa(v, Int) && (m.inv_match[v] = i)
+        if isa(v, Int)
+            for vv in (length(m.inv_match) + 1):v
+                push!(m.inv_match, unassigned)
+            end
+            m.inv_match[v] = i
+        end
     end
     return m.match[i] = v
 end
 
 function Base.push!(m::Matching, v)
     push!(m.match, v)
     if v isa Integer && m.inv_match !== nothing
+        for vv in (length(m.inv_match) + 1):v
+            push!(m.inv_match, unassigned)
+        end
         m.inv_match[v] = length(m.match)
     end
 end

src/structural_transformation/partial_state_selection.jl

@@ -214,6 +214,7 @@ function dummy_derivative_graph!(
     eqcolor = falses(nsrcs(graph))
     dummy_derivatives = Int[]
     col_order = Int[]
+    neqs = nsrcs(graph)
     nvars = ndsts(graph)
     eqs = Int[]
     vars = Int[]
@@ -236,7 +237,7 @@ function dummy_derivative_graph!(
         end
         isempty(eqs) && continue
 
-        rank_matching = Matching(nvars)
+        rank_matching = Matching(max(nvars, neqs))
         isfirst = true
         if jac === nothing
             J = nothing
@@ -389,11 +390,13 @@ function tearing_with_dummy_derivatives(structure, dummy_derivatives)
         Base.Fix1(isdiffed, (structure, dummy_derivatives)),
         Union{Unassigned, SelectedState};
         varfilter = Base.Fix1(getindex, can_eliminate))
+
     for v in 𝑑vertices(structure.graph)
         is_present(structure, v) || continue
         dv = var_to_diff[v]
         (dv === nothing || !is_some_diff(structure, dummy_derivatives, dv)) && continue
         var_eq_matching[v] = SelectedState()
     end
+
     return var_eq_matching, full_var_eq_matching, var_sccs, can_eliminate
 end

src/structural_transformation/utils.jl

@@ -91,8 +91,10 @@ function check_consistency(state::TransformationState, orig_inputs)
                                                         map(collect, edges(var_to_diff))])
     extended_var_eq_matching = maximal_matching(extended_graph)
 
+    nvars = ndsts(graph)
     unassigned_var = []
     for (vj, eq) in enumerate(extended_var_eq_matching)
+        vj > nvars && break
         if eq === unassigned && !isempty(𝑑neighbors(graph, vj))
             push!(unassigned_var, fullvars[vj])
         end