Network Flow

A flow network is a directed graph G = (V, E) with two distinguished vertices: a source s and a sink t. Every edge (u, v) ∈ E has a non-negative integer capacity c(u, v) ≥ 0. For edges not in E we define c(u, v) = 0.

A flow f is a function on ordered pairs (u, v) that satisfies two constraints:

1. Capacity constraint: 0 ≤ f(u, v) ≤ c(u, v) for all u, v ∈ V.
2. Flow conservation: For every vertex u ∈ V \ {s, t}, the total flow into u equals the total flow out of u: Σ_v f(v, u) = Σ_v f(u, v).

The value of the flow |f| is the total flow leaving the source: |f| = Σ_v f(s, v) − Σ_v f(v, s). The maximum flow problem asks: find a flow of maximum value.

Skew symmetry is a useful notational trick: define f(u, v) = −f(v, u) for all u, v. This lets you write flow conservation as Σ_v f(u, v) = 0 for all u ≠ s, t, which is cleaner for proofs. The capacity constraint becomes f(u, v) ≤ c(u, v).

Analogy: think of a city water network. Pipes (edges) have maximum throughput (capacity). Junctions (internal vertices) can neither create nor destroy water (conservation). The pumping station is s, the reservoir is t. How much water per second can you pump from station to reservoir?

# Representing a flow network as an adjacency structure
# Each edge stores: to, capacity, flow, and reverse_edge_index
from collections import defaultdict

class FlowNetwork:
    def __init__(self, n):
        self.n = n
        self.graph = defaultdict(list)  # adjacency list
    
    def add_edge(self, u, v, cap):
        # Forward edge
        self.graph[u].append([v, cap, 0, len(self.graph[v])])
        # Reverse edge (capacity 0, used for residual)
        self.graph[v].append([u, 0, 0, len(self.graph[u]) - 1])

# Example: a 4-node network
# s=0, t=3
G = FlowNetwork(4)
G.add_edge(0, 1, 10)  # s -> 1, cap 10
G.add_edge(0, 2, 10)  # s -> 2, cap 10
G.add_edge(1, 3, 10)  # 1 -> t, cap 10
G.add_edge(2, 3, 10)  # 2 -> t, cap 10
G.add_edge(1, 2, 1)   # 1 -> 2, cap 1 (the bottleneck bridge)

The residual graph G_f = (V, E_f) is the key technical device that makes flow algorithms work. Given a flow f on G, the residual graph captures how much additional flow can be pushed along each edge — in either direction.

For each edge (u, v) ∈ E, the residual capacity is:

• Forward residual: c_f(u, v) = c(u, v) − f(u, v) — how much more can go from u to v.
• Backward residual: c_f(v, u) = f(u, v) — you can "cancel" f(u, v) units of flow by sending flow backwards.

An augmenting path is any simple path from s to t in the residual graph G_f where every edge has positive residual capacity. The bottleneck of an augmenting path P is δ = min_{(u,v) ∈ P} c_f(u, v).

If an augmenting path P with bottleneck δ exists, you can increase the flow value by δ: for each (u,v) in P, set f(u,v) += δ and f(v,u) −= δ (the backward cancellation). This preserves both the capacity constraint and conservation.

Crucial insight: the backward edges in the residual graph represent rerouting, not literal backward flow. If you previously sent flow from A to B, a backward residual edge B→A means you can reroute that flow elsewhere — effectively undoing a prior decision. This is what makes greedy-then-correct algorithms possible.

Lemma (augmenting path theorem): A flow f is maximum if and only if the residual graph G_f contains no augmenting path from s to t. This is one of three equivalent characterizations given by the max-flow min-cut theorem.

def bfs_residual(graph, s, t, parent):
    """BFS in residual graph. Returns True if augmenting path found.
    parent[v] = u means we reached v from u."""
    visited = set([s])
    queue = [s]
    while queue:
        u = queue.pop(0)
        for edge in graph[u]:
            v, cap, flow, _ = edge
            residual = cap - flow
            if v not in visited and residual > 0:
                visited.add(v)
                parent[v] = (u, graph[u].index(edge))
                if v == t:
                    return True
                queue.append(v)
    return False

