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();
124template <
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) {
162template <
typename Graph>
166 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
167 if (residual_arc_capacity_[arc] < 0) {
174template <
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);
208template <
typename Graph>
219 residual_arc_capacity_.Set(Opposite(arc), -new_flow);
220 residual_arc_capacity_.Set(arc,
capacity - new_flow);
221 status_ = NOT_SOLVED;
224template <
typename Graph>
226 std::vector<NodeIndex>* result) {
227 ComputeReachableNodes<false>(source_, result);
230template <
typename Graph>
232 ComputeReachableNodes<true>(sink_, result);
235template <
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";
277template <
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_];
306template <
typename Graph>
312 DCHECK(!IsAdmissible(arc)) << DebugString(
"CheckRelabelPrecondition:", arc);
317template <
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]);
333template <
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;
378template <
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;
409template <
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_]);
564template <
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);
674template <
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);
717template <
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;
739template <
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);
754template <
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();
793template <
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();
846template <
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;
878template <
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;
911template <
typename Graph>
916template <
typename Graph>
918 return IsArcValid(arc) && arc >= 0;
921template <
typename Graph>
926template <
typename Graph>
930template <
typename Graph>
931template <
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_;
965template <
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));
#define DCHECK_LE(val1, val2)
#define DCHECK_NE(val1, val2)
#define DCHECK_GE(val1, val2)
#define DCHECK_GT(val1, val2)
#define DCHECK(condition)
#define DCHECK_EQ(val1, val2)
#define VLOG(verboselevel)
Graph::OutgoingArcIterator OutgoingArcIterator
bool CheckInputConsistency() const
void Relabel(NodeIndex node)
std::vector< NodeIndex > active_nodes_
void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity)
std::string DebugString(const std::string &context, ArcIndex arc) const
bool SaturateOutgoingArcsFromSource()
Graph::OutgoingOrOppositeIncomingArcIterator OutgoingOrOppositeIncomingArcIterator
bool CheckRelabelPrecondition(NodeIndex node) const
std::vector< NodeIndex > bfs_queue_
ArcIndexArray first_admissible_arc_
void SetArcFlow(ArcIndex arc, FlowQuantity new_flow)
void PushFlowExcessBackToSource()
Graph::NodeIndex NodeIndex
void GetSourceSideMinCut(std::vector< NodeIndex > *result)
NodeHeightArray node_potential_
void InitializeActiveNodeContainer()
const Graph * graph() const
QuantityArray residual_arc_capacity_
void PushFlow(FlowQuantity flow, ArcIndex arc)
void ComputeReachableNodes(NodeIndex start, std::vector< NodeIndex > *result)
bool IsArcValid(ArcIndex arc) const
GenericMaxFlow(const Graph *graph, NodeIndex source, NodeIndex sink)
void RefineWithGlobalUpdate()
bool AugmentingPathExists() const
void GetSinkSideMinCut(std::vector< NodeIndex > *result)
ArcIndex Opposite(ArcIndex arc) const
bool IsArcDirect(ArcIndex arc) const
void Discharge(NodeIndex node)
QuantityArray node_excess_
FlowModel CreateFlowModel()
FlowQuantity Flow(ArcIndex arc) const
NodeIndex NumNodes() const
FlowModel CreateFlowModelOfLastSolve()
void GetSourceSideMinCut(std::vector< NodeIndex > *result)
Status Solve(NodeIndex source, NodeIndex sink)
FlowQuantity OptimalFlow() const
ArcIndex AddArcWithCapacity(NodeIndex tail, NodeIndex head, FlowQuantity capacity)
FlowQuantity Capacity(ArcIndex arc) const
NodeIndex Head(ArcIndex arc) const
NodeIndex Tail(ArcIndex arc) const
void SetArcCapacity(ArcIndex arc, FlowQuantity capacity)
void GetSinkSideMinCut(std::vector< NodeIndex > *result)
bool Reserve(int64 new_min_index, int64 new_max_index)
GurobiMPCallbackContext * context
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
std::pair< int64, int64 > Arc
#define IF_STATS_ENABLED(instructions)
#define SCOPED_TIME_STAT(stats)
static ArcIndex ArcReservation(const Graph &graph)
static NodeIndex NodeReservation(const Graph &graph)
static bool IsArcValid(const Graph &graph, ArcIndex arc)
static ArcIndex OppositeArc(const Graph &graph, ArcIndex arc)