Edmonds-Karp Max Flow/Min Cut

Compatibility with v0.3

doc/source/algorithms.rst

@@ -9,7 +9,7 @@ Graph Algorithms
 - topological sorting
 - shortest paths: Dijkstra, Floyd-Warshall, A*
 - minimum spanning trees: Prim, Kruskal
-- flow: Minimum Cut
+- flow: Minimum s-t Cut, Maximum Flow, Simple Minimum Cut
 - random graph generation
 - more algorithms are being implemented
 
@@ -342,10 +342,33 @@ Kruskal's algorithm finds a minimum spanning tree (or forest) by gradually uniti
 Flow
 -----------------------
 
-This package implements Simple Minimum Cut
+This package implements Minimum s-t Cut, Maximum Flow, and Simple Minimum Cut
+
+
+Minimum s-t Cut
+~~~~~~~~~~~~~~~
+
+The minimum cut that separates vertex s and vertex t. 
+
+.. py:function:: min_st_cut(graph, capacity)
+
+    :param graph:       the input graph
+    :param capacity:    the edge capacities (vector of floats)
+
+    :returns:   ``(parity, bestcut)``, where ``parity`` is a vector of boolean values that determines the partition and ``bestcut`` is the weight of the cut that makes this partition.
+
+Maximum Flow
+~~~~~~~~~~~~
+
+.. py:function:: max_flow(graph, capacity)
+
+    :param graph:       the input graph
+    :param capacity:    the edge capacities (vector of floats)
+
+    :returns:   ``maxflow``, where ``maxflow`` is the maximum flow of the network. 
 
 Simple Minimum Cut
-~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~
 
 Stoer's simple minimum cut gets the minimum cut of an undirected graph.
 

src/Graphs.jl

@@ -71,6 +71,9 @@ export
     # maximum_adjacency_visit
     MaximumAdjacency, AbstractMASVisitor, min_cut, maximum_adjacency_visit,
 
+    # edmonds_karp
+    min_st_cut, max_flow,
+
     # connected_components
     connected_components, strongly_connected_components,
 
@@ -136,6 +139,7 @@ include("breadth_first_visit.jl")
 include("depth_first_visit.jl")
 include("maximum_adjacency_visit.jl")
 
+include("edmonds_karp.jl")
 include("connected_components.jl")
 include("dijkstra_spath.jl")
 include("bellmanford.jl")

src/edmonds_karp.jl