The Ford-Fulkerson method is the foundational max-flow algorithm, proposed by L.R. Ford Jr. and D.R. Fulkerson in 1956. It is a method (not a specific algorithm) because it doesn't specify how to find augmenting paths.

Ford-Fulkerson Method:
  Initialize f(u,v) = 0 for all edges
  while there exists an augmenting path P in G_f:
      δ = bottleneck(P)
      augment flow along P by δ
  return f

Correctness follows from the max-flow min-cut theorem (proved in the next concept). When no augmenting path exists, the flow is maximum.

Termination and complexity — here is where it gets subtle. If capacities are integers, each augmentation increases |f| by at least 1, so the algorithm terminates in at most |f*| iterations where f* is the max flow value. With BFS or DFS path-finding, each iteration costs O(E), giving O(E · |f|)* total.

This is a pseudo-polynomial runtime — polynomial in |f*|, which can be exponentially large in the input size. Worse, with irrational capacities, Ford-Fulkerson can actually fail to terminate (the flow values converge to a non-maximum value). The classic bad example is a network where each DFS augmentation picks a path that increases flow by only 1, even when the max flow is huge.

When it's fine: If capacities are small integers, Ford-Fulkerson with DFS is simple and fast. For general inputs, you need a smarter path selection strategy — enter Edmonds-Karp.

def ford_fulkerson_dfs(graph, s, t):
    """Ford-Fulkerson with DFS path finding.
    graph: FlowNetwork adjacency list (see above)
    Returns max flow value."""
    max_flow = 0
    
    def dfs(u, t, pushed, visited):
        if u == t:
            return pushed
        visited.add(u)
        for edge in graph[u]:
            v, cap, flow, rev = edge
            if v not in visited and cap - flow > 0:
                d = dfs(v, t, min(pushed, cap - flow), visited)
                if d > 0:
                    edge[2] += d          # increase forward flow
                    graph[v][rev][2] -= d  # decrease reverse flow
                    return d
        return 0
    
    while True:
        pushed = dfs(s, t, float('inf'), set())
        if pushed == 0:
            break
        max_flow += pushed
    
    return max_flow

This is the crown jewel of network flow theory. Let us state and prove it carefully.

Definition: An (s, t)-cut (S, T) is a partition of V into sets S (containing s) and T (containing t). The capacity of cut (S, T) is the sum of capacities of edges going from S to T: cap(S, T) = Σ_{u∈S, v∈T} c(u, v).

Note: edges from T to S do not count toward cut capacity. A cut is like severing all the pipes that cross from the s-side to the t-side.

Lemma (weak duality): For any flow f and any (s,t)-cut (S, T), |f| ≤ cap(S, T).

Proof: |f| = net flow across any cut separating s from t (this takes a small algebraic argument using flow conservation). Since flow on each edge is bounded by capacity, the net flow ≤ Σ_{u∈S,v∈T} c(u,v) = cap(S,T). ∎

Max-Flow Min-Cut Theorem: The following three conditions are equivalent:

1. f is a maximum flow.
2. There is no augmenting path from s to t in G_f.
3. There exists an (s,t)-cut (S, T) with cap(S, T) = |f|.

Proof sketch: (1 ⟹ 2) trivially (augmenting would increase the flow). (2 ⟹ 3) Let S = vertices reachable from s in G_f. Since no augmenting path exists, t ∉ S. For every edge (u,v) with u∈S, v∈T: c_f(u,v) = 0 (else v would be reachable), so f(u,v) = c(u,v). For edges (v,u) with v∈T, u∈S: c_f(v,u) = f(u,v) ≥ 0 so these contribute 0 to the cut. Hence |f| = cap(S,T). (3 ⟹ 1) By weak duality, no flow can exceed cap(S,T), and |f| = cap(S,T), so f is maximum. ∎

The constructive min-cut: when Ford-Fulkerson terminates, run BFS from s in G_f. The reachable set S and unreachable set T form the minimum cut — you get both the max flow and a witness proving its optimality.

Duality perspective: This is strong LP duality made combinatorial. The primal LP maximizes flow; the dual LP minimizes cut capacity. The theorem says the gap is zero — primal and dual optima coincide exactly, with integer solutions whenever capacities are integers.

