Part Ⅲ

16 Minimum Spanning Trees

16.1 Introduction to MSTs

Spanning tree: A spanning tree is a subgraph that is both a tree (connected and acyclic) and spanning (includes all of the vertices).

Application:

Dithering
Cluster analysis
Max bottleneck paths
Real-time face verification
LDPC codes for error correction
Image registration with Renyi entropy
Find road networks in satellite and aerial imagery
Reducing data storage in sequencing amino acids in a protein
Model locality of particle interactions in turbulent fluid flows
Autoconfig protocol for Ethernet bridging to avoid cycles in a network
Approximation algorithms for NP-hard problems (e.g., TSP, Steiner tree)
Network design (communication, electrical, hydraulic, computer, road)

16.2 Greedy Algorithm

Definitions

Cut: A cut in a graph is a partition of its vertices into two (nonempty) sets.
Crossing edge: A crossing edge connects a vertex in one set with a vertex in the other.

Greedy Algorithm for MST

Start with all edges colored gray.
Find cut with no black crossing edges; color its min-weight edge black.
Repeat until edges are colored black.

Correctness Proof

Given any cut, the crossing edge of min weight is in MST.
Proof
Suppose min-weight crossing edge is not in the MST.
- Adding to the MST creates a cycle.
- Some other edge in cycle must be a crossing edge.
- Removing and adding is also a spanning edge.
- Since weight of is less than the weight of , that spanning tree is lower height.
- Contradiction.
The greedy algorithm computes the MST.
Proof
- Any edge colored black is in the MST (via cut property).
- Fewer than black edges => cut with no black crossing edges (consider cut whose vertices are one connected component).

16.3 Edge-weighted Graph Implementation

Edge

public class Edge implements Comparable<Edge> { private final int source; private final int destination; private final double weight; public Edge(int source, int destination, double weight) { this.source = source; this.destination = destination; this.weight = weight; } public int getEitherVertex() { return source; } public int getOtherVertex(int vertex) { if (vertex == source) { return destination; } else if (vertex == destination) { return source; } else { throw new IllegalArgumentException("Invalid vertex"); } } public double getWeight() { return weight; } @Override public String toString() { return "(" + source + " - " + destination + " : " + weight + ")"; } @Override public int compareTo(Edge other) { return Double.compare(this.weight, other.weight); } }

#ifndef EDGE_H #define EDGE_H #include <string> class Edge { private: int source; int destination; double weight; public: Edge(int source, int destination, double weight); [[nodiscard]] int getEitherVertex() const; [[nodiscard]] int getOtherVertex(int vertex) const; [[nodiscard]] double getWeight() const; bool operator<(const Edge& other) const; [[nodiscard]] std::string toString() const; }; #endif // EDGE_H

#include "Edge.h" #include <stdexcept> #include <sstream> #include <iomanip> Edge::Edge(const int source, const int destination, const double weight) : source(source), destination(destination), weight(weight) {} int Edge::getEitherVertex() const { return source; } int Edge::getOtherVertex(const int vertex) const { if (vertex == source) { return destination; } else if (vertex == destination) { return source; } else { throw std::invalid_argument("Invalid vertex"); } } double Edge::getWeight() const { return weight; } bool Edge::operator<(const Edge& other) const { return weight < other.weight; } std::string Edge::toString() const { std::ostringstream oss; oss << "(" << source << " - " << destination << " : " << std::fixed << std::setprecision(2) << weight << ")"; return oss.str(); }

class Edge: def __init__(self, source: int, destination: int, weight: float): self.source = source self.destination = destination self.weight = weight def get_either_vertex(self) -> int: return self.source def get_other_vertex(self, vertex: int) -> int: if vertex == self.source: return self.destination elif vertex == self.destination: return self.source else: raise ValueError("Invalid vertex") def get_weight(self) -> float: return self.weight def __str__(self) -> str: return f"({self.source} - {self.destination} : {self.weight})" def __lt__(self, other: 'Edge') -> bool: return self.weight < other.weight

Edge-weighted Graph

import java.util.ArrayList; import java.util.List; public class EdgeWeightedGraph { private final int vertices; private final List<Edge>[] adjacencyList; public EdgeWeightedGraph(int vertices) { this.vertices = vertices; this.adjacencyList = new ArrayList[vertices]; for (int i = 0; i < vertices; i++) { adjacencyList[i] = new ArrayList<>(); } } public void addEdge(int source, int destination, double weight) { adjacencyList[source].add(new Edge(source, destination, weight)); adjacencyList[destination].add(new Edge(destination, source, weight)); } public int getVertices() { return vertices; } public List<Edge> getAdjacencyList(int vertex) { return adjacencyList[vertex]; } public void printGraph() { for (int i = 0; i < vertices; i++) { List<Edge> edges = adjacencyList[i]; System.out.print("Vertex " + i + ":"); for (Edge edge : edges) { System.out.print(" " + edge); } System.out.println(); } } }

#ifndef EDGEWEIGHTEDGRAPH_H #define EDGEWEIGHTEDGRAPH_H #include <vector7gt; #include "Edge.h" class EdgeWeightedGraph { private: int vertices; std::vector<std::vector<Edge>> adjacencyList; public: explicit EdgeWeightedGraph(int vertices); void addEdge(int source, int destination, double weight); [[nodiscard]] int getVertices() const; [[nodiscard]] const std::vector<Edge>& getAdjacencyList(int vertex) const; void printGraph() const; }; #endif // EDGEWEIGHTEDGRAPH_H

#include "EdgeWeightedGraph.h" #include <iostream> EdgeWeightedGraph::EdgeWeightedGraph(const int vertices) : vertices(vertices), adjacencyList(vertices) {} void EdgeWeightedGraph::addEdge(const int source, const int destination, const double weight) { const Edge edge(source, destination, weight); adjacencyList[source].push_back(edge); adjacencyList[destination].emplace_back(destination, source, weight); } int EdgeWeightedGraph::getVertices() const { return vertices; } const std::vector<Edge>& EdgeWeightedGraph::getAdjacencyList(int vertex) const { return adjacencyList[vertex]; } void EdgeWeightedGraph::printGraph() const { for (int i = 0; i < vertices; ++i) { const std::vector<Edge>& edges = adjacencyList[i]; std::cout << "Vertex " << i << ":"; for (const Edge& edge : edges) { std::cout << " " << edge.toString(); } std::cout << std::endl; } }

from Edge import Edge class EdgeWeightedGraph: def __init__(self, vertices: int): self.vertices = vertices self.adjacency_list: list[list[Edge]] = [[] for _ in range(vertices)] def add_edge(self, source: int, destination: int, weight: float): self.adjacency_list[source].append(Edge(source, destination, weight)) self.adjacency_list[destination].append(Edge(destination, source, weight)) def get_vertices(self) -> int: return self.vertices def get_adjacency_list(self, vertex: int) -> list[Edge]: return self.adjacency_list[vertex] def print_graph(self): for i in range(self.vertices): print(f"Vertex {i}:", end="") for edge in self.adjacency_list[i]: print(f" {edge}", end="") print()

16.4 Kruskal's Algorithm

Kruskal's Algorithm

Consider edges in ascending order of weight.
Add next edge to tree unless doing so would create a cycle.

Union-Find for Cycle Challenge

Maintain a set for each connected component in
If and are in same set, then adding would create a cycle.
To add to , merge sets containing and .

Correctness Proof

Kruskal's Algorithm is a special case of the greedy MST algorithm.

Suppose Kruskal's algorithm colors the edge black.
Cut = set of vertices connected to in tree .
No crossing edge is black.
No crossing edge has lower weight.

Property: Kruskal's algorithm computes MST in time proportional to (in the worst case).

Proof

Operation	Frequency	Time per op
Build pq
Delete-min
Build pq
Connected

import java.util.ArrayList; import java.util.List; import java.util.PriorityQueue; public class KruskalsAlgorithm { public static List<Edge> findMinimumSpanningTree(EdgeWeightedGraph graph) { int vertices = graph.getVertices(); List<Edge> minimumSpanningTree = new ArrayList<>(); PriorityQueue<Edge> minHeap = new PriorityQueue<>(graph.getVertices()); UnionFind unionFind = new UnionFind(vertices); for (int i = 0; i < vertices; i++) { for (Edge edge : graph.getAdjacencyList(i)) { if (edge.getEitherVertex() < i) { minHeap.offer(edge); } } } while (!minHeap.isEmpty() && minimumSpanningTree.size() < vertices - 1) { Edge edge = minHeap.poll(); int source = edge.getEitherVertex(); int destination = edge.getOtherVertex(source); int sourceRoot = unionFind.find(source); int destinationRoot = unionFind.find(destination); if (sourceRoot != destinationRoot) { minimumSpanningTree.add(edge); unionFind.union(sourceRoot, destinationRoot); } } return minimumSpanningTree; } }

#include "EdgeWeightedGraph.h" #include <queue> #include <vector> class UnionFind { public: explicit UnionFind(int size); int find(int element); void unionSets(int element1, int element2); private: std::vector<int> parent; std::vector<int> rank; }; std::vector<Edge> findMinimumSpanningTree(const EdgeWeightedGraph& graph) { const int vertices = graph.getVertices(); std::vector<Edge> minimumSpanningTree; std::priority_queue<Edge, std::vector<Edge>, std::greater<>> minHeap; UnionFind unionFind(vertices); for (int i = 0; i < vertices; ++i) { for (const Edge& edge : graph.getAdjacencyList(i)) { minHeap.push(edge); } } // Build the minimum spanning tree while (!minHeap.empty() && minimumSpanningTree.size() < vertices - 1) { Edge edge = minHeap.top(); minHeap.pop(); const int source = edge.getEitherVertex(); int destination = edge.getOtherVertex(source); // Check if connecting these vertices creates a cycle if (unionFind.find(source) != unionFind.find(destination)) { minimumSpanningTree.push_back(edge); unionFind.unionSets(source, destination); } } return minimumSpanningTree; } UnionFind::UnionFind(const int size) : parent(size), rank(size, 1) { for (int i = 0; i < size; ++i) { parent[i] = i; } } int UnionFind::find(const int element) { if (parent[element] != element) { parent[element] = find(parent[element]); } return parent[element]; } void UnionFind::unionSets(const int element1, const int element2) { const int root1 = find(element1); const int root2 = find(element2); if (root1 != root2) { if (rank[root1] > rank[root2]) { parent[root2] = root1; } else if (rank[root1] < rank[root2]) { parent[root1] = root2; } else { parent[root2] = root1; rank[root1]++; } } }

from typing import List import heapq from EdgeWeightedGraph import EdgeWeightedGraph, Edge class KruskalsAlgorithm: @staticmethod def find_minimum_spanning_tree(graph: EdgeWeightedGraph) -> List[Edge]: vertices: int = graph.get_vertices() minimum_spanning_tree: List[Edge] = [] min_heap: list[Edge] = [] union_find = UnionFind(vertices) # Add edges to the min-heap, ensuring no duplicates for i in range(vertices): for edge in graph.get_adjacency_list(i): # Add edge only if its source vertex is smaller than its destination if edge.source < edge.destination: heapq.heappush(min_heap, edge) while min_heap and len(minimum_spanning_tree) < vertices - 1: edge: Edge = heapq.heappop(min_heap) source: int = edge.get_either_vertex() destination: int = edge.get_other_vertex(source) source_root: int = union_find.find(source) destination_root: int = union_find.find(destination) if source_root != destination_root: minimum_spanning_tree.append(edge) union_find.union(source_root, destination_root) return minimum_spanning_tree class UnionFind: def __init__(self, size: int): self.parent: List[int] = [i for i in range(size)] self.rank: List[int] = [1] * size def find(self, element: int) -> int: if self.parent[element] != element: self.parent[element] = self.find(self.parent[element]) return self.parent[element] def union(self, element1: int, element2: int): root1: int = self.find(element1) root2: int = self.find(element2) if root1 != root2: if self.rank[root1] > self.rank[root2]: self.parent[root2] = root1 elif self.rank[root1] < self.rank[root2]: self.parent[root1] = root2 else: self.parent[root2] = root1 self.rank[root1] += 1