@@ -0,0 +1,122 @@
+function min_st_cut{V,E}(g::AbstractGraph{V,E},s::V,t::V,capacity::Vector{Float64})
+  @graph_requires g incidence_list vertex_list
+  @assert is_directed(g)
+
+  r = residual_graph(g,s,capacity)
+  flow = edmonds_karp_max_flow!(r,s,t)
+  parity = dfs(r,s)
+  return parity, flow
+end
+
+function max_flow{V,E}(g::AbstractGraph{V,E},s::V,t::V,capacity::Vector{Float64})
+  @graph_requires g incidence_list vertex_list
+  @assert is_directed(g)
+
+  r = residual_graph(g,s,capacity)
+  return edmonds_karp_max_flow!(r,s,t)
+end
+
+function residual_graph{V,E}(g::AbstractGraph{V,E},s::V,capacity::Vector{Float64})
+  visited = fill(false,num_vertices(g))
+  res_g = inclist(vertices(g),ExEdge,is_directed=true)
+  residual_graph_sub!(g,res_g,s,visited,capacity)
+  return res_g
+end
+
+function residual_graph_sub!{V,E1,E2}(g::AbstractGraph{V,E1},r::AbstractGraph{V,E2},
+                                      s::V, visited::Vector{Bool}, capacity::Vector{Float64})
+  visited[vertex_index(s,g)] = true
+  for edge in out_edges(s,g)
+    i = edge_index(edge); u = edge.source; v = edge.target
+    d1 = Dict{UTF8String,Any}(); d1["capacity"] = capacity[i]; d1["flow"] = 0
+    d2 = Dict{UTF8String,Any}(); d2["capacity"] = capacity[i]; d2["flow"] = capacity[i]
+    edge     = ExEdge(i, u, v, d1)
+    rev_edge = ExEdge(i, v, u, d2)
+    add_edge!(r,edge)
+    add_edge!(r,rev_edge)
+    if !visited[vertex_index(v,g)]
+      residual_graph_sub!(g,r,v,visited,capacity)
+    end
+  end
+end
+
+function edmonds_karp_max_flow!{V,E}(g::AbstractGraph{V,E},s::V,t::V)
+  flow = 0
+
+  while true
+
+    #run BFS to find shortest s-t path
+    #store edges taken to get to each vertex in 'pred'
+    pred = bfs(g,s,t)
+
+    #stop if we weren't able to find a path from s to t
+    if !haskey(pred,t)
+      break
+    end
+
+    #Otherwise see how much flow we can send
+    df = Inf
+    edge = pred[t]
+    while true
+      df = min(df, edge.attributes["capacity"] - edge.attributes["flow"])
+      if haskey(pred,edge.source) 
+        edge = pred[edge.source]
+      else
+        break
+      end
+    end
+
+    #and update edges by that amount
+    edge = pred[t]
+    while true
+      #find rev edge
+      t_edges = out_edges(edge.target,g)
+      idx = find(x-> x.target==edge.source,t_edges) #there should be only one!
+      rev_edge = t_edges[idx[1]]
+
+      edge.attributes["flow"] += df
+      rev_edge.attributes["flow"] -= df
+
+      if haskey(pred,edge.source) 
+        edge = pred[edge.source]
+      else
+        break
+      end
+    end
+    flow += df
+  end
+  return flow
+end
+
+function bfs{V,E}(g::AbstractGraph{V,E},s::V,t::V)
+  q = DataStructures.Queue(V)
+  enqueue!(q,s)
+
+  pred = Dict{V,E}()
+
+  while length(q) > 0
+    cur = dequeue!(q)
+    for edge in out_edges(cur,g)
+      if !haskey(pred,edge.target) && edge.target != s && edge.attributes["capacity"] > edge.attributes["flow"]
+        pred[edge.target] = edge
+        enqueue!(q,edge.target)
+      end
+    end
+  end
+  return pred
+end
+
+function dfs{V,E}(g::AbstractGraph{V,E},s::V)
+  colormap = fill(false,num_vertices(g))
+  return dfs!(g,s,colormap)
+end
+
+function dfs!{V,E}(g::AbstractGraph{V,E},s::V, colormap::Vector{Bool})
+  colormap[vertex_index(s,g)] = true
+  for edge in out_edges(s,g)
+    if edge.attributes["capacity"] > edge.attributes["flow"] && !colormap[vertex_index(edge.target,g)]
+      dfs!(g,edge.target,colormap)
+    end
+  end
+  return colormap
+end

test/edmonds_karp.jl

@@ -0,0 +1,40 @@
+using Graphs
+using Base.Test
+
+#example from wikipedia
+
+g = inclist(collect(1:7),is_directed=true)
+
+#(u, v, c) edge from u to v with capacity c
+inputs = [
+(1, 2, 3.),
+(1, 4, 3.),
+(2, 3, 4.),
+(3, 1, 3.),
+(3, 4, 1.),
+(3, 5, 2.),
+(4, 5, 2.),
+(4, 6, 6.),
+(5, 2, 1.),
+(5, 7, 1.),
+(6, 7, 9.)]
+
+
+m = length(inputs)
+c = zeros(m)
+for i = 1 : m
+  add_edge!(g, inputs[i][1],inputs[i][2])
+  c[i] = inputs[i][3]
+end
+
+@assert num_vertices(g) == 7
+@assert num_edges(g) == 11
+
+parity, f = min_st_cut(g,1,7,c)
+
+@test length(parity) == 7
+@test parity == Bool[true,true,true,false,true,false,false]
+@test f == 5.0
+
+f = max_flow(g,1,7,c)
+@test f == 5.0

test/runtests.jl

@@ -19,7 +19,8 @@ tests = [
     "cliques",
     "random",
     "generators",
-    "maximum_adjacency_visit" ]
+    "maximum_adjacency_visit",
+    "edmonds_karp"]
 
 
 for t in tests