def min_cut(graph, s, t):
    """After running max-flow, find the min cut.
    Returns (S, T) partition and cut edges."""
    # BFS/DFS in residual graph from s
    from collections import deque
    visited = set()
    queue = deque([s])
    visited.add(s)
    
    while queue:
        u = queue.popleft()
        for edge in graph[u]:
            v, cap, flow, _ = edge
            residual = cap - flow
            if v not in visited and residual > 0:
                visited.add(v)
                queue.append(v)
    
    S = visited  # reachable from s in residual graph
    T = set(graph.keys()) - S
    
    # Cut edges: from S to T at full capacity
    cut_edges = []
    for u in S:
        for edge in graph[u]:
            v, cap, flow, _ = edge
            if v in T and cap > 0:  # original edge (not reverse)
                cut_edges.append((u, v, cap))
    
    return S, T, cut_edges

Edmonds and Karp (1972) gave the first strongly polynomial implementation of Ford-Fulkerson: always augment along the shortest augmenting path (fewest edges), found by BFS. This innocent-sounding change yields a O(VE²) time bound — independent of capacity values.

Key lemma: Under BFS augmentation, the shortest-path distance d_f(s, v) in the residual graph is monotonically non-decreasing over all augmentations (it never decreases as we push flow).

Intuition: When you augment along a shortest path, you saturate at least one edge (the bottleneck). Saturating (u,v) adds a reverse edge (v,u) to the residual graph. This reverse edge can only create paths that are longer than the original shortest path. Therefore distances can only increase.

Theorem: Edmonds-Karp runs in O(VE²) time.

Proof sketch: An edge (u,v) can become the bottleneck (be saturated) only after d_f(s,u) increases by at least 2 since the last time it was saturated (because the reverse edge is used first, increasing u's distance). Since distances are bounded by V, each edge can be a bottleneck at most O(V) times, giving O(VE) total augmentations. Each BFS is O(E), so the total is O(VE²). ∎

In practice: For dense graphs Edmonds-Karp is excellent. For unit-capacity graphs (common in matching), you can do better. For highly capacitated graphs, Dinic's algorithm (which pushes flow along all shortest paths simultaneously using blocking flows) achieves O(V²E) and is faster in practice. But Edmonds-Karp is the go-to for correctness-first implementations and competitive programming.

from collections import deque

def edmonds_karp(graph, s, t, n):
    """
    Edmonds-Karp: BFS-augmenting Ford-Fulkerson.
    graph: adjacency list as list of [to, cap, flow, rev_idx]
    Returns maximum flow value.
    Time: O(V * E^2)
    """
    max_flow = 0
    
    while True:
        # BFS to find shortest augmenting path
        parent = [-1] * n       # parent[v] = index of edge in graph[parent_node]
        parent_node = [-1] * n
        parent_node[s] = s
        queue = deque([s])
        
        while queue and parent_node[t] == -1:
            u = queue.popleft()
            for i, (v, cap, flow, _) in enumerate(graph[u]):
                if parent_node[v] == -1 and cap - flow > 0:
                    parent_node[v] = u
                    parent[v] = i
                    queue.append(v)
        
        if parent_node[t] == -1:
            break  # No augmenting path found
        
        # Find bottleneck along the path
        path_flow = float('inf')
        v = t
        while v != s:
            u = parent_node[v]
            idx = parent[v]
            path_flow = min(path_flow, graph[u][idx][1] - graph[u][idx][2])
            v = u
        
        # Update flows along the path
        v = t
        while v != s:
            u = parent_node[v]
            idx = parent[v]
            graph[u][idx][2] += path_flow          # forward: increase flow
            rev = graph[u][idx][3]
            graph[v][rev][2] -= path_flow           # reverse: decrease flow
            v = u
        
        max_flow += path_flow
    
    return max_flow

One of the most elegant reductions in algorithms: maximum bipartite matching reduces to maximum flow.

Setup: Given a bipartite graph G = (L ∪ R, E) where every edge goes from L to R, find the largest set of edges with no vertex repeated (a matching).

Reduction: Build a flow network:

1. Add a source s with edges s → l for every l ∈ L (capacity 1).
2. Add a sink t with edges r → t for every r ∈ R (capacity 1).
3. For every edge (l, r) ∈ E, add edge l → r (capacity 1).

Now: a max flow of value k in this network corresponds to a maximum matching of size k. Each unit of flow uses one left vertex, one right vertex, and crosses one bipartite edge — exactly a matching edge.

Correctness: Integer capacities → Ford-Fulkerson finds an integer flow → the flow is 0/1 on each edge → the set of l→r edges carrying flow 1 is a valid matching. By max-flow min-cut, this matching is maximum.

Runtime: O(E · √V) using Hopcroft-Karp (the bipartite-specific algorithm), or O(VE) using Edmonds-Karp. For sparse graphs Hopcroft-Karp is preferred, but the flow reduction is simpler to implement and still fast enough for competitive programming (V, E ≤ 10^4).

König's theorem falls out immediately: in a bipartite graph, the size of the maximum matching equals the size of the minimum vertex cover. This is just max-flow min-cut applied to the bipartite flow network — the min cut translates directly to the min vertex cover.

Applications: Assigning workers to tasks, pairing students with projects, scheduling jobs on machines, image pixel labeling — any time you have two sides and want to maximize pairings subject to capacity-1 constraints.

def bipartite_max_matching(left_n, right_n, edges):
    """
    Reduce bipartite matching to max flow.
    Nodes: 0=source, 1..left_n = left, left_n+1..left_n+right_n = right,
           left_n+right_n+1 = sink
    edges: list of (l, r) pairs (0-indexed)
    """
    s = 0
    t = left_n + right_n + 1
    n = t + 1
    
    # Build adjacency list for flow network
    # Each entry: [to, cap, flow, rev_idx]
    graph = [[] for _ in range(n)]
    
    def add_edge(u, v, cap):
        graph[u].append([v, cap, 0, len(graph[v])])
        graph[v].append([u, 0,   0, len(graph[u]) - 1])
    
    # Source to left vertices
    for l in range(1, left_n + 1):
        add_edge(s, l, 1)
    
    # Right vertices to sink
    for r in range(1, right_n + 1):
        add_edge(left_n + r, t, 1)
    
    # Bipartite edges
    for (l, r) in edges:
        add_edge(l + 1, left_n + r + 1, 1)  # l is 0-indexed
    
    return edmonds_karp(graph, s, t, n)  # reuse from above

# Example: 3 left, 3 right, edges (0,0),(0,1),(1,1),(2,2)
# matching = 3 (perfect matching)
print(bipartite_max_matching(3, 3, [(0,0),(0,1),(1,1),(2,2)]))

One of the most powerful consequences of the max-flow min-cut theorem is the integrality theorem:

Theorem: If all capacities in a flow network are integers, then there exists a maximum flow in which every f(u,v) is an integer.

Proof: Ford-Fulkerson (with any augmenting path finder) starts with integer flow 0 and increases it by the integer bottleneck δ at each step. Induction: if all current flows are integers, then all residual capacities are integers, so δ is an integer, so the updated flows are integers. ∎

This is enormously useful: it means that combinatorial problems with 0/1 structure (matchings, packings, coverings) that reduce to flow will automatically get integer optimal solutions, without any rounding or branch-and-bound.

LP duality perspective: The maximum flow problem can be written as a linear program:

• Maximize Σ_v f(s,v) − Σ_v f(v,s)
• Subject to: 0 ≤ f(u,v) ≤ c(u,v) for all (u,v), flow conservation.

Its dual LP is the minimum cut problem. LP strong duality says primal optimum = dual optimum — and the max-flow min-cut theorem is exactly this statement, proved combinatorially. When the constraint matrix is totally unimodular (which flow networks' constraint matrices are), integer programs automatically have integer LP relaxation optima. This is why integrality comes for free in flow problems but not in general integer programming.

Hall's theorem (bipartite matching): A bipartite graph has a perfect matching iff for every subset S ⊆ L, |N(S)| ≥ |S|. This is proved by applying max-flow min-cut to the bipartite flow network and interpreting the min cut condition. Hall's theorem becomes a special case of the general flow duality machinery.

Edmonds-Karp is great, but for large networks you want Dinic's algorithm (Yefim Dinitz, 1970), which achieves O(V²E) — and for unit-capacity networks, O(E√V), the same as Hopcroft-Karp.

Level graph: In the residual graph G_f, define the level of each vertex as its BFS distance from s. The level graph L_f contains only edges (u,v) from G_f where level(v) = level(u) + 1. Only forward-BFS edges are kept.

Blocking flow: A flow f' in the level graph is a blocking flow if every s-t path in the level graph has at least one saturated edge (i.e., there is no augmenting path in the level graph after adding f').

Dinic's algorithm:

Initialize f = 0
Repeat:
  Build level graph L_f via BFS
  If t is unreachable: return f
  Find a blocking flow f' in L_f
  Augment f by f'

Why faster? After each phase, the shortest-path distance from s to t strictly increases (because the blocking flow saturates every short path). Since distances are bounded by V, there are at most V phases. A blocking flow can be found in O(VE) using DFS with "dead-end pruning" (advance-retreat algorithm). Total: O(V²E).

For unit-capacity graphs, distances increase by at least 1 per phase and after O(√E) phases the residual graph has no more augmenting paths (because the remaining flow is small). This gives O(E√V) = O(E√E) — matching Hopcroft-Karp.

In competitive programming, Dinic's is the standard max-flow implementation: it handles all cases well, runs fast in practice (much faster than the worst-case bound), and the code is clean once you understand blocking flows with the pointer optimization (current arc optimization).

from collections import deque

def dinic(graph, s, t, n):
    """Dinic's algorithm. O(V^2 * E).
    graph: list of [to, cap, rev_idx] per node."""
    
    def bfs_level():
        level[:] = [-1] * n
        level[s] = 0
        q = deque([s])
        while q:
            u = q.popleft()
            for v, cap, _ in graph[u]:
                if cap > 0 and level[v] == -1:
                    level[v] = level[u] + 1
                    q.append(v)
        return level[t] != -1
    
    def dfs_blocking(u, pushed):
        if u == t:
            return pushed
        while iter_[u] < len(graph[u]):
            e = graph[u][iter_[u]]
            v, cap, rev = e
            if cap > 0 and level[v] == level[u] + 1:
                d = dfs_blocking(v, min(pushed, cap))
                if d > 0:
                    e[1] -= d               # reduce forward cap
                    graph[v][rev][1] += d   # increase reverse cap
                    return d
            iter_[u] += 1  # dead-end pruning: advance pointer
        return 0
    
    level = [-1] * n
    total = 0
    while bfs_level():
        iter_ = [0] * n  # current arc pointer per node
        while True:
            f = dfs_blocking(s, float('inf'))
            if f == 0:
                break
            total += f
    return total

The true power of max-flow reveals itself through reductions — recognizing that a hard-looking problem is secretly asking for a max flow.

1. Vertex-disjoint paths: How many edge-disjoint (or vertex-disjoint) s-t paths exist? By Menger's theorem, this equals the minimum edge cut (or vertex cut). Reduction: edge-disjoint paths → set each edge capacity to 1 and find max flow. Vertex-disjoint → split each vertex v into v_in and v_out with a capacity-1 edge between them.

2. Baseball elimination: Team x is eliminated if even the best possible scenario can't help them tie for first. Build a flow network where a game between teams i and j is a "game node" with supply equal to remaining games, feeding into team nodes i and j (with constraints that team x cannot win). If max flow saturates all game edges, x is not eliminated. This problem became famous when Ken Steiglitz and others showed it reduces beautifully to flow.

3. Project selection / closed-set problem: You have projects with profits/costs and dependencies (doing project p requires doing prerequisite q). Find the most profitable subset. Reduction: build source → project node (profit), cost node → sink (cost), and dependency edges with infinite capacity. Min cut = max profit lost = optimal selection via s-t cut.

4. Image segmentation (Graph Cuts): In computer vision, each pixel is a node. Source = foreground label, sink = background label. Edge weights encode unary potentials (how likely each pixel is foreground/background) and pairwise smoothness terms. Min s-t cut gives a minimum-energy segmentation. This is the foundation of GrabCut and other interactive segmentation tools (Boykov & Kolmogorov, 2004).

5. Circulation with lower bounds: Sometimes edges have both lower and upper bounds on flow (L(u,v) ≤ f(u,v) ≤ U(u,v)). This generalizes max-flow. The standard reduction creates a new super-source and super-sink, converts lower bounds to demands, and solves two max-flow instances.

Recognition pattern: Look for problems where you're routing something (flow, paths, assignments) through a network with capacity constraints. If the structure is a DAG with capacity constraints and conservation-like equations, it's likely a flow problem.

# Example: edge-disjoint s-t paths
# Each edge capacity = 1; max flow = max edge-disjoint paths (Menger's theorem)

# Example: vertex-disjoint paths via vertex splitting
def vertex_split_network(n, s, t, adj_edges):
    """
    For vertex-disjoint paths:
    Node v -> split into v_in (v) and v_out (v + n)
    Internal edge: v_in -> v_out with capacity 1
    Original edge (u,v) -> u_out -> v_in with capacity infinity (or 1)
    """
    # 2n nodes: [0..n-1] = v_in, [n..2n-1] = v_out
    from collections import defaultdict
    graph = defaultdict(list)
    
    def add_edge(u, v, cap):
        graph[u].append([v, cap, 0, len(graph[v])])
        graph[v].append([u,   0, 0, len(graph[u]) - 1])
    
    # Internal vertex edges
    for v in range(n):
        if v == s or v == t:
            add_edge(v, v + n, 10**9)  # source/sink unconstrained
        else:
            add_edge(v, v + n, 1)      # capacity-1 internal edge
    
    # Graph edges: u_out -> v_in
    for (u, v) in adj_edges:
        add_edge(u + n, v, 1)  # u_out -> v_in
        add_edge(v + n, u, 1)  # v_out -> u_in (undirected)
    
    return graph, s, t + n, 2 * n  # new source=s_in, new sink=t_out

Is O(V²E) the best possible? This is an open question for general networks. In practice, push-relabel algorithms (Goldberg-Tarjan, 1988) achieve O(V²√E) and are often faster than Dinic's on dense graphs. The highest-label push-relabel variant is the practical choice for the hardest instances.

Push-relabel intuition: Instead of augmenting along entire paths, push-relabel maintains a "preflow" (where vertices can temporarily have excess inflow) and uses a height function to label vertices. Flow is pushed downhill (from higher to lower labels). When a vertex has excess but no valid downhill neighbor, it is relabeled (raised). The invariant ensures progress. This local, edge-by-edge update avoids the global BFS of Dinic's.

Minimum cost flow: Each edge has both a capacity and a cost per unit of flow. Goal: find a max flow of minimum total cost (or route a given amount of flow at minimum cost). Applications include optimal transportation, scheduling with costs, and min-cost bipartite matching. Algorithms: Successive Shortest Paths (repeatedly find min-cost augmenting path using Bellman-Ford or SPFA), Minimum Mean Cycle Canceling, or Hungarian algorithm as a special case.

Multi-commodity flow: Multiple source-sink pairs each need to route their own commodity through a shared network with shared capacities. Unlike single-commodity flow, the LP relaxation may have fractional optima, and the integer version is NP-hard in general. However, for two commodities in planar graphs it remains polynomial. This is a deep and active research area.

Parametric max-flow: How does the max flow change as a single capacity parameter varies? Gallo-Grigoriadis-Tarjan (1989) solve all parametric instances in the same time as a single max-flow, using a clever sensitivity structure. Used heavily in computer vision (graph cuts with parameter sweep for regularization strength).

# Pseudocode for push-relabel (FIFO variant)
#
# height[v]: label (integer, represents height in residual graph)
# excess[v]: excess flow at vertex v
# Invariant: if f(u,v) < c(u,v), then height[u] <= height[v] + 1
#
# PUSH(u, v):
#   delta = min(excess[u], c_f(u,v))
#   send delta units from u to v
#   excess[u] -= delta; excess[v] += delta
#
# RELABEL(u):
#   height[u] = 1 + min(height[v] for v in residual neighbors of u)
#
# GENERIC-PUSH-RELABEL:
#   height[s] = n; height[v] = 0 for v != s
#   Saturate all edges out of s (push s's initial excess)
#   while there exists u != s,t with excess[u] > 0:
#       if there exists admissible edge (u,v) [height[u] = height[v]+1]:
#           PUSH(u, v)
#       else:
#           RELABEL(u)