ISAP, in depth: distance labels, advance and retreat, the gap heuristic and a tested max-flow implementation

Residual graphs and shortest augmenting paths

A flow network is a directed graph with a capacity on each arc, a source s and a sink t. A flow assigns each arc a value between zero and its capacity so that every node other than the source and sink has inflow equal to outflow. The maximum flow is the largest total amount that can leave the source.

Every augmenting-path algorithm works on the residual graph. For an arc u to v with capacity c carrying flow f, the residual graph has a forward arc with capacity c minus f (room to push more) and a backward arc v to u with capacity f (room to cancel flow already sent). A path from s to t in the residual graph is an augmenting path; pushing its bottleneck capacity increases the flow. When no augmenting path exists, the flow is maximum and the nodes reachable from s form one side of a minimum cut. The Ford-Fulkerson method and the max-flow min-cut theorem articles cover this foundation in detail.

The choice of which augmenting path to push is what separates algorithms. Edmonds-Karp always takes a shortest path, found by a fresh BFS for every augmentation, giving O(VE²). Dinic finds all shortest paths of the current length in one phase with one BFS plus blocking-flow DFS, giving O(V²E). ISAP also only uses shortest paths, but it never redoes the BFS: it repairs distance information locally as the residual graph changes.

Distance labels and admissible arcs

A distance labelling assigns each node an integer d(v) with two properties: d(t) = 0, and for every residual arc u to v, d(u) is at most d(v) + 1. Such a labelling is called valid, and it has a useful consequence: d(v) is a lower bound on the number of arcs on any residual path from v to t. Labels can only underestimate.

An arc u to v is admissible when it has residual capacity and d(u) = d(v) + 1. A path made only of admissible arcs descends one level per step, so its length equals d(s), which is the smallest possible. That is the whole trick: if you only walk admissible arcs, every augmenting path you find is a shortest one, which is what gives the shortest-augmenting-path bound, and you never need to search for it.

Because any simple path has at most n minus 1 arcs, where n is the number of nodes, a label of n or more on the source proves that no augmenting path exists. ISAP stops there.

Advance, augment, retreat

ISAP maintains a current node u, starting at the source, and the arcs that led to it. At each step it does exactly one of three things.

1. Advance. If u has an admissible arc u to v, record that arc as v's parent and move to v.
2. Augment. If u is the sink, walk the parent arcs back to the source, find the bottleneck residual capacity, subtract it along the path and add it to the reverse arcs, then restart at the source.
3. Retreat. If u has no admissible arc, relabel it: set d(u) to one plus the minimum label among nodes reachable by a residual arc from u (or n if there are none), then step back to u's parent. Relabelling keeps the labelling valid and makes at least one arc admissible if any residual arc exists.

Two optimisations turn this from correct into fast. The current-arc pointer remembers, for each node, where its last scan for an admissible arc stopped. An arc that is not admissible cannot become admissible again until u is relabelled, so the scan resumes from the pointer and only resets on relabel. The gap heuristic keeps a count of nodes at each label. If relabelling u empties its level k, no node above level k can reach the sink, because labels drop by at most one per residual arc and some node would have to sit at level k. The source is above that level, so the algorithm stops immediately instead of relabelling its way up to n one node at a time. In practice the gap check is often the difference between ISAP finishing quickly and spending most of its time proving there is nothing left.

A complete implementation

The implementation below stores arcs in flat arrays with a linked list per node, the layout most fast max-flow codes use. Arcs are added in pairs so the residual twin of arc e is e ^ 1. It is iterative, and it was checked against Edmonds-Karp on thousands of random graphs.

from collections import deque

INF = float("inf")


class ISAP:
    def __init__(self, n):
        self.n = n
        self.head = [-1] * n
        self.to, self.cap, self.nxt = [], [], []

    def add_edge(self, u, v, c):
        # Arcs are stored in pairs: e is u->v, e ^ 1 is its residual twin v->u.
        for a, b, cc in ((u, v, c), (v, u, 0)):
            self.to.append(b)
            self.cap.append(cc)
            self.nxt.append(self.head[a])
            self.head[a] = len(self.to) - 1

    def max_flow(self, s, t):
        if s == t:
            raise ValueError("source and sink must differ")
        n, to, cap, nxt, head = self.n, self.to, self.cap, self.nxt, self.head
        # 1. Exact labels: reverse BFS from the sink over residual arcs.
        d = [n] * n
        d[t] = 0
        q = deque([t])
        while q:
            v = q.popleft()
            e = head[v]
            while e != -1:
                u = to[e]
                if cap[e ^ 1] > 0 and d[u] == n:  # residual arc u->v
                    d[u] = d[v] + 1
                    q.append(u)
                e = nxt[e]
        if d[s] == n:
            return 0
        cnt = [0] * (n + 1)
        for x in d:
            cnt[x] += 1
        cur = head[:]            # current-arc pointer per node
        pe = [-1] * n            # arc used to enter each node on the path
        flow, u = 0, s
        while d[s] < n:
            if u == t:           # 2. augment along the stored path
                f, v = INF, t
                while v != s:
                    f = min(f, cap[pe[v]])
                    v = to[pe[v] ^ 1]
                v = t
                while v != s:
                    cap[pe[v]] -= f
                    cap[pe[v] ^ 1] += f
                    v = to[pe[v] ^ 1]
                flow += f
                u = s
                continue
            e = cur[u]
            while e != -1 and not (cap[e] > 0 and d[u] == d[to[e]] + 1):
                e = nxt[e]
            if e != -1:          # 3. advance on an admissible arc
                cur[u] = e
                pe[to[e]] = e
                u = to[e]
                continue
            m = n - 1            # 4. retreat: relabel u
            e = head[u]
            while e != -1:
                if cap[e] > 0:
                    m = min(m, d[to[e]])
                e = nxt[e]
            cnt[d[u]] -= 1
            if cnt[d[u]] == 0:   # gap: nothing above this level reaches t
                break
            d[u] = m + 1
            cnt[d[u]] += 1
            cur[u] = head[u]
            if u != s:
                u = to[pe[u] ^ 1]
        return flow

