Advanced graphs & flows

1 · Floyd-Warshall — all-pairs shortest paths

Picture the little map your phone draws before you've even typed a destination. A ride-hailing app doesn't just want "how far is the driver from this passenger" — it wants a table: every driver's travel time to every waiting rider, all at once, so it can assign the whole fleet in one shot. Run Dijkstra once and you get a single row of that table. You need the whole grid. That "distance between everybody and everybody" problem is exactly what Floyd-Warshall was built for, and it solves it with three nested loops so plain you could write them from memory.

Here's the everyday analogy that makes the algorithm click. Imagine you know the direct driving time between certain pairs of towns, but not all — some journeys only work if you connect through a hub. Floyd-Warshall asks one patient question, over and over: "If I'm allowed to pass through town k on my way, does that give me a shorter trip between i and j than the best I've found so far?" It starts with no hubs allowed (direct edges only), then permits town 0 as a stopover, then towns 0 and 1, and so on — each round unlocking one more possible connection point. After every town has had its turn as a permitted middleman, every entry dist[i][j] holds the true cheapest route. That is precisely what "intermediate vertices numbered < k" means: the set of stopovers you've unlocked so far. Because the answer for the current k is built on the answers already banked for smaller k, this is really dynamic programming wearing a matrix.

The common mistake is reaching for it too eagerly. Floyd-Warshall is O(V³) — cubic in the number of vertices — so on a graph of 10,000 nodes that's a trillion operations and it simply will not finish in any reasonable time. It shines on dense or small graphs (a few hundred vertices), or when you genuinely need the full matrix. If you have a big sparse graph and only a handful of sources, run Dijkstra from each of those instead. The second thing people forget is its two quiet superpowers: it happily eats negative edge weights (which Dijkstra chokes on), and it doubles as a negative-cycle detector. If, after the algorithm runs, any vertex can improve on its own distance-to-itself of zero, there must be a loop you can walk around forever shedding cost — that's why the check is simply "did any dist[i][i] drop below zero?"

Now the precise version. Dijkstra answers "shortest path from one source." When you need the distance between every pair of vertices — routing tables, a distance matrix for clustering, transitive closure — Floyd-Warshall does it in one triple loop. The idea: dist[i][j] using only intermediate vertices numbered < k, grown one k at a time. It handles negative edges (unlike Dijkstra) and detects a negative cycle when any dist[i][i] goes below zero.

INF = float("inf")

def floyd_warshall(n, edges):
    """All-pairs shortest paths. edges: list of (u, v, w), directed.
    Returns dist matrix; dist[i][j] = cheapest i->j (INF if unreachable)."""
    dist = [[INF] * n for _ in range(n)]
    for i in range(n):
        dist[i][i] = 0
    for u, v, w in edges:
        dist[u][v] = min(dist[u][v], w)      # keep the cheapest parallel edge
    # k is the "highest intermediate vertex allowed" — grow it one at a time
    for k in range(n):
        dk = dist[k]
        for i in range(n):
            dik = dist[i][k]
            if dik == INF:                   # no path i..k, skip the whole row
                continue
            di = dist[i]
            for j in range(n):
                through_k = dik + dk[j]
                if through_k < di[j]:
                    di[j] = through_k
    return dist

def negative_cycle(dist):
    """A negative cycle exists iff some vertex can improve its own 0 distance."""
    return any(dist[i][i] < 0 for i in range(len(dist)))

# 4 nodes; note the cheaper 0->1->2 (=3) beats direct 0->2 (=10)
edges = [(0, 1, 1), (1, 2, 2), (0, 2, 10), (2, 3, 1), (3, 1, 4)]
d = floyd_warshall(4, edges)
print("0->2 :", d[0][2])          # 3, via node 1
print("0->3 :", d[0][3])          # 4
print("negative cycle?", negative_cycle(d))

0->2 : 3
0->3 : 4
negative cycle? False

2 · A* — heuristic (informed) search

Think about how you actually walk to a café you've never been to. You don't fan out equally in all directions like a spilled drink spreading across a table — you head roughly toward where you know the café is, and only detour when a street dead-ends. That instinct, "prefer the direction that seems to close the gap," is the entire idea behind A*. Dijkstra, by contrast, is the spilled-drink search: it expands outward by distance-from-start alone, exploring plenty of ground in the wrong direction because it has no sense of where the goal is. On a game map or a road network that difference can be the gap between checking a few hundred cells and checking a few hundred thousand.

