Merge e9b2fe4 into 744b362

Project.toml

@@ -3,6 +3,10 @@ uuid = "664ebbd3-01ea-5ff0-9621-55aaaf3546e3"
 authors = ["Mathieu Besançon <mathieu.besancon@gmail.com>"]
 version = "0.1.0"
 
+[deps]
+LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
+SimpleTraits = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
+
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

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.

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
+

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 qh<qt
+        v = aux.queue[qh]
+        qh = qh + one(nvg) 
+        for v_neighbor in neighbors(g, v)
+            if (aux.base[v] == aux.base[v_neighbor]) || (aux.mate[v] == v_neighbor)
+                continue
+            end
+
+            # The edge closes the cycle of odd length, i.e. a flower is found. 
+            # An odd-length cycle is detected when the following condition is met.
+            # In this case, we need to compress the flower.
+            if (v_neighbor == root) || (aux.mate[v_neighbor] != zero(nvg)) && (aux.parent[aux.mate[v_neighbor]] != zero(nvg))
+                #######Code for compressing the flower######
+                curbase = lowest_common_ancestor!(aux, v, v_neighbor, nvg)
+                for i in one(nvg):nvg
+                    aux.blossom[i] = 0
+                end
+
+                mark_path!(aux, v, curbase, v_neighbor, nvg)
+                mark_path!(aux, v_neighbor, curbase, v, nvg)
+
+                for i in one(nvg):nvg
+                    if aux.blossom[aux.base[i]] == one(nvg) 
+                        aux.base[i] = curbase
+                        if  aux.used[i] == 0
+                            aux.used[i] = one(nvg)
+                            aux.queue[qt] = i
+                            qt = qt + one(nvg) 
+                        end
+                    end
+                end
+                ############################################
+
+            # Otherwise, it is an "usual" edge, we act as in a normal wide search. 
+            # The only subtlety - when checking that we have not visited this vertex yet, we must not look 
+            # into the array aux.used, but instead into the array aux.parent as it is filled for the odd visited vertices.
+            # If we did not go to the top yet, and it turned out to be unsaturated, then we found an 
+            # augmenting path ending at the top v_neighbor, return it.
+
+            elseif aux.parent[v_neighbor] == zero(nvg)
+                aux.parent[v_neighbor] = v
+                if aux.mate[v_neighbor] == zero(nvg) 
+                    return v_neighbor
+                end
+                v_neighbor = aux.mate[v_neighbor]
+                aux.used[v_neighbor] = one(nvg)   
+                aux.queue[qt] = v_neighbor
+                qt = qt + one(nvg)
+            end
+        end
+    end
+
+    return zero(nvg)
+end
+
+
+"""
+    max_cardinality_matching(g)
+Compute the maximum cardinality matching in an undirected (weighted/unweighted) graph 
+`g` using Edmonds Blossom algorithm.
+(https://en.wikipedia.org/wiki/Blossom_algorithm).
+ 
+Returns a namedtuple(named 'solution') containing 2 vectors. 
+solution.mate contains the matched vertex of a vertex(also called as mate of vertex). 
+solution.matched_edges is a list of all matched edges represented by 2 end vertices. 
+ 
+### Performance
+Time Complexity : O(n^3)
+There are a total of n iterations, each of which performs a wide pass for O(m), in addition,
+there can be flower squeezing operations - which are O(n).  
+Thus the general asymptotics of the algorithm will be O(n(m+n^2)) = O(n^3).
+
+Complexity : O(n) {Excluding the memory required for storing graph}
+n = Number of vertices
+m = Number of edges
+
+max_vertices = 1000. O(n^3) algorithm shouldn't work in short amount(<=10sec) of time if n>10^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<aux.mate[u]
+            temp = Vector{T}()
+            push!(temp,u)
+            push!(temp,aux.mate[u])
+            push!(matched_edges,temp)       
+        end
+    end
+
+    result = (mate = aux.mate, matched_edges = matched_edges)
+
+    return result
+end

test/EdmondsMaxCardinality.jl

@@ -0,0 +1,87 @@
+@testset "EdmondsMaxCardinality" begin
+
+    g = PathGraph(5)
+    match = @inferred(max_cardinality_matching(g))
+    @test match.mate == [2, 1, 4, 3, 0]
+    @test match.matched_edges == [[1, 2], [3, 4]]
+
+    g = CycleGraph(11)
+    match = @inferred(max_cardinality_matching(g))
+    @test match.mate == [2, 1, 4, 3, 6, 5, 8, 7, 10, 9, 0]
+    @test match.matched_edges == [[1, 2], [3, 4], [5, 6], [7, 8], [9, 10]]
+
+    g = CompleteGraph(11)
+    match = @inferred(max_cardinality_matching(g))
+    @test match.mate == [2, 1, 4, 3, 6, 5, 8, 7, 10, 9, 0]
+    @test match.matched_edges == [[1, 2], [3, 4], [5, 6], [7, 8], [9, 10]]
+
+    g = SimpleGraph(11)
+    add_edge!(g,1,2)
+    add_edge!(g,2,3)
+    add_edge!(g,3,4)
+    add_edge!(g,4,1)   
+    add_edge!(g,3,5)
+    add_edge!(g,5,6)
+    add_edge!(g,6,7)
+    add_edge!(g,7,5)
+    add_edge!(g,5,8)
+    add_edge!(g,8,9)
+    add_edge!(g,10,11)
+    match = @inferred(max_cardinality_matching(g))
+    @test match.mate == [2, 1, 4, 3, 6, 5, 0, 9, 8, 11, 10]
+    @test match.matched_edges == [[1, 2], [3, 4], [5, 6], [8, 9], [10, 11]] 
+
+    g = SimpleGraph(9)
+    add_edge!(g,1,2)
+    add_edge!(g,2,3)
+    add_edge!(g,3,1)
+    add_edge!(g,1,8)
+    add_edge!(g,7,8)
+    add_edge!(g,7,3)
+    add_edge!(g,3,6)
+    add_edge!(g,3,9)
+    add_edge!(g,3,4)
+    add_edge!(g,5,6)
+    add_edge!(g,5,4)
+    match = @inferred(max_cardinality_matching(g))
+    @test match.mate == [2, 1, 4, 3, 6, 5, 8, 7, 0]
+    @test match.matched_edges == [[1, 2], [3, 4], [5, 6], [7, 8]]
+
+    g = SimpleGraph(12)
+    add_edge!(g,1,2)
+    add_edge!(g,1,3)
+    add_edge!(g,3,2)
+    add_edge!(g,5,2)
+    add_edge!(g,6,2)
+    add_edge!(g,4,2)
+    add_edge!(g,6,7)
+    add_edge!(g,9,7)
+    add_edge!(g,8,7)
+    add_edge!(g,8,11)
+    add_edge!(g,9,10)
+    add_edge!(g,9,8)
+    add_edge!(g,10,11)
+    add_edge!(g,10,12)
+    match = @inferred(max_cardinality_matching(g))
+    @test match.mate == [3, 4, 1, 2, 0, 7, 6, 9, 8, 11, 10, 0]
+    @test match.matched_edges == [[1, 3], [2, 4], [6, 7], [8, 9], [10, 11]]
+
+    g = SimpleGraph(12)
+    add_edge!(g,1,2)
+    add_edge!(g,1,3)
+    add_edge!(g,3,2)
+    add_edge!(g,5,2)
+    add_edge!(g,6,2)
+    add_edge!(g,4,2)
+    add_edge!(g,9,7)
+    add_edge!(g,8,7)
+    add_edge!(g,8,11)
+    add_edge!(g,9,10)
+    add_edge!(g,9,8)
+    add_edge!(g,10,11)
+    add_edge!(g,10,12)
+    match = @inferred(max_cardinality_matching(g))
+    @test match.mate == [3, 4, 1, 2, 0, 0, 9, 11, 7, 12, 8, 10]
+    @test match.matched_edges == [[1, 3], [2, 4], [7, 9], [8, 11], [10, 12]] 
+
+end