3.5. Ford-Fulkerson in Python#
We will now use Python to help us run the Ford-Fulkerson algorithm.
We start by defining the input in Python. We will use the package “NetworkX” to manipulate and display graphs. We will consider a new maximum flow instance. The nodes in this instance are numbered 1, 2, 3, 4, 5 and 6, and we want to find the maximum flow from node 1 to node 6.
Read the code below, and run it.
# To use a package, we use the command "import" in Python
import networkx as nx
# G is directed graph
G = nx.DiGraph()
# We define G based on the edges and the capacity on each edge.
# One way to do this is to specify a list of edges with each endpoint and the edge capacity.
# For instance, edge (1,2) has capacity 2, (1,3) has capacity 5, etc.
edgeList = [(1,2,2),(1,3,5),(2,4,1),(3,2,1),(3,5,3),(4,6,6),(5,2,1),(5,4,2),(5,6,1)]
# We create the graph based on this list.
G.add_weighted_edges_from(edgeList, 'cap')
# Each edge will have a flow associated with it. At the start of the algorithm, the flow on each edge is 0.
# We initialize this by iterating through the edges in a for-loop and setting the attribute 'flow' to 0.
for i, j in G.edges:
G.edges[i,j]['flow'] = 0
Note
Pay attention to the white space in the code. In Python, white space indentation is used to signal that certain commands are part of a group; for example, the last command above is part of the for-loop.
The code above initializes the graph, but it does not display anything. Let’s run some code that shows what we are doing — run the code below.
pos = ((0,0),(0,10),(10,20),(10,0),(30,20),(30,0),(40,10)) # node positions
nx.draw_networkx(G,pos,node_size=500,node_color='darkorange') # draw graph
As you can see, the flows and capacities are not displayed on the graph. To do this, you can run the following code.
# Define edge labels
label = {}
for i, j in G.edges:
label[(i,j)] = str( G.edges[i,j]["flow"] ) + " / " + str( G.edges[i,j]["cap"] )
nx.draw_networkx(G,pos,node_size=500,node_color='darkorange') # draw graph
nx.draw_networkx_edge_labels(G,pos,edge_labels=label); # draw labels
So this is the input that we want to run the Ford-Fulkerson algorithm on. As a reminder, every iteration of the algorithm works as follows:
Create the residual graph \(G_f\) based on the current flow
Find a path in \(G_f\) from \(s\) to \(t\) (here from 1 to 6)
If a path exists, update the flow
Otherwise, you have found a maximum flow
So first we want to create the residual graph \(G_f\). Recall that we create this graph by adding forward edges with residual capacity equal to the remaining capacity of the edge on which flow could be pushed and backward edges with residual capacity equal to the flow on that edge.
Q: Assume we have a feasible flow. Consider an edge with capacity \(c\) and flow \(f\). What inequality must hold true for \(c\) and \(f\)? Why?
Q: Now, assume \(f > 0\). What residual edges will be created? For each residual edge created, state whether the edge is a forward or backward edge, and give its residual capacity in terms of \(c\) and \(f\).
The following code implements the logic you just described.
residualGraph = nx.DiGraph() # blank graph
for i, j in G.edges:
c = G.edges[ i,j ]['cap'] # capacity on edge (i,j)
f = G.edges[ i,j ]['flow'] # flow on edge (i,j)
# If there is leftover capacity
if( c > f ):
residualGraph.add_edge(i,j)
residualGraph.edges[i,j]['residual capacity'] = c-f
residualGraph.edges[i,j]['forward edge'] = True
# If there is nonzero flow
if ( f > 0 ):
residualGraph.add_edge(j,i)
residualGraph.edges[j,i]['residual capacity'] = f
residualGraph.edges[j,i]['forward edge'] = False
Now, the code below plots the residual graph.
nx.draw_networkx(residualGraph,pos,node_size=500,node_color='limegreen',connectionstyle='arc3, rad=0.1')
residualcap = nx.get_edge_attributes(residualGraph, 'residual capacity')
nx.draw_networkx_edge_labels(residualGraph,pos,edge_labels=residualcap,label_pos=0.66);
The nodes are drawn in green as a visual reminder that we are looking at the residual graph. Also, the edges are drawn with a slight bend because we know there will often be both an arc \((i,j)\) and \((j,i)\) in the residual graph. The residual capacities are given closer to the tail of the edges.
Let’s now do an iteration of Ford-Fulkerson.
Q: Find a path from 1 to 6 in the residual graph, and define delta to be the minimum residual capacity of the edges in this path.
You can define a path in Python by listing the edges in the path as follows:
flowpath = [(1,2),(2,3)]
defines the variable “flowpath” as the path with edges (1,2) and (2,3).
# TODO
flowpath = []
delta =
Q: How will the flow change on the edges in the original graph with a corresponding forward edge in the s-t path in the residual graph? What about those with a corresponding backward edge in the s-t path in the residual graph?
The following code implements the logic you just described.
for i, j in flowpath:
if residualGraph.edges[i,j]['forward edge']:
G.edges[i,j]['flow'] = G.edges[i,j]['flow'] + delta
else:
G.edges[j,i]['flow'] = G.edges[j,i]['flow'] - delta
Now, we plot the new flow!
# Define edge labels
label = {}
for i, j in G.edges:
label[(i,j)] = str( G.edges[i,j]["flow"] ) + " / " + str( G.edges[i,j]["cap"] )
nx.draw_networkx(G,pos,node_size=500,node_color='darkorange') # draw graph
nx.draw_networkx_edge_labels(G,pos,edge_labels=label); # draw labels
This concludes the first iteration of the Ford-Fulkerson algorithm! Let’s do this again! Did you notice we used the same code to plot the flow graph? When we have code we want to reuse, we should make a “function”. You can do this in Python using the keyword def. As an example, let’s define a function to plot the graph with the current flow, called plot_flow.
def plot_flow():
# Define edge labels
label = {}
for i, j in G.edges:
label[(i,j)] = str( G.edges[i,j]["flow"] ) + " / " + str( G.edges[i,j]["cap"] )
nx.draw_networkx(G,pos,node_size=500,node_color='darkorange') # draw graph
nx.draw_networkx_edge_labels(G,pos,edge_labels=label); # draw labels
After running the code in the cell above, it now suffices to just type plot_flow() to see the current flow.
plot_flow()
Great! Let’s do the same thing for plotting the residual graph. (The global keyword ensures the variable residualGraph has a scope outside this function. Speaking of scope, this is outside the scope of this course).
def plot_residual_graph():
global residualGraph
residualGraph = nx.DiGraph() # blank graph
for i, j in G.edges:
c = G.edges[ i,j ]['cap'] # capacity on edge (i,j)
f = G.edges[ i,j ]['flow'] # flow on edge (i,j)
# If there is leftover capacity
if( c > f ):
residualGraph.add_edge(i,j)
residualGraph.edges[i,j]['residual capacity'] = c-f
residualGraph.edges[i,j]['forward edge'] = True
# If there is nonzero flow
if ( f > 0 ):
residualGraph.add_edge(j,i)
residualGraph.edges[j,i]['residual capacity'] = f
residualGraph.edges[j,i]['forward edge'] = False
nx.draw_networkx(residualGraph,pos,node_size=500,node_color='limegreen',connectionstyle='arc3, rad=0.1')
residualcap = nx.get_edge_attributes(residualGraph, 'residual capacity')
nx.draw_networkx_edge_labels(residualGraph,pos,edge_labels=residualcap,label_pos=0.66);
Again, we can now plot the residual graph using plot_residual_graph().
plot_residual_graph()
Note
The residual capacities are given closer to the tail of the edges than the head, again, because there will often be two edges between a pair of nodes.
Q: Find another 1 to 6 path in the residual graph, and find the minimum residual capacity of the edges on this path.
# TODO
flowpath = []
delta =
Let’s update the flow again. But first, we will reduce repeated code by introducing a function.
def update_flow():
for i, j in flowpath:
if residualGraph.edges[i,j]['forward edge']:
G.edges[i,j]['flow'] = G.edges[i,j]['flow'] + delta
else:
G.edges[j,i]['flow'] = G.edges[j,i]['flow'] - delta
Run this cell to update the flow.
update_flow()
And run this cell to view the new flow graph.
plot_flow()
Let’s do another iteration! Here is the updated residual graph:
plot_residual_graph()
Q: Find another path.
# TODO
flowpath = []
delta =
update_flow()
plot_flow()
plot_residual_graph()
Did we find an optimal flow? If you did the previous iterations correctly, there is no longer a path from 1 to 6 in the residual graph. Recall, this means that the flow is optimal, and you can show this by providing a cut that has capacity equal to the flow value. This cut can be found using a “labeling algorithm” that maintains a list of checked nodes that it still needs to process.
Labeling Algorithm#
Labeling algorithm
Mark \(s\) with a \(\checkmark\). Let the list initially contain only node \(s\).
Until the list is empty:
Take any node (e.g., the first one) off of the list. Suppose that \(i\) is the name of this node. For each arc leaving \(i\), that is, each \((i,j) \in E_f\), if \(j\) is unchecked, check it and add it to the list.
The nodes that have been checked are those nodes that are reachable from \(s\) by a path in \(G_f\). Note that the labeling algorithm is extremely similar to Dijkstra’s shortest path algorithm. However, it does not care about distance, only if a node is reachable.
We implement the labeling algorithm below.
# Uncheck all nodes except for 1
for i in G.nodes:
G.nodes[i]["check"] = False
G.nodes[1]["check"] = True
# Initialize the list of checked nodes that we still need to process
list = { 1 }
# Until the list is empty
while (len( list ) > 0):
# Take any node off of the list
i = list.pop()
# For each arc in the residual graph leaving i
for j in residualGraph.neighbors(i):
# If j is unchecked
if not G.nodes[j]["check"]:
G.nodes[j]["check"] = True # check it
list.add( j ) # and add it to the list
# Output the checked nodes
for i in G.nodes:
if G.nodes[i]["check"]:
print( i )
Q: What is the cut that the labeling algorithm finds? What is the capacity of the cut?
We can modify the code to find a path in the residual graph if there is one by keeping track of the node \(i\) we came from when we first check a node \(j\). We can do this by adding another attribute to the nodes, called “prev”, and defining this in the algorithm. As soon as we encounter node 6 we can backtrack using “prev” to find the path.
We implement this idea below. Look for the slight change in implementation from the previous one.
def labeling_algorithm():
# Uncheck all nodes except for 1
for i in G.nodes:
G.nodes[i]["check"] = False
G.nodes[1]["check"] = True
# Initialize the list of checked nodes that we still need to process
list = { 1 }
# Until the list is empty
while ( len( list ) > 0 ):
# Take any node off of the list
i = list.pop()
# For each arc in the residual graph leaving i
for j in residualGraph.neighbors(i):
# If j is unchecked
if not G.nodes[j]["check"]:
G.nodes[j]["check"] = True # check it
G.nodes[j]["prev"] = i # set prev
list.add( j ) # and add it to the list
Let’s restart with the same graph again.
G = nx.DiGraph() # blank directed graph
edgeList = [(1,2,2),(1,3,5),(2,4,1),(3,2,1),(3,5,3),(4,6,6),(5,2,1),(5,4,2),(5,6,1)] # create edges / capacities
G.add_weighted_edges_from( edgeList, 'cap' ) # add edges / capacities to graph
for i, j in G.edges:
G.edges[i,j]['flow'] = 0 # set flows to zero
Let’s plot the flow and the residual graphs.
plot_flow()
plot_residual_graph()
Before using our function labeling_algorithm to run the labeling algorithm on this instance, we will use a visual and interactive implementation.
import urllib.request
urllib.request.urlretrieve('https://engri-1101.github.io/textbook/modules/max_flow.py', "max_flow.py")
from max_flow import *
pos = ((0,0),(0,10),(10,20),(10,0),(30,20),(30,0),(40,10))
for i in range(1,7):
G.nodes[i]['pos'] = pos[i]
ex = max_flow(G)
ex.label(s=1,auto=False,show=True)
Now let’s use our implementation!
labeling_algorithm()
We can use the attribute “check” to see whether node 6 is reachable from node 1 in the residual graph.
print(G.nodes[6]["check"])
Q: Is node 6 reachable from node 1 in the residual graph? Why?
Let’s use the “prev” attribute to find a path from 1 to 6 in the residual graph. Run the cell to generate the augmenting path in the residual graph.
j = 6
flowpath = []
while j != 1:
i = G.nodes[j]["prev"]
flowpath.insert( 0, (i,j) )
j = i
print(flowpath)
And let’s find the minimum residual capacity on this path. First, we can get a list of the residual capacities for every edge in the path. Then, we use min to find the smallest which we call delta.
path_capacities = [residualGraph.edges[i,j]['residual capacity'] for i,j in flowpath]
print(path_capacities)
delta = min(path_capacities)
print(delta)
Now we can use the function defined earlier to update the flow!
update_flow()
plot_flow()
plot_residual_graph()
Let’s run the labeling algorithm again.
labeling_algorithm()
Let’s define a new function to find and print an \(s-t\) path in the residual graph and its corresponding delta.
def print_path_and_delta():
global flowpath, delta
j = 6
flowpath = []
while j != 1:
i = G.nodes[j]["prev"]
flowpath.insert( 0, (i,j) )
j = i
print(flowpath)
path_capacities = [residualGraph.edges[i,j]['residual capacity'] for i,j in flowpath]
delta = min(path_capacities)
print(delta)
Run this cell to make sure your function works as expected.
print_path_and_delta()
Let’s update the flow again.
update_flow()
plot_flow()
print_path_and_delta()
Putting Everything Together#
You’re probably getting a bit sick from running each cell every time. Let’s automate the procedure!
import matplotlib.pyplot as plt
def ford_fulkerson(G):
for i, j in G.edges:
G.edges[i,j]['flow'] = 0
plot_residual_graph() # plot the residual graph
labeling_algorithm() # run the labeling algorithm to find a s-t path in the residual graph
while G.nodes[6]["check"]: # while there is an s-t path in the residual graph
plt.figure()
print_path_and_delta() # show it
update_flow() # update the flow
plot_residual_graph() # plot the new residual graph
labeling_algorithm() # run the labeling algorithm to find a s-t path in the new residual graph
# Define the input graph
G = nx.DiGraph()
edgeList = [(1,2,2),(1,3,5),(2,4,1),(3,2,1),(3,5,3),(4,6,6),(5,2,1),(5,4,2),(5,6,1)]
G.add_weighted_edges_from(edgeList, 'cap')
# Run Ford-Fulkerson
ford_fulkerson(G)
Let’s have a look at the maximum flow.
plot_flow()
Use the visual/interactive labeling algorithm to identify the checked nodes.
pos = ((0,0),(0,10),(10,20),(10,0),(30,20),(30,0),(40,10))
for i in range(1,7):
G.nodes[i]['pos'] = pos[i]
ex = max_flow(G)
ex.ford_fulkerson(s=1, t=6, show=False)
ex.label(s=1, auto=False, show=True)
Q: Which nodes were checked in the last run of the labeling algorithm? Use this to give a minimum cut. What is its capacity? How does it relate to the value of the optimal flow?
We have now implemented the Ford-Fulkerson algorithm in Python! With a few more tweaks, you could use this on any graph.
Bonus: You have seen all the ingredients necessary to also implement Dijkstra’s algorithm for the shortest path problem. Implement this algorithm in Python.