diff --git a/Project.toml b/Project.toml index cfffc35..d697870 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,10 @@ uuid = "664ebbd3-01ea-5ff0-9621-55aaaf3546e3" authors = ["Mathieu Besançon "] version = "0.1.0" +[deps] +LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" +SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d" + [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index 70b10fd..a4132ea 100644 --- a/README.md +++ b/README.md @@ -3,3 +3,6 @@ [![Build Status](https://travis-ci.org/JuliaGraphs/BlossomMatching.jl.svg?branch=master)](https://travis-ci.org/JuliaGraphs/BlossomMatching.jl) [![codecov.io](http://codecov.io/github/JuliaGraphs/BlossomMatching.jl/coverage.svg?branch=master)](http://codecov.io/github/JuliaGraphs/BlossomMatching.jl?branch=master) [![Coverage Status](https://coveralls.io/repos/github/JuliaGraphs/BlossomMatching.jl/badge.svg?branch=master)](https://coveralls.io/github/JuliaGraphs/BlossomMatching.jl?branch=master) + +This repository contains (under developement) the code for Edmonds Blossom Algorithm and its different variants. This repository is mainly focussed on +the implementation of Blossom V algorithm for computing min cost perfect matching in a general graph. \ No newline at end of file diff --git a/REQUIRE b/REQUIRE deleted file mode 100644 index 05b5ab4..0000000 --- a/REQUIRE +++ /dev/null @@ -1 +0,0 @@ -julia 1.0 diff --git a/src/BlossomMatching.jl b/src/BlossomMatching.jl index 75687df..9eb4b33 100644 --- a/src/BlossomMatching.jl +++ b/src/BlossomMatching.jl @@ -1,5 +1,11 @@ module BlossomMatching -greet() = print("Hello World!") +using LightGraphs +using SimpleTraits + +export max_cardinality_matching + +include("EdmondsMaxCardinality.jl") end # module + diff --git a/src/EdmondsMaxCardinality.jl b/src/EdmondsMaxCardinality.jl new file mode 100644 index 0000000..5c7c579 --- /dev/null +++ b/src/EdmondsMaxCardinality.jl @@ -0,0 +1,275 @@ +mutable struct AuxStruct{T} + mate :: Vector{T} #mate[i] stores the vertex which is matched to vertex i + parent :: Vector{T} #parent[i] stores the parent of a vertex in the tree + base :: Vector{T} #base[i] stores the base of the flower to which vertex i belongs. It equals i if vertex i doesn't belong to any flower. + queue :: Vector{T} #array used in traversing the tree + used :: Vector{Bool} #boolean array to mark used vertex + blossom :: Vector{Bool} #boolean array to mark vertices in a blossom +end + + +""" +Function to calculate lowest common ancestor of 2 vertices in a tree. +The algorithm has O(N) time complexity. +LCA algorithms with better time complexity e.g. O(log(N)) can be used but not used here as this is not of primary importance +to the algorithm. PRs welcome for this. +""" +function lowest_common_ancestor!(aux, a, b, nvg) + fill!(aux.used,false) + + ##Rise from vertex a to root, marking all even vertices + while true + a = aux.base[a] + aux.used[a] = one(nvg) + if aux.mate[a] == zero(nvg) + break + end + a = aux.parent[aux.mate[a]] + end + + ##Rise from the vertex b until we find the labeled vertex + while true + b = aux.base[b] + if aux.used[b] == one(nvg) + return b + end + b = aux.parent[aux.mate[b]] + end +end + + +""" +This function on the way from the top to the aux.base of the flower, marks true for the vertices +in the array aux.blossom[] and stores the ancestors for the even nodes in the tree. +The parameter child is the son for the vertex v itself(with the help of this parameter we close +the loop in the ancestors). +""" +function mark_path!(aux, v, b, child, nvg) + while aux.base[v] != b + aux.blossom[aux.base[v]] = one(nvg) + aux.blossom[aux.base[aux.mate[v]]] = one(nvg) + aux.parent[v] = child + child = aux.mate[v] + v = aux.parent[aux.mate[v]] + end + return nothing +end + + +""" +This function looks for augmenting path from the exposed vertex root +and returns the last vertex of this path, or zero(nvg) if the augmenting path is not found. +""" +function find_path!(aux, root, g) + nvg = nv(g) + + fill!(aux.used, 0) + + fill!(aux.parent, zero(nvg)) + + @inbounds for i in one(nvg):nvg + aux.base[i] = i + end + + aux.used[root] = one(nvg) + + qh = one(nvg) + qt = one(nvg) + + aux.queue[qt] = root + qt = qt + one(nvg) + + while qh10^3. But because of +the greedy initialization of the matching in the beginning, the algorithm works in about 2 second, even for +a CompleteGraph(10001). However this might not be the case for other such huge graphs because the greedy +initialization of the matching might not work for all type of graphs. Infact it is possible to create counter +test cases to fail the greedy intiliaztion. + +### Examples +```jldoctest + +julia> g = SimpleGraph([0 1 0 ; 1 0 1; 0 1 0]) +{3, 2} undirected simple Int64 graph + +julia> max_cardinality_matching(g) +(mate = [2, 1, 0], matched_edges = Array{Int64,1}[[1, 2]]) + + +julia> g=SimpleGraph(11) +{11, 0} undirected simple Int64 graph + +julia> add_edge!(g,1,2); + +julia> add_edge!(g,2,3); + +julia> add_edge!(g,3,4); + +julia> add_edge!(g,4,1); + +julia> add_edge!(g,3,5); + +julia> add_edge!(g,5,6); + +julia> add_edge!(g,6,7); + +julia> add_edge!(g,7,5); + +julia> add_edge!(g,5,8); + +julia> add_edge!(g,8,9); + +julia> add_edge!(g,10,11); + +julia> max_cardinality_matching(g) +(mate = [2, 1, 4, 3, 6, 5, 0, 9, 8, 11, 10], matched_edges = Array{Int64,1}[[1, 2], [3, 4], [5, 6], [8, 9], [10, 11]]) + +``` +""" +function max_cardinality_matching end +# see https://github.com/mauro3/SimpleTraits.jl/issues/47#issuecomment-327880153 for syntax +@traitfn function max_cardinality_matching(g::AG::(!IsDirected)) where {T<:Integer, AG <:AbstractGraph{T}} + + nvg = nv(g) + neg = ne(g) + + aux = AuxStruct{T}(T[],T[],T[],T[],Bool[],Bool[]) + + #Initializing vectors + aux.mate = fill(zero(nvg), nvg) + + aux.parent = fill(zero(nvg), nvg) + + sizehint!(aux.base, nvg) + @inbounds for i in one(nvg):nvg + push!(aux.base, i) + end + + aux.queue = fill(0, 2*nvg) + + sizehint!(aux.used, nvg) + aux.used = fill(0, nvg) + + sizehint!(aux.blossom, nvg) + aux.blossom = fill(0, nvg) + + #Using a simple greedy algorithm to mark the preliminary matching to begin with the algorithm. This speeds up + #the matching algorithm by several times. + #Better heuristic based algorithm can be used which start with near optimal matching thus reducing the number + #of augmentating paths and thus speeding up the algorithm. PRs are welcome for this. + @inbounds for v in one(nvg):nvg + if aux.mate[v]==zero(nvg) + for v_neighbor in neighbors(g, v) + if aux.mate[v_neighbor]==zero(nvg) + aux.mate[v_neighbor] = v + aux.mate[v] = v_neighbor + break + end + end + end + end + + #Iteratively going through all the vertices to find an unmarked vertex + @inbounds for u in one(nvg):nvg + if aux.mate[u]==zero(nvg) # If vertex u is not in matching. + v = find_path!(aux, u, g) # Find augmenting path starting with u as one of end points. + while v!=zero(nvg) # Alternate along the path from i to last_v (whole while loop is for that). + pv = aux.parent[v] # Increasing the number of matched edges by alternating the edges in + ppv = aux.mate[pv] # augmenting path. + aux.mate[v] = pv + aux.mate[pv] = v + v = ppv + end + end + end + + matched_edges = Vector{Vector{T}}() # Result vector containing end points of matched edges. + + @inbounds for u in one(nvg):nvg + if u