Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
generic_max_flow.h
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
14// An implementation of a push-relabel algorithm for the max flow problem.
15//
16// In the following, we consider a graph G = (V,E,s,t) where V denotes the set
17// of nodes (vertices) in the graph, E denotes the set of arcs (edges). s and t
18// denote distinguished nodes in G called source and target. n = |V| denotes the
19// number of nodes in the graph, and m = |E| denotes the number of arcs in the
20// graph.
21//
22// Each arc (v,w) is associated a capacity c(v,w).
23//
24// A flow is a function from E to R such that:
25//
26// a) f(v,w) <= c(v,w) for all (v,w) in E (capacity constraint.)
27//
28// b) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint.)
29//
30// c) sum on v f(v,w) = 0 (flow conservation.)
31//
32// The goal of this algorithm is to find the maximum flow from s to t, i.e.
33// for example to maximize sum v f(s,v).
34//
35// The starting reference for this class of algorithms is:
36// A.V. Goldberg and R.E. Tarjan. A new approach to the maximum flow problem.
37// ACM Symposium on Theory of Computing, pp. 136-146.
38// http://portal.acm.org/citation.cfm?id=12144.
39//
40// The basic idea of the algorithm is to handle preflows instead of flows,
41// and to refine preflows until a maximum flow is obtained.
42// A preflow is like a flow, except that the inflow can be larger than the
43// outflow. If it is the case at a given node v, it is said that there is an
44// excess at node v, and inflow = outflow + excess.
45//
46// More formally, a preflow is a function f such that:
47//
48// 1) f(v,w) <= c(v,w) for all (v,w) in E (capacity constraint). c(v,w) is a
49// value representing the maximum capacity for arc (v,w).
50//
51// 2) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint)
52//
53// 3) excess(v) = sum on u f(u,v) >= 0 is the excess at node v, the
54// algebraic sum of all the incoming preflows at this node.
55//
56// Each node has an associated "height", in addition to its excess. The
57// height of the source is defined to be equal to n, and cannot change. The
58// height of the target is defined to be zero, and cannot change either. The
59// height of all the other nodes is initialized at zero and is updated during
60// the algorithm (see below). For those who want to know the details, the height
61// of a node, corresponds to a reduced cost, and this enables one to prove that
62// the algorithm actually computes the max flow. Note that the height of a node
63// can be initialized to the distance to the target node in terms of number of
64// nodes. This has not been tried in this implementation.
65//
66// A node v is said to be *active* if excess(v) > 0.
67//
68// In this case the following operations can be applied to it:
69//
70// - if there are *admissible* incident arcs, i.e. arcs which are not saturated,
71// and whose head's height is lower than the height of the active node
72// considered, a PushFlow operation can be applied. It consists in sending as
73// much flow as both the excess at the node and the capacity of the arc
74// permit.
75// - if there are no admissible arcs, the active node considered is relabeled,
76// i.e. its height is increased to 1 + the minimum height of its neighboring
77// nodes on admissible arcs.
78// This is implemented in Discharge, which itself calls PushFlow and Relabel.
79//
80// Before running Discharge, it is necessary to initialize the algorithm with a
81// preflow. This is done in InitializePreflow, which saturates all the arcs
82// leaving the source node, and sets the excess at the heads of those arcs
83// accordingly.
84//
85// The algorithm terminates when there are no remaining active nodes, i.e. all
86// the excesses at all nodes are equal to zero. In this case, a maximum flow is
87// obtained.
88//
89// The complexity of this algorithm depends amongst other things on the choice
90// of the next active node. It has been shown, for example in:
91// L. Tuncel, "On the Complexity of Preflow-Push Algorithms for Maximum-Flow
92// Problems", Algorithmica 11(4): 353-359 (1994).
93// and
94// J. Cheriyan and K. Mehlhorn, "An analysis of the highest-level selection rule
95// in the preflow-push max-flow algorithm", Information processing letters,
96// 69(5):239-242 (1999).
97// http://www.math.uwaterloo.ca/~jcheriya/PS_files/me3.0.ps
98//
99// ...that choosing the active node with the highest level yields a
100// complexity of O(n^2 * sqrt(m)).
101//
102// TODO(user): implement the above active node choice rule.
103//
104// This has been validated experimentally in:
105// R.K. Ahuja, M. Kodialam, A.K. Mishra, and J.B. Orlin, "Computational
106// Investigations of Maximum Flow Algorithms", EJOR 97:509-542(1997).
107// http://jorlin.scripts.mit.edu/docs/publications/58-comput%20investigations%20of.pdf.
108//
109//
110// TODO(user): an alternative would be to evaluate:
111// A.V. Goldberg, "The Partial Augment-Relabel Algorithm for the Maximum Flow
112// Problem.” In Proceedings of Algorithms ESA, LNCS 5193:466-477, Springer 2008.
113// http://www.springerlink.com/index/5535k2j1mt646338.pdf
114//
115// An interesting general reference on network flows is:
116// R. K. Ahuja, T. L. Magnanti, J. B. Orlin, "Network Flows: Theory, Algorithms,
117// and Applications," Prentice Hall, 1993, ISBN: 978-0136175490,
118// http://www.amazon.com/dp/013617549X
119//
120// Keywords: Push-relabel, max-flow, network, graph, Goldberg, Tarjan, Dinic,
121// Dinitz.
122
123#ifndef OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_
124#define OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_
125
126#include <algorithm>
127#include <cstdint>
128#include <limits>
129#include <memory>
130#include <string>
131#include <utility>
132#include <vector>
133
134#include "absl/strings/string_view.h"
135#include "ortools/base/logging.h"
136#include "ortools/graph/flow_problem.pb.h"
137#include "ortools/util/stats.h"
138#include "ortools/util/zvector.h"
139
140namespace operations_research {
141
142// Specific but efficient priority queue implementation. The priority type must
143// be an integer. The queue allows to retrieve the element with highest priority
144// but only allows pushes with a priority greater or equal to the highest
145// priority in the queue minus one. All operations are in O(1) and the memory is
146// in O(num elements in the queue). Elements with the same priority are
147// retrieved with LIFO order.
148//
149// Note(user): As far as I know, this is an original idea and is the only code
150// that use this in the Maximum Flow context. Papers usually refer to an
151// height-indexed array of simple linked lists of active node with the same
152// height. Even worse, sometimes they use double-linked list to allow arbitrary
153// height update in order to detect missing height (used for the Gap heuristic).
154// But this can actually be implemented a lot more efficiently by just
155// maintaining the height distribution of all the node in the graph.
156template <typename Element, typename IntegerPriority>
158 public:
159 PriorityQueueWithRestrictedPush() : even_queue_(), odd_queue_() {}
160
161#ifndef SWIG
162 // This type is neither copyable nor movable.
164 delete;
166 const PriorityQueueWithRestrictedPush&) = delete;
167#endif
168
169 // Is the queue empty?
170 bool IsEmpty() const;
171
172 // Clears the queue.
173 void Clear();
174
175 // Push a new element in the queue. Its priority must be greater or equal to
176 // the highest priority present in the queue, minus one. This condition is
177 // DCHECKed, and violating it yields erroneous queue behavior in NDEBUG mode.
178 void Push(Element element, IntegerPriority priority);
179
180 // Returns the element with highest priority and remove it from the queue.
181 // IsEmpty() must be false, this condition is DCHECKed.
182 Element Pop();
183
184 private:
185 // Helper function to get the last element of a vector and pop it.
186 Element PopBack(std::vector<std::pair<Element, IntegerPriority> >* queue);
187
188 // This is the heart of the algorithm. basically we split the elements by
189 // parity of their priority and the precondition on the Push() ensures that
190 // both vectors are always sorted by increasing priority.
191 std::vector<std::pair<Element, IntegerPriority> > even_queue_;
192 std::vector<std::pair<Element, IntegerPriority> > odd_queue_;
193};
194
195// We want an enum for the Status of a max flow run, and we want this
196// enum to be scoped under GenericMaxFlow<>. Unfortunately, swig
197// doesn't handle templated enums very well, so we need a base,
198// untemplated class to hold it.
200 public:
201 enum Status {
202 NOT_SOLVED, // The problem was not solved, or its data were edited.
203 OPTIMAL, // Solve() was called and found an optimal solution.
204 INT_OVERFLOW, // There is a feasible flow > max possible flow.
205
206 // TODO(user): These are no longer used. Remove.
209 };
210};
211
212// Generic MaxFlow that works on graphs with the notion of "reverse" arcs.
213// That is, the graph is directed, and each arc 'tail -> head' must be
214// associated to an unique reverse arc going in the opposite direction
215// 'head -> tail'. We must also have reverse[reverse[arc]] = arc.
216//
217// This works with all the reverse arc graphs from 'util/graph/graph.h' and
218// uses the API defined there.
219//
220// We actually support two kind of graphs with "reverse" arcs depending on the
221// value of Graph::kHasNegativeReverseArcs:
222// - If it is true, the arcs from the "user graph" are called "direct" arcs and
223// are assumed to be indexed in [0, num_arcs). These are the only ones that
224// can have a positive capacity. All the "reverse/opposite" of these arcs
225// will always have negative indices in [-num_arcs, 0) and a capacity of
226// zero.
227// - If it is false, then we only have "direct" arcs in [0, num_arcs), the
228// reverse of each arc must lives in the same space, and both an arc and its
229// reverse can have a positive initial capacity. This can lead to twice fewer
230// arcs and a faster algo if the "user graph" as a lot of reverse arcs
231// already.
232//
233// Flow types and overflow: To optimize memory usage and speed we have two
234// integer types representing flow quantities.
235// - ArcFlowType (can be unsigned) is used to represent the flow per arc. The
236// type should be big enough to hold the initial capacity of an arc + its
237// reverse. So if all your initial capacity <=
238// std::numeric_limit<ArcFlowType>::max() / 2, you are safe.
239// - FlowSumType (must be signed) is used to hold the sum of a flow going
240// through a node or through the full network. Ideally, the sum of all the
241// capacity of the arcs leaving the source should fit in a FlowSumType. If
242// this is the case, you will get the better performance, and always an
243// OPTIMAL result. Otherwise, the code might be slower, and you will get
244// OPTIMAL only if the max_flow that can be sent through the network fit a
245// FlowSumType, you will get the INT_OVERFLOW status if more flow than that
246// can be sent.
247template <typename Graph, typename ArcFlowT = int64_t,
248 typename FlowSumT = int64_t>
250 public:
251 typedef typename Graph::NodeIndex NodeIndex;
252 typedef typename Graph::ArcIndex ArcIndex;
253 typedef ArcFlowT ArcFlowType;
254 typedef FlowSumT FlowSumType;
255 static_assert(sizeof(FlowSumType) >= sizeof(ArcFlowType));
256 static_assert(std::is_signed_v<FlowSumType>);
257
258 // The height of a node never excess 2 times the number of node, so we
259 // use the same type as a Node index.
261
262 // Initialize a MaxFlow instance on the given graph. The graph does not need
263 // to be fully built yet, but its capacity reservation are used to initialize
264 // the memory of this class. source and sink must also be valid node of
265 // graph.
267
268#ifndef SWIG
269 // This type is neither copyable nor movable.
272#endif
273
274 // Returns the graph associated to the current object.
275 const Graph* graph() const { return graph_; }
276
277 // Returns the status of last call to Solve(). NOT_SOLVED is returned if
278 // Solve() has never been called or if the problem has been modified in such a
279 // way that the previous solution becomes invalid.
280 Status status() const { return status_; }
281
282 // Returns the index of the node corresponding to the source of the network.
284
285 // Returns the index of the node corresponding to the sink of the network.
286 NodeIndex GetSinkNodeIndex() const { return sink_; }
287
288 // Sets the capacity for arc to new_capacity.
289 //
290 // Note that this will be ignored for self-arc, so do not be surprised if you
291 // get zero when reading the capacity of a self-arc back.
292 void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity);
293
294 // Returns true if a maximum flow was solved.
295 bool Solve();
296
297 // Returns the total flow found by the algorithm.
299
300 // Returns the current flow on the given arc.
301 //
302 // Note that on a couple (arc, opposite_arc) the flow only goes in one
303 // direction (where it is positive) and the other direction will have the
304 // negation of that flow.
305 //
306 // Tricky: This returns a FlowSumType in order to be able to distinguish
307 // negative flow from positive flow when ArcFlowType is unsigned.
308 // TODO(user): We could change the API to avoid that if needed.
310 if constexpr (Graph::kHasNegativeReverseArcs) {
311 if (IsArcDirect(arc)) {
312 return static_cast<FlowSumType>(residual_arc_capacity_[Opposite(arc)]);
313 }
314 return -static_cast<FlowSumType>(residual_arc_capacity_[arc]);
315 }
316 return static_cast<FlowSumType>(initial_capacity_[arc]) -
317 static_cast<FlowSumType>(residual_arc_capacity_[arc]);
318 }
319
320 // Returns the initial capacity of an arc.
322 if constexpr (Graph::kHasNegativeReverseArcs) {
323 if (!IsArcDirect(arc)) return 0;
324 return residual_arc_capacity_[arc] +
326 }
327 return initial_capacity_[arc];
328 }
329
330 // Returns the nodes reachable from the source in the residual graph, the
331 // outgoing arcs of this set form a minimum cut.
332 void GetSourceSideMinCut(std::vector<NodeIndex>* result);
333
334 // Returns the nodes that can reach the sink in the residual graph, the
335 // outgoing arcs of this set form a minimum cut. Note that if this is the
336 // complement of GetNodeReachableFromSource(), then the min-cut is unique.
337 //
338 // TODO(user): In the two-phases algorithm, we can get this minimum cut
339 // without doing the second phase. Add an option for this if there is a need
340 // to, note that the second phase is pretty fast so the gain will be small.
341 void GetSinkSideMinCut(std::vector<NodeIndex>* result);
342
343 // Returns true if there exists a path from the source to the sink with
344 // remaining capacity. This allows us to easily check at the end that the flow
345 // we computed is indeed optimal (provided that all the conditions tested by
346 // CheckResult() also hold).
348
349 // Returns the protocol buffer representation of the current problem.
350 FlowModelProto CreateFlowModel();
351
352 protected:
353 // Checks whether the result is valid, i.e. that node excesses are all equal
354 // to zero (we have a flow) and that residual capacities are all non-negative
355 // or zero.
356 bool CheckResult() const;
357
358 // Returns true if arc is admissible.
360 const NodeHeight* potentials) const {
361 DCHECK_EQ(tail, Tail(arc));
362 return residual_arc_capacity_[arc] > 0 &&
363 potentials[tail] == potentials[Head(arc)] + 1;
364 }
365
366 // Returns true if node is active, i.e. if its excess is positive and it
367 // is neither the source or the sink of the graph.
368 bool IsActive(NodeIndex node) const {
369 return (node != source_) && (node != sink_) && (node_excess_[node] > 0);
370 }
371
372 // Returns true if a precondition for Relabel is met, i.e. the outgoing arcs
373 // of node are all either saturated or the heights of their heads are greater
374 // or equal to the height of node.
376
377 // Returns context concatenated with information about arc
378 // in a human-friendly way.
379 std::string DebugString(absl::string_view context, ArcIndex arc) const;
380
381 // Initializes the container active_nodes_.
383
384 // Get the first element from the active node container.
388
389 // Push element to the active node container.
390 void PushActiveNode(const NodeIndex& node) {
391 active_node_by_height_.Push(node, node_potential_[node]);
392 }
393
394 // Check the emptiness of the container.
396
397 // Performs optimization step.
398 void Refine();
400
401 // Discharges an active node node by saturating its admissible adjacent arcs,
402 // if any, and by relabelling it when it becomes inactive.
404
405 // Initializes the preflow to a state that enables to run Refine.
407
408 // Clears the flow excess at each node by pushing the flow back to the source:
409 // - Do a depth-first search from the source in the direct graph to cancel
410 // flow cycles.
411 // - Then, return flow excess along the depth-first search tree (by pushing
412 // the flow in the reverse dfs topological order).
413 // The theoretical complexity is O(mn), but it is a lot faster in practice.
415
416 // Computes the best possible node potential given the current flow using a
417 // reverse breadth-first search from the sink in the reverse residual graph.
418 // This is an implementation of the global update heuristic mentioned in many
419 // max-flow papers. See for instance: B.V. Cherkassky, A.V. Goldberg, "On
420 // implementing push-relabel methods for the maximum flow problem",
421 // Algorithmica, 19:390-410, 1997.
422 // ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/94/1523/CS-TR-94-1523.pdf
424
425 // Tries to saturate all the outgoing arcs from the source that can reach the
426 // sink. Most of the time, we can do that in one go, except when more flow
427 // than kMaxFlowSum can be pushed out of the source in which case we have
428 // to be careful. Returns true if some flow was pushed.
430
431 // Pushes flow on arc, i.e. consumes flow on residual_arc_capacity_[arc],
432 // and consumes -flow on residual_arc_capacity_[Opposite(arc)]. Updates
433 // node_excess_ at the tail and head of arc accordingly.
434 void PushFlow(ArcFlowType flow, NodeIndex tail, ArcIndex arc);
435
436 // Similar to PushFlow(), but this will always push
437 // std::min(node_excess_[tail], residual_arc_capacity_[arc]).
438 //
439 // We also use a cached node_excess_.data() as this is used in hot loops and
440 // allow to remove bound checking and fetching the pointer which is constant.
442 FlowSumType* node_excess);
443
444 // Relabels a node, i.e. increases its height by the minimum necessary amount.
445 // This version of Relabel is relaxed in a way such that if an admissible arc
446 // exists at the current node height, then the node is not relabeled. This
447 // enables us to deal with wrong values of first_admissible_arc_[node] when
448 // updating it is too costly.
449 void Relabel(NodeIndex node);
450
451 // Handy member functions to make the code more compact.
452 NodeIndex Head(ArcIndex arc) const { return graph_->Head(arc); }
453 NodeIndex Tail(ArcIndex arc) const { return graph_->Tail(arc); }
455 bool IsArcDirect(ArcIndex arc) const;
456 bool IsArcValid(ArcIndex arc) const;
457
458 // Returns the set of nodes reachable from start in the residual graph or in
459 // the reverse residual graph (if reverse is true).
460 template <bool reverse>
461 void ComputeReachableNodes(NodeIndex start, std::vector<NodeIndex>* result);
462
463 // Maximum manageable flow.
464 static constexpr FlowSumType kMaxFlowSum =
465 std::numeric_limits<FlowSumType>::max();
466
467 // A pointer to the graph passed as argument.
468 const Graph* graph_;
469
470 // An array representing the excess for each node in graph_.
471 // With the current algorithm this is always positive except at the source.
472 // TODO(user): If needed we could tweak that to allow unsigned FlowSumType.
473 std::unique_ptr<FlowSumType[]> node_excess_;
474
475 // An array representing the height function for each node in graph_. For a
476 // given node, this is a lower bound on the shortest path length from this
477 // node to the sink in the residual network. The height of a node always goes
478 // up during the course of a Solve().
479 //
480 // Since initially we saturate all the outgoing arcs of the source, we can
481 // never reach the sink from the source in the residual graph. Initially we
482 // set the height of the source to n (the number of node of the graph) and it
483 // never changes. If a node as an height >= n, then this node can't reach the
484 // sink and its height minus n is a lower bound on the shortest path length
485 // from this node to the source in the residual graph.
486 std::unique_ptr<NodeHeight[]> node_potential_;
487
488 // An array representing the residual_capacity for each arc in graph_.
489 // Residual capacities enable one to represent the capacity and flow for all
490 // arcs in the graph in the following manner.
491 // For all arc, residual_arc_capacity_[arc] = capacity[arc] - flow[arc]
492 // Moreover, for reverse arcs, capacity[arc] = 0 by definition,
493 // Also flow[Opposite(arc)] = -flow[arc] by definition.
494 // Therefore:
495 // - for a direct arc:
496 // flow[arc] = 0 - flow[Opposite(arc)]
497 // = capacity[Opposite(arc)] - flow[Opposite(arc)]
498 // = residual_arc_capacity_[Opposite(arc)]
499 // - for a reverse arc:
500 // flow[arc] = -residual_arc_capacity_[arc]
501 // Using these facts enables one to only maintain residual_arc_capacity_,
502 // instead of both capacity and flow, for each direct and indirect arc. This
503 // reduces the amount of memory for this information by a factor 2.
505
506 // The initial capacity as set by SetArcCapacity(), unused if
507 // `Graph::kNegativeReverseArcs`, as we can always recover the initial
508 // capacity from the residual capacities of an arc and its reverse.
509 std::unique_ptr<ArcFlowType[]> initial_capacity_;
510
511 // An array representing the first admissible arc for each node in graph_.
512 std::unique_ptr<ArcIndex[]> first_admissible_arc_;
513
514 // A priority queue used for managing active nodes in the algorithm. It allows
515 // to select the active node with highest height before each Discharge().
516 // Moreover, since all pushes from this node will be to nodes with height
517 // greater or equal to the initial discharged node height minus one, the
518 // PriorityQueueWithRestrictedPush is a perfect fit.
520
521 // The index of the source node in graph_.
523
524 // The index of the sink node in graph_.
526
527 // The status of the problem.
529
530 // BFS queue used by the GlobalUpdate() function. We do not use a C++ queue
531 // because we need access to the vector for different optimizations.
532 std::vector<bool> node_in_bfs_queue_;
533 std::vector<NodeIndex> bfs_queue_;
534
535 // Statistics about this class.
537};
538
539template <typename Element, typename IntegerPriority>
541 const {
542 return even_queue_.empty() && odd_queue_.empty();
543}
544
545template <typename Element, typename IntegerPriority>
547 even_queue_.clear();
548 odd_queue_.clear();
549}
550
551template <typename Element, typename IntegerPriority>
553 Element element, IntegerPriority priority) {
554 // Since users may rely on it, we DCHECK the exact condition.
555 DCHECK(even_queue_.empty() || priority >= even_queue_.back().second - 1);
556 DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second - 1);
557
558 // Note that the DCHECK() below are less restrictive than the ones above but
559 // check a necessary and sufficient condition for the priority queue to behave
560 // as expected.
561 if (priority & 1) {
562 DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second);
563 odd_queue_.push_back(std::make_pair(element, priority));
564 } else {
565 DCHECK(even_queue_.empty() || priority >= even_queue_.back().second);
566 even_queue_.push_back(std::make_pair(element, priority));
567 }
568}
569
570template <typename Element, typename IntegerPriority>
572 DCHECK(!IsEmpty());
573 if (even_queue_.empty()) return PopBack(&odd_queue_);
574 if (odd_queue_.empty()) return PopBack(&even_queue_);
575 if (odd_queue_.back().second > even_queue_.back().second) {
576 return PopBack(&odd_queue_);
577 } else {
578 return PopBack(&even_queue_);
579 }
580}
581
582template <typename Element, typename IntegerPriority>
583Element PriorityQueueWithRestrictedPush<Element, IntegerPriority>::PopBack(
584 std::vector<std::pair<Element, IntegerPriority> >* queue) {
585 DCHECK(!queue->empty());
586 Element element = queue->back().first;
587 queue->pop_back();
588 return element;
589}
590
591template <typename Graph, typename ArcFlowT, typename FlowSumT>
593 NodeIndex source,
594 NodeIndex sink)
595 : graph_(graph),
597 source_(source),
598 sink_(sink),
599 stats_("MaxFlow") {
601 DCHECK(graph->IsNodeValid(source));
602 DCHECK(graph->IsNodeValid(sink));
603 const NodeIndex max_num_nodes = graph_->node_capacity();
604 if (max_num_nodes > 0) {
605 // We will initialize them in InitializePreflow(), so no need for memset.
606 //
607 // TODO(user): Unfortunately std/absl::make_unique_for_overwrite is not
608 // yet available in open-source.
609 node_excess_ = std::make_unique<FlowSumType[]>(max_num_nodes);
610 node_potential_ = std::make_unique<NodeHeight[]>(max_num_nodes);
611 first_admissible_arc_ = std::make_unique<ArcIndex[]>(max_num_nodes);
612 bfs_queue_.reserve(max_num_nodes);
613 }
614 const ArcIndex max_num_arcs = graph_->arc_capacity();
615 if (max_num_arcs > 0) {
616 if constexpr (Graph::kHasNegativeReverseArcs) {
617 residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1);
618 } else {
619 // We will need to store the initial capacity in this case.
620 initial_capacity_ = std::make_unique<ArcFlowType[]>(max_num_arcs);
621 residual_arc_capacity_.Reserve(0, max_num_arcs - 1);
622 }
623 residual_arc_capacity_.SetAll(0);
624 }
625}
626
627template <typename Graph, typename ArcFlowT, typename FlowSumT>
629 ArcIndex arc, ArcFlowType new_capacity) {
631 DCHECK_LE(0, new_capacity);
632 DCHECK(IsArcDirect(arc));
633
634 // This serves no purpose from a max-flow point of view, so it is safer to
635 // just leave the capacity of all self-arc at zero.
636 if (Head(arc) == Tail(arc)) return;
637
639 residual_arc_capacity_[arc] = new_capacity;
640
641 // Since the class might have already be used, we need to make sure we clear
642 // any previous state on this arc and its reverse.
643 if constexpr (Graph::kHasNegativeReverseArcs) {
645 } else {
646 initial_capacity_[arc] = new_capacity;
647 DCHECK_LE(initial_capacity_[arc], std::numeric_limits<ArcFlowType>::max() -
649 }
650}
651
652template <typename Graph, typename ArcFlowT, typename FlowSumT>
657
658template <typename Graph, typename ArcFlowT, typename FlowSumT>
660 std::vector<NodeIndex>* result) {
662}
663
664template <typename Graph, typename ArcFlowT, typename FlowSumT>
668 LOG(DFATAL) << "-node_excess_[source_] = " << -node_excess_[source_]
669 << " != node_excess_[sink_] = " << node_excess_[sink_];
670 return false;
671 }
672 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
673 if (node != source_ && node != sink_) {
674 if (node_excess_[node] != 0) {
675 LOG(DFATAL) << "node_excess_[" << node << "] = " << node_excess_[node]
676 << " != 0";
677 return false;
678 }
679 }
680 }
681 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
682 const ArcIndex opposite = Opposite(arc);
683 const ArcFlowType direct_capacity = residual_arc_capacity_[arc];
684 const ArcFlowType opposite_capacity = residual_arc_capacity_[opposite];
685 if (direct_capacity < 0) {
686 LOG(DFATAL) << "residual_arc_capacity_[" << arc
687 << "] = " << direct_capacity << " < 0";
688 return false;
689 }
690 if (opposite_capacity < 0) {
691 LOG(DFATAL) << "residual_arc_capacity_[" << opposite
692 << "] = " << opposite_capacity << " < 0";
693 return false;
694 }
695 // The initial capacity of the direct arcs is non-negative.
696 if (direct_capacity + opposite_capacity < 0) {
697 LOG(DFATAL) << "initial capacity [" << arc
698 << "] = " << direct_capacity + opposite_capacity << " < 0";
699 return false;
700 }
701 }
702
704 LOG(DFATAL) << "The algorithm terminated, but the flow is not maximal!";
705 return false;
706 }
707
708 return true;
709}
710
711template <typename Graph, typename ArcFlowT, typename FlowSumT>
714
715 // We simply compute the reachability from the source in the residual graph.
716 const NodeIndex num_nodes = graph_->num_nodes();
717 std::vector<bool> is_reached(num_nodes, false);
718 std::vector<NodeIndex> to_process;
719
720 to_process.push_back(source_);
721 is_reached[source_] = true;
722 while (!to_process.empty()) {
723 const NodeIndex node = to_process.back();
724 to_process.pop_back();
725 for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
726 if (residual_arc_capacity_[arc] > 0) {
727 const NodeIndex head = graph_->Head(arc);
728 if (!is_reached[head]) {
729 is_reached[head] = true;
730 to_process.push_back(head);
731 }
732 }
733 }
734 }
735 return is_reached[sink_];
736}
737
738template <typename Graph, typename ArcFlowT, typename FlowSumT>
740 NodeIndex node) const {
741 DCHECK(IsActive(node));
742 for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
743 DCHECK(!IsAdmissible(node, arc, node_potential_.data()))
744 << DebugString("CheckRelabelPrecondition:", arc);
745 }
746 return true;
747}
748
749template <typename Graph, typename ArcFlowT, typename FlowSumT>
751 absl::string_view context, ArcIndex arc) const {
752 const NodeIndex tail = Tail(arc);
753 const NodeIndex head = Head(arc);
754 return absl::StrFormat(
755 "%s Arc %d, from %d to %d, "
756 "Residual capacity = %d, "
757 "Residual capacity for reverse arc = %d, "
758 "Height(tail) = %d, Height(head) = %d, "
759 "Excess(tail) = %d, Excess(head) = %d",
760 context, arc, tail, head, residual_arc_capacity_[arc],
762 node_potential_[head], node_excess_[tail], node_excess_[head]);
763}
764
765template <typename Graph, typename ArcFlowT, typename FlowSumT>
769
770 // Deal with the case when source_ or sink_ is not inside graph_.
771 // Since they are both specified independently of the graph, we do need to
772 // take care of this corner case.
773 const NodeIndex num_nodes = graph_->num_nodes();
774 if (sink_ >= num_nodes || source_ >= num_nodes) {
775 // Behave like a normal graph where source_ and sink_ are disconnected.
776 // Note that the arc flow is set to 0 by InitializePreflow().
778 return true;
779 }
780
782
784 DCHECK(CheckResult());
785
787 // In this case, we are sure that the flow is > kMaxFlowSum.
789 }
790 IF_STATS_ENABLED(VLOG(1) << stats_.StatString());
791 return true;
792}
793
794template <typename Graph, typename ArcFlowT, typename FlowSumT>
797
798 // TODO(user): Ebert graph has an issue with nodes with no arcs, so we
799 // use max_num_nodes here to resize vectors.
800 const NodeIndex num_nodes = graph_->num_nodes();
801 const NodeIndex max_num_nodes = graph_->node_capacity();
802
803 // InitializePreflow() clears the whole flow that could have been computed
804 // by a previous Solve(). This is not optimal in terms of complexity.
805 //
806 // TODO(user): find a way to make the re-solving incremental (not an obvious
807 // task, and there has not been a lot of literature on the subject.)
808 std::fill(node_excess_.get(), node_excess_.get() + max_num_nodes, 0);
809
810 // Restart from a clear state with no flow, and an initial arc capacity.
811 const ArcIndex num_arcs = graph_->num_arcs();
812 if constexpr (Graph::kHasNegativeReverseArcs) {
813 for (ArcIndex arc = 0; arc < num_arcs; ++arc) {
814 const ArcIndex opposite_arc = Opposite(arc);
816 residual_arc_capacity_[opposite_arc] = 0;
817 }
818 } else {
819 for (ArcIndex arc = 0; arc < num_arcs; ++arc) {
821 }
822 }
823
824 // All the initial heights are zero except for the source whose height is
825 // equal to the number of nodes and will never change during the algorithm.
826 std::fill(node_potential_.get(), node_potential_.get() + max_num_nodes, 0);
827 node_potential_[source_] = num_nodes;
828
829 // Initially we set first_admissible_arc_ to the first arc in the iteration.
830 for (NodeIndex node = 0; node < num_nodes; ++node) {
832 for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
833 first_admissible_arc_[node] = arc;
834 break;
835 }
836 }
837}
838
839// Note(user): Calling this function will break the property on the node
840// potentials because of the way we cancel flow on cycle. However, we only call
841// that at the end of the algorithm, or just before a GlobalUpdate() that will
842// restore the precondition on the node potentials.
843template <typename Graph, typename ArcFlowT, typename FlowSumT>
846 const NodeIndex num_nodes = graph_->num_nodes();
847
848 // We implement a variation of Tarjan's strongly connected component algorithm
849 // to detect cycles published in: Tarjan, R. E. (1972), "Depth-first search
850 // and linear graph algorithms", SIAM Journal on Computing. A description can
851 // also be found in wikipedia.
852
853 // Stored nodes are settled nodes already stored in the
854 // reverse_topological_order (except the sink_ that we do not actually store).
855 std::vector<bool> stored(num_nodes, false);
856 stored[sink_] = true;
857
858 // The visited nodes that are not yet stored are all the nodes from the
859 // source_ to the current node in the current dfs branch.
860 std::vector<bool> visited(num_nodes, false);
861 visited[sink_] = true;
862
863 // Stack of arcs to explore in the dfs search.
864 // The current node is Head(arc_stack.back()).
865 std::vector<ArcIndex> arc_stack;
866
867 // Increasing list of indices into the arc_stack that correspond to the list
868 // of arcs in the current dfs branch from the source_ to the current node.
869 std::vector<int> index_branch;
870
871 // Node in reverse_topological_order in the final dfs tree.
872 std::vector<NodeIndex> reverse_topological_order;
873
874 // We start by pushing all the outgoing arcs from the source on the stack to
875 // avoid special conditions in the code. As a result, source_ will not be
876 // stored in reverse_topological_order, and this is what we want.
877 for (const ArcIndex arc : graph_->OutgoingArcs(source_)) {
878 if (Flow(arc) > 0) {
879 arc_stack.push_back(arc);
880 }
881 }
882 visited[source_] = true;
883
884 // Start the dfs on the subgraph formed by the direct arcs with positive flow.
885 while (!arc_stack.empty()) {
886 DCHECK_GT(Flow(arc_stack.back()), 0)
887 << arc_stack.size() - 1 << " arc " << arc_stack.back() << " "
888 << Tail(arc_stack.back()) << "->" << Head(arc_stack.back());
889 const NodeIndex node = Head(arc_stack.back());
890
891 // If the node is visited, it means we have explored all its arcs and we
892 // have just backtracked in the dfs. Store it if it is not already stored
893 // and process the next arc on the stack.
894 if (visited[node]) {
895 if (!stored[node]) {
896 stored[node] = true;
897 reverse_topological_order.push_back(node);
898 DCHECK(!index_branch.empty());
899 index_branch.pop_back();
900 }
901 arc_stack.pop_back();
902 continue;
903 }
904
905 // The node is a new unexplored node, add all its outgoing arcs with
906 // positive flow to the stack and go deeper in the dfs.
907 DCHECK(!stored[node]);
908 DCHECK(index_branch.empty() ||
909 (arc_stack.size() - 1 > index_branch.back()));
910 visited[node] = true;
911 index_branch.push_back(arc_stack.size() - 1);
912 for (const ArcIndex arc : graph_->OutgoingArcs(node)) {
913 const FlowSumType flow = Flow(arc);
914 const NodeIndex head = Head(arc);
915 if (flow > 0 && !stored[head]) {
916 if (!visited[head]) {
917 arc_stack.push_back(arc);
918 } else {
919 // There is a cycle.
920 // Find the first index to consider,
921 // arc_stack[index_branch[cycle_begin]] will be the first arc on the
922 // cycle.
923 int cycle_begin = index_branch.size();
924 while (cycle_begin > 0 &&
925 Head(arc_stack[index_branch[cycle_begin - 1]]) != head) {
926 --cycle_begin;
927 }
928
929 // Compute the maximum flow that can be canceled on the cycle and the
930 // min index such that arc_stack[index_branch[i]] will be saturated.
931 ArcFlowType flow_on_cycle = static_cast<ArcFlowType>(flow);
932 int first_saturated_index = index_branch.size();
933 for (int i = index_branch.size() - 1; i >= cycle_begin; --i) {
934 const ArcIndex arc_on_cycle = arc_stack[index_branch[i]];
935 if (const ArcFlowType arc_flow = Flow(arc_on_cycle);
936 arc_flow <= flow_on_cycle) {
937 flow_on_cycle = arc_flow;
938 first_saturated_index = i;
939 }
940 }
941
942 // This is just here for a DCHECK() below.
943 const FlowSumType excess = node_excess_[head];
944
945 // Cancel the flow on the cycle, and set visited[node] = false for
946 // the node that will be backtracked over.
947 PushFlow(-flow_on_cycle, node, arc);
948 for (int i = index_branch.size() - 1; i >= cycle_begin; --i) {
949 const ArcIndex arc_on_cycle = arc_stack[index_branch[i]];
950 PushFlow(-flow_on_cycle, Tail(arc_on_cycle), arc_on_cycle);
951 if (i >= first_saturated_index) {
952 DCHECK(visited[Head(arc_on_cycle)]);
953 visited[Head(arc_on_cycle)] = false;
954 } else {
955 DCHECK_GT(Flow(arc_on_cycle), 0);
956 }
957 }
958
959 // This is a simple check that the flow was pushed properly.
960 DCHECK_EQ(excess, node_excess_[head]);
961
962 // Backtrack the dfs just before index_branch[first_saturated_index].
963 // If the current node is still active, there is nothing to do.
964 if (first_saturated_index < index_branch.size()) {
965 arc_stack.resize(index_branch[first_saturated_index]);
966 index_branch.resize(first_saturated_index);
967
968 // We backtracked over the current node, so there is no need to
969 // continue looping over its arcs.
970 break;
971 }
972 }
973 }
974 }
975 }
976 DCHECK(arc_stack.empty());
977 DCHECK(index_branch.empty());
978
979 // Return the flow to the source. Note that the sink_ and the source_ are not
980 // stored in reverse_topological_order.
981 for (int i = 0; i < reverse_topological_order.size(); i++) {
982 const NodeIndex node = reverse_topological_order[i];
983 if (node_excess_[node] == 0) continue;
984 for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
985 const FlowSumType flow = Flow(arc);
986 if (flow < 0) {
987 DCHECK_GT(residual_arc_capacity_[arc], 0);
988 const ArcFlowType to_push =
989 static_cast<ArcFlowType>(std::min(node_excess_[node], -flow));
990 PushFlow(to_push, node, arc);
991 if (node_excess_[node] == 0) break;
992 }
993 }
994 DCHECK_EQ(0, node_excess_[node]);
995 }
996 DCHECK_EQ(-node_excess_[source_], node_excess_[sink_]);
997}
998
999template <typename Graph, typename ArcFlowT, typename FlowSumT>
1002 bfs_queue_.clear();
1003 int queue_index = 0;
1004 const NodeIndex num_nodes = graph_->num_nodes();
1005 node_in_bfs_queue_.assign(num_nodes, false);
1006 node_in_bfs_queue_[sink_] = true;
1007
1008 // We cache this as this is a hot-loop and we don't want bound checking.
1009 FlowSumType* const node_excess = node_excess_.get();
1010
1011 // We do a BFS in the reverse residual graph, starting from the sink.
1012 // Because all the arcs from the source are saturated (except in
1013 // presence of integer overflow), the source cannot reach the sink in the
1014 // residual graph. If there is possible overflow and the source is reachable,
1015 // we still do not want to relabel it, so we start with the source marked.
1017 bfs_queue_.push_back(sink_);
1018
1019 while (queue_index != bfs_queue_.size()) {
1020 const NodeIndex node = bfs_queue_[queue_index];
1021 ++queue_index;
1022 const NodeIndex candidate_distance = node_potential_[node] + 1;
1023 for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
1024 const NodeIndex head = Head(arc);
1025
1026 // Skip the arc if the height of head was already set to the correct
1027 // value (Remember we are doing reverse BFS).
1028 if (node_in_bfs_queue_[head]) continue;
1029
1030 // TODO(user): By using more memory we can speed this up quite a bit by
1031 // avoiding to take the opposite arc here, too options:
1032 // - if (residual_arc_capacity_[arc] != arc_capacity_[arc])
1033 // - if (opposite_arc_is_admissible_[arc]) // need updates.
1034 // Experiment with the first option shows more than 10% gain on this
1035 // function running time, which is the bottleneck on many instances.
1036 const ArcIndex opposite_arc = Opposite(arc);
1037 if (residual_arc_capacity_[opposite_arc] > 0) {
1038 // Note(user): We used to have a DCHECK_GE(candidate_distance,
1039 // node_potential_[head]); which is always true except in the case
1040 // where we can push more than kMaxFlowSum out of the source. The
1041 // problem comes from the fact that in this case, we call
1042 // PushFlowExcessBackToSource() in the middle of the algorithm. The
1043 // later call will break the properties of the node potential. Note
1044 // however, that this function will recompute a good node potential
1045 // for all the nodes and thus fix the issue.
1046
1047 // If head is active, we can steal some or all of its excess.
1048 // This brings a huge gain on some problems.
1049 // Note(user): I haven't seen this anywhere in the literature.
1050 // TODO(user): Investigate more and maybe write a publication :)
1051 if (node_excess[head] > 0) {
1052 PushAsMuchFlowAsPossible(head, opposite_arc, node_excess);
1053
1054 // If the arc became saturated, it is no longer in the residual
1055 // graph, so we do not need to consider head at this time.
1056 if (residual_arc_capacity_[opposite_arc] == 0) continue;
1057 }
1058
1059 // Note that there is no need to touch first_admissible_arc_[node]
1060 // because of the relaxed Relabel() we use.
1061 node_potential_[head] = candidate_distance;
1062 node_in_bfs_queue_[head] = true;
1063 bfs_queue_.push_back(head);
1064 }
1065 }
1066 }
1067
1068 // At the end of the search, some nodes may not be in the bfs_queue_. Such
1069 // nodes cannot reach the sink_ or source_ in the residual graph, so there is
1070 // no point trying to push flow toward them. We obtain this effect by setting
1071 // their height to something unreachable.
1072 //
1073 // Note that this also prevents cycling due to our anti-overflow procedure.
1074 // For instance, suppose there is an edge s -> n outgoing from the source. If
1075 // node n has no other connection and some excess, we will push the flow back
1076 // to the source, but if we don't update the height of n
1077 // SaturateOutgoingArcsFromSource() will push the flow to n again.
1078 // TODO(user): This is another argument for another anti-overflow algorithm.
1079 for (NodeIndex node = 0; node < num_nodes; ++node) {
1080 if (!node_in_bfs_queue_[node]) {
1081 node_potential_[node] = 2 * num_nodes - 1;
1082 }
1083 }
1084
1085 // Reset the active nodes. Doing it like this pushes the nodes in increasing
1086 // order of height. Note that bfs_queue_[0] is the sink_ so we skip it.
1088 for (int i = 1; i < bfs_queue_.size(); ++i) {
1089 const NodeIndex node = bfs_queue_[i];
1090 if (node_excess_[node] > 0) {
1091 DCHECK(IsActive(node));
1092 PushActiveNode(node);
1093 }
1094 }
1095}
1096
1097template <typename Graph, typename ArcFlowT, typename FlowSumT>
1098bool GenericMaxFlow<Graph, ArcFlowT,
1101 const NodeIndex num_nodes = graph_->num_nodes();
1102
1103 // If sink_ or source_ already have kMaxFlowSum, then there is no
1104 // point pushing more flow since it will cause an integer overflow.
1105 if (node_excess_[sink_] == kMaxFlowSum) return false;
1106 if (node_excess_[source_] == -kMaxFlowSum) return false;
1107
1108 bool flow_pushed = false;
1109 for (const ArcIndex arc : graph_->OutgoingArcs(source_)) {
1110 const ArcFlowType flow = residual_arc_capacity_[arc];
1111
1112 // This is a special IsAdmissible() condition for the source.
1113 if (flow == 0 || node_potential_[Head(arc)] >= num_nodes) continue;
1114
1115 // We are careful in case the sum of the flow out of the source is greater
1116 // than kMaxFlowSum to avoid overflow.
1117 const FlowSumType current_flow_out_of_source = -node_excess_[source_];
1118 DCHECK_GE(flow, 0) << flow;
1119 DCHECK_GE(current_flow_out_of_source, 0) << current_flow_out_of_source;
1120 const FlowSumType capped_flow = kMaxFlowSum - current_flow_out_of_source;
1121 if (capped_flow < static_cast<FlowSumType>(flow)) {
1122 // We push as much flow as we can so the current flow on the network will
1123 // be kMaxFlowSum.
1124
1125 // Since at the beginning of the function, current_flow_out_of_source
1126 // was different from kMaxFlowSum, we are sure to have pushed some
1127 // flow before if capped_flow is 0.
1128 if (capped_flow == 0) return true;
1129 PushFlow(static_cast<ArcFlowType>(capped_flow), source_, arc);
1130 return true;
1131 }
1132 PushFlow(flow, source_, arc);
1133 flow_pushed = true;
1134 }
1135 DCHECK_LE(node_excess_[source_], 0);
1136 return flow_pushed;
1137}
1138
1139template <typename Graph, typename ArcFlowT, typename FlowSumT>
1141 NodeIndex tail,
1142 ArcIndex arc) {
1144
1145 // Update the residual capacity of the arc and its opposite arc.
1146 DCHECK_NE(flow, 0);
1147 residual_arc_capacity_[arc] -= flow;
1148 residual_arc_capacity_[Opposite(arc)] += flow;
1149 DCHECK_GE(residual_arc_capacity_[arc], 0);
1150 DCHECK_GE(residual_arc_capacity_[Opposite(arc)], 0);
1151
1152 // Update the excesses at the tail and head of the arc.
1153 //
1154 // Note(user): node_excess_ should be always greater than or equal to 0 except
1155 // for the source where it should always be smaller than or equal to 0. Note
1156 // however that we cannot check this because when we cancel the flow on a
1157 // cycle in PushFlowExcessBackToSource(), we may break this invariant during
1158 // the operation even if it is still valid at the end.
1159 node_excess_[tail] -= static_cast<FlowSumType>(flow);
1160 node_excess_[Head(arc)] += static_cast<FlowSumType>(flow);
1161}
1162
1163template <typename Graph, typename ArcFlowT, typename FlowSumT>
1165 NodeIndex tail, ArcIndex arc, FlowSumType* node_excess) {
1167
1168 // We either push all node_excess or saturate the arc by consuming all its
1169 // residual capacity.
1170 const ArcFlowType flow = static_cast<ArcFlowType>(
1171 std::min(static_cast<FlowSumType>(residual_arc_capacity_[arc]),
1172 node_excess[tail]));
1173
1174 DCHECK_NE(flow, 0);
1175 residual_arc_capacity_[arc] -= flow;
1176 residual_arc_capacity_[Opposite(arc)] += flow;
1177 DCHECK_GE(residual_arc_capacity_[arc], 0);
1178 DCHECK_GE(residual_arc_capacity_[Opposite(arc)], 0);
1179
1180 node_excess[tail] -= static_cast<FlowSumType>(flow);
1181 node_excess[Head(arc)] += static_cast<FlowSumType>(flow);
1182}
1183
1184template <typename Graph, typename ArcFlowT, typename FlowSumT>
1185void GenericMaxFlow<Graph, ArcFlowT,
1189 const NodeIndex num_nodes = graph_->num_nodes();
1190 for (NodeIndex node = 0; node < num_nodes; ++node) {
1191 if (IsActive(node)) {
1192 // This node cannot reach the sink in the residual graph, we don't need to
1193 // consider it anymore, and we will just push its potential excess back to
1194 // the source in PushFlowExcessBackToSource() at the end.
1195 if (node_potential_[node] >= num_nodes) continue;
1196 PushActiveNode(node);
1197 }
1198 }
1199}
1200
1201template <typename Graph, typename ArcFlowT, typename FlowSumT>
1204
1205 const NodeIndex num_nodes = graph_->num_nodes();
1206 std::vector<int> skip_active_node;
1207
1208 // Usually SaturateOutgoingArcsFromSource() will saturate all the arcs from
1209 // the source in one go, and we will loop just once. But in case we can push
1210 // more than kMaxFlowSum out of the source the loop is as follow:
1211 // - Push up to kMaxFlowSum out of the source on the admissible outgoing
1212 // arcs. Stop if no flow was pushed.
1213 // - Compute the current max-flow. This will push some flow back to the
1214 // source and render more outgoing arcs from the source not admissible.
1215 //
1216 // TODO(user): This may not be the most efficient algorithm if we need to loop
1217 // many times. An alternative may be to handle the source like the other nodes
1218 // in the algorithm, initially putting an excess of kMaxFlowSum on it,
1219 // and making the source active like any other node with positive excess. To
1220 // investigate.
1222 int num_skipped;
1223 do {
1224 num_skipped = 0;
1225 skip_active_node.assign(num_nodes, 0);
1226 skip_active_node[sink_] = 2;
1227 skip_active_node[source_] = 2;
1228 GlobalUpdate();
1229 while (!IsEmptyActiveNodeContainer()) {
1231 if (skip_active_node[node] > 1) {
1232 if (node != sink_ && node != source_) ++num_skipped;
1233 continue;
1234 }
1235 const NodeIndex old_height = node_potential_[node];
1236 Discharge(node);
1237
1238 // The idea behind this is that if a node height augments by more than
1239 // one, then it is likely to push flow back the way it came. This can
1240 // lead to very costly loops. A bad case is: source -> n1 -> n2 and n2
1241 // just recently isolated from the sink. Then n2 will push flow back to
1242 // n1, and n1 to n2 and so on. The height of each node will increase by
1243 // steps of two until the height of the source is reached, which can
1244 // take a long time. If the chain is longer, the situation is even
1245 // worse. The behavior of this heuristic is related to the Gap
1246 // heuristic.
1247 //
1248 // Note that the global update will fix all such cases efficiently. So
1249 // the idea is to discharge the active node as much as possible, and
1250 // then do a global update.
1251 //
1252 // We skip a node when this condition was true 2 times to avoid doing a
1253 // global update too frequently.
1254 if (node_potential_[node] > old_height + 1) {
1255 ++skip_active_node[node];
1256 }
1257 }
1258 } while (num_skipped > 0);
1259
1260 // We use a two-phase algorithm:
1261 // 1/ Only deal with nodes that can reach the sink. At the end we know the
1262 // value of the maximum flow and we have a min-cut.
1263 // 2/ Call PushFlowExcessBackToSource() to obtain a max-flow. This is
1264 // usually a lot faster than the first phase.
1266 }
1267}
1268
1269template <typename Graph, typename ArcFlowT, typename FlowSumT>
1271 const NodeIndex node) {
1273 const NodeIndex num_nodes = graph_->num_nodes();
1274
1275 // We cache this as this is a hot-loop and we don't want bound checking.
1276 FlowSumType* const node_excess = node_excess_.get();
1277 NodeHeight* const node_potentials = node_potential_.get();
1278 ArcIndex* const first_admissible_arc = first_admissible_arc_.get();
1279
1280 while (true) {
1281 DCHECK(IsActive(node));
1282 for (const ArcIndex arc :
1283 graph_->OutgoingOrOppositeIncomingArcsStartingFrom(
1284 node, first_admissible_arc_[node])) {
1285 if (IsAdmissible(node, arc, node_potentials)) {
1286 DCHECK(IsActive(node));
1287 const NodeIndex head = Head(arc);
1288 if (node_excess[head] == 0) {
1289 // The push below will make the node active for sure. Note that we may
1290 // push the sink_, but that is handled properly in Refine().
1291 PushActiveNode(head);
1292 }
1293 PushAsMuchFlowAsPossible(node, arc, node_excess);
1294 if (node_excess[node] == 0) {
1295 first_admissible_arc[node] = arc; // arc may still be admissible.
1296 return;
1297 }
1298 }
1299 }
1300 Relabel(node);
1301
1302 // This node can no longer reach the sink, skip until
1303 // PushFlowExcessBackToSource().
1304 if (node_potentials[node] >= num_nodes) break;
1305 }
1306}
1307
1308template <typename Graph, typename ArcFlowT, typename FlowSumT>
1311 // Because we use a relaxed version, this is no longer true if the
1312 // first_admissible_arc_[node] was not actually the first arc!
1313 // DCHECK(CheckRelabelPrecondition(node));
1314 NodeHeight min_height = std::numeric_limits<NodeHeight>::max();
1315 ArcIndex first_admissible_arc = Graph::kNilArc;
1316 for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
1317 if (residual_arc_capacity_[arc] > 0) {
1318 // Update min_height only for arcs with available capacity.
1319 NodeHeight head_height = node_potential_[Head(arc)];
1320 if (head_height < min_height) {
1321 min_height = head_height;
1322 first_admissible_arc = arc;
1323
1324 // We found an admissible arc at the current height, just stop there.
1325 // This is the true first_admissible_arc_[node].
1326 if (min_height + 1 == node_potential_[node]) break;
1327 }
1328 }
1329 }
1330 DCHECK_NE(first_admissible_arc, Graph::kNilArc);
1331 node_potential_[node] = min_height + 1;
1332
1333 // Note that after a Relabel(), the loop will continue in Discharge(), and
1334 // we are sure that all the arcs before first_admissible_arc are not
1335 // admissible since their height is > min_height.
1336 first_admissible_arc_[node] = first_admissible_arc;
1337}
1338
1339template <typename Graph, typename ArcFlowT, typename FlowSumT>
1341 ArcIndex arc) const {
1342 return graph_->OppositeArc(arc);
1343}
1344
1345template <typename Graph, typename ArcFlowT, typename FlowSumT>
1347 ArcIndex arc) const {
1348 return IsArcValid(arc) && arc >= 0;
1349}
1350
1351template <typename Graph, typename ArcFlowT, typename FlowSumT>
1353 return graph_->IsArcValid(arc);
1354}
1355
1356template <typename Graph, typename ArcFlowT, typename FlowSumT>
1357template <bool reverse>
1359 NodeIndex start, std::vector<NodeIndex>* result) {
1360 // If start is not a valid node index, it can reach only itself.
1361 // Note(user): This is needed because source and sink are given independently
1362 // of the graph and sometimes before it is even constructed.
1363 const NodeIndex num_nodes = graph_->num_nodes();
1364 if (start >= num_nodes) {
1365 result->clear();
1366 result->push_back(start);
1367 return;
1368 }
1369 bfs_queue_.clear();
1370 node_in_bfs_queue_.assign(num_nodes, false);
1371
1372 int queue_index = 0;
1373 bfs_queue_.push_back(start);
1374 node_in_bfs_queue_[start] = true;
1375 while (queue_index != bfs_queue_.size()) {
1376 const NodeIndex node = bfs_queue_[queue_index];
1377 ++queue_index;
1378 for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) {
1379 const NodeIndex head = Head(arc);
1380 if (node_in_bfs_queue_[head]) continue;
1381 if (residual_arc_capacity_[reverse ? Opposite(arc) : arc] == 0) continue;
1382 node_in_bfs_queue_[head] = true;
1383 bfs_queue_.push_back(head);
1384 }
1385 }
1386 *result = bfs_queue_;
1387}
1388
1389template <typename Graph, typename ArcFlowT, typename FlowSumT>
1391 FlowModelProto model;
1392 model.set_problem_type(FlowModelProto::MAX_FLOW);
1393 for (int n = 0; n < graph_->num_nodes(); ++n) {
1394 FlowNodeProto* node = model.add_nodes();
1395 node->set_id(n);
1396 if (n == source_) node->set_supply(1);
1397 if (n == sink_) node->set_supply(-1);
1398 }
1399 for (int a = 0; a < graph_->num_arcs(); ++a) {
1400 FlowArcProto* arc = model.add_arcs();
1401 arc->set_tail(graph_->Tail(a));
1402 arc->set_head(graph_->Head(a));
1403 arc->set_capacity(Capacity(a));
1404 }
1405 return model;
1406}
1407
1408} // namespace operations_research
1409
1410#endif // OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_
void InitializePreflow()
Initializes the preflow to a state that enables to run Refine.
const Graph * graph_
A pointer to the graph passed as argument.
std::unique_ptr< NodeHeight[]> node_potential_
PriorityQueueWithRestrictedPush< NodeIndex, NodeHeight > active_node_by_height_
bool IsAdmissible(NodeIndex tail, ArcIndex arc, const NodeHeight *potentials) const
Returns true if arc is admissible.
std::unique_ptr< ArcIndex[]> first_admissible_arc_
An array representing the first admissible arc for each node in graph_.
void GetSourceSideMinCut(std::vector< NodeIndex > *result)
void PushFlow(ArcFlowType flow, NodeIndex tail, ArcIndex arc)
NodeIndex GetSourceNodeIndex() const
Returns the index of the node corresponding to the source of the network.
ArcIndex Opposite(ArcIndex arc) const
bool CheckRelabelPrecondition(NodeIndex node) const
StatsGroup stats_
Statistics about this class.
NodeIndex sink_
The index of the sink node in graph_.
void InitializeActiveNodeContainer()
Initializes the container active_nodes_.
std::unique_ptr< ArcFlowType[]> initial_capacity_
GenericMaxFlow & operator=(const GenericMaxFlow &)=delete
const util::ReverseArcStaticGraph< NodeIndex, ArcIndex > * graph() const
FlowSumType Flow(ArcIndex arc) const
NodeIndex Tail(ArcIndex arc) const
std::vector< NodeIndex > bfs_queue_
FlowModelProto CreateFlowModel()
Returns the protocol buffer representation of the current problem.
std::unique_ptr< FlowSumType[]> node_excess_
void PushActiveNode(const NodeIndex &node)
Push element to the active node container.
bool Solve()
Returns true if a maximum flow was solved.
void GetSinkSideMinCut(std::vector< NodeIndex > *result)
NodeIndex source_
The index of the source node in graph_.
void ComputeReachableNodes(NodeIndex start, std::vector< NodeIndex > *result)
void PushAsMuchFlowAsPossible(NodeIndex tail, ArcIndex arc, FlowSumType *node_excess)
NodeIndex Head(ArcIndex arc) const
Handy member functions to make the code more compact.
ArcFlowType Capacity(ArcIndex arc) const
Returns the initial capacity of an arc.
GenericMaxFlow(const Graph *graph, NodeIndex source, NodeIndex sink)
void Refine()
Performs optimization step.
bool IsActive(NodeIndex node) const
static constexpr FlowSumType kMaxFlowSum
Maximum manageable flow.
std::string DebugString(absl::string_view context, ArcIndex arc) const
NodeIndex GetSinkNodeIndex() const
Returns the index of the node corresponding to the sink of the network.
ZVector< ArcFlowType > residual_arc_capacity_
NodeIndex GetAndRemoveFirstActiveNode()
Get the first element from the active node container.
FlowSumType GetOptimalFlow() const
Returns the total flow found by the algorithm.
bool IsEmptyActiveNodeContainer()
Check the emptiness of the container.
Status status_
The status of the problem.
void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity)
GenericMaxFlow(const GenericMaxFlow &)=delete
This type is neither copyable nor movable.
void Push(Element element, IntegerPriority priority)
PriorityQueueWithRestrictedPush & operator=(const PriorityQueueWithRestrictedPush &)=delete
PriorityQueueWithRestrictedPush(const PriorityQueueWithRestrictedPush &)=delete
This type is neither copyable nor movable.
Base class to print a nice summary of a group of statistics.
Definition stats.h:131
static constexpr int32_t kNilArc
Definition graph.h:293
static constexpr bool kHasNegativeReverseArcs
Definition graph.h:215
In SWIG mode, we don't want anything besides these top-level includes.
util::ReverseArcStaticGraph Graph
Type of graph to use.
#define IF_STATS_ENABLED(instructions)
Definition stats.h:414
#define SCOPED_TIME_STAT(stats)
Definition stats.h:421