initial transposition

Author: matbesancon

REQUIRE

@@ -1 +1,6 @@
 julia 0.6
+LightGraphs 0.9.0
+JuMP 0.13.2
+MatrixDepot
+BlossomV 0.1
+GLPKMathProgInterface 0.3.2

src/LightGraphsMatching.jl

@@ -1,5 +1,34 @@
+__precompile__(true)
 module LightGraphsMatching
+using LightGraphs
+export MatchingResult, maximum_weight_matching, maximum_weight_maximal_matching, minimum_weight_perfect_matching
 
-# package code goes here
+"""
+    type MatchingResult{T}
+        weight::T
+        mate::Vector{Int}
+    end
+
+A type representing the result of a matching algorithm.
+
+    weight: total weight of the matching
+
+    mate:    `mate[i] = j` if vertex `i` is matched to vertex `j`.
+             `mate[i] = -1` for unmatched vertices.
+"""
+struct MatchingResult{T<:Real}
+    weight::T
+    mate::Vector{Int}
+end
+
+import BlossomV
+include("blossomv.jl")
+
+using GLPKMathProgInterface: GLPKSolverMIP
+
+using JuMP
+include("lp.jl")
+include("maximum_weight_matching.jl")
 
 end # module
+

src/blossomv.jl

@@ -0,0 +1,70 @@
+"""
+    minimum_weight_perfect_matching{T<:Real}(g, w::Dict{Edge,T})
+    minimum_weight_perfect_matching{T<:Real}(g, w::Dict{Edge,T}, cutoff)
+
+Given a graph `g` and an edgemap `w` containing weights associated to edges,
+returns a matching with the mimimum total weight among the ones containing
+exactly `nv(g)/2` edges.
+
+Edges in `g` not present in `w` will not be considered for the matching.
+
+This function relies on the BlossomV.jl package, a julia wrapper
+around Kolmogorov's BlossomV algorithm.
+
+Eventually a `cutoff` argument can be given, to the reduce computational time
+excluding edges with weights higher than the cutoff.
+
+The returned object is of type `MatchingResult`.
+
+In case of error try to change the optional argument `tmaxscale` (default is `tmaxscale=10`).
+"""
+function minimum_weight_perfect_matching end
+
+function minimum_weight_perfect_matching(g::Graph, w::Dict{E,T}, cutoff, kws...) where {T<:Real, E<:Edge}
+    wnew = Dict{E, T}()
+    for (e, c) in w
+        if c <= cutoff
+            wnew[e] = c
+        end
+    end
+    return minimum_weight_perfect_matching(g, wnew; kws...)
+end
+
+function minimum_weight_perfect_matching(g::Graph, w::Dict{E,T}; tmaxscale=10.) where {T<:AbstractFloat, E<:Edge}
+    wnew = Dict{E, Int32}()
+    cmax = maximum(values(w))
+    cmin = minimum(values(w))
+    tmax = typemax(Int32)  / tmaxscale # /10 is kinda arbitrary,
+                                # hopefully high enough to not incurr in overflow problems
+    for (e, c) in w
+        wnew[e] = round(Int32, (c-cmin) / (cmax-cmin) * tmax)
+    end
+    match = minimum_weight_perfect_matching(g, wnew)
+    weight = T(0)
+    for i=1:nv(g)
+        j = match.mate[i]
+        if j > i
+            weight += w[E(i,j)]
+        end
+    end
+    return MatchingResult(weight, match.mate)
+end
+
+function minimum_weight_perfect_matching(g::Graph, w::Dict{E,T}) where {T<:Integer, E<:Edge}
+    m = BlossomV.Matching(nv(g))
+    for (e, c) in w
+        BlossomV.add_edge(m, src(e)-1, dst(e)-1, c)
+    end
+    BlossomV.solve(m)
+
+    mate = fill(-1, nv(g))
+    totweight = T(0)
+    for i=1:nv(g)
+        j = BlossomV.get_match(m, i-1) + 1
+        mate[i] = j <= 0 ? -1 : j
+        if i < j
+            totweight += w[Edge(i,j)]
+        end
+    end
+    return MatchingResult(totweight, mate)
+end

src/lp.jl

