OR-Tools  8.2
max_flow.cc
Go to the documentation of this file.
1// Copyright 2010-2018 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
15
16#include <algorithm>
17
18#include "absl/memory/memory.h"
19#include "absl/strings/str_format.h"
20#include "ortools/graph/graph.h"
22
23namespace operations_research {
24
25SimpleMaxFlow::SimpleMaxFlow() : num_nodes_(0) {}
26
29 const ArcIndex num_arcs = arc_tail_.size();
30 num_nodes_ = std::max(num_nodes_, tail + 1);
31 num_nodes_ = std::max(num_nodes_, head + 1);
32 arc_tail_.push_back(tail);
33 arc_head_.push_back(head);
34 arc_capacity_.push_back(capacity);
35 return num_arcs;
36}
37
38NodeIndex SimpleMaxFlow::NumNodes() const { return num_nodes_; }
39
40ArcIndex SimpleMaxFlow::NumArcs() const { return arc_tail_.size(); }
41
42NodeIndex SimpleMaxFlow::Tail(ArcIndex arc) const { return arc_tail_[arc]; }
43
44NodeIndex SimpleMaxFlow::Head(ArcIndex arc) const { return arc_head_[arc]; }
45
47 return arc_capacity_[arc];
48}
49
51 arc_capacity_[arc] = capacity;
52}
53
55 const ArcIndex num_arcs = arc_capacity_.size();
56 arc_flow_.assign(num_arcs, 0);
57 underlying_max_flow_.reset();
58 underlying_graph_.reset();
59 optimal_flow_ = 0;
60 if (source == sink || source < 0 || sink < 0) {
61 return BAD_INPUT;
62 }
63 if (source >= num_nodes_ || sink >= num_nodes_) {
64 return OPTIMAL;
65 }
66 underlying_graph_ = absl::make_unique<Graph>(num_nodes_, num_arcs);
67 underlying_graph_->AddNode(source);
68 underlying_graph_->AddNode(sink);
69 for (int arc = 0; arc < num_arcs; ++arc) {
70 underlying_graph_->AddArc(arc_tail_[arc], arc_head_[arc]);
71 }
72 underlying_graph_->Build(&arc_permutation_);
73 underlying_max_flow_ = absl::make_unique<GenericMaxFlow<Graph>>(
74 underlying_graph_.get(), source, sink);
75 for (ArcIndex arc = 0; arc < num_arcs; ++arc) {
76 ArcIndex permuted_arc =
77 arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
78 underlying_max_flow_->SetArcCapacity(permuted_arc, arc_capacity_[arc]);
79 }
80 if (underlying_max_flow_->Solve()) {
81 optimal_flow_ = underlying_max_flow_->GetOptimalFlow();
82 for (ArcIndex arc = 0; arc < num_arcs; ++arc) {
83 ArcIndex permuted_arc =
84 arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
85 arc_flow_[arc] = underlying_max_flow_->Flow(permuted_arc);
86 }
87 }
88 // Translate the GenericMaxFlow::Status. It is different because NOT_SOLVED
89 // does not make sense in the simple api.
90 switch (underlying_max_flow_->status()) {
92 return BAD_RESULT;
94 return OPTIMAL;
96 return POSSIBLE_OVERFLOW;
98 return BAD_INPUT;
100 return BAD_RESULT;
101 }
102 return BAD_RESULT;
103}
104
105FlowQuantity SimpleMaxFlow::OptimalFlow() const { return optimal_flow_; }
106
107FlowQuantity SimpleMaxFlow::Flow(ArcIndex arc) const { return arc_flow_[arc]; }
108
109void SimpleMaxFlow::GetSourceSideMinCut(std::vector<NodeIndex>* result) {
110 if (underlying_max_flow_ == nullptr) return;
111 underlying_max_flow_->GetSourceSideMinCut(result);
112}
113
114void SimpleMaxFlow::GetSinkSideMinCut(std::vector<NodeIndex>* result) {
115 if (underlying_max_flow_ == nullptr) return;
116 underlying_max_flow_->GetSinkSideMinCut(result);
117}
118
120 if (underlying_max_flow_ == nullptr) return FlowModel();
121 return underlying_max_flow_->CreateFlowModel();
122}
123
124template <typename Graph>
126 NodeIndex sink)
127 : graph_(graph),
128 node_excess_(),
129 node_potential_(),
130 residual_arc_capacity_(),
131 first_admissible_arc_(),
132 active_nodes_(),
133 source_(source),
134 sink_(sink),
135 use_global_update_(true),
136 use_two_phase_algorithm_(true),
137 process_node_by_height_(true),
138 check_input_(true),
139 check_result_(true),
140 stats_("MaxFlow") {
142 DCHECK(graph->IsNodeValid(source));
143 DCHECK(graph->IsNodeValid(sink));
144 const NodeIndex max_num_nodes = Graphs<Graph>::NodeReservation(*graph_);
145 if (max_num_nodes > 0) {
146 node_excess_.Reserve(0, max_num_nodes - 1);
148 node_potential_.Reserve(0, max_num_nodes - 1);
150 first_admissible_arc_.Reserve(0, max_num_nodes - 1);
151 first_admissible_arc_.SetAll(Graph::kNilArc);
152 bfs_queue_.reserve(max_num_nodes);
153 active_nodes_.reserve(max_num_nodes);
154 }
155 const ArcIndex max_num_arcs = Graphs<Graph>::ArcReservation(*graph_);
156 if (max_num_arcs > 0) {
157 residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1);
159 }
160}
161
162template <typename Graph>
164 SCOPED_TIME_STAT(&stats_);
165 bool ok = true;
166 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
167 if (residual_arc_capacity_[arc] < 0) {
168 ok = false;
169 }
170 }
171 return ok;
172}
173
174template <typename Graph>
176 FlowQuantity new_capacity) {
177 SCOPED_TIME_STAT(&stats_);
178 DCHECK_LE(0, new_capacity);
179 DCHECK(IsArcDirect(arc));
180 const FlowQuantity free_capacity = residual_arc_capacity_[arc];
181 const FlowQuantity capacity_delta = new_capacity - Capacity(arc);
182 if (capacity_delta == 0) {
183 return; // Nothing to do.
184 }
185 status_ = NOT_SOLVED;
186 if (free_capacity + capacity_delta >= 0) {
187 // The above condition is true if one of the two conditions is true:
188 // 1/ (capacity_delta > 0), meaning we are increasing the capacity
189 // 2/ (capacity_delta < 0 && free_capacity + capacity_delta >= 0)
190 // meaning we are reducing the capacity, but that the capacity
191 // reduction is not larger than the free capacity.
192 DCHECK((capacity_delta > 0) ||
193 (capacity_delta < 0 && free_capacity + capacity_delta >= 0));
194 residual_arc_capacity_.Set(arc, free_capacity + capacity_delta);
195 DCHECK_LE(0, residual_arc_capacity_[arc]);
196 } else {
197 // Note that this breaks the preflow invariants but it is currently not an
198 // issue since we restart from scratch on each Solve() and we set the status
199 // to NOT_SOLVED.
200 //
201 // TODO(user): The easiest is probably to allow negative node excess in
202 // other places than the source, but the current implementation does not
203 // deal with this.
204 SetCapacityAndClearFlow(arc, new_capacity);
205 }
206}
207
208template <typename Graph>
210 SCOPED_TIME_STAT(&stats_);
211 DCHECK(IsArcValid(arc));
212 DCHECK_GE(new_flow, 0);
213 const FlowQuantity capacity = Capacity(arc);
214 DCHECK_GE(capacity, new_flow);
215
216 // Note that this breaks the preflow invariants but it is currently not an
217 // issue since we restart from scratch on each Solve() and we set the status
218 // to NOT_SOLVED.
219 residual_arc_capacity_.Set(Opposite(arc), -new_flow);
220 residual_arc_capacity_.Set(arc, capacity - new_flow);
221 status_ = NOT_SOLVED;
222}
223
224template <typename Graph>
226 std::vector<NodeIndex>* result) {
227 ComputeReachableNodes<false>(source_, result);
228}
229
230template <typename Graph>
231void GenericMaxFlow<Graph>::GetSinkSideMinCut(std::vector<NodeIndex>* result) {
232 ComputeReachableNodes<true>(sink_, result);
233}
234
235template <typename Graph>
237 SCOPED_TIME_STAT(&stats_);
238 bool ok = true;
239 if (node_excess_[source_] != -node_excess_[sink_]) {
240 LOG(DFATAL) << "-node_excess_[source_] = " << -node_excess_[source_]
241 << " != node_excess_[sink_] = " << node_excess_[sink_];
242 ok = false;
243 }
244 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
245 if (node != source_ && node != sink_) {
246 if (node_excess_[node] != 0) {
247 LOG(DFATAL) << "node_excess_[" << node << "] = " << node_excess_[node]
248 << " != 0";
249 ok = false;
250 }
251 }
252 }
253 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
254 const ArcIndex opposite = Opposite(arc);
255 const FlowQuantity direct_capacity = residual_arc_capacity_[arc];
256 const FlowQuantity opposite_capacity = residual_arc_capacity_[opposite];
257 if (direct_capacity < 0) {
258 LOG(DFATAL) << "residual_arc_capacity_[" << arc
259 << "] = " << direct_capacity << " < 0";
260 ok = false;
261 }
262 if (opposite_capacity < 0) {
263 LOG(DFATAL) << "residual_arc_capacity_[" << opposite
264 << "] = " << opposite_capacity << " < 0";
265 ok = false;
266 }
267 // The initial capacity of the direct arcs is non-negative.
268 if (direct_capacity + opposite_capacity < 0) {
269 LOG(DFATAL) << "initial capacity [" << arc
270 << "] = " << direct_capacity + opposite_capacity << " < 0";
271 ok = false;
272 }
273 }
274 return ok;
275}
276
277template <typename Graph>
279 SCOPED_TIME_STAT(&stats_);
280
281 // We simply compute the reachability from the source in the residual graph.
282 const NodeIndex num_nodes = graph_->num_nodes();
283 std::vector<bool> is_reached(num_nodes, false);
284 std::vector<NodeIndex> to_process;
285
286 to_process.push_back(source_);
287 is_reached[source_] = true;
288 while (!to_process.empty()) {
289 const NodeIndex node = to_process.back();
290 to_process.pop_back();
291 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
292 it.Next()) {
293 const ArcIndex arc = it.Index();
294 if (residual_arc_capacity_[arc] > 0) {
295 const NodeIndex head = graph_->Head(arc);
296 if (!is_reached[head]) {
297 is_reached[head] = true;
298 to_process.push_back(head);
299 }
300 }
301 }
302 }
303 return is_reached[sink_];
304}
305
306template <typename Graph>
308 DCHECK(IsActive(node));
309 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
310 it.Next()) {
311 const ArcIndex arc = it.Index();
312 DCHECK(!IsAdmissible(arc)) << DebugString("CheckRelabelPrecondition:", arc);
313 }
314 return true;
315}
316
317template <typename Graph>
318std::string GenericMaxFlow<Graph>::DebugString(const std::string& context,
319 ArcIndex arc) const {
320 const NodeIndex tail = Tail(arc);
321 const NodeIndex head = Head(arc);
322 return absl::StrFormat(
323 "%s Arc %d, from %d to %d, "
324 "Capacity = %d, Residual capacity = %d, "
325 "Flow = residual capacity for reverse arc = %d, "
326 "Height(tail) = %d, Height(head) = %d, "
327 "Excess(tail) = %d, Excess(head) = %d",
328 context, arc, tail, head, Capacity(arc), residual_arc_capacity_[arc],
329 Flow(arc), node_potential_[tail], node_potential_[head],
330 node_excess_[tail], node_excess_[head]);
331}
332
333template <typename Graph>
335 status_ = NOT_SOLVED;
336 if (check_input_ && !CheckInputConsistency()) {
337 status_ = BAD_INPUT;
338 return false;
339 }
340 InitializePreflow();
341
342 // Deal with the case when source_ or sink_ is not inside graph_.
343 // Since they are both specified independently of the graph, we do need to
344 // take care of this corner case.
345 const NodeIndex num_nodes = graph_->num_nodes();
346 if (sink_ >= num_nodes || source_ >= num_nodes) {
347 // Behave like a normal graph where source_ and sink_ are disconnected.
348 // Note that the arc flow is set to 0 by InitializePreflow().
349 status_ = OPTIMAL;
350 return true;
351 }
352 if (use_global_update_) {
353 RefineWithGlobalUpdate();
354 } else {
355 Refine();
356 }
357 if (check_result_) {
358 if (!CheckResult()) {
359 status_ = BAD_RESULT;
360 return false;
361 }
362 if (GetOptimalFlow() < kMaxFlowQuantity && AugmentingPathExists()) {
363 LOG(ERROR) << "The algorithm terminated, but the flow is not maximal!";
364 status_ = BAD_RESULT;
365 return false;
366 }
367 }
368 DCHECK_EQ(node_excess_[sink_], -node_excess_[source_]);
369 status_ = OPTIMAL;
370 if (GetOptimalFlow() == kMaxFlowQuantity && AugmentingPathExists()) {
371 // In this case, we are sure that the flow is > kMaxFlowQuantity.
372 status_ = INT_OVERFLOW;
373 }
374 IF_STATS_ENABLED(VLOG(1) << stats_.StatString());
375 return true;
376}
377
378template <typename Graph>
380 SCOPED_TIME_STAT(&stats_);
381 // InitializePreflow() clears the whole flow that could have been computed
382 // by a previous Solve(). This is not optimal in terms of complexity.
383 // TODO(user): find a way to make the re-solving incremental (not an obvious
384 // task, and there has not been a lot of literature on the subject.)
385 node_excess_.SetAll(0);
386 const ArcIndex num_arcs = graph_->num_arcs();
387 for (ArcIndex arc = 0; arc < num_arcs; ++arc) {
388 SetCapacityAndClearFlow(arc, Capacity(arc));
389 }
390
391 // All the initial heights are zero except for the source whose height is
392 // equal to the number of nodes and will never change during the algorithm.
393 node_potential_.SetAll(0);
394 node_potential_.Set(source_, graph_->num_nodes());
396 // Initially no arcs are admissible except maybe the one leaving the source,
397 // but we treat the source in a special way, see
398 // SaturateOutgoingArcsFromSource().
399 const NodeIndex num_nodes = graph_->num_nodes();
400 for (NodeIndex node = 0; node < num_nodes; ++node) {
401 first_admissible_arc_[node] = Graph::kNilArc;
402 }
403}
405// Note(user): Calling this function will break the property on the node
406// potentials because of the way we cancel flow on cycle. However, we only call
407// that at the end of the algorithm, or just before a GlobalUpdate() that will
408// restore the precondition on the node potentials.
409template <typename Graph>
411 SCOPED_TIME_STAT(&stats_);
412 const NodeIndex num_nodes = graph_->num_nodes();
413
414 // We implement a variation of Tarjan's strongly connected component algorithm
415 // to detect cycles published in: Tarjan, R. E. (1972), "Depth-first search
416 // and linear graph algorithms", SIAM Journal on Computing. A description can
417 // also be found in wikipedia.
418
419 // Stored nodes are settled nodes already stored in the
420 // reverse_topological_order (except the sink_ that we do not actually store).
421 std::vector<bool> stored(num_nodes, false);
422 stored[sink_] = true;
423
424 // The visited nodes that are not yet stored are all the nodes from the
425 // source_ to the current node in the current dfs branch.
426 std::vector<bool> visited(num_nodes, false);
427 visited[sink_] = true;
428
429 // Stack of arcs to explore in the dfs search.
430 // The current node is Head(arc_stack.back()).
431 std::vector<ArcIndex> arc_stack;
432
433 // Increasing list of indices into the arc_stack that correspond to the list
434 // of arcs in the current dfs branch from the source_ to the current node.
435 std::vector<int> index_branch;
436
437 // Node in reverse_topological_order in the final dfs tree.
438 std::vector<NodeIndex> reverse_topological_order;
439
440 // We start by pushing all the outgoing arcs from the source on the stack to
441 // avoid special conditions in the code. As a result, source_ will not be
442 // stored in reverse_topological_order, and this is what we want.
443 for (OutgoingArcIterator it(*graph_, source_); it.Ok(); it.Next()) {
444 const ArcIndex arc = it.Index();
445 const FlowQuantity flow = Flow(arc);
446 if (flow > 0) {
447 arc_stack.push_back(arc);
448 }
449 }
450 visited[source_] = true;
451
452 // Start the dfs on the subgraph formed by the direct arcs with positive flow.
453 while (!arc_stack.empty()) {
454 const NodeIndex node = Head(arc_stack.back());
455
456 // If the node is visited, it means we have explored all its arcs and we
457 // have just backtracked in the dfs. Store it if it is not already stored
458 // and process the next arc on the stack.
459 if (visited[node]) {
460 if (!stored[node]) {
461 stored[node] = true;
462 reverse_topological_order.push_back(node);
463 DCHECK(!index_branch.empty());
464 index_branch.pop_back();
465 }
466 arc_stack.pop_back();
467 continue;
468 }
469
470 // The node is a new unexplored node, add all its outgoing arcs with
471 // positive flow to the stack and go deeper in the dfs.
472 DCHECK(!stored[node]);
473 DCHECK(index_branch.empty() ||
474 (arc_stack.size() - 1 > index_branch.back()));
475 visited[node] = true;
476 index_branch.push_back(arc_stack.size() - 1);
477
478 for (OutgoingArcIterator it(*graph_, node); it.Ok(); it.Next()) {
479 const ArcIndex arc = it.Index();
480 const FlowQuantity flow = Flow(arc);
481 const NodeIndex head = Head(arc);
482 if (flow > 0 && !stored[head]) {
483 if (!visited[head]) {
484 arc_stack.push_back(arc);
485 } else {
486 // There is a cycle.
487 // Find the first index to consider,
488 // arc_stack[index_branch[cycle_begin]] will be the first arc on the
489 // cycle.
490 int cycle_begin = index_branch.size();
491 while (cycle_begin > 0 &&
492 Head(arc_stack[index_branch[cycle_begin - 1]]) != head) {
493 --cycle_begin;
495
496 // Compute the maximum flow that can be canceled on the cycle and the
497 // min index such that arc_stack[index_branch[i]] will be saturated.
498 FlowQuantity max_flow = flow;
499 int first_saturated_index = index_branch.size();
500 for (int i = index_branch.size() - 1; i >= cycle_begin; --i) {
501 const ArcIndex arc_on_cycle = arc_stack[index_branch[i]];
502 if (Flow(arc_on_cycle) <= max_flow) {
503 max_flow = Flow(arc_on_cycle);
504 first_saturated_index = i;
505 }
506 }
507
508 // This is just here for a DCHECK() below.
509 const FlowQuantity excess = node_excess_[head];
510
511 // Cancel the flow on the cycle, and set visited[node] = false for
512 // the node that will be backtracked over.
513 PushFlow(-max_flow, arc);
514 for (int i = index_branch.size() - 1; i >= cycle_begin; --i) {
515 const ArcIndex arc_on_cycle = arc_stack[index_branch[i]];
516 PushFlow(-max_flow, arc_on_cycle);
517 if (i >= first_saturated_index) {
518 DCHECK(visited[Head(arc_on_cycle)]);
519 visited[Head(arc_on_cycle)] = false;
520 } else {
521 DCHECK_GT(Flow(arc_on_cycle), 0);
523 }
524
525 // This is a simple check that the flow was pushed properly.
526 DCHECK_EQ(excess, node_excess_[head]);
527
528 // Backtrack the dfs just before index_branch[first_saturated_index].
529 // If the current node is still active, there is nothing to do.
530 if (first_saturated_index < index_branch.size()) {
531 arc_stack.resize(index_branch[first_saturated_index]);
532 index_branch.resize(first_saturated_index);
533
534 // We backtracked over the current node, so there is no need to
535 // continue looping over its arcs.
536 break;
537 }
538 }
539 }
540 }
542 DCHECK(arc_stack.empty());
543 DCHECK(index_branch.empty());
544
545 // Return the flow to the sink. Note that the sink_ and the source_ are not
546 // stored in reverse_topological_order.
547 for (int i = 0; i < reverse_topological_order.size(); i++) {
548 const NodeIndex node = reverse_topological_order[i];
549 if (node_excess_[node] == 0) continue;
550 for (IncomingArcIterator it(*graph_, node); it.Ok(); it.Next()) {
551 const ArcIndex opposite_arc = Opposite(it.Index());
552 if (residual_arc_capacity_[opposite_arc] > 0) {
553 const FlowQuantity flow =
554 std::min(node_excess_[node], residual_arc_capacity_[opposite_arc]);
555 PushFlow(flow, opposite_arc);
556 if (node_excess_[node] == 0) break;
557 }
558 }
559 DCHECK_EQ(0, node_excess_[node]);
560 }
561 DCHECK_EQ(-node_excess_[source_], node_excess_[sink_]);
562}
563
564template <typename Graph>
566 SCOPED_TIME_STAT(&stats_);
567 bfs_queue_.clear();
568 int queue_index = 0;
569 const NodeIndex num_nodes = graph_->num_nodes();
570 node_in_bfs_queue_.assign(num_nodes, false);
571 node_in_bfs_queue_[sink_] = true;
572 node_in_bfs_queue_[source_] = true;
573
574 // We do two BFS in the reverse residual graph, one from the sink and one from
575 // the source. Because all the arcs from the source are saturated (except in
576 // presence of integer overflow), the source cannot reach the sink in the
577 // residual graph. However, we still want to relabel all the nodes that cannot
578 // reach the sink but can reach the source (because if they have excess, we
579 // need to push it back to the source).
580 //
581 // Note that the second pass is not needed here if we use a two-pass algorithm
582 // to return the flow to the source after we found the min cut.
583 const int num_passes = use_two_phase_algorithm_ ? 1 : 2;
584 for (int pass = 0; pass < num_passes; ++pass) {
585 if (pass == 0) {
586 bfs_queue_.push_back(sink_);
587 } else {
588 bfs_queue_.push_back(source_);
589 }
590
591 while (queue_index != bfs_queue_.size()) {
592 const NodeIndex node = bfs_queue_[queue_index];
593 ++queue_index;
594 const NodeIndex candidate_distance = node_potential_[node] + 1;
595 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
596 it.Next()) {
597 const ArcIndex arc = it.Index();
598 const NodeIndex head = Head(arc);
599
600 // Skip the arc if the height of head was already set to the correct
601 // value (Remember we are doing reverse BFS).
602 if (node_in_bfs_queue_[head]) continue;
603
604 // TODO(user): By using more memory we can speed this up quite a bit by
605 // avoiding to take the opposite arc here, too options:
606 // - if (residual_arc_capacity_[arc] != arc_capacity_[arc])
607 // - if (opposite_arc_is_admissible_[arc]) // need updates.
608 // Experiment with the first option shows more than 10% gain on this
609 // function running time, which is the bottleneck on many instances.
610 const ArcIndex opposite_arc = Opposite(arc);
611 if (residual_arc_capacity_[opposite_arc] > 0) {
612 // Note(user): We used to have a DCHECK_GE(candidate_distance,
613 // node_potential_[head]); which is always true except in the case
614 // where we can push more than kMaxFlowQuantity out of the source. The
615 // problem comes from the fact that in this case, we call
616 // PushFlowExcessBackToSource() in the middle of the algorithm. The
617 // later call will break the properties of the node potential. Note
618 // however, that this function will recompute a good node potential
619 // for all the nodes and thus fix the issue.
620
621 // If head is active, we can steal some or all of its excess.
622 // This brings a huge gain on some problems.
623 // Note(user): I haven't seen this anywhere in the literature.
624 // TODO(user): Investigate more and maybe write a publication :)
625 if (node_excess_[head] > 0) {
626 const FlowQuantity flow = std::min(
627 node_excess_[head], residual_arc_capacity_[opposite_arc]);
628 PushFlow(flow, opposite_arc);
629
630 // If the arc became saturated, it is no longer in the residual
631 // graph, so we do not need to consider head at this time.
632 if (residual_arc_capacity_[opposite_arc] == 0) continue;
633 }
634
635 // Note that there is no need to touch first_admissible_arc_[node]
636 // because of the relaxed Relabel() we use.
637 node_potential_[head] = candidate_distance;
638 node_in_bfs_queue_[head] = true;
639 bfs_queue_.push_back(head);
640 }
641 }
642 }
643 }
644
645 // At the end of the search, some nodes may not be in the bfs_queue_. Such
646 // nodes cannot reach the sink_ or source_ in the residual graph, so there is
647 // no point trying to push flow toward them. We obtain this effect by setting
648 // their height to something unreachable.
649 //
650 // Note that this also prevents cycling due to our anti-overflow procedure.
651 // For instance, suppose there is an edge s -> n outgoing from the source. If
652 // node n has no other connection and some excess, we will push the flow back
653 // to the source, but if we don't update the height of n
654 // SaturateOutgoingArcsFromSource() will push the flow to n again.
655 // TODO(user): This is another argument for another anti-overflow algorithm.
656 for (NodeIndex node = 0; node < num_nodes; ++node) {
657 if (!node_in_bfs_queue_[node]) {
658 node_potential_[node] = 2 * num_nodes - 1;
659 }
660 }
661
662 // Reset the active nodes. Doing it like this pushes the nodes in increasing
663 // order of height. Note that bfs_queue_[0] is the sink_ so we skip it.
664 DCHECK(IsEmptyActiveNodeContainer());
665 for (int i = 1; i < bfs_queue_.size(); ++i) {
666 const NodeIndex node = bfs_queue_[i];
667 if (node_excess_[node] > 0) {
668 DCHECK(IsActive(node));
669 PushActiveNode(node);
670 }
671 }
672}
673
674template <typename Graph>
676 SCOPED_TIME_STAT(&stats_);
677 const NodeIndex num_nodes = graph_->num_nodes();
678
679 // If sink_ or source_ already have kMaxFlowQuantity, then there is no
680 // point pushing more flow since it will cause an integer overflow.
681 if (node_excess_[sink_] == kMaxFlowQuantity) return false;
682 if (node_excess_[source_] == -kMaxFlowQuantity) return false;
683
684 bool flow_pushed = false;
685 for (OutgoingArcIterator it(*graph_, source_); it.Ok(); it.Next()) {
686 const ArcIndex arc = it.Index();
687 const FlowQuantity flow = residual_arc_capacity_[arc];
688
689 // This is a special IsAdmissible() condition for the source.
690 if (flow == 0 || node_potential_[Head(arc)] >= num_nodes) continue;
691
692 // We are careful in case the sum of the flow out of the source is greater
693 // than kMaxFlowQuantity to avoid overflow.
694 const FlowQuantity current_flow_out_of_source = -node_excess_[source_];
695 DCHECK_GE(flow, 0) << flow;
696 DCHECK_GE(current_flow_out_of_source, 0) << current_flow_out_of_source;
697 const FlowQuantity capped_flow =
698 kMaxFlowQuantity - current_flow_out_of_source;
699 if (capped_flow < flow) {
700 // We push as much flow as we can so the current flow on the network will
701 // be kMaxFlowQuantity.
702
703 // Since at the beginning of the function, current_flow_out_of_source
704 // was different from kMaxFlowQuantity, we are sure to have pushed some
705 // flow before if capped_flow is 0.
706 if (capped_flow == 0) return true;
707 PushFlow(capped_flow, arc);
708 return true;
709 }
710 PushFlow(flow, arc);
711 flow_pushed = true;
712 }
713 DCHECK_LE(node_excess_[source_], 0);
714 return flow_pushed;
715}
716
717template <typename Graph>
719 SCOPED_TIME_STAT(&stats_);
720 // TODO(user): Do not allow a zero flow after fixing the UniformMaxFlow code.
721 DCHECK_GE(residual_arc_capacity_[Opposite(arc)] + flow, 0);
722 DCHECK_GE(residual_arc_capacity_[arc] - flow, 0);
723
724 // node_excess_ should be always greater than or equal to 0 except for the
725 // source where it should always be smaller than or equal to 0. Note however
726 // that we cannot check this because when we cancel the flow on a cycle in
727 // PushFlowExcessBackToSource(), we may break this invariant during the
728 // operation even if it is still valid at the end.
729
730 // Update the residual capacity of the arc and its opposite arc.
731 residual_arc_capacity_[arc] -= flow;
732 residual_arc_capacity_[Opposite(arc)] += flow;
733
734 // Update the excesses at the tail and head of the arc.
735 node_excess_[Tail(arc)] -= flow;
736 node_excess_[Head(arc)] += flow;
737}
738
739template <typename Graph>
741 SCOPED_TIME_STAT(&stats_);
742 DCHECK(IsEmptyActiveNodeContainer());
743 const NodeIndex num_nodes = graph_->num_nodes();
744 for (NodeIndex node = 0; node < num_nodes; ++node) {
745 if (IsActive(node)) {
746 if (use_two_phase_algorithm_ && node_potential_[node] >= num_nodes) {
747 continue;
748 }
749 PushActiveNode(node);
750 }
751 }
752}
753
754template <typename Graph>
756 SCOPED_TIME_STAT(&stats_);
757 // Usually SaturateOutgoingArcsFromSource() will saturate all the arcs from
758 // the source in one go, and we will loop just once. But in case we can push
759 // more than kMaxFlowQuantity out of the source the loop is as follow:
760 // - Push up to kMaxFlowQuantity out of the source on the admissible outgoing
761 // arcs. Stop if no flow was pushed.
762 // - Compute the current max-flow. This will push some flow back to the
763 // source and render more outgoing arcs from the source not admissible.
764 //
765 // TODO(user): This may not be the most efficient algorithm if we need to loop
766 // many times. An alternative may be to handle the source like the other nodes
767 // in the algorithm, initially putting an excess of kMaxFlowQuantity on it,
768 // and making the source active like any other node with positive excess. To
769 // investigate.
770 //
771 // TODO(user): The code below is buggy when more than kMaxFlowQuantity can be
772 // pushed out of the source (i.e. when we loop more than once in the while()).
773 // This is not critical, since this code is not used in the default algorithm
774 // computation. The issue is twofold:
775 // - InitializeActiveNodeContainer() doesn't push the nodes in
776 // the correct order.
777 // - PushFlowExcessBackToSource() may break the node potential properties, and
778 // we will need a call to GlobalUpdate() to fix that.
779 while (SaturateOutgoingArcsFromSource()) {
780 DCHECK(IsEmptyActiveNodeContainer());
781 InitializeActiveNodeContainer();
782 while (!IsEmptyActiveNodeContainer()) {
783 const NodeIndex node = GetAndRemoveFirstActiveNode();
784 if (node == source_ || node == sink_) continue;
785 Discharge(node);
786 }
787 if (use_two_phase_algorithm_) {
788 PushFlowExcessBackToSource();
789 }
790 }
791}
792
793template <typename Graph>
795 SCOPED_TIME_STAT(&stats_);
796
797 // TODO(user): This should be graph_->num_nodes(), but ebert graph does not
798 // have a correct size if the highest index nodes have no arcs.
799 const NodeIndex num_nodes = Graphs<Graph>::NodeReservation(*graph_);
800 std::vector<int> skip_active_node;
801
802 while (SaturateOutgoingArcsFromSource()) {
803 int num_skipped;
804 do {
805 num_skipped = 0;
806 skip_active_node.assign(num_nodes, 0);
807 skip_active_node[sink_] = 2;
808 skip_active_node[source_] = 2;
809 GlobalUpdate();
810 while (!IsEmptyActiveNodeContainer()) {
811 const NodeIndex node = GetAndRemoveFirstActiveNode();
812 if (skip_active_node[node] > 1) {
813 if (node != sink_ && node != source_) ++num_skipped;
814 continue;
815 }
816 const NodeIndex old_height = node_potential_[node];
817 Discharge(node);
818
819 // The idea behind this is that if a node height augments by more than
820 // one, then it is likely to push flow back the way it came. This can
821 // lead to very costly loops. A bad case is: source -> n1 -> n2 and n2
822 // just recently isolated from the sink. Then n2 will push flow back to
823 // n1, and n1 to n2 and so on. The height of each node will increase by
824 // steps of two until the height of the source is reached, which can
825 // take a long time. If the chain is longer, the situation is even
826 // worse. The behavior of this heuristic is related to the Gap
827 // heuristic.
828 //
829 // Note that the global update will fix all such cases efficiently. So
830 // the idea is to discharge the active node as much as possible, and
831 // then do a global update.
832 //
833 // We skip a node when this condition was true 2 times to avoid doing a
834 // global update too frequently.
835 if (node_potential_[node] > old_height + 1) {
836 ++skip_active_node[node];
837 }
838 }
839 } while (num_skipped > 0);
840 if (use_two_phase_algorithm_) {
841 PushFlowExcessBackToSource();
842 }
843 }
844}
845
846template <typename Graph>
848 SCOPED_TIME_STAT(&stats_);
849 const NodeIndex num_nodes = graph_->num_nodes();
850 while (true) {
851 DCHECK(IsActive(node));
852 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
853 first_admissible_arc_[node]);
854 it.Ok(); it.Next()) {
855 const ArcIndex arc = it.Index();
856 if (IsAdmissible(arc)) {
857 DCHECK(IsActive(node));
858 const NodeIndex head = Head(arc);
859 if (node_excess_[head] == 0) {
860 // The push below will make the node active for sure. Note that we may
861 // push the sink_, but that is handled properly in Refine().
862 PushActiveNode(head);
863 }
864 const FlowQuantity delta =
865 std::min(node_excess_[node], residual_arc_capacity_[arc]);
866 PushFlow(delta, arc);
867 if (node_excess_[node] == 0) {
868 first_admissible_arc_[node] = arc; // arc may still be admissible.
869 return;
870 }
871 }
872 }
873 Relabel(node);
874 if (use_two_phase_algorithm_ && node_potential_[node] >= num_nodes) break;
875 }
876}
877
878template <typename Graph>
880 SCOPED_TIME_STAT(&stats_);
881 // Because we use a relaxed version, this is no longer true if the
882 // first_admissible_arc_[node] was not actually the first arc!
883 // DCHECK(CheckRelabelPrecondition(node));
885 ArcIndex first_admissible_arc = Graph::kNilArc;
886 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
887 it.Next()) {
888 const ArcIndex arc = it.Index();
889 if (residual_arc_capacity_[arc] > 0) {
890 // Update min_height only for arcs with available capacity.
891 NodeHeight head_height = node_potential_[Head(arc)];
892 if (head_height < min_height) {
893 min_height = head_height;
894 first_admissible_arc = arc;
895
896 // We found an admissible arc at the current height, just stop there.
897 // This is the true first_admissible_arc_[node].
898 if (min_height + 1 == node_potential_[node]) break;
899 }
900 }
901 }
902 DCHECK_NE(first_admissible_arc, Graph::kNilArc);
903 node_potential_[node] = min_height + 1;
904
905 // Note that after a Relabel(), the loop will continue in Discharge(), and
906 // we are sure that all the arcs before first_admissible_arc are not
907 // admissible since their height is > min_height.
908 first_admissible_arc_[node] = first_admissible_arc;
909}
910
911template <typename Graph>
913 return Graphs<Graph>::OppositeArc(*graph_, arc);
914}
915
916template <typename Graph>
918 return IsArcValid(arc) && arc >= 0;
919}
920
921template <typename Graph>
923 return Graphs<Graph>::IsArcValid(*graph_, arc);
924}
925
926template <typename Graph>
929
930template <typename Graph>
931template <bool reverse>
933 NodeIndex start, std::vector<NodeIndex>* result) {
934 // If start is not a valid node index, it can reach only itself.
935 // Note(user): This is needed because source and sink are given independently
936 // of the graph and sometimes before it is even constructed.
937 const NodeIndex num_nodes = graph_->num_nodes();
938 if (start >= num_nodes) {
939 result->clear();
940 result->push_back(start);
941 return;
942 }
943 bfs_queue_.clear();
944 node_in_bfs_queue_.assign(num_nodes, false);
945
946 int queue_index = 0;
947 bfs_queue_.push_back(start);
948 node_in_bfs_queue_[start] = true;
949 while (queue_index != bfs_queue_.size()) {
950 const NodeIndex node = bfs_queue_[queue_index];
951 ++queue_index;
952 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
953 it.Next()) {
954 const ArcIndex arc = it.Index();
955 const NodeIndex head = Head(arc);
956 if (node_in_bfs_queue_[head]) continue;
957 if (residual_arc_capacity_[reverse ? Opposite(arc) : arc] == 0) continue;
958 node_in_bfs_queue_[head] = true;
959 bfs_queue_.push_back(head);
960 }
961 }
962 *result = bfs_queue_;
963}
964
965template <typename Graph>
967 FlowModel model;
968 model.set_problem_type(FlowModel::MAX_FLOW);
969 for (int n = 0; n < graph_->num_nodes(); ++n) {
970 Node* node = model.add_node();
971 node->set_id(n);
972 if (n == source_) node->set_supply(1);
973 if (n == sink_) node->set_supply(-1);
974 }
975 for (int a = 0; a < graph_->num_arcs(); ++a) {
976 Arc* arc = model.add_arc();
977 arc->set_tail_node_id(graph_->Tail(a));
978 arc->set_head_node_id(graph_->Head(a));
979 arc->set_capacity(Capacity(a));
980 }
981 return model;
982}
983
984// Explicit instanciations that can be used by a client.
985//
986// TODO(user): moves this code out of a .cc file and include it at the end of
987// the header so it can work with any graph implementation ?
988template class GenericMaxFlow<StarGraph>;
992
993} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK_GT(val1, val2)
Definition: base/logging.h:890
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
#define VLOG(verboselevel)
Definition: base/logging.h:978
Graph::OutgoingArcIterator OutgoingArcIterator
Definition: max_flow.h:319
void Relabel(NodeIndex node)
Definition: max_flow.cc:879
std::vector< NodeIndex > active_nodes_
Definition: max_flow.h:590
void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity)
Definition: max_flow.cc:175
std::string DebugString(const std::string &context, ArcIndex arc) const
Definition: max_flow.cc:318
Graph::OutgoingOrOppositeIncomingArcIterator OutgoingOrOppositeIncomingArcIterator
Definition: max_flow.h:321
bool CheckRelabelPrecondition(NodeIndex node) const
Definition: max_flow.cc:307
std::vector< NodeIndex > bfs_queue_
Definition: max_flow.h:611
void SetArcFlow(ArcIndex arc, FlowQuantity new_flow)
Definition: max_flow.cc:209
void GetSourceSideMinCut(std::vector< NodeIndex > *result)
Definition: max_flow.cc:225
const Graph * graph() const
Definition: max_flow.h:338
void PushFlow(FlowQuantity flow, ArcIndex arc)
Definition: max_flow.cc:718
void ComputeReachableNodes(NodeIndex start, std::vector< NodeIndex > *result)
Definition: max_flow.cc:932
bool IsArcValid(ArcIndex arc) const
Definition: max_flow.cc:922
GenericMaxFlow(const Graph *graph, NodeIndex source, NodeIndex sink)
Definition: max_flow.cc:125
void GetSinkSideMinCut(std::vector< NodeIndex > *result)
Definition: max_flow.cc:231
ArcIndex Opposite(ArcIndex arc) const
Definition: max_flow.cc:912
bool IsArcDirect(ArcIndex arc) const
Definition: max_flow.cc:917
void Discharge(NodeIndex node)
Definition: max_flow.cc:847
FlowQuantity Flow(ArcIndex arc) const
Definition: max_flow.cc:107
void GetSourceSideMinCut(std::vector< NodeIndex > *result)
Definition: max_flow.cc:109
Status Solve(NodeIndex source, NodeIndex sink)
Definition: max_flow.cc:54
FlowQuantity OptimalFlow() const
Definition: max_flow.cc:105
ArcIndex AddArcWithCapacity(NodeIndex tail, NodeIndex head, FlowQuantity capacity)
Definition: max_flow.cc:27
FlowQuantity Capacity(ArcIndex arc) const
Definition: max_flow.cc:46
NodeIndex Head(ArcIndex arc) const
Definition: max_flow.cc:44
NodeIndex Tail(ArcIndex arc) const
Definition: max_flow.cc:42
void SetArcCapacity(ArcIndex arc, FlowQuantity capacity)
Definition: max_flow.cc:50
void GetSinkSideMinCut(std::vector< NodeIndex > *result)
Definition: max_flow.cc:114
bool Reserve(int64 new_min_index, int64 new_max_index)
Definition: zvector.h:98
GRBmodel * model
GurobiMPCallbackContext * context
const int ERROR
Definition: log_severity.h:32
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
std::pair< int64, int64 > Arc
Definition: search.cc:3361
ListGraph Graph
Definition: graph.h:2360
int64 delta
Definition: resource.cc:1684
int64 tail
int64 head
int64 capacity
#define IF_STATS_ENABLED(instructions)
Definition: stats.h:437
#define SCOPED_TIME_STAT(stats)
Definition: stats.h:438
static ArcIndex ArcReservation(const Graph &graph)
Definition: graphs.h:39
static NodeIndex NodeReservation(const Graph &graph)
Definition: graphs.h:36
static bool IsArcValid(const Graph &graph, ArcIndex arc)
Definition: graphs.h:33
static ArcIndex OppositeArc(const Graph &graph, ArcIndex arc)
Definition: graphs.h:30