From 10e26931228de940a0fec6f01279c8c80c86b023 Mon Sep 17 00:00:00 2001 From: Kevin Lim Date: Tue, 8 Mar 2022 14:17:33 +0800 Subject: [PATCH 1/3] optimise java max flow template and add java min cost max flow template --- ch8/maxflow.java | 226 ++++++++++++++++++++++------------------ ch9/MinCostMaxFlow.java | 127 ++++++++++++++++++++++ 2 files changed, 250 insertions(+), 103 deletions(-) create mode 100644 ch9/MinCostMaxFlow.java diff --git a/ch8/maxflow.java b/ch8/maxflow.java index 58ed20f..8bf55d6 100644 --- a/ch8/maxflow.java +++ b/ch8/maxflow.java @@ -1,55 +1,54 @@ import java.io.*; import java.util.*; -// (Nearly) 1-1 translation of maxflow.cpp with min-cut augmentation -class Pair { - int first, second; - - Pair(int first, int second) { - this.first = first; - this.second = second; - } -} - class Edge { - int v; - long w, f; + int u, v; + long f; - Edge(int v, long w, long f) { + Edge(int u, int v, long f) { + this.u = u; this.v = v; - this.w = w; this.f = f; } + + @Override + public String toString() { + return u + " " + v + " " + f; + } } -public class maxflow { - static long INF = 1000000000000000000l; - - int V; - Edge[] EL; - int ELSize; - List> AL; - int[] d, last; - Pair[] p; - - maxflow(int initialV, int numEdges) { - V = initialV; - EL = new Edge[2 * numEdges]; // 2 * to account for back edges - ELSize = 0; - AL = new ArrayList<>(); - for (int i = 0; i < V; i++) { - AL.add(new ArrayList<>()); - } +// Partially adapted from https://github.com/lewin/Algorithms/blob/master/src/Graph/Algorithms/NetworkFlow/Dinic.java +// Contains additional functions from https://github.com/stevenhalim/cpbook-code/blob/master/ch8/maxflow.cpp +class maxflow { + static long INF = 1l << 62; // 4E18 + + int numVertices, numEdges; + int edgeIndex; + int[] adjNodes, prevEdges, lastEdges, d, last, bfsVertices, bfsEdges; + long[] flows, caps; + + maxflow(int numVertices, int numEdges) { + this.numVertices = numVertices; + numEdges *= 2; // * 2 to account for back edges + this.numEdges = numEdges; + this.adjNodes = new int[numEdges]; + this.prevEdges = new int[numEdges]; + this.lastEdges = new int[numVertices]; + this.flows = new long[numEdges]; + this.caps = new long[numEdges]; + this.edgeIndex = 0; + this.last = new int[numVertices]; + Arrays.fill(lastEdges, -1); } void addEdge(int u, int v, long w, boolean directed) { if (u == v) return; // safeguard: no self loop - EL[ELSize] = new Edge(v, w, 0); // u->v, cap w, flow 0 - ELSize++; - AL.get(u).add(ELSize - 1); // remember this index - EL[ELSize] = new Edge(u, directed ? 0 : w, 0); // back edge - ELSize++; - AL.get(v).add(ELSize - 1); // remember this index + + adjNodes[edgeIndex] = v; caps[edgeIndex] = w; flows[edgeIndex] = 0; // u->v, cap w, flow 0 + prevEdges[edgeIndex] = lastEdges[u]; lastEdges[u] = edgeIndex++; + + adjNodes[edgeIndex] = u; caps[edgeIndex] = directed ? 0 : w; flows[edgeIndex] = 0; // back edge + prevEdges[edgeIndex] = lastEdges[v]; lastEdges[v] = edgeIndex++; } long edmondsKarp(int s, int t) { @@ -65,7 +64,7 @@ long edmondsKarp(int s, int t) { long dinic(int s, int t) { long mf = 0; // mf stands for max_flow while (BFS(s, t)) { // an O(V^2*E) algorithm - last = new int[V]; // important speedup + System.arraycopy(lastEdges, 0, last, 0, numVertices); // important speedup long f; while ((f = DFS(s, t, INF)) > 0) { // exhaust blocking flow mf += f; @@ -75,98 +74,119 @@ long dinic(int s, int t) { } boolean BFS(int s, int t) { // find augmenting path - d = new int[V]; - p = new Pair[V]; // record BFS sp tree - for (int i = 0; i < V; i++) { + d = new int[numVertices]; + bfsVertices = new int[numVertices]; // record BFS sp tree + bfsEdges = new int[numVertices]; + for (int i = 0; i < numVertices; i++) { d[i] = -1; - p[i] = new Pair(-1, -1); + bfsVertices[i] = bfsEdges[i] = -1; } + d[s] = 0; - Queue q = new LinkedList<>(); - q.offer(s); - while (!q.isEmpty()) { - int u = q.poll(); + int[] queue = new int[numVertices]; + int front = 0, back = 0; + queue[back++] = s; + + while (front < back) { + int u = queue[front++]; if (u == t) break; // stop as sink t reached - for (int idx : AL.get(u)) { // explore neighbors of u - Edge e = EL[idx]; // stored in EL[idx] - int v = e.v; - long cap = e.w; - long flow = e.f; - if ((cap - flow > 0) && (d[v] == -1)) { // positive residual edge + + for (int idx = lastEdges[u]; idx != -1; idx = prevEdges[idx]) { // explore neighbors of u + int v = adjNodes[idx]; // stored in adjNodes[idx] + + if (d[v] == -1 & flows[idx] < caps[idx]) { d[v] = d[u] + 1; - q.offer(v); - p[v] = new Pair(u, idx); + queue[back++] = v; + bfsVertices[v] = u; + bfsEdges[v] = idx; } } } + return d[t] != -1; } - // Run after performing Dinics/Edmonds Karp to get nodes in min-cut - // Basically performs BFS on the flow graph one more time + long sendOneFlow(int s, int t, long f) { // send one flow from s->t + if (s == t) return f; // bottleneck edge f found + + int u = bfsVertices[t]; + int idx = bfsEdges[t]; + + long pushed = sendOneFlow(s, u, Math.min(f, caps[idx] - flows[idx])); + flows[idx] += pushed; + flows[idx ^ 1] -= pushed; // back flow + return pushed; + } + + long DFS(int u, int t, long f) { // traverse from u->t + if ((u == t) || (f == 0)) return f; + + for (int idx = last[u]; idx != -1; last[u] = idx = prevEdges[idx]) { // from last edge + int v = adjNodes[idx]; + if (d[v] != d[u] + 1) continue; // not part of layer graph + + long pushed = DFS(v, t, Math.min(f, caps[idx] - flows[idx])); + + if (pushed > 0) { + flows[idx] += pushed; + flows[idx ^ 1] -= pushed; // back flow + return pushed; + } + } + + return 0; + } + + // Run after performing Dinics/Edmonds Karp to get + // nodes in min cut List minCut(int s, int t) { List result = new ArrayList<>(); - d = new int[V]; - p = new Pair[V]; // record BFS sp tree - for (int i = 0; i < V; i++) { + d = new int[numVertices]; + bfsVertices = new int[numVertices]; // record BFS sp tree + bfsEdges = new int[numVertices]; + for (int i = 0; i < numVertices; i++) { d[i] = -1; - p[i] = new Pair(-1, -1); + bfsVertices[i] = bfsEdges[i] = -1; } + d[s] = 0; - Queue q = new LinkedList<>(); - q.offer(s); + int[] queue = new int[numVertices]; + int front = 0, back = 0; + queue[back++] = s; result.add(s); - while (!q.isEmpty()) { - int u = q.poll(); + + while (front < back) { + int u = queue[front++]; if (u == t) break; // stop as sink t reached - for (int idx : AL.get(u)) { // explore neighbors of u - Edge e = EL[idx]; // stored in EL[idx] - int v = e.v; - long cap = e.w; - long flow = e.f; - if ((cap - flow > 0) && (d[v] == -1)) { // positive residual edge + + for (int idx = lastEdges[u]; idx != -1; idx = prevEdges[idx]) { // explore neighbors of u + int v = adjNodes[idx]; // stored in adjNodes[idx] + + if (d[v] == -1 & flows[idx] < caps[idx]) { d[v] = d[u] + 1; - q.offer(v); - p[v] = new Pair(u, idx); + queue[back++] = v; + bfsVertices[v] = u; + bfsEdges[v] = idx; result.add(v); } } } - return result; - } - long sendOneFlow(int s, int t, long f) { // send one flow from s->t - if (s == t) return f; // bottleneck edge f found - Pair pair = p[t]; - int u = pair.first; - int idx = pair.second; - Edge e = EL[idx]; - long cap = e.w; - long flow = e.f; - long pushed = sendOneFlow(s, u, Math.min(f, cap - flow)); - e.f += pushed; - EL[idx ^ 1].f -= pushed; // back flow - return pushed; + return result; } + + // Run after performing Dinics/Edmonds Karp to get + // edges in flow graph + List getFlowGraphEdges() { + List result = new ArrayList<>(); - long DFS(int u, int t, long f) { // traverse from s->t - if ((u == t) || (f == 0)) return f; - int start = last[u]; - int stop = AL.get(u).size(); - for (int i = start; i < stop; i++) { // from last edge - Edge e = EL[AL.get(u).get(i)]; - int v = e.v; - long cap = e.w; - long flow = e.f; - if (d[v] != d[u] + 1) continue; // not part of layer graph - long pushed; - if ((pushed = DFS(v, t, Math.min(f, cap - flow))) > 0) { - e.f += pushed; - EL[AL.get(u).get(i) ^ 1].f -= pushed; // back flow - return pushed; + for (int i = 0; i < numEdges; i++) { + if (flows[i] > 0) { + result.add(new Edge(adjNodes[i ^ 1], adjNodes[i], flows[i])); } } - return 0; + + return result; } } diff --git a/ch9/MinCostMaxFlow.java b/ch9/MinCostMaxFlow.java new file mode 100644 index 0000000..ff137ff --- /dev/null +++ b/ch9/MinCostMaxFlow.java @@ -0,0 +1,127 @@ +import java.io.*; +import java.util.*; + +class Pair { + long flow, cost; + + Pair(long flow, long cost) { + this.flow = flow; + this.cost = cost; + } + + @Override + public String toString() { + return flow + " " + cost; + } +} + +// Partially adapted from https://github.com/lewin/Algorithms/blob/master/src/Graph/Algorithms/NetworkFlow/Dinic.java +// Contains additional functions from https://github.com/stevenhalim/cpbook-code/blob/master/ch8/maxflow.cpp +class MinCostMaxFlow { + static long INF = (long) 1e18; // INF = 1e18, not 2^63-1 to avoid overflow + + int numVertices, numEdges; + int edgeIndex; + int[] adjNodes, prevEdges, lastEdges, last; + long[] flows, caps, costs, d; + long totalCost; + boolean[] vis; + + MinCostMaxFlow(int numVertices, int numEdges) { + this.numVertices = numVertices; + numEdges *= 2; // * 2 to account for back edges + this.numEdges = numEdges; + this.adjNodes = new int[numEdges]; + this.prevEdges = new int[numEdges]; + this.lastEdges = new int[numVertices]; + this.flows = new long[numEdges]; + this.caps = new long[numEdges]; + this.costs = new long[numEdges]; + this.edgeIndex = 0; + this.last = new int[numVertices]; + this.vis = new boolean[numVertices]; + Arrays.fill(lastEdges, -1); + } + + void addEdge(int u, int v, long w, long c, boolean directed) { + if (u == v) return; // safeguard: no self loop + + adjNodes[edgeIndex] = v; caps[edgeIndex] = w; costs[edgeIndex] = c; flows[edgeIndex] = 0; // u->v, cap w, cost c, flow 0 + prevEdges[edgeIndex] = lastEdges[u]; lastEdges[u] = edgeIndex++; + + adjNodes[edgeIndex] = u; caps[edgeIndex] = directed ? 0 : w; costs[edgeIndex] = -c; flows[edgeIndex] = 0; // back edge + prevEdges[edgeIndex] = lastEdges[v]; lastEdges[v] = edgeIndex++; + } + + Pair mcmf(int s, int t) { + long mf = 0; // mf stands for max_flow + + while (SPFA(s, t)) { // an O(V^2*E) algorithm + System.arraycopy(lastEdges, 0, last, 0, numVertices); // important speedup + long f; + while ((f = DFS(s, t, INF)) > 0) { // exhaust blocking flow + mf += f; + } + } + + return new Pair(mf, totalCost); + } + + boolean SPFA(int s, int t) { // SPFA to find augmenting path in residual graph + d = new long[numVertices]; + Arrays.fill(d, INF); + + d[s] = 0; + vis[s] = true; + Queue queue = new LinkedList<>(); + queue.offer(s); + + while (!queue.isEmpty()) { + int u = queue.poll(); + vis[u] = false; + + for (int idx = lastEdges[u]; idx != -1; idx = prevEdges[idx]) { // explore neighbors of u + int v = adjNodes[idx]; // stored in adjNodes[idx] + + if (caps[idx] - flows[idx] > 0 && d[v] > d[u] + costs[idx]) { // positive residual edge + d[v] = d[u] + costs[idx]; + if (!vis[v]) { + queue.offer(v); + vis[v] = true; + } + } + } + } + + return d[t] != INF; // has an augmenting path + } + + long DFS(int u, int t, long f) { // traverse from u->t + if ((u == t) || (f == 0)) return f; + + vis[u] = true; + + for (int idx = last[u]; idx != -1; last[u] = idx = prevEdges[idx]) { // from last edge + int v = adjNodes[idx]; + long cap = caps[idx]; + long flow = flows[idx]; + long cost = costs[idx]; + + if (!vis[v] && d[v] == d[u] + cost) { // in current layer graph + long pushed; + + if ((pushed = DFS(v, t, Math.min(f, cap - flow))) > 0) { + totalCost += pushed * cost; + flows[idx] += pushed; + flows[idx ^ 1] -= pushed; // back edge + vis[u] = false; + return pushed; + } + } + } + + vis[u] = false; + return 0; + } +} + From 5efe9758976f45448e1e3e6dedb9a8999e9e8087 Mon Sep 17 00:00:00 2001 From: Kevin Lim Date: Tue, 8 Mar 2022 14:20:08 +0800 Subject: [PATCH 2/3] update min cost max flow java template --- ch9/MinCostMaxFlow.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ch9/MinCostMaxFlow.java b/ch9/MinCostMaxFlow.java index ff137ff..39c1858 100644 --- a/ch9/MinCostMaxFlow.java +++ b/ch9/MinCostMaxFlow.java @@ -16,7 +16,7 @@ public String toString() { } // Partially adapted from https://github.com/lewin/Algorithms/blob/master/src/Graph/Algorithms/NetworkFlow/Dinic.java -// Contains additional functions from https://github.com/stevenhalim/cpbook-code/blob/master/ch8/maxflow.cpp +// Contains additional functions from https://github.com/stevenhalim/cpbook-code/blob/master/ch9/mcmf.cpp class MinCostMaxFlow { static long INF = (long) 1e18; // INF = 1e18, not 2^63-1 to avoid overflow From 628fb31ffc7dea98b44c45090f9206c97311c4ad Mon Sep 17 00:00:00 2001 From: Kevin Lim Date: Tue, 12 Apr 2022 14:45:49 +0800 Subject: [PATCH 3/3] add java hungrarian algo --- ch9/Hungrarian.java | 168 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 168 insertions(+) create mode 100644 ch9/Hungrarian.java diff --git a/ch9/Hungrarian.java b/ch9/Hungrarian.java new file mode 100644 index 0000000..43a94b8 --- /dev/null +++ b/ch9/Hungrarian.java @@ -0,0 +1,168 @@ +import java.util.*; + +/* +Adapted from: https://github.com/ShahjalalShohag/code-library/blob/master/Graph%20Theory/Hungarian%20Algorithm.cpp +Complexity: O(n^3) but optimized +It finds minimum cost maximum matching. +For finding maximum cost maximum matching +add -cost and return -matching() +1-indexed +*/ +class Hungrarian { + static int N = 2005; + static long inf = (long) 1e18; + + long[][] c; + long[] fx, fy, d; + int[] l, r, arg, trace; + Queue q; + int start, finish, n; + + Hungrarian() { + c = new long[N][N]; + fx = new long[N]; + fy = new long[N]; + d = new long[N]; + l = new int[N]; + r = new int[N]; + arg = new int[N]; + trace = new int[N]; + } + + Hungrarian(int n1, int n2) { + this(); + n = Math.max(n1, n2); + + for (int i = 1; i <= n; ++i) { + fy[i] = l[i] = r[i] = 0; + for (int j = 1; j <= n; ++j) c[i][j] = inf; + } + } + + void addEdge(int u, int v, long cost) { + c[u][v] = Math.min(c[u][v], cost); + } + + long getC(int u, int v) { + return c[u][v] - fx[u] - fy[v]; + } + + void initBFS() { + q = new LinkedList<>(); + q.offer(start); + + for (int i = 0; i <= n; ++i) trace[i] = 0; + + for (int v = 1; v <= n; ++v) { + d[v] = getC(start, v); + arg[v] = start; + } + + finish = 0; + } + + void findAugPath() { + while (!q.isEmpty()) { + int u = q.poll(); + + for (int v = 1; v <= n; ++v) { + if (trace[v] == 0) { + long w = getC(u, v); + + if (w == 0) { + trace[v] = u; + + if (r[v] == 0) { + finish = v; + return; + } + + q.offer(r[v]); + } + + if (d[v] > w) { + d[v] = w; + arg[v] = u; + } + } + } + } + } + + void subXAddY() { + long delta = inf; + + for (int v = 1; v <= n; ++v) { + if (trace[v] == 0 && d[v] < delta) { + delta = d[v]; + } + } + + // Rotate + fx[start] += delta; + + for (int v = 1; v <= n; ++v) { + if(trace[v] != 0) { + int u = r[v]; + fy[v] -= delta; + fx[u] += delta; + } else { + d[v] -= delta; + } + } + + for (int v = 1; v <= n; ++v) if (trace[v] == 0 && d[v] == 0) { + trace[v] = arg[v]; + + if (r[v] == 0) { + finish = v; + return; + } + + q.offer(r[v]); + } + } + + void Enlarge() { + do { + int u = trace[finish]; + int nxt = l[u]; + l[u] = finish; + r[finish] = u; + finish = nxt; + } while (finish != 0); + } + + long maximumMatching() { + for (int u = 1; u <= n; ++u) { + fx[u] = c[u][1]; + for (int v = 1; v <= n; ++v) { + fx[u] = Math.min(fx[u], c[u][v]); + } + } + + for (int v = 1; v <= n; ++v) { + fy[v] = c[1][v] - fx[1]; + for (int u = 1; u <= n; ++u) { + fy[v] = Math.min(fy[v], c[u][v] - fx[u]); + } + } + + for (int u = 1; u <= n; ++u) { + start = u; + initBFS(); + while (finish == 0) { + findAugPath(); + if (finish == 0) subXAddY(); + } + Enlarge(); + } + + long ans = 0; + for (int i = 1; i <= n; ++i) { + if (c[i][l[i]] != inf) ans += c[i][l[i]]; + else l[i] = 0; + } + return ans; + } +}