Edge weights #5: maximum flow

doc/weights.xml

@@ -219,3 +219,55 @@ fail]]></Example>
   </Description>
 </ManSection>
 <#/GAPDoc>
+
+<#GAPDoc Label="DigraphMaximumFlow">
+<ManSection>
+  <Attr Name="DigraphMaximumFlow" Arg="digraph, start, destination"/>
+  <Returns>A list of lists of integers.</Returns>
+  <Description>
+    If <A>digraph</A> is an edge-weighted digraph with vertices <A>start</A> and
+    <A>destination</A>, this returns a record representing the maximum flow from
+    <A>start</A> to <A>destination</A> in the digraph. <P/>
+
+    A <E>flow</E> is a function from the weighted edges of <A>digraph</A> to the
+    positive real numbers, such that:
+    <List>
+      <Item>
+        Each edge's flow is no more than its weight;
+      </Item>
+      <Item>
+        For each vertex other than <A>start</A> and <A>destination</A>, the sum
+        of flows for all incoming edges is equal to the sum of flows for all
+        outgoing edges;
+      </Item>
+      <Item>
+        The sum of flows of edges leaving <A>start</A> is equal to the sum of
+        flows of edges entering <A>destination</A> (this sum is denoted
+        <M>M</M>).
+      </Item>
+    </List>
+    A <E>maximum flow</E> is a flow that maximises the value of <M>M</M>. <P/>
+
+    The flow is represented as a list of lists where each entry is a number
+    representing the flow on the edge in the corresponding position in
+    <C>OutNeighbours(<A>digraph</A>)</C>.
+    Note that the value <M>M</M> of the flow can be found with
+    <C>Sum(DigraphMaximumFlow(<A>digraph</A>, <A>start</A>,
+    <A>destination</A>)[<A>start</A>])</C>. <P/>
+
+    This attribute is computed by an implementation of the push–relabel maximum
+    flow algorithm, which has time complexity <M>O(v^2 e)</M> where <M>v</M> is
+    the number of vertices of the digraph, and <M>e</M> is the number of
+    edges. <P/>
+
+    See <Ref Attr="EdgeWeights" Func="EdgeWeightedDigraph"/>.
+    <Example><![CDATA[
+gap> g := EdgeWeightedDigraph([[2, 2], [3], []], [[3, 2], [1], []]);
+<immutable multidigraph with 3 vertices, 3 edges>
+gap> flow := DigraphMaximumFlow(g, 1, 3);
+[ [ 1, 0 ], [ 1 ], [  ] ]
+gap> Sum(flow[1]);
+1]]></Example>
+  </Description>
+</ManSection>
+<#/GAPDoc>

doc/z-chap5.xml

@@ -31,6 +31,7 @@
     <#Include Label="EdgeWeightedDigraphMinimumSpanningTree">
     <#Include Label="EdgeWeightedDigraphShortestPaths">
     <#Include Label="EdgeWeightedDigraphShortestPath">
+    <#Include Label="DigraphMaximumFlow">
   </Section>
 
   <Section><Heading>Orders</Heading>

gap/weights.gd

@@ -34,3 +34,7 @@ DeclareGlobalFunction("DIGRAPHS_Edge_Weighted_Johnson");
 DeclareGlobalFunction("DIGRAPHS_Edge_Weighted_FloydWarshall");
 DeclareGlobalFunction("DIGRAPHS_Edge_Weighted_Bellman_Ford");
 DeclareGlobalFunction("DIGRAPHS_Edge_Weighted_Dijkstra");
+
+# 5. Maximum Flow
+DeclareOperation("DigraphMaximumFlow",
+                 [IsDigraph and HasEdgeWeights, IsPosInt, IsPosInt]);

gap/weights.gi

@@ -645,3 +645,138 @@ function(digraph, source)
 
   return rec(distances := distances, parents := parents, edges := edges);
 end);
