andres::graph
max-flow.hxx
1 #pragma once
2 #ifndef ANDRES_GRAPH_MAX_FLOW_HXX
3 #define ANDRES_GRAPH_MAX_FLOW_HXX
4 
5 #include <cassert>
6 #include <cstddef>
7 #include <vector>
8 #include <queue>
9 #include <limits> // std::numeric_limit
10 #include <algorithm> // std::max, std::min
11 
12 #include "andres/graph/digraph.hxx"
13 #include "andres/graph/shortest-paths.hxx"
14 
15 namespace andres {
16 namespace graph {
17 
27 template<class GRAPH, class FLOW>
29 public:
30  typedef GRAPH GraphType;
31  typedef FLOW Flow;
32  typedef typename GraphType::VertexIterator VertexIterator;
33  typedef typename GraphType::EdgeIterator EdgeIterator;
34  typedef typename GraphType::AdjacencyIterator AdjacencyIterator;
35 
37  void clear();
38  Flow maxFlow() const;
39  Flow flow(const std::size_t) const;
40  std::size_t numberOfPushes() const;
41  std::size_t numberOfRelabels() const;
42  template<class EDGE_WEIGHT_ITERATOR>
43  MaxFlowPushRelabel(const GraphType&, EDGE_WEIGHT_ITERATOR, const std::size_t, const std::size_t);
44  template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
45  MaxFlowPushRelabel(const GraphType&, const SUBGRAPH_MASK&, EDGE_WEIGHT_ITERATOR, const std::size_t, const std::size_t);
46  template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
47  Flow operator()(const GraphType&, const SUBGRAPH_MASK&, EDGE_WEIGHT_ITERATOR, const std::size_t, const std::size_t);
48 
49 private:
50  template<class EDGE_WEIGHT_ITERATOR>
51  void push(const GraphType&, EDGE_WEIGHT_ITERATOR, const std::size_t);
52  template<class EDGE_WEIGHT_ITERATOR>
53  void pushBack(const GraphType&, EDGE_WEIGHT_ITERATOR, const std::size_t);
54  template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
55  void relabel(const GraphType&, const SUBGRAPH_MASK&, EDGE_WEIGHT_ITERATOR, const std::size_t);
56  template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
57  void discharge(const GraphType&, const SUBGRAPH_MASK&, EDGE_WEIGHT_ITERATOR, const std::size_t);
58  void gapRelabel(const GraphType&, const std::size_t);
59 
60  std::vector<std::size_t> height_;
61  std::vector<std::size_t> labelCount_;
62  std::vector<Flow> excess_;
63  std::vector<Flow> flow_;
64  std::vector<bool> active_;
65  std::queue<std::size_t> queue_;
66  std::size_t sourceVertexIndex_;
67  std::size_t sinkVertexIndex_;
68  std::size_t pushCount_;
69  std::size_t relabelCount_;
70 };
71 
74 template<class GRAPH, class FLOW>
75 inline
77 : height_(),
78  labelCount_(),
79  excess_(),
80  flow_(),
81  active_(),
82  queue_(),
83  sourceVertexIndex_(),
84  sinkVertexIndex_(),
85  pushCount_(),
86  relabelCount_()
87 {}
88 
96 template<class GRAPH, class FLOW>
97 template<class EDGE_WEIGHT_ITERATOR>
98 inline
100  const GraphType& graph,
101  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
102  const std::size_t sourceVertexIndex,
103  const std::size_t sinkVertexIndex
104 )
105 : height_(),
106  labelCount_(),
107  excess_(),
108  flow_(),
109  active_(),
110  queue_(),
111  sourceVertexIndex_(),
112  sinkVertexIndex_(),
113  pushCount_(),
114  relabelCount_()
115 {
116  (*this)(graph, DefaultSubgraphMask<>(), edgeWeightIterator, sourceVertexIndex, sinkVertexIndex);
117 }
118 
127 template<class GRAPH, class FLOW>
128 template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
129 inline
131  const GraphType& graph,
132  const SUBGRAPH_MASK& mask,
133  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
134  const std::size_t sourceVertexIndex,
135  const std::size_t sinkVertexIndex
136 )
137 : height_(),
138  labelCount_(),
139  excess_(),
140  flow_(),
141  active_(),
142  queue_(),
143  sourceVertexIndex_(),
144  sinkVertexIndex_(),
145  pushCount_(),
146  relabelCount_()
147 {
148  (*this)(graph, mask, edgeWeightIterator, sourceVertexIndex, sinkVertexIndex);
149 }
150 
153 template<class GRAPH, class FLOW>
154 inline void
156  height_.clear();
157  labelCount_.clear();
158  excess_.clear();
159  flow_.clear();
160  active_.clear();
161  sourceVertexIndex_ = 0;
162  sinkVertexIndex_ = 0;
163  pushCount_ = 0;
164  relabelCount_ = 0;
165 
166  assert(queue_.empty());
167 }
168 
173 template<class GRAPH, class FLOW>
176  return excess_[sinkVertexIndex_];
177 }
178 
184 template<class GRAPH, class FLOW>
187  const std::size_t edgeIndex
188 ) const {
189  assert(edgeIndex < flow_.size());
190  return flow_[edgeIndex];
191 }
192 
197 template<class GRAPH, class FLOW>
198 inline std::size_t
200  return pushCount_;
201 }
202 
207 template<class GRAPH, class FLOW>
208 inline std::size_t
210  return relabelCount_;
211 }
212 
221 template<class GRAPH, class FLOW>
222 template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
225  const GraphType& graph,
226  const SUBGRAPH_MASK& mask,
227  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
228  const std::size_t sourceVertexIndex,
229  const std::size_t sinkVertexIndex
230 ) {
231  assert(mask.vertex(sourceVertexIndex) && mask.vertex(sinkVertexIndex));
232 
233  const std::size_t numberOfVertices = graph.numberOfVertices();
234  const std::size_t numberOfEdges = graph.numberOfEdges();
235 
236  // initialization
237  sourceVertexIndex_ = sourceVertexIndex;
238  sinkVertexIndex_ = sinkVertexIndex;
239  pushCount_ = 0;
240  relabelCount_ = 0;
241  height_.resize(numberOfVertices);
242  excess_.resize(numberOfVertices);
243  active_.resize(numberOfVertices);
244  std::fill(height_.begin(), height_.end(), std::size_t());
245  std::fill(excess_.begin(), excess_.end(), Flow());
246  std::fill(active_.begin(), active_.end(), false);
247  height_[sourceVertexIndex] = numberOfVertices;
248  excess_[sourceVertexIndex] = std::numeric_limits<Flow>::max(); // this is supposed to be infinite
249  active_[sourceVertexIndex] = true;
250  active_[sinkVertexIndex] = true;
251  labelCount_.resize((2 * numberOfVertices) + 2); // 2n + 1 is the maximum possible height for a vertex
252  std::fill(labelCount_.begin(), labelCount_.end(), std::size_t());
253  labelCount_[0] = numberOfVertices - 1;
254  labelCount_[numberOfVertices] = 1;
255  flow_.resize(numberOfEdges);
256  std::fill(flow_.begin(), flow_.end(), Flow());
257 
258  // first, push as much flow as possible from the source to all adjacent vertices
259  for(EdgeIterator it = graph.edgesFromVertexBegin(sourceVertexIndex); it != graph.edgesFromVertexEnd(sourceVertexIndex); ++it) {
260  const std::size_t edgeIndex = *it;
261  if (mask.edge(edgeIndex)) {
262  const std::size_t v = graph.vertexOfEdge(edgeIndex, 1);
263  if (mask.vertex(v)) {
264  push(graph, edgeWeightIterator, edgeIndex);
265  }
266  }
267  }
268 
269  while(!queue_.empty()) { // main loop
270  const std::size_t u = queue_.front();
271  queue_.pop();
272  active_[u] = false;
273  discharge(graph, mask, edgeWeightIterator, u);
274  active_[u] = false; // TODO: why does putting active_[u] = false twice decrease the number of pushes?
275  }
276 
277  return maxFlow();
278 }
279 
286 template<class GRAPH, class FLOW>
287 template<class EDGE_WEIGHT_ITERATOR>
288 inline void
290  const GRAPH& graph,
291  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
292  const std::size_t edgeIndex
293 ) {
294  const std::size_t u = graph.vertexOfEdge(edgeIndex, 0);
295  const std::size_t v = graph.vertexOfEdge(edgeIndex, 1);
296  // push from u to v
297  Flow pushAmount = std::min(excess_[u], edgeWeightIterator[edgeIndex] - flow_[edgeIndex]);
298  flow_[edgeIndex] += pushAmount;
299  excess_[u] -= pushAmount;
300  excess_[v] += pushAmount;
301  if (!active_[v] && excess_[v] > 0) {
302  active_[v] = true;
303  queue_.push(v);
304  }
305  pushCount_++;
306 }
307 
314 template<class GRAPH, class FLOW>
315 template<class EDGE_WEIGHT_ITERATOR>
316 inline void
317 MaxFlowPushRelabel<GRAPH, FLOW>::pushBack(
318  const GRAPH& graph,
319  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
320  const std::size_t edgeIndex
321 ) {
322  const std::size_t u = graph.vertexOfEdge(edgeIndex, 0);
323  const std::size_t v = graph.vertexOfEdge(edgeIndex, 1);
324  // push back from v to u
325  Flow pushAmount = std::min(excess_[v], flow_[edgeIndex]);
326  flow_[edgeIndex] -= pushAmount;
327  excess_[u] += pushAmount;
328  excess_[v] -= pushAmount;
329  if (!active_[u] && excess_[u] > 0) {
330  active_[u] = true;
331  queue_.push(u);
332  }
333  pushCount_++;
334 }
335 
343 template<class GRAPH, class FLOW>
344 template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
345 inline void
346 MaxFlowPushRelabel<GRAPH, FLOW>::relabel(
347  const GRAPH& graph,
348  const SUBGRAPH_MASK& mask,
349  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
350  const std::size_t u
351 ) {
352  std::size_t minHeight = 2 * graph.numberOfVertices();
353  const std::size_t oldHeight = height_[u];
354  for(EdgeIterator it = graph.edgesFromVertexBegin(u); it != graph.edgesFromVertexEnd(u); ++it) {
355  const std::size_t edgeIndex = *it;
356  if (mask.edge(edgeIndex)) {
357  const std::size_t v = graph.vertexOfEdge(edgeIndex, 1); // edge is (u, v)
358  if (mask.vertex(v)) {
359  if(edgeWeightIterator[edgeIndex] - flow_[edgeIndex] > 0) {
360  minHeight = std::min(minHeight, height_[v]);
361  }
362  }
363  }
364  }
365  for(EdgeIterator it = graph.edgesToVertexBegin(u); it != graph.edgesToVertexEnd(u); ++it) {
366  const std::size_t edgeIndex = *it;
367  if (mask.edge(edgeIndex)) {
368  const std::size_t v = graph.vertexOfEdge(edgeIndex, 0); // edge is (v, u)
369  if (mask.vertex(v)) {
370  if(flow_[edgeIndex] > 0) {
371  minHeight = std::min(minHeight, height_[v]);
372  }
373  }
374  }
375  }
376  height_[u] = minHeight + 1;
377  if (!active_[u] && excess_[u] > 0) {
378  active_[u] = true;
379  queue_.push(u);
380  }
381  relabelCount_++;
382 
383  // gap relabeling
384  labelCount_[oldHeight]--;
385  labelCount_[height_[u]]++;
386  if (labelCount_[oldHeight] == 0) {
387  gapRelabel(graph, oldHeight);
388  }
389 }
390 
398 template<class GRAPH, class FLOW>
399 template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
400 inline void
401 MaxFlowPushRelabel<GRAPH, FLOW>::discharge(
402  const GRAPH& graph,
403  const SUBGRAPH_MASK& mask,
404  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
405  const std::size_t u
406 ) {
407  while(excess_[u] > 0) {
408  for(EdgeIterator it = graph.edgesFromVertexBegin(u); it != graph.edgesFromVertexEnd(u); ++it) {
409  const std::size_t edgeIndex = *it;
410  if (mask.edge(edgeIndex)) {
411  const std::size_t v = graph.vertexOfEdge(edgeIndex, 1); // edge is the pair (u, v)
412  if (mask.vertex(v)) {
413  if(edgeWeightIterator[edgeIndex] - flow_[edgeIndex] > 0 && height_[u] > height_[v]) {
414  push(graph, edgeWeightIterator, edgeIndex);
415  }
416  }
417  }
418  }
419  for(EdgeIterator it = graph.edgesToVertexBegin(u); it != graph.edgesToVertexEnd(u); ++it) {
420  const std::size_t edgeIndex = *it;
421  if (mask.edge(edgeIndex)) {
422  const std::size_t v = graph.vertexOfEdge(edgeIndex, 0); // edge is the pair (v, u)
423  if (mask.vertex(v)) {
424  if(flow_[edgeIndex] > 0 && height_[u] > height_[v]) {
425  pushBack(graph, edgeWeightIterator, edgeIndex);
426  }
427  }
428  }
429  }
430  relabel(graph, mask, edgeWeightIterator, u);
431  }
432 }
433 
439 template<class GRAPH, class FLOW>
440 inline void
441 MaxFlowPushRelabel<GRAPH, FLOW>::gapRelabel(
442  const GRAPH& graph,
443  const std::size_t threshold
444 ) {
445  for(std::size_t i = 0; i < graph.numberOfVertices(); i++) {
446  if(height_[i] > threshold) {
447  height_[i] = std::max(height_[i], graph.numberOfVertices() + 1);;
448  }
449  }
450 }
451 
459 template<class GRAPH, class FLOW>
461 public:
462  typedef GRAPH GraphType;
463  typedef FLOW Flow;
464 
466  template <class EDGE_WEIGHT_ITERATOR>
467  MaxFlowEdmondsKarp(const GraphType&, EDGE_WEIGHT_ITERATOR, const std::size_t, const std::size_t);
468  template <class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
469  MaxFlowEdmondsKarp(const GraphType&, const SUBGRAPH_MASK&, EDGE_WEIGHT_ITERATOR, const std::size_t, const std::size_t);
470  void clear();
471  Flow maxFlow() const;
472  Flow flow(const std::size_t) const;
473  template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
474  Flow operator()(const GraphType&, const SUBGRAPH_MASK&, const EDGE_WEIGHT_ITERATOR, const std::size_t, const std::size_t);
475 
476 private:
477  std::deque<std::size_t> augmentingPath_;
478  Digraph<> rgraph_;
479  std::vector<Flow> flow_;
480  std::size_t sourceVertexIndex_;
481  std::size_t sinkVertexIndex_;
482  Flow maxFlow_;
483 
484  template <class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
485  class ResidualMask {
486  public:
487  ResidualMask(
488  const std::size_t numberOfEdges,
489  const std::vector<Flow>& flow,
490  const EDGE_WEIGHT_ITERATOR& edgeWeightIterator,
491  const SUBGRAPH_MASK& subgraphMask
492  )
493  : numberOfEdges_(numberOfEdges),
494  flow_(flow),
495  edgeWeightIterator_(edgeWeightIterator),
496  subgraphMask_(subgraphMask)
497  {}
498  bool vertex(const std::size_t v) const
499  { return subgraphMask_.vertex(v); }
500  bool edge(const std::size_t e) const
501  { return capacity(e) > Flow() && inMask(e); }
502  Flow capacity(const std::size_t e) const
503  {
504  if(e >= numberOfEdges_) {
505  return flow_[e - numberOfEdges_];
506  }
507  else {
508  return edgeWeightIterator_[e] - flow_[e];
509  }
510  }
511 
512  private:
513  const std::size_t numberOfEdges_;
514  const std::vector<Flow>& flow_;
515  const EDGE_WEIGHT_ITERATOR& edgeWeightIterator_;
516  const SUBGRAPH_MASK& subgraphMask_;
517  bool inMask(const std::size_t e) const
518  {
519  if(e >= numberOfEdges_) {
520  return subgraphMask_.edge(e - numberOfEdges_);
521  }
522  else {
523  return subgraphMask_.edge(e);
524  }
525  }
526  };
527 };
528 
531 template<class GRAPH, class FLOW>
532 inline
534 : augmentingPath_(),
535  rgraph_(),
536  flow_(),
537  sourceVertexIndex_(),
538  sinkVertexIndex_(),
539  maxFlow_()
540 {}
541 
549 template<class GRAPH, class FLOW>
550 template<class EDGE_WEIGHT_ITERATOR>
551 inline
553  const GraphType& graph,
554  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
555  const std::size_t sourceVertexIndex,
556  const std::size_t sinkVertexIndex
557 )
558 : augmentingPath_(),
559  rgraph_(),
560  flow_(),
561  sourceVertexIndex_(),
562  sinkVertexIndex_(),
563  maxFlow_()
564 {
565  (*this)(graph, DefaultSubgraphMask<>(), edgeWeightIterator, sourceVertexIndex, sinkVertexIndex);
566 }
567 
576 template<class GRAPH, class FLOW>
577 template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
578 inline
580  const GraphType& graph,
581  const SUBGRAPH_MASK& subgraphMask,
582  EDGE_WEIGHT_ITERATOR edgeWeightIterator,
583  const std::size_t sourceVertexIndex,
584  const std::size_t sinkVertexIndex
585 )
586 : augmentingPath_(),
587  rgraph_(),
588  flow_(),
589  sourceVertexIndex_(),
590  sinkVertexIndex_(),
591  maxFlow_()
592 {
593  (*this)(graph, subgraphMask, edgeWeightIterator, sourceVertexIndex, sinkVertexIndex);
594 }
595 
598 template<class GRAPH, class FLOW>
599 inline void
601  augmentingPath_.clear();
602  rgraph_.assign();
603  flow_.clear();
604  sourceVertexIndex_ = 0;
605  sinkVertexIndex_ = 0;
606  maxFlow_ = Flow();
607 }
608 
613 template<class GRAPH, class FLOW>
616  return maxFlow_;
617 }
618 
624 template<class GRAPH, class FLOW>
627  const std::size_t edgeIndex
628 ) const {
629  assert(edgeIndex < flow_.size());
630  return flow_[edgeIndex];
631 }
632 
641 template<class GRAPH, class FLOW>
642 template<class EDGE_WEIGHT_ITERATOR, class SUBGRAPH_MASK>
645  const GraphType& graph,
646  const SUBGRAPH_MASK& subgraphMask,
647  const EDGE_WEIGHT_ITERATOR edgeWeightIterator,
648  const std::size_t sourceVertexIndex,
649  const std::size_t sinkVertexIndex
650 ) {
651  const std::size_t numberOfEdges = graph.numberOfEdges();
652 
653  // initialization
654  sourceVertexIndex_ = sourceVertexIndex;
655  sinkVertexIndex_ = sinkVertexIndex;
656  flow_.resize(numberOfEdges);
657  std::fill(flow_.begin(), flow_.end(), Flow());
658  augmentingPath_.clear();
659  rgraph_.assign(graph.numberOfVertices());
660  for(std::size_t edge = 0; edge < numberOfEdges; ++edge) {
661  rgraph_.insertEdge(graph.vertexOfEdge(edge, 0), graph.vertexOfEdge(edge, 1));
662  }
663  for(std::size_t edge = 0; edge < numberOfEdges; ++edge) {
664  rgraph_.insertEdge(graph.vertexOfEdge(edge, 1), graph.vertexOfEdge(edge, 0));
665  }
666  ResidualMask<EDGE_WEIGHT_ITERATOR, SUBGRAPH_MASK> residualMask(numberOfEdges, flow_, edgeWeightIterator, subgraphMask);
667  maxFlow_ = Flow();
668 
669  // while there is an augmenting path, augment flow along it.
670  // choose the shortest augmenting path by number of edges
671  Flow minResidualCapacity;
672  while(spspEdges(rgraph_, residualMask, sourceVertexIndex_, sinkVertexIndex_, augmentingPath_)) {
673  minResidualCapacity = residualMask.capacity(augmentingPath_[0]);
674  for(std::deque<std::size_t>::iterator it = augmentingPath_.begin(); it < augmentingPath_.end(); ++it) {
675  if(residualMask.capacity(*it) < minResidualCapacity) {
676  minResidualCapacity = residualMask.capacity(*it);
677  }
678  }
679  for(std::deque<std::size_t>::iterator it = augmentingPath_.begin(); it < augmentingPath_.end(); ++it) {
680  if(*it < numberOfEdges) {
681  flow_[*it] += minResidualCapacity; // this call updates also the residualMask which holds a reference to flow_
682  }
683  else {
684  flow_[*it - numberOfEdges] -= minResidualCapacity; // this call updates also the residualMask which holds a reference to flow_
685  }
686  }
687  }
688 
689  // sum flow out of source to get the max flow
690  for(std::size_t edge = 0; edge < numberOfEdges; ++edge) {
691  if(graph.vertexOfEdge(edge, 0) == sourceVertexIndex_) {
692  maxFlow_ += flow_[edge];
693  }
694  }
695 
696  return maxFlow_;
697 }
698 
699 } // namespace graph
700 } // namespace andres
701 
702 #endif // #ifndef ANDRES_GRAPH_MAX_FLOW_HXX