The trick A* adds is a heuristic — a cheap guess, h(n), of how much cost still remains from node n to the goal. Instead of ordering its to-do list by "cost so far" (g) like Dijkstra, it orders by f = g + h: the real cost already spent plus the estimated cost still to go. That single addition biases exploration toward the goal. The guess doesn't have to be exact — on a grid, the straight-line "Manhattan" distance (how many steps if there were no walls) is a fine estimate you can compute in one line.

But there's a rule you cannot break, and it's where nearly everyone stumbles: the heuristic must be admissible, meaning it never overestimates the true remaining cost. Think of it as an optimistic-but-honest advisor: it may say "you're closer than you think," never "you're farther than you are." The moment your heuristic exaggerates the distance, A* may confidently prune away a route that was actually cheaper — because it believed that route led somewhere far — and hand you a suboptimal path while looking fast and smug about it. The classic mistake is to scale the heuristic up ("multiply Manhattan by 2 so it searches harder") precisely because it feels more aggressive; it is, and that's the bug. Two useful anchors: with h = 0 the guess is uselessly conservative and A* degrades exactly into Dijkstra (correct, just slower); with a perfect heuristic it walks straight to the goal expanding almost nothing.

Now the precise statement. Dijkstra expands nodes purely by distance-from-start. A* adds a heuristic h(n) — an estimate of the remaining cost to the goal — and orders the frontier by f(n) = g(n) + h(n). If h is admissible (never overestimates) and consistent, A* is optimal and expands no more nodes than Dijkstra — usually far fewer. With h = 0 it is Dijkstra; with a perfect heuristic it walks straight to the goal.

import heapq

def astar(graph, start, goal, h):
    """A* on a weighted graph. graph: {node: [(nbr, cost), ...]}.
    h(n) = admissible heuristic (never overestimates the true remaining cost).
    Returns (path, cost) or (None, inf)."""
    # frontier holds (f = g + h, g, node); g is the real cost from start
    frontier = [(h(start), 0, start)]
    best_g = {start: 0}
    parent = {start: None}
    while frontier:
        f, g, node = heapq.heappop(frontier)
        if node == goal:
            path = []
            while node is not None:
                path.append(node)
                node = parent[node]
            return path[::-1], g
        if g > best_g.get(node, float("inf")):
            continue                          # stale heap entry, skip
        for nbr, cost in graph[node]:
            ng = g + cost
            if ng < best_g.get(nbr, float("inf")):
                best_g[nbr] = ng
                parent[nbr] = node
                heapq.heappush(frontier, (ng + h(nbr), ng, nbr))
    return None, float("inf")

# A 3x3 grid; move cost 1; heuristic = Manhattan distance (admissible on a grid)
def grid_graph(n):
    g = {}
    for r in range(n):
        for c in range(n):
            g[(r, c)] = [((r+dr, c+dc), 1) for dr, dc in ((1,0),(-1,0),(0,1),(0,-1))
                         if 0 <= r+dr < n and 0 <= c+dc < n]
    return g

goal = (2, 2)
manhattan = lambda p: abs(p[0]-goal[0]) + abs(p[1]-goal[1])
path, cost = astar(grid_graph(3), (0, 0), goal, manhattan)
print("path :", path)
print("cost :", cost)             # 4 (Manhattan distance corner to corner)

path : [(0, 0), (0, 1), (0, 2), (1, 2), (2, 2)]
cost : 4

3 · Strongly connected components — Kosaraju

Here's a bug that ships in real systems more often than anyone admits: the circular import. Module A imports B, B imports C, and C — three files away — imports A. The build seems fine until one day it deadlocks or throws a baffling half-initialized error, and you spend an afternoon drawing arrows on a whiteboard trying to find the loop. What you're hunting for, without knowing its name, is a strongly connected component: a cluster of things that can all reach each other, round and round. Package managers, deadlock detectors, and "which services are mutually entangled" dashboards all run this exact analysis for you.