16.5 Prim's Algorithm

Prim's Algorithm

Start with vertex and greedily grow tree
Add to the min weight edge with exactly one endpoint in .
Repeat until edges.

Correctness Proof

Prim's Algorithm is a special case of the greedy MST algorithm.

Suppose edge e = min weight edge connecting a vertex on the tree to a vertex not on the tree.
Cut = set of vertices connected on tree.
No crossing edge is black.
No crossing edge has lower weight.

16.5.1 Lazy Implementation

Lazy Implementation

Maintain a PQ of edges with (at least) one endpoint in T.
Key = edge; priority = weight of edge.
Delete-min to determine next edge to add to .
Disregard if both endpoints and are in .
Otherwise, let be the vertex not in .
Add to PQ any edge incident to (assuming other endpoint not in T)
Add to and mark .

Property: Lazy Prim's algorithm computes the MST in time proportional to and extra space proportional to (in the worst case).

Operation	Frequency	Binary Heap
Delete min
Insert

import java.util.ArrayList; import java.util.List; import java.util.PriorityQueue; public class PrimMSTLazy { private final boolean[] marked; private final PriorityQueue<Edge> pq; private final List<Edge> mst; private double weight; public PrimMSTLazy(EdgeWeightedGraph graph) { marked = new boolean[graph.getVertices()]; pq = new PriorityQueue<>(); mst = new ArrayList<>(); weight = 0.0; visit(graph, 0); while (!pq.isEmpty()) { Edge e = pq.poll(); int v = e.getEitherVertex(); int w = e.getOtherVertex(v); if (marked[v] && marked[w]) continue; mst.add(e); weight += e.getWeight(); if (!marked[v]) visit(graph, v); if (!marked[w]) visit(graph, w); } } private void visit(EdgeWeightedGraph graph, int v) { marked[v] = true; for (Edge e : graph.getAdjacencyList(v)) { if (!marked[e.getOtherVertex(v)]) { pq.offer(e); } } } public Iterable<Edge> edges() { return mst; } public double weight() { return weight; } }

#include <iostream> #include <vector> #include <queue> #include "EdgeWeightedGraph.h" class PrimMSTLazy { private: std::vector<bool> marked; std::priority_queue<Edge, std::vector<Edge>, std::greater<>> pq; std::vector<Edge> mst; double weight; void visit(const EdgeWeightedGraph& graph, int v) { marked[v] = true; for (const Edge& e : graph.getAdjacencyList(v)) { if (!marked[e.getOtherVertex(v)]) { pq.push(e); } } } public: explicit PrimMSTLazy(const EdgeWeightedGraph& graph) : marked(graph.getVertices(), false), weight(0.0) { visit(graph, 0); while (!pq.empty()) { Edge e = pq.top(); pq.pop(); int v = e.getEitherVertex(); int w = e.getOtherVertex(v); if (marked[v] && marked[w]) continue; mst.push_back(e); weight += e.getWeight(); if (!marked[v]) visit(graph, v); if (!marked[w]) visit(graph, w); } } [[nodiscard]] const std::vector<Edge>& edges() const { return mst; } [[nodiscard]] double getWeight() const { return weight; } };

from typing import List, Iterable import heapq from EdgeWeightedGraph import EdgeWeightedGraph, Edge class PrimMSTLazy: def __init__(self, graph: EdgeWeightedGraph): self.marked: List[bool] = [False] * graph.get_vertices() self.pq: List[Edge] = [] # Min-heap for edges self.mst: List[Edge] = [] # Stores the MST edges self.weight: float = 0.0 self._visit(graph, 0) # Start from vertex 0 while self.pq: edge: Edge = heapq.heappop(self.pq) v: int = edge.get_either_vertex() w: int = edge.get_other_vertex(v) if self.marked[v] and self.marked[w]: continue # Ignore if both vertices are already in the MST self.mst.append(edge) self.weight += edge.get_weight() if not self.marked[v]: self._visit(graph, v) if not self.marked[w]: self._visit(graph, w) def _visit(self, graph: EdgeWeightedGraph, v: int): """Adds edges connected to vertex v to the priority queue.""" self.marked[v] = True for edge in graph.get_adjacency_list(v): if not self.marked[edge.get_other_vertex(v)]: heapq.heappush(self.pq, edge) def edges(self) -> Iterable[Edge]: return self.mst def weight(self) -> float: return self.weight

16.5.2 Eager Implementation

Property

Running time depends on PQ implementation: insert, delete-min, decrease-key.

PQ Implementation	Insert	Delete-Min	Decrease-Key	Total
Array
Binary Heap
d-way Heap (Johnson 1975)
Fibonacci Heap (Fredman-Tarjan 1984)

*: amortized

Bottom Line

Array implementation optimal for dense graph.
Binary heap much faster for sparse graphs.
4-way heap worth the trouble in performance-critical situations.
Fibonacci heap best in theory, but not worth implementing.

import java.util.ArrayList; import java.util.List; public class PrimMST { private final boolean[] marked; private final Edge[] edgeTo; private final double[] distTo; private final IndexedPriorityQueue pq; private final List<Edge> mst; public PrimMST(EdgeWeightedGraph graph) { marked = new boolean[graph.getVertices()]; edgeTo = new Edge[graph.getVertices()]; distTo = new double[graph.getVertices()]; pq = new IndexedPriorityQueue(graph.getVertices()); mst = new ArrayList<>(); for (int v = 0; v < graph.getVertices(); v++) { distTo[v] = Double.POSITIVE_INFINITY; } distTo[0] = 0.0; pq.insert(0, 0.0); while (!pq.isEmpty()) { visit(graph, pq.delMin()); } } private void visit(EdgeWeightedGraph graph, int vertex) { marked[vertex] = true; for (Edge edge : graph.getAdjacencyList(vertex)) { int w = edge.getOtherVertex(vertex); if (marked[w]) continue; if (edge.getWeight() < distTo[w]) { edgeTo[w] = edge; distTo[w] = edge.getWeight(); if (pq.contains(w)) { pq.decreaseKey(w, distTo[w]); } else { pq.insert(w, distTo[w]); } } } } public Iterable<Edge> edges() { for (int v = 1; v < edgeTo.length; v++) { if (edgeTo[v] != null) { mst.add(edgeTo[v]); } } return mst; } public double weight() { double weight = 0.0; for (Edge edge : mst) { weight += edge.getWeight(); } return weight; } }

#include "IndexedPriorityQueue.h" #include "EdgeWeightedGraph.h" #include <vector> #include <limits> class PrimMST { private: std::vector<bool> marked; std::vector<Edge> edgeTo; std::vector<double> distTo; IndexedPriorityQueue pq; std::vector<Edge> mst; void visit(const EdgeWeightedGraph& graph, int vertex) { marked[vertex] = true; for (const Edge& edge : graph.getAdjacencyList(vertex)) { const int w = edge.getOtherVertex(vertex); if (marked[w]) continue; if (edge.getWeight() < distTo[w]) { edgeTo[w] = edge; distTo[w] = edge.getWeight(); if (pq.contains(w)) { pq.decreaseKey(w, distTo[w]); } else { pq.insert(w, distTo[w]); } } } } public: explicit PrimMST(const EdgeWeightedGraph& graph) : marked(graph.getVertices(), false), edgeTo(graph.getVertices()), distTo(graph.getVertices(), std::numeric_limits<double>::infinity()), pq(graph.getVertices()) { distTo[0] = 0.0; pq.insert(0, 0.0); while (!pq.isEmpty()) { visit(graph, pq.delMin()); } } const std::vector<Edge>& edges() { mst.clear(); for (int v = 1; v < edgeTo.size(); v++) { if (edgeTo[v].getWeight() != 0.0) { mst.push_back(edgeTo[v]); } } return mst; } [[nodiscard]] double weight() const { double weight = 0.0; for (const Edge& edge : mst) { weight += edge.getWeight(); } return weight; } };

from typing import List, Iterable from EdgeWeightedGraph import EdgeWeightedGraph, Edge from IndexedPriorityQueue import IndexedPriorityQueue class PrimMSTEager: def __init__(self, graph: EdgeWeightedGraph): self.marked: List[bool] = [False] * graph.get_vertices() self.edge_to: List[Edge] = [None] * graph.get_vertices() self.dist_to: List[float] = [float('inf')] * graph.get_vertices() self.pq: IndexedPriorityQueue = IndexedPriorityQueue(graph.get_vertices()) self.mst: List[Edge] = [] self.dist_to[0] = 0.0 self.pq.insert(0, 0.0) while not self.pq.is_empty(): self._visit(graph, self.pq.del_min()) def _visit(self, graph: EdgeWeightedGraph, vertex: int): self.marked[vertex] = True for edge in graph.get_adjacency_list(vertex): w: int = edge.get_other_vertex(vertex) if self.marked[w]: continue if edge.get_weight() < self.dist_to[w]: self.dist_to[w] = edge.get_weight() self.edge_to[w] = edge if self.pq.contains(w): self.pq.decrease_key(w, self.dist_to[w]) else: self.pq.insert(w, self.dist_to[w]) def edges(self) -> Iterable[Edge]: """Returns an iterable of edges in the MST.""" for v in range(1, len(self.edge_to)): if self.edge_to[v] is not None: self.mst.append(self.edge_to[v]) return self.mst def weight(self) -> float: """Returns the total weight of the MST.""" return sum(edge.get_weight() for edge in self.mst)

16.6 MST Context

16.6.1 Euclidean MST

Euclidean MST: Given points in the plane, find MST connecting them, where the distances between point pairs are their Euclidean distances.

Methods: Exploit geometry and do it in

16.6.2 Single Link Clustering

Definitions

k-clustering: Divide a set of objects calssify into coherent groups.
Distance Function: Numeric value specifying "closeness" of two objects.
Single link: Distance between two clusters equals the distance between the two closest objects (one in each cluster).
Single-link clustering: Given an integer k, find a k-clustering that maximizes the distance between two closest clusters.

"Well-known" algorithm in science literature for single-link clustering:

Form clusters of one object each.
Find the closest pair of objects such that each object is in a different cluster, and merge the two clusters.
Repeat until there are exactly clusters.

Applications

Routing in mobile ad hoc networks.
Document categorization for web search.
Similarity searching in medical image databases.
Skycat: cluster sky objects into stars, quasars, galaxies.

16.7 Important Questions

Q: Bottleneck minimum spanning tree: Given a connected edge-weighted graph, design an efficient algorithm to find a minimum bottleneck spanning tree. The bottleneck capacity of a spanning tree is the weights of its largest edge. A minimum bottleneck spanning tree is a spanning tree of minimum bottleneck capacity.
A: Prove that an MST is a minimum bottleneck spanning tree.
Q: Is an edge in a MST: Given an edge-weighted graph and an edge , design a linear-time algorithm to determine whether appears in some MST of .
Note: Since your algorithm must take linear time in the worst case, you cannot afford to compute the MST itself.
A: Consider the subgraph of containing only those edges whose weight is strictly less than that of .
Q: Minimum-weight feedback edge set: A feedback edge set of a graph is a subset of edges that contains at least one edge from every cycle in the graph. If the edges of a feedback edge set are removed, the resulting graph is acyclic. Given an edge-weighted graph, design an efficient algorithm to find a feedback edge set of minimum weight. Assume the edge weights are positive.
A: Complement of an MST.

17 Shortest Paths

17.1 Shortest Paths APIs

Goal: Given an edge-weighted digraph, find the shortest path from to .

Navigation
PERT/CPM
Map routing
Seam carving
Robot navigation
Texture mapping
Typesetting in TeX
Urban traffic planning
Optimal pipelining of VLSI chip
Telemarketer operator scheduling
Routing of telecommunications messages
Network routing protocols (OSPF, BGP, RIP)
Exploiting arbitrage opportunities in currency exchange
Optimal truck routing through given traffic congestion pattern

Directed Edge

public class DirectedEdge { private final int v; private final int w; private final double weight; public DirectedEdge(int v, int w, double weight) { this.v = v; this.w = w; this.weight = weight; } public int from() { return v; } public int to() { return w; } public double weight() { return weight; } @Override public String toString() { return v + "->" + w + " " + String.format("%5.2f", weight); } }

