Alex's Anthology of Algorithms Common Code for Contests in Concise C++
Graphs / Flows and Cuts

4.5.5 Flow with Lower Bounds

4-Graphs/4.5.5_Flow_with_Lower_Bounds.cpp

Given a directed network where each edge has lower and upper capacity bounds, find whether the lower bounds can be satisfied, and optionally optimize an $s$-$t$ flow value. Reserve each edge's lower bound, leaving the difference between its bounds as residual capacity, then record the forced imbalance: for an edge from $u$ to $v$, its lower-bound flow must leave $u$ and enter $v$. A super-source supplies every node with net demand, and every node with net surplus sends that surplus to a super-sink. The bounds are feasible exactly when all such auxiliary edges can be saturated.

For $s$-$t$ flow, add infinite auxiliary edges in both directions between $s$ and $t$ before checking feasibility. The difference between their flows is one feasible signed value. After the feasibility check, augmenting from $s$ to $t$ maximizes the value; augmenting from $t$ to $s$ minimizes it.

  • BoundedFlow(n) constructs a directed lower-bound flow network with nodes in $[0, {\htmlClass{math-inline-code}{\texttt{n}}})$.
  • add_edge(u, v, lo, hi) adds an edge with lower capacity lo and upper capacity hi.
  • feasible_circulation() returns whether all edge bounds can be satisfied with flow conserved at every node.
  • max_flow(source, sink) returns the maximum feasible flow from source to sink, or std::nullopt if no feasible flow exists.
  • min_flow(source, sink) returns the minimum feasible flow from source to sink, or std::nullopt if no feasible flow exists.
  • edge_flows() returns one feasible flow value for each original edge after a successful call.

The flow value is the net flow leaving source and may be negative when the bounds force flow in the opposite direction. All capacities should be nonnegative integers and lo $\leq$ hi. Choose INF larger than the absolute value of any possible finite flow. Each feasibility or optimization call starts from the original network, so different variants may be solved successively on the same instance.

Time Complexity
$O(n^{2} \cdot (n + m))$ per feasibility check or optimized flow call, where $n$ is the number of nodes and $m$ is the number of edges.
Space Complexity
$O(\max\left(n, m\right))$ for network storage and auxiliary arrays.

Implementation

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

class BoundedFlow {
  struct Edge {
    int v, rev;
    int64_t cap;
  };

  struct OriginalEdge {
    int u, v, index;
    int64_t lo, hi;
  };

  static constexpr int64_t INF = (1LL << 60);
  int nodes;
  std::vector<std::vector<Edge>> adj;
  std::vector<OriginalEdge> original;
  std::vector<int> level, ptr;

  void add_residual_edge(int u, int v, int64_t cap) {
    adj[u].push_back(Edge{v, static_cast<int>(adj[v].size()), cap});
    adj[v].push_back(Edge{u, static_cast<int>(adj[u].size()) - 1, 0});
  }

  bool bfs(int source, int sink) {
    level.assign(static_cast<int>(adj.size()), -1);
    level[source] = 0;
    std::queue<int> q;
    q.push(source);
    while (!q.empty()) {
      int u = q.front();
      q.pop();
      for (const Edge &e : adj[u]) {
        if (e.cap > 0 && level[e.v] == -1) {
          level[e.v] = level[u] + 1;
          q.push(e.v);
        }
      }
    }
    return level[sink] != -1;
  }

  int64_t dfs(int u, int sink, int64_t pushed) {
    if (u == sink || pushed == 0) {
      return pushed;
    }
    for (; ptr[u] < static_cast<int>(adj[u].size()); ptr[u]++) {
      Edge &e = adj[u][ptr[u]];
      if (e.cap > 0 && level[e.v] == level[u] + 1) {
        int64_t next = dfs(e.v, sink, std::min(pushed, e.cap));
        if (next > 0) {
          e.cap -= next;
          adj[e.v][e.rev].cap += next;
          return next;
        }
      }
    }
    return 0;
  }

  int64_t dinic(int source, int sink) {
    int64_t flow = 0;
    while (bfs(source, sink)) {
      ptr.assign(static_cast<int>(adj.size()), 0);
      while (int64_t pushed = dfs(source, sink, INF)) {
        flow += pushed;
      }
    }
    return flow;
  }

  bool satisfy_demands(
      int source, int sink, int &backward_index, int &forward_index, int64_t &forced_flow
  ) {
    int super_source = static_cast<int>(adj.size());
    int super_sink = super_source + 1;
    adj.resize(super_sink + 1);
    level.resize(super_sink + 1);
    ptr.resize(super_sink + 1);
    std::vector<int64_t> balance(nodes);
    for (const OriginalEdge &e : original) {
      balance[e.u] -= e.lo;
      balance[e.v] += e.lo;
    }
    if (source != -1) {
      backward_index = static_cast<int>(adj[sink].size());
      add_residual_edge(sink, source, INF);
      forward_index = static_cast<int>(adj[source].size());
      add_residual_edge(source, sink, INF);
    }
    int64_t need = 0;
    for (int u = 0; u < nodes; u++) {
      if (balance[u] > 0) {
        add_residual_edge(super_source, u, balance[u]);
        need += balance[u];
      } else if (balance[u] < 0) {
        add_residual_edge(u, super_sink, -balance[u]);
      }
    }
    bool ok = dinic(super_source, super_sink) == need;
    forced_flow = 0;
    if (source != -1) {
      forced_flow = adj[source][adj[sink][backward_index].rev].cap -
                    adj[sink][adj[source][forward_index].rev].cap;
    }
    return ok;
  }

  std::optional<int64_t> optimize_flow(int source, int sink, bool maximize) {
    int backward_index = -1, forward_index = -1;
    int64_t value = 0;
    if (!satisfy_demands(source, sink, backward_index, forward_index, value)) {
      return std::nullopt;
    }
    Edge &backward = adj[sink][backward_index];
    Edge &backward_reverse = adj[source][backward.rev];
    Edge &forward = adj[source][forward_index];
    Edge &forward_reverse = adj[sink][forward.rev];
    backward.cap = backward_reverse.cap = 0;
    forward.cap = forward_reverse.cap = 0;
    return maximize ? value + dinic(source, sink) : value - dinic(sink, source);
  }

  void reset() {
    adj.assign(nodes, {});
    level.assign(nodes, 0);
    ptr.assign(nodes, 0);
    for (OriginalEdge &e : original) {
      e.index = static_cast<int>(adj[e.u].size());
      add_residual_edge(e.u, e.v, e.hi - e.lo);
    }
  }

 public:
  explicit BoundedFlow(int n) : nodes(n), adj(n), level(n), ptr(n) {}

  void add_edge(int u, int v, int64_t lo, int64_t hi) {
    assert(0 <= u && u < nodes && 0 <= v && v < nodes);
    assert(0 <= lo && lo <= hi);
    original.push_back(OriginalEdge{u, v, static_cast<int>(adj[u].size()), lo, hi});
    add_residual_edge(u, v, hi - lo);
  }

  bool feasible_circulation() {
    reset();
    int backward_index = -1, forward_index = -1;
    int64_t value = 0;
    return satisfy_demands(-1, -1, backward_index, forward_index, value);
  }

  std::optional<int64_t> max_flow(int source, int sink) {
    assert(0 <= source && source < nodes && 0 <= sink && sink < nodes && source != sink);
    reset();
    return optimize_flow(source, sink, true);
  }

  std::optional<int64_t> min_flow(int source, int sink) {
    assert(0 <= source && source < nodes && 0 <= sink && sink < nodes && source != sink);
    reset();
    return optimize_flow(source, sink, false);
  }

  std::vector<int64_t> edge_flows() const {
    std::vector<int64_t> flow;
    flow.reserve(original.size());
    for (const OriginalEdge &e : original) {
      const Edge &forward = adj[e.u][e.index];
      flow.push_back(e.hi - forward.cap);
    }
    return flow;
  }
};

Example Usage

using namespace std;

int main() {
  //           [1,4]       [1,3]
  //       0 --------> 1 --------> 3
  //       |                       ^
  // [0,2] |                       | [0,2]
  //       v                       |
  //       2 ----------------------+
  BoundedFlow g(4);
  g.add_edge(0, 1, 1, 4);
  g.add_edge(0, 2, 0, 2);
  g.add_edge(1, 3, 1, 3);
  g.add_edge(2, 3, 0, 2);
  assert(g.max_flow(0, 3) == 5);
  vector<int64_t> max_flow = g.edge_flows();
  assert(max_flow[0] + max_flow[1] == 5);

  assert(g.min_flow(0, 3) == 1);

  BoundedFlow rev(2);
  rev.add_edge(1, 0, 1, 1);
  assert(rev.max_flow(0, 1) == -1);
  assert(rev.min_flow(0, 1) == -1);
  assert((rev.edge_flows() == vector<int64_t>{1}));

  //            [2,3]
  //       0 ----------> 1
  //       ^             |
  // [2,3] |             | [2,3]
  //       |             |
  //       2 <-----------+
  BoundedFlow circulation(3);
  circulation.add_edge(0, 1, 2, 3);
  circulation.add_edge(1, 2, 2, 3);
  circulation.add_edge(2, 0, 2, 3);
  assert(circulation.feasible_circulation());
  assert((circulation.edge_flows() == vector<int64_t>{2, 2, 2}));
  return 0;
}
Back to top