Merge pull request #738 from SciML/myb/bg

Add more methods on BipartiteGraph and add SystemStructure

Author: YingboMa

Project.toml

@@ -24,6 +24,7 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -50,9 +51,10 @@ RecursiveArrayTools = "2.3"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.4, 0.5"
 SafeTestsets = "0.0.1"
+Setfield = "0.7"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicUtils = "0.7.3"
+SymbolicUtils = "0.7.4"
 TreeViews = "0.3"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"

src/ModelingToolkit.jl

@@ -168,6 +168,8 @@ Get the set of parameters variables for the given system.
 """
 function parameters end
 
+include("bipartite_graph.jl")
+
 include("variables.jl")
 include("context_dsl.jl")
 include("differentials.jl")
@@ -181,6 +183,7 @@ include("domains.jl")
 include("register_function.jl")
 
 include("systems/abstractsystem.jl")
+include("systems/systemstructure.jl")
 
 include("systems/diffeqs/odesystem.jl")
 include("systems/diffeqs/sdesystem.jl")
@@ -211,6 +214,7 @@ include("extra_functions.jl")
 
 export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
+export SystemStructure
 export JumpSystem
 export ODEProblem, SDEProblem
 export NonlinearProblem, NonlinearProblemExpr

src/bipartite_graph.jl

@@ -0,0 +1,231 @@
+using UnPack
+using SparseArrays
+using LightGraphs
+using Setfield
+
+###
+### Edges & Vertex
+###
+@enum VertType SRC DST ALL
+
+struct BipartiteEdge{I<:Integer} <: LightGraphs.AbstractEdge{I}
+    src::I
+    dst::I
+    function BipartiteEdge(src::I, dst::V) where {I,V}
+        T = promote_type(I, V)
+        new{T}(T(src), T(dst))
+    end
+end
+
+LightGraphs.src(edge::BipartiteEdge) = edge.src
+LightGraphs.dst(edge::BipartiteEdge) = edge.dst
+
+function Base.show(io::IO, edge::BipartiteEdge)
+    @unpack src, dst = edge
+    print(io, "[src: ", src, "] => [dst: ", dst, "]")
+end
+
+Base.:(==)(a::BipartiteEdge, b::BipartiteEdge) = src(a) == src(b) && dst(a) == dst(b)
+
+###
+### Graph
+###
+"""
+$(TYPEDEF)
+
+A bipartite graph representation between two, possibly distinct, sets of vertices
+(source and dependencies). Maps source vertices, labelled `1:N₁`, to vertices
+on which they depend (labelled `1:N₂`).
+
+# Fields
+$(FIELDS)
+
+# Example
+```julia
+using ModelingToolkit
+
+ne = 4
+srcverts = 1:4
+depverts = 1:2
+
+# six source vertices
+fadjlist = [[1],[1],[2],[2],[1],[1,2]]
+
+# two vertices they depend on
+badjlist = [[1,2,5,6],[3,4,6]]
+
+bg = BipartiteGraph(7, fadjlist, badjlist)
+```
+"""
+mutable struct BipartiteGraph{I<:Integer} <: LightGraphs.AbstractGraph{I}
+    ne::Int
+    fadjlist::Vector{Vector{I}} # `fadjlist[src] => dsts`
+    badjlist::Vector{Vector{I}} # `badjlist[dst] => srcs`
+end
+
+"""
+```julia
+Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}
+```
+
+Test whether two [`BipartiteGraph`](@ref)s are equal.
+"""
+function Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}
+    iseq = (bg1.ne == bg2.ne)
+    iseq &= (bg1.fadjlist == bg2.fadjlist)
+    iseq &= (bg1.badjlist == bg2.badjlist)
+    iseq
+end
+
+"""
+$(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))
+
+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.length(::BipartiteGraph) = error("length is not well defined! Use `ne` or `nv`.")
+
+if isdefined(LightGraphs, :has_contiguous_vertices)
+    LightGraphs.has_contiguous_vertices(::Type{<:BipartiteGraph}) = false
+end
+LightGraphs.is_directed(::Type{<:BipartiteGraph}) = false
+LightGraphs.vertices(g::BipartiteGraph) = (𝑠vertices(g), 𝑑vertices(g))
+𝑠vertices(g::BipartiteGraph) = axes(g.fadjlist, 1)
+𝑑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]
+LightGraphs.ne(g::BipartiteGraph) = g.ne
+LightGraphs.nv(g::BipartiteGraph) = sum(length, vertices(g))
+LightGraphs.edgetype(g::BipartiteGraph{I}) where I = BipartiteEdge{I}
+
+nsrcs(g::BipartiteGraph) = length(𝑠vertices(g))
+ndsts(g::BipartiteGraph) = length(𝑑vertices(g))
+
+function LightGraphs.has_edge(g::BipartiteGraph, edge::BipartiteEdge)
+    @unpack src, dst = edge
+    (src in 𝑠vertices(g) && dst in 𝑑vertices(g)) || return false  # edge out of bounds
+    insorted(𝑠neighbors(src), dst)
+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)
+    @unpack fadjlist, badjlist = g
+    verts = vertices(g)
+    s, d = src(edge), dst(edge)
+    (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.")
+    @inbounds list = fadjlist[s]
+    index = searchsortedfirst(list, d)
+    @inbounds (index <= length(list) && list[index] == d) && return false  # edge already in graph
+    insert!(list, index, d)
+
+    g.ne += 1
+    @inbounds list = badjlist[d]
+    index = searchsortedfirst(list, s)
+    insert!(list, index, s)
+    return true  # edge successfully added
+end
+
+function LightGraphs.add_vertex!(g::BipartiteGraph{T}, type::VertType) where T
+    if type === DST
+        push!(g.badjlist, T[])
+    elseif type === SRC
+        push!(g.fadjlist, T[])
+    else
+        error("type ($type) must be either `DST` or `SRC`")
+    end
+    return true  # vertex successfully added
+end
+
+###
+### Edges iteration
+###
+LightGraphs.edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(ALL))
+𝑠edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(SRC))
+𝑑edges(g::BipartiteGraph) = BipartiteEdgeIter(g, Val(DST))
+
+struct BipartiteEdgeIter{T,G} <: LightGraphs.AbstractEdgeIter
+    g::G
+    type::Val{T}
+end
+
+Base.length(it::BipartiteEdgeIter) = ne(it.g)
+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
+    @unpack g = it
+    neqs = nsrcs(g)
+    neqs == 0 && return nothing
+    eq, jvar = state
+
+    while eq <= neqs
+        eq′ = eq
+        vars = 𝑠neighbors(g, eq′)
+        if jvar > length(vars)
+            eq += 1
+            jvar = 1
+            continue
+        end
+        edge = BipartiteEdge(eq′, vars[jvar])
+        state = (eq, jvar + 1, SRC)
+        return edge, state
+    end
+    return nothing
+end
+
+function Base.iterate(it::BipartiteEdgeIter{DST,BipartiteGraph{T}}, state=(1, 1, DST)) where T
+    @unpack g = it
+    nvars = ndsts(g)
+    nvars == 0 && return nothing
+    ieq, jvar = state
+
+    while jvar <= nvars
+        eqs = 𝑑neighbors(g, jvar)
+        if ieq > length(eqs)
+            ieq = 1
+            jvar += 1
+            continue
+        end
+        edge = BipartiteEdge(eqs[ieq], jvar)
+        state = (ieq + 1, jvar, DST)
+        return edge, state
+    end
+    return nothing
+end
+
+function Base.iterate(it::BipartiteEdgeIter{ALL,<:BipartiteGraph}, state=nothing)
+    if state === nothing
+        ss = iterate((@set it.type = Val(SRC)))
+    elseif state[3] === SRC
+        ss = iterate((@set it.type = Val(SRC)), state)
+    elseif state[3] == DST
+        ss = iterate((@set it.type = Val(DST)), state)
+    end
+    if ss === nothing && state[3] == SRC
+        return iterate((@set it.type = Val(DST)))
+    else
+        return ss
+    end
+end
+
+###
+### Utils
+###
+function LightGraphs.incidence_matrix(g::BipartiteGraph, val=true)
+    I = Int[]
+    J = Int[]
+    for i in 𝑠vertices(g), n in 𝑠neighbors(g, i)
+        push!(I, i)
+        push!(J, n)
+    end
+    S = sparse(I, J, val, nsrcs(g), ndsts(g))
+end

src/context_dsl.jl

@@ -10,7 +10,7 @@ SymbolicUtils.@number_methods(Sym{Parameter{Real}},
                               term(f, a, b), skipbasics)
 
 SymbolicUtils.symtype(s::Symbolic{Parameter{T}}) where T = T
-SymbolicUtils.similarterm(t::Term{T}, f, args) where {T<:Parameter} = Term{T}(f, args)
+SymbolicUtils.similarterm(t::Term{<:Parameter}, f, args) = Term(f, args)
 
 Base.convert(::Type{Num}, x::Symbolic{Parameter{T}}) where {T<:Number} = Num(x)
 

src/differentials.jl

@@ -29,6 +29,7 @@ struct Differential <: Function
 end
 (D::Differential)(x) = Term{symtype(x)}(D, [x])
 (D::Differential)(x::Num) = Num(D(value(x)))
+SymbolicUtils.promote_symtype(::Differential, x) = x
 
 Base.show(io::IO, D::Differential) = print(io, "(D'~", D.x, ")")
 

src/systems/abstractsystem.jl

@@ -297,4 +297,3 @@ function (f::AbstractSysToExpr)(O)
     end
     return build_expr(:call, Any[operation(O); f.(arguments(O))])
 end
-