Usage is two lines: g = ISAP(n), then g.add_edge(u, v, c) for each arc and g.max_flow(s, t). For an undirected edge, add one pair and give both directions the capacity rather than adding two separate pairs.

Worked example: tracing a four-node network

Take the network in the diagram: s to a (3), s to b (2), a to b (1), a to t (2), b to t (3). The reverse BFS gives d(t) = 0, d(a) = d(b) = 1, d(s) = 2, so the level counts are one node at 0, two at 1 and one at 2. The trace below scans arcs in the order listed; the linked-list code visits them in reverse insertion order, which changes which path is found first but not the answer.

The answer is 5, and the residual graph confirms it: s can reach nothing, so the cut around s alone has capacity 3 + 2 = 5. Notice that the augmenting path lengths went 2, 2, 3, never decreasing, which is exactly the shortest-path guarantee at work.

Why it is O(V²E)

The bound follows the standard shortest-augmenting-path argument. Labels never decrease and never exceed n, so there are at most n relabels per node and O(n²) relabels in total; each relabel scans the node's arcs, so all relabels together cost O(nm), where m is the number of arcs. Each augmentation saturates at least one arc. After arc u to v is saturated, it can only carry flow again once flow is pushed back along v to u, which needs d(v) = d(u) + 1, so d(u) must rise by at least 2 between consecutive saturations. Each arc is therefore saturated O(n) times, giving O(nm) augmentations, each costing O(n) to walk the path. Advances are paid for by augmentations and retreats, and the current-arc pointer keeps scans to O(nm) overall. Total: O(n²m), the same worst case as Dinic.

Worst-case bounds are pessimistic for both; on real networks constant factors decide.

ISAP versus Dinic and push-relabel

ISAP and push-relabel share the distance-label idea; the difference is that ISAP moves flow only along complete paths, so it never has excess stranded on interior nodes, which makes it easier to reason about and to extract a flow decomposition from. Dinic recomputes every label with a BFS each phase, while ISAP updates only the nodes it gets stuck at; when many labels change at once, Dinic's full rebuild can be cheaper. Benchmark both on your own graph family before choosing.

Operating it in real code

• Memory layout. Flat integer arrays for to, cap and nxt are far more cache-friendly than per-node lists of objects. In C++ or Rust, use 32-bit indices and 64-bit capacities.
• Capacity types. Sums of capacities can overflow 32-bit integers long before any single capacity does. With floating-point capacities, treat residuals below a small epsilon as zero or the algorithm can loop on dust; better, scale to integers.
• Reading the min cut. After the run, BFS from s over arcs with positive residual capacity. Reached nodes are the source side; original arcs from that side to the rest form the minimum cut. Do not use ISAP's labels for this directly: after a gap break, the labels are stale above the gap.
• Reuse. For many queries on one graph, keep a copy of the original cap array and slice-assign it back before each run.

Failure modes

• No current-arc reset on relabel. The scan skips arcs that just became admissible, so the algorithm hangs or stops early with too little flow. It hides on small tests.
• Gap counts out of sync. Forgetting to decrement the old level or increment the new one makes the gap fire too early (wrong answer) or never (slow).
• Missing initial BFS. Starting every label at zero is still a valid labelling, so the algorithm stays correct, but it spends a long warm-up relabelling nodes one level at a time.
• Arc pairing broken. If arcs are not added in exact pairs, e ^ 1 points at an unrelated arc and the flow becomes nonsense without any crash.
• Source equals sink. Without the guard, the sink is reached at once with an infinite bottleneck and the loop never ends.

What to do next

1. Type the implementation from memory and test it against a brute-force Edmonds-Karp on a few thousand random small graphs, including disconnected ones.
2. Remove the gap check and the current-arc pointer one at a time and measure relabel and augmentation counts, so you know what each optimisation buys.
3. Solve a bipartite matching instance with ISAP and with Dinic and compare times.
4. Add min-cut extraction with a residual BFS and verify the cut capacity equals the flow.
5. Port it to your production language with 64-bit capacities and flat arrays.
6. Keep a benchmark of your real graph family so future algorithm swaps are measured, not guessed.