+
+#############################################################################
+# 5. Maximum Flow
+#############################################################################
+
+InstallMethod(DigraphMaximumFlow, "for an edge weighted digraph",
+[IsDigraph and HasEdgeWeights, IsPosInt, IsPosInt],
+function(digraph, start, destination)
+  local push, relabel, discharge, weights, vertices, nrVertices, outs, capacity,
+        flowMatrix, seen, height, excess, queue, u, outNeighbours, i, v, w,
+        flows, delta;
+
+  # Push flow from u to v
+  push := function(u, v)
+    local delta;
+
+    delta := Minimum(excess[u], capacity[u][v] - flowMatrix[u][v]);
+
+    flowMatrix[u][v] := flowMatrix[u][v] + delta;
+    flowMatrix[v][u] := flowMatrix[v][u] - delta;
+    excess[u]        := excess[u] - delta;
+    excess[v]        := excess[v] + delta;
+
+    if excess[v] = delta then
+      PlistDequePushBack(queue, v);
+    fi;
+  end;
+
+  # Decide height of u, which must be above any vertex it can flow to
+  relabel := function(u)
+    local d, v;
+    d := infinity;
+    for v in vertices do
+      if capacity[u][v] - flowMatrix[u][v] > 0 then
+        d := Minimum(d, height[v]);
+      fi;
+    od;
+    if d < infinity then
+      height[u] := d + 1;
+    fi;
+  end;
+
+  # Push all excess flow from u to other vertices
+  discharge := function(u)
+    local v;
+    while excess[u] > 0 do
+      if seen[u] <= nrVertices then
+        v := seen[u];
+        if capacity[u][v] - flowMatrix[u][v] > 0
+            and height[u] > height[v] then
+          push(u, v);
+        else
+          seen[u] := seen[u] + 1;
+        fi;
+      else
+        relabel(u);
+        seen[u] := 1;
+      fi;
+    od;
+  end;
+
+  # Extract important data
+  weights    := EdgeWeights(digraph);
+  vertices   := DigraphVertices(digraph);
+  nrVertices := Length(vertices);
+  outs       := OutNeighbors(digraph);
+
+  # Check input
+  if not start in vertices then
+    ErrorNoReturn("<start> must be a vertex of <digraph>,");
+  elif not destination in vertices then
+    ErrorNoReturn("<destination> must be a vertex of <digraph>,");
+  fi;
+
+  # Setup data structures
+  capacity   := EmptyPlist(nrVertices);
+  flowMatrix := EmptyPlist(nrVertices);
+  seen       := EmptyPlist(nrVertices);
+  height     := EmptyPlist(nrVertices);
+  excess     := EmptyPlist(nrVertices);
+  queue      := PlistDeque();
+  for u in vertices do
+    capacity[u]   := ListWithIdenticalEntries(nrVertices, 0);
+    flowMatrix[u] := ListWithIdenticalEntries(nrVertices, 0);
+    seen[u]       := 1;  # highest vertex to which we've pushed flow from u
+    height[u]     := 0;  # proximity to start (further is lower)
+    excess[u]     := 0;  # flow coming into u that isn't yet leaving u
+    if u <> start and u <> destination then
+      PlistDequePushBack(queue, u);
+    fi;
+  od;
+
+  # Compute total capacity from each u to v
+  for u in vertices do
+    outNeighbours := outs[u];
+    for i in [1 .. Size(outNeighbours)] do
+      v := outNeighbours[i];  # the out neighbour
+      w := weights[u][i];     # the weight to the out neighbour
+      capacity[u][v] := capacity[u][v] + w;
+    od;
+  od;
+
+  # start vertex is "top", with unlimited flow coming out
+  height[start] := nrVertices;
+  excess[start] := infinity;
+
+  # Push flow from start to the other vertices
+  for v in vertices do
+    if v <> start then
+      push(start, v);
+    fi;
+  od;
+
+  # Discharge from vertices until there are none left to do
+  while not IsEmpty(queue) do
+    u := PlistDequePopFront(queue);
+    if u <> start and u <> destination then
+      discharge(u);
+    fi;
+  od;
+
+  # Convert to output format: a list of lists with same shape as weights
+  flows := List(weights, l -> EmptyPlist(Length(l)));
+  for u in vertices do
+    for i in [1 .. Length(outs[u])] do
+      v := outs[u][i];
+      delta := Minimum(flowMatrix[u][v], weights[u][i]);
+      delta := Maximum(0, delta);
+      flowMatrix[u][v] := flowMatrix[u][v] - delta;
+      flows[u][i] := delta;
+    od;
+  od;
+
+  return flows;
+end);

tst/standard/weights.tst

@@ -1,4 +1,3 @@
-
 #############################################################################
 ##
 #W  standard/weights.tst
@@ -15,10 +14,14 @@ gap> LoadPackage("digraphs", false);;
 
 #
 gap> DIGRAPHS_StartTest();
-gap> d := EdgeWeightedDigraph([[2], []], [[5], []]);
-<immutable digraph with 2 vertices, 1 edge>
+
+#############################################################################
+# 1. Edge Weights
+#############################################################################
 
 # create edge weighted digraph
+gap> d := EdgeWeightedDigraph([[2], []], [[5], []]);
+<immutable digraph with 2 vertices, 1 edge>
 gap> d := EdgeWeightedDigraph(Digraph([[2], []]), [[5], []]);
 <immutable digraph with 2 vertices, 1 edge>
 gap> EdgeWeightedDigraphTotalWeight(d);
@@ -59,6 +62,10 @@ Error, the sizes of the out neighbours and weights for vertex 1 must be equal,
 gap> d := EdgeWeightedDigraph([[2], []], [[5, 10], []]);
 Error, the sizes of the out neighbours and weights for vertex 1 must be equal,
 
+#############################################################################
+# 2. Copies of edge weights
+#############################################################################
+
 # changing edge weights mutable copy
 gap> d := EdgeWeightedDigraph([[2], [1]], [[5], [10]]);
 <immutable digraph with 2 vertices, 2 edges>
@@ -87,6 +94,10 @@ gap> d := EdgeWeightedDigraph([[2], [1]], [[-5], [10]]);
 gap> IsNegativeEdgeWeightedDigraph(d);
 true
 
+#############################################################################
+# 3. Minimum Spanning Trees
+#############################################################################
+
 # not connnected digraph
 gap> d := EdgeWeightedDigraph([[1], [2]], [[5], [10]]);
 <immutable digraph with 2 vertices, 2 edges>
