378 lines
13 KiB
Python
378 lines
13 KiB
Python
|
# boykovkolmogorov.py - Boykov Kolmogorov algorithm for maximum flow problems.
|
||
|
|
#
|
||
|
|
# Copyright 2016-2019 NetworkX developers.
|
||
|
|
#
|
||
|
|
# This file is part of NetworkX.
|
||
|
|
#
|
||
|
|
# NetworkX is distributed under a BSD license; see LICENSE.txt for more
|
||
|
|
# information.
|
||
|
|
#
|
||
|
|
# Author: Jordi Torrents <jordi.t21@gmail.com>
|
||
|
|
"""
|
||
|
|
Boykov-Kolmogorov algorithm for maximum flow problems.
|
||
|
|
"""
|
||
|
|
from collections import deque
|
||
|
|
from operator import itemgetter
|
||
|
|
|
||
|
|
import networkx as nx
|
||
|
|
from networkx.algorithms.flow.utils import build_residual_network
|
||
|
|
|
||
|
|
__all__ = ['boykov_kolmogorov']
|
||
|
|
|
||
|
|
|
||
|
|
def boykov_kolmogorov(G, s, t, capacity='capacity', residual=None,
|
||
|
|
value_only=False, cutoff=None):
|
||
|
|
r"""Find a maximum single-commodity flow using Boykov-Kolmogorov algorithm.
|
||
|
|
|
||
|
|
This function returns the residual network resulting after computing
|
||
|
|
the maximum flow. See below for details about the conventions
|
||
|
|
NetworkX uses for defining residual networks.
|
||
|
|
|
||
|
|
This algorithm has worse case complexity $O(n^2 m |C|)$ for $n$ nodes, $m$
|
||
|
|
edges, and $|C|$ the cost of the minimum cut [1]_. This implementation
|
||
|
|
uses the marking heuristic defined in [2]_ which improves its running
|
||
|
|
time in many practical problems.
|
||
|
|
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
G : NetworkX graph
|
||
|
|
Edges of the graph are expected to have an attribute called
|
||
|
|
'capacity'. If this attribute is not present, the edge is
|
||
|
|
considered to have infinite capacity.
|
||
|
|
|
||
|
|
s : node
|
||
|
|
Source node for the flow.
|
||
|
|
|
||
|
|
t : node
|
||
|
|
Sink node for the flow.
|
||
|
|
|
||
|
|
capacity : string
|
||
|
|
Edges of the graph G are expected to have an attribute capacity
|
||
|
|
that indicates how much flow the edge can support. If this
|
||
|
|
attribute is not present, the edge is considered to have
|
||
|
|
infinite capacity. Default value: 'capacity'.
|
||
|
|
|
||
|
|
residual : NetworkX graph
|
||
|
|
Residual network on which the algorithm is to be executed. If None, a
|
||
|
|
new residual network is created. Default value: None.
|
||
|
|
|
||
|
|
value_only : bool
|
||
|
|
If True compute only the value of the maximum flow. This parameter
|
||
|
|
will be ignored by this algorithm because it is not applicable.
|
||
|
|
|
||
|
|
cutoff : integer, float
|
||
|
|
If specified, the algorithm will terminate when the flow value reaches
|
||
|
|
or exceeds the cutoff. In this case, it may be unable to immediately
|
||
|
|
determine a minimum cut. Default value: None.
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
R : NetworkX DiGraph
|
||
|
|
Residual network after computing the maximum flow.
|
||
|
|
|
||
|
|
Raises
|
||
|
|
------
|
||
|
|
NetworkXError
|
||
|
|
The algorithm does not support MultiGraph and MultiDiGraph. If
|
||
|
|
the input graph is an instance of one of these two classes, a
|
||
|
|
NetworkXError is raised.
|
||
|
|
|
||
|
|
NetworkXUnbounded
|
||
|
|
If the graph has a path of infinite capacity, the value of a
|
||
|
|
feasible flow on the graph is unbounded above and the function
|
||
|
|
raises a NetworkXUnbounded.
|
||
|
|
|
||
|
|
See also
|
||
|
|
--------
|
||
|
|
:meth:`maximum_flow`
|
||
|
|
:meth:`minimum_cut`
|
||
|
|
:meth:`preflow_push`
|
||
|
|
:meth:`shortest_augmenting_path`
|
||
|
|
|
||
|
|
Notes
|
||
|
|
-----
|
||
|
|
The residual network :samp:`R` from an input graph :samp:`G` has the
|
||
|
|
same nodes as :samp:`G`. :samp:`R` is a DiGraph that contains a pair
|
||
|
|
of edges :samp:`(u, v)` and :samp:`(v, u)` iff :samp:`(u, v)` is not a
|
||
|
|
self-loop, and at least one of :samp:`(u, v)` and :samp:`(v, u)` exists
|
||
|
|
in :samp:`G`.
|
||
|
|
|
||
|
|
For each edge :samp:`(u, v)` in :samp:`R`, :samp:`R[u][v]['capacity']`
|
||
|
|
is equal to the capacity of :samp:`(u, v)` in :samp:`G` if it exists
|
||
|
|
in :samp:`G` or zero otherwise. If the capacity is infinite,
|
||
|
|
:samp:`R[u][v]['capacity']` will have a high arbitrary finite value
|
||
|
|
that does not affect the solution of the problem. This value is stored in
|
||
|
|
:samp:`R.graph['inf']`. For each edge :samp:`(u, v)` in :samp:`R`,
|
||
|
|
:samp:`R[u][v]['flow']` represents the flow function of :samp:`(u, v)` and
|
||
|
|
satisfies :samp:`R[u][v]['flow'] == -R[v][u]['flow']`.
|
||
|
|
|
||
|
|
The flow value, defined as the total flow into :samp:`t`, the sink, is
|
||
|
|
stored in :samp:`R.graph['flow_value']`. If :samp:`cutoff` is not
|
||
|
|
specified, reachability to :samp:`t` using only edges :samp:`(u, v)` such
|
||
|
|
that :samp:`R[u][v]['flow'] < R[u][v]['capacity']` induces a minimum
|
||
|
|
:samp:`s`-:samp:`t` cut.
|
||
|
|
|
||
|
|
Examples
|
||
|
|
--------
|
||
|
|
>>> import networkx as nx
|
||
|
|
>>> from networkx.algorithms.flow import boykov_kolmogorov
|
||
|
|
|
||
|
|
The functions that implement flow algorithms and output a residual
|
||
|
|
network, such as this one, are not imported to the base NetworkX
|
||
|
|
namespace, so you have to explicitly import them from the flow package.
|
||
|
|
|
||
|
|
>>> G = nx.DiGraph()
|
||
|
|
>>> G.add_edge('x','a', capacity=3.0)
|
||
|
|
>>> G.add_edge('x','b', capacity=1.0)
|
||
|
|
>>> G.add_edge('a','c', capacity=3.0)
|
||
|
|
>>> G.add_edge('b','c', capacity=5.0)
|
||
|
|
>>> G.add_edge('b','d', capacity=4.0)
|
||
|
|
>>> G.add_edge('d','e', capacity=2.0)
|
||
|
|
>>> G.add_edge('c','y', capacity=2.0)
|
||
|
|
>>> G.add_edge('e','y', capacity=3.0)
|
||
|
|
>>> R = boykov_kolmogorov(G, 'x', 'y')
|
||
|
|
>>> flow_value = nx.maximum_flow_value(G, 'x', 'y')
|
||
|
|
>>> flow_value
|
||
|
|
3.0
|
||
|
|
>>> flow_value == R.graph['flow_value']
|
||
|
|
True
|
||
|
|
|
||
|
|
A nice feature of the Boykov-Kolmogorov algorithm is that a partition
|
||
|
|
of the nodes that defines a minimum cut can be easily computed based
|
||
|
|
on the search trees used during the algorithm. These trees are stored
|
||
|
|
in the graph attribute `trees` of the residual network.
|
||
|
|
|
||
|
|
>>> source_tree, target_tree = R.graph['trees']
|
||
|
|
>>> partition = (set(source_tree), set(G) - set(source_tree))
|
||
|
|
|
||
|
|
Or equivalently:
|
||
|
|
|
||
|
|
>>> partition = (set(G) - set(target_tree), set(target_tree))
|
||
|
|
|
||
|
|
References
|
||
|
|
----------
|
||
|
|
.. [1] Boykov, Y., & Kolmogorov, V. (2004). An experimental comparison
|
||
|
|
of min-cut/max-flow algorithms for energy minimization in vision.
|
||
|
|
Pattern Analysis and Machine Intelligence, IEEE Transactions on,
|
||
|
|
26(9), 1124-1137.
|
||
|
|
http://www.csd.uwo.ca/~yuri/Papers/pami04.pdf
|
||
|
|
|
||
|
|
.. [2] Vladimir Kolmogorov. Graph-based Algorithms for Multi-camera
|
||
|
|
Reconstruction Problem. PhD thesis, Cornell University, CS Department,
|
||
|
|
2003. pp. 109-114.
|
||
|
|
https://pub.ist.ac.at/~vnk/papers/thesis.pdf
|
||
|
|
|
||
|
|
"""
|
||
|
|
R = boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff)
|
||
|
|
R.graph['algorithm'] = 'boykov_kolmogorov'
|
||
|
|
return R
|
||
|
|
|
||
|
|
|
||
|
|
def boykov_kolmogorov_impl(G, s, t, capacity, residual, cutoff):
|
||
|
|
if s not in G:
|
||
|
|
raise nx.NetworkXError('node %s not in graph' % str(s))
|
||
|
|
if t not in G:
|
||
|
|
raise nx.NetworkXError('node %s not in graph' % str(t))
|
||
|
|
if s == t:
|
||
|
|
raise nx.NetworkXError('source and sink are the same node')
|
||
|
|
|
||
|
|
if residual is None:
|
||
|
|
R = build_residual_network(G, capacity)
|
||
|
|
else:
|
||
|
|
R = residual
|
||
|
|
|
||
|
|
# Initialize/reset the residual network.
|
||
|
|
# This is way too slow
|
||
|
|
#nx.set_edge_attributes(R, 0, 'flow')
|
||
|
|
for u in R:
|
||
|
|
for e in R[u].values():
|
||
|
|
e['flow'] = 0
|
||
|
|
|
||
|
|
# Use an arbitrary high value as infinite. It is computed
|
||
|
|
# when building the residual network.
|
||
|
|
INF = R.graph['inf']
|
||
|
|
|
||
|
|
if cutoff is None:
|
||
|
|
cutoff = INF
|
||
|
|
|
||
|
|
R_succ = R.succ
|
||
|
|
R_pred = R.pred
|
||
|
|
|
||
|
|
def grow():
|
||
|
|
"""Bidirectional breadth-first search for the growth stage.
|
||
|
|
|
||
|
|
Returns a connecting edge, that is and edge that connects
|
||
|
|
a node from the source search tree with a node from the
|
||
|
|
target search tree.
|
||
|
|
The first node in the connecting edge is always from the
|
||
|
|
source tree and the last node from the target tree.
|
||
|
|
"""
|
||
|
|
while active:
|
||
|
|
u = active[0]
|
||
|
|
if u in source_tree:
|
||
|
|
this_tree = source_tree
|
||
|
|
other_tree = target_tree
|
||
|
|
neighbors = R_succ
|
||
|
|
else:
|
||
|
|
this_tree = target_tree
|
||
|
|
other_tree = source_tree
|
||
|
|
neighbors = R_pred
|
||
|
|
for v, attr in neighbors[u].items():
|
||
|
|
if attr['capacity'] - attr['flow'] > 0:
|
||
|
|
if v not in this_tree:
|
||
|
|
if v in other_tree:
|
||
|
|
return (u, v) if this_tree is source_tree else (v, u)
|
||
|
|
this_tree[v] = u
|
||
|
|
dist[v] = dist[u] + 1
|
||
|
|
timestamp[v] = timestamp[u]
|
||
|
|
active.append(v)
|
||
|
|
elif v in this_tree and _is_closer(u, v):
|
||
|
|
this_tree[v] = u
|
||
|
|
dist[v] = dist[u] + 1
|
||
|
|
timestamp[v] = timestamp[u]
|
||
|
|
_ = active.popleft()
|
||
|
|
return None, None
|
||
|
|
|
||
|
|
def augment(u, v):
|
||
|
|
"""Augmentation stage.
|
||
|
|
|
||
|
|
Reconstruct path and determine its residual capacity.
|
||
|
|
We start from a connecting edge, which links a node
|
||
|
|
from the source tree to a node from the target tree.
|
||
|
|
The connecting edge is the output of the grow function
|
||
|
|
and the input of this function.
|
||
|
|
"""
|
||
|
|
attr = R_succ[u][v]
|
||
|
|
flow = min(INF, attr['capacity'] - attr['flow'])
|
||
|
|
path = [u]
|
||
|
|
# Trace a path from u to s in source_tree.
|
||
|
|
w = u
|
||
|
|
while w != s:
|
||
|
|
n = w
|
||
|
|
w = source_tree[n]
|
||
|
|
attr = R_pred[n][w]
|
||
|
|
flow = min(flow, attr['capacity'] - attr['flow'])
|
||
|
|
path.append(w)
|
||
|
|
path.reverse()
|
||
|
|
# Trace a path from v to t in target_tree.
|
||
|
|
path.append(v)
|
||
|
|
w = v
|
||
|
|
while w != t:
|
||
|
|
n = w
|
||
|
|
w = target_tree[n]
|
||
|
|
attr = R_succ[n][w]
|
||
|
|
flow = min(flow, attr['capacity'] - attr['flow'])
|
||
|
|
path.append(w)
|
||
|
|
# Augment flow along the path and check for saturated edges.
|
||
|
|
it = iter(path)
|
||
|
|
u = next(it)
|
||
|
|
these_orphans = []
|
||
|
|
for v in it:
|
||
|
|
R_succ[u][v]['flow'] += flow
|
||
|
|
R_succ[v][u]['flow'] -= flow
|
||
|
|
if R_succ[u][v]['flow'] == R_succ[u][v]['capacity']:
|
||
|
|
if v in source_tree:
|
||
|
|
source_tree[v] = None
|
||
|
|
these_orphans.append(v)
|
||
|
|
if u in target_tree:
|
||
|
|
target_tree[u] = None
|
||
|
|
these_orphans.append(u)
|
||
|
|
u = v
|
||
|
|
orphans.extend(sorted(these_orphans, key=dist.get))
|
||
|
|
return flow
|
||
|
|
|
||
|
|
def adopt():
|
||
|
|
"""Adoption stage.
|
||
|
|
|
||
|
|
Reconstruct search trees by adopting or discarding orphans.
|
||
|
|
During augmentation stage some edges got saturated and thus
|
||
|
|
the source and target search trees broke down to forests, with
|
||
|
|
orphans as roots of some of its trees. We have to reconstruct
|
||
|
|
the search trees rooted to source and target before we can grow
|
||
|
|
them again.
|
||
|
|
"""
|
||
|
|
while orphans:
|
||
|
|
u = orphans.popleft()
|
||
|
|
if u in source_tree:
|
||
|
|
tree = source_tree
|
||
|
|
neighbors = R_pred
|
||
|
|
else:
|
||
|
|
tree = target_tree
|
||
|
|
neighbors = R_succ
|
||
|
|
nbrs = ((n, attr, dist[n]) for n, attr in neighbors[u].items()
|
||
|
|
if n in tree)
|
||
|
|
for v, attr, d in sorted(nbrs, key=itemgetter(2)):
|
||
|
|
if attr['capacity'] - attr['flow'] > 0:
|
||
|
|
if _has_valid_root(v, tree):
|
||
|
|
tree[u] = v
|
||
|
|
dist[u] = dist[v] + 1
|
||
|
|
timestamp[u] = time
|
||
|
|
break
|
||
|
|
else:
|
||
|
|
nbrs = ((n, attr, dist[n]) for n, attr in neighbors[u].items()
|
||
|
|
if n in tree)
|
||
|
|
for v, attr, d in sorted(nbrs, key=itemgetter(2)):
|
||
|
|
if attr['capacity'] - attr['flow'] > 0:
|
||
|
|
if v not in active:
|
||
|
|
active.append(v)
|
||
|
|
if tree[v] == u:
|
||
|
|
tree[v] = None
|
||
|
|
orphans.appendleft(v)
|
||
|
|
if u in active:
|
||
|
|
active.remove(u)
|
||
|
|
del tree[u]
|
||
|
|
|
||
|
|
def _has_valid_root(n, tree):
|
||
|
|
path = []
|
||
|
|
v = n
|
||
|
|
while v is not None:
|
||
|
|
path.append(v)
|
||
|
|
if v == s or v == t:
|
||
|
|
base_dist = 0
|
||
|
|
break
|
||
|
|
elif timestamp[v] == time:
|
||
|
|
base_dist = dist[v]
|
||
|
|
break
|
||
|
|
v = tree[v]
|
||
|
|
else:
|
||
|
|
return False
|
||
|
|
length = len(path)
|
||
|
|
for i, u in enumerate(path, 1):
|
||
|
|
dist[u] = base_dist + length - i
|
||
|
|
timestamp[u] = time
|
||
|
|
return True
|
||
|
|
|
||
|
|
def _is_closer(u, v):
|
||
|
|
return timestamp[v] <= timestamp[u] and dist[v] > dist[u] + 1
|
||
|
|
|
||
|
|
source_tree = {s: None}
|
||
|
|
target_tree = {t: None}
|
||
|
|
active = deque([s, t])
|
||
|
|
orphans = deque()
|
||
|
|
flow_value = 0
|
||
|
|
# data structures for the marking heuristic
|
||
|
|
time = 1
|
||
|
|
timestamp = {s: time, t: time}
|
||
|
|
dist = {s: 0, t: 0}
|
||
|
|
while flow_value < cutoff:
|
||
|
|
# Growth stage
|
||
|
|
u, v = grow()
|
||
|
|
if u is None:
|
||
|
|
break
|
||
|
|
time += 1
|
||
|
|
# Augmentation stage
|
||
|
|
flow_value += augment(u, v)
|
||
|
|
# Adoption stage
|
||
|
|
adopt()
|
||
|
|
|
||
|
|
if flow_value * 2 > INF:
|
||
|
|
raise nx.NetworkXUnbounded('Infinite capacity path, flow unbounded above.')
|
||
|
|
|
||
|
|
# Add source and target tree in a graph attribute.
|
||
|
|
# A partition that defines a minimum cut can be directly
|
||
|
|
# computed from the search trees as explained in the docstrings.
|
||
|
|
R.graph['trees'] = (source_tree, target_tree)
|
||
|
|
# Add the standard flow_value graph attribute.
|
||
|
|
R.graph['flow_value'] = flow_value
|
||
|
|
return R
|