@@ -0,0 +1,104 @@
+"""
+    maximum_weight_maximal_matching{T<:Real}(g, w::Dict{Edge,T})
+    maximum_weight_maximal_matching{T<:Real}(g, w::Dict{Edge,T}, cutoff)
+
+Given a bipartite graph `g` and an edgemap `w` containing weights associated to edges,
+returns a matching with the maximum total weight among the ones containing the
+greatest number of edges.
+
+Edges in `g` not present in `w` will not be considered for the matching.
+
+The algorithm relies on a linear relaxation on of the matching problem, which is
+guaranteed to have integer solution on bipartite graps.
+
+Eventually a `cutoff` argument can be given, to reduce computational times
+excluding edges with weights lower than the cutoff.
+
+The package JuMP.jl and one of its supported solvers is required.
+
+The returned object is of type `MatchingResult`.
+"""
+function maximum_weight_maximal_matching end
+
+function maximum_weight_maximal_matching(g::Graph, w::Dict{Edge,T}, cutoff; solver = GLPKSolverMIP) where {T<:Real}
+    wnew = Dict{Edge,T}()
+    for (e,x) in w
+        if x >= cutoff
+            wnew[e] = x
+        end
+    end
+
+    return maximum_weight_maximal_matching(g, wnew)
+end
+
+
+function maximum_weight_maximal_matching(g::Graph, w::Dict{Edge,T}; solver = GLPKSolverMIP) where {T<:Real}
+# TODO support for graphs with zero degree nodes
+# TODO apply separately on each connected component
+    bpmap = bipartite_map(g)
+    length(bpmap) != nv(g) && error("Graph is not bipartite")
+    v1 = findin(bpmap, 1)
+    v2 = findin(bpmap, 2)
+    if length(v1) > length(v2)
+        v1, v2 = v2, v1
+    end
+
+    nedg = 0
+    edgemap = Dict{Edge,Int}()
+    for (e,_) in w
+        nedg += 1
+        edgemap[e] = nedg
+        edgemap[reverse(e)] = nedg
+    end
+
+    model = Model(solver=solver())
+    @variable(model, x[1:length(w)] >= 0)
+
+    for i in v1
+        idx = Vector{Int}()
+        for j in neighbors(g, i)
+            if haskey(edgemap, Edge(i,j))
+                push!(idx, edgemap[Edge(i,j)])
+            end
+        end
+        if length(idx) > 0
+            @constraint(model, sum{x[id], id=idx} == 1)
+        end
+    end
+
+    for j in v2
+        idx = Vector{Int}()
+        for i in neighbors(g, j)
+            if haskey(edgemap, Edge(i,j))
+                push!(idx, edgemap[Edge(i,j)])
+            end
+        end
+
+        if length(idx) > 0
+            @constraint(model, sum{x[id], id=idx} <= 1)
+        end
+    end
+
+    @objective(model, Max, sum{c * x[edgemap[e]], (e,c)=w})
+
+    status = solve(model)
+    status != :Optimal && error("JuMP solver failed to find optimal solution.")
+    sol = getvalue(x)
+
+    all(Bool[s == 1 || s == 0 for s in sol]) || error("Found non-integer solution.")
+
+    cost = getobjectivevalue(model)
+
+    mate = fill(-1, nv(g))
+    for e in edges(g)
+        if haskey(w, e)
+            inmatch = convert(Bool, sol[edgemap[e]])
+            if inmatch
+                mate[src(e)] = dst(e)
+                mate[dst(e)] = src(e)
+            end
+        end
+    end
+
+    return MatchingResult(cost, mate)
+end

src/maximum_weight_matching.jl

@@ -0,0 +1,73 @@
+"""
+maximum_weight_matching{T <:Real}(g::Graph, w::Dict{Edge,T} = Dict{Edge,Int64}())
+
+Given a graph `g` and an edgemap `w` containing weights associated to edges,
+returns a matching with the maximum total weight.
+`w` is a dictionary that maps edges i => j to weights.
+If no weight parameter is given, all edges will be considered to have weight 1
+(results in max cardinality matching)
+
+The efficiency of the algorithm depends on the input graph:
+  - If the graph is bipartite, then the LP relaxation is integral.
+  - If the graph is not bipartite, then it requires a MIP solver and
+  the computation time may grow exponentially.
+
+The package JuMP.jl and one of its supported solvers is required.
+
+Returns MatchingResult containing:
+  - a solve status (indicating whether the problem was solved to optimality)
+  - the optimal cost
+  - a list of each vertex's match (or -1 for unmatched vertices)
+"""
+function maximum_weight_matching end
+
+function maximum_weight_matching(g::Graph,
+          w::Dict{Edge,T} = Dict{Edge,Int64}(i => 1 for i in collect(edges(g)));
+          solver = GLPKSolverMIP) where {T <:Real}
+
+    model = Model(solver = solver())
+    n = nv(g)
+    edge_list = collect(edges(g))
+
+    # put the edge weights in w in the right order to be compatible with edge_list
+    for edge in keys(w)
+      redge = reverse(edge)
+      if !is_ordered(edge) && !haskey(w, redge) # replace i=>j by j=>i if necessary.
+        w[redge] = w[edge]
+      end
+    end
+
+    if setdiff(edge_list, keys(w)) != [] # If some edges do not have a key in w.
+      error("Some edge weights are missing, check that keys i => j in w satisfy i <= j")
+    end
+
+    if is_bipartite(g)
+      @variable(model, x[edge_list] >= 0) # no need to enforce integrality
+    else
+      @variable(model, x[edge_list] >= 0, Int) # requires MIP solver
+    end
+    @objective(model, Max, sum(x[edge]*w[edge] for edge in edge_list))
+    @constraint(model, c1[i=1:n],
+                sum(x[Edge(i,j)] for j=filter(l -> l > i, neighbors(g,i))) +
+                sum(x[Edge(j,i)] for j=filter(l -> l <= i, neighbors(g,i)))
+                <= 1)
+
+    status = solve(model)
+    solution = getvalue(x)
+    cost = getobjectivevalue(model)
+    ## TODO: add the option of returning the solve status as part of the MatchingResult type.
+    return MatchingResult(cost, dict_to_arr(n, solution))
+end
+
+""" Returns an array of mates from a dictionary that maps edges to {0,1} """
+function dict_to_arr(n::Int64, solution::JuMP.JuMPArray{T,1,Tuple{Array{E,1}}}) where {T<: Real, E<: Edge}
+  mate = fill(-1,n)
+  for i in keys(solution)
+    key = i[1] # i is a tuple with 1 element.
+    if solution[key] >= 1 - 1e-5 # Some tolerance to numerical approximations by the solver.
+        mate[src(key)] = dst(key)
+        mate[dst(key)] = src(key)
+    end
+  end
+  return mate
+end