@@ -133,6 +144,10 @@ gap> d := EdgeWeightedDigraph([[2, 2, 2], [1]], [[10, 5, 15], [7]]);
 gap> EdgeWeightedDigraphMinimumSpanningTree(d);
 <immutable digraph with 2 vertices, 1 edge>
 
+#############################################################################
+# 4. Shortest Path
+#############################################################################
+
 # Shortest paths: one node
 gap> d := EdgeWeightedDigraph([[]], [[]]);
 <immutable empty digraph with 1 vertex>
@@ -275,6 +290,84 @@ Error, the 2nd argument <source> must be a vertex of the 1st argument <digraph\
 gap> EdgeWeightedDigraphShortestPath(d, 1, 3);
 [ [ 1, 2, 3 ], [ 1, 1 ] ]
 
+#############################################################################
+# 5. Maximum Flow
+#############################################################################
+
+# Maximum flow: empty digraphs
+gap> d := EdgeWeightedDigraph([], []);
+<immutable empty digraph with 0 vertices>
+gap> DigraphMaximumFlow(d, 1, 1);
+Error, <start> must be a vertex of <digraph>,
+
+# Maximum flow: single vertex (also empty digraphs)
+gap> d := EdgeWeightedDigraph([[]], [[]]);
+<immutable empty digraph with 1 vertex>
+gap> DigraphMaximumFlow(d, 1, 1);
+[ [  ] ]
+
+# Maximum flow: source = dest
+gap> d := EdgeWeightedDigraph([[2], []], [[5], []]);
+<immutable digraph with 2 vertices, 1 edge>
+gap> DigraphMaximumFlow(d, 1, 1);
+[ [ 0 ], [  ] ]
+
+# Maximum flow: has loop
+gap> d := EdgeWeightedDigraph([[1, 2], []], [[5, 10], []]);
+<immutable digraph with 2 vertices, 2 edges>
+gap> DigraphMaximumFlow(d, 1, 2);
+[ [ 0, 10 ], [  ] ]
+
+# Maximum flow: invalid source
+gap> d := EdgeWeightedDigraph([[1, 2], []], [[5, 10], []]);
+<immutable digraph with 2 vertices, 2 edges>
+gap> DigraphMaximumFlow(d, 5, 2);
+Error, <start> must be a vertex of <digraph>,
+
+# Maximum flow: invalid sink
+gap> d := EdgeWeightedDigraph([[1, 2], []], [[5, 10], []]);
+<immutable digraph with 2 vertices, 2 edges>
+gap> DigraphMaximumFlow(d, 1, 5);
+Error, <destination> must be a vertex of <digraph>,
+
+# Maximum flow: sink not reachable
+gap> d := EdgeWeightedDigraph([[1], []], [[5], []]);
+<immutable digraph with 2 vertices, 1 edge>
+gap> DigraphMaximumFlow(d, 1, 2);
+[ [ 0 ], [  ] ]
+
+# Maximum flow: source has in neighbours
+gap> d := EdgeWeightedDigraph([[2], [3], []], [[5], [10], []]);
+<immutable digraph with 3 vertices, 2 edges>
+gap> DigraphMaximumFlow(d, 2, 3);
+[ [ 0 ], [ 10 ], [  ] ]
+
+# Maximum flow: sink has out neighbours
+gap> d := EdgeWeightedDigraph([[2], [3], [2]], [[5], [10], [7]]);
+<immutable digraph with 3 vertices, 3 edges>
+gap> DigraphMaximumFlow(d, 2, 3);                           
+[ [ 0 ], [ 10 ], [ 0 ] ]
+
+# Maximum flow: cycle
+gap> d := EdgeWeightedDigraph([[2], [3], [1]], [[5], [10], [7]]);
+<immutable digraph with 3 vertices, 3 edges>
+gap> DigraphMaximumFlow(d, 1, 3);
+[ [ 5 ], [ 5 ], [ 0 ] ]
+
+# Maximum flow: example from Wikipedia
+gap> gr := EdgeWeightedDigraph([[3, 4], [], [2, 4], [2]],
+>                              [[10, 5], [], [5, 15], [10]]);;
+gap> DigraphMaximumFlow(gr, 1, 2);
+[ [ 10, 5 ], [  ], [ 5, 5 ], [ 10 ] ]
+gap> DigraphMaximumFlow(gr, 3, 2);
+[ [ 0, 0 ], [  ], [ 5, 10 ], [ 10 ] ]
+
+# Maximum flow: example from Wikipedia article on Push-label maximum flow
+gap> gr := EdgeWeightedDigraph([[2], [3, 6], [4], [1, 6], [1, 3], []],
+>                              [[12], [3, 7], [10], [5, 10], [15, 4], []]);;
+gap> DigraphMaximumFlow(gr, 5, 6);
+[ [ 10 ], [ 3, 7 ], [ 7 ], [ 0, 7 ], [ 10, 4 ], [  ] ]
+
 #  DIGRAPHS_UnbindVariables
 gap> Unbind(d);
 gap> Unbind(tree);