#ifndef DIRECTEDEDGE_H #define DIRECTEDEDGE_H #include <iostream> class DirectedEdge { private: int v; int w; double weight; public: explicit DirectedEdge(int v = -1, int w = -1, double weight = 0.0); [[nodiscard]] int from() const; [[nodiscard]] int to() const; [[nodiscard]] double getWeight() const; friend std::ostream& operator<<(std::ostream& out, const DirectedEdge& e); }; #endif // DIRECTEDEDGE_H

#include "DirectedEdge.h" #include <iostream> DirectedEdge::DirectedEdge(const int v, const int w, const double weight) : v(v), w(w), weight(weight) {} int DirectedEdge::from() const { return v; } int DirectedEdge::to() const { return w; } double DirectedEdge::getWeight() const { return weight; } std::ostream& operator<<(std::ostream& out, const DirectedEdge& e) { out << e.v << "->" << e.w << " " << e.weight; return out; }

class DirectedEdge: def __init__(self, v, w, weight): self.v = v self.w = w self.weight = weight def from_vertex(self): return self.v def to_vertex(self): return self.w def get_weight(self): return self.weight def __str__(self): return f"{self.v}->{self.w} ({self.weight})"

Edge-Weighted Digraph

import java.util.ArrayList; import java.util.List; public class EdgeWeightedDigraph { private final int V; private final List<DirectedEdge>[] adj; public EdgeWeightedDigraph(int V) { this.V = V; adj = (List<DirectedEdge>[]) new ArrayList[V]; for (int v = 0; v < V; v++) adj[v] = new ArrayList<DirectedEdge>(); } public void addEdge(int source, int destination, double weight) { DirectedEdge e = new DirectedEdge(source, destination, weight); adj[source].add(e); } public Iterable<DirectedEdge> adj(int v) { return adj[v]; } public int V() { return V; } public int E() { int count = 0; for (int v = 0; v < V; v++) count += adj[v].size(); return count; } public String toString() { StringBuilder s = new StringBuilder(); s.append(V).append(" vertices, ").append(E()).append(" edges ").append("\n"); for (int v = 0; v < V; v++) { s.append(v).append(": "); for (DirectedEdge e : adj[v]) { s.append(e).append(" "); } s.append("\n"); } return s.toString(); } }

#ifndef EDGEWEIGHTEDDIGRAPH_H #define EDGEWEIGHTEDDIGRAPH_H #include "DirectedEdge.h" #include <vector> #include <iostream> class EdgeWeightedDigraph { private: int V; std::vector<std::vector<DirectedEdge>> adj; public: explicit EdgeWeightedDigraph(int V); void addEdge(int source, int destination, double weight); [[nodiscard]] std::vector<DirectedEdge> getAdj(int v) const; [[nodiscard]] int getV() const; [[nodiscard]] int getE() const; friend std::ostream& operator<<(std::ostream& out, const EdgeWeightedDigraph& G); }; #endif // EDGEWEIGHTEDDIGRAPH_H

#include "EdgeWeightedDigraph.h" EdgeWeightedDigraph::EdgeWeightedDigraph(const int V) : V(V), adj(V) {} void EdgeWeightedDigraph::addEdge(const int source, const int destination, const double weight) { const DirectedEdge e(source, destination, weight); adj[source].push_back(e); } std::vector<DirectedEdge> EdgeWeightedDigraph::getAdj(const int v) const { return adj[v]; } int EdgeWeightedDigraph::getV() const { return V; } int EdgeWeightedDigraph::getE() const { std::size_t count = 0; for (int v = 0; v < V; ++v) { count += adj[v].size(); } return static_cast<int>(count); } std::ostream& operator<<(std::ostream& out, const EdgeWeightedDigraph& G) { out << G.V << " vertices, " << G.getE() << " edges\n"; for (int v = 0; v < G.V; ++v) { out << v << ": "; for (const auto& e : G.adj[v]) { out << e << " "; } out << "\n"; } return out; }

from DirecteEdge import DirectedEdge class EdgeWeightedDigraph: def __init__(self, V): self.V = V self.adj = [[] for _ in range(V)] def add_edge(self, source, destination, weight): e = DirectedEdge(source, destination, weight) self.adj[source].append(e) def get_adj(self, v): return self.adj[v] def get_V(self): return self.V def get_E(self): count = 0 for v in range(self.V): count += len(self.adj[v]) return count def __str__(self): s = f"{self.V} vertices, {self.get_E()} edges\n" for v in range(self.V): s += f"{v}: " for e in self.adj[v]: s += f"{e} " s += "\n" return s

17.2 Shortest Path Properties

Cost of undirected shortest paths:

Proof: Undirected shortest paths (with nonnegative weights) reduces to directed shortest path.

cost of shortest paths in digraph
cost of reduction

Edge Relaxation

distTo[v] is length of shortest known path from to .
distTo[w] is length of shortest known path from to .
edgeTo[w] is last edge on shortest known path from to .
If gives shorter path to through , update both distTo[w] and edgeTo[w].

Correctness Proof: Shortest-paths optimality conditions

Let be an edge-weighted digraph.

Then distTo[] are the shortest path distances from iff:

distTo[s] = 0.
For each vertex v, distTo[v] is the length of some path from to .
For each edge , distTo[w] <= distTo[v] + e.weight().

Proof

Suppose that is a shortest path from to .
Then,
= edge on shortest path from to .
Add inequalities; simplify; and substitute
is the weight of shortest path from to .
Thus, distTo[w] is the weight of shortest path to .

Different Implementations

Dijkstra's algorithm (nonnegative weights).
Topological sort algorithm (no directed cycles).
Bellman-Ford algorithm (no negative cycles).

17.3 Dijkstra's Algorithm

Dijkstra's Algorithm

Consider vertices in increasing order of distance from s.
(non-tree vertex with the lowest distTo[] value)
Add vertex to tree and relaw all edges pointing from that vertex.

Correctness Proof: Dijkstra's algorithm computes a SPT in any edge-weighted digraph with nonnegative weights.

Each edge is relaxed exactly once (when v is relaxed), leaving .
Inequality holds until algorithm terminates because:
- distTo[w] cannot increase => distTo[] values are monotone decreasing.
- distTo[v] will not change => we choose lowest distTo[] value at each step and edge weights are nonnegative)
Thus, upon termination, shortest-paths optimality conditions hold.

Prim’s algorithm is essentially the same algorithm as Dijkstra's algorithm

Both are in a family of algorithms that compute a graph's spanning tree.
Prim's: Closest vertex to the tree (via an undirected edge).
Dijkstra's: Closest vertex to the source (via a directed path).

import java.util.ArrayList; import java.util.Comparator; import java.util.List; import java.util.PriorityQueue; public class Dijkstra { private final double[] distTo; private final DirectedEdge[] edgeTo; private final boolean[] marked; private final PriorityQueue<Integer> pq; public Dijkstra(EdgeWeightedDigraph G, int s) { distTo = new double[G.V()]; edgeTo = new DirectedEdge[G.V()]; marked = new boolean[G.V()]; for (int v = 0; v < G.V(); v++) distTo[v] = Double.POSITIVE_INFINITY; distTo[s] = 0.0; pq = new PriorityQueue<>(Comparator.comparingDouble(v -> distTo[v])); pq.offer(s); while (!pq.isEmpty()) { int v = pq.poll(); marked[v] = true; for (DirectedEdge e : G.adj(v)) { relax(e); } } } private void relax(DirectedEdge e) { int v = e.from(); int w = e.to(); if (distTo[w] > distTo[v] + e.weight()) { distTo[w] = distTo[v] + e.weight(); edgeTo[w] = e; if (!marked[w]) { pq.offer(w); } } } public double distTo(int v) { return distTo[v]; } public boolean hasPathTo(int v) { return distTo[v] < Double.POSITIVE_INFINITY; } public Iterable<DirectedEdge&gt pathTo(int v) { if (!hasPathTo(v)) return null; List<DirectedEdge> path = new ArrayList<>(); for (DirectedEdge e = edgeTo[v]; e != null; e = edgeTo[e.from()]) { path.add(e); } return path; } }

#ifndef DIJKSTRA_H #define DIJKSTRA_H #include "EdgeWeightedDigraph.h" #include <vector> #include <queue> class Dijkstra { private: std::vector<double> distTo; std::vector<DirectedEdge> edgeTo; std::vector<bool> marked; std::priority_queue<std::pair<double, int>, std::vector<std::pair<double, int>>, std::greater<>> pq; void relax(const DirectedEdge& e); public: explicit Dijkstra(const EdgeWeightedDigraph& G, int s); [[nodiscard]] double getdistTo(int v) const; [[nodiscard]] bool hasPathTo(int v) const; [[nodiscard]] std::vector<DirectedEdge> pathTo(int v) const; }; #endif // DIJKSTRA_H

#include "dijkstra.h" #include <limits> Dijkstra::Dijkstra(const EdgeWeightedDigraph& G, int s) : distTo(G.getV(), std::numeric_limits<double>::infinity()), edgeTo(G.getV(), DirectedEdge()), marked(G.getV(), false) { distTo[s] = 0.0; pq.emplace(0.0, s); while (!pq.empty()) { int v = pq.top().second; pq.pop(); if (marked[v]) continue; marked[v] = true; for (const auto& e : G.getAdj(v)) { relax(e); } } } void Dijkstra::relax(const DirectedEdge& e) { int v = e.from(); int w = e.to(); if (distTo[w] > distTo[v] + e.getWeight()) { distTo[w] = distTo[v] + e.getWeight(); edgeTo[w] = e; pq.emplace(distTo[w], w); } } double Dijkstra::getdistTo(int v) const { return distTo[v]; } bool Dijkstra::hasPathTo(int v) const { return distTo[v] < std::numeric_limits<double>::infinity(); } std::vector<DirectedEdge> Dijkstra::pathTo(int v) const { if (!hasPathTo(v)) return {}; std::vector<DirectedEdge> path; for (DirectedEdge e = edgeTo[v]; e.from() != -1; e = edgeTo[e.from()]) { path.push_back(e); } return path; }

from EdgeWeightedDigraph import EdgeWeightedDigraph import heapq class Dijkstra: def __init__(self, G, s): self.distTo = [float('inf')] * G.get_V() self.edgeTo = [None] * G.get_V() self.marked = [False] * G.get_V() self.pq = [] self.distTo[s] = 0.0 heapq.heappush(self.pq, (0.0, s)) while self.pq: _, v = heapq.heappop(self.pq) if self.marked[v]: continue self.marked[v] = True for e in G.get_adj(v): self.relax(e) def relax(self, e): v = e.from_vertex() w = e.to_vertex() if self.distTo[w] > self.distTo[v] + e.get_weight(): self.distTo[w] = self.distTo[v] + e.get_weight() self.edgeTo[w] = e heapq.heappush(self.pq, (self.distTo[w], w)) def dist_to(self, v): return self.distTo[v] def has_path_to(self, v): return self.distTo[v] < float('inf') def path_to(self, v): if not self.has_path_to(v): return None path = [] for e in reversed(self.edgeTo[v:v + 1]): if e is not None: path.append(e) return path

Property

Running time depends on PQ implementation: insert, delete-min, decrease-key.

PQ Implementation	Insert	Delete-Min	Decrease-Key	Total
Array
Binary Heap
d-way Heap (Johnson 1975)
Fibonacci Heap (Fredman-Tarjan 1984)

*: amortized

Bottom Line

Array implementation optimal for dense graph.
Binary heap much faster for sparse graphs.
4-way heap worth the trouble in performance-critical situations.
Fibonacci heap best in theory, but not worth implementing.

17.4 Edge-Weighted DAGs

Topological Sort Algorithm for Shortest Path

Consider all vertices in topological order.
Relax all edges pointing from that vertex.

Property: Topological sort algorithm computes SPT in any edgeweighted DAG in time proportional to .