The key intuition is the difference between "connected" and strongly connected in a directed graph, where arrows point one way. Two towns joined by a one-way street are connected on a map but not mutually reachable — you can drive A→B but not B→A. An SCC is a set of vertices where you genuinely can get from any member to any other and back: a knot of two-way reachability. Every such knot, collapsed to a single dot, turns a tangled directed graph into a clean acyclic one — which is why SCC-finding is the standard first step before you can topologically schedule anything cyclic (you'll build that condensation in the exercises).

Kosaraju's method sounds almost too cute to work, so here's the plain-language why. First you DFS the whole graph and write down the order in which vertices finish (the moment all their descendants are explored). Then — the clever move — you reverse every edge and DFS again, but starting from the vertex that finished last, then the next-latest, and so on. Reversing the edges means you can now only travel "backwards"; a vertex you reach in this second pass must be one that could reach you in the original — and combined with the finish-order trick, each tree the second DFS grows is exactly one SCC. The common mistake is confusing SCCs (directed, two-way reachability) with plain connected components (undirected) — using the wrong one silently gives wrong clusters. Everything runs in two linear passes.

Formally: in a directed graph, a strongly connected component (SCC) is a maximal set of vertices where every vertex can reach every other. SCCs reveal cyclic dependency clusters (mutually-recursive modules, deadlock cycles). Kosaraju is the easiest to reason about: DFS the graph recording finish order, reverse every edge, then DFS in decreasing finish order — each tree is one SCC. Two linear passes, O(V+E).

from collections import defaultdict

def kosaraju(n, adj):
    """Strongly connected components via Kosaraju's two-pass DFS.
    adj: {u: [v, ...]} directed. Returns list of components (lists of nodes)."""
    visited = [False] * n
    order = []                                # finish-time order (post-order)

    def dfs1(u):                              # iterative to dodge recursion limits
        stack = [(u, iter(adj[u]))]
        visited[u] = True
        while stack:
            node, it = stack[-1]
            for v in it:
                if not visited[v]:
                    visited[v] = True
                    stack.append((v, iter(adj[v])))
                    break
            else:
                order.append(node)            # all children done -> record finish
                stack.pop()

    for u in range(n):
        if not visited[u]:
            dfs1(u)

    radj = defaultdict(list)                   # reverse every edge
    for u in range(n):
        for v in adj[u]:
            radj[v].append(u)

    comp = [-1] * n
    def dfs2(u, cid):
        stack = [u]
        comp[u] = cid
        while stack:
            node = stack.pop()
            for v in radj[node]:
                if comp[v] == -1:
                    comp[v] = cid
                    stack.append(v)

    cid = 0
    for u in reversed(order):                  # decreasing finish time
        if comp[u] == -1:
            dfs2(u, cid); cid += 1
    groups = defaultdict(list)
    for node, c in enumerate(comp):
        groups[c].append(node)
    return list(groups.values())

# 0->1->2->0 is one cycle (an SCC); 3 hangs off it alone
adj = {0: [1], 1: [2], 2: [0, 3], 3: []}
print(sorted(sorted(c) for c in kosaraju(4, adj)))   # [[0,1,2],[3]]

[[0, 1, 2], [3]]

4 · Strongly connected components — Tarjan

Kosaraju works, but it walks the graph twice and builds a whole reversed copy — and on a graph with tens of millions of edges, "twice" and "a whole copy" are exactly the words that make a memory dashboard turn red. Tarjan's algorithm answers the same question in a single pass with no reversed graph, which is why it's the version interviewers love and the version production code usually ships. The catch is that it's cleverer, so it earns a slower read.

The idea is a piece of bookkeeping you can picture. As the DFS visits vertices it stamps each with a index — a discovery number, "you're the 5th vertex I met." It also keeps a running low-link for each vertex: the smallest discovery number it can still climb back up to by following edges within its own subtree, including at most one edge that loops back to an ancestor. Think of the index as "when I was born" and the low-link as "the oldest ancestor I can still phone." If a vertex's low-link equals its own index, that means nobody in its subtree found a way back to anything older than it — so it's the root of a self-contained knot, and everything piled above it on the DFS stack since it was discovered is precisely its SCC. You pop that group off and record it.

The mistake to guard against is mismanaging the stack: only vertices currently on the stack (still "open") count when you update a low-link through a back-edge. A vertex that already got assigned to a finished SCC is off-limits — folding its index in would wrongly glue two separate knots together. That single on_stack flag is the whole difference between correct and subtly-wrong Tarjan. Everything else is the same linear budget as Kosaraju; you just pay it once instead of twice.

Precisely: Tarjan finds SCCs in a single DFS using a discovery index and a low-link (the lowest index reachable from a vertex’s subtree via at most one back-edge). When a vertex’s low == index, it is the root of an SCC and everything above it on the stack forms the component. Same O(V+E) as Kosaraju but one pass and no reversed graph.

def tarjan_scc(n, adj):
    """Tarjan's single-pass SCC. Returns list of components.
    index = DFS discovery time; low = lowest index reachable from the subtree."""
    index = [-1] * n
    low = [0] * n
    on_stack = [False] * n
    stack, comps = [], []
    counter = [0]

    def strongconnect(v):
        # explicit-stack DFS carrying an iterator per frame
        work = [(v, iter(adj[v]))]
        index[v] = low[v] = counter[0]; counter[0] += 1
        stack.append(v); on_stack[v] = True
        while work:
            node, it = work[-1]
            advanced = False
            for w in it:
                if index[w] == -1:
                    index[w] = low[w] = counter[0]; counter[0] += 1
                    stack.append(w); on_stack[w] = True
                    work.append((w, iter(adj[w])))
                    advanced = True
                    break
                elif on_stack[w]:
                    low[node] = min(low[node], index[w])   # back-edge
            if advanced:
                continue
            work.pop()
            if work:                                        # fold child low into parent
                low[work[-1][0]] = min(low[work[-1][0]], low[node])
            if low[node] == index[node]:                    # node is an SCC root
                comp = []
                while True:
                    w = stack.pop(); on_stack[w] = False
                    comp.append(w)
                    if w == node:
                        break
                comps.append(comp)

    for v in range(n):
        if index[v] == -1:
            strongconnect(v)
    return comps

adj = {0: [1], 1: [2], 2: [0, 3], 3: []}
print(sorted(sorted(c) for c in tarjan_scc(4, adj)))   # [[0,1,2],[3]]

[[0, 1, 2], [3]]

5 · Max-flow and the min-cut duality

Turn on every tap in your house at once and watch what happens. The water doesn't gush infinitely — it's throttled by the narrowest pipe in the path from the street main to your faucet. That single skinniest segment sets the ceiling on how much can ever flow, no matter how wide everything else is. Networks of all kinds have this shape: data pushed through routers with bandwidth limits, cargo routed over roads with capacity caps, requests fanned across servers with connection ceilings. The question "how much can I push from source to sink?" is the max-flow problem, and the answer is always pinned by a bottleneck.

The algorithm's core move is disarmingly simple and one detail is genuinely surprising. Ford-Fulkerson keeps finding a path from source to sink that still has spare room — an augmenting path — and pushes as much as the tightest edge on it allows, then repeats until no path has any room left. The surprise is the residual back-edge: every time you push, say, 3 units along an edge, you also create a phantom reverse edge of capacity 3 pointing backward. Why? Because a greedy early choice can wedge flow into a route that blocks a better global arrangement, and the back-edge is a built-in "undo" — a later path can shove flow backward along it, effectively rerouting an earlier commitment. Without that undo, the method would be a mere heuristic that sometimes gets stuck; with it, Ford-Fulkerson is provably optimal. The everyday analogy: it's like letting yourself trade a seat you already claimed so two people can both end up happier.

Then comes one of the most beautiful facts in all of algorithms — the max-flow min-cut theorem. The maximum you can push equals the cheapest set of pipes you'd have to sever to fully disconnect source from sink. Flow and cut are two views of the same bottleneck: the most you can send is exactly the least you'd have to destroy to send nothing. That duality is why so many seemingly unrelated problems — image segmentation, project selection, "which team is the single point of failure" — quietly reduce to a flow computation. The common mistake is picking augmenting paths carelessly with plain DFS, which can crawl when capacities are huge; use BFS to always take the shortest augmenting path and the running time stops depending on the capacity numbers at all.

To state it exactly: a flow network has a source, a sink, and edges with capacities. Ford-Fulkerson repeatedly finds an augmenting path with spare capacity and pushes flow along it, adding a residual back-edge so later paths can "undo" earlier choices. Picking the path with BFS (shortest augmenting path) is Edmonds-Karp, which runs in O(V·E²) — polynomial and independent of capacity magnitudes. The max-flow min-cut theorem: the maximum flow equals the minimum total capacity you must cut to disconnect source from sink.

from collections import defaultdict, deque

class MaxFlow:
    """Edmonds-Karp: Ford-Fulkerson using BFS to find augmenting paths.
    Runs in O(V * E^2) — polynomial, unlike naive DFS Ford-Fulkerson."""
    def __init__(self, n):
        self.n = n
        self.cap = defaultdict(int)          # (u, v) -> residual capacity
        self.adj = defaultdict(set)          # neighbours in the residual graph

    def add_edge(self, u, v, c):
        self.cap[(u, v)] += c
        self.adj[u].add(v); self.adj[v].add(u)   # reverse edge starts at 0

    def _bfs(self, s, t, parent):
        parent.clear(); parent[s] = s
        q = deque([s])
        while q:
            u = q.popleft()
            for v in self.adj[u]:
                if v not in parent and self.cap[(u, v)] > 0:
                    parent[v] = u
                    if v == t:
                        return True
                    q.append(v)
        return False

    def max_flow(self, s, t):
        flow, parent = 0, {}
        while self._bfs(s, t, parent):
            # bottleneck = smallest residual capacity along the found path
            v, bottleneck = t, float("inf")
            while v != s:
                u = parent[v]
                bottleneck = min(bottleneck, self.cap[(u, v)])
                v = u
            v = t
            while v != s:                    # push flow, add residual back-edge
                u = parent[v]
                self.cap[(u, v)] -= bottleneck
                self.cap[(v, u)] += bottleneck
                v = u
            flow += bottleneck
        return flow

    def min_cut(self, s):
        """After max_flow, vertices reachable from s in the residual graph are
        the source side of a minimum cut (max-flow min-cut theorem)."""
        seen = {s}; q = deque([s])
        while q:
            u = q.popleft()
            for v in self.adj[u]:
                if v not in seen and self.cap[(u, v)] > 0:
                    seen.add(v); q.append(v)
        return seen

mf = MaxFlow(4)
for u, v, c in [(0,1,3),(0,2,2),(1,2,1),(1,3,2),(2,3,3)]:
    mf.add_edge(u, v, c)
print("max flow 0->3 :", mf.max_flow(0, 3))    # 5 (both source edges saturate)
print("source side of min cut :", sorted(mf.min_cut(0)))

max flow 0->3 : 5
source side of min cut : [0]

6 · Bipartite matching

You've almost certainly felt this problem as a manager staring at a whiteboard. Five engineers, five open tasks, and everyone can only do some of the tasks based on skill. Assign greedily — hand each person the first thing they're able to do — and you routinely paint yourself into a corner: the last engineer is qualified only for a task you already gave away, and now someone sits idle for no good reason. Sound familiar? It's the same shape as pairing drivers to deliveries, students to project slots, or ad impressions to advertisers. The goal is a maximum matching: the largest possible set of pairings where nobody is double-booked.

The structure has a name — a bipartite graph — and it just means the world splits into two sides (people on the left, tasks on the right) with connections only crossing the gap, never within a side. Two engineers never "match" each other; two tasks never do. What makes maximum matching more than greedy is the idea of an augmenting path, and the analogy is exactly the manager's fix: when an engineer's only option is already taken, you don't give up — you ask whether that task's current owner could be shuffled to a different task they're also qualified for, freeing the slot. If such a chain of "you move here so I can take your old spot" exists and ends on a genuinely free task, following it makes one more pairing than before. Repeat until no augmenting path remains and you're provably optimal.

The common mistake is trusting that greedy assignment is "good enough" — it can leave people unmatched who could have been placed, and the loss is invisible unless you know to look for the reroute. Notice this is the same rerouting idea as the residual back-edge from max-flow, and that's no accident: maximum bipartite matching is literally a max-flow problem with every edge given capacity 1. Kuhn's algorithm just expresses that flow directly in terms of augmenting paths, which is simpler to code when all you need is the matching.

Precisely: a bipartite graph splits into two sides (workers/tasks, applicants/jobs) with edges only across. Maximum matching pairs as many as possible with no vertex used twice. It is a special case of max-flow (unit capacities), but Kuhn’s algorithm expresses it directly: for each left vertex, try to find an augmenting path that either lands on a free right vertex or bumps an existing match to a different slot. O(V·E).

def hopcroft_style_matching(left, adj):
    """Maximum bipartite matching via repeated augmenting paths (Hungarian /
    Kuhn's algorithm). left = list of left-side nodes; adj[l] = allowed rights.
    Returns dict right -> left. O(V * E)."""
    match_r = {}                              # right node -> matched left node

    def try_augment(u, seen):
        for v in adj[u]:
            if v in seen:
                continue
            seen.add(v)
            # v is free, or its current partner can be rematched elsewhere
            if v not in match_r or try_augment(match_r[v], seen):
                match_r[v] = u
                return True
        return False

    for u in left:
        try_augment(u, set())
    return match_r

# workers -> tasks they can do; find the largest set of disjoint assignments
adj = {"w1": ["t1", "t2"], "w2": ["t1"], "w3": ["t2", "t3"]}
m = hopcroft_style_matching(["w1", "w2", "w3"], adj)
pairs = sorted((l, r) for r, l in m.items())
print("assignments :", pairs)
print("matching size :", len(m))              # 3 — everyone gets a task

assignments : [('w1', 't2'), ('w2', 't1'), ('w3', 't3')]
matching size : 3

7 · Articulation points and bridges

Every network diagram hides a nightmare question: which single thing, if it dies, splits us in two? One switch whose failure isolates a whole wing of the building. One undersea cable whose cut severs a country from the internet. One router that, when it reboots, strands half your services from the other half. Postmortems are full of these — a component nobody flagged as special turns out to have been the only bridge between two islands of the system. Finding those single points of failure before the outage is exactly what this section's two definitions give you.

The two shapes are worth separating clearly. An articulation point (or cut vertex) is a node whose removal disconnects the graph — pull it out and the network falls into pieces. A bridge is an edge with the same property — cut that one link and connectivity shatters. The intuition for spotting them: a piece of the graph is safely attached only if it has a second way home. Picture a subtree hanging off vertex u during a DFS. If anything down in that subtree has a back-edge climbing to u or higher — an alternate route around u — then removing u doesn't strand it; there's another path. But if the entire subtree can reach no higher than u, then u is the only door, and u is an articulation point.

That "can it reach higher?" test is measured with the same low-link machinery from Tarjan: each vertex tracks the oldest ancestor its subtree can still climb back to. The comparison is a hair's-breadth apart for the two cases and it's the classic place people slip: low[child] ≥ disc[u] (child can't get above u) makes u a cut vertex, while the strict low[child] > disc[u] (child can't even reach u itself by another route) makes the connecting edge a bridge. The other subtlety is the DFS root: it has no parent, so the "some child is trapped" rule doesn't apply — the root is an articulation point only if it has two or more separate children in the DFS tree. Mishandle that special case and you'll wrongly flag or miss the root every time. One DFS finds all of it.

Formally: an articulation point (cut vertex) is a node whose removal disconnects the graph; a bridge is an edge whose removal does. These are the single points of failure in a network topology. A single DFS with discovery times and low-links finds both: a non-root vertex u is an articulation point if some child’s subtree cannot reach above u (low[child] ≥ disc[u]); the edge is a bridge if it strictly cannot (low[child] > disc[u]). O(V+E).

def articulation_and_bridges(n, adj):
    """Find articulation points (cut vertices) and bridges (cut edges) in one
    DFS using discovery times and low-links. Undirected graph. O(V + E)."""
    disc = [-1] * n
    low = [0] * n
    timer = [0]
    aps = set()
    bridges = []

    def dfs(root):
        stack = [(root, -1, iter(adj[root]))]
        disc[root] = low[root] = timer[0]; timer[0] += 1
        child_count = 0                       # children of the DFS root
        while stack:
            u, parent, it = stack[-1]
            for w in it:
                if disc[w] == -1:
                    if parent == -1:
                        child_count += 1
                    disc[w] = low[w] = timer[0]; timer[0] += 1
                    stack.append((w, u, iter(adj[w])))
                    break
                elif w != parent:
                    low[u] = min(low[u], disc[w])   # back-edge
            else:
                stack.pop()
                if stack:
                    pu = stack[-1][0]
                    low[pu] = min(low[pu], low[u])
                    # non-root u's parent is an AP if u can't reach above it
                    if stack[-1][1] != -1 and low[u] >= disc[pu]:
                        aps.add(pu)
                    if low[u] > disc[pu]:           # nothing below reaches pu -> bridge
                        bridges.append(tuple(sorted((pu, u))))
        if child_count > 1:                          # root is an AP iff 2+ children
            aps.add(root)

    for v in range(n):
        if disc[v] == -1:
            dfs(v)
    return sorted(aps), sorted(bridges)

# path 0-1-2 with a triangle 2-3-4-2; node 2 and edge (1,2) are cut points
adj = {0: [1], 1: [0, 2], 2: [1, 3, 4], 3: [2, 4], 4: [2, 3]}
aps, bridges = articulation_and_bridges(5, adj)
print("articulation points :", aps)     # [1, 2]
print("bridges :", bridges)             # [(0, 1), (1, 2)]

articulation points : [1, 2]
bridges : [(0, 1), (1, 2)]

def transitive_closure(n, edges):
    """Reachability matrix: reach[i][j] = can i reach j? Floyd-Warshall variant."""
    reach = [[False] * n for _ in range(n)]
    for i in range(n):
        reach[i][i] = True
    for u, v in edges:
        reach[u][v] = True
    for k in range(n):
        for i in range(n):
            if reach[i][k]:
                for j in range(n):
                    if reach[k][j]:
                        reach[i][j] = True
    return reach

edges = [(0, 1), (1, 2), (3, 3)]
r = transitive_closure(4, edges)
print("0 reaches 2 :", r[0][2])    # True (0->1->2)
print("2 reaches 0 :", r[2][0])    # False
print("3 reaches 1 :", r[3][1])    # False

INF = float("inf")

def floyd_with_paths(n, edges):
    """All-pairs distances plus path reconstruction via a next-hop matrix."""
    dist = [[INF] * n for _ in range(n)]
    nxt = [[None] * n for _ in range(n)]
    for i in range(n):
        dist[i][i] = 0
    for u, v, w in edges:
        if w < dist[u][v]:
            dist[u][v] = w
            nxt[u][v] = v
    for k in range(n):
        for i in range(n):
            if dist[i][k] == INF:
                continue
            for j in range(n):
                if dist[i][k] + dist[k][j] < dist[i][j]:
                    dist[i][j] = dist[i][k] + dist[k][j]
                    nxt[i][j] = nxt[i][k]      # first hop toward j goes via k's path
    return dist, nxt

def reconstruct(nxt, u, v):
    if nxt[u][v] is None:
        return []
    path = [u]
    while u != v:
        u = nxt[u][v]
        path.append(u)
    return path

edges = [(0, 1, 1), (1, 2, 2), (0, 2, 10), (2, 3, 1)]
dist, nxt = floyd_with_paths(4, edges)
print("0->3 cost :", dist[0][3])            # 4
print("0->3 path :", reconstruct(nxt, 0, 3))   # [0, 1, 2, 3]

import heapq

def astar_vs_dijkstra(graph, start, goal, h):
    """Return (nodes_expanded_astar, nodes_expanded_dijkstra) to show A* explores
    fewer nodes when the heuristic is informative. h=0 makes A* == Dijkstra."""
    def search(heur):
        frontier = [(heur(start), 0, start)]
        best = {start: 0}
        expanded = 0
        while frontier:
            f, g, node = heapq.heappop(frontier)
            if g > best.get(node, float("inf")):
                continue
            expanded += 1
            if node == goal:
                return expanded
            for nbr, c in graph[node]:
                ng = g + c
                if ng < best.get(nbr, float("inf")):
                    best[nbr] = ng
                    heapq.heappush(frontier, (ng + heur(nbr), ng, nbr))
        return expanded
    return search(h), search(lambda _n: 0)

def grid(n):
    g = {}
    for r in range(n):
        for c in range(n):
            g[(r, c)] = [((r+dr, c+dc), 1) for dr, dc in ((1,0),(-1,0),(0,1),(0,-1))
                         if 0 <= r+dr < n and 0 <= c+dc < n]
    return g

goal = (4, 4)
man = lambda p: abs(p[0]-goal[0]) + abs(p[1]-goal[1])
a, d = astar_vs_dijkstra(grid(5), (0, 0), goal, man)
print("A* expanded    :", a)
print("Dijkstra expanded:", d)
print("A* explored fewer:", a < d)

from collections import defaultdict

def condensation(n, adj):
    """Collapse each SCC into one super-node; the result is always a DAG.
    Uses Tarjan-style SCC ids, then dedupes cross-component edges."""
    index = [-1] * n; low = [0] * n; on = [False] * n
    stack, comp_id = [], [-1] * n
    counter = [0]; cid = [0]

    def sc(v):
        work = [(v, iter(adj[v]))]
        index[v] = low[v] = counter[0]; counter[0] += 1
        stack.append(v); on[v] = True
        while work:
            node, it = work[-1]
            adv = False
            for w in it:
                if index[w] == -1:
                    index[w] = low[w] = counter[0]; counter[0] += 1
                    stack.append(w); on[w] = True
                    work.append((w, iter(adj[w]))); adv = True; break
                elif on[w]:
                    low[node] = min(low[node], index[w])
            if adv:
                continue
            work.pop()
            if work:
                low[work[-1][0]] = min(low[work[-1][0]], low[node])
            if low[node] == index[node]:
                while True:
                    w = stack.pop(); on[w] = False; comp_id[w] = cid[0]
                    if w == node:
                        break
                cid[0] += 1

    for v in range(n):
        if index[v] == -1:
            sc(v)
    dag = defaultdict(set)
    for u in range(n):
        for v in adj[u]:
            if comp_id[u] != comp_id[v]:
                dag[comp_id[u]].add(comp_id[v])
    return comp_id, {k: sorted(v) for k, v in dag.items()}

adj = {0: [1], 1: [2], 2: [0, 3], 3: [4], 4: [3]}   # {0,1,2} and {3,4} are SCCs
comp_id, dag = condensation(5, adj)
print("component of each node :", comp_id)
print("condensed DAG edges    :", dict(dag))

from collections import defaultdict, deque

def project_selection(profits, prereqs):
    """Maximum-weight closure via min-cut: pick a subset of projects to maximise
    profit, where selecting a project forces selecting its prerequisites.
    Classic reduction of a scheduling problem to max-flow/min-cut."""
    n = len(profits)
    S, T = n, n + 1                            # super source / sink
    cap = defaultdict(int); adj = defaultdict(set)
    def edge(u, v, c):
        cap[(u, v)] += c; adj[u].add(v); adj[v].add(u)

    total_gain = 0
    for i, pr in enumerate(profits):
        if pr > 0:
            edge(S, i, pr); total_gain += pr   # profit if kept on source side
        elif pr < 0:
            edge(i, T, -pr)                     # cost if kept on source side
    INF = 10 ** 9
    for i, deps in enumerate(prereqs):
        for d in deps:
            edge(i, d, INF)                     # can't take i without prerequisite d

    def bfs(parent):
        parent.clear(); parent[S] = S; q = deque([S])
        while q:
            u = q.popleft()
            for v in adj[u]:
                if v not in parent and cap[(u, v)] > 0:
                    parent[v] = u
                    if v == T:
                        return True
                    q.append(v)
        return False

    flow, parent = 0, {}
    while bfs(parent):
        v, b = T, INF
        while v != S:
            b = min(b, cap[(parent[v], v)]); v = parent[v]
        v = T
        while v != S:
            u = parent[v]; cap[(u, v)] -= b; cap[(v, u)] += b; v = u
        flow += b
    return total_gain - flow                    # max profit = gains - min cut

# project 0 (+10) needs 1 (-4); project 2 (+3) needs 1 too -> take 0,1,2 = 9
profit = project_selection([10, -4, 3], [[1], [], [1]])
print("max profit :", profit)                  # 9

from collections import defaultdict, deque

class FlowNetwork:
    """Edmonds-Karp with capacity scaling awareness kept simple. Used here to
    solve a real assignment: route on-call requests to engineers with capacity
    and skill constraints, maximising covered shifts."""
    def __init__(self, n):
        self.n = n
        self.cap = defaultdict(int); self.adj = defaultdict(set)
    def add_edge(self, u, v, c):
        self.cap[(u, v)] += c; self.adj[u].add(v); self.adj[v].add(u)
    def max_flow(self, s, t):
        flow, parent = 0, {}
        while self._bfs(s, t, parent):
            v, b = t, float("inf")
            while v != s:
                b = min(b, self.cap[(parent[v], v)]); v = parent[v]
            v = t
            while v != s:
                u = parent[v]; self.cap[(u, v)] -= b; self.cap[(v, u)] += b; v = u
            flow += b
        return flow
    def _bfs(self, s, t, parent):
        parent.clear(); parent[s] = s; q = deque([s])
        while q:
            u = q.popleft()
            for v in self.adj[u]:
                if v not in parent and self.cap[(u, v)] > 0:
                    parent[v] = u
                    if v == t:
                        return True
                    q.append(v)
        return False

# 2 engineers (cap 2 shifts each) cover 3 shifts they are skilled for.
# Nodes: S=0, engineers 1..2, shifts 3..5, T=6
net = FlowNetwork(7)
S, T = 0, 6
net.add_edge(S, 1, 2); net.add_edge(S, 2, 2)     # per-engineer capacity
skills = {1: [3, 4], 2: [4, 5]}                   # which shifts each can take
for eng, shifts in skills.items():
    for sh in shifts:
        net.add_edge(eng, sh, 1)
for sh in (3, 4, 5):
    net.add_edge(sh, T, 1)                        # each shift needs 1 engineer
print("shifts covered :", net.max_flow(S, T))    # 3 — full coverage possible