Merge pull request #809 from SciML/myb/alias

WIP: Alias elimination with singularity removal

Author: YingboMa

src/ModelingToolkit.jl

@@ -248,7 +248,7 @@ include("systems/dependency_graphs.jl")
 include("systems/systemstructure.jl")
 using .SystemStructures
 
-include("systems/reduction.jl")
+include("systems/alias_elimination.jl")
 
 include("latexify_recipes.jl")
 include("build_function.jl")
@@ -266,6 +266,7 @@ export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem
 export ControlSystem
+export alias_elimination, flatten
 export ode_order_lowering, liouville_transform
 export runge_kutta_discretize
 export PDESystem

src/bipartite_graph.jl

@@ -66,11 +66,13 @@ badjlist = [[1,2,5,6],[3,4,6]]
 bg = BipartiteGraph(7, fadjlist, badjlist)
 ```
 """
-mutable struct BipartiteGraph{I<:Integer} <: LightGraphs.AbstractGraph{I}
+mutable struct BipartiteGraph{I<:Integer,M} <: LightGraphs.AbstractGraph{I}
     ne::Int
     fadjlist::Vector{Vector{I}} # `fadjlist[src] => dsts`
     badjlist::Vector{Vector{I}} # `badjlist[dst] => srcs`
+    metadata::M
 end
+BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::AbstractVector) = BipartiteGraph(ne, fadj, badj, nothing)
 
 """
 ```julia
@@ -91,10 +93,22 @@ $(SIGNATURES)
 
 Build an empty `BipartiteGraph` with `nsrcs` sources and `ndsts` destinations.
 """
-BipartiteGraph(nsrcs::T, ndsts::T) where T = BipartiteGraph(0, map(_->T[], 1:nsrcs), map(_->T[], 1:ndsts))
+function BipartiteGraph(nsrcs::T, ndsts::T; metadata=nothing) where T
+    fadjlist = map(_->T[], 1:nsrcs)
+    badjlist = map(_->T[], 1:ndsts)
+    BipartiteGraph(0, fadjlist, badjlist, metadata)
+end
 
-Base.eltype(::Type{BipartiteGraph{I}}) where I = I
-Base.empty!(g::BipartiteGraph) = (foreach(empty!, g.fadjlist); foreach(empty!, g.badjlist); g.ne = 0; g)
+Base.eltype(::Type{<:BipartiteGraph{I}}) where I = I
+function Base.empty!(g::BipartiteGraph)
+    foreach(empty!, g.fadjlist)
+    foreach(empty!, g.badjlist)
+    g.ne = 0
+    if g.metadata !== nothing
+        foreach(empty!, g.metadata)
+    end
+    g
+end
 Base.length(::BipartiteGraph) = error("length is not well defined! Use `ne` or `nv`.")
 
 if isdefined(LightGraphs, :has_contiguous_vertices)
@@ -106,8 +120,8 @@ LightGraphs.vertices(g::BipartiteGraph) = (𝑠vertices(g), 𝑑vertices(g))
 𝑑vertices(g::BipartiteGraph) = axes(g.badjlist, 1)
 has_𝑠vertex(g::BipartiteGraph, v::Integer) = v in 𝑠vertices(g)
 has_𝑑vertex(g::BipartiteGraph, v::Integer) = v in 𝑑vertices(g)
-𝑠neighbors(g::BipartiteGraph, i::Integer) = g.fadjlist[i]
-𝑑neighbors(g::BipartiteGraph, i::Integer) = g.badjlist[i]
+𝑠neighbors(g::BipartiteGraph, i::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.fadjlist[i], g.metadata[i]) : g.fadjlist[i]
+𝑑neighbors(g::BipartiteGraph, j::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j]
 LightGraphs.ne(g::BipartiteGraph) = g.ne
 LightGraphs.nv(g::BipartiteGraph) = sum(length, vertices(g))
 LightGraphs.edgetype(g::BipartiteGraph{I}) where I = BipartiteEdge{I}
@@ -124,8 +138,12 @@ end
 ###
 ### Populate
 ###
-LightGraphs.add_edge!(g::BipartiteGraph, i::Integer, j::Integer) = add_edge!(g, BipartiteEdge(i, j))
-function LightGraphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge)
+struct NoMetadata
+end
+const NO_METADATA = NoMetadata()
+
+LightGraphs.add_edge!(g::BipartiteGraph, i::Integer, j::Integer, md=NO_METADATA) = add_edge!(g, BipartiteEdge(i, j), md)
+function LightGraphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md=NO_METADATA)
     @unpack fadjlist, badjlist = g
     verts = vertices(g)
     s, d = src(edge), dst(edge)
@@ -134,6 +152,9 @@ function LightGraphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge)
     index = searchsortedfirst(list, d)
     @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
     insert!(list, index, d)
+    if md !== NO_METADATA
+        insert!(g.metadata[s], index, md)
+    end
 
     g.ne += 1
     @inbounds list = badjlist[d]
@@ -170,7 +191,7 @@ Base.length(it::BipartiteEdgeIter{ALL}) = 2ne(it.g)
 
 Base.eltype(it::BipartiteEdgeIter) = edgetype(it.g)
 
-function Base.iterate(it::BipartiteEdgeIter{SRC,BipartiteGraph{T}}, state=(1, 1, SRC)) where T
+function Base.iterate(it::BipartiteEdgeIter{SRC,<:BipartiteGraph{T}}, state=(1, 1, SRC)) where T
     @unpack g = it
     neqs = nsrcs(g)
     neqs == 0 && return nothing
@@ -191,7 +212,7 @@ function Base.iterate(it::BipartiteEdgeIter{SRC,BipartiteGraph{T}}, state=(1, 1,
     return nothing
 end
 
-function Base.iterate(it::BipartiteEdgeIter{DST,BipartiteGraph{T}}, state=(1, 1, DST)) where T
+function Base.iterate(it::BipartiteEdgeIter{DST,<:BipartiteGraph{T}}, state=(1, 1, DST)) where T
     @unpack g = it
     nvars = ndsts(g)
     nvars == 0 && return nothing