Each edge is relaxed exactly once (when v is relaxed), leaving .
Inequality holds until algorithm terminates because:
- distTo[w] cannot increase => distTo[] values are monotone decreasing.
- distTo[v] will not change => we choose lowest distTo[] value at each step and edge weights are nonnegative)
Thus, upon termination, shortest-paths optimality conditions hold.

Longest paths in edge-weighted DAGs:

Formulate as a shortest paths problem in edge-weighted DAGs.

Negate all weights.
Find shortest paths.
Negate weights in result.

import java.util.*; public class ShortestPathTopological { private final EdgeWeightedDigraph graph; private final int source; private final double[] distTo; private final DirectedEdge[] edgeTo; public ShortestPathTopological(EdgeWeightedDigraph graph, int source) { this.graph = graph; this.source = source; distTo = new double[graph.V()]; edgeTo = new DirectedEdge[graph.V()]; for (int v = 0; v < graph.V(); v++) { distTo[v] = Double.POSITIVE_INFINITY; } distTo[source] = 0.0; TopologicalSort topologicalSort = new TopologicalSort(); List<Integer> sorted = topologicalSort.sort(graph); for (int v : sorted) { relax(v); } } private void relax(int v) { for (DirectedEdge edge : graph.adj(v)) { int w = edge.to(); if (distTo[w] > distTo[v] + edge.weight()) { distTo[w] = distTo[v] + edge.weight(); edgeTo[w] = edge; } } } public double distTo(int v) { return distTo[v]; } public boolean hasPathTo(int v) { return distTo[v] < Double.POSITIVE_INFINITY; } public Iterable<DirectedEdge> pathTo(int v) { if (!hasPathTo(v)) return null; List<DirectedEdge> path = new ArrayList<>(); for (DirectedEdge e = edgeTo[v]; e != null; e = edgeTo[e.from()]) { path.addFirst(e); } return path; } private static class TopologicalSort { public List<Integer> sort(EdgeWeightedDigraph graph) { int V = graph.V(); List<Integer> sorted = new ArrayList<>(); int[] inDegree = new int[V]; Queue<Integer> queue = new LinkedList<>(); for (int v = 0; v < V; v++) { for (DirectedEdge edge : graph.adj(v)) { inDegree[edge.to()]++; } } for (int v = 0; v < V; v++) { if (inDegree[v] == 0) { queue.offer(v); } } while (!queue.isEmpty()) { int u = queue.poll(); sorted.add(u); for (DirectedEdge edge : graph.adj(u)) { int v = edge.to(); inDegree[v]--; if (inDegree[v] == 0) { queue.offer(v); } } } if (sorted.size() != V) { throw new IllegalArgumentException("Graph contains a cycle."); } return sorted; } } }

#ifndef SHORTESTPATHTOPOLOGICAL_H #define SHORTESTPATHTOPOLOGICAL_H #include "EdgeWeightedDigraph.h" #include <vector> class ShortestPathTopological { private: const EdgeWeightedDigraph& graph; int source; std::vector<double> distTo; std::vector<DirectedEdge> edgeTo; class TopologicalSort { public: explicit TopologicalSort(const EdgeWeightedDigraph& graph); [[nodiscard]] std::vector<int> sort() const; private: const EdgeWeightedDigraph& graph; }; void relax(int v); public: explicit ShortestPathTopological(const EdgeWeightedDigraph& graph, int source); [[nodiscard]] double getdistTo(int v) const; [[nodiscard]] bool hasPathTo(int v) const; [[nodiscard]] std::vector<DirectedEdge> pathTo(int v) const; }; #endif // SHORTESTPATHTOPOLOGICAL_H

#include "ShortestPathTopological.h" #include <algorithm> #include <iostream> #include <queue> #include <limits> ShortestPathTopological::ShortestPathTopological(const EdgeWeightedDigraph &graph, int source) : graph(graph), source(source), distTo(graph.getV(), std::numeric_limits<double>::infinity()), edgeTo(graph.getV()) { distTo[source] = 0.0; TopologicalSort topologicalSort(graph); std::vector<int> sorted = topologicalSort.sort(); for (int v : sorted) { relax(v); } } void ShortestPathTopological::relax(const int v) { for (const DirectedEdge& edge : graph.getAdj(v)) { int w = edge.to(); if (distTo[w] > distTo[v] + edge.getWeight()) { distTo[w] = distTo[v] + edge.getWeight(); edgeTo[w] = edge; } } } double ShortestPathTopological::getdistTo(const int v) const { return distTo[v]; } bool ShortestPathTopological::hasPathTo(const int v) const { return distTo[v] < std::numeric_limits<double>::infinity(); } std::vector<DirectedEdge> ShortestPathTopological::pathTo(const int v) const { std::vector<DirectedEdge> path; if (!hasPathTo(v)) { return path; } for (DirectedEdge e = edgeTo[v]; e.from() != -1; e = edgeTo[e.from()]) { path.push_back(e); } std::ranges::reverse(path); return path; } ShortestPathTopological::TopologicalSort::TopologicalSort(const EdgeWeightedDigraph &graph) : graph(graph) {} std::vector<int> ShortestPathTopological::TopologicalSort::sort() const { const int V = graph.getV(); std::vector<int> sorted; std::vector<int> inDegree(V, 0); std::queue<int> queue; for (int v = 0; v < V; ++v) { for (const DirectedEdge& e : graph.getAdj(v)) { inDegree[e.to()]++; } } for (int v = 0; v < V; ++v) { if (inDegree[v] == 0) { queue.push(v); } } while (!queue.empty()) { int u = queue.front(); queue.pop(); sorted.push_back(u); for (const DirectedEdge& e : graph.getAdj(u)) { int w = e.to(); if (--inDegree[w] == 0) { queue.push(w); } } } if (sorted.size() != static_cast<size_t>(V)) { throw std::runtime_error("Graph contains a cycle!"); } return sorted; }

from collections import deque from typing import List, Optional from DirectedEdge import DirectedEdge from EdgeWeightedDigraph import EdgeWeightedDigraph class ShortestPathTopological: def __init__(self, graph: EdgeWeightedDigraph, source: int): self.graph = graph self.source = source self.dist_to = [float('inf')] * graph.get_V() self.edge_to: List[Optional[DirectedEdge]] = [None] * graph.get_V() self.dist_to[source] = 0.0 topological_order = self._topological_sort() for v in topological_order: self._relax(v) def _relax(self, v: int): for edge in self.graph.get_adj(v): w = edge.to_vertex() if self.dist_to[w] > self.dist_to[v] + edge.get_weight(): self.dist_to[w] = self.dist_to[v] + edge.get_weight() self.edge_to[w] = edge def get_dist_to(self, v: int) -> float: return self.dist_to[v] def has_path_to(self, v: int) -> bool: return self.dist_to[v] < float('inf') def path_to(self, v: int) -> Optional[List[str]]: if not self.has_path_to(v): return None path: List[DirectedEdge] = [] e = self.edge_to[v] while e is not None: path.append(e) e = self.edge_to[e.from_vertex()] return [str(edge) for edge in path[::-1]] def _topological_sort(self) -> List[int]: V = self.graph.get_V() in_degree = [0] * V for v in range(V): for edge in self.graph.get_adj(v): in_degree[edge.to_vertex()] += 1 queue = deque([v for v in range(V) if in_degree[v] == 0]) sorted_order = [] while queue: u = queue.popleft() sorted_order.append(u) for edge in self.graph.get_adj(u): v = edge.to_vertex() in_degree[v] -= 1 if in_degree[v] == 0: queue.append(v) if len(sorted_order) != V: raise ValueError("Graph contains a cycle.") return sorted_order

Application Ⅰ - Content-Aware Resizing

Seam Carving: Resize an image without distortion for display on cell phones and web browsers.

Grid DAG: vertex = pixel; edge = from pixel to 3 downward neighbors.
Weight of pixel = energy function of 8 neighboring pixels.
Seam = shortest path (sum of vertex weights) from top to bottom.

Application Ⅱ - Parallel Job Scheduling

Parallel Job Scheduling: Given a set of jobs with durations and precedence constraints, schedule the jobs (by finding a start time for each) so as to achieve the minimum completion time, while respecting the constraints.

To solve a parallel job-scheduling problem, create edge-weighted DAG, use longest path from the source to schedule each job:

Source and sink vertices.
Two vertices (begin and end) for each job.
Three edges for each job.
- begin to end (weighted by duration)
- source to begin (0 weight)
- end to sink (0 weight)
One edge for each precedence constraint (0 weight).

17.5 Negative Weights

Negative Cycle: A negative cycle is a directed cycle whose sum of edge weights is negative.

Bellman-Ford Algorithm

Initialize distTo[s] = 0 and distTo[v] = ∞ for all other vertices.
Repeat V times, relax each edge.

Practical Improvement: If distTo[v] does not change during pass , no need to relax any edge pointing from v in pass => maintain queue of vertices whose distTo[] changed.

Algorithm	Restriction	Typical Case	Worst Case	Extra Space
Topological Sort	No Directed Cycles
Dijkstra (Binary Heap)	No Negative Weights
Bellman-Ford	No Negative Cycles
Bellman-Ford (queue-based)	No Negative Cycles

import java.util.ArrayList; import java.util.LinkedList; import java.util.List; import java.util.Queue; public class BellmanFordSP { private final double[] distTo; private final DirectedEdge[] edgeTo; private final boolean[] onQueue; private final int[] cost; private final int s; private boolean hasNegativeCycle; private final Queue<Integer> q; public BellmanFordSP(EdgeWeightedDigraph G, int s) { this.s = s; distTo = new double[G.V()]; edgeTo = new DirectedEdge[G.V()]; onQueue = new boolean[G.V()]; cost = new int[G.V()]; for (int v = 0; v < G.V(); v++) distTo[v] = Double.POSITIVE_INFINITY; distTo[s] = 0.0; q = new LinkedList<>(); q.add(s); onQueue[s] = true; while (!q.isEmpty()) { int v = q.remove(); onQueue[v] = false; relax(G, v); } } private void relax(EdgeWeightedDigraph G, int v) { for (DirectedEdge e : G.adj(v)) { int w = e.to(); if (distTo[w] > distTo[v] + e.weight()) { distTo[w] = distTo[v] + e.weight(); edgeTo[w] = e; cost[w]++; if (!onQueue[w]) { q.add(w); onQueue[w] = true; } if (cost[w] >= G.V()) { hasNegativeCycle = true; return; } } } } public double distTo(int v) { return distTo[v]; } public boolean hasPathTo(int v) { return distTo[v] < Double.POSITIVE_INFINITY; } public Iterable<DirectedEdge> pathTo(int v) { if (!hasPathTo(v)) return null; List<DirectedEdge> path = new ArrayList<>(); for (DirectedEdge e = edgeTo[v]; e != null; e = edgeTo[e.from()]) { path.add(e); } return path; } public boolean hasNegativeCycle() { return hasNegativeCycle; } }

#ifndef BELLMANFORDSP_H #define BELLMANFORDSP_H #include "EdgeWeightedDigraph.h" #include <vector> #include <queue> class BellmanFordSP { private: std::vector<double> distTo; std::vector<DirectedEdge> edgeTo; std::vector<bool> onQueue; std::vector<int> cost; int s; bool hasNegativeCycle; std::queue<int> q; public: BellmanFordSP(const EdgeWeightedDigraph& G, int s); [[nodiscard]] double getdistTo(int v) const; [[nodiscard]] bool hasPathTo(int v) const; [[nodiscard]] std::vector<DirectedEdge> pathTo(int v) const; [[nodiscard]] bool NegativeCycle() const; private: void relax(const EdgeWeightedDigraph& G, int v); }; #endif // BELLMANFORDSP_H

#include "BellmanFordSP.h" #include <limits> BellmanFordSP::BellmanFordSP(const EdgeWeightedDigraph& G, const int s) : distTo(G.getV(), std::numeric_limits<double>::infinity()), edgeTo(G.getV(), DirectedEdge(-1, -1, 0.0)), onQueue(G.getV(), false), cost(G.getV(), 0), s(s), hasNegativeCycle(false) { distTo[s] = 0.0; q.push(s); onQueue[s] = true; while (!q.empty()) { int v = q.front(); q.pop(); onQueue[v] = false; relax(G, v); } } double BellmanFordSP::getdistTo(const int v) const { return distTo[v]; } bool BellmanFordSP::hasPathTo(const int v) const { return distTo[v] < std::numeric_limits<double>::infinity(); } std::vector<DirectedEdge> BellmanFordSP::pathTo(const int v) const { if (!hasPathTo(v)) { return {}; } std::vector<DirectedEdge> path; for (DirectedEdge e = edgeTo[v]; e.from() != -1; e = edgeTo[e.from()]) { path.push_back(e); } return path; } bool BellmanFordSP::NegativeCycle() const { return hasNegativeCycle; } void BellmanFordSP::relax(const EdgeWeightedDigraph& G, const int v) { for (const DirectedEdge& e : G.getAdj(v)) { int w = e.to(); if (distTo[w] > distTo[v] + e.getWeight()) { distTo[w] = distTo[v] + e.getWeight(); edgeTo[w] = e; cost[w]++; if (!onQueue[w]) { q.push(w); onQueue[w] = true; } if (cost[w] >= G.getV()) { hasNegativeCycle = true; return; } } } }

from EdgeWeightedDigraph import EdgeWeightedDigraph class BellmanFordSP: def __init__(self, G, s): self.distTo = [float("inf") for _ in range(G.get_V())] self.edgeTo = [None for _ in range(G.get_V())] self.onQueue = [False for _ in range(G.get_V())] self.cost = [0 for _ in range(G.get_V())] self.s = s self.hasNegativeCycle = False self.distTo[s] = 0.0 self.q = [s] self.onQueue[s] = True while self.q: v = self.q.pop(0) self.onQueue[v] = False self.relax(G, v) def distTo(self, v): return self.distTo[v] def hasPathTo(self, v): return self.distTo[v] != float("inf") def pathTo(self, v): if not self.hasPathTo(v): return None path = [] e = self.edgeTo[v] while e is not None: path.append(e) e = self.edgeTo[e.from_vertex()] return path def hasNegativeCycle(self): return self.hasNegativeCycle def relax(self, G, v): for e in G.get_adj(v): w = e.to_vertex() if self.distTo[w] > self.distTo[v] + e.get_weight(): self.distTo[w] = self.distTo[v] + e.get_weight() self.edgeTo[w] = e self.cost[w] += 1 if not self.onQueue[w]: self.q.append(w) self.onQueue[w] = True if self.cost[w] >= G.get_V(): self.hasNegativeCycle = True return

Find A Negative Cycle

If there is a negative cycle, Bellman-Ford gets stuck in loop, updating distTo[] and edgeTo[] entries of vertices in the cycle.

If any vertex v is updated in phase V, there exists a negative cycle (and can trace back edgeTo[v] entries to find it).

Application - Arbitrage Detection

Currency exchange graph.

Vertex = currency.
Edge = transaction, with weight equal to exchange rate.
Find a directed cycle whose product of edge weights is > 1.

Arbitrage Detection

Let weight of edge be (exchange rate from currency to ).
Multiplication turns to addition; turns to
Find a directed cycle whose sum of edge weights is (negative cycle).

18 Maximum Flow and Minimum Cut

18.1 Introduction

Definitions

-cut: A -cut (cut) is a partition of the vertices into two disjoint sets, with in one set and in the other set .

-cut capacity: Its capacity is the sum of the capacities of the edges from to .

Minimum cut problem: Find a cut of minimum capacity.

Definitions

-flow: An -flow (flow) is an assignment of values to the edges such that:

Capacity constraint: 0 ≤ edge's flow ≤ edge's capacity.
Local equilibrium: inflow = outflow at every vertex (except and ).

Value of a flow: The value of a flow is the inflow at (assuming no edge points to or from ).

Maximum st-flow (maxflow) problem: Find a flow of maximum value.

18.2 Ford-Fulkerson Algorithm

Ford-Fulkerson Algorithm

Start with 0 flow.
Find an undirected path from s to t such that:
1. Can increase flow on forward edges (not full).
2. Can decrease flow on backward edges (not empty).
Terminates when all paths from s to t are blocked by either a full forward edge or an empty backward edge.

18.3 Maxflow-Mincut Theorem

Definition

Net Flow: The net flow across a cut (, ) is the sum of the flows on its edges from to minus the sum of the flows on its edges from from to ).

Flow-value lemma: Let be any flow and let (, ) be any cut. Then, the net flow across (, ) equals the value of .

Proof: By induction on the size of .

Base case:
Induction step: remains true by local equilibrium when moving any vertex from to .

Weak duality: Let be any flow and let be any cut. Then, the value of the flow ≤ the capacity of the cut.

Proof

Value of flow = net flow across cut ≤ capacity of cut .

Augmenting path theorem: A flow f is a maxflow iff no augmenting paths.

Maxflow-mincut theorem: Value of the maxflow = capacity of mincut.

Proof: The following three conditions are equivalent for any flow .

There exists a cut whose capacity equals the value of the flow .
is a maxflow.
There is no augmenting path with respect to .

1 -> 2

Suppose that is a cut with capacity equal to the value of .
Then, the value of any flow ≤ capacity of = value of .
Thus, is a maxflow.

2 -> 3: We prove contrapositive: ~3 -> ~2

Suppose that there is an augmenting path with respect to .
Can improve flow by sending flow along this path.
Thus, f is not a maxflow.

3 -> 1: Suppose that there is no augmenting path with respect to .

Let be a cut where is the set of vertices connected to by an undirected path with no full forward or empty backward edges.
By definition, is in ; since no augmenting path, is in .
Capacity of cut = net flow across cut (forward edges full; backward edges empty) = value of flow (flow-value lemma).

To compute mincut from maxflow :

By augmenting path theorem, no augmenting paths with respect to .
Compute = set of vertices connected to by an undirected path with \ no full forward or empty backward edges.

18.4 Running Time Analysis

Properties

The flow is integer-valued throughout Ford-Fulkerson.
Proof
- Bottleneck capacity is an integer.
- Flow on an edge increases/decreases by bottleneck capacity.
Number of augmentations ≤ the value of the maxflow.
Proof: Each augmentation increases the value by at least 1.
Integrality theorem: There exists an integer-valued maxflow.
Proof: Ford-Fulkerson terminates and maxflow that it finds is integer-valued.

Running time: FF performance depends on choice of augmenting paths.

Digraph with vertices, edges, and integer capacities between 1 and

Augmenting Path	Number of Paths	Implementation
Shortest Path		Queue (BFS)
Fattest path		Priority Queue
Random Path		Randomized Queue
DFS Path		Stack (DFS)

18.5 Implementation

18.5.1 Flow Edge

Implementation

Use residual capcity:

Forward edge: residual capacity .
Backward edge: residual capacity .

public class FlowEdge { private static final double FLOATING_POINT_EPSILON = 1.0E-10; private final int v; private final int w; private final double capacity; private double flow; public FlowEdge(int v, int w, double capacity) { if (v < 0) throw new IllegalArgumentException("vertex index must be a non-negative integer"); if (w < 0) throw new IllegalArgumentException("vertex index must be a non-negative integer"); if (!(capacity >= 0.0)) throw new IllegalArgumentException("Edge capacity must be non-negative"); this.v = v; this.w = w; this.capacity = capacity; this.flow = 0.0; } public FlowEdge(int v, int w, double capacity, double flow) { if (v < 0) throw new IllegalArgumentException("vertex index must be a non-negative integer"); if (w < 0) throw new IllegalArgumentException("vertex index must be a non-negative integer"); if (!(capacity >= 0.0)) throw new IllegalArgumentException("edge capacity must be non-negative"); if (!(flow <= capacity)) throw new IllegalArgumentException("flow exceeds capacity"); if (!(flow >= 0.0)) throw new IllegalArgumentException("flow must be non-negative"); this.v = v; this.w = w; this.capacity = capacity; this.flow = flow; } public FlowEdge(FlowEdge e) { this.v = e.v; this.w = e.w; this.capacity = e.capacity; this.flow = e.flow; } public int from() { return v; } public int to() { return w; } public double capacity() { return capacity; } public double flow() { return flow; } public int other(int vertex) { if (vertex == v) return w; else if (vertex == w) return v; else throw new IllegalArgumentException("invalid endpoint"); } public double residualCapacityTo(int vertex) { if (vertex == v) return flow; else if (vertex == w) return capacity - flow; else throw new IllegalArgumentException("invalid endpoint"); } public void addResidualFlowTo(int vertex, double delta) { if (!(delta >= 0.0)) throw new IllegalArgumentException("Delta must be non-negative"); if (vertex == v) flow -= delta; else if (vertex == w) flow += delta; else throw new IllegalArgumentException("invalid endpoint"); if (Math.abs(flow) <= FLOATING_POINT_EPSILON) flow = 0; if (Math.abs(flow - capacity) <= FLOATING_POINT_EPSILON) flow = capacity; if (!(flow >= 0.0)) throw new IllegalArgumentException("Flow is negative"); if (!(flow <= capacity)) throw new IllegalArgumentException("Flow exceeds capacity"); } public String toString() { return v + "->" + w + " " + flow + "/" + capacity; } }

#ifndef FLOWEDGE_H #define FLOWEDGE_H #include <iostream> class FlowEdge { private: static constexpr double FLOATING_POINT_EPSILON = 1.0E-10; int v; int w; double capacity; double flow; public: FlowEdge(int v, int w, double capacity); FlowEdge(int v, int w, double capacity, double flow); FlowEdge(const FlowEdge& e); [[nodiscard]] int from() const; [[nodiscard]] int to() const; [[nodiscard]] double getcapacity() const; [[nodiscard]] double getflow() const; [[nodiscard]] int other(int vertex) const; [[nodiscard]] double residualCapacityTo(int vertex) const; void addResidualFlowTo(int vertex, double delta); friend std::ostream& operator<<(std::ostream& os, const FlowEdge& e); }; #endif // FLOWEDGE_H

#include "FlowEdge.h" #include <stdexcept> #include <cmath> FlowEdge::FlowEdge(const int v, const int w, const double capacity) : v(v), w(w), capacity(capacity), flow(0.0) { if (v < 0) throw std::invalid_argument("vertex index must be a non-negative integer"); if (w < 0) throw std::invalid_argument("vertex index must be a non-negative integer"); if (!(capacity >= 0.0)) throw std::invalid_argument("Edge capacity must be non-negative"); } FlowEdge::FlowEdge(const int v, const int w, const double capacity, const double flow) : v(v), w(w), capacity(capacity), flow(flow) { if (v < 0) throw std::invalid_argument("vertex index must be a non-negative integer"); if (w < 0) throw std::invalid_argument("vertex index must be a non-negative integer"); if (!(capacity >= 0.0)) throw std::invalid_argument("edge capacity must be non-negative"); if (!(flow <= capacity)) throw std::invalid_argument("flow exceeds capacity"); if (!(flow >= 0.0)) throw std::invalid_argument("flow must be non-negative"); } FlowEdge::FlowEdge(const FlowEdge& e) = default; int FlowEdge::from() const { return v; } int FlowEdge::to() const { return w; } double FlowEdge::getcapacity() const { return capacity; } double FlowEdge::getflow() const { return flow; } int FlowEdge::other(const int vertex) const { if (vertex == v) return w; else if (vertex == w) return v; else throw std::invalid_argument("invalid endpoint"); } double FlowEdge::residualCapacityTo(const int vertex) const { if (vertex == v) return flow; else if (vertex == w) return capacity - flow; else throw std::invalid_argument("invalid endpoint"); } void FlowEdge::addResidualFlowTo(const int vertex, const double delta) { if (!(delta >= 0.0)) throw std::invalid_argument("Delta must be non-negative"); if (vertex == v) flow -= delta; else if (vertex == w) flow += delta; else throw std::invalid_argument("invalid endpoint"); if (std::abs(flow) <= FLOATING_POINT_EPSILON) flow = 0; if (std::abs(flow - capacity) <= FLOATING_POINT_EPSILON) flow = capacity; if (!(flow >= 0.0)) throw std::invalid_argument("Flow is negative"); if (!(flow <= capacity)) throw std::invalid_argument("Flow exceeds capacity"); } std::ostream& operator<<(std::ostream& os, const FlowEdge& e) { os << e.v << "->" << e.w << " " << e.flow << "/" << e.capacity; return os; }

class FlowEdge: FLOATING_POINT_EPSILON = 1e-10 def __init__(self, v, w, capacity, flow=0.0): if v < 0: raise ValueError("vertex index must be a non-negative integer") if w < 0: raise ValueError("vertex index must be a non-negative integer") if capacity < 0.0: raise ValueError("Edge capacity must be non-negative") if flow > capacity: raise ValueError("flow exceeds capacity") if flow < 0.0: raise ValueError("flow must be non-negative") self._v = v self._w = w self._capacity = capacity self._flow = flow def from_(self): return self._v def to(self): return self._w def capacity(self): return self._capacity def flow(self): return self._flow def other(self, vertex): if vertex == self._v: return self._w elif vertex == self._w: return self._v else: raise ValueError("invalid endpoint") def residualCapacityTo(self, vertex): if vertex == self._v: return self._flow elif vertex == self._w: return self._capacity - self._flow else: raise ValueError("invalid endpoint") def addResidualFlowTo(self, vertex, delta): if delta < 0.0: raise ValueError("Delta must be non-negative") if vertex == self._v: self._flow -= delta elif vertex == self._w: self._flow += delta else: raise ValueError("invalid endpoint") if abs(self._flow) <= self.FLOATING_POINT_EPSILON: self._flow = 0 if abs(self._flow - self._capacity) <= self.FLOATING_POINT_EPSILON: self._flow = self._capacity if self._flow < 0.0: raise ValueError("Flow is negative") if self._flow > self._capacity: raise ValueError("Flow exceeds capacity") def __str__(self): return f"{self._v}->{self._w} {self._flow}/{self._capacity}"

18.5.2 Flow Network

import java.util.ArrayList; import java.util.List; public class FlowNetwork { private final int V; private int E; private final List<FlowEdge>[] adj; public FlowNetwork(int V, int E, List<int[]> edges) { if (V < 0) throw new IllegalArgumentException("Number of vertices in a Graph must be non-negative"); if (E < 0) throw new IllegalArgumentException("Number of edges must be non-negative"); this.V = V; this.E = 0; adj = (List<FlowEdge>[]) new List[V]; for (int v = 0; v < V; v++) adj[v] = new ArrayList<>(); for (int[] edge : edges) { int v = edge[0]; int w = edge[1]; double capacity = edge[2]; validateVertex(v); validateVertex(w); addEdge(new FlowEdge(v, w, capacity)); } } public int V() { return V; } public int E() { return E; } private void validateVertex(int v) { if (v < 0 || v >= V) throw new IllegalArgumentException("vertex " + v + " is not between 0 and " + (V-1)); } public void addEdge(FlowEdge e) { int v = e.from(); int w = e.to(); validateVertex(v); validateVertex(w); adj[v].add(e); adj[w].add(e); E++; } public Iterable<FlowEdge> adj(int v) { validateVertex(v); return adj[v]; } public Iterable<FlowEdge> edges() { List<FlowEdge> list = new ArrayList<>(); for (int v = 0; v < V; v++) for (FlowEdge e : adj(v)) { if (e.to() != v) list.add(e); } return list; } public String toString() { StringBuilder s = new StringBuilder(); s.append(V).append(" ").append(E).append(System.lineSeparator()); for (int v = 0; v < V; v++) { s.append(v).append(": "); for (FlowEdge e : adj[v]) { if (e.to() != v) s.append(e).append(" "); } s.append(System.lineSeparator()); } return s.toString(); } }

#ifndef FLOWNETWORK_H #define FLOWNETWORK_H #include <vector> #include "FlowEdge.h" class FlowNetwork { private: int V; int E; std::vector<FlowEdge> adj; void validateVertex(int v) const; public: FlowNetwork(int V, int E, const std::vector<std::vector<int>>& edges); [[nodiscard]] int getV() const; [[nodiscard]] int getE() const; void addEdge(const FlowEdge& e); [[nodiscard]] std::vector<FlowEdge> getadj(int v) const; [[nodiscard]] std::vector<FlowEdge> edges() const; friend std::ostream& operator<<(std::ostream& os, const FlowNetwork& network); }; #endif // FLOWNETWORK_H

#include <iostream> #include "FlowNetwork.h" FlowNetwork::FlowNetwork(int V, int E, const std::vector<std::vector<int>>& edges) : V(V), E(0) { if (V < 0) throw std::invalid_argument("Number of vertices in a Graph must be non-negative"); if (E > 0) throw std::invalid_argument("Number of edges must be non-negative"); adj = new std::vector<FlowEdge>[V]; for (const auto& edge : edges) { int v = edge[0]; int w = edge[1]; double capacity = edge[2]; validateVertex(v); validateVertex(w); addEdge(FlowEdge(v, w, capacity)); } } int FlowNetwork::V() const { return V; } int FlowNetwork::E() const { return E; } void FlowNetwork::validateVertex(int v) const { if (v < 0 || v >= V) throw std::invalid_argument("vertex " + std::to_string(v) + " is not between 0 and " + std::to_string(V - 1)); } void FlowNetwork::addEdge(const FlowEdge& e) { int v = e.from(); int w = e.to(); validateVertex(v); validateVertex(w); adj[v].push_back(e); adj[w].push_back(e); E++; } std::vector<FlowEdge> FlowNetwork::adj(int v) const { validateVertex(v); return adj[v]; } std::vector<FlowEdge> FlowNetwork::edges() const { std::vector<FlowEdge> list; for (int v = 0; v < V; v++) { for (const FlowEdge& e : adj(v)) { if (e.to() != v) list.push_back(e); } } return list; } std::ostream& operator<<(std::ostream& os, const FlowNetwork& network) { os << network.V << " " << network.E << std::endl; for (int v = 0; v < network.V; v++) { os << v << ": "; for (const FlowEdge& e : network.adj[v]) { if (e.to() != v) os << e << " "; } os << std::endl; } return os; }

from FlowEdge import FlowEdge class FlowNetwork: def __init__(self, V, E, edges): if V < 0: raise ValueError("Number of vertices in a Graph must be non-negative") if E < 0: raise ValueError("Number of edges must be non-negative") self._V = V self._E = 0 self._adj = [[] for _ in range(V)] for edge in edges: v, w, capacity = edge self._validate_vertex(v) self._validate_vertex(w) self._add_edge(FlowEdge(v, w, capacity)) def V(self): return self._V def E(self): return self._E def _validate_vertex(self, v): if v < 0 or v >= self._V: raise ValueError(f"vertex {v} is not between 0 and {self._V - 1}") def _add_edge(self, e): v = e.from_() w = e.to() self._validate_vertex(v) self._validate_vertex(w) self._adj[v].append(e) self._adj[w].append(e) self._E += 1 def adj(self, v): self._validate_vertex(v) return self._adj[v] def edges(self): all_edges = [] for v in range(self._V): for edge in self.adj(v): if edge.to() != v: all_edges.append(edge) return all_edges def __str__(self): s = f"{self._V} {self._E}\n" for v in range(self._V): s += f"{v}: " for edge in self._adj[v]: if edge.to() != v: s += str(edge) + " " s += "\n" return s

18.5.3 Ford-Fulkerson Algorithm

import java.util.LinkedList; import java.util.Queue; public class FordFulkerson { private static final double FLOATING_POINT_EPSILON = 1.0E-11; private final int V; private boolean[] marked; private FlowEdge[] edgeTo; private double value; public FordFulkerson(FlowNetwork G, int s, int t) { V = G.V(); validate(s); validate(t); if (s == t) throw new IllegalArgumentException("Source equals sink"); if (!isFeasible(G, s, t)) throw new IllegalArgumentException("Initial flow is infeasible"); value = excess(G, t); while (hasAugmentingPath(G, s, t)) { double bottle = Double.POSITIVE_INFINITY; for (int v = t; v != s; v = edgeTo[v].other(v)) { bottle = Math.min(bottle, edgeTo[v].residualCapacityTo(v)); } for (int v = t; v != s; v = edgeTo[v].other(v)) { edgeTo[v].addResidualFlowTo(v, bottle); } value += bottle; } assert check(G, s, t); } public double value() { return value; } public boolean inCut(int v) { validate(v); return marked[v]; } private void validate(int v) { if (v < 0 || v >= V) throw new IllegalArgumentException("vertex " + v + " is not between 0 and " + (V - 1)); } private boolean hasAugmentingPath(FlowNetwork G, int s, int t) { edgeTo = new FlowEdge[G.V()]; marked = new boolean[G.V()]; Queue<Integer> queue = new LinkedList<>(); queue.add(s); marked[s] = true; while (!queue.isEmpty() && !marked[t]) { int v = queue.remove(); for (FlowEdge e : G.adj(v)) { int w = e.other(v); if (e.residualCapacityTo(w) > 0) { if (!marked[w]) { edgeTo[w] = e; marked[w] = true; queue.add(w); } } } } return marked[t]; } private double excess(FlowNetwork G, int v) { double excess = 0.0; for (FlowEdge e : G.adj(v)) { if (v == e.from()) excess -= e.flow(); else excess += e.flow(); } return excess; } private boolean isFeasible(FlowNetwork G, int s, int t) { for (int v = 0; v < G.V(); v++) { for (FlowEdge e : G.adj(v)) { if (e.flow() < -FLOATING_POINT_EPSILON || e.flow() > e.capacity() + FLOATING_POINT_EPSILON) { System.err.println("Edge does not satisfy capacity constraints: " + e); return false; } } } if (Math.abs(value + excess(G, s)) > FLOATING_POINT_EPSILON) { System.err.println("Excess at source = " + excess(G, s)); System.err.println("Max flow = " + value); return false; } if (Math.abs(value - excess(G, t)) > FLOATING_POINT_EPSILON) { System.err.println("Excess at sink = " + excess(G, t)); System.err.println("Max flow = " + value); return false; } for (int v = 0; v < G.V(); v++) { if (v == s || v == t) continue; else if (Math.abs(excess(G, v)) > FLOATING_POINT_EPSILON) { System.err.println("Net flow out of " + v + " doesn't equal zero"); return false; } } return true; } private boolean check(FlowNetwork G, int s, int t) { if (!isFeasible(G, s, t)) { System.err.println("Flow is infeasible"); return false; } if (!inCut(s)) { System.err.println("source " + s + " is not on source side of min cut"); return false; } if (inCut(t)) { System.err.println("sink " + t + " is on source side of min cut"); return false; } double mincutValue = 0.0; for (int v = 0; v < G.V(); v++) { for (FlowEdge e : G.adj(v)) { if ((v == e.from()) && inCut(e.from()) && !inCut(e.to())) mincutValue += e.capacity(); } } if (Math.abs(mincutValue - value) > FLOATING_POINT_EPSILON) { System.err.println("Max flow value = " + value + ", min cut value = " + mincutValue); return false; } return true; } }

#ifndef FORDFULKERSON_H #define FORDFULKERSON_H #include <vector> #include "FlowEdge.h" #include "FlowNetwork.h" class FordFulkerson { private: static constexpr double FLOATING_POINT_EPSILON = 1.0E-11; int V; std::vector<bool> marked; std::vector<FlowEdge> edgeTo; double value; void validate(int v) const; bool hasAugmentingPath(const FlowNetwork& G, int s, int t); static double excess(const FlowNetwork& G, int v) ; [[nodiscard]] bool isFeasible(const FlowNetwork& G, int s, int t) const; [[nodiscard]] bool check(const FlowNetwork& G, int s, int t) const; public: FordFulkerson(const FlowNetwork& G, int s, int t); [[nodiscard]] double getvalue() const; [[nodiscard]] bool inCut(int v) const; }; #endif // FORDFULKERSON_H

#include <cassert> #include <iostream> #include <limits> #include <queue> #include <vector> #include "FordFulkerson.h" FordFulkerson::FordFulkerson(const FlowNetwork& G, const int s, const int t) : V(G.getV()), value(0.0) { validate(s); validate(t); if (s == t) throw std::invalid_argument("Source equals sink"); if (!isFeasible(G, s, t)) throw std::invalid_argument("Initial flow is infeasible"); value = excess(G, t); while (hasAugmentingPath(G, s, t)) { double bottle = std::numeric_limits<double>::infinity(); // Use numeric_limits for infinity for (int v = t; v != s; v = edgeTo[v].other(v)) { bottle = std::min(bottle, edgeTo[v].residualCapacityTo(v)); } for (int v = t; v != s; v = edgeTo[v].other(v)) { edgeTo[v].addResidualFlowTo(v, bottle); } value += bottle; } assert(check(G, s, t)); } double FordFulkerson::getvalue() const { return value; } bool FordFulkerson::inCut(int v) const { validate(v); return marked[v]; } void FordFulkerson::validate(int v) const { if (v < 0 || v >= V) throw std::invalid_argument("vertex " + std::to_string(v) + " is not between 0 and " + std::to_string(V - 1)); } bool FordFulkerson::hasAugmentingPath(const FlowNetwork& G, const int s, const int t) { edgeTo.assign(G.getV(), FlowEdge(0, 0, 0.0)); marked.assign(G.getV(), false); std::queue<int> queue; queue.push(s); marked[s] = true; while (!queue.empty() && !marked[t]) { const int v = queue.front(); queue.pop(); for (const FlowEdge& e : G.getadj(v)) { const int w = e.other(v); if (e.residualCapacityTo(w) > 0) { if (!marked[w]) { edgeTo[w] = e; marked[w] = true; queue.push(w); } } } } return marked[t]; } double FordFulkerson::excess(const FlowNetwork& G, const int v) { double excess = 0.0; for (const FlowEdge& e : G.getadj(v)) { if (v == e.from()) excess -= e.getflow(); else excess += e.getflow(); } return excess; } bool FordFulkerson::isFeasible(const FlowNetwork& G, int s, int t) const { for (int v = 0; v < G.getV(); v++) { for (const FlowEdge& e : G.getadj(v)) { if (e.getflow() < -FLOATING_POINT_EPSILON || e.getflow() > e.getcapacity() + FLOATING_POINT_EPSILON) { std::cerr << "Edge does not satisfy capacity constraints: " << e << std::endl; return false; } } } if (std::abs(value + excess(G, s)) > FLOATING_POINT_EPSILON) { std::cerr << "Excess at source = " << excess(G, s) << std::endl; std::cerr << "Max flow = " << value << std::endl; return false; } if (std::abs(value - excess(G, t)) > FLOATING_POINT_EPSILON) { std::cerr << "Excess at sink = " << excess(G, t) << std::endl; std::cerr << "Max flow = " << value << std::endl; return false; } for (int v = 0; v < G.getV(); v++) { if (v == s || v == t) continue; else if (std::abs(excess(G, v)) > FLOATING_POINT_EPSILON) { std::cerr << "Net flow out of " << v << " doesn't equal zero" << std::endl; return false; } } return true; } bool FordFulkerson::check(const FlowNetwork& G, int s, int t) const { if (!isFeasible(G, s, t)) { std::cerr << "Flow is infeasible" << std::endl; return false; } if (!inCut(s)) { std::cerr << "source " << s << " is not on source side of min cut" << std::endl; return false; } if (inCut(t)) { std::cerr << "sink " << t << " is on source side of min cut" << std::endl; return false; } double mincutValue = 0.0; for (int v = 0; v < G.getV(); v++) { for (const FlowEdge& e : G.getadj(v)) { if ((v == e.from()) && inCut(e.from()) && !inCut(e.to())) mincutValue += e.getcapacity(); } } if (std::abs(mincutValue - value) > FLOATING_POINT_EPSILON) { std::cerr << "Max flow value = " << value << ", min cut value = " << mincutValue << std::endl; return false; } return true; }

from collections import deque class FordFulkerson: FLOATING_POINT_EPSILON = 1e-11 def __init__(self, G, s, t): self._V = G.V() self._validate(s) self._validate(t) if s == t: raise ValueError("Source equals sink") if not self._is_feasible(G, s, t): raise ValueError("Initial flow is infeasible") self._value = self._excess(G, t) while self._has_augmenting_path(G, s, t): # compute bottleneck capacity bottle = float('inf') for v in range(t, s -1, -1): if v != s: bottle = min(bottle, self._edgeTo[v].residualCapacityTo(v)) v = self._edgeTo[v].other(v) # augment flow for v in range(t, s-1, -1): if v != s: self._edgeTo[v].addResidualFlowTo(v, bottle) v = self._edgeTo[v].other(v) self._value += bottle # check optimality conditions assert self._check(G, s, t) def value(self): return self._value def in_cut(self, v): self._validate(v) return self._marked[v] def _validate(self, v): if v < 0 or v >= self._V: raise ValueError(f"vertex {v} is not between 0 and {self._V - 1}") def _has_augmenting_path(self, G, s, t): self._edgeTo = [None] * G.V() self._marked = [False] * G.V() queue = deque() queue.append(s) self._marked[s] = True while queue and not self._marked[t]: v = queue.popleft() for e in G.adj(v): w = e.other(v) if e.residualCapacityTo(w) > 0: if not self._marked[w]: self._edgeTo[w] = e self._marked[w] = True queue.append(w) return self._marked[t] def _excess(self, G, v): excess = 0.0 for e in G.adj(v): if v == e.from_(): excess -= e.flow() else: excess += e.flow() return excess def _is_feasible(self, G, s, t): for v in range(G.V()): for e in G.adj(v): if e.flow() < -self.FLOATING_POINT_EPSILON or e.flow() > e.capacity() + self.FLOATING_POINT_EPSILON: print(f"Edge does not satisfy capacity constraints: {e}") return False if abs(self._value + self._excess(G, s)) > self.FLOATING_POINT_EPSILON: print(f"Excess at source = {self._excess(G, s)}") print(f"Max flow = {self._value}") return False if abs(self._value - self._excess(G, t)) > self.FLOATING_POINT_EPSILON: print(f"Excess at sink = {self._excess(G, t)}") print(f"Max flow = {self._value}") return False for v in range(G.V()): if v == s or v == t: continue elif abs(self._excess(G, v)) > self.FLOATING_POINT_EPSILON: print(f"Net flow out of {v} doesn't equal zero") return False return True def _check(self, G, s, t): if not self._is_feasible(G, s, t): print("Flow is infeasible") return False if not self.in_cut(s): print(f"source {s} is not on source side of min cut") return False if self.in_cut(t): print(f"sink {t} is on source side of min cut") return False mincut_value = 0.0 for v in range(G.V()): for e in G.adj(v): if v == e.from_() and self.in_cut(e.from_()) and not self.in_cut(e.to()): mincut_value += e.capacity() if abs(mincut_value - self._value) > self.FLOATING_POINT_EPSILON: print(f"Max flow value = {self._value}, min cut value = {mincut_value}") return False return True

18.6 Maxflow Applications

Applications

Data mining.
Open-pit mining.
Bipartite matching.
Network reliability.
Baseball elimination.
Image segmentation.
Network connectivity.
Distributed computing.
Security of statistical data.
Egalitarian stable matching.
Multi-camera scene reconstruction.
Sensor placement for homeland security.
Many, many, more.

18.6.1 Bipartite Matching

N students apply for N jobs. Given a bipartite graph, find a perfect matching.

Bipartite Matching

Create , , one vertex for each student, and one vertex for each job.
Add edge from to each student (capacity 1).
Add edge from each job to (capacity 1).
Add edge from student to each job offered (infinite capacity).
1-1 correspondence between perfect matchings in bipartite graph and integer-valued maxflows of value .

When no perfect matching, mincut explains why

Consider mincut (, ):

Let = students on side of cut.
Let = companies on side of cut.
Fact: ; students in can be matched only to companies in .

18.6.2 Baseball Elimination

Problem: Which teams have a chance of finishing the season with the most wins?

Team 4 not eliminated iff all edges pointing from s are full in maxflow.

19 Radix Sort

19.1 Strings in Java

String: Sequence of characters.

The char data type

C char data type: Typically an 8-bit integer.
- Supports 7-bit ASCII.
- Can represent only 256 characters.
Java char data type: A 16-bit unsigned integer.
- Supports original 16-bit Unicode.
- Supports 21-bit Unicode 3.0 (awkwardly).

Examples

public class StringTest { public static void main(String[] args) { String s1 = "Hello"; System.out.println(s1.length()); // 5 System.out.println(s1.charAt(0)); // H System.out.println(s1.substring(0, 1)); // H System.out.println(s1.concat(" World")); // Hello World } }

		String		StringBuilder
Data Type in Java		Sequence of characters (immutable)		Sequence of characters (mutable)
Underlying Implementation		Immutable `char[]` array, offset, and length		Resizing `char[]` array and length.
		Guarantee	Extra Space	Guarantee	Extra Space
Operation	`length()`
	`charAt()`
	`substring()`
	`concat()`

19.2 Key-Indexed Counting

Assumption: Keys are integers between and . Can use key as an array index.

Goal: Sort an array of integers between and .

Key-Indexed Counting

Count frequencies of each letter using key as index.
Compute frequency cumulates which specify destinations.
Access cumulates using key as index to move items.
Copy back into original array.

Properties

Key-indexed counting uses array accesses to sort items whose keys are integers between and .
Key-indexed counting uses extra space proportional to .
Key-indexed counting is stable.

public class KeyIndexedSorting { public static void sort(Squirrel[] students) { int N = students.length; int R = 5; // Assuming grades are from 0 to 4 int[] count = new int[R + 1]; for (Squirrel student : students) { count[student.grade + 1]++; } for (int r = 0; r < R; r++) { count[r + 1] += count[r]; } Squirrel[] aux = new Squirrel[N]; for (Squirrel student : students) { aux[count[student.grade]++] = student; } System.arraycopy(aux, 0, students, 0, N); } static class Squirrel { String name; int grade; public Squirrel(String name, int grade) { this.name = name; this.grade = grade; } @Override public String toString() { return name + " (Grade: " + grade + ")"; } } }

#include <iostream> #include <utility> #include <vector> #include <string> struct Squirrel { std::string name; int grade; Squirrel(std::string n, const int g) : name(std::move(n)), grade(g) {} friend std::ostream& operator<<(std::ostream& os, const Squirrel& s) { os << s.name << " (Grade: " << s.grade << ")"; return os; } }; void sort(std::vector<Squirrel>& students) { const int N = static_cast<int>(students.size()); constexpr int R = 5; // Assuming grades are from 0 to 4 std::vector<int> count(R + 1, 0); for (int i = 0; i < N; i++) { count[students[i].grade + 1]++; } for (int r = 0; r < R; r++) { count[r + 1] += count[r]; } std::vector<Squirrel> aux(N); for (int i = 0; i < N; i++) { aux[count[students[i].grade]++] = students[i]; } students = aux; }

class Squirrel: def __init__(self, name, grade): self.name = name self.grade = grade def __str__(self): return f"{self.name} (Grade: {self.grade})" def sort(students): N = len(students) R = 5 # Assuming grades are from 0 to 4 count = [0] * (R + 1) for student in students: count[student.grade + 1] += 1 for r in range(R): count[r + 1] += count[r] aux = [None] * N for student in students: aux[count[student.grade]] = student count[student.grade] += 1 students[:] = aux

19.3 LSD Radix Sort

LSD Radix Sort

Consider characters from right to left.
Stably sort using dth character as the key (using key-indexed counting).

Correctness Proof: LSD sorts fixed-length strings in ascending order.

Proof

After pass , strings are sorted by last characters.

If two strings differ on sort key, key-indexed sort puts them in proper relative order.
If two strings agree on sort key, stability keeps them in proper relative order.

Property: LSD sort is stable.

Google (or presidential) Interview Question: Sort one million 32-bit integers using LSD radix sort.

public class LSDStringSort { public static void sort(String[] a, int W) { int N = a.length; int R = 256; // extended ASCII alphabet size String[] aux = new String[N]; for (int d = W - 1; d >= 0; d--) { int[] count = new int[R + 1]; for (String string : a) { count[string.charAt(d) + 1]++; } for (int r = 0; r < R; r++) { count[r + 1] += count[r]; } for (String s : a) { aux[count[s.charAt(d)]++] = s; } System.arraycopy(aux, 0, a, 0, N); } } }

#include <iostream> #include <string> #include <vector> void lsdSort(std::vector<std::string>& a, const int w) { const int n = static_cast<int>(a.size()); int R = 256; std::vector<std::string> aux(n); for (int d = w - 1; d >= 0; d--) { std::vector<int> count(R + 1, 0); for (int i = 0; i < n; i++) { count[a[i][d] + 1]++; } for (int r = 0; r < R; r++) { count[r + 1] += count[r]; } for (int i = 0; i < n; i++) { aux[count[a[i][d]]++] = a[i]; } for (int i = 0; i < n; i++) { a[i] = aux[i]; } } }

def lsd_sort(a, w): n = len(a) R = 256 aux = [""] * n for d in range(w - 1, -1, -1): count = [0] * (R + 1) for i in range(n): count[ord(a[i][d]) + 1] += 1 for r in range(R): count[r + 1] += count[r] for i in range(n): aux[count[ord(a[i][d])]] = a[i] count[ord(a[i][d])] += 1 a[:] = aux

19.4 MSD Radix Sort

MSD Radix Sort

Partition array into pieces according to first character (use key-indexed counting ).
Recursively sort all strings that start with each character (key-indexed counts delineate subarrays to sort).

Variable-length strings: Treat strings as if they had an extra char at end (smaller than any char)

Potential for Disastrous Performance

Much too slow for small subarrays.
- Each function call needs its own count[] array.
- ASCII (256 counts): 100x slower than copy pass for .
- Unicode (65,536 counts): 32,000x slower for .
Huge number of small subarrays because of recursion.

Improvements

Cutoff to insertion sort for small subarrays.

Insertion sort, but start at character.
Implement less() so that it compares starting at character.

Performance

Number of characters examined.

MSD examines just enough characters to sort the keys.
Number of characters examined depends on keys.
Can be sublinear in input size!

MSD String Sort vs. Quicksort for Strings

Disadvantages of MSD string sort:
- Extra space for aux[] arrays.
- Extra space for count[] arrays.
- Inner loop has a lot of instructions.
- Accesses memory "randomly" (cache inefficient).
Disadvantage of quicksort
- Linearithmic number of string compares (not linear).
- Has to rescan many characters in keys with long prefix matches .

public class MSDStringSort { private static final int R = 256; private static final int CUTOFF = 15; public static void sort(String[] a) { String[] aux = new String[a.length]; sort(a, aux, 0, a.length - 1, 0); } private static void sort(String[] a, String[] aux, int low, int high, int d) { if (high <= low + CUTOFF) { insertionSort(a, low, high, d); return; } int[] count = new int[R + 2]; for (int i = low; i <= high; i++) { int c = charAt(a[i], d); count[c + 2]++; } for (int r = 0; r < R + 1; r++) { count[r + 1] += count[r]; } for (int i = low; i <= high; i++) { int c = charAt(a[i], d); aux[count[c + 1]++] = a[i]; } if (high + 1 - low >= 0) System.arraycopy(aux, 0, a, low, high + 1 - low); for (int r = 0; r < R; r++) { sort(a, aux, low + count[r], low + count[r + 1] - 1, d + 1); } } private static int charAt(String s, int d) { if (d < s.length()) { return s.charAt(d); } else { return -1; } } private static void insertionSort(String[] a, int low, int high, int d) { for (int i = low; i <= high; i++) { for (int j = i; j > low && less(a[j], a[j - 1], d); j--) { swap(a, j, j - 1); } } } private static boolean less(String v, String w, int d) { return v.substring(d).compareTo(w.substring(d)) < 0; } private static void swap(String[] a, int i, int j) { String temp = a[i]; a[i] = a[j]; a[j] = temp; } }

#include <vector> #include <string> #include <algorithm> #include <iostream> class MSDStringSort { private: static constexpr int R = 256; static constexpr int CUTOFF = 15; public: static void sort(std::vector<std::string>& a) { std::vector<std::string> aux(a.size()); sort(a, aux, 0, static_cast<int>(a.size()) - 1, 0); } private: static void sort(std::vector<std::string>& a, std::vector<std::string>& aux, const int low, const int high, const int d) { if (high <= low + CUTOFF) { insertionSort(a, low, high, d); return; } std::vector<int> count(R + 2, 0); for (int i = low; i <= high; i++) { const int c = charAt(a[i], d); count[c + 2]++; } for (int r = 0; r < R + 1; r++) { count[r + 1] += count[r]; } for (int i = low; i <= high; i++) { const int c = charAt(a[i], d); aux[count[c + 1]++] = a[i]; } std::copy_n(aux.begin(), (high + 1 - low), a.begin() + low); for (int r = 0; r < R; r++) { sort(a, aux, low + count[r], low + count[r + 1] - 1, d + 1); } } static int charAt(const std::string& s, const int d) { if (d < s.length()) { return s[d]; } else { return -1; } } static void insertionSort(std::vector<std::string>& a, const int low, const int high, const int d) { for (int i = low; i <= high; i++) { for (int j = i; j > low && less(a[j], a[j - 1], d); j--) { std::swap(a[j], a[j - 1]); } } } static bool less(const std::string& v, const std::string& w, const int d) { return v.substr(d) < w.substr(d); } };

def char_at(s, d): if d < len(s): return ord(s[d]) else: return -1 def insertion_sort(arr, low, high, d): for i in range(low, high + 1): for j in range(i, low, -1): if arr[j][d:] < arr[j - 1][d:]: arr[j], arr[j - 1] = arr[j - 1], arr[j] else: break def msd_string_sort(arr): CUTOFF = 15 aux = [None] * len(arr) sort(arr, 0, len(arr) - 1, 0, aux, CUTOFF) def sort(arr, low, high, d, aux, CUTOFF): if high <= low + CUTOFF: insertion_sort(arr, low, high, d) return R = 256 count = [0] * (R + 2) for i in range(low, high + 1): c = char_at(arr[i], d) count[c + 2] += 1 for r in range(R + 1): count[r + 1] += count[r] for i in range(low, high + 1): c = char_at(arr[i], d) aux[count[c + 1]] = arr[i] count[c + 1] += 1 for i in range(low, high + 1): arr[i] = aux[i - low] for r in range(R): sort(arr, low + count[r], low + count[r + 1] - 1, d + 1, aux, CUTOFF)

19.5 3-Way Radix Quicksort

Do 3-way partitioning on the character.

Properties

Less overhead than R-way partitioning in MSD string sort.
Does not re-examine characters equal to the partitioning char.
(but does re-examine characters not equal to the partitioning char)

3-Way String Quicksort vs. Standard Quicksort

Standard quicksort
- Uses string compares on average.
- Costly for keys with long common prefixes (and this is a common case!)
3-way string (radix) quicksort
- Uses character compares on average for random strings.
- Avoids re-comparing long common prefixes.

3-way String Quicksort vs. MSD String Sort

MSD string sort
- Is cache-inefficient.
- Too much memory storing count[].
- Too much overhead reinitializing count[] and aux[].
3-way string quicksort
- Has a short inner loop.
- Is cache-friendly.
- Is in-place.

public class ThreeWayRadixQuicksortStrings { private static final int CUTOFF = 15; private static int charAt(String s, int d) { if (d < s.length()) return s.charAt(d); else return -1; } public static void sort(String[] a) { sort(a, 0, a.length - 1, 0); } private static void sort(String[] a, int lo, int hi, int d) { if (hi <= lo + CUTOFF) { insertionSort(a, lo, hi, d); return; } int lt = lo, gt = hi; int v = charAt(a[lo], d); int i = lo + 1; while (i <= gt) { int t = charAt(a[i], d); if (t < v) exch(a, lt++, i++); else if (t > v) exch(a, i, gt--); else i++; } sort(a, lo, lt - 1, d); if (v >= 0) sort(a, lt, gt, d + 1); sort(a, gt + 1, hi, d); } private static void insertionSort(String[] a, int lo, int hi, int d) { for (int i = lo; i <= hi; i++) { for (int j = i; j > lo && less(a[j], a[j - 1], d); j--) { exch(a, j, j - 1); } } } private static boolean less(String v, String w, int d) { for (int i = d; i < Math.min(v.length(), w.length()); i++) { if (v.charAt(i) < w.charAt(i)) return true; if (v.charAt(i) > w.charAt(i)) return false; } return v.length() < w.length(); } private static void exch(String[] a, int i, int j) { String temp = a[i]; a[i] = a[j]; a[j] = temp; } }

#include <iostream> #include <vector> #include <string> class ThreeWayRadixQuicksortStrings { private: static constexpr int CUTOFF = 15; static int charAt(const std::string& s, const int d) { if (d < s.length()) return s[d]; else return -1; } static void sort(std::vector<std::string>& a, const int lo, const int hi, const int d) { if (hi <= lo + CUTOFF) { insertionSort(a, lo, hi, d); return; } int lt = lo, gt = hi; const int v = charAt(a[lo], d); int i = lo + 1; while (i <= gt) { int t = charAt(a[i], d); if (t < v) std::swap(a[lt++], a[i++]); else if (t > v) std::swap(a[i], a[gt--]); else i++; } sort(a, lo, lt - 1, d); if (v >= 0) sort(a, lt, gt, d + 1); sort(a, gt + 1, hi, d); } static void insertionSort(std::vector<std::string>& a, const int lo, const int hi, const int d) { for (int i = lo; i <= hi; i++) { for (int j = i; j > lo && less(a[j], a[j - 1], d); j--) { std::swap(a[j], a[j - 1]); } } } static bool less(const std::string& v, const std::string& w, const int d) { for (int i = d; i < std::min(v.length(), w.length()); i++) { if (v[i] < w[i]) return true; if (v[i] > w[i]) return false; } return v.length() < w.length(); } public: static void sort(std::vector<std::string>& a) { sort(a, 0, static_cast<int>(a.size()) - 1, 0); } };

CUTOFF = 15 def char_at(s, d): if d < len(s): return ord(s[d]) else: return -1 def insertion_sort(arr, lo, hi, d): for i in range(lo, hi + 1): for j in range(i, lo, -1): if arr[j][d:] < arr[j - 1][d:]: arr[j], arr[j - 1] = arr[j - 1], arr[j] else: break def three_way_radix_quicksort(arr): def sort(arr, lo, hi, d): if hi <= lo + CUTOFF: insertion_sort(arr, lo, hi, d) return lt, gt = lo, hi v = char_at(arr[lo], d) i = lo + 1 while i <= gt: t = char_at(arr[i], d) if t < v: arr[lt], arr[i] = arr[i], arr[lt] lt += 1 i += 1 elif t > v: arr[gt], arr[i] = arr[i], arr[gt] gt -= 1 else: i += 1 sort(arr, lo, lt - 1, d) if v >= 0: sort(arr, lt, gt, d + 1) sort(arr, gt + 1, hi, d) sort(arr, 0, len(arr) - 1, 0)

Summary of the Performance of Sorting Algorithms

Algorithm	Guarantee	Stable?	Operations on keys
Insertion Sort		Yes	`compareTo()`
Mergesort		Yes	`compareTo()`
Quicksort	*	No	`compareTo()`
Heapsort		No	`compareTo()`
LSD☆		Yes	`charAt()`
MSD★		Yes	`charAt()`
3-Way String Quicksort	*	No	`charAt()`

*: probabilistic

☆: fixed-length W keys

★: average-length W keys

19.6 Suffix Arrays

Keyword in context search: Given a text of N characters, preprocess it to enable fast substring search (find all occurrences of query string context).

Applications: Linguistics, databases, web search, word processing, ...

Keyword-in-context search: suffix-sorting solution.

Preprocess: suffix sort the text.
Query: binary search for query; scan until mismatch.

Last modified: 25 November 2024