Network Flow

from collections import deque

def edmonds_karp(n, source, sink, capacity):
    """
    Edmonds-Karp max flow (Ford-Fulkerson + BFS shortest augmenting path).
    n: number of vertices (0-indexed).
    capacity: 2D list, capacity[u][v] = c(u,v).
    Returns (max_flow_value, flow_matrix).
    Time: O(V * E^2).
    """
    flow = [[0]*n for _ in range(n)]   # current flow on each edge

    def bfs_find_path():
        """Find shortest augmenting path s->t in residual graph via BFS."""
        parent = [-1] * n
        parent[source] = source
        queue  = deque([source])
        while queue and parent[sink] == -1:
            u = queue.popleft()
            for v in range(n):
                # Residual capacity: c_f(u,v) = capacity[u][v] - flow[u][v]
                if parent[v] == -1 and capacity[u][v] - flow[u][v] > 0:
                    parent[v] = u
                    queue.append(v)
        return parent if parent[sink] != -1 else None

    max_flow = 0
    while True:
        parent = bfs_find_path()
        if parent is None:
            break                          # no augmenting path: done

        # Find bottleneck capacity along the path
        path_flow = float('inf')
        v = sink
        while v != source:
            u = parent[v]
            path_flow = min(path_flow, capacity[u][v] - flow[u][v])
            v = u

        # Augment flow along the path
        v = sink
        while v != source:
            u = parent[v]
            flow[u][v] += path_flow
            flow[v][u] -= path_flow        # update reverse edge (skew symmetry)
            v = u

        max_flow += path_flow

    return max_flow, flow


def min_cut_edges(n, source, capacity, flow):
    """
    Find the minimum cut edges from a max flow solution.
    BFS on residual graph: S = vertices reachable from source.
    Min-cut edges: (u,v) with u in S, v not in S, capacity[u][v] > 0.
    """
    visited = [False] * n
    queue   = deque([source])
    visited[source] = True
    while queue:
        u = queue.popleft()
        for v in range(n):
            if not visited[v] and capacity[u][v] - flow[u][v] > 0:
                visited[v] = True
                queue.append(v)
    # Min-cut edges cross from S (visited) to T (not visited)
    return [(u,v) for u in range(n) for v in range(n)
            if visited[u] and not visited[v] and capacity[u][v] > 0]

def bipartite_max_matching(left_n, right_n, edges):
    """
    Maximum bipartite matching via max flow (Edmonds-Karp).
    left_n, right_n: sizes of left and right sets.
    edges: list of (l, r) pairs (0-indexed within each side).
    Returns (matching_size, matched_pairs).

    Network: s -> L_i (cap 1), L_i -> R_j (cap inf), R_j -> t (cap 1).
    Max flow = maximum matching size (by flow integrality).
    O(V*E) time: O(V) augmentations (one per matched pair), O(E) per BFS.
    """
    # Node numbering: 0=source, 1..left_n = left, left_n+1..left_n+right_n = right, last=sink
    n = left_n + right_n + 2
    s = 0
    t = n - 1
    INF = float('inf')
    cap = [[0]*n for _ in range(n)]

    # s -> each left vertex (cap 1)
    for i in range(left_n):
        cap[s][1 + i] = 1

    # left -> right edges (cap inf — won't be bottleneck)
    for l, r in edges:
        cap[1 + l][left_n + 1 + r] = INF

    # each right vertex -> t (cap 1)
    for j in range(right_n):
        cap[left_n + 1 + j][t] = 1

    match_size, flow = edmonds_karp(n, s, t, cap)

    # Extract matching from flow
    matched = []
    for l in range(left_n):
        for r in range(right_n):
            if flow[1+l][left_n+1+r] > 0:
                matched.append((l, r))

    return match_size, matched


def min_vertex_cover(left_n, right_n, edges):
    """
    Minimum vertex cover in bipartite graph = max matching (König's theorem).
    Returns (cover_size, left_in_cover, right_in_cover).
    Uses the König construction: from max flow, find min cut,
    then identify cover as (T ∩ L) ∪ (S ∩ R).
    """
    n = left_n + right_n + 2
    s, t = 0, n - 1
    INF = float('inf')
    cap = [[0]*n for _ in range(n)]
    for i in range(left_n):     cap[s][1+i] = 1
    for l, r in edges:          cap[1+l][left_n+1+r] = INF
    for j in range(right_n):    cap[left_n+1+j][t] = 1

    _, flow = edmonds_karp(n, s, t, cap)
    cut_edges = min_cut_edges(n, s, cap, flow)

    # König's theorem: cover = vertices incident to min-cut edges
    # S = reachable from s in residual, T = rest
    from collections import deque
    visited = [False]*n
    q = deque([s]); visited[s] = True
    while q:
        u = q.popleft()
        for v in range(n):
            if not visited[v] and cap[u][v]-flow[u][v] > 0:
                visited[v]=True; q.append(v)
    # Left in cover: left vertices in T (not reachable)
    left_cover  = [i for i in range(left_n) if not visited[1+i]]
    # Right in cover: right vertices in S (reachable)
    right_cover = [j for j in range(right_n) if visited[left_n+1+j]]
    return len(left_cover)+len(right_cover), left_cover, right_cover

def ford_fulkerson_dfs(n, source, sink, capacity):
    """
    Ford-Fulkerson with DFS augmentation (can be very slow!).
    Demonstrates why DFS path choice is dangerous.
    Same network as Edmonds-Karp but uses DFS to find augmenting paths.
    On the adversarial 4-node graph with capacities 10^9, this
    takes 2*10^9 iterations where Edmonds-Karp takes 2.
    """
    flow = [[0]*n for _ in range(n)]

    def dfs_find_path():
        visited = [False]*n
        parent  = [-1]*n
        stack   = [source]
        visited[source] = True
        parent[source]  = source
        while stack:
            u = stack.pop()
            if u == sink: return parent
            for v in range(n):
                if not visited[v] and capacity[u][v]-flow[u][v] > 0:
                    visited[v] = True
                    parent[v]  = u
                    stack.append(v)
        return None

    total, iterations = 0, 0
    while True:
        parent = dfs_find_path()
        if parent is None: break
        # Bottleneck
        pf, v = float('inf'), sink
        while v != source:
            u = parent[v]
            pf = min(pf, capacity[u][v] - flow[u][v])
            v = u
        # Augment
        v = sink
        while v != source:
            u = parent[v]
            flow[u][v] += pf; flow[v][u] -= pf
            v = u
        total += pf; iterations += 1
    return total, iterations

# Adversarial graph: 4 nodes, capacities that cause 2*10^9 iterations
# with DFS but only 2 with BFS (Edmonds-Karp)
cap_adversarial = [
    [0, 10**9, 10**9, 0],
    [0, 0,     1,     10**9],
    [0, 1,     0,     10**9],
    [0, 0,     0,     0],
]
# Edmonds-Karp handles this in 2 augmentations
ek_flow, _ = edmonds_karp(4, 0, 3, cap_adversarial)
print(f"Edmonds-Karp max flow: {ek_flow}")  # 2*10^9, 2 iterations