18 #include "absl/memory/memory.h"
19 #include "absl/strings/str_format.h"
29 const ArcIndex num_arcs = arc_tail_.size();
32 arc_tail_.push_back(
tail);
33 arc_head_.push_back(
head);
47 return arc_capacity_[arc];
55 const ArcIndex num_arcs = arc_capacity_.size();
56 arc_flow_.assign(num_arcs, 0);
57 underlying_max_flow_.reset();
58 underlying_graph_.reset();
60 if (source == sink || source < 0 || sink < 0) {
63 if (source >= num_nodes_ || sink >= num_nodes_) {
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]);
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) {
77 arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
78 underlying_max_flow_->SetArcCapacity(permuted_arc, arc_capacity_[arc]);
80 if (underlying_max_flow_->Solve()) {
81 optimal_flow_ = underlying_max_flow_->GetOptimalFlow();
82 for (
ArcIndex arc = 0; arc < num_arcs; ++arc) {
84 arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
85 arc_flow_[arc] = underlying_max_flow_->Flow(permuted_arc);
90 switch (underlying_max_flow_->status()) {
110 if (underlying_max_flow_ ==
nullptr)
return;
111 underlying_max_flow_->GetSourceSideMinCut(result);
115 if (underlying_max_flow_ ==
nullptr)
return;
116 underlying_max_flow_->GetSinkSideMinCut(result);
120 if (underlying_max_flow_ ==
nullptr)
return FlowModel();
121 return underlying_max_flow_->CreateFlowModel();
124 template <
typename Graph>
130 residual_arc_capacity_(),
131 first_admissible_arc_(),
135 use_global_update_(true),
136 use_two_phase_algorithm_(true),
137 process_node_by_height_(true),
145 if (max_num_nodes > 0) {
156 if (max_num_arcs > 0) {
162 template <
typename Graph>
166 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
167 if (residual_arc_capacity_[arc] < 0) {
174 template <
typename Graph>
180 const FlowQuantity free_capacity = residual_arc_capacity_[arc];
181 const FlowQuantity capacity_delta = new_capacity - Capacity(arc);
182 if (capacity_delta == 0) {
185 status_ = NOT_SOLVED;
186 if (free_capacity + capacity_delta >= 0) {
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]);
204 SetCapacityAndClearFlow(arc, new_capacity);
208 template <
typename Graph>
219 residual_arc_capacity_.Set(Opposite(arc), -new_flow);
220 residual_arc_capacity_.Set(arc,
capacity - new_flow);
221 status_ = NOT_SOLVED;
224 template <
typename Graph>
226 std::vector<NodeIndex>* result) {
227 ComputeReachableNodes<false>(source_, result);
230 template <
typename Graph>
232 ComputeReachableNodes<true>(sink_, result);
235 template <
typename Graph>
239 if (node_excess_[source_] != -node_excess_[sink_]) {
240 LOG(DFATAL) <<
"-node_excess_[source_] = " << -node_excess_[source_]
241 <<
" != node_excess_[sink_] = " << node_excess_[sink_];
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]
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";
262 if (opposite_capacity < 0) {
263 LOG(DFATAL) <<
"residual_arc_capacity_[" << opposite
264 <<
"] = " << opposite_capacity <<
" < 0";
268 if (direct_capacity + opposite_capacity < 0) {
269 LOG(DFATAL) <<
"initial capacity [" << arc
270 <<
"] = " << direct_capacity + opposite_capacity <<
" < 0";
277 template <
typename Graph>
282 const NodeIndex num_nodes = graph_->num_nodes();
283 std::vector<bool> is_reached(num_nodes,
false);
284 std::vector<NodeIndex> to_process;
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();
294 if (residual_arc_capacity_[arc] > 0) {
296 if (!is_reached[
head]) {
297 is_reached[
head] =
true;
298 to_process.push_back(
head);
303 return is_reached[sink_];
306 template <
typename Graph>
312 DCHECK(!IsAdmissible(arc)) << DebugString(
"CheckRelabelPrecondition:", arc);
317 template <
typename Graph>
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",
329 Flow(arc), node_potential_[
tail], node_potential_[
head],
330 node_excess_[
tail], node_excess_[
head]);
333 template <
typename Graph>
335 status_ = NOT_SOLVED;
336 if (check_input_ && !CheckInputConsistency()) {
345 const NodeIndex num_nodes = graph_->num_nodes();
346 if (sink_ >= num_nodes || source_ >= num_nodes) {
352 if (use_global_update_) {
353 RefineWithGlobalUpdate();
358 if (!CheckResult()) {
359 status_ = BAD_RESULT;
362 if (GetOptimalFlow() < kMaxFlowQuantity && AugmentingPathExists()) {
363 LOG(
ERROR) <<
"The algorithm terminated, but the flow is not maximal!";
364 status_ = BAD_RESULT;
368 DCHECK_EQ(node_excess_[sink_], -node_excess_[source_]);
370 if (GetOptimalFlow() == kMaxFlowQuantity && AugmentingPathExists()) {
372 status_ = INT_OVERFLOW;
378 template <
typename Graph>
385 node_excess_.SetAll(0);
387 for (
ArcIndex arc = 0; arc < num_arcs; ++arc) {
388 SetCapacityAndClearFlow(arc, Capacity(arc));
393 node_potential_.SetAll(0);
394 node_potential_.Set(source_, graph_->num_nodes());
400 for (
NodeIndex node = 0; node < num_nodes; ++node) {
401 first_admissible_arc_[node] = Graph::kNilArc;
409 template <
typename Graph>
412 const NodeIndex num_nodes = graph_->num_nodes();
421 std::vector<bool> stored(num_nodes,
false);
422 stored[sink_] =
true;
426 std::vector<bool> visited(num_nodes,
false);
427 visited[sink_] =
true;
431 std::vector<ArcIndex> arc_stack;
435 std::vector<int> index_branch;
438 std::vector<NodeIndex> reverse_topological_order;
447 arc_stack.push_back(arc);
450 visited[source_] =
true;
453 while (!arc_stack.empty()) {
462 reverse_topological_order.push_back(node);
463 DCHECK(!index_branch.empty());
464 index_branch.pop_back();
466 arc_stack.pop_back();
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);
478 for (OutgoingArcIterator it(*graph_, node); it.Ok(); it.Next()) {
482 if (flow > 0 && !stored[
head]) {
483 if (!visited[
head]) {
484 arc_stack.push_back(arc);
490 int cycle_begin = index_branch.size();
491 while (cycle_begin > 0 &&
492 Head(arc_stack[index_branch[cycle_begin - 1]]) !=
head) {
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;
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;
530 if (first_saturated_index < index_branch.size()) {
531 arc_stack.resize(index_branch[first_saturated_index]);
532 index_branch.resize(first_saturated_index);
542 DCHECK(arc_stack.empty());
543 DCHECK(index_branch.empty());
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) {
554 std::min(node_excess_[node], residual_arc_capacity_[opposite_arc]);
555 PushFlow(flow, opposite_arc);
556 if (node_excess_[node] == 0)
break;
561 DCHECK_EQ(-node_excess_[source_], node_excess_[sink_]);
564 template <
typename Graph>
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;
583 const int num_passes = use_two_phase_algorithm_ ? 1 : 2;
584 for (
int pass = 0; pass < num_passes; ++pass) {
586 bfs_queue_.push_back(sink_);
588 bfs_queue_.push_back(source_);
591 while (queue_index != bfs_queue_.size()) {
592 const NodeIndex node = bfs_queue_[queue_index];
594 const NodeIndex candidate_distance = node_potential_[node] + 1;
602 if (node_in_bfs_queue_[
head])
continue;
610 const ArcIndex opposite_arc = Opposite(arc);
611 if (residual_arc_capacity_[opposite_arc] > 0) {
625 if (node_excess_[
head] > 0) {
627 node_excess_[
head], residual_arc_capacity_[opposite_arc]);
628 PushFlow(flow, opposite_arc);
632 if (residual_arc_capacity_[opposite_arc] == 0)
continue;
637 node_potential_[
head] = candidate_distance;
638 node_in_bfs_queue_[
head] =
true;
639 bfs_queue_.push_back(
head);
656 for (
NodeIndex node = 0; node < num_nodes; ++node) {
657 if (!node_in_bfs_queue_[node]) {
658 node_potential_[node] = 2 * num_nodes - 1;
664 DCHECK(IsEmptyActiveNodeContainer());
665 for (
int i = 1; i < bfs_queue_.size(); ++i) {
667 if (node_excess_[node] > 0) {
669 PushActiveNode(node);
674 template <
typename Graph>
677 const NodeIndex num_nodes = graph_->num_nodes();
681 if (node_excess_[sink_] == kMaxFlowQuantity)
return false;
682 if (node_excess_[source_] == -kMaxFlowQuantity)
return false;
684 bool flow_pushed =
false;
690 if (flow == 0 || node_potential_[Head(arc)] >= num_nodes)
continue;
694 const FlowQuantity current_flow_out_of_source = -node_excess_[source_];
696 DCHECK_GE(current_flow_out_of_source, 0) << current_flow_out_of_source;
698 kMaxFlowQuantity - current_flow_out_of_source;
699 if (capped_flow < flow) {
706 if (capped_flow == 0)
return true;
707 PushFlow(capped_flow, arc);
717 template <
typename Graph>
721 DCHECK_GE(residual_arc_capacity_[Opposite(arc)] + flow, 0);
722 DCHECK_GE(residual_arc_capacity_[arc] - flow, 0);
731 residual_arc_capacity_[arc] -= flow;
732 residual_arc_capacity_[Opposite(arc)] += flow;
735 node_excess_[Tail(arc)] -= flow;
736 node_excess_[Head(arc)] += flow;
739 template <
typename Graph>
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) {
749 PushActiveNode(node);
754 template <
typename Graph>
779 while (SaturateOutgoingArcsFromSource()) {
780 DCHECK(IsEmptyActiveNodeContainer());
781 InitializeActiveNodeContainer();
782 while (!IsEmptyActiveNodeContainer()) {
783 const NodeIndex node = GetAndRemoveFirstActiveNode();
784 if (node == source_ || node == sink_)
continue;
787 if (use_two_phase_algorithm_) {
788 PushFlowExcessBackToSource();
793 template <
typename Graph>
800 std::vector<int> skip_active_node;
802 while (SaturateOutgoingArcsFromSource()) {
806 skip_active_node.assign(num_nodes, 0);
807 skip_active_node[sink_] = 2;
808 skip_active_node[source_] = 2;
810 while (!IsEmptyActiveNodeContainer()) {
811 const NodeIndex node = GetAndRemoveFirstActiveNode();
812 if (skip_active_node[node] > 1) {
813 if (node != sink_ && node != source_) ++num_skipped;
816 const NodeIndex old_height = node_potential_[node];
835 if (node_potential_[node] > old_height + 1) {
836 ++skip_active_node[node];
839 }
while (num_skipped > 0);
840 if (use_two_phase_algorithm_) {
841 PushFlowExcessBackToSource();
846 template <
typename Graph>
849 const NodeIndex num_nodes = graph_->num_nodes();
853 first_admissible_arc_[node]);
854 it.Ok(); it.Next()) {
856 if (IsAdmissible(arc)) {
859 if (node_excess_[
head] == 0) {
862 PushActiveNode(
head);
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;
874 if (use_two_phase_algorithm_ && node_potential_[node] >= num_nodes)
break;
878 template <
typename Graph>
885 ArcIndex first_admissible_arc = Graph::kNilArc;
889 if (residual_arc_capacity_[arc] > 0) {
891 NodeHeight head_height = node_potential_[Head(arc)];
892 if (head_height < min_height) {
893 min_height = head_height;
894 first_admissible_arc = arc;
898 if (min_height + 1 == node_potential_[node])
break;
902 DCHECK_NE(first_admissible_arc, Graph::kNilArc);
903 node_potential_[node] = min_height + 1;
908 first_admissible_arc_[node] = first_admissible_arc;
911 template <
typename Graph>
916 template <
typename Graph>
918 return IsArcValid(arc) && arc >= 0;
921 template <
typename Graph>
926 template <
typename Graph>
930 template <
typename Graph>
931 template <
bool reverse>
933 NodeIndex start, std::vector<NodeIndex>* result) {
937 const NodeIndex num_nodes = graph_->num_nodes();
938 if (start >= num_nodes) {
940 result->push_back(start);
944 node_in_bfs_queue_.assign(num_nodes,
false);
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];
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);
962 *result = bfs_queue_;
965 template <
typename Graph>
968 model.set_problem_type(FlowModel::MAX_FLOW);
969 for (
int n = 0; n < graph_->num_nodes(); ++n) {
970 Node* node =
model.add_node();
972 if (n == source_) node->set_supply(1);
973 if (n == sink_) node->set_supply(-1);
975 for (
int a = 0;
a < graph_->num_arcs(); ++
a) {
977 arc->set_tail_node_id(graph_->Tail(
a));
978 arc->set_head_node_id(graph_->Head(
a));
979 arc->set_capacity(Capacity(
a));