4.5.1 Maximum Flow (Ford-Fulkerson)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

/*

Given a flow network with integer capacities, find the maximum flow from a given source node to a
given sink node. The flow along each edge may not exceed its capacity, and flow is conserved at
every node other than the source and sink.

The Ford-Fulkerson method repeatedly finds any augmenting path from source to sink in the residual
graph (here by depth-first search) and pushes as much flow as the path allows, until no augmenting
path remains. This algorithm should only be used on graphs with integer capacities, because some
real-valued inputs can make it fail to terminate. The Edmonds-Karp algorithm instead uses
breadth-first search and always terminates in polynomial time.

- `max_flow(source, sink)` modifies the global residual capacity matrix `cap` and returns the
  maximum flow. Nodes are numbered $[0, n)$, where $n$ is `cap.size()`.

Time Complexity:
- O(n^2*f) per call, where $n$ is the number of nodes and $f$ is the maximum flow.

Space Complexity:
- O(n^2) for storage of the flow network, where $n$ is the number of nodes.
- O(n) auxiliary stack space.

*/

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <vector>

const int64_t INF = INT64_MAX / 4;
std::vector<std::vector<int64_t>> cap;
std::vector<char> visit;

int64_t dfs(int u, int64_t f, int sink) {
  int n = static_cast<int>(cap.size());
  if (u == sink) {
    return f;
  }
  visit[u] = true;
  for (int v = 0; v < n; v++) {
    if (!visit[v] && cap[u][v] > 0) {
      int64_t flow = dfs(v, std::min(f, cap[u][v]), sink);
      if (flow > 0) {
        cap[u][v] -= flow;
        cap[v][u] += flow;
        return flow;
      }
    }
  }
  return 0;
}

int64_t max_flow(int source, int sink) {
  assert(source != sink);
  int n = static_cast<int>(cap.size());
  int64_t total = 0;
  while (true) {
    visit.assign(n, false);
    int64_t flow = dfs(source, INF, sink);
    if (flow == 0) {
      break;
    }
    total += flow;
  }
  return total;
}

/*** Example Usage ***/

#include <cassert>
using namespace std;

int main() {
  // Example graph after max flow, with each edge labeled flow/capacity:
  //            2/2
  //       1 --------> 3
  //      / \          |
  // 3/4 /   \ 1/1     | 2/2
  //    /     v        v
  //   0       4 ----> 5
  //    \     ^   3/3
  // 2/3 \   / 2/2
  //      v /
  //       2
  int nodes = 6;
  cap.assign(nodes, vector<int64_t>(nodes));
  cap[0][1] = 4;
  cap[0][2] = 3;
  cap[1][3] = 2;
  cap[1][4] = 1;
  cap[2][4] = 2;
  cap[3][5] = 2;
  cap[4][5] = 3;
  assert(max_flow(0, 5) == 5);
  return 0;
}