123#ifndef OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_
124#define OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_
134#include "absl/strings/string_view.h"
136#include "ortools/graph/flow_problem.pb.h"
156template <
typename Element,
typename IntegerPriority>
178 void Push(Element element, IntegerPriority priority);
186 Element PopBack(std::vector<std::pair<Element, IntegerPriority> >* queue);
191 std::vector<std::pair<Element, IntegerPriority> > even_queue_;
192 std::vector<std::pair<Element, IntegerPriority> > odd_queue_;
247template <
typename Graph,
typename ArcFlowT = int64_t,
248 typename FlowSumT = int64_t>
256 static_assert(std::is_signed_v<FlowSumType>);
361 DCHECK_EQ(tail,
Tail(arc));
363 potentials[tail] == potentials[
Head(arc)] + 1;
460 template <
bool reverse>
465 std::numeric_limits<FlowSumType>::max();
539template <
typename Element,
typename IntegerPriority>
542 return even_queue_.empty() && odd_queue_.empty();
545template <
typename Element,
typename IntegerPriority>
551template <
typename Element,
typename IntegerPriority>
553 Element element, IntegerPriority priority) {
555 DCHECK(even_queue_.empty() || priority >= even_queue_.back().second - 1);
556 DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second - 1);
562 DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second);
563 odd_queue_.push_back(std::make_pair(element, priority));
565 DCHECK(even_queue_.empty() || priority >= even_queue_.back().second);
566 even_queue_.push_back(std::make_pair(element, priority));
570template <
typename Element,
typename IntegerPriority>
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_);
578 return PopBack(&even_queue_);
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;
591template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
601 DCHECK(
graph->IsNodeValid(source));
602 DCHECK(
graph->IsNodeValid(sink));
604 if (max_num_nodes > 0) {
609 node_excess_ = std::make_unique<FlowSumType[]>(max_num_nodes);
615 if (max_num_arcs > 0) {
627template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
631 DCHECK_LE(0, new_capacity);
652template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
654 std::vector<NodeIndex>* result) {
658template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
660 std::vector<NodeIndex>* result) {
664template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
675 LOG(DFATAL) <<
"node_excess_[" << node <<
"] = " <<
node_excess_[node]
685 if (direct_capacity < 0) {
686 LOG(DFATAL) <<
"residual_arc_capacity_[" << arc
687 <<
"] = " << direct_capacity <<
" < 0";
690 if (opposite_capacity < 0) {
691 LOG(DFATAL) <<
"residual_arc_capacity_[" << opposite
692 <<
"] = " << opposite_capacity <<
" < 0";
696 if (direct_capacity + opposite_capacity < 0) {
697 LOG(DFATAL) <<
"initial capacity [" << arc
698 <<
"] = " << direct_capacity + opposite_capacity <<
" < 0";
704 LOG(DFATAL) <<
"The algorithm terminated, but the flow is not maximal!";
711template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
717 std::vector<bool> is_reached(num_nodes,
false);
718 std::vector<NodeIndex> to_process;
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)) {
728 if (!is_reached[head]) {
729 is_reached[head] =
true;
730 to_process.push_back(head);
735 return is_reached[
sink_];
738template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
742 for (
const ArcIndex arc :
graph_->OutgoingOrOppositeIncomingArcs(node)) {
749template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
751 absl::string_view context,
ArcIndex arc)
const {
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",
765template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
794template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
813 for (
ArcIndex arc = 0; arc < num_arcs; ++arc) {
819 for (
ArcIndex arc = 0; arc < num_arcs; ++arc) {
830 for (
NodeIndex node = 0; node < num_nodes; ++node) {
832 for (
const ArcIndex arc :
graph_->OutgoingOrOppositeIncomingArcs(node)) {
843template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
855 std::vector<bool> stored(num_nodes,
false);
856 stored[
sink_] =
true;
860 std::vector<bool> visited(num_nodes,
false);
861 visited[
sink_] =
true;
865 std::vector<ArcIndex> arc_stack;
869 std::vector<int> index_branch;
872 std::vector<NodeIndex> reverse_topological_order;
879 arc_stack.push_back(arc);
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());
897 reverse_topological_order.push_back(node);
898 DCHECK(!index_branch.empty());
899 index_branch.pop_back();
901 arc_stack.pop_back();
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);
915 if (flow > 0 && !stored[head]) {
916 if (!visited[head]) {
917 arc_stack.push_back(arc);
923 int cycle_begin = index_branch.size();
924 while (cycle_begin > 0 &&
925 Head(arc_stack[index_branch[cycle_begin - 1]]) != head) {
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]];
936 arc_flow <= flow_on_cycle) {
937 flow_on_cycle = arc_flow;
938 first_saturated_index = i;
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;
955 DCHECK_GT(
Flow(arc_on_cycle), 0);
964 if (first_saturated_index < index_branch.size()) {
965 arc_stack.resize(index_branch[first_saturated_index]);
966 index_branch.resize(first_saturated_index);
976 DCHECK(arc_stack.empty());
977 DCHECK(index_branch.empty());
981 for (
int i = 0; i < reverse_topological_order.size(); i++) {
982 const NodeIndex node = reverse_topological_order[i];
984 for (
const ArcIndex arc :
graph_->OutgoingOrOppositeIncomingArcs(node)) {
999template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1003 int queue_index = 0;
1023 for (
const ArcIndex arc :
graph_->OutgoingOrOppositeIncomingArcs(node)) {
1051 if (node_excess[head] > 0) {
1079 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1088 for (
int i = 1; i <
bfs_queue_.size(); ++i) {
1097template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1108 bool flow_pushed =
false;
1118 DCHECK_GE(flow, 0) << flow;
1119 DCHECK_GE(current_flow_out_of_source, 0) << current_flow_out_of_source;
1121 if (capped_flow <
static_cast<FlowSumType>(flow)) {
1128 if (capped_flow == 0)
return true;
1139template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1163template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1172 node_excess[tail]));
1180 node_excess[tail] -=
static_cast<FlowSumType>(flow);
1184template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1190 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1201template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1206 std::vector<int> skip_active_node;
1225 skip_active_node.assign(num_nodes, 0);
1226 skip_active_node[
sink_] = 2;
1227 skip_active_node[
source_] = 2;
1231 if (skip_active_node[node] > 1) {
1255 ++skip_active_node[node];
1258 }
while (num_skipped > 0);
1269template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1283 graph_->OutgoingOrOppositeIncomingArcsStartingFrom(
1288 if (node_excess[head] == 0) {
1294 if (node_excess[node] == 0) {
1295 first_admissible_arc[node] = arc;
1304 if (node_potentials[node] >= num_nodes)
break;
1308template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1314 NodeHeight min_height = std::numeric_limits<NodeHeight>::max();
1316 for (
const ArcIndex arc :
graph_->OutgoingOrOppositeIncomingArcs(node)) {
1320 if (head_height < min_height) {
1321 min_height = head_height;
1322 first_admissible_arc = arc;
1339template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1342 return graph_->OppositeArc(arc);
1345template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1351template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1353 return graph_->IsArcValid(arc);
1356template <
typename Graph,
typename ArcFlowT,
typename FlowSumT>
1357template <
bool reverse>
1359 NodeIndex start, std::vector<NodeIndex>* result) {
1364 if (start >= num_nodes) {
1366 result->push_back(start);
1372 int queue_index = 0;
1378 for (
const ArcIndex arc :
graph_->OutgoingOrOppositeIncomingArcs(node)) {
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();
1396 if (n ==
source_) node->set_supply(1);
1397 if (n ==
sink_) node->set_supply(-1);
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));
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_
void Discharge(NodeIndex node)
void RefineWithGlobalUpdate()
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
bool SaturateOutgoingArcsFromSource()
bool IsArcValid(ArcIndex arc) const
const util::ReverseArcStaticGraph< NodeIndex, ArcIndex > * graph() const
FlowSumType Flow(ArcIndex arc) const
bool AugmentingPathExists() 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_
bool IsArcDirect(ArcIndex arc) const
void PushActiveNode(const NodeIndex &node)
Push element to the active node container.
void PushFlowExcessBackToSource()
bool Solve()
Returns true if a maximum flow was solved.
void Relabel(NodeIndex node)
Graph::NodeIndex NodeIndex
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.
std::vector< bool > node_in_bfs_queue_
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.
bool IsEmpty() const
Is the queue empty?
void Push(Element element, IntegerPriority priority)
PriorityQueueWithRestrictedPush & operator=(const PriorityQueueWithRestrictedPush &)=delete
PriorityQueueWithRestrictedPush(const PriorityQueueWithRestrictedPush &)=delete
This type is neither copyable nor movable.
PriorityQueueWithRestrictedPush()
void Clear()
Clears the queue.
Base class to print a nice summary of a group of statistics.
static constexpr int32_t kNilArc
static constexpr bool kHasNegativeReverseArcs
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)
#define SCOPED_TIME_STAT(stats)