diff --git a/ortools/graph/BUILD.bazel b/ortools/graph/BUILD.bazel index 6e0d2a41e6..79918ea430 100644 --- a/ortools/graph/BUILD.bazel +++ b/ortools/graph/BUILD.bazel @@ -96,6 +96,27 @@ cc_test( ], ) +cc_library( + name = "minimum_vertex_cover", + srcs = ["minimum_vertex_cover.cc"], + hdrs = ["minimum_vertex_cover.h"], + deps = [ + ":ebert_graph", + ":max_flow", + "@com_google_absl//absl/log:check", + ], +) + +cc_test( + name = "minimum_vertex_cover_test", + srcs = ["minimum_vertex_cover_test.cc"], + deps = [ + ":minimum_vertex_cover", + "//ortools/base:gmock_main", + "@com_google_absl//absl/algorithm:container", + ], +) + cc_library( name = "multi_dijkstra", hdrs = ["multi_dijkstra.h"], @@ -443,40 +464,57 @@ cc_library( srcs = ["max_flow.cc"], hdrs = ["max_flow.h"], deps = [ - ":ebert_graph", ":flow_problem_cc_proto", + ":generic_max_flow", ":graph", - ":graphs", - "//ortools/base", - "//ortools/util:stats", - "//ortools/util:zvector", - "@com_google_absl//absl/log:check", - "@com_google_absl//absl/memory", - "@com_google_absl//absl/strings", - "@com_google_absl//absl/strings:str_format", ], ) cc_test( name = "max_flow_test", - size = "medium", srcs = ["max_flow_test.cc"], data = ["//ortools/graph/testdata:max_flow_test1.pb.txt"], deps = [ - ":graph", - ":graphs", + ":ebert_graph", + ":flow_problem_cc_proto", ":max_flow", - "//ortools/base", "//ortools/base:gmock_main", "//ortools/base:path", - "//ortools/linear_solver", "//ortools/util:file_util", - "@com_google_absl//absl/algorithm:container", + "@com_google_protobuf//:protobuf", + ], +) + +cc_library( + name = "generic_max_flow", + hdrs = ["generic_max_flow.h"], + deps = [ + ":ebert_graph", + ":flow_problem_cc_proto", + ":graphs", + "//ortools/base", + "//ortools/util:stats", + "//ortools/util:zvector", + "@com_google_absl//absl/strings", + ], +) + +cc_test( + name = "generic_max_flow_test", + size = "medium", + srcs = ["generic_max_flow_test.cc"], + deps = [ + ":ebert_graph", + ":generic_max_flow", + ":graph", + ":graphs", + "//ortools/base", + "//ortools/base:gmock_main", + "//ortools/linear_solver", "@com_google_absl//absl/random", "@com_google_absl//absl/strings:str_format", "@com_google_absl//absl/types:span", "@com_google_benchmark//:benchmark", - "@com_google_protobuf//:protobuf", ], ) @@ -496,14 +534,13 @@ cc_library( ":graph", ":graphs", ":max_flow", - "//ortools/base", - "//ortools/base:dump_vars", "//ortools/base:mathutil", - "//ortools/base:types", "//ortools/util:saturated_arithmetic", "//ortools/util:stats", "//ortools/util:zvector", "@com_google_absl//absl/flags:flag", + "@com_google_absl//absl/log", + "@com_google_absl//absl/log:check", "@com_google_absl//absl/strings", "@com_google_absl//absl/strings:str_format", ], diff --git a/ortools/graph/generic_max_flow.h b/ortools/graph/generic_max_flow.h new file mode 100644 index 0000000000..8976696f13 --- /dev/null +++ b/ortools/graph/generic_max_flow.h @@ -0,0 +1,1364 @@ +// Copyright 2010-2024 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// An implementation of a push-relabel algorithm for the max flow problem. +// +// In the following, we consider a graph G = (V,E,s,t) where V denotes the set +// of nodes (vertices) in the graph, E denotes the set of arcs (edges). s and t +// denote distinguished nodes in G called source and target. n = |V| denotes the +// number of nodes in the graph, and m = |E| denotes the number of arcs in the +// graph. +// +// Each arc (v,w) is associated a capacity c(v,w). +// +// A flow is a function from E to R such that: +// +// a) f(v,w) <= c(v,w) for all (v,w) in E (capacity constraint.) +// +// b) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint.) +// +// c) sum on v f(v,w) = 0 (flow conservation.) +// +// The goal of this algorithm is to find the maximum flow from s to t, i.e. +// for example to maximize sum v f(s,v). +// +// The starting reference for this class of algorithms is: +// A.V. Goldberg and R.E. Tarjan. A new approach to the maximum flow problem. +// ACM Symposium on Theory of Computing, pp. 136-146. +// http://portal.acm.org/citation.cfm?id=12144. +// +// The basic idea of the algorithm is to handle preflows instead of flows, +// and to refine preflows until a maximum flow is obtained. +// A preflow is like a flow, except that the inflow can be larger than the +// outflow. If it is the case at a given node v, it is said that there is an +// excess at node v, and inflow = outflow + excess. +// +// More formally, a preflow is a function f such that: +// +// 1) f(v,w) <= c(v,w) for all (v,w) in E (capacity constraint). c(v,w) is a +// value representing the maximum capacity for arc (v,w). +// +// 2) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint) +// +// 3) excess(v) = sum on u f(u,v) >= 0 is the excess at node v, the +// algebraic sum of all the incoming preflows at this node. +// +// Each node has an associated "height", in addition to its excess. The +// height of the source is defined to be equal to n, and cannot change. The +// height of the target is defined to be zero, and cannot change either. The +// height of all the other nodes is initialized at zero and is updated during +// the algorithm (see below). For those who want to know the details, the height +// of a node, corresponds to a reduced cost, and this enables one to prove that +// the algorithm actually computes the max flow. Note that the height of a node +// can be initialized to the distance to the target node in terms of number of +// nodes. This has not been tried in this implementation. +// +// A node v is said to be *active* if excess(v) > 0. +// +// In this case the following operations can be applied to it: +// +// - if there are *admissible* incident arcs, i.e. arcs which are not saturated, +// and whose head's height is lower than the height of the active node +// considered, a PushFlow operation can be applied. It consists in sending as +// much flow as both the excess at the node and the capacity of the arc +// permit. +// - if there are no admissible arcs, the active node considered is relabeled, +// i.e. its height is increased to 1 + the minimum height of its neighboring +// nodes on admissible arcs. +// This is implemented in Discharge, which itself calls PushFlow and Relabel. +// +// Before running Discharge, it is necessary to initialize the algorithm with a +// preflow. This is done in InitializePreflow, which saturates all the arcs +// leaving the source node, and sets the excess at the heads of those arcs +// accordingly. +// +// The algorithm terminates when there are no remaining active nodes, i.e. all +// the excesses at all nodes are equal to zero. In this case, a maximum flow is +// obtained. +// +// The complexity of this algorithm depends amongst other things on the choice +// of the next active node. It has been shown, for example in: +// L. Tuncel, "On the Complexity of Preflow-Push Algorithms for Maximum-Flow +// Problems", Algorithmica 11(4): 353-359 (1994). +// and +// J. Cheriyan and K. Mehlhorn, "An analysis of the highest-level selection rule +// in the preflow-push max-flow algorithm", Information processing letters, +// 69(5):239-242 (1999). +// http://www.math.uwaterloo.ca/~jcheriya/PS_files/me3.0.ps +// +// ...that choosing the active node with the highest level yields a +// complexity of O(n^2 * sqrt(m)). +// +// TODO(user): implement the above active node choice rule. +// +// This has been validated experimentally in: +// R.K. Ahuja, M. Kodialam, A.K. Mishra, and J.B. Orlin, "Computational +// Investigations of Maximum Flow Algorithms", EJOR 97:509-542(1997). +// http://jorlin.scripts.mit.edu/docs/publications/58-comput%20investigations%20of.pdf. +// +// +// TODO(user): an alternative would be to evaluate: +// A.V. Goldberg, "The Partial Augment-Relabel Algorithm for the Maximum Flow +// Problem.” In Proceedings of Algorithms ESA, LNCS 5193:466-477, Springer 2008. +// http://www.springerlink.com/index/5535k2j1mt646338.pdf +// +// An interesting general reference on network flows is: +// R. K. Ahuja, T. L. Magnanti, J. B. Orlin, "Network Flows: Theory, Algorithms, +// and Applications," Prentice Hall, 1993, ISBN: 978-0136175490, +// http://www.amazon.com/dp/013617549X +// +// Keywords: Push-relabel, max-flow, network, graph, Goldberg, Tarjan, Dinic, +// Dinitz. + +#ifndef OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_ +#define OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_ + +#include +#include +#include +#include + +#include "absl/strings/string_view.h" +#include "ortools/base/logging.h" +#include "ortools/graph/ebert_graph.h" +#include "ortools/graph/flow_problem.pb.h" +#include "ortools/graph/graphs.h" +#include "ortools/util/stats.h" +#include "ortools/util/zvector.h" + +namespace operations_research { + +// Specific but efficient priority queue implementation. The priority type must +// be an integer. The queue allows to retrieve the element with highest priority +// but only allows pushes with a priority greater or equal to the highest +// priority in the queue minus one. All operations are in O(1) and the memory is +// in O(num elements in the queue). Elements with the same priority are +// retrieved with LIFO order. +// +// Note(user): As far as I know, this is an original idea and is the only code +// that use this in the Maximum Flow context. Papers usually refer to an +// height-indexed array of simple linked lists of active node with the same +// height. Even worse, sometimes they use double-linked list to allow arbitrary +// height update in order to detect missing height (used for the Gap heuristic). +// But this can actually be implemented a lot more efficiently by just +// maintaining the height distribution of all the node in the graph. +template +class PriorityQueueWithRestrictedPush { + public: + PriorityQueueWithRestrictedPush() : even_queue_(), odd_queue_() {} + +#ifndef SWIG + // This type is neither copyable nor movable. + PriorityQueueWithRestrictedPush(const PriorityQueueWithRestrictedPush&) = + delete; + PriorityQueueWithRestrictedPush& operator=( + const PriorityQueueWithRestrictedPush&) = delete; +#endif + + // Is the queue empty? + bool IsEmpty() const; + + // Clears the queue. + void Clear(); + + // Push a new element in the queue. Its priority must be greater or equal to + // the highest priority present in the queue, minus one. This condition is + // DCHECKed, and violating it yields erroneous queue behavior in NDEBUG mode. + void Push(Element element, IntegerPriority priority); + + // Returns the element with highest priority and remove it from the queue. + // IsEmpty() must be false, this condition is DCHECKed. + Element Pop(); + + private: + // Helper function to get the last element of a vector and pop it. + Element PopBack(std::vector >* queue); + + // This is the heart of the algorithm. basically we split the elements by + // parity of their priority and the precondition on the Push() ensures that + // both vectors are always sorted by increasing priority. + std::vector > even_queue_; + std::vector > odd_queue_; +}; + +// We want an enum for the Status of a max flow run, and we want this +// enum to be scoped under GenericMaxFlow<>. Unfortunately, swig +// doesn't handle templated enums very well, so we need a base, +// untemplated class to hold it. +class MaxFlowStatusClass { + public: + enum Status { + NOT_SOLVED, // The problem was not solved, or its data were edited. + OPTIMAL, // Solve() was called and found an optimal solution. + INT_OVERFLOW, // There is a feasible flow > max possible flow. + + // TODO(user): These are no longer used. Remove. + BAD_INPUT, + BAD_RESULT + }; +}; + +// Generic MaxFlow (there is a default MaxFlow specialization defined below) +// that works with StarGraph and all the reverse arc graphs from graph.h, see +// the end of max_flow.cc for the exact types this class is compiled for. +template +class GenericMaxFlow : public MaxFlowStatusClass { + public: + typedef typename Graph::NodeIndex NodeIndex; + typedef typename Graph::ArcIndex ArcIndex; + typedef typename Graph::OutgoingArcIterator OutgoingArcIterator; + typedef typename Graph::OutgoingOrOppositeIncomingArcIterator + OutgoingOrOppositeIncomingArcIterator; + typedef typename Graph::IncomingArcIterator IncomingArcIterator; + + // The height of a node never excess 2 times the number of node, so we + // use the same type as a Node index. + typedef NodeIndex NodeHeight; + + // Initialize a MaxFlow instance on the given graph. The graph does not need + // to be fully built yet, but its capacity reservation are used to initialize + // the memory of this class. source and sink must also be valid node of + // graph. + GenericMaxFlow(const Graph* graph, NodeIndex source, NodeIndex sink); + +#ifndef SWIG + // This type is neither copyable nor movable. + GenericMaxFlow(const GenericMaxFlow&) = delete; + GenericMaxFlow& operator=(const GenericMaxFlow&) = delete; +#endif + + virtual ~GenericMaxFlow() {} + + // Returns the graph associated to the current object. + const Graph* graph() const { return graph_; } + + // Returns the status of last call to Solve(). NOT_SOLVED is returned if + // Solve() has never been called or if the problem has been modified in such a + // way that the previous solution becomes invalid. + Status status() const { return status_; } + + // Returns the index of the node corresponding to the source of the network. + NodeIndex GetSourceNodeIndex() const { return source_; } + + // Returns the index of the node corresponding to the sink of the network. + NodeIndex GetSinkNodeIndex() const { return sink_; } + + // Sets the capacity for arc to new_capacity. + void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity); + + // Returns true if a maximum flow was solved. + bool Solve(); + + // Returns the total flow found by the algorithm. + FlowQuantity GetOptimalFlow() const { return node_excess_[sink_]; } + + // Returns the flow on arc using the equations given in the comment on + // residual_arc_capacity_. + FlowQuantity Flow(ArcIndex arc) const { + if (IsArcDirect(arc)) { + return residual_arc_capacity_[Opposite(arc)]; + } else { + return -residual_arc_capacity_[arc]; + } + } + + // Returns the capacity of arc using the equations given in the comment on + // residual_arc_capacity_. + FlowQuantity Capacity(ArcIndex arc) const { + if (IsArcDirect(arc)) { + return residual_arc_capacity_[arc] + + residual_arc_capacity_[Opposite(arc)]; + } else { + return 0; + } + } + + // Returns the nodes reachable from the source in the residual graph, the + // outgoing arcs of this set form a minimum cut. + void GetSourceSideMinCut(std::vector* result); + + // Returns the nodes that can reach the sink in the residual graph, the + // outgoing arcs of this set form a minimum cut. Note that if this is the + // complement of GetNodeReachableFromSource(), then the min-cut is unique. + // + // TODO(user): In the two-phases algorithm, we can get this minimum cut + // without doing the second phase. Add an option for this if there is a need + // to, note that the second phase is pretty fast so the gain will be small. + void GetSinkSideMinCut(std::vector* result); + + // Returns true if there exists a path from the source to the sink with + // remaining capacity. This allows us to easily check at the end that the flow + // we computed is indeed optimal (provided that all the conditions tested by + // CheckResult() also hold). + bool AugmentingPathExists() const; + + // Returns the protocol buffer representation of the current problem. + FlowModelProto CreateFlowModel(); + + protected: + // Checks whether the result is valid, i.e. that node excesses are all equal + // to zero (we have a flow) and that residual capacities are all non-negative + // or zero. + bool CheckResult() const; + + // Returns true if arc is admissible. + bool IsAdmissible(NodeIndex tail, ArcIndex arc, + const NodeHeight* potentials) const { + DCHECK_EQ(tail, Tail(arc)); + return residual_arc_capacity_[arc] > 0 && + potentials[tail] == potentials[Head(arc)] + 1; + } + + // Returns true if node is active, i.e. if its excess is positive and it + // is neither the source or the sink of the graph. + bool IsActive(NodeIndex node) const { + return (node != source_) && (node != sink_) && (node_excess_[node] > 0); + } + + // Sets the capacity of arc to 'capacity' and clears the flow on arc. + void SetCapacityAndClearFlow(ArcIndex arc, FlowQuantity capacity) { + residual_arc_capacity_.Set(arc, capacity); + residual_arc_capacity_.Set(Opposite(arc), 0); + } + + // Returns true if a precondition for Relabel is met, i.e. the outgoing arcs + // of node are all either saturated or the heights of their heads are greater + // or equal to the height of node. + bool CheckRelabelPrecondition(NodeIndex node) const; + + // Returns context concatenated with information about arc + // in a human-friendly way. + std::string DebugString(absl::string_view context, ArcIndex arc) const; + + // Initializes the container active_nodes_. + void InitializeActiveNodeContainer(); + + // Get the first element from the active node container. + NodeIndex GetAndRemoveFirstActiveNode() { + return active_node_by_height_.Pop(); + } + + // Push element to the active node container. + void PushActiveNode(const NodeIndex& node) { + active_node_by_height_.Push(node, node_potential_[node]); + } + + // Check the emptiness of the container. + bool IsEmptyActiveNodeContainer() { return active_node_by_height_.IsEmpty(); } + + // Performs optimization step. + void Refine(); + void RefineWithGlobalUpdate(); + + // Discharges an active node node by saturating its admissible adjacent arcs, + // if any, and by relabelling it when it becomes inactive. + void Discharge(NodeIndex node); + + // Initializes the preflow to a state that enables to run Refine. + void InitializePreflow(); + + // Clears the flow excess at each node by pushing the flow back to the source: + // - Do a depth-first search from the source in the direct graph to cancel + // flow cycles. + // - Then, return flow excess along the depth-first search tree (by pushing + // the flow in the reverse dfs topological order). + // The theoretical complexity is O(mn), but it is a lot faster in practice. + void PushFlowExcessBackToSource(); + + // Computes the best possible node potential given the current flow using a + // reverse breadth-first search from the sink in the reverse residual graph. + // This is an implementation of the global update heuristic mentioned in many + // max-flow papers. See for instance: B.V. Cherkassky, A.V. Goldberg, "On + // implementing push-relabel methods for the maximum flow problem", + // Algorithmica, 19:390-410, 1997. + // ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/94/1523/CS-TR-94-1523.pdf + void GlobalUpdate(); + + // Tries to saturate all the outgoing arcs from the source that can reach the + // sink. Most of the time, we can do that in one go, except when more flow + // than kMaxFlowQuantity can be pushed out of the source in which case we + // have to be careful. Returns true if some flow was pushed. + bool SaturateOutgoingArcsFromSource(); + + // Pushes flow on arc, i.e. consumes flow on residual_arc_capacity_[arc], + // and consumes -flow on residual_arc_capacity_[Opposite(arc)]. Updates + // node_excess_ at the tail and head of arc accordingly. + void PushFlow(FlowQuantity flow, NodeIndex tail, ArcIndex arc); + + // Same as PushFlow() but with a cached node_excess_.data(), this is faster + // in hot loop as we remove bound checking and the pointer is constant. + void FastPushFlow(FlowQuantity flow, NodeIndex tail, ArcIndex arc, + FlowQuantity* node_excess); + + // Relabels a node, i.e. increases its height by the minimum necessary amount. + // This version of Relabel is relaxed in a way such that if an admissible arc + // exists at the current node height, then the node is not relabeled. This + // enables us to deal with wrong values of first_admissible_arc_[node] when + // updating it is too costly. + void Relabel(NodeIndex node); + + // Handy member functions to make the code more compact. + NodeIndex Head(ArcIndex arc) const { return graph_->Head(arc); } + NodeIndex Tail(ArcIndex arc) const { return graph_->Tail(arc); } + ArcIndex Opposite(ArcIndex arc) const; + bool IsArcDirect(ArcIndex arc) const; + bool IsArcValid(ArcIndex arc) const; + + // Returns the set of nodes reachable from start in the residual graph or in + // the reverse residual graph (if reverse is true). + template + void ComputeReachableNodes(NodeIndex start, std::vector* result); + + // Maximum manageable flow. + static constexpr FlowQuantity kMaxFlowQuantity = + std::numeric_limits::max(); + + // A pointer to the graph passed as argument. + const Graph* graph_; + + // An array representing the excess for each node in graph_. + std::vector node_excess_; + + // An array representing the height function for each node in graph_. For a + // given node, this is a lower bound on the shortest path length from this + // node to the sink in the residual network. The height of a node always goes + // up during the course of a Solve(). + // + // Since initially we saturate all the outgoing arcs of the source, we can + // never reach the sink from the source in the residual graph. Initially we + // set the height of the source to n (the number of node of the graph) and it + // never changes. If a node as an height >= n, then this node can't reach the + // sink and its height minus n is a lower bound on the shortest path length + // from this node to the source in the residual graph. + std::vector node_potential_; + + // An array representing the residual_capacity for each arc in graph_. + // Residual capacities enable one to represent the capacity and flow for all + // arcs in the graph in the following manner. + // For all arc, residual_arc_capacity_[arc] = capacity[arc] - flow[arc] + // Moreover, for reverse arcs, capacity[arc] = 0 by definition, + // Also flow[Opposite(arc)] = -flow[arc] by definition. + // Therefore: + // - for a direct arc: + // flow[arc] = 0 - flow[Opposite(arc)] + // = capacity[Opposite(arc)] - flow[Opposite(arc)] + // = residual_arc_capacity_[Opposite(arc)] + // - for a reverse arc: + // flow[arc] = -residual_arc_capacity_[arc] + // Using these facts enables one to only maintain residual_arc_capacity_, + // instead of both capacity and flow, for each direct and indirect arc. This + // reduces the amount of memory for this information by a factor 2. + ZVector residual_arc_capacity_; + + // An array representing the first admissible arc for each node in graph_. + std::vector first_admissible_arc_; + + // A priority queue used for managing active nodes in the algorithm. It allows + // to select the active node with highest height before each Discharge(). + // Moreover, since all pushes from this node will be to nodes with height + // greater or equal to the initial discharged node height minus one, the + // PriorityQueueWithRestrictedPush is a perfect fit. + PriorityQueueWithRestrictedPush active_node_by_height_; + + // The index of the source node in graph_. + NodeIndex source_; + + // The index of the sink node in graph_. + NodeIndex sink_; + + // The status of the problem. + Status status_; + + // BFS queue used by the GlobalUpdate() function. We do not use a C++ queue + // because we need access to the vector for different optimizations. + std::vector node_in_bfs_queue_; + std::vector bfs_queue_; + + // Statistics about this class. + mutable StatsGroup stats_; +}; + +#if !SWIG + +// Default instance MaxFlow that uses StarGraph. Note that we cannot just use a +// typedef because of dependent code expecting MaxFlow to be a real class. +// TODO(user): Modify this code and remove it. +class MaxFlow : public GenericMaxFlow { + public: + MaxFlow(const StarGraph* graph, NodeIndex source, NodeIndex target) + : GenericMaxFlow(graph, source, target) {} +}; + +#endif // SWIG + +template +bool PriorityQueueWithRestrictedPush::IsEmpty() + const { + return even_queue_.empty() && odd_queue_.empty(); +} + +template +void PriorityQueueWithRestrictedPush::Clear() { + even_queue_.clear(); + odd_queue_.clear(); +} + +template +void PriorityQueueWithRestrictedPush::Push( + Element element, IntegerPriority priority) { + // Since users may rely on it, we DCHECK the exact condition. + DCHECK(even_queue_.empty() || priority >= even_queue_.back().second - 1); + DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second - 1); + + // Note that the DCHECK() below are less restrictive than the ones above but + // check a necessary and sufficient condition for the priority queue to behave + // as expected. + if (priority & 1) { + DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second); + odd_queue_.push_back(std::make_pair(element, priority)); + } else { + DCHECK(even_queue_.empty() || priority >= even_queue_.back().second); + even_queue_.push_back(std::make_pair(element, priority)); + } +} + +template +Element PriorityQueueWithRestrictedPush::Pop() { + DCHECK(!IsEmpty()); + if (even_queue_.empty()) return PopBack(&odd_queue_); + if (odd_queue_.empty()) return PopBack(&even_queue_); + if (odd_queue_.back().second > even_queue_.back().second) { + return PopBack(&odd_queue_); + } else { + return PopBack(&even_queue_); + } +} + +template +Element PriorityQueueWithRestrictedPush::PopBack( + std::vector >* queue) { + DCHECK(!queue->empty()); + Element element = queue->back().first; + queue->pop_back(); + return element; +} + +template +GenericMaxFlow::GenericMaxFlow(const Graph* graph, NodeIndex source, + NodeIndex sink) + : graph_(graph), + residual_arc_capacity_(), + source_(source), + sink_(sink), + stats_("MaxFlow") { + SCOPED_TIME_STAT(&stats_); + DCHECK(graph->IsNodeValid(source)); + DCHECK(graph->IsNodeValid(sink)); + const NodeIndex max_num_nodes = Graphs::NodeReservation(*graph_); + if (max_num_nodes > 0) { + node_excess_.assign(max_num_nodes, 0); + node_potential_.assign(max_num_nodes, 0); + first_admissible_arc_.assign(max_num_nodes, Graph::kNilArc); + bfs_queue_.reserve(max_num_nodes); + } + const ArcIndex max_num_arcs = Graphs::ArcReservation(*graph_); + if (max_num_arcs > 0) { + residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1); + residual_arc_capacity_.SetAll(0); + } +} + +template +void GenericMaxFlow::SetArcCapacity(ArcIndex arc, + FlowQuantity new_capacity) { + SCOPED_TIME_STAT(&stats_); + DCHECK_LE(0, new_capacity); + DCHECK(IsArcDirect(arc)); + const FlowQuantity free_capacity = residual_arc_capacity_[arc]; + const FlowQuantity capacity_delta = new_capacity - Capacity(arc); + if (capacity_delta == 0) { + return; // Nothing to do. + } + status_ = NOT_SOLVED; + if (free_capacity + capacity_delta >= 0) { + // The above condition is true if one of the two conditions is true: + // 1/ (capacity_delta > 0), meaning we are increasing the capacity + // 2/ (capacity_delta < 0 && free_capacity + capacity_delta >= 0) + // meaning we are reducing the capacity, but that the capacity + // reduction is not larger than the free capacity. + DCHECK((capacity_delta > 0) || + (capacity_delta < 0 && free_capacity + capacity_delta >= 0)); + residual_arc_capacity_.Set(arc, free_capacity + capacity_delta); + DCHECK_LE(0, residual_arc_capacity_[arc]); + } else { + // Note that this breaks the preflow invariants but it is currently not an + // issue since we restart from scratch on each Solve() and we set the status + // to NOT_SOLVED. + // + // TODO(user): The easiest is probably to allow negative node excess in + // other places than the source, but the current implementation does not + // deal with this. + SetCapacityAndClearFlow(arc, new_capacity); + } +} + +template +void GenericMaxFlow::GetSourceSideMinCut( + std::vector* result) { + ComputeReachableNodes(source_, result); +} + +template +void GenericMaxFlow::GetSinkSideMinCut(std::vector* result) { + ComputeReachableNodes(sink_, result); +} + +template +bool GenericMaxFlow::CheckResult() const { + SCOPED_TIME_STAT(&stats_); + if (node_excess_[source_] != -node_excess_[sink_]) { + LOG(DFATAL) << "-node_excess_[source_] = " << -node_excess_[source_] + << " != node_excess_[sink_] = " << node_excess_[sink_]; + return false; + } + for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) { + if (node != source_ && node != sink_) { + if (node_excess_[node] != 0) { + LOG(DFATAL) << "node_excess_[" << node << "] = " << node_excess_[node] + << " != 0"; + return false; + } + } + } + for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) { + const ArcIndex opposite = Opposite(arc); + const FlowQuantity direct_capacity = residual_arc_capacity_[arc]; + const FlowQuantity opposite_capacity = residual_arc_capacity_[opposite]; + if (direct_capacity < 0) { + LOG(DFATAL) << "residual_arc_capacity_[" << arc + << "] = " << direct_capacity << " < 0"; + return false; + } + if (opposite_capacity < 0) { + LOG(DFATAL) << "residual_arc_capacity_[" << opposite + << "] = " << opposite_capacity << " < 0"; + return false; + } + // The initial capacity of the direct arcs is non-negative. + if (direct_capacity + opposite_capacity < 0) { + LOG(DFATAL) << "initial capacity [" << arc + << "] = " << direct_capacity + opposite_capacity << " < 0"; + return false; + } + } + + if (GetOptimalFlow() < kMaxFlowQuantity && AugmentingPathExists()) { + LOG(DFATAL) << "The algorithm terminated, but the flow is not maximal!"; + return false; + } + + return true; +} + +template +bool GenericMaxFlow::AugmentingPathExists() const { + SCOPED_TIME_STAT(&stats_); + + // We simply compute the reachability from the source in the residual graph. + const NodeIndex num_nodes = graph_->num_nodes(); + std::vector is_reached(num_nodes, false); + std::vector to_process; + + to_process.push_back(source_); + is_reached[source_] = true; + while (!to_process.empty()) { + const NodeIndex node = to_process.back(); + to_process.pop_back(); + for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); + it.Next()) { + const ArcIndex arc = it.Index(); + if (residual_arc_capacity_[arc] > 0) { + const NodeIndex head = graph_->Head(arc); + if (!is_reached[head]) { + is_reached[head] = true; + to_process.push_back(head); + } + } + } + } + return is_reached[sink_]; +} + +template +bool GenericMaxFlow::CheckRelabelPrecondition(NodeIndex node) const { + DCHECK(IsActive(node)); + for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); + it.Next()) { + const ArcIndex arc = it.Index(); + DCHECK(!IsAdmissible(node, arc, node_potential_.data())) + << DebugString("CheckRelabelPrecondition:", arc); + } + return true; +} + +template +std::string GenericMaxFlow::DebugString(absl::string_view context, + ArcIndex arc) const { + const NodeIndex tail = Tail(arc); + const NodeIndex head = Head(arc); + return absl::StrFormat( + "%s Arc %d, from %d to %d, " + "Capacity = %d, Residual capacity = %d, " + "Flow = residual capacity for reverse arc = %d, " + "Height(tail) = %d, Height(head) = %d, " + "Excess(tail) = %d, Excess(head) = %d", + context, arc, tail, head, Capacity(arc), residual_arc_capacity_[arc], + Flow(arc), node_potential_[tail], node_potential_[head], + node_excess_[tail], node_excess_[head]); +} + +template +bool GenericMaxFlow::Solve() { + status_ = NOT_SOLVED; + InitializePreflow(); + + // Deal with the case when source_ or sink_ is not inside graph_. + // Since they are both specified independently of the graph, we do need to + // take care of this corner case. + const NodeIndex num_nodes = graph_->num_nodes(); + if (sink_ >= num_nodes || source_ >= num_nodes) { + // Behave like a normal graph where source_ and sink_ are disconnected. + // Note that the arc flow is set to 0 by InitializePreflow(). + status_ = OPTIMAL; + return true; + } + + RefineWithGlobalUpdate(); + + status_ = OPTIMAL; + DCHECK(CheckResult()); + + if (GetOptimalFlow() == kMaxFlowQuantity && AugmentingPathExists()) { + // In this case, we are sure that the flow is > kMaxFlowQuantity. + status_ = INT_OVERFLOW; + } + IF_STATS_ENABLED(VLOG(1) << stats_.StatString()); + return true; +} + +template +void GenericMaxFlow::InitializePreflow() { + SCOPED_TIME_STAT(&stats_); + + // TODO(user): Ebert graph has an issue with nodes with no arcs, so we + // use max_num_nodes here to resize vectors. + const NodeIndex num_nodes = graph_->num_nodes(); + const NodeIndex max_num_nodes = Graphs::NodeReservation(*graph_); + const ArcIndex num_arcs = graph_->num_arcs(); + + // InitializePreflow() clears the whole flow that could have been computed + // by a previous Solve(). This is not optimal in terms of complexity. + // + // TODO(user): find a way to make the re-solving incremental (not an obvious + // task, and there has not been a lot of literature on the subject.) + node_excess_.assign(max_num_nodes, 0); + for (ArcIndex arc = 0; arc < num_arcs; ++arc) { + SetCapacityAndClearFlow(arc, Capacity(arc)); + } + + // All the initial heights are zero except for the source whose height is + // equal to the number of nodes and will never change during the algorithm. + node_potential_.assign(max_num_nodes, 0); + node_potential_[source_] = num_nodes; + + // Initially no arcs are admissible except maybe the one leaving the source, + // but we treat the source in a special way, see + // SaturateOutgoingArcsFromSource(). + first_admissible_arc_.assign(max_num_nodes, Graph::kNilArc); +} + +// Note(user): Calling this function will break the property on the node +// potentials because of the way we cancel flow on cycle. However, we only call +// that at the end of the algorithm, or just before a GlobalUpdate() that will +// restore the precondition on the node potentials. +template +void GenericMaxFlow::PushFlowExcessBackToSource() { + SCOPED_TIME_STAT(&stats_); + const NodeIndex num_nodes = graph_->num_nodes(); + + // We implement a variation of Tarjan's strongly connected component algorithm + // to detect cycles published in: Tarjan, R. E. (1972), "Depth-first search + // and linear graph algorithms", SIAM Journal on Computing. A description can + // also be found in wikipedia. + + // Stored nodes are settled nodes already stored in the + // reverse_topological_order (except the sink_ that we do not actually store). + std::vector stored(num_nodes, false); + stored[sink_] = true; + + // The visited nodes that are not yet stored are all the nodes from the + // source_ to the current node in the current dfs branch. + std::vector visited(num_nodes, false); + visited[sink_] = true; + + // Stack of arcs to explore in the dfs search. + // The current node is Head(arc_stack.back()). + std::vector arc_stack; + + // Increasing list of indices into the arc_stack that correspond to the list + // of arcs in the current dfs branch from the source_ to the current node. + std::vector index_branch; + + // Node in reverse_topological_order in the final dfs tree. + std::vector reverse_topological_order; + + // We start by pushing all the outgoing arcs from the source on the stack to + // avoid special conditions in the code. As a result, source_ will not be + // stored in reverse_topological_order, and this is what we want. + for (OutgoingArcIterator it(*graph_, source_); it.Ok(); it.Next()) { + const ArcIndex arc = it.Index(); + const FlowQuantity flow = Flow(arc); + if (flow > 0) { + arc_stack.push_back(arc); + } + } + visited[source_] = true; + + // Start the dfs on the subgraph formed by the direct arcs with positive flow. + while (!arc_stack.empty()) { + const NodeIndex node = Head(arc_stack.back()); + + // If the node is visited, it means we have explored all its arcs and we + // have just backtracked in the dfs. Store it if it is not already stored + // and process the next arc on the stack. + if (visited[node]) { + if (!stored[node]) { + stored[node] = true; + reverse_topological_order.push_back(node); + DCHECK(!index_branch.empty()); + index_branch.pop_back(); + } + arc_stack.pop_back(); + continue; + } + + // The node is a new unexplored node, add all its outgoing arcs with + // positive flow to the stack and go deeper in the dfs. + DCHECK(!stored[node]); + DCHECK(index_branch.empty() || + (arc_stack.size() - 1 > index_branch.back())); + visited[node] = true; + index_branch.push_back(arc_stack.size() - 1); + + for (OutgoingArcIterator it(*graph_, node); it.Ok(); it.Next()) { + const ArcIndex arc = it.Index(); + const FlowQuantity flow = Flow(arc); + const NodeIndex head = Head(arc); + if (flow > 0 && !stored[head]) { + if (!visited[head]) { + arc_stack.push_back(arc); + } else { + // There is a cycle. + // Find the first index to consider, + // arc_stack[index_branch[cycle_begin]] will be the first arc on the + // cycle. + int cycle_begin = index_branch.size(); + while (cycle_begin > 0 && + Head(arc_stack[index_branch[cycle_begin - 1]]) != head) { + --cycle_begin; + } + + // Compute the maximum flow that can be canceled on the cycle and the + // min index such that arc_stack[index_branch[i]] will be saturated. + FlowQuantity max_flow = flow; + int first_saturated_index = index_branch.size(); + for (int i = index_branch.size() - 1; i >= cycle_begin; --i) { + const ArcIndex arc_on_cycle = arc_stack[index_branch[i]]; + if (Flow(arc_on_cycle) <= max_flow) { + max_flow = Flow(arc_on_cycle); + first_saturated_index = i; + } + } + + // This is just here for a DCHECK() below. + const FlowQuantity excess = node_excess_[head]; + + // Cancel the flow on the cycle, and set visited[node] = false for + // the node that will be backtracked over. + PushFlow(-max_flow, node, arc); + for (int i = index_branch.size() - 1; i >= cycle_begin; --i) { + const ArcIndex arc_on_cycle = arc_stack[index_branch[i]]; + PushFlow(-max_flow, Tail(arc_on_cycle), arc_on_cycle); + if (i >= first_saturated_index) { + DCHECK(visited[Head(arc_on_cycle)]); + visited[Head(arc_on_cycle)] = false; + } else { + DCHECK_GT(Flow(arc_on_cycle), 0); + } + } + + // This is a simple check that the flow was pushed properly. + DCHECK_EQ(excess, node_excess_[head]); + + // Backtrack the dfs just before index_branch[first_saturated_index]. + // If the current node is still active, there is nothing to do. + if (first_saturated_index < index_branch.size()) { + arc_stack.resize(index_branch[first_saturated_index]); + index_branch.resize(first_saturated_index); + + // We backtracked over the current node, so there is no need to + // continue looping over its arcs. + break; + } + } + } + } + } + DCHECK(arc_stack.empty()); + DCHECK(index_branch.empty()); + + // Return the flow to the sink. Note that the sink_ and the source_ are not + // stored in reverse_topological_order. + for (int i = 0; i < reverse_topological_order.size(); i++) { + const NodeIndex node = reverse_topological_order[i]; + if (node_excess_[node] == 0) continue; + for (IncomingArcIterator it(*graph_, node); it.Ok(); it.Next()) { + const ArcIndex opposite_arc = Opposite(it.Index()); + if (residual_arc_capacity_[opposite_arc] > 0) { + const FlowQuantity flow = + std::min(node_excess_[node], residual_arc_capacity_[opposite_arc]); + PushFlow(flow, node, opposite_arc); + if (node_excess_[node] == 0) break; + } + } + DCHECK_EQ(0, node_excess_[node]); + } + DCHECK_EQ(-node_excess_[source_], node_excess_[sink_]); +} + +template +void GenericMaxFlow::GlobalUpdate() { + SCOPED_TIME_STAT(&stats_); + bfs_queue_.clear(); + int queue_index = 0; + const NodeIndex num_nodes = graph_->num_nodes(); + node_in_bfs_queue_.assign(num_nodes, false); + node_in_bfs_queue_[sink_] = true; + + // We cache this as this is a hot-loop and we don't want bound checking. + FlowQuantity* const node_excess = node_excess_.data(); + + // We do a BFS in the reverse residual graph, starting from the sink. + // Because all the arcs from the source are saturated (except in + // presence of integer overflow), the source cannot reach the sink in the + // residual graph. If there is possible overflow and the source is reachable, + // we still do not want to relabel it, so we start with the source marked. + node_in_bfs_queue_[source_] = true; + bfs_queue_.push_back(sink_); + + while (queue_index != bfs_queue_.size()) { + const NodeIndex node = bfs_queue_[queue_index]; + ++queue_index; + const NodeIndex candidate_distance = node_potential_[node] + 1; + for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); + it.Next()) { + const ArcIndex arc = it.Index(); + const NodeIndex head = Head(arc); + + // Skip the arc if the height of head was already set to the correct + // value (Remember we are doing reverse BFS). + if (node_in_bfs_queue_[head]) continue; + + // TODO(user): By using more memory we can speed this up quite a bit by + // avoiding to take the opposite arc here, too options: + // - if (residual_arc_capacity_[arc] != arc_capacity_[arc]) + // - if (opposite_arc_is_admissible_[arc]) // need updates. + // Experiment with the first option shows more than 10% gain on this + // function running time, which is the bottleneck on many instances. + const ArcIndex opposite_arc = Opposite(arc); + if (residual_arc_capacity_[opposite_arc] > 0) { + // Note(user): We used to have a DCHECK_GE(candidate_distance, + // node_potential_[head]); which is always true except in the case + // where we can push more than kMaxFlowQuantity out of the source. The + // problem comes from the fact that in this case, we call + // PushFlowExcessBackToSource() in the middle of the algorithm. The + // later call will break the properties of the node potential. Note + // however, that this function will recompute a good node potential + // for all the nodes and thus fix the issue. + + // If head is active, we can steal some or all of its excess. + // This brings a huge gain on some problems. + // Note(user): I haven't seen this anywhere in the literature. + // TODO(user): Investigate more and maybe write a publication :) + if (node_excess[head] > 0) { + const FlowQuantity flow = + std::min(node_excess[head], residual_arc_capacity_[opposite_arc]); + FastPushFlow(flow, head, opposite_arc, node_excess); + + // If the arc became saturated, it is no longer in the residual + // graph, so we do not need to consider head at this time. + if (residual_arc_capacity_[opposite_arc] == 0) continue; + } + + // Note that there is no need to touch first_admissible_arc_[node] + // because of the relaxed Relabel() we use. + node_potential_[head] = candidate_distance; + node_in_bfs_queue_[head] = true; + bfs_queue_.push_back(head); + } + } + } + + // At the end of the search, some nodes may not be in the bfs_queue_. Such + // nodes cannot reach the sink_ or source_ in the residual graph, so there is + // no point trying to push flow toward them. We obtain this effect by setting + // their height to something unreachable. + // + // Note that this also prevents cycling due to our anti-overflow procedure. + // For instance, suppose there is an edge s -> n outgoing from the source. If + // node n has no other connection and some excess, we will push the flow back + // to the source, but if we don't update the height of n + // SaturateOutgoingArcsFromSource() will push the flow to n again. + // TODO(user): This is another argument for another anti-overflow algorithm. + for (NodeIndex node = 0; node < num_nodes; ++node) { + if (!node_in_bfs_queue_[node]) { + node_potential_[node] = 2 * num_nodes - 1; + } + } + + // Reset the active nodes. Doing it like this pushes the nodes in increasing + // order of height. Note that bfs_queue_[0] is the sink_ so we skip it. + DCHECK(IsEmptyActiveNodeContainer()); + for (int i = 1; i < bfs_queue_.size(); ++i) { + const NodeIndex node = bfs_queue_[i]; + if (node_excess_[node] > 0) { + DCHECK(IsActive(node)); + PushActiveNode(node); + } + } +} + +template +bool GenericMaxFlow::SaturateOutgoingArcsFromSource() { + SCOPED_TIME_STAT(&stats_); + const NodeIndex num_nodes = graph_->num_nodes(); + + // If sink_ or source_ already have kMaxFlowQuantity, then there is no + // point pushing more flow since it will cause an integer overflow. + if (node_excess_[sink_] == kMaxFlowQuantity) return false; + if (node_excess_[source_] == -kMaxFlowQuantity) return false; + + bool flow_pushed = false; + for (OutgoingArcIterator it(*graph_, source_); it.Ok(); it.Next()) { + const ArcIndex arc = it.Index(); + const FlowQuantity flow = residual_arc_capacity_[arc]; + + // This is a special IsAdmissible() condition for the source. + if (flow == 0 || node_potential_[Head(arc)] >= num_nodes) continue; + + // We are careful in case the sum of the flow out of the source is greater + // than kMaxFlowQuantity to avoid overflow. + const FlowQuantity current_flow_out_of_source = -node_excess_[source_]; + DCHECK_GE(flow, 0) << flow; + DCHECK_GE(current_flow_out_of_source, 0) << current_flow_out_of_source; + const FlowQuantity capped_flow = + kMaxFlowQuantity - current_flow_out_of_source; + if (capped_flow < flow) { + // We push as much flow as we can so the current flow on the network will + // be kMaxFlowQuantity. + + // Since at the beginning of the function, current_flow_out_of_source + // was different from kMaxFlowQuantity, we are sure to have pushed some + // flow before if capped_flow is 0. + if (capped_flow == 0) return true; + PushFlow(capped_flow, source_, arc); + return true; + } + PushFlow(flow, source_, arc); + flow_pushed = true; + } + DCHECK_LE(node_excess_[source_], 0); + return flow_pushed; +} + +template +void GenericMaxFlow::PushFlow(FlowQuantity flow, NodeIndex tail, + ArcIndex arc) { + SCOPED_TIME_STAT(&stats_); + + // Update the residual capacity of the arc and its opposite arc. + DCHECK_NE(flow, 0); + residual_arc_capacity_[arc] -= flow; + residual_arc_capacity_[Opposite(arc)] += flow; + DCHECK_GE(residual_arc_capacity_[arc], 0); + DCHECK_GE(residual_arc_capacity_[Opposite(arc)], 0); + + // Update the excesses at the tail and head of the arc. + // + // Note(user): node_excess_ should be always greater than or equal to 0 except + // for the source where it should always be smaller than or equal to 0. Note + // however that we cannot check this because when we cancel the flow on a + // cycle in PushFlowExcessBackToSource(), we may break this invariant during + // the operation even if it is still valid at the end. + node_excess_[tail] -= flow; + node_excess_[Head(arc)] += flow; +} + +template +void GenericMaxFlow::FastPushFlow(FlowQuantity flow, NodeIndex tail, + ArcIndex arc, + FlowQuantity* node_excess) { + SCOPED_TIME_STAT(&stats_); + + DCHECK_NE(flow, 0); + residual_arc_capacity_[arc] -= flow; + residual_arc_capacity_[Opposite(arc)] += flow; + DCHECK_GE(residual_arc_capacity_[arc], 0); + DCHECK_GE(residual_arc_capacity_[Opposite(arc)], 0); + + node_excess[tail] -= flow; + node_excess[Head(arc)] += flow; +} + +template +void GenericMaxFlow::InitializeActiveNodeContainer() { + SCOPED_TIME_STAT(&stats_); + DCHECK(IsEmptyActiveNodeContainer()); + const NodeIndex num_nodes = graph_->num_nodes(); + for (NodeIndex node = 0; node < num_nodes; ++node) { + if (IsActive(node)) { + // This node cannot reach the sink in the residual graph, we don't need to + // consider it anymore, and we will just push its potential excess back to + // the source in PushFlowExcessBackToSource() at the end. + if (node_potential_[node] >= num_nodes) continue; + PushActiveNode(node); + } + } +} + +template +void GenericMaxFlow::RefineWithGlobalUpdate() { + SCOPED_TIME_STAT(&stats_); + + // TODO(user): This should be graph_->num_nodes(), but ebert graph does not + // have a correct size if the highest index nodes have no arcs. + const NodeIndex num_nodes = Graphs::NodeReservation(*graph_); + std::vector skip_active_node; + + // Usually SaturateOutgoingArcsFromSource() will saturate all the arcs from + // the source in one go, and we will loop just once. But in case we can push + // more than kMaxFlowQuantity out of the source the loop is as follow: + // - Push up to kMaxFlowQuantity out of the source on the admissible outgoing + // arcs. Stop if no flow was pushed. + // - Compute the current max-flow. This will push some flow back to the + // source and render more outgoing arcs from the source not admissible. + // + // TODO(user): This may not be the most efficient algorithm if we need to loop + // many times. An alternative may be to handle the source like the other nodes + // in the algorithm, initially putting an excess of kMaxFlowQuantity on it, + // and making the source active like any other node with positive excess. To + // investigate. + while (SaturateOutgoingArcsFromSource()) { + int num_skipped; + do { + num_skipped = 0; + skip_active_node.assign(num_nodes, 0); + skip_active_node[sink_] = 2; + skip_active_node[source_] = 2; + GlobalUpdate(); + while (!IsEmptyActiveNodeContainer()) { + const NodeIndex node = GetAndRemoveFirstActiveNode(); + if (skip_active_node[node] > 1) { + if (node != sink_ && node != source_) ++num_skipped; + continue; + } + const NodeIndex old_height = node_potential_[node]; + Discharge(node); + + // The idea behind this is that if a node height augments by more than + // one, then it is likely to push flow back the way it came. This can + // lead to very costly loops. A bad case is: source -> n1 -> n2 and n2 + // just recently isolated from the sink. Then n2 will push flow back to + // n1, and n1 to n2 and so on. The height of each node will increase by + // steps of two until the height of the source is reached, which can + // take a long time. If the chain is longer, the situation is even + // worse. The behavior of this heuristic is related to the Gap + // heuristic. + // + // Note that the global update will fix all such cases efficiently. So + // the idea is to discharge the active node as much as possible, and + // then do a global update. + // + // We skip a node when this condition was true 2 times to avoid doing a + // global update too frequently. + if (node_potential_[node] > old_height + 1) { + ++skip_active_node[node]; + } + } + } while (num_skipped > 0); + + // We use a two-phase algorithm: + // 1/ Only deal with nodes that can reach the sink. At the end we know the + // value of the maximum flow and we have a min-cut. + // 2/ Call PushFlowExcessBackToSource() to obtain a max-flow. This is + // usually a lot faster than the first phase. + PushFlowExcessBackToSource(); + } +} + +template +void GenericMaxFlow::Discharge(const NodeIndex node) { + SCOPED_TIME_STAT(&stats_); + const NodeIndex num_nodes = graph_->num_nodes(); + + // We cache this as this is a hot-loop and we don't want bound checking. + FlowQuantity* const node_excess = node_excess_.data(); + NodeHeight* const node_potentials = node_potential_.data(); + ArcIndex* const first_admissible_arc = first_admissible_arc_.data(); + + while (true) { + DCHECK(IsActive(node)); + for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node, + first_admissible_arc[node]); + it.Ok(); it.Next()) { + const ArcIndex arc = it.Index(); + if (IsAdmissible(node, arc, node_potentials)) { + DCHECK(IsActive(node)); + const NodeIndex head = Head(arc); + if (node_excess[head] == 0) { + // The push below will make the node active for sure. Note that we may + // push the sink_, but that is handled properly in Refine(). + PushActiveNode(head); + } + const FlowQuantity delta = + std::min(node_excess[node], residual_arc_capacity_[arc]); + FastPushFlow(delta, node, arc, node_excess); + if (node_excess[node] == 0) { + first_admissible_arc[node] = arc; // arc may still be admissible. + return; + } + } + } + Relabel(node); + + // This node can no longer reach the sink, skip until + // PushFlowExcessBackToSource(). + if (node_potentials[node] >= num_nodes) break; + } +} + +template +void GenericMaxFlow::Relabel(NodeIndex node) { + SCOPED_TIME_STAT(&stats_); + // Because we use a relaxed version, this is no longer true if the + // first_admissible_arc_[node] was not actually the first arc! + // DCHECK(CheckRelabelPrecondition(node)); + NodeHeight min_height = std::numeric_limits::max(); + ArcIndex first_admissible_arc = Graph::kNilArc; + for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); + it.Next()) { + const ArcIndex arc = it.Index(); + if (residual_arc_capacity_[arc] > 0) { + // Update min_height only for arcs with available capacity. + NodeHeight head_height = node_potential_[Head(arc)]; + if (head_height < min_height) { + min_height = head_height; + first_admissible_arc = arc; + + // We found an admissible arc at the current height, just stop there. + // This is the true first_admissible_arc_[node]. + if (min_height + 1 == node_potential_[node]) break; + } + } + } + DCHECK_NE(first_admissible_arc, Graph::kNilArc); + node_potential_[node] = min_height + 1; + + // Note that after a Relabel(), the loop will continue in Discharge(), and + // we are sure that all the arcs before first_admissible_arc are not + // admissible since their height is > min_height. + first_admissible_arc_[node] = first_admissible_arc; +} + +template +typename Graph::ArcIndex GenericMaxFlow::Opposite(ArcIndex arc) const { + return Graphs::OppositeArc(*graph_, arc); +} + +template +bool GenericMaxFlow::IsArcDirect(ArcIndex arc) const { + return IsArcValid(arc) && arc >= 0; +} + +template +bool GenericMaxFlow::IsArcValid(ArcIndex arc) const { + return Graphs::IsArcValid(*graph_, arc); +} + +template +template +void GenericMaxFlow::ComputeReachableNodes( + NodeIndex start, std::vector* result) { + // If start is not a valid node index, it can reach only itself. + // Note(user): This is needed because source and sink are given independently + // of the graph and sometimes before it is even constructed. + const NodeIndex num_nodes = graph_->num_nodes(); + if (start >= num_nodes) { + result->clear(); + result->push_back(start); + return; + } + bfs_queue_.clear(); + node_in_bfs_queue_.assign(num_nodes, false); + + int queue_index = 0; + bfs_queue_.push_back(start); + node_in_bfs_queue_[start] = true; + while (queue_index != bfs_queue_.size()) { + const NodeIndex node = bfs_queue_[queue_index]; + ++queue_index; + for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); + it.Next()) { + const ArcIndex arc = it.Index(); + const NodeIndex head = Head(arc); + if (node_in_bfs_queue_[head]) continue; + if (residual_arc_capacity_[reverse ? Opposite(arc) : arc] == 0) continue; + node_in_bfs_queue_[head] = true; + bfs_queue_.push_back(head); + } + } + *result = bfs_queue_; +} + +template +FlowModelProto GenericMaxFlow::CreateFlowModel() { + FlowModelProto model; + model.set_problem_type(FlowModelProto::MAX_FLOW); + for (int n = 0; n < graph_->num_nodes(); ++n) { + FlowNodeProto* node = model.add_nodes(); + node->set_id(n); + if (n == source_) node->set_supply(1); + if (n == sink_) node->set_supply(-1); + } + for (int a = 0; a < graph_->num_arcs(); ++a) { + FlowArcProto* arc = model.add_arcs(); + arc->set_tail(graph_->Tail(a)); + arc->set_head(graph_->Head(a)); + arc->set_capacity(Capacity(a)); + } + return model; +} + +} // namespace operations_research + +#endif // OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_ diff --git a/ortools/graph/generic_max_flow_test.cc b/ortools/graph/generic_max_flow_test.cc new file mode 100644 index 0000000000..8f3fb5982c --- /dev/null +++ b/ortools/graph/generic_max_flow_test.cc @@ -0,0 +1,685 @@ +// Copyright 2010-2024 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/graph/generic_max_flow.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "absl/random/random.h" +#include "absl/strings/str_format.h" +#include "absl/types/span.h" +#include "benchmark/benchmark.h" +#include "gtest/gtest.h" +#include "ortools/base/gmock.h" +#include "ortools/base/logging.h" +#include "ortools/graph/ebert_graph.h" +#include "ortools/graph/graph.h" +#include "ortools/graph/graphs.h" +#include "ortools/linear_solver/linear_solver.h" + +namespace operations_research { +namespace { + +using ::testing::ContainerEq; +using ::testing::WhenSorted; + +template +typename GenericMaxFlow::Status MaxFlowTester( + const typename Graph::NodeIndex num_nodes, + const typename Graph::ArcIndex num_arcs, + const typename Graph::NodeIndex* tail, + const typename Graph::NodeIndex* head, const FlowQuantity* capacity, + const FlowQuantity* expected_flow, const FlowQuantity expected_total_flow, + const std::vector* expected_source_min_cut = + nullptr, + const std::vector* expected_sink_min_cut = + nullptr) { + Graph graph(num_nodes, num_arcs); + for (int i = 0; i < num_arcs; ++i) { + graph.AddArc(tail[i], head[i]); + } + std::vector permutation; + Graphs::Build(&graph, &permutation); + EXPECT_TRUE(permutation.empty()); + + GenericMaxFlow max_flow(&graph, 0, num_nodes - 1); + for (typename Graph::ArcIndex arc = 0; arc < num_arcs; ++arc) { + max_flow.SetArcCapacity(arc, capacity[arc]); + EXPECT_EQ(max_flow.Capacity(arc), capacity[arc]); + } + EXPECT_TRUE(max_flow.Solve()); + if (max_flow.status() == GenericMaxFlow::OPTIMAL) { + const FlowQuantity total_flow = max_flow.GetOptimalFlow(); + EXPECT_EQ(expected_total_flow, total_flow); + for (int i = 0; i < num_arcs; ++i) { + EXPECT_EQ(expected_flow[i], max_flow.Flow(i)) << " i = " << i; + } + } + + // Tests the min-cut functions. + if (expected_source_min_cut != nullptr) { + std::vector cut; + max_flow.GetSourceSideMinCut(&cut); + std::sort(cut.begin(), cut.end()); + EXPECT_THAT(*expected_source_min_cut, WhenSorted(ContainerEq(cut))); + } + if (expected_sink_min_cut != nullptr) { + std::vector cut; + max_flow.GetSinkSideMinCut(&cut); + std::sort(cut.begin(), cut.end()); + EXPECT_THAT(*expected_sink_min_cut, WhenSorted(ContainerEq(cut))); + } + + return max_flow.status(); +} + +template +class GenericMaxFlowTest : public ::testing::Test {}; + +typedef ::testing::Types, + util::ReverseArcStaticGraph<>, + util::ReverseArcMixedGraph<> > + GraphTypes; + +TYPED_TEST_SUITE(GenericMaxFlowTest, GraphTypes); + +TYPED_TEST(GenericMaxFlowTest, FeasibleFlow1) { + const int kNumNodes = 4; + const int kNumArcs = 3; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 1, 2}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3}; + const FlowQuantity kCapacity[kNumArcs] = {8, 10, 8}; + const FlowQuantity kExpectedFlow[kNumArcs] = {8, 8, 8}; + const FlowQuantity kExpectedTotalFlow = 8; + std::vector source_cut({0}); + std::vector sink_cut({3}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +TYPED_TEST(GenericMaxFlowTest, FeasibleFlow2) { + const int kNumNodes = 6; + const int kNumArcs = 9; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 0, 0, 1, + 2, 3, 3, 4}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 4, 3, + 4, 4, 5, 5}; + const FlowQuantity kCapacity[kNumArcs] = {6, 8, 5, 0, 1, 4, 0, 6, 4}; + const FlowQuantity kExpectedFlow[kNumArcs] = {1, 4, 5, 0, 1, 4, 0, 6, 4}; + const FlowQuantity kExpectedTotalFlow = 10; + std::vector source_cut({0, 1, 2}); + std::vector sink_cut({5}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +TYPED_TEST(GenericMaxFlowTest, FeasibleFlowWithMultipleArcs) { + const int kNumNodes = 5; + const int kNumArcs = 8; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 1, + 2, 2, 3, 3}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 1, 2, 2, + 3, 3, 4, 4}; + const FlowQuantity kCapacity[kNumArcs] = {5, 3, 5, 3, 4, 4, 4, 4}; + const FlowQuantity kExpectedFlow[kNumArcs] = {5, 3, 5, 3, 4, 4, 4, 4}; + const FlowQuantity kExpectedTotalFlow = 8; + std::vector source_cut({0}); + std::vector sink_cut({4}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +TYPED_TEST(GenericMaxFlowTest, HugeCapacity) { + const FlowQuantity kCapacityMax = std::numeric_limits::max(); + const int kNumNodes = 5; + const int kNumArcs = 5; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 2, 3}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 3, 4}; + const FlowQuantity kCapacity[kNumArcs] = {kCapacityMax, kCapacityMax, 5, 3, + kCapacityMax}; + const FlowQuantity kExpectedFlow[kNumArcs] = {5, 3, 5, 3, 8}; + const FlowQuantity kExpectedTotalFlow = 8; + std::vector source_cut({0, 1, 2}); + std::vector sink_cut({4, 3}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +TYPED_TEST(GenericMaxFlowTest, FlowQuantityOverflowLimitCase) { + const FlowQuantity kCapacityMax = std::numeric_limits::max(); + const FlowQuantity kHalfLow = kCapacityMax / 2; + const FlowQuantity kHalfHigh = kCapacityMax - kHalfLow; + const int kNumNodes = 5; + const int kNumArcs = 5; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 2, 3}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 3, 4}; + const FlowQuantity kCapacity[kNumArcs] = {kCapacityMax, kCapacityMax, + kHalfLow, kHalfHigh, kCapacityMax}; + const FlowQuantity kExpectedFlow[kNumArcs] = {kHalfLow, kHalfHigh, kHalfLow, + kHalfHigh, kCapacityMax}; + const FlowQuantity kExpectedTotalFlow = kCapacityMax; + std::vector source_cut({0, 1, 2}); + std::vector sink_cut({4}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +TYPED_TEST(GenericMaxFlowTest, FlowQuantityOverflow) { + const FlowQuantity kCapacityMax = std::numeric_limits::max(); + const int kNumNodes = 4; + const int kNumArcs = 4; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 2}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 3}; + const FlowQuantity kCapacity[kNumArcs] = {kCapacityMax, kCapacityMax, + kCapacityMax, kCapacityMax}; + const FlowQuantity kExpectedFlow[kNumArcs] = {kCapacityMax, kCapacityMax, + kCapacityMax, kCapacityMax}; + const FlowQuantity kExpectedTotalFlow = kCapacityMax; + EXPECT_EQ( + GenericMaxFlow::INT_OVERFLOW, + MaxFlowTester(kNumNodes, kNumArcs, kTail, kHead, kCapacity, + kExpectedFlow, kExpectedTotalFlow)); +} + +TYPED_TEST(GenericMaxFlowTest, DirectArcFromSourceToSink) { + const int kNumNodes = 4; + const int kNumArcs = 5; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 0, 1, 2}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 3, 2, 3, 3}; + const FlowQuantity kCapacity[kNumArcs] = {5, 8, 5, 2, 2}; + const FlowQuantity kExpectedFlow[kNumArcs] = {2, 8, 2, 2, 2}; + const FlowQuantity kExpectedTotalFlow = 12; + std::vector source_cut({0, 1, 2}); + std::vector sink_cut({3}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +TYPED_TEST(GenericMaxFlowTest, FlowOnDisconnectedGraph1) { + const int kNumNodes = 6; + const int kNumArcs = 7; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 0, 0, 1, 2, 3}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 4, 3, 4, 4}; + const FlowQuantity kCapacity[kNumArcs] = {5, 8, 5, 3, 4, 5, 6}; + const FlowQuantity kExpectedFlow[kNumArcs] = {0, 0, 0, 0, 0, 0, 0}; + const FlowQuantity kExpectedTotalFlow = 0; + std::vector source_cut({0, 1, 2, 3, 4}); + std::vector sink_cut({5}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +TYPED_TEST(GenericMaxFlowTest, FlowOnDisconnectedGraph2) { + const int kNumNodes = 6; + const int kNumArcs = 5; + const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 3, 3, 4}; + const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 4, 5, 5}; + const FlowQuantity kCapacity[kNumArcs] = {5, 8, 6, 6, 4}; + const FlowQuantity kExpectedFlow[kNumArcs] = {0, 0, 0, 0, 0}; + const FlowQuantity kExpectedTotalFlow = 0; + std::vector source_cut({0, 1, 2}); + std::vector sink_cut({3, 4, 5}); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, + MaxFlowTester( + kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, + kExpectedTotalFlow, &source_cut, &sink_cut)); +} + +template +void AddSourceAndSink(const typename Graph::NodeIndex num_tails, + const typename Graph::NodeIndex num_heads, Graph* graph) { + const typename Graph::NodeIndex source = num_tails + num_heads; + const typename Graph::NodeIndex sink = num_tails + num_heads + 1; + for (typename Graph::NodeIndex tail = 0; tail < num_tails; ++tail) { + graph->AddArc(source, tail); + } + for (typename Graph::NodeIndex head = 0; head < num_heads; ++head) { + graph->AddArc(num_tails + head, sink); + } +} + +template +void GenerateCompleteGraph(const typename Graph::NodeIndex num_tails, + const typename Graph::NodeIndex num_heads, + Graph* graph) { + const typename Graph::NodeIndex num_nodes = num_tails + num_heads + 2; + const typename Graph::ArcIndex num_arcs = + num_tails * num_heads + num_tails + num_heads; + graph->Reserve(num_nodes, num_arcs); + for (typename Graph::NodeIndex tail = 0; tail < num_tails; ++tail) { + for (typename Graph::NodeIndex head = 0; head < num_heads; ++head) { + graph->AddArc(tail, head + num_tails); + } + } + AddSourceAndSink(num_tails, num_heads, graph); +} + +template +void GeneratePartialRandomGraph(const typename Graph::NodeIndex num_tails, + const typename Graph::NodeIndex num_heads, + const typename Graph::NodeIndex degree, + Graph* graph) { + const typename Graph::NodeIndex num_nodes = num_tails + num_heads + 2; + const typename Graph::ArcIndex num_arcs = + num_tails * degree + num_tails + num_heads; + graph->Reserve(num_nodes, num_arcs); + std::mt19937 randomizer(0); + for (typename Graph::NodeIndex tail = 0; tail < num_tails; ++tail) { + for (typename Graph::NodeIndex d = 0; d < degree; ++d) { + const typename Graph::NodeIndex head = + absl::Uniform(randomizer, 0, num_heads); + graph->AddArc(tail, head + num_tails); + } + } + AddSourceAndSink(num_tails, num_heads, graph); +} + +template +void GenerateRandomArcValuations(const Graph& graph, const int64_t max_range, + std::vector* arc_valuation) { + const typename Graph::ArcIndex num_arcs = graph.num_arcs(); + arc_valuation->resize(num_arcs); + std::mt19937 randomizer(0); + for (typename Graph::ArcIndex arc = 0; arc < graph.num_arcs(); ++arc) { + (*arc_valuation)[arc] = absl::Uniform(randomizer, 0, max_range); + } +} + +template +void SetUpNetworkData(absl::Span arc_capacity, + GenericMaxFlow* max_flow) { + const Graph* graph = max_flow->graph(); + for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { + max_flow->SetArcCapacity(arc, arc_capacity[arc]); + } +} + +template +FlowQuantity SolveMaxFlow(GenericMaxFlow* max_flow) { + EXPECT_TRUE(max_flow->Solve()); + EXPECT_EQ(GenericMaxFlow::OPTIMAL, max_flow->status()); + const Graph* graph = max_flow->graph(); + for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { + EXPECT_LE(max_flow->Flow(Graphs::OppositeArc(*graph, arc)), 0) + << arc; + EXPECT_EQ(0, max_flow->Capacity(Graphs::OppositeArc(*graph, arc))) + << arc; + EXPECT_LE(0, max_flow->Flow(arc)) << arc; + EXPECT_LE(max_flow->Flow(arc), max_flow->Capacity(arc)) << arc; + } + return max_flow->GetOptimalFlow(); +} + +template +FlowQuantity SolveMaxFlowWithLP(GenericMaxFlow* max_flow) { + MPSolver solver("LPSolver", MPSolver::GLOP_LINEAR_PROGRAMMING); + const double infinity = solver.infinity(); + const Graph* graph = max_flow->graph(); + const typename Graph::NodeIndex num_nodes = graph->num_nodes(); + const typename Graph::ArcIndex num_arcs = graph->num_arcs(); + const typename Graph::NodeIndex source_index = num_nodes - 2; + std::unique_ptr constraint(new MPConstraint*[num_nodes]); + for (typename Graph::NodeIndex node = 0; node < graph->num_nodes(); ++node) { + constraint[node] = solver.MakeRowConstraint(); + if (node < source_index) { // Node is neither source nor sink. + constraint[node]->SetBounds(0.0, 0.0); // Flow is conserved. + } else { + constraint[node]->SetBounds(-infinity, +infinity); + } + } + std::unique_ptr var(new MPVariable*[num_arcs]); + for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { + var[arc] = solver.MakeNumVar(0.0, max_flow->Capacity(arc), + absl::StrFormat("v%d", arc)); + constraint[graph->Tail(arc)]->SetCoefficient(var[arc], 1.0); + constraint[graph->Head(arc)]->SetCoefficient(var[arc], -1.0); + } + MPObjective* const objective = solver.MutableObjective(); + for (typename Graph::OutgoingArcIterator arc_it(*graph, source_index); + arc_it.Ok(); arc_it.Next()) { + const typename Graph::ArcIndex arc = arc_it.Index(); + objective->SetCoefficient(var[arc], -1.0); + } + solver.Solve(); + return static_cast(-objective->Value() + .5); +} + +template +struct MaxFlowSolver { + typedef FlowQuantity (*Solver)(GenericMaxFlow* max_flow); +}; + +template +void FullRandomAssignment(typename MaxFlowSolver::Solver f, + typename Graph::NodeIndex num_tails, + typename Graph::NodeIndex num_heads, + FlowQuantity expected_flow1, + FlowQuantity expected_flow2) { + Graph graph; + GenerateCompleteGraph(num_tails, num_heads, &graph); + Graphs::Build(&graph); + std::vector arc_capacity(graph.num_arcs(), 1); + std::unique_ptr > max_flow(new GenericMaxFlow( + &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); + SetUpNetworkData(arc_capacity, max_flow.get()); + FlowQuantity flow = f(max_flow.get()); + EXPECT_EQ(expected_flow1, flow); +} + +template +void PartialRandomAssignment(typename MaxFlowSolver::Solver f, + typename Graph::NodeIndex num_tails, + typename Graph::NodeIndex num_heads, + FlowQuantity expected_flow1, + FlowQuantity expected_flow2) { + const typename Graph::NodeIndex kDegree = 10; + Graph graph; + GeneratePartialRandomGraph(num_tails, num_heads, kDegree, &graph); + Graphs::Build(&graph); + CHECK_EQ(graph.num_arcs(), num_tails * kDegree + num_tails + num_heads); + std::vector arc_capacity(graph.num_arcs(), 1); + std::unique_ptr > max_flow(new GenericMaxFlow( + &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); + SetUpNetworkData(arc_capacity, max_flow.get()); + FlowQuantity flow = f(max_flow.get()); + EXPECT_EQ(expected_flow1, flow); +} + +template +void ChangeCapacities(absl::Span arc_capacity, + FlowQuantity delta, GenericMaxFlow* max_flow) { + const Graph* graph = max_flow->graph(); + for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { + max_flow->SetArcCapacity(arc, + std::max(arc_capacity[arc] - delta, int64_t{0})); + } +} + +template +void PartialRandomFlow(typename MaxFlowSolver::Solver f, + typename Graph::NodeIndex num_tails, + typename Graph::NodeIndex num_heads, + FlowQuantity expected_flow1, + FlowQuantity expected_flow2) { + const typename Graph::NodeIndex kDegree = 10; + const FlowQuantity kCapacityRange = 10000; + const FlowQuantity kCapacityDelta = 1000; + Graph graph; + GeneratePartialRandomGraph(num_tails, num_heads, kDegree, &graph); + std::vector arc_capacity(graph.num_arcs()); + GenerateRandomArcValuations(graph, kCapacityRange, &arc_capacity); + + std::vector permutation; + Graphs::Build(&graph, &permutation); + util::Permute(permutation, &arc_capacity); + + std::unique_ptr > max_flow(new GenericMaxFlow( + &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); + SetUpNetworkData(arc_capacity, max_flow.get()); + FlowQuantity flow = f(max_flow.get()); + EXPECT_EQ(expected_flow1, flow); + ChangeCapacities(arc_capacity, kCapacityDelta, max_flow.get()); + flow = f(max_flow.get()); + EXPECT_EQ(expected_flow2, flow); + ChangeCapacities(arc_capacity, 0, max_flow.get()); + flow = f(max_flow.get()); + EXPECT_EQ(expected_flow1, flow); +} + +template +void FullRandomFlow(typename MaxFlowSolver::Solver f, + typename Graph::NodeIndex num_tails, + typename Graph::NodeIndex num_heads, + FlowQuantity expected_flow1, FlowQuantity expected_flow2) { + const FlowQuantity kCapacityRange = 10000; + const FlowQuantity kCapacityDelta = 1000; + Graph graph; + GenerateCompleteGraph(num_tails, num_heads, &graph); + std::vector arc_capacity(graph.num_arcs()); + GenerateRandomArcValuations(graph, kCapacityRange, &arc_capacity); + + std::vector permutation; + Graphs::Build(&graph, &permutation); + util::Permute(permutation, &arc_capacity); + + std::unique_ptr > max_flow(new GenericMaxFlow( + &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); + SetUpNetworkData(arc_capacity, max_flow.get()); + FlowQuantity flow = f(max_flow.get()); + EXPECT_EQ(expected_flow1, flow); + ChangeCapacities(arc_capacity, kCapacityDelta, max_flow.get()); + flow = f(max_flow.get()); + EXPECT_EQ(expected_flow2, flow); + ChangeCapacities(arc_capacity, 0, max_flow.get()); + flow = f(max_flow.get()); + EXPECT_EQ(expected_flow1, flow); +} + +#define LP_AND_FLOW_TEST(test_name, size, expected_flow1, expected_flow2) \ + LP_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ + FLOW_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ + FLOW_ONLY_TEST_SG(test_name, size, expected_flow1, expected_flow2) + +#define LP_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ + TEST(LPMaxFlowTest, test_name##size) { \ + test_name(SolveMaxFlowWithLP, size, size, \ + expected_flow1, expected_flow2); \ + } + +#define FLOW_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ + TEST(MaxFlowTest, test_name##size) { \ + test_name(SolveMaxFlow, size, size, expected_flow1, \ + expected_flow2); \ + } + +#define FLOW_ONLY_TEST_SG(test_name, size, expected_flow1, expected_flow2) \ + TEST(MaxFlowTestStaticGraph, test_name##size) { \ + test_name >(SolveMaxFlow, size, size, \ + expected_flow1, expected_flow2); \ + } + +LP_AND_FLOW_TEST(FullRandomAssignment, 300, 300, 300); +LP_AND_FLOW_TEST(PartialRandomAssignment, 100, 100, 100); +LP_AND_FLOW_TEST(PartialRandomAssignment, 1000, 1000, 1000); +LP_AND_FLOW_TEST(PartialRandomFlow, 400, 1898664, 1515203); +LP_AND_FLOW_TEST(FullRandomFlow, 100, 482391, 386587); + +// LARGE must be defined from the build command line to test larger instances. +#ifdef LARGE +LP_AND_FLOW_TEST(PartialRandomAssignment, 10000, 10000, 10000); +#endif + +template +static void BM_FullRandomAssignment(benchmark::State& state) { + const int kSize = 3000; + for (auto _ : state) { + FullRandomAssignment(SolveMaxFlow, kSize, kSize, kSize, kSize); + } + state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); +} + +template +static void BM_PartialRandomAssignment(benchmark::State& state) { + const int kSize = 10100; + for (auto _ : state) { + PartialRandomAssignment(SolveMaxFlow, kSize, kSize, kSize, kSize); + } + state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); +} + +template +static void BM_PartialRandomFlow(benchmark::State& state) { + const int kSize = 800; + for (auto _ : state) { + PartialRandomFlow(SolveMaxFlow, kSize, kSize, 3884850, 3112123); + } + state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); +} + +template +static void BM_FullRandomFlow(benchmark::State& state) { + const int kSize = 800; + for (auto _ : state) { + FullRandomFlow(SolveMaxFlow, kSize, kSize, 4000549, 3239512); + } + state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); +} + +BENCHMARK_TEMPLATE(BM_FullRandomAssignment, StarGraph); +BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcListGraph<>); +BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcStaticGraph<>); +BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcMixedGraph<>); +BENCHMARK_TEMPLATE(BM_PartialRandomFlow, StarGraph); +BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcListGraph<>); +BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcStaticGraph<>); +BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcMixedGraph<>); + +// One iteration of each of the following tests is slow. +BENCHMARK_TEMPLATE(BM_FullRandomFlow, StarGraph); +BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcListGraph<>); +BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcStaticGraph<>); +BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcMixedGraph<>); +BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, StarGraph); +BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcListGraph<>); +BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcStaticGraph<>); +BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcMixedGraph<>); + +#undef LP_AND_FLOW_TEST +#undef LP_ONLY_TEST +#undef FLOW_ONLY_TEST +#undef FLOW_ONLY_TEST_SG + +// ---------------------------------------------------------- +// PriorityQueueWithRestrictedPush tests. +// ---------------------------------------------------------- + +TEST(PriorityQueueWithRestrictedPushTest, BasicBehavior) { + PriorityQueueWithRestrictedPush queue; + EXPECT_TRUE(queue.IsEmpty()); + queue.Push("A", 1); + queue.Push("B", 0); + queue.Push("C", 2); + queue.Push("D", 10); + queue.Push("E", 9); + EXPECT_EQ("D", queue.Pop()); + EXPECT_EQ("E", queue.Pop()); + EXPECT_EQ("C", queue.Pop()); + EXPECT_EQ("A", queue.Pop()); + EXPECT_EQ("B", queue.Pop()); + EXPECT_TRUE(queue.IsEmpty()); + queue.Push("A", 1); + queue.Push("B", 0); + EXPECT_FALSE(queue.IsEmpty()); + queue.Clear(); + EXPECT_TRUE(queue.IsEmpty()); +} + +TEST(PriorityQueueWithRestrictedPushTest, BasicBehaviorWithMixedPushPop) { + PriorityQueueWithRestrictedPush queue; + EXPECT_TRUE(queue.IsEmpty()); + queue.Push("A", 1); + queue.Push("B", 0); + queue.Push("C", 2); + EXPECT_EQ("C", queue.Pop()); + EXPECT_EQ("A", queue.Pop()); + queue.Push("D", 1); + queue.Push("E", 0); + EXPECT_EQ("D", queue.Pop()); + EXPECT_EQ("E", queue.Pop()); + EXPECT_EQ("B", queue.Pop()); + EXPECT_TRUE(queue.IsEmpty()); + queue.Push("E", 1); + EXPECT_FALSE(queue.IsEmpty()); + EXPECT_EQ("E", queue.Pop()); + EXPECT_TRUE(queue.IsEmpty()); +} + +TEST(PriorityQueueWithRestrictedPushTest, RandomPushPop) { + struct ElementWithPriority { + ElementWithPriority(int _element, int _priority) + : element(_element), priority(_priority) {} + bool operator<(const ElementWithPriority& other) const { + return priority < other.priority; + } + int element; + int priority; + }; + std::vector pairs; + std::mt19937 randomizer(1); + const int kNumElement = 10000; + const int kMaxPriority = 10000; // We want duplicates and gaps. + for (int i = 0; i < kNumElement; ++i) { + pairs.push_back( + ElementWithPriority(i, absl::Uniform(randomizer, 0, kMaxPriority))); + } + std::sort(pairs.begin(), pairs.end()); + + // Randomly add +1 and push to the queue. + PriorityQueueWithRestrictedPush queue; + for (int i = 0; i < pairs.size(); ++i) { + pairs[i].priority += absl::Bernoulli(randomizer, 1.0 / 2) ? 1 : 0; + queue.Push(pairs[i].element, pairs[i].priority); + } + + // Stable sort the pairs for checking (the queue order is stable). + std::stable_sort(pairs.begin(), pairs.end()); + + // Random Push() and Pop() with more Pop(). + int current = pairs.size(); + while (!queue.IsEmpty()) { + EXPECT_GT(current, 0); + if (absl::Bernoulli(randomizer, 1.0 / 4) && current < pairs.size()) { + queue.Push(pairs[current].element, pairs[current].priority); + ++current; + } else { + --current; + EXPECT_EQ(pairs[current].element, queue.Pop()); + } + } +} + +TEST(PriorityQueueWithRestrictedPushDeathTest, DCHECK) { + // Don't run this test in opt mode. + if (!DEBUG_MODE) GTEST_SKIP(); + + PriorityQueueWithRestrictedPush queue; + EXPECT_TRUE(queue.IsEmpty()); + ASSERT_DEATH(queue.Pop(), ""); + queue.Push("A", 10); + queue.Push("B", 9); + ASSERT_DEATH(queue.Push("C", 4), ""); + ASSERT_DEATH(queue.Push("C", 8), ""); +} + +} // namespace +} // namespace operations_research diff --git a/ortools/graph/max_flow.cc b/ortools/graph/max_flow.cc index cbe1a725b4..630129d60e 100644 --- a/ortools/graph/max_flow.cc +++ b/ortools/graph/max_flow.cc @@ -14,25 +14,18 @@ #include "ortools/graph/max_flow.h" #include -#include #include -#include #include -#include "absl/log/check.h" -#include "absl/memory/memory.h" -#include "absl/strings/str_format.h" -#include "absl/strings/string_view.h" -#include "ortools/graph/ebert_graph.h" +#include "ortools/graph/generic_max_flow.h" #include "ortools/graph/graph.h" -#include "ortools/graph/graphs.h" namespace operations_research { SimpleMaxFlow::SimpleMaxFlow() : num_nodes_(0) {} -ArcIndex SimpleMaxFlow::AddArcWithCapacity(NodeIndex tail, NodeIndex head, - FlowQuantity capacity) { +SimpleMaxFlow::ArcIndex SimpleMaxFlow::AddArcWithCapacity( + NodeIndex tail, NodeIndex head, FlowQuantity capacity) { const ArcIndex num_arcs = arc_tail_.size(); num_nodes_ = std::max(num_nodes_, tail + 1); num_nodes_ = std::max(num_nodes_, head + 1); @@ -42,15 +35,21 @@ ArcIndex SimpleMaxFlow::AddArcWithCapacity(NodeIndex tail, NodeIndex head, return num_arcs; } -NodeIndex SimpleMaxFlow::NumNodes() const { return num_nodes_; } +SimpleMaxFlow::NodeIndex SimpleMaxFlow::NumNodes() const { return num_nodes_; } -ArcIndex SimpleMaxFlow::NumArcs() const { return arc_tail_.size(); } +SimpleMaxFlow::ArcIndex SimpleMaxFlow::NumArcs() const { + return arc_tail_.size(); +} -NodeIndex SimpleMaxFlow::Tail(ArcIndex arc) const { return arc_tail_[arc]; } +SimpleMaxFlow::NodeIndex SimpleMaxFlow::Tail(ArcIndex arc) const { + return arc_tail_[arc]; +} -NodeIndex SimpleMaxFlow::Head(ArcIndex arc) const { return arc_head_[arc]; } +SimpleMaxFlow::NodeIndex SimpleMaxFlow::Head(ArcIndex arc) const { + return arc_head_[arc]; +} -FlowQuantity SimpleMaxFlow::Capacity(ArcIndex arc) const { +SimpleMaxFlow::FlowQuantity SimpleMaxFlow::Capacity(ArcIndex arc) const { return arc_capacity_[arc]; } @@ -109,9 +108,13 @@ SimpleMaxFlow::Status SimpleMaxFlow::Solve(NodeIndex source, NodeIndex sink) { return BAD_RESULT; } -FlowQuantity SimpleMaxFlow::OptimalFlow() const { return optimal_flow_; } +SimpleMaxFlow::FlowQuantity SimpleMaxFlow::OptimalFlow() const { + return optimal_flow_; +} -FlowQuantity SimpleMaxFlow::Flow(ArcIndex arc) const { return arc_flow_[arc]; } +SimpleMaxFlow::FlowQuantity SimpleMaxFlow::Flow(ArcIndex arc) const { + return arc_flow_[arc]; +} void SimpleMaxFlow::GetSourceSideMinCut(std::vector* result) { if (underlying_max_flow_ == nullptr) return; @@ -142,968 +145,4 @@ FlowModelProto SimpleMaxFlow::CreateFlowModelProto(NodeIndex source, return model; } -template -GenericMaxFlow::GenericMaxFlow(const Graph* graph, NodeIndex source, - NodeIndex sink) - : graph_(graph), - node_excess_(), - node_potential_(), - residual_arc_capacity_(), - first_admissible_arc_(), - active_nodes_(), - source_(source), - sink_(sink), - use_global_update_(true), - use_two_phase_algorithm_(true), - process_node_by_height_(true), - check_input_(true), - check_result_(true), - stats_("MaxFlow") { - SCOPED_TIME_STAT(&stats_); - DCHECK(graph->IsNodeValid(source)); - DCHECK(graph->IsNodeValid(sink)); - const NodeIndex max_num_nodes = Graphs::NodeReservation(*graph_); - if (max_num_nodes > 0) { - node_excess_.Reserve(0, max_num_nodes - 1); - node_excess_.SetAll(0); - node_potential_.Reserve(0, max_num_nodes - 1); - node_potential_.SetAll(0); - first_admissible_arc_.Reserve(0, max_num_nodes - 1); - first_admissible_arc_.SetAll(Graph::kNilArc); - bfs_queue_.reserve(max_num_nodes); - active_nodes_.reserve(max_num_nodes); - } - const ArcIndex max_num_arcs = Graphs::ArcReservation(*graph_); - if (max_num_arcs > 0) { - residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1); - residual_arc_capacity_.SetAll(0); - } -} - -template -bool GenericMaxFlow::CheckInputConsistency() const { - SCOPED_TIME_STAT(&stats_); - bool ok = true; - for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) { - if (residual_arc_capacity_[arc] < 0) { - ok = false; - } - } - return ok; -} - -template -void GenericMaxFlow::SetArcCapacity(ArcIndex arc, - FlowQuantity new_capacity) { - SCOPED_TIME_STAT(&stats_); - DCHECK_LE(0, new_capacity); - DCHECK(IsArcDirect(arc)); - const FlowQuantity free_capacity = residual_arc_capacity_[arc]; - const FlowQuantity capacity_delta = new_capacity - Capacity(arc); - if (capacity_delta == 0) { - return; // Nothing to do. - } - status_ = NOT_SOLVED; - if (free_capacity + capacity_delta >= 0) { - // The above condition is true if one of the two conditions is true: - // 1/ (capacity_delta > 0), meaning we are increasing the capacity - // 2/ (capacity_delta < 0 && free_capacity + capacity_delta >= 0) - // meaning we are reducing the capacity, but that the capacity - // reduction is not larger than the free capacity. - DCHECK((capacity_delta > 0) || - (capacity_delta < 0 && free_capacity + capacity_delta >= 0)); - residual_arc_capacity_.Set(arc, free_capacity + capacity_delta); - DCHECK_LE(0, residual_arc_capacity_[arc]); - } else { - // Note that this breaks the preflow invariants but it is currently not an - // issue since we restart from scratch on each Solve() and we set the status - // to NOT_SOLVED. - // - // TODO(user): The easiest is probably to allow negative node excess in - // other places than the source, but the current implementation does not - // deal with this. - SetCapacityAndClearFlow(arc, new_capacity); - } -} - -template -void GenericMaxFlow::SetArcFlow(ArcIndex arc, FlowQuantity new_flow) { - SCOPED_TIME_STAT(&stats_); - DCHECK(IsArcValid(arc)); - DCHECK_GE(new_flow, 0); - const FlowQuantity capacity = Capacity(arc); - DCHECK_GE(capacity, new_flow); - - // Note that this breaks the preflow invariants but it is currently not an - // issue since we restart from scratch on each Solve() and we set the status - // to NOT_SOLVED. - residual_arc_capacity_.Set(Opposite(arc), -new_flow); - residual_arc_capacity_.Set(arc, capacity - new_flow); - status_ = NOT_SOLVED; -} - -template -void GenericMaxFlow::GetSourceSideMinCut( - std::vector* result) { - ComputeReachableNodes(source_, result); -} - -template -void GenericMaxFlow::GetSinkSideMinCut(std::vector* result) { - ComputeReachableNodes(sink_, result); -} - -template -bool GenericMaxFlow::CheckResult() const { - SCOPED_TIME_STAT(&stats_); - bool ok = true; - if (node_excess_[source_] != -node_excess_[sink_]) { - LOG(DFATAL) << "-node_excess_[source_] = " << -node_excess_[source_] - << " != node_excess_[sink_] = " << node_excess_[sink_]; - ok = false; - } - for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) { - if (node != source_ && node != sink_) { - if (node_excess_[node] != 0) { - LOG(DFATAL) << "node_excess_[" << node << "] = " << node_excess_[node] - << " != 0"; - ok = false; - } - } - } - for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) { - const ArcIndex opposite = Opposite(arc); - const FlowQuantity direct_capacity = residual_arc_capacity_[arc]; - const FlowQuantity opposite_capacity = residual_arc_capacity_[opposite]; - if (direct_capacity < 0) { - LOG(DFATAL) << "residual_arc_capacity_[" << arc - << "] = " << direct_capacity << " < 0"; - ok = false; - } - if (opposite_capacity < 0) { - LOG(DFATAL) << "residual_arc_capacity_[" << opposite - << "] = " << opposite_capacity << " < 0"; - ok = false; - } - // The initial capacity of the direct arcs is non-negative. - if (direct_capacity + opposite_capacity < 0) { - LOG(DFATAL) << "initial capacity [" << arc - << "] = " << direct_capacity + opposite_capacity << " < 0"; - ok = false; - } - } - return ok; -} - -template -bool GenericMaxFlow::AugmentingPathExists() const { - SCOPED_TIME_STAT(&stats_); - - // We simply compute the reachability from the source in the residual graph. - const NodeIndex num_nodes = graph_->num_nodes(); - std::vector is_reached(num_nodes, false); - std::vector to_process; - - to_process.push_back(source_); - is_reached[source_] = true; - while (!to_process.empty()) { - const NodeIndex node = to_process.back(); - to_process.pop_back(); - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); - if (residual_arc_capacity_[arc] > 0) { - const NodeIndex head = graph_->Head(arc); - if (!is_reached[head]) { - is_reached[head] = true; - to_process.push_back(head); - } - } - } - } - return is_reached[sink_]; -} - -template -bool GenericMaxFlow::CheckRelabelPrecondition(NodeIndex node) const { - DCHECK(IsActive(node)); - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); - DCHECK(!IsAdmissible(arc)) << DebugString("CheckRelabelPrecondition:", arc); - } - return true; -} - -template -std::string GenericMaxFlow::DebugString(absl::string_view context, - ArcIndex arc) const { - const NodeIndex tail = Tail(arc); - const NodeIndex head = Head(arc); - return absl::StrFormat( - "%s Arc %d, from %d to %d, " - "Capacity = %d, Residual capacity = %d, " - "Flow = residual capacity for reverse arc = %d, " - "Height(tail) = %d, Height(head) = %d, " - "Excess(tail) = %d, Excess(head) = %d", - context, arc, tail, head, Capacity(arc), residual_arc_capacity_[arc], - Flow(arc), node_potential_[tail], node_potential_[head], - node_excess_[tail], node_excess_[head]); -} - -template -bool GenericMaxFlow::Solve() { - status_ = NOT_SOLVED; - if (check_input_ && !CheckInputConsistency()) { - status_ = BAD_INPUT; - return false; - } - InitializePreflow(); - - // Deal with the case when source_ or sink_ is not inside graph_. - // Since they are both specified independently of the graph, we do need to - // take care of this corner case. - const NodeIndex num_nodes = graph_->num_nodes(); - if (sink_ >= num_nodes || source_ >= num_nodes) { - // Behave like a normal graph where source_ and sink_ are disconnected. - // Note that the arc flow is set to 0 by InitializePreflow(). - status_ = OPTIMAL; - return true; - } - if (use_global_update_) { - RefineWithGlobalUpdate(); - } else { - Refine(); - } - if (check_result_) { - if (!CheckResult()) { - status_ = BAD_RESULT; - return false; - } - if (GetOptimalFlow() < kMaxFlowQuantity && AugmentingPathExists()) { - LOG(ERROR) << "The algorithm terminated, but the flow is not maximal!"; - status_ = BAD_RESULT; - return false; - } - } - DCHECK_EQ(node_excess_[sink_], -node_excess_[source_]); - status_ = OPTIMAL; - if (GetOptimalFlow() == kMaxFlowQuantity && AugmentingPathExists()) { - // In this case, we are sure that the flow is > kMaxFlowQuantity. - status_ = INT_OVERFLOW; - } - IF_STATS_ENABLED(VLOG(1) << stats_.StatString()); - return true; -} - -template -void GenericMaxFlow::InitializePreflow() { - SCOPED_TIME_STAT(&stats_); - // InitializePreflow() clears the whole flow that could have been computed - // by a previous Solve(). This is not optimal in terms of complexity. - // TODO(user): find a way to make the re-solving incremental (not an obvious - // task, and there has not been a lot of literature on the subject.) - node_excess_.SetAll(0); - const ArcIndex num_arcs = graph_->num_arcs(); - for (ArcIndex arc = 0; arc < num_arcs; ++arc) { - SetCapacityAndClearFlow(arc, Capacity(arc)); - } - - // All the initial heights are zero except for the source whose height is - // equal to the number of nodes and will never change during the algorithm. - node_potential_.SetAll(0); - node_potential_.Set(source_, graph_->num_nodes()); - - // Initially no arcs are admissible except maybe the one leaving the source, - // but we treat the source in a special way, see - // SaturateOutgoingArcsFromSource(). - const NodeIndex num_nodes = graph_->num_nodes(); - for (NodeIndex node = 0; node < num_nodes; ++node) { - first_admissible_arc_[node] = Graph::kNilArc; - } -} - -// Note(user): Calling this function will break the property on the node -// potentials because of the way we cancel flow on cycle. However, we only call -// that at the end of the algorithm, or just before a GlobalUpdate() that will -// restore the precondition on the node potentials. -template -void GenericMaxFlow::PushFlowExcessBackToSource() { - SCOPED_TIME_STAT(&stats_); - const NodeIndex num_nodes = graph_->num_nodes(); - - // We implement a variation of Tarjan's strongly connected component algorithm - // to detect cycles published in: Tarjan, R. E. (1972), "Depth-first search - // and linear graph algorithms", SIAM Journal on Computing. A description can - // also be found in wikipedia. - - // Stored nodes are settled nodes already stored in the - // reverse_topological_order (except the sink_ that we do not actually store). - std::vector stored(num_nodes, false); - stored[sink_] = true; - - // The visited nodes that are not yet stored are all the nodes from the - // source_ to the current node in the current dfs branch. - std::vector visited(num_nodes, false); - visited[sink_] = true; - - // Stack of arcs to explore in the dfs search. - // The current node is Head(arc_stack.back()). - std::vector arc_stack; - - // Increasing list of indices into the arc_stack that correspond to the list - // of arcs in the current dfs branch from the source_ to the current node. - std::vector index_branch; - - // Node in reverse_topological_order in the final dfs tree. - std::vector reverse_topological_order; - - // We start by pushing all the outgoing arcs from the source on the stack to - // avoid special conditions in the code. As a result, source_ will not be - // stored in reverse_topological_order, and this is what we want. - for (OutgoingArcIterator it(*graph_, source_); it.Ok(); it.Next()) { - const ArcIndex arc = it.Index(); - const FlowQuantity flow = Flow(arc); - if (flow > 0) { - arc_stack.push_back(arc); - } - } - visited[source_] = true; - - // Start the dfs on the subgraph formed by the direct arcs with positive flow. - while (!arc_stack.empty()) { - const NodeIndex node = Head(arc_stack.back()); - - // If the node is visited, it means we have explored all its arcs and we - // have just backtracked in the dfs. Store it if it is not already stored - // and process the next arc on the stack. - if (visited[node]) { - if (!stored[node]) { - stored[node] = true; - reverse_topological_order.push_back(node); - DCHECK(!index_branch.empty()); - index_branch.pop_back(); - } - arc_stack.pop_back(); - continue; - } - - // The node is a new unexplored node, add all its outgoing arcs with - // positive flow to the stack and go deeper in the dfs. - DCHECK(!stored[node]); - DCHECK(index_branch.empty() || - (arc_stack.size() - 1 > index_branch.back())); - visited[node] = true; - index_branch.push_back(arc_stack.size() - 1); - - for (OutgoingArcIterator it(*graph_, node); it.Ok(); it.Next()) { - const ArcIndex arc = it.Index(); - const FlowQuantity flow = Flow(arc); - const NodeIndex head = Head(arc); - if (flow > 0 && !stored[head]) { - if (!visited[head]) { - arc_stack.push_back(arc); - } else { - // There is a cycle. - // Find the first index to consider, - // arc_stack[index_branch[cycle_begin]] will be the first arc on the - // cycle. - int cycle_begin = index_branch.size(); - while (cycle_begin > 0 && - Head(arc_stack[index_branch[cycle_begin - 1]]) != head) { - --cycle_begin; - } - - // Compute the maximum flow that can be canceled on the cycle and the - // min index such that arc_stack[index_branch[i]] will be saturated. - FlowQuantity max_flow = flow; - int first_saturated_index = index_branch.size(); - for (int i = index_branch.size() - 1; i >= cycle_begin; --i) { - const ArcIndex arc_on_cycle = arc_stack[index_branch[i]]; - if (Flow(arc_on_cycle) <= max_flow) { - max_flow = Flow(arc_on_cycle); - first_saturated_index = i; - } - } - - // This is just here for a DCHECK() below. - const FlowQuantity excess = node_excess_[head]; - - // Cancel the flow on the cycle, and set visited[node] = false for - // the node that will be backtracked over. - PushFlow(-max_flow, arc); - for (int i = index_branch.size() - 1; i >= cycle_begin; --i) { - const ArcIndex arc_on_cycle = arc_stack[index_branch[i]]; - PushFlow(-max_flow, arc_on_cycle); - if (i >= first_saturated_index) { - DCHECK(visited[Head(arc_on_cycle)]); - visited[Head(arc_on_cycle)] = false; - } else { - DCHECK_GT(Flow(arc_on_cycle), 0); - } - } - - // This is a simple check that the flow was pushed properly. - DCHECK_EQ(excess, node_excess_[head]); - - // Backtrack the dfs just before index_branch[first_saturated_index]. - // If the current node is still active, there is nothing to do. - if (first_saturated_index < index_branch.size()) { - arc_stack.resize(index_branch[first_saturated_index]); - index_branch.resize(first_saturated_index); - - // We backtracked over the current node, so there is no need to - // continue looping over its arcs. - break; - } - } - } - } - } - DCHECK(arc_stack.empty()); - DCHECK(index_branch.empty()); - - // Return the flow to the sink. Note that the sink_ and the source_ are not - // stored in reverse_topological_order. - for (int i = 0; i < reverse_topological_order.size(); i++) { - const NodeIndex node = reverse_topological_order[i]; - if (node_excess_[node] == 0) continue; - for (IncomingArcIterator it(*graph_, node); it.Ok(); it.Next()) { - const ArcIndex opposite_arc = Opposite(it.Index()); - if (residual_arc_capacity_[opposite_arc] > 0) { - const FlowQuantity flow = - std::min(node_excess_[node], residual_arc_capacity_[opposite_arc]); - PushFlow(flow, opposite_arc); - if (node_excess_[node] == 0) break; - } - } - DCHECK_EQ(0, node_excess_[node]); - } - DCHECK_EQ(-node_excess_[source_], node_excess_[sink_]); -} - -template -void GenericMaxFlow::GlobalUpdate() { - SCOPED_TIME_STAT(&stats_); - bfs_queue_.clear(); - int queue_index = 0; - const NodeIndex num_nodes = graph_->num_nodes(); - node_in_bfs_queue_.assign(num_nodes, false); - node_in_bfs_queue_[sink_] = true; - node_in_bfs_queue_[source_] = true; - - // We do two BFS in the reverse residual graph, one from the sink and one from - // the source. Because all the arcs from the source are saturated (except in - // presence of integer overflow), the source cannot reach the sink in the - // residual graph. However, we still want to relabel all the nodes that cannot - // reach the sink but can reach the source (because if they have excess, we - // need to push it back to the source). - // - // Note that the second pass is not needed here if we use a two-pass algorithm - // to return the flow to the source after we found the min cut. - const int num_passes = use_two_phase_algorithm_ ? 1 : 2; - for (int pass = 0; pass < num_passes; ++pass) { - if (pass == 0) { - bfs_queue_.push_back(sink_); - } else { - bfs_queue_.push_back(source_); - } - - while (queue_index != bfs_queue_.size()) { - const NodeIndex node = bfs_queue_[queue_index]; - ++queue_index; - const NodeIndex candidate_distance = node_potential_[node] + 1; - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); - const NodeIndex head = Head(arc); - - // Skip the arc if the height of head was already set to the correct - // value (Remember we are doing reverse BFS). - if (node_in_bfs_queue_[head]) continue; - - // TODO(user): By using more memory we can speed this up quite a bit by - // avoiding to take the opposite arc here, too options: - // - if (residual_arc_capacity_[arc] != arc_capacity_[arc]) - // - if (opposite_arc_is_admissible_[arc]) // need updates. - // Experiment with the first option shows more than 10% gain on this - // function running time, which is the bottleneck on many instances. - const ArcIndex opposite_arc = Opposite(arc); - if (residual_arc_capacity_[opposite_arc] > 0) { - // Note(user): We used to have a DCHECK_GE(candidate_distance, - // node_potential_[head]); which is always true except in the case - // where we can push more than kMaxFlowQuantity out of the source. The - // problem comes from the fact that in this case, we call - // PushFlowExcessBackToSource() in the middle of the algorithm. The - // later call will break the properties of the node potential. Note - // however, that this function will recompute a good node potential - // for all the nodes and thus fix the issue. - - // If head is active, we can steal some or all of its excess. - // This brings a huge gain on some problems. - // Note(user): I haven't seen this anywhere in the literature. - // TODO(user): Investigate more and maybe write a publication :) - if (node_excess_[head] > 0) { - const FlowQuantity flow = std::min( - node_excess_[head], residual_arc_capacity_[opposite_arc]); - PushFlow(flow, opposite_arc); - - // If the arc became saturated, it is no longer in the residual - // graph, so we do not need to consider head at this time. - if (residual_arc_capacity_[opposite_arc] == 0) continue; - } - - // Note that there is no need to touch first_admissible_arc_[node] - // because of the relaxed Relabel() we use. - node_potential_[head] = candidate_distance; - node_in_bfs_queue_[head] = true; - bfs_queue_.push_back(head); - } - } - } - } - - // At the end of the search, some nodes may not be in the bfs_queue_. Such - // nodes cannot reach the sink_ or source_ in the residual graph, so there is - // no point trying to push flow toward them. We obtain this effect by setting - // their height to something unreachable. - // - // Note that this also prevents cycling due to our anti-overflow procedure. - // For instance, suppose there is an edge s -> n outgoing from the source. If - // node n has no other connection and some excess, we will push the flow back - // to the source, but if we don't update the height of n - // SaturateOutgoingArcsFromSource() will push the flow to n again. - // TODO(user): This is another argument for another anti-overflow algorithm. - for (NodeIndex node = 0; node < num_nodes; ++node) { - if (!node_in_bfs_queue_[node]) { - node_potential_[node] = 2 * num_nodes - 1; - } - } - - // Reset the active nodes. Doing it like this pushes the nodes in increasing - // order of height. Note that bfs_queue_[0] is the sink_ so we skip it. - DCHECK(IsEmptyActiveNodeContainer()); - for (int i = 1; i < bfs_queue_.size(); ++i) { - const NodeIndex node = bfs_queue_[i]; - if (node_excess_[node] > 0) { - DCHECK(IsActive(node)); - PushActiveNode(node); - } - } -} - -template -bool GenericMaxFlow::SaturateOutgoingArcsFromSource() { - SCOPED_TIME_STAT(&stats_); - const NodeIndex num_nodes = graph_->num_nodes(); - - // If sink_ or source_ already have kMaxFlowQuantity, then there is no - // point pushing more flow since it will cause an integer overflow. - if (node_excess_[sink_] == kMaxFlowQuantity) return false; - if (node_excess_[source_] == -kMaxFlowQuantity) return false; - - bool flow_pushed = false; - for (OutgoingArcIterator it(*graph_, source_); it.Ok(); it.Next()) { - const ArcIndex arc = it.Index(); - const FlowQuantity flow = residual_arc_capacity_[arc]; - - // This is a special IsAdmissible() condition for the source. - if (flow == 0 || node_potential_[Head(arc)] >= num_nodes) continue; - - // We are careful in case the sum of the flow out of the source is greater - // than kMaxFlowQuantity to avoid overflow. - const FlowQuantity current_flow_out_of_source = -node_excess_[source_]; - DCHECK_GE(flow, 0) << flow; - DCHECK_GE(current_flow_out_of_source, 0) << current_flow_out_of_source; - const FlowQuantity capped_flow = - kMaxFlowQuantity - current_flow_out_of_source; - if (capped_flow < flow) { - // We push as much flow as we can so the current flow on the network will - // be kMaxFlowQuantity. - - // Since at the beginning of the function, current_flow_out_of_source - // was different from kMaxFlowQuantity, we are sure to have pushed some - // flow before if capped_flow is 0. - if (capped_flow == 0) return true; - PushFlow(capped_flow, arc); - return true; - } - PushFlow(flow, arc); - flow_pushed = true; - } - DCHECK_LE(node_excess_[source_], 0); - return flow_pushed; -} - -template -void GenericMaxFlow::PushFlow(FlowQuantity flow, ArcIndex arc) { - SCOPED_TIME_STAT(&stats_); - // TODO(user): Do not allow a zero flow after fixing the UniformMaxFlow code. - DCHECK_GE(residual_arc_capacity_[Opposite(arc)] + flow, 0); - DCHECK_GE(residual_arc_capacity_[arc] - flow, 0); - - // node_excess_ should be always greater than or equal to 0 except for the - // source where it should always be smaller than or equal to 0. Note however - // that we cannot check this because when we cancel the flow on a cycle in - // PushFlowExcessBackToSource(), we may break this invariant during the - // operation even if it is still valid at the end. - - // Update the residual capacity of the arc and its opposite arc. - residual_arc_capacity_[arc] -= flow; - residual_arc_capacity_[Opposite(arc)] += flow; - - // Update the excesses at the tail and head of the arc. - node_excess_[Tail(arc)] -= flow; - node_excess_[Head(arc)] += flow; -} - -template -void GenericMaxFlow::InitializeActiveNodeContainer() { - SCOPED_TIME_STAT(&stats_); - DCHECK(IsEmptyActiveNodeContainer()); - const NodeIndex num_nodes = graph_->num_nodes(); - for (NodeIndex node = 0; node < num_nodes; ++node) { - if (IsActive(node)) { - if (use_two_phase_algorithm_ && node_potential_[node] >= num_nodes) { - continue; - } - PushActiveNode(node); - } - } -} - -template -void GenericMaxFlow::Refine() { - SCOPED_TIME_STAT(&stats_); - // Usually SaturateOutgoingArcsFromSource() will saturate all the arcs from - // the source in one go, and we will loop just once. But in case we can push - // more than kMaxFlowQuantity out of the source the loop is as follow: - // - Push up to kMaxFlowQuantity out of the source on the admissible outgoing - // arcs. Stop if no flow was pushed. - // - Compute the current max-flow. This will push some flow back to the - // source and render more outgoing arcs from the source not admissible. - // - // TODO(user): This may not be the most efficient algorithm if we need to loop - // many times. An alternative may be to handle the source like the other nodes - // in the algorithm, initially putting an excess of kMaxFlowQuantity on it, - // and making the source active like any other node with positive excess. To - // investigate. - // - // TODO(user): The code below is buggy when more than kMaxFlowQuantity can be - // pushed out of the source (i.e. when we loop more than once in the while()). - // This is not critical, since this code is not used in the default algorithm - // computation. The issue is twofold: - // - InitializeActiveNodeContainer() doesn't push the nodes in - // the correct order. - // - PushFlowExcessBackToSource() may break the node potential properties, and - // we will need a call to GlobalUpdate() to fix that. - while (SaturateOutgoingArcsFromSource()) { - DCHECK(IsEmptyActiveNodeContainer()); - InitializeActiveNodeContainer(); - while (!IsEmptyActiveNodeContainer()) { - const NodeIndex node = GetAndRemoveFirstActiveNode(); - if (node == source_ || node == sink_) continue; - Discharge(node); - } - if (use_two_phase_algorithm_) { - PushFlowExcessBackToSource(); - } - } -} - -template -void GenericMaxFlow::RefineWithGlobalUpdate() { - SCOPED_TIME_STAT(&stats_); - - // TODO(user): This should be graph_->num_nodes(), but ebert graph does not - // have a correct size if the highest index nodes have no arcs. - const NodeIndex num_nodes = Graphs::NodeReservation(*graph_); - std::vector skip_active_node; - - while (SaturateOutgoingArcsFromSource()) { - int num_skipped; - do { - num_skipped = 0; - skip_active_node.assign(num_nodes, 0); - skip_active_node[sink_] = 2; - skip_active_node[source_] = 2; - GlobalUpdate(); - while (!IsEmptyActiveNodeContainer()) { - const NodeIndex node = GetAndRemoveFirstActiveNode(); - if (skip_active_node[node] > 1) { - if (node != sink_ && node != source_) ++num_skipped; - continue; - } - const NodeIndex old_height = node_potential_[node]; - Discharge(node); - - // The idea behind this is that if a node height augments by more than - // one, then it is likely to push flow back the way it came. This can - // lead to very costly loops. A bad case is: source -> n1 -> n2 and n2 - // just recently isolated from the sink. Then n2 will push flow back to - // n1, and n1 to n2 and so on. The height of each node will increase by - // steps of two until the height of the source is reached, which can - // take a long time. If the chain is longer, the situation is even - // worse. The behavior of this heuristic is related to the Gap - // heuristic. - // - // Note that the global update will fix all such cases efficiently. So - // the idea is to discharge the active node as much as possible, and - // then do a global update. - // - // We skip a node when this condition was true 2 times to avoid doing a - // global update too frequently. - if (node_potential_[node] > old_height + 1) { - ++skip_active_node[node]; - } - } - } while (num_skipped > 0); - if (use_two_phase_algorithm_) { - PushFlowExcessBackToSource(); - } - } -} - -template -void GenericMaxFlow::Discharge(NodeIndex node) { - SCOPED_TIME_STAT(&stats_); - const NodeIndex num_nodes = graph_->num_nodes(); - while (true) { - DCHECK(IsActive(node)); - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node, - first_admissible_arc_[node]); - it.Ok(); it.Next()) { - const ArcIndex arc = it.Index(); - if (IsAdmissible(arc)) { - DCHECK(IsActive(node)); - const NodeIndex head = Head(arc); - if (node_excess_[head] == 0) { - // The push below will make the node active for sure. Note that we may - // push the sink_, but that is handled properly in Refine(). - PushActiveNode(head); - } - const FlowQuantity delta = - std::min(node_excess_[node], residual_arc_capacity_[arc]); - PushFlow(delta, arc); - if (node_excess_[node] == 0) { - first_admissible_arc_[node] = arc; // arc may still be admissible. - return; - } - } - } - Relabel(node); - if (use_two_phase_algorithm_ && node_potential_[node] >= num_nodes) break; - } -} - -template -void GenericMaxFlow::Relabel(NodeIndex node) { - SCOPED_TIME_STAT(&stats_); - // Because we use a relaxed version, this is no longer true if the - // first_admissible_arc_[node] was not actually the first arc! - // DCHECK(CheckRelabelPrecondition(node)); - NodeHeight min_height = std::numeric_limits::max(); - ArcIndex first_admissible_arc = Graph::kNilArc; - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); - if (residual_arc_capacity_[arc] > 0) { - // Update min_height only for arcs with available capacity. - NodeHeight head_height = node_potential_[Head(arc)]; - if (head_height < min_height) { - min_height = head_height; - first_admissible_arc = arc; - - // We found an admissible arc at the current height, just stop there. - // This is the true first_admissible_arc_[node]. - if (min_height + 1 == node_potential_[node]) break; - } - } - } - DCHECK_NE(first_admissible_arc, Graph::kNilArc); - node_potential_[node] = min_height + 1; - - // Note that after a Relabel(), the loop will continue in Discharge(), and - // we are sure that all the arcs before first_admissible_arc are not - // admissible since their height is > min_height. - first_admissible_arc_[node] = first_admissible_arc; -} - -template -typename Graph::ArcIndex GenericMaxFlow::Opposite(ArcIndex arc) const { - return Graphs::OppositeArc(*graph_, arc); -} - -template -bool GenericMaxFlow::IsArcDirect(ArcIndex arc) const { - return IsArcValid(arc) && arc >= 0; -} - -template -bool GenericMaxFlow::IsArcValid(ArcIndex arc) const { - return Graphs::IsArcValid(*graph_, arc); -} - -template -const FlowQuantity GenericMaxFlow::kMaxFlowQuantity = - std::numeric_limits::max(); - -template -template -void GenericMaxFlow::ComputeReachableNodes( - NodeIndex start, std::vector* result) { - // If start is not a valid node index, it can reach only itself. - // Note(user): This is needed because source and sink are given independently - // of the graph and sometimes before it is even constructed. - const NodeIndex num_nodes = graph_->num_nodes(); - if (start >= num_nodes) { - result->clear(); - result->push_back(start); - return; - } - bfs_queue_.clear(); - node_in_bfs_queue_.assign(num_nodes, false); - - int queue_index = 0; - bfs_queue_.push_back(start); - node_in_bfs_queue_[start] = true; - while (queue_index != bfs_queue_.size()) { - const NodeIndex node = bfs_queue_[queue_index]; - ++queue_index; - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); - const NodeIndex head = Head(arc); - if (node_in_bfs_queue_[head]) continue; - if (residual_arc_capacity_[reverse ? Opposite(arc) : arc] == 0) continue; - node_in_bfs_queue_[head] = true; - bfs_queue_.push_back(head); - } - } - *result = bfs_queue_; -} - -template -FlowModelProto GenericMaxFlow::CreateFlowModel() { - FlowModelProto model; - model.set_problem_type(FlowModelProto::MAX_FLOW); - for (int n = 0; n < graph_->num_nodes(); ++n) { - FlowNodeProto* node = model.add_nodes(); - node->set_id(n); - if (n == source_) node->set_supply(1); - if (n == sink_) node->set_supply(-1); - } - for (int a = 0; a < graph_->num_arcs(); ++a) { - FlowArcProto* arc = model.add_arcs(); - arc->set_tail(graph_->Tail(a)); - arc->set_head(graph_->Head(a)); - arc->set_capacity(Capacity(a)); - } - return model; -} - -// Explicit instantiations that can be used by a client. -// -// TODO(user): moves this code out of a .cc file and include it at the end of -// the header so it can work with any graph implementation ? -template <> -const FlowQuantity GenericMaxFlow::kMaxFlowQuantity = - std::numeric_limits::max(); -template <> -const FlowQuantity - GenericMaxFlow<::util::ReverseArcListGraph<>>::kMaxFlowQuantity = - std::numeric_limits::max(); -template <> -const FlowQuantity - GenericMaxFlow<::util::ReverseArcStaticGraph<>>::kMaxFlowQuantity = - std::numeric_limits::max(); -template <> -const FlowQuantity - GenericMaxFlow<::util::ReverseArcMixedGraph<>>::kMaxFlowQuantity = - std::numeric_limits::max(); - -template class GenericMaxFlow; -template class GenericMaxFlow<::util::ReverseArcListGraph<>>; -template class GenericMaxFlow<::util::ReverseArcStaticGraph<>>; -template class GenericMaxFlow<::util::ReverseArcMixedGraph<>>; - -std::vector BipartiteMinimumVertexCover( - const std::vector>& left_to_right_arcs, int num_right) { - // This algorithm first uses the maximum flow to find a maximum matching. Then - // it uses the same method outlined in the proof of Konig's theorem to - // transform the maximum matching into a minimum vertex cover. - // - // More concretely, it uses a DFS starting with unmatched nodes and - // alternating matched/unmatched edges to find a minimum vertex cover. - SimpleMaxFlow max_flow; - const int num_left = left_to_right_arcs.size(); - std::vector arcs; - for (int i = 0; i < num_left; ++i) { - for (const int right_node : left_to_right_arcs[i]) { - DCHECK_GE(right_node, num_left); - DCHECK_LT(right_node, num_right + num_left); - arcs.push_back(max_flow.AddArcWithCapacity(i, right_node, 1)); - } - } - std::vector> adj_list = left_to_right_arcs; - adj_list.resize(num_left + num_right); - for (int i = 0; i < num_left; ++i) { - for (const int right_node : left_to_right_arcs[i]) { - adj_list[right_node].push_back(i); - } - } - const int sink = num_left + num_right; - const int source = num_left + num_right + 1; - for (int i = 0; i < num_left; ++i) { - max_flow.AddArcWithCapacity(source, i, 1); - } - for (int i = 0; i < num_right; ++i) { - max_flow.AddArcWithCapacity(i + num_left, sink, 1); - } - CHECK(max_flow.Solve(source, sink) == SimpleMaxFlow::OPTIMAL); - std::vector maximum_matching(num_left + num_right, -1); - for (const ArcIndex arc : arcs) { - if (max_flow.Flow(arc) > 0) { - maximum_matching[max_flow.Tail(arc)] = max_flow.Head(arc); - maximum_matching[max_flow.Head(arc)] = max_flow.Tail(arc); - } - } - // We do a DFS starting with unmatched nodes and alternating matched/unmatched - // edges. - std::vector in_alternating_path(num_left + num_right, false); - std::vector to_visit; - for (int i = 0; i < num_left; ++i) { - if (maximum_matching[i] == -1) { - to_visit.push_back(i); - } - } - while (!to_visit.empty()) { - const int current = to_visit.back(); - to_visit.pop_back(); - if (in_alternating_path[current]) { - continue; - } - in_alternating_path[current] = true; - for (const int j : adj_list[current]) { - if (current < num_left && maximum_matching[current] != j) { - to_visit.push_back(j); - } else if (current >= num_left && maximum_matching[current] == j) { - to_visit.push_back(j); - } - } - } - std::vector minimum_vertex_cover(num_left + num_right, false); - for (int i = 0; i < num_left; ++i) { - if (!in_alternating_path[i]) { - minimum_vertex_cover[i] = true; - } - } - for (int i = num_left; i < num_left + num_right; ++i) { - if (in_alternating_path[i]) { - minimum_vertex_cover[i] = true; - } - } - return minimum_vertex_cover; -} - } // namespace operations_research diff --git a/ortools/graph/max_flow.h b/ortools/graph/max_flow.h index b9cfa38d69..13c158f1b6 100644 --- a/ortools/graph/max_flow.h +++ b/ortools/graph/max_flow.h @@ -11,137 +11,19 @@ // See the License for the specific language governing permissions and // limitations under the License. -// An implementation of a push-relabel algorithm for the max flow problem. -// -// In the following, we consider a graph G = (V,E,s,t) where V denotes the set -// of nodes (vertices) in the graph, E denotes the set of arcs (edges). s and t -// denote distinguished nodes in G called source and target. n = |V| denotes the -// number of nodes in the graph, and m = |E| denotes the number of arcs in the -// graph. -// -// Each arc (v,w) is associated a capacity c(v,w). -// -// A flow is a function from E to R such that: -// -// a) f(v,w) <= c(v,w) for all (v,w) in E (capacity constraint.) -// -// b) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint.) -// -// c) sum on v f(v,w) = 0 (flow conservation.) -// -// The goal of this algorithm is to find the maximum flow from s to t, i.e. -// for example to maximize sum v f(s,v). -// -// The starting reference for this class of algorithms is: -// A.V. Goldberg and R.E. Tarjan. A new approach to the maximum flow problem. -// ACM Symposium on Theory of Computing, pp. 136-146. -// http://portal.acm.org/citation.cfm?id=12144. -// -// The basic idea of the algorithm is to handle preflows instead of flows, -// and to refine preflows until a maximum flow is obtained. -// A preflow is like a flow, except that the inflow can be larger than the -// outflow. If it is the case at a given node v, it is said that there is an -// excess at node v, and inflow = outflow + excess. -// -// More formally, a preflow is a function f such that: -// -// 1) f(v,w) <= c(v,w) for all (v,w) in E (capacity constraint). c(v,w) is a -// value representing the maximum capacity for arc (v,w). -// -// 2) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint) -// -// 3) excess(v) = sum on u f(u,v) >= 0 is the excess at node v, the -// algebraic sum of all the incoming preflows at this node. -// -// Each node has an associated "height", in addition to its excess. The -// height of the source is defined to be equal to n, and cannot change. The -// height of the target is defined to be zero, and cannot change either. The -// height of all the other nodes is initialized at zero and is updated during -// the algorithm (see below). For those who want to know the details, the height -// of a node, corresponds to a reduced cost, and this enables one to prove that -// the algorithm actually computes the max flow. Note that the height of a node -// can be initialized to the distance to the target node in terms of number of -// nodes. This has not been tried in this implementation. -// -// A node v is said to be *active* if excess(v) > 0. -// -// In this case the following operations can be applied to it: -// -// - if there are *admissible* incident arcs, i.e. arcs which are not saturated, -// and whose head's height is lower than the height of the active node -// considered, a PushFlow operation can be applied. It consists in sending as -// much flow as both the excess at the node and the capacity of the arc -// permit. -// - if there are no admissible arcs, the active node considered is relabeled, -// i.e. its height is increased to 1 + the minimum height of its neighboring -// nodes on admissible arcs. -// This is implemented in Discharge, which itself calls PushFlow and Relabel. -// -// Before running Discharge, it is necessary to initialize the algorithm with a -// preflow. This is done in InitializePreflow, which saturates all the arcs -// leaving the source node, and sets the excess at the heads of those arcs -// accordingly. -// -// The algorithm terminates when there are no remaining active nodes, i.e. all -// the excesses at all nodes are equal to zero. In this case, a maximum flow is -// obtained. -// -// The complexity of this algorithm depends amongst other things on the choice -// of the next active node. It has been shown, for example in: -// L. Tuncel, "On the Complexity of Preflow-Push Algorithms for Maximum-Flow -// Problems", Algorithmica 11(4): 353-359 (1994). -// and -// J. Cheriyan and K. Mehlhorn, "An analysis of the highest-level selection rule -// in the preflow-push max-flow algorithm", Information processing letters, -// 69(5):239-242 (1999). -// http://www.math.uwaterloo.ca/~jcheriya/PS_files/me3.0.ps -// -// ...that choosing the active node with the highest level yields a -// complexity of O(n^2 * sqrt(m)). -// -// TODO(user): implement the above active node choice rule. -// -// This has been validated experimentally in: -// R.K. Ahuja, M. Kodialam, A.K. Mishra, and J.B. Orlin, "Computational -// Investigations of Maximum Flow Algorithms", EJOR 97:509-542(1997). -// http://jorlin.scripts.mit.edu/docs/publications/58-comput%20investigations%20of.pdf. -// -// -// TODO(user): an alternative would be to evaluate: -// A.V. Goldberg, "The Partial Augment-Relabel Algorithm for the Maximum Flow -// Problem.” In Proceedings of Algorithms ESA, LNCS 5193:466-477, Springer 2008. -// http://www.springerlink.com/index/5535k2j1mt646338.pdf -// -// An interesting general reference on network flows is: -// R. K. Ahuja, T. L. Magnanti, J. B. Orlin, "Network Flows: Theory, Algorithms, -// and Applications," Prentice Hall, 1993, ISBN: 978-0136175490, -// http://www.amazon.com/dp/013617549X -// -// Keywords: Push-relabel, max-flow, network, graph, Goldberg, Tarjan, Dinic, -// Dinitz. - #ifndef OR_TOOLS_GRAPH_MAX_FLOW_H_ #define OR_TOOLS_GRAPH_MAX_FLOW_H_ +#include #include -#include -#include #include -#include "absl/strings/string_view.h" -#include "ortools/base/logging.h" -#include "ortools/graph/ebert_graph.h" #include "ortools/graph/flow_problem.pb.h" +#include "ortools/graph/generic_max_flow.h" #include "ortools/graph/graph.h" -#include "ortools/util/stats.h" -#include "ortools/util/zvector.h" namespace operations_research { -// Forward declaration. -template -class GenericMaxFlow; - // A simple and efficient max-cost flow interface. This is as fast as // GenericMaxFlow, which is the fastest, but uses // more memory in order to hide the somewhat involved construction of the @@ -150,6 +32,10 @@ class GenericMaxFlow; // TODO(user): If the need arises, extend this interface to support warm start. class SimpleMaxFlow { public: + typedef int32_t NodeIndex; + typedef int32_t ArcIndex; + typedef int64_t FlowQuantity; + // The constructor takes no size. // New node indices will be created lazily by AddArcWithCapacity(). SimpleMaxFlow(); @@ -246,514 +132,9 @@ class SimpleMaxFlow { // instance that uses it. typedef ::util::ReverseArcStaticGraph Graph; std::unique_ptr underlying_graph_; - std::unique_ptr> underlying_max_flow_; + std::unique_ptr > underlying_max_flow_; }; -// Specific but efficient priority queue implementation. The priority type must -// be an integer. The queue allows to retrieve the element with highest priority -// but only allows pushes with a priority greater or equal to the highest -// priority in the queue minus one. All operations are in O(1) and the memory is -// in O(num elements in the queue). Elements with the same priority are -// retrieved with LIFO order. -// -// Note(user): As far as I know, this is an original idea and is the only code -// that use this in the Maximum Flow context. Papers usually refer to an -// height-indexed array of simple linked lists of active node with the same -// height. Even worse, sometimes they use double-linked list to allow arbitrary -// height update in order to detect missing height (used for the Gap heuristic). -// But this can actually be implemented a lot more efficiently by just -// maintaining the height distribution of all the node in the graph. -template -class PriorityQueueWithRestrictedPush { - public: - PriorityQueueWithRestrictedPush() : even_queue_(), odd_queue_() {} - -#ifndef SWIG - // This type is neither copyable nor movable. - PriorityQueueWithRestrictedPush(const PriorityQueueWithRestrictedPush&) = - delete; - PriorityQueueWithRestrictedPush& operator=( - const PriorityQueueWithRestrictedPush&) = delete; -#endif - - // Is the queue empty? - bool IsEmpty() const; - - // Clears the queue. - void Clear(); - - // Push a new element in the queue. Its priority must be greater or equal to - // the highest priority present in the queue, minus one. This condition is - // DCHECKed, and violating it yields erroneous queue behavior in NDEBUG mode. - void Push(Element element, IntegerPriority priority); - - // Returns the element with highest priority and remove it from the queue. - // IsEmpty() must be false, this condition is DCHECKed. - Element Pop(); - - private: - // Helper function to get the last element of a vector and pop it. - Element PopBack(std::vector>* queue); - - // This is the heart of the algorithm. basically we split the elements by - // parity of their priority and the precondition on the Push() ensures that - // both vectors are always sorted by increasing priority. - std::vector> even_queue_; - std::vector> odd_queue_; -}; - -// We want an enum for the Status of a max flow run, and we want this -// enum to be scoped under GenericMaxFlow<>. Unfortunately, swig -// doesn't handle templated enums very well, so we need a base, -// untemplated class to hold it. -class MaxFlowStatusClass { - public: - enum Status { - NOT_SOLVED, // The problem was not solved, or its data were edited. - OPTIMAL, // Solve() was called and found an optimal solution. - INT_OVERFLOW, // There is a feasible flow > max possible flow. - BAD_INPUT, // The input is inconsistent. - BAD_RESULT // There was an error. - }; -}; - -// Generic MaxFlow (there is a default MaxFlow specialization defined below) -// that works with StarGraph and all the reverse arc graphs from graph.h, see -// the end of max_flow.cc for the exact types this class is compiled for. -template -class GenericMaxFlow : public MaxFlowStatusClass { - public: - typedef typename Graph::NodeIndex NodeIndex; - typedef typename Graph::ArcIndex ArcIndex; - typedef typename Graph::OutgoingArcIterator OutgoingArcIterator; - typedef typename Graph::OutgoingOrOppositeIncomingArcIterator - OutgoingOrOppositeIncomingArcIterator; - typedef typename Graph::IncomingArcIterator IncomingArcIterator; - typedef ZVector ArcIndexArray; - - // The height of a node never excess 2 times the number of node, so we - // use the same type as a Node index. - typedef NodeIndex NodeHeight; - typedef ZVector NodeHeightArray; - - // Initialize a MaxFlow instance on the given graph. The graph does not need - // to be fully built yet, but its capacity reservation are used to initialize - // the memory of this class. source and sink must also be valid node of - // graph. - GenericMaxFlow(const Graph* graph, NodeIndex source, NodeIndex sink); - -#ifndef SWIG - // This type is neither copyable nor movable. - GenericMaxFlow(const GenericMaxFlow&) = delete; - GenericMaxFlow& operator=(const GenericMaxFlow&) = delete; -#endif - - virtual ~GenericMaxFlow() {} - - // Returns the graph associated to the current object. - const Graph* graph() const { return graph_; } - - // Returns the status of last call to Solve(). NOT_SOLVED is returned if - // Solve() has never been called or if the problem has been modified in such a - // way that the previous solution becomes invalid. - Status status() const { return status_; } - - // Returns the index of the node corresponding to the source of the network. - NodeIndex GetSourceNodeIndex() const { return source_; } - - // Returns the index of the node corresponding to the sink of the network. - NodeIndex GetSinkNodeIndex() const { return sink_; } - - // Sets the capacity for arc to new_capacity. - void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity); - - // Sets the flow for arc. - void SetArcFlow(ArcIndex arc, FlowQuantity new_flow); - - // Returns true if a maximum flow was solved. - bool Solve(); - - // Returns the total flow found by the algorithm. - FlowQuantity GetOptimalFlow() const { return node_excess_[sink_]; } - - // Returns the flow on arc using the equations given in the comment on - // residual_arc_capacity_. - FlowQuantity Flow(ArcIndex arc) const { - if (IsArcDirect(arc)) { - return residual_arc_capacity_[Opposite(arc)]; - } else { - return -residual_arc_capacity_[arc]; - } - } - - // Returns the capacity of arc using the equations given in the comment on - // residual_arc_capacity_. - FlowQuantity Capacity(ArcIndex arc) const { - if (IsArcDirect(arc)) { - return residual_arc_capacity_[arc] + - residual_arc_capacity_[Opposite(arc)]; - } else { - return 0; - } - } - - // Returns the nodes reachable from the source in the residual graph, the - // outgoing arcs of this set form a minimum cut. - void GetSourceSideMinCut(std::vector* result); - - // Returns the nodes that can reach the sink in the residual graph, the - // outgoing arcs of this set form a minimum cut. Note that if this is the - // complement of GetNodeReachableFromSource(), then the min-cut is unique. - // - // TODO(user): In the two-phases algorithm, we can get this minimum cut - // without doing the second phase. Add an option for this if there is a need - // to, note that the second phase is pretty fast so the gain will be small. - void GetSinkSideMinCut(std::vector* result); - - // Checks the consistency of the input, i.e. that capacities on the arcs are - // non-negative or null. - bool CheckInputConsistency() const; - - // Checks whether the result is valid, i.e. that node excesses are all equal - // to zero (we have a flow) and that residual capacities are all non-negative - // or zero. - bool CheckResult() const; - - // Returns true if there exists a path from the source to the sink with - // remaining capacity. This allows us to easily check at the end that the flow - // we computed is indeed optimal (provided that all the conditions tested by - // CheckResult() also hold). - bool AugmentingPathExists() const; - - // Sets the different algorithm options. All default to true. - // See the corresponding variable declaration below for more details. - void SetUseGlobalUpdate(bool value) { - use_global_update_ = value; - if (!use_global_update_) process_node_by_height_ = false; - } - void SetUseTwoPhaseAlgorithm(bool value) { use_two_phase_algorithm_ = value; } - void SetCheckInput(bool value) { check_input_ = value; } - void SetCheckResult(bool value) { check_result_ = value; } - void ProcessNodeByHeight(bool value) { - process_node_by_height_ = value && use_global_update_; - } - - // Returns the protocol buffer representation of the current problem. - FlowModelProto CreateFlowModel(); - - protected: - // Returns true if arc is admissible. - bool IsAdmissible(ArcIndex arc) const { - return residual_arc_capacity_[arc] > 0 && - node_potential_[Tail(arc)] == node_potential_[Head(arc)] + 1; - } - - // Returns true if node is active, i.e. if its excess is positive and it - // is neither the source or the sink of the graph. - bool IsActive(NodeIndex node) const { - return (node != source_) && (node != sink_) && (node_excess_[node] > 0); - } - - // Sets the capacity of arc to 'capacity' and clears the flow on arc. - void SetCapacityAndClearFlow(ArcIndex arc, FlowQuantity capacity) { - residual_arc_capacity_.Set(arc, capacity); - residual_arc_capacity_.Set(Opposite(arc), 0); - } - - // Returns true if a precondition for Relabel is met, i.e. the outgoing arcs - // of node are all either saturated or the heights of their heads are greater - // or equal to the height of node. - bool CheckRelabelPrecondition(NodeIndex node) const; - - // Returns context concatenated with information about arc - // in a human-friendly way. - std::string DebugString(absl::string_view context, ArcIndex arc) const; - - // Initializes the container active_nodes_. - void InitializeActiveNodeContainer(); - - // Get the first element from the active node container. - NodeIndex GetAndRemoveFirstActiveNode() { - if (process_node_by_height_) return active_node_by_height_.Pop(); - const NodeIndex node = active_nodes_.back(); - active_nodes_.pop_back(); - return node; - } - - // Push element to the active node container. - void PushActiveNode(const NodeIndex& node) { - if (process_node_by_height_) { - active_node_by_height_.Push(node, node_potential_[node]); - } else { - active_nodes_.push_back(node); - } - } - - // Check the emptiness of the container. - bool IsEmptyActiveNodeContainer() { - if (process_node_by_height_) { - return active_node_by_height_.IsEmpty(); - } else { - return active_nodes_.empty(); - } - } - - // Performs optimization step. - void Refine(); - void RefineWithGlobalUpdate(); - - // Discharges an active node node by saturating its admissible adjacent arcs, - // if any, and by relabelling it when it becomes inactive. - void Discharge(NodeIndex node); - - // Initializes the preflow to a state that enables to run Refine. - void InitializePreflow(); - - // Clears the flow excess at each node by pushing the flow back to the source: - // - Do a depth-first search from the source in the direct graph to cancel - // flow cycles. - // - Then, return flow excess along the depth-first search tree (by pushing - // the flow in the reverse dfs topological order). - // The theoretical complexity is O(mn), but it is a lot faster in practice. - void PushFlowExcessBackToSource(); - - // Computes the best possible node potential given the current flow using a - // reverse breadth-first search from the sink in the reverse residual graph. - // This is an implementation of the global update heuristic mentioned in many - // max-flow papers. See for instance: B.V. Cherkassky, A.V. Goldberg, "On - // implementing push-relabel methods for the maximum flow problem", - // Algorithmica, 19:390-410, 1997. - // ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/94/1523/CS-TR-94-1523.pdf - void GlobalUpdate(); - - // Tries to saturate all the outgoing arcs from the source that can reach the - // sink. Most of the time, we can do that in one go, except when more flow - // than kMaxFlowQuantity can be pushed out of the source in which case we - // have to be careful. Returns true if some flow was pushed. - bool SaturateOutgoingArcsFromSource(); - - // Pushes flow on arc, i.e. consumes flow on residual_arc_capacity_[arc], - // and consumes -flow on residual_arc_capacity_[Opposite(arc)]. Updates - // node_excess_ at the tail and head of arc accordingly. - void PushFlow(FlowQuantity flow, ArcIndex arc); - - // Relabels a node, i.e. increases its height by the minimum necessary amount. - // This version of Relabel is relaxed in a way such that if an admissible arc - // exists at the current node height, then the node is not relabeled. This - // enables us to deal with wrong values of first_admissible_arc_[node] when - // updating it is too costly. - void Relabel(NodeIndex node); - - // Handy member functions to make the code more compact. - NodeIndex Head(ArcIndex arc) const { return graph_->Head(arc); } - NodeIndex Tail(ArcIndex arc) const { return graph_->Tail(arc); } - ArcIndex Opposite(ArcIndex arc) const; - bool IsArcDirect(ArcIndex arc) const; - bool IsArcValid(ArcIndex arc) const; - - // Returns the set of nodes reachable from start in the residual graph or in - // the reverse residual graph (if reverse is true). - template - void ComputeReachableNodes(NodeIndex start, std::vector* result); - - // Maximum manageable flow. - static const FlowQuantity kMaxFlowQuantity; - - // A pointer to the graph passed as argument. - const Graph* graph_; - - // An array representing the excess for each node in graph_. - QuantityArray node_excess_; - - // An array representing the height function for each node in graph_. For a - // given node, this is a lower bound on the shortest path length from this - // node to the sink in the residual network. The height of a node always goes - // up during the course of a Solve(). - // - // Since initially we saturate all the outgoing arcs of the source, we can - // never reach the sink from the source in the residual graph. Initially we - // set the height of the source to n (the number of node of the graph) and it - // never changes. If a node as an height >= n, then this node can't reach the - // sink and its height minus n is a lower bound on the shortest path length - // from this node to the source in the residual graph. - NodeHeightArray node_potential_; - - // An array representing the residual_capacity for each arc in graph_. - // Residual capacities enable one to represent the capacity and flow for all - // arcs in the graph in the following manner. - // For all arc, residual_arc_capacity_[arc] = capacity[arc] - flow[arc] - // Moreover, for reverse arcs, capacity[arc] = 0 by definition, - // Also flow[Opposite(arc)] = -flow[arc] by definition. - // Therefore: - // - for a direct arc: - // flow[arc] = 0 - flow[Opposite(arc)] - // = capacity[Opposite(arc)] - flow[Opposite(arc)] - // = residual_arc_capacity_[Opposite(arc)] - // - for a reverse arc: - // flow[arc] = -residual_arc_capacity_[arc] - // Using these facts enables one to only maintain residual_arc_capacity_, - // instead of both capacity and flow, for each direct and indirect arc. This - // reduces the amount of memory for this information by a factor 2. - QuantityArray residual_arc_capacity_; - - // An array representing the first admissible arc for each node in graph_. - ArcIndexArray first_admissible_arc_; - - // A stack used for managing active nodes in the algorithm. - // Note that the papers cited above recommend the use of a queue, but - // benchmarking so far has not proved it is better. In particular, processing - // nodes in LIFO order has better cache locality. - std::vector active_nodes_; - - // A priority queue used for managing active nodes in the algorithm. It allows - // to select the active node with highest height before each Discharge(). - // Moreover, since all pushes from this node will be to nodes with height - // greater or equal to the initial discharged node height minus one, the - // PriorityQueueWithRestrictedPush is a perfect fit. - PriorityQueueWithRestrictedPush active_node_by_height_; - - // The index of the source node in graph_. - NodeIndex source_; - - // The index of the sink node in graph_. - NodeIndex sink_; - - // The status of the problem. - Status status_; - - // BFS queue used by the GlobalUpdate() function. We do not use a C++ queue - // because we need access to the vector for different optimizations. - std::vector node_in_bfs_queue_; - std::vector bfs_queue_; - - // Whether or not to use GlobalUpdate(). - bool use_global_update_; - - // Whether or not we use a two-phase algorithm: - // 1/ Only deal with nodes that can reach the sink. At the end we know the - // value of the maximum flow and we have a min-cut. - // 2/ Call PushFlowExcessBackToSource() to obtain a max-flow. This is usually - // a lot faster than the first phase. - bool use_two_phase_algorithm_; - - // Whether or not we use the PriorityQueueWithRestrictedPush to process the - // active nodes rather than a simple queue. This can only be true if - // use_global_update_ is true. - // - // Note(user): using a template will be slightly faster, but since we test - // this in a non-critical path, this only has a minor impact. - bool process_node_by_height_; - - // Whether or not we check the input, this is a small price to pay for - // robustness. Disable only if you know the input is valid because an invalid - // input can cause the algorithm to run into an infinite loop! - bool check_input_; - - // Whether or not we check the result. - // TODO(user): Make the check more exhaustive by checking the optimality? - bool check_result_; - - // Statistics about this class. - mutable StatsGroup stats_; -}; - -#if !SWIG - -// Note: SWIG does not seem to understand explicit template specialization and -// instantiation declarations. - -template <> -const FlowQuantity GenericMaxFlow::kMaxFlowQuantity; -template <> -const FlowQuantity - GenericMaxFlow<::util::ReverseArcListGraph<>>::kMaxFlowQuantity; -template <> -const FlowQuantity - GenericMaxFlow<::util::ReverseArcStaticGraph<>>::kMaxFlowQuantity; -template <> -const FlowQuantity - GenericMaxFlow<::util::ReverseArcMixedGraph<>>::kMaxFlowQuantity; - -extern template class GenericMaxFlow; -extern template class GenericMaxFlow<::util::ReverseArcListGraph<>>; -extern template class GenericMaxFlow<::util::ReverseArcStaticGraph<>>; -extern template class GenericMaxFlow<::util::ReverseArcMixedGraph<>>; - -// This method computes a minimum vertex cover for the bipartite graph. -// -// If we define num_left=left_to_right_arcs.size(), the "left" nodes are -// integers in [0, num_left), and the "right" nodes are integers in [num_left, -// num_left + num_right). -// -// Returns a vector of size num_left+num_right, such that element #l is true if -// it is part of the minimum vertex cover and false if it is part of the maximum -// independent set (one is the complement of the other). -std::vector BipartiteMinimumVertexCover( - const std::vector>& left_to_right_arcs, int num_right); - -// Default instance MaxFlow that uses StarGraph. Note that we cannot just use a -// typedef because of dependent code expecting MaxFlow to be a real class. -// TODO(user): Modify this code and remove it. -class MaxFlow : public GenericMaxFlow { - public: - MaxFlow(const StarGraph* graph, NodeIndex source, NodeIndex target) - : GenericMaxFlow(graph, source, target) {} -}; - -#endif // SWIG - -template -bool PriorityQueueWithRestrictedPush::IsEmpty() - const { - return even_queue_.empty() && odd_queue_.empty(); -} - -template -void PriorityQueueWithRestrictedPush::Clear() { - even_queue_.clear(); - odd_queue_.clear(); -} - -template -void PriorityQueueWithRestrictedPush::Push( - Element element, IntegerPriority priority) { - // Since users may rely on it, we DCHECK the exact condition. - DCHECK(even_queue_.empty() || priority >= even_queue_.back().second - 1); - DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second - 1); - - // Note that the DCHECK() below are less restrictive than the ones above but - // check a necessary and sufficient condition for the priority queue to behave - // as expected. - if (priority & 1) { - DCHECK(odd_queue_.empty() || priority >= odd_queue_.back().second); - odd_queue_.push_back(std::make_pair(element, priority)); - } else { - DCHECK(even_queue_.empty() || priority >= even_queue_.back().second); - even_queue_.push_back(std::make_pair(element, priority)); - } -} - -template -Element PriorityQueueWithRestrictedPush::Pop() { - DCHECK(!IsEmpty()); - if (even_queue_.empty()) return PopBack(&odd_queue_); - if (odd_queue_.empty()) return PopBack(&even_queue_); - if (odd_queue_.back().second > even_queue_.back().second) { - return PopBack(&odd_queue_); - } else { - return PopBack(&even_queue_); - } -} - -template -Element PriorityQueueWithRestrictedPush::PopBack( - std::vector>* queue) { - DCHECK(!queue->empty()); - Element element = queue->back().first; - queue->pop_back(); - return element; -} - } // namespace operations_research #endif // OR_TOOLS_GRAPH_MAX_FLOW_H_ diff --git a/ortools/graph/max_flow_test.cc b/ortools/graph/max_flow_test.cc index 8c7d6b5b0a..7930155a9f 100644 --- a/ortools/graph/max_flow_test.cc +++ b/ortools/graph/max_flow_test.cc @@ -13,38 +13,22 @@ #include "ortools/graph/max_flow.h" -#include -#include #include -#include -#include -#include -#include -#include "absl/algorithm/container.h" -#include "absl/random/random.h" -#include "absl/strings/str_format.h" -#include "absl/types/span.h" -#include "benchmark/benchmark.h" #include "google/protobuf/text_format.h" #include "gtest/gtest.h" #include "ortools/base/gmock.h" -#include "ortools/base/logging.h" #include "ortools/base/path.h" -#include "ortools/graph/graph.h" -#include "ortools/graph/graphs.h" -#include "ortools/linear_solver/linear_solver.h" +#include "ortools/graph/ebert_graph.h" +#include "ortools/graph/flow_problem.pb.h" #include "ortools/util/file_util.h" -#if not defined(ROOT_DIR) -#define ROOT_DIR "_main/" -#endif - namespace operations_research { namespace { -using ::testing::ContainerEq; -using ::testing::WhenSorted; +#if not defined(ROOT_DIR) +#define ROOT_DIR "_main/" +#endif TEST(SimpleMaxFlowTest, EmptyWithValidSourceAndSink) { SimpleMaxFlow max_flow; @@ -198,681 +182,5 @@ TEST(SimpleMaxFlowTest, ProblematicProblemWithMaxCapacity) { EXPECT_EQ(10290243, solver.OptimalFlow()); } -template -typename GenericMaxFlow::Status MaxFlowTester( - const typename Graph::NodeIndex num_nodes, - const typename Graph::ArcIndex num_arcs, - const typename Graph::NodeIndex* tail, - const typename Graph::NodeIndex* head, const FlowQuantity* capacity, - const FlowQuantity* expected_flow, const FlowQuantity expected_total_flow, - const std::vector* expected_source_min_cut = - nullptr, - const std::vector* expected_sink_min_cut = - nullptr) { - Graph graph(num_nodes, num_arcs); - for (int i = 0; i < num_arcs; ++i) { - graph.AddArc(tail[i], head[i]); - } - std::vector permutation; - Graphs::Build(&graph, &permutation); - EXPECT_TRUE(permutation.empty()); - - typename GenericMaxFlow::Status result = - GenericMaxFlow::OPTIMAL; - for (int options = 0; options < 8; ++options) { - GenericMaxFlow max_flow(&graph, 0, num_nodes - 1); - max_flow.SetUseGlobalUpdate(options & 1); - max_flow.SetUseTwoPhaseAlgorithm(options & 2); - max_flow.ProcessNodeByHeight(options & 3); - for (typename Graph::ArcIndex arc = 0; arc < num_arcs; ++arc) { - max_flow.SetArcCapacity(arc, capacity[arc]); - EXPECT_EQ(max_flow.Capacity(arc), capacity[arc]); - } - EXPECT_TRUE(max_flow.CheckInputConsistency()); - EXPECT_TRUE(max_flow.Solve()); - if (max_flow.status() == GenericMaxFlow::OPTIMAL) { - const FlowQuantity total_flow = max_flow.GetOptimalFlow(); - EXPECT_EQ(expected_total_flow, total_flow); - for (int i = 0; i < num_arcs; ++i) { - EXPECT_EQ(expected_flow[i], max_flow.Flow(i)) << " i = " << i; - } - } - EXPECT_TRUE(max_flow.CheckResult()); - if (options == 0) { - result = max_flow.status(); - } else { - EXPECT_EQ(result, max_flow.status()); - } - - // Tests the min-cut functions. - if (expected_source_min_cut != nullptr) { - std::vector cut; - max_flow.GetSourceSideMinCut(&cut); - std::sort(cut.begin(), cut.end()); - EXPECT_THAT(*expected_source_min_cut, WhenSorted(ContainerEq(cut))); - } - if (expected_sink_min_cut != nullptr) { - std::vector cut; - max_flow.GetSinkSideMinCut(&cut); - std::sort(cut.begin(), cut.end()); - EXPECT_THAT(*expected_sink_min_cut, WhenSorted(ContainerEq(cut))); - } - } - - return result; -} - -template -class GenericMaxFlowTest : public ::testing::Test {}; - -typedef ::testing::Types, - util::ReverseArcStaticGraph<>, - util::ReverseArcMixedGraph<>> - GraphTypes; - -TYPED_TEST_SUITE(GenericMaxFlowTest, GraphTypes); - -TYPED_TEST(GenericMaxFlowTest, FeasibleFlow1) { - const int kNumNodes = 4; - const int kNumArcs = 3; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 1, 2}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3}; - const FlowQuantity kCapacity[kNumArcs] = {8, 10, 8}; - const FlowQuantity kExpectedFlow[kNumArcs] = {8, 8, 8}; - const FlowQuantity kExpectedTotalFlow = 8; - std::vector source_cut({0}); - std::vector sink_cut({3}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -TYPED_TEST(GenericMaxFlowTest, FeasibleFlow2) { - const int kNumNodes = 6; - const int kNumArcs = 9; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 0, 0, 1, - 2, 3, 3, 4}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 4, 3, - 4, 4, 5, 5}; - const FlowQuantity kCapacity[kNumArcs] = {6, 8, 5, 0, 1, 4, 0, 6, 4}; - const FlowQuantity kExpectedFlow[kNumArcs] = {1, 4, 5, 0, 1, 4, 0, 6, 4}; - const FlowQuantity kExpectedTotalFlow = 10; - std::vector source_cut({0, 1, 2}); - std::vector sink_cut({5}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -TYPED_TEST(GenericMaxFlowTest, FeasibleFlowWithMultipleArcs) { - const int kNumNodes = 5; - const int kNumArcs = 8; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 1, - 2, 2, 3, 3}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 1, 2, 2, - 3, 3, 4, 4}; - const FlowQuantity kCapacity[kNumArcs] = {5, 3, 5, 3, 4, 4, 4, 4}; - const FlowQuantity kExpectedFlow[kNumArcs] = {5, 3, 5, 3, 4, 4, 4, 4}; - const FlowQuantity kExpectedTotalFlow = 8; - std::vector source_cut({0}); - std::vector sink_cut({4}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -TYPED_TEST(GenericMaxFlowTest, HugeCapacity) { - const FlowQuantity kCapacityMax = std::numeric_limits::max(); - const int kNumNodes = 5; - const int kNumArcs = 5; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 2, 3}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 3, 4}; - const FlowQuantity kCapacity[kNumArcs] = {kCapacityMax, kCapacityMax, 5, 3, - kCapacityMax}; - const FlowQuantity kExpectedFlow[kNumArcs] = {5, 3, 5, 3, 8}; - const FlowQuantity kExpectedTotalFlow = 8; - std::vector source_cut({0, 1, 2}); - std::vector sink_cut({4, 3}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -TYPED_TEST(GenericMaxFlowTest, FlowQuantityOverflowLimitCase) { - const FlowQuantity kCapacityMax = std::numeric_limits::max(); - const FlowQuantity kHalfLow = kCapacityMax / 2; - const FlowQuantity kHalfHigh = kCapacityMax - kHalfLow; - const int kNumNodes = 5; - const int kNumArcs = 5; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 2, 3}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 3, 4}; - const FlowQuantity kCapacity[kNumArcs] = {kCapacityMax, kCapacityMax, - kHalfLow, kHalfHigh, kCapacityMax}; - const FlowQuantity kExpectedFlow[kNumArcs] = {kHalfLow, kHalfHigh, kHalfLow, - kHalfHigh, kCapacityMax}; - const FlowQuantity kExpectedTotalFlow = kCapacityMax; - std::vector source_cut({0, 1, 2}); - std::vector sink_cut({4}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -TYPED_TEST(GenericMaxFlowTest, FlowQuantityOverflow) { - const FlowQuantity kCapacityMax = std::numeric_limits::max(); - const int kNumNodes = 4; - const int kNumArcs = 4; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 1, 2}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 3}; - const FlowQuantity kCapacity[kNumArcs] = {kCapacityMax, kCapacityMax, - kCapacityMax, kCapacityMax}; - const FlowQuantity kExpectedFlow[kNumArcs] = {kCapacityMax, kCapacityMax, - kCapacityMax, kCapacityMax}; - const FlowQuantity kExpectedTotalFlow = kCapacityMax; - EXPECT_EQ( - GenericMaxFlow::INT_OVERFLOW, - MaxFlowTester(kNumNodes, kNumArcs, kTail, kHead, kCapacity, - kExpectedFlow, kExpectedTotalFlow)); -} - -TYPED_TEST(GenericMaxFlowTest, DirectArcFromSourceToSink) { - const int kNumNodes = 4; - const int kNumArcs = 5; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 0, 1, 2}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 3, 2, 3, 3}; - const FlowQuantity kCapacity[kNumArcs] = {5, 8, 5, 2, 2}; - const FlowQuantity kExpectedFlow[kNumArcs] = {2, 8, 2, 2, 2}; - const FlowQuantity kExpectedTotalFlow = 12; - std::vector source_cut({0, 1, 2}); - std::vector sink_cut({3}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -TYPED_TEST(GenericMaxFlowTest, FlowOnDisconnectedGraph1) { - const int kNumNodes = 6; - const int kNumArcs = 7; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 0, 0, 1, 2, 3}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 3, 4, 3, 4, 4}; - const FlowQuantity kCapacity[kNumArcs] = {5, 8, 5, 3, 4, 5, 6}; - const FlowQuantity kExpectedFlow[kNumArcs] = {0, 0, 0, 0, 0, 0, 0}; - const FlowQuantity kExpectedTotalFlow = 0; - std::vector source_cut({0, 1, 2, 3, 4}); - std::vector sink_cut({5}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -TYPED_TEST(GenericMaxFlowTest, FlowOnDisconnectedGraph2) { - const int kNumNodes = 6; - const int kNumArcs = 5; - const typename TypeParam::NodeIndex kTail[kNumArcs] = {0, 0, 3, 3, 4}; - const typename TypeParam::NodeIndex kHead[kNumArcs] = {1, 2, 4, 5, 5}; - const FlowQuantity kCapacity[kNumArcs] = {5, 8, 6, 6, 4}; - const FlowQuantity kExpectedFlow[kNumArcs] = {0, 0, 0, 0, 0}; - const FlowQuantity kExpectedTotalFlow = 0; - std::vector source_cut({0, 1, 2}); - std::vector sink_cut({3, 4, 5}); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, - MaxFlowTester( - kNumNodes, kNumArcs, kTail, kHead, kCapacity, kExpectedFlow, - kExpectedTotalFlow, &source_cut, &sink_cut)); -} - -template -void AddSourceAndSink(const typename Graph::NodeIndex num_tails, - const typename Graph::NodeIndex num_heads, Graph* graph) { - const typename Graph::NodeIndex source = num_tails + num_heads; - const typename Graph::NodeIndex sink = num_tails + num_heads + 1; - for (typename Graph::NodeIndex tail = 0; tail < num_tails; ++tail) { - graph->AddArc(source, tail); - } - for (typename Graph::NodeIndex head = 0; head < num_heads; ++head) { - graph->AddArc(num_tails + head, sink); - } -} - -template -void GenerateCompleteGraph(const typename Graph::NodeIndex num_tails, - const typename Graph::NodeIndex num_heads, - Graph* graph) { - const typename Graph::NodeIndex num_nodes = num_tails + num_heads + 2; - const typename Graph::ArcIndex num_arcs = - num_tails * num_heads + num_tails + num_heads; - graph->Reserve(num_nodes, num_arcs); - for (typename Graph::NodeIndex tail = 0; tail < num_tails; ++tail) { - for (typename Graph::NodeIndex head = 0; head < num_heads; ++head) { - graph->AddArc(tail, head + num_tails); - } - } - AddSourceAndSink(num_tails, num_heads, graph); -} - -template -void GeneratePartialRandomGraph(const typename Graph::NodeIndex num_tails, - const typename Graph::NodeIndex num_heads, - const typename Graph::NodeIndex degree, - Graph* graph) { - const typename Graph::NodeIndex num_nodes = num_tails + num_heads + 2; - const typename Graph::ArcIndex num_arcs = - num_tails * degree + num_tails + num_heads; - graph->Reserve(num_nodes, num_arcs); - std::mt19937 randomizer(0); - for (typename Graph::NodeIndex tail = 0; tail < num_tails; ++tail) { - for (typename Graph::NodeIndex d = 0; d < degree; ++d) { - const typename Graph::NodeIndex head = - absl::Uniform(randomizer, 0, num_heads); - graph->AddArc(tail, head + num_tails); - } - } - AddSourceAndSink(num_tails, num_heads, graph); -} - -template -void GenerateRandomArcValuations(const Graph& graph, const int64_t max_range, - std::vector* arc_valuation) { - const typename Graph::ArcIndex num_arcs = graph.num_arcs(); - arc_valuation->resize(num_arcs); - std::mt19937 randomizer(0); - for (typename Graph::ArcIndex arc = 0; arc < graph.num_arcs(); ++arc) { - (*arc_valuation)[arc] = absl::Uniform(randomizer, 0, max_range); - } -} - -template -void SetUpNetworkData(absl::Span arc_capacity, - GenericMaxFlow* max_flow) { - const Graph* graph = max_flow->graph(); - for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { - max_flow->SetArcCapacity(arc, arc_capacity[arc]); - } -} - -template -FlowQuantity SolveMaxFlow(GenericMaxFlow* max_flow) { - EXPECT_TRUE(max_flow->Solve()); - EXPECT_EQ(GenericMaxFlow::OPTIMAL, max_flow->status()); - const Graph* graph = max_flow->graph(); - for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { - EXPECT_LE(max_flow->Flow(Graphs::OppositeArc(*graph, arc)), 0) - << arc; - EXPECT_EQ(0, max_flow->Capacity(Graphs::OppositeArc(*graph, arc))) - << arc; - EXPECT_LE(0, max_flow->Flow(arc)) << arc; - EXPECT_LE(max_flow->Flow(arc), max_flow->Capacity(arc)) << arc; - } - return max_flow->GetOptimalFlow(); -} - -template -FlowQuantity SolveMaxFlowWithLP(GenericMaxFlow* max_flow) { - MPSolver solver("LPSolver", MPSolver::GLOP_LINEAR_PROGRAMMING); - const double infinity = solver.infinity(); - const Graph* graph = max_flow->graph(); - const typename Graph::NodeIndex num_nodes = graph->num_nodes(); - const typename Graph::ArcIndex num_arcs = graph->num_arcs(); - const typename Graph::NodeIndex source_index = num_nodes - 2; - std::unique_ptr constraint(new MPConstraint*[num_nodes]); - for (typename Graph::NodeIndex node = 0; node < graph->num_nodes(); ++node) { - constraint[node] = solver.MakeRowConstraint(); - if (node < source_index) { // Node is neither source nor sink. - constraint[node]->SetBounds(0.0, 0.0); // Flow is conserved. - } else { - constraint[node]->SetBounds(-infinity, +infinity); - } - } - std::unique_ptr var(new MPVariable*[num_arcs]); - for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { - var[arc] = solver.MakeNumVar(0.0, max_flow->Capacity(arc), - absl::StrFormat("v%d", arc)); - constraint[graph->Tail(arc)]->SetCoefficient(var[arc], 1.0); - constraint[graph->Head(arc)]->SetCoefficient(var[arc], -1.0); - } - MPObjective* const objective = solver.MutableObjective(); - for (typename Graph::OutgoingArcIterator arc_it(*graph, source_index); - arc_it.Ok(); arc_it.Next()) { - const typename Graph::ArcIndex arc = arc_it.Index(); - objective->SetCoefficient(var[arc], -1.0); - } - solver.Solve(); - return static_cast(-objective->Value() + .5); -} - -template -struct MaxFlowSolver { - typedef FlowQuantity (*Solver)(GenericMaxFlow* max_flow); -}; - -template -void FullRandomAssignment(typename MaxFlowSolver::Solver f, - typename Graph::NodeIndex num_tails, - typename Graph::NodeIndex num_heads, - FlowQuantity expected_flow1, - FlowQuantity expected_flow2) { - Graph graph; - GenerateCompleteGraph(num_tails, num_heads, &graph); - Graphs::Build(&graph); - std::vector arc_capacity(graph.num_arcs(), 1); - std::unique_ptr> max_flow(new GenericMaxFlow( - &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); - SetUpNetworkData(arc_capacity, max_flow.get()); - FlowQuantity flow = f(max_flow.get()); - EXPECT_EQ(expected_flow1, flow); -} - -template -void PartialRandomAssignment(typename MaxFlowSolver::Solver f, - typename Graph::NodeIndex num_tails, - typename Graph::NodeIndex num_heads, - FlowQuantity expected_flow1, - FlowQuantity expected_flow2) { - const typename Graph::NodeIndex kDegree = 10; - Graph graph; - GeneratePartialRandomGraph(num_tails, num_heads, kDegree, &graph); - Graphs::Build(&graph); - CHECK_EQ(graph.num_arcs(), num_tails * kDegree + num_tails + num_heads); - std::vector arc_capacity(graph.num_arcs(), 1); - std::unique_ptr> max_flow(new GenericMaxFlow( - &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); - SetUpNetworkData(arc_capacity, max_flow.get()); - FlowQuantity flow = f(max_flow.get()); - EXPECT_EQ(expected_flow1, flow); -} - -template -void ChangeCapacities(absl::Span arc_capacity, - FlowQuantity delta, GenericMaxFlow* max_flow) { - const Graph* graph = max_flow->graph(); - for (typename Graph::ArcIndex arc = 0; arc < graph->num_arcs(); ++arc) { - max_flow->SetArcCapacity(arc, - std::max(arc_capacity[arc] - delta, int64_t{0})); - } -} - -template -void PartialRandomFlow(typename MaxFlowSolver::Solver f, - typename Graph::NodeIndex num_tails, - typename Graph::NodeIndex num_heads, - FlowQuantity expected_flow1, - FlowQuantity expected_flow2) { - const typename Graph::NodeIndex kDegree = 10; - const FlowQuantity kCapacityRange = 10000; - const FlowQuantity kCapacityDelta = 1000; - Graph graph; - GeneratePartialRandomGraph(num_tails, num_heads, kDegree, &graph); - std::vector arc_capacity(graph.num_arcs()); - GenerateRandomArcValuations(graph, kCapacityRange, &arc_capacity); - - std::vector permutation; - Graphs::Build(&graph, &permutation); - util::Permute(permutation, &arc_capacity); - - std::unique_ptr> max_flow(new GenericMaxFlow( - &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); - SetUpNetworkData(arc_capacity, max_flow.get()); - FlowQuantity flow = f(max_flow.get()); - EXPECT_EQ(expected_flow1, flow); - ChangeCapacities(arc_capacity, kCapacityDelta, max_flow.get()); - flow = f(max_flow.get()); - EXPECT_EQ(expected_flow2, flow); - ChangeCapacities(arc_capacity, 0, max_flow.get()); - flow = f(max_flow.get()); - EXPECT_EQ(expected_flow1, flow); -} - -template -void FullRandomFlow(typename MaxFlowSolver::Solver f, - typename Graph::NodeIndex num_tails, - typename Graph::NodeIndex num_heads, - FlowQuantity expected_flow1, FlowQuantity expected_flow2) { - const FlowQuantity kCapacityRange = 10000; - const FlowQuantity kCapacityDelta = 1000; - Graph graph; - GenerateCompleteGraph(num_tails, num_heads, &graph); - std::vector arc_capacity(graph.num_arcs()); - GenerateRandomArcValuations(graph, kCapacityRange, &arc_capacity); - - std::vector permutation; - Graphs::Build(&graph, &permutation); - util::Permute(permutation, &arc_capacity); - - std::unique_ptr> max_flow(new GenericMaxFlow( - &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); - SetUpNetworkData(arc_capacity, max_flow.get()); - FlowQuantity flow = f(max_flow.get()); - EXPECT_EQ(expected_flow1, flow); - ChangeCapacities(arc_capacity, kCapacityDelta, max_flow.get()); - flow = f(max_flow.get()); - EXPECT_EQ(expected_flow2, flow); - ChangeCapacities(arc_capacity, 0, max_flow.get()); - flow = f(max_flow.get()); - EXPECT_EQ(expected_flow1, flow); -} - -#define LP_AND_FLOW_TEST(test_name, size, expected_flow1, expected_flow2) \ - LP_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ - FLOW_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ - FLOW_ONLY_TEST_SG(test_name, size, expected_flow1, expected_flow2) - -#define LP_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ - TEST(LPMaxFlowTest, test_name##size) { \ - test_name(SolveMaxFlowWithLP, size, size, \ - expected_flow1, expected_flow2); \ - } - -#define FLOW_ONLY_TEST(test_name, size, expected_flow1, expected_flow2) \ - TEST(MaxFlowTest, test_name##size) { \ - test_name(SolveMaxFlow, size, size, expected_flow1, \ - expected_flow2); \ - } - -#define FLOW_ONLY_TEST_SG(test_name, size, expected_flow1, expected_flow2) \ - TEST(MaxFlowTestStaticGraph, test_name##size) { \ - test_name>(SolveMaxFlow, size, size, \ - expected_flow1, expected_flow2); \ - } - -LP_AND_FLOW_TEST(FullRandomAssignment, 300, 300, 300); -LP_AND_FLOW_TEST(PartialRandomAssignment, 100, 100, 100); -LP_AND_FLOW_TEST(PartialRandomAssignment, 1000, 1000, 1000); -LP_AND_FLOW_TEST(PartialRandomFlow, 400, 1898664, 1515203); -LP_AND_FLOW_TEST(FullRandomFlow, 100, 482391, 386587); - -// LARGE must be defined from the build command line to test larger instances. -#ifdef LARGE -LP_AND_FLOW_TEST(PartialRandomAssignment, 10000, 10000, 10000); -#endif - -template -static void BM_FullRandomAssignment(benchmark::State& state) { - const int kSize = 3000; - for (auto _ : state) { - FullRandomAssignment(SolveMaxFlow, kSize, kSize, kSize, kSize); - } - state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); -} - -template -static void BM_PartialRandomAssignment(benchmark::State& state) { - const int kSize = 10100; - for (auto _ : state) { - PartialRandomAssignment(SolveMaxFlow, kSize, kSize, kSize, kSize); - } - state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); -} - -template -static void BM_PartialRandomFlow(benchmark::State& state) { - const int kSize = 800; - for (auto _ : state) { - PartialRandomFlow(SolveMaxFlow, kSize, kSize, 3884850, 3112123); - } - state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); -} - -template -static void BM_FullRandomFlow(benchmark::State& state) { - const int kSize = 800; - for (auto _ : state) { - FullRandomFlow(SolveMaxFlow, kSize, kSize, 4000549, 3239512); - } - state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); -} - -BENCHMARK_TEMPLATE(BM_FullRandomAssignment, StarGraph); -BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcListGraph<>); -BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcStaticGraph<>); -BENCHMARK_TEMPLATE(BM_FullRandomAssignment, util::ReverseArcMixedGraph<>); -BENCHMARK_TEMPLATE(BM_PartialRandomFlow, StarGraph); -BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcListGraph<>); -BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcStaticGraph<>); -BENCHMARK_TEMPLATE(BM_PartialRandomFlow, util::ReverseArcMixedGraph<>); - -// One iteration of each of the following tests is slow. -BENCHMARK_TEMPLATE(BM_FullRandomFlow, StarGraph); -BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcListGraph<>); -BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcStaticGraph<>); -BENCHMARK_TEMPLATE(BM_FullRandomFlow, util::ReverseArcMixedGraph<>); -BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, StarGraph); -BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcListGraph<>); -BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcStaticGraph<>); -BENCHMARK_TEMPLATE(BM_PartialRandomAssignment, util::ReverseArcMixedGraph<>); - -#undef LP_AND_FLOW_TEST -#undef LP_ONLY_TEST -#undef FLOW_ONLY_TEST -#undef FLOW_ONLY_TEST_SG - -// ---------------------------------------------------------- -// PriorityQueueWithRestrictedPush tests. -// ---------------------------------------------------------- - -TEST(PriorityQueueWithRestrictedPushTest, BasicBehavior) { - PriorityQueueWithRestrictedPush queue; - EXPECT_TRUE(queue.IsEmpty()); - queue.Push("A", 1); - queue.Push("B", 0); - queue.Push("C", 2); - queue.Push("D", 10); - queue.Push("E", 9); - EXPECT_EQ("D", queue.Pop()); - EXPECT_EQ("E", queue.Pop()); - EXPECT_EQ("C", queue.Pop()); - EXPECT_EQ("A", queue.Pop()); - EXPECT_EQ("B", queue.Pop()); - EXPECT_TRUE(queue.IsEmpty()); - queue.Push("A", 1); - queue.Push("B", 0); - EXPECT_FALSE(queue.IsEmpty()); - queue.Clear(); - EXPECT_TRUE(queue.IsEmpty()); -} - -TEST(PriorityQueueWithRestrictedPushTest, BasicBehaviorWithMixedPushPop) { - PriorityQueueWithRestrictedPush queue; - EXPECT_TRUE(queue.IsEmpty()); - queue.Push("A", 1); - queue.Push("B", 0); - queue.Push("C", 2); - EXPECT_EQ("C", queue.Pop()); - EXPECT_EQ("A", queue.Pop()); - queue.Push("D", 1); - queue.Push("E", 0); - EXPECT_EQ("D", queue.Pop()); - EXPECT_EQ("E", queue.Pop()); - EXPECT_EQ("B", queue.Pop()); - EXPECT_TRUE(queue.IsEmpty()); - queue.Push("E", 1); - EXPECT_FALSE(queue.IsEmpty()); - EXPECT_EQ("E", queue.Pop()); - EXPECT_TRUE(queue.IsEmpty()); -} - -TEST(PriorityQueueWithRestrictedPushTest, RandomPushPop) { - struct ElementWithPriority { - ElementWithPriority(int _element, int _priority) - : element(_element), priority(_priority) {} - bool operator<(const ElementWithPriority& other) const { - return priority < other.priority; - } - int element; - int priority; - }; - std::vector pairs; - std::mt19937 randomizer(1); - const int kNumElement = 10000; - const int kMaxPriority = 10000; // We want duplicates and gaps. - for (int i = 0; i < kNumElement; ++i) { - pairs.push_back( - ElementWithPriority(i, absl::Uniform(randomizer, 0, kMaxPriority))); - } - std::sort(pairs.begin(), pairs.end()); - - // Randomly add +1 and push to the queue. - PriorityQueueWithRestrictedPush queue; - for (int i = 0; i < pairs.size(); ++i) { - pairs[i].priority += absl::Bernoulli(randomizer, 1.0 / 2) ? 1 : 0; - queue.Push(pairs[i].element, pairs[i].priority); - } - - // Stable sort the pairs for checking (the queue order is stable). - std::stable_sort(pairs.begin(), pairs.end()); - - // Random Push() and Pop() with more Pop(). - int current = pairs.size(); - while (!queue.IsEmpty()) { - EXPECT_GT(current, 0); - if (absl::Bernoulli(randomizer, 1.0 / 4) && current < pairs.size()) { - queue.Push(pairs[current].element, pairs[current].priority); - ++current; - } else { - --current; - EXPECT_EQ(pairs[current].element, queue.Pop()); - } - } -} - -TEST(BipartiteMinimumVertexCoverTest, BasicBehavior) { - const int num_right = 4; - const std::vector> left_to_right = { - {5}, {4, 5, 6}, {5}, {5, 6, 7}}; - EXPECT_EQ(absl::c_count(BipartiteMinimumVertexCover(left_to_right, num_right), - true), - 3); - EXPECT_EQ(absl::c_count(BipartiteMinimumVertexCover(left_to_right, num_right), - false), - 5); -} - -TEST(BipartiteMinimumVertexCoverTest, Empty) { - const int num_right = 4; - const std::vector> left_to_right = {{}, {}}; - EXPECT_EQ(absl::c_count(BipartiteMinimumVertexCover(left_to_right, num_right), - false), - 6); -} - -TEST(PriorityQueueWithRestrictedPushDeathTest, DCHECK) { - // Don't run this test in opt mode. - if (!DEBUG_MODE) GTEST_SKIP(); - - PriorityQueueWithRestrictedPush queue; - EXPECT_TRUE(queue.IsEmpty()); - ASSERT_DEATH(queue.Pop(), ""); - queue.Push("A", 10); - queue.Push("B", 9); - ASSERT_DEATH(queue.Push("C", 4), ""); - ASSERT_DEATH(queue.Push("C", 8), ""); -} - } // namespace } // namespace operations_research diff --git a/ortools/graph/min_cost_flow.cc b/ortools/graph/min_cost_flow.cc index 7ddfed5e28..2cb528d041 100644 --- a/ortools/graph/min_cost_flow.cc +++ b/ortools/graph/min_cost_flow.cc @@ -22,9 +22,10 @@ #include #include "absl/flags/flag.h" +#include "absl/log/check.h" #include "absl/strings/str_format.h" #include "absl/strings/string_view.h" -#include "ortools/base/dump_vars.h" +#include "ortools/base/logging.h" #include "ortools/base/mathutil.h" #include "ortools/graph/graph.h" #include "ortools/graph/graphs.h" @@ -39,8 +40,6 @@ ABSL_FLAG(bool, min_cost_flow_check_feasibility, true, "Check that the graph has enough capacity to send all supplies " "and serve all demands. Also check that the sum of supplies " "is equal to the sum of demands."); -ABSL_FLAG(bool, min_cost_flow_check_balance, true, - "Check that the sum of supplies is equal to the sum of demands."); ABSL_FLAG(bool, min_cost_flow_check_result, true, "Check that the result is valid."); @@ -59,7 +58,6 @@ GenericMinCostFlow::GenericMinCostFlow( node_potential_.assign(max_num_nodes, 0); node_excess_.assign(max_num_nodes, 0); initial_node_excess_.assign(max_num_nodes, 0); - feasible_node_excess_.assign(max_num_nodes, 0); } const ArcIndex max_num_arcs = Graphs::ArcReservation(*graph_); if (max_num_arcs > 0) { @@ -127,49 +125,97 @@ void GenericMinCostFlow::SetArcCapacity( } } -template -void GenericMinCostFlow::SetArcFlow( - ArcIndex arc, ArcFlowType new_flow) { - DCHECK(IsArcValid(arc)); - const FlowQuantity capacity = Capacity(arc); - DCHECK_GE(capacity, new_flow); - residual_arc_capacity_.Set(Opposite(arc), new_flow); - residual_arc_capacity_.Set(arc, capacity - new_flow); - status_ = NOT_SOLVED; - feasibility_checked_ = false; -} - template bool GenericMinCostFlow::CheckInputConsistency() { - FlowQuantity total_supply = 0; - uint64_t max_capacity = 0; // uint64_t because it is positive and will be - // used to check against FlowQuantity overflows. - for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) { - const uint64_t capacity = - static_cast(residual_arc_capacity_[arc]); - max_capacity = std::max(capacity, max_capacity); - } - uint64_t total_flow = 0; // uint64_t for the same reason as max_capacity. + const FlowQuantity kMaxFlow = std::numeric_limits::max(); + + // First lets make sure supply == demand and the total supply/demand do not + // overflow. + FlowQuantity sum_of_supplies = 0; + FlowQuantity sum_of_demands = 0; for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) { const FlowQuantity excess = node_excess_[node]; - total_supply += excess; if (excess > 0) { - total_flow += excess; - if (std::numeric_limits::max() < - max_capacity + total_flow) { - status_ = BAD_COST_RANGE; - LOG(ERROR) << "Input consistency error: max capacity + flow exceed " - << "precision"; + sum_of_supplies = CapAdd(sum_of_supplies, excess); + if (sum_of_supplies >= kMaxFlow) { + status_ = BAD_CAPACITY_RANGE; + LOG(ERROR) << "Input consistency error: sum of supplies overflow"; + return false; + } + } else if (excess < 0) { + sum_of_demands = CapAdd(sum_of_demands, -excess); + if (sum_of_demands >= kMaxFlow) { + status_ = BAD_CAPACITY_RANGE; + LOG(ERROR) << "Input consistency error: sum of demands overflow"; return false; } } } - if (total_supply != 0) { + if (sum_of_supplies != sum_of_demands) { status_ = UNBALANCED; LOG(ERROR) << "Input consistency error: unbalanced problem"; return false; } + + std::vector max_node_excess = node_excess_; + std::vector min_node_excess = node_excess_; + for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) { + const FlowQuantity capacity = residual_arc_capacity_[arc]; + CHECK_GE(capacity, 0); + if (capacity == 0) continue; + const int tail = graph_->Tail(arc); + const int head = graph_->Head(arc); + min_node_excess[tail] = CapSub(min_node_excess[tail], capacity); + max_node_excess[head] = CapAdd(max_node_excess[head], capacity); + } + + const int num_nodes = graph_->num_nodes(); + for (NodeIndex node = 0; node < num_nodes; ++node) { + if (max_node_excess[node] >= kMaxFlow || + min_node_excess[node] <= -kMaxFlow) { + // Try to fix it. + // Some user just use arc with infinite capacity out of the source/sink. + // This is why we use CappedCapacity() for a name. + if (max_node_excess[node] < std::numeric_limits::max()) { + // There is no point having an outgoing arc with more than this. + const ArcFlowType upper_bound = + std::max(0, max_node_excess[node]); + + // Adjust and recompute min_node_excess[node]. + min_node_excess[node] = node_excess_[node]; + for (OutgoingArcIterator it(*graph_, node); it.Ok(); it.Next()) { + const int arc = it.Index(); + residual_arc_capacity_[arc] = + std::min(residual_arc_capacity_[arc], upper_bound); + min_node_excess[node] = + CapSub(min_node_excess[node], residual_arc_capacity_[arc]); + } + if (min_node_excess[node] > -kMaxFlow) continue; + } + if (min_node_excess[node] > -std::numeric_limits::max()) { + // There is no point having an incoming arc with more than this. + const ArcFlowType upper_bound = + std::max(0, -min_node_excess[node]); + + // Adjust and recompute max_node_excess[node]. + max_node_excess[node] = node_excess_[node]; + for (IncomingArcIterator it(*graph_, node); it.Ok(); it.Next()) { + const int arc = it.Index(); + residual_arc_capacity_[arc] = + std::min(residual_arc_capacity_[arc], upper_bound); + max_node_excess[node] = + CapAdd(max_node_excess[node], residual_arc_capacity_[arc]); + } + if (max_node_excess[node] < kMaxFlow) continue; + } + + status_ = BAD_CAPACITY_RANGE; + LOG(ERROR) << "Maximum in or out flow of node + excess " << node + << " overflow the FlowQuantity type (int64_t)."; + return false; + } + } return true; } @@ -248,9 +294,8 @@ GenericMinCostFlow::DebugString( } template -bool GenericMinCostFlow:: - CheckFeasibility(std::vector* const infeasible_supply_node, - std::vector* const infeasible_demand_node) { +bool GenericMinCostFlow::CheckFeasibility() { SCOPED_TIME_STAT(&stats_); // Create a new graph, which is a copy of graph_, with the following // modifications: @@ -276,8 +321,6 @@ bool GenericMinCostFlow:: const NodeIndex sink = num_nodes_in_max_flow - 1; StarGraph checker_graph(num_nodes_in_max_flow, num_arcs_in_max_flow); MaxFlow checker(&checker_graph, source, sink); - checker.SetCheckInput(false); - checker.SetCheckResult(false); // Copy graph_ to checker_graph. for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) { const ArcIndex new_arc = @@ -310,44 +353,10 @@ bool GenericMinCostFlow:: return false; } const FlowQuantity optimal_max_flow = checker.GetOptimalFlow(); - feasible_node_excess_.assign(checker_graph.num_nodes(), 0); - for (StarGraph::OutgoingArcIterator it(checker_graph, source); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); - const NodeIndex node = checker_graph.Head(arc); - const FlowQuantity flow = checker.Flow(arc); - feasible_node_excess_[node] = flow; - if (infeasible_supply_node != nullptr) { - infeasible_supply_node->push_back(node); - } - } - for (StarGraph::IncomingArcIterator it(checker_graph, sink); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); - const NodeIndex node = checker_graph.Tail(arc); - const FlowQuantity flow = checker.Flow(arc); - feasible_node_excess_[node] = -flow; - if (infeasible_demand_node != nullptr) { - infeasible_demand_node->push_back(node); - } - } feasibility_checked_ = true; return optimal_max_flow == total_supply; } -template -bool GenericMinCostFlow::MakeFeasible() { - if (!feasibility_checked_) { - return false; - } - for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) { - const FlowQuantity excess = feasible_node_excess_[node]; - node_excess_[node] = excess; - initial_node_excess_[node] = excess; - } - return true; -} - template FlowQuantity GenericMinCostFlow::Flow( ArcIndex arc) const { @@ -385,20 +394,6 @@ FlowQuantity GenericMinCostFlow::Supply( return node_excess_[node]; } -template -FlowQuantity -GenericMinCostFlow::InitialSupply( - NodeIndex node) const { - return initial_node_excess_[node]; -} - -template -FlowQuantity -GenericMinCostFlow::FeasibleSupply( - NodeIndex node) const { - return feasible_node_excess_[node]; -} - template bool GenericMinCostFlow::IsAdmissible( ArcIndex arc) const { @@ -452,17 +447,17 @@ GenericMinCostFlow:: template bool GenericMinCostFlow::Solve() { - if (absl::GetFlag(FLAGS_min_cost_flow_check_balance) && - !CheckInputConsistency()) { + if (!CheckInputConsistency()) { return false; } - if (check_feasibility_ && !CheckFeasibility(nullptr, nullptr)) { + if (check_feasibility_ && !CheckFeasibility()) { status_ = INFEASIBLE; return false; } status_ = NOT_SOLVED; node_potential_.assign(node_potential_.size(), 0); + ResetFirstAdmissibleArcs(); if (!ScaleCosts()) return false; if (!Optimize()) return false; @@ -1170,6 +1165,7 @@ SimpleMinCostFlow::Status SimpleMinCostFlow::SolveWithPossibleAdjustment( arc_flow_[arc] = min_cost_flow.Flow(PermutedArc(arc)); } } + return min_cost_flow.status(); } diff --git a/ortools/graph/min_cost_flow.h b/ortools/graph/min_cost_flow.h index 270d2c6157..3ac55df0c8 100644 --- a/ortools/graph/min_cost_flow.h +++ b/ortools/graph/min_cost_flow.h @@ -201,9 +201,6 @@ class MinCostFlowBase { // execution. // // Some details on how to deal with this: - // - The sum of all incoming/outgoing capacity at each node should not - // overflow. TODO(user): this is not always properly checked and probably - // deserve a different return status. // - Since we scale cost, each arc unit cost times (num_nodes + 1) should // not overflow. We detect that at the beginning of the Solve(). // - This is however not sufficient as the node potential depends on the @@ -225,7 +222,28 @@ class MinCostFlowBase { // deal with since we will still have the proper flow on each arc. It is // thus possible to recompute the total cost in double or using // absl::int128 if the need arise. - BAD_COST_RANGE + BAD_COST_RANGE, + + // This is returned in our initial check if the arc capacity are too large. + // For each node these quantity should not overflow the FlowQuantity type + // which is int64_t by default. + // - Its initial excess (+supply or -demand) + sum incoming arc capacity + // - Its initial excess (+supply or -demand) - sum outgoing arc capacity. + // + // Note that these are upper bounds that guarantee that no overflow will + // take place during the algorithm execution. It is possible to go over that + // and still encounter no overflow, but since we cannot guarantee this, we + // rather check early and return a proper status. + // + // Note that we might cap the capacity of the arcs if we can detect it + // is too large in order to avoid returning this status for simple case, + // like if a client used int64_t max for all arcs out of the source. + // + // TODO(user): Not sure this is a good idea, probably better to make sure + // client use reasonable capacities. Also we should template by FlowQuantity + // and allow use of absl::int128 so we never have issue if the input is + // int64_t. + BAD_CAPACITY_RANGE, }; }; @@ -379,6 +397,7 @@ class GenericMinCostFlow : public MinCostFlowBase { public: typedef typename Graph::NodeIndex NodeIndex; typedef typename Graph::ArcIndex ArcIndex; + typedef typename Graph::IncomingArcIterator IncomingArcIterator; typedef typename Graph::OutgoingArcIterator OutgoingArcIterator; typedef typename Graph::OutgoingOrOppositeIncomingArcIterator OutgoingOrOppositeIncomingArcIterator; @@ -413,31 +432,9 @@ class GenericMinCostFlow : public MinCostFlowBase { // Sets the capacity for the given arc. void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity); - // Sets the flow for the given arc. Note that new_flow must be smaller than - // the capacity of the arc. - void SetArcFlow(ArcIndex arc, ArcFlowType new_flow); - // Solves the problem, returning true if a min-cost flow could be found. bool Solve(); - // Checks for feasibility, i.e., that all the supplies and demands can be - // matched without exceeding bottlenecks in the network. - // If infeasible_supply_node (resp. infeasible_demand_node) are not NULL, - // they are populated with the indices of the nodes where the initial supplies - // (resp. demands) are too large. Feasible values for the supplies and - // demands are accessible through FeasibleSupply. - // Note that CheckFeasibility is called by Solve() when the flag - // min_cost_flow_check_feasibility is set to true (which is the default.) - bool CheckFeasibility(std::vector* infeasible_supply_node, - std::vector* infeasible_demand_node); - - // Makes the min-cost flow problem solvable by truncating supplies and - // demands to a level acceptable by the network. There may be several ways to - // do it. In our case, the levels are computed from the result of the max-flow - // algorithm run in CheckFeasibility(). - // MakeFeasible returns false if CheckFeasibility() was not called before. - bool MakeFeasible(); - // Returns the cost of the minimum-cost flow found by the algorithm. This // works in O(num_arcs). This will only work if the last Solve() call was // successful and returned true, otherwise it will return 0. Note that the @@ -450,24 +447,18 @@ class GenericMinCostFlow : public MinCostFlowBase { FlowQuantity Flow(ArcIndex arc) const; // Returns the capacity of the given arc. + // + // Warning: If the capacity were close to ArcFlowType::max() we might have + // adjusted them in order to avoid overflow. FlowQuantity Capacity(ArcIndex arc) const; // Returns the unscaled cost for the given arc. CostValue UnitCost(ArcIndex arc) const; - // Returns the supply at a given node. Demands are modelled as negative - // supplies. + // Returns the supply at a given node. + // Demands are modelled as negative supplies. FlowQuantity Supply(NodeIndex node) const; - // Returns the initial supply at a given node. - FlowQuantity InitialSupply(NodeIndex node) const; - - // Returns the largest supply (if > 0) or largest demand in absolute value - // (if < 0) admissible at node. If the problem is not feasible, some of these - // values will be smaller (in absolute value) than the initial supplies - // and demand given as input. - FlowQuantity FeasibleSupply(NodeIndex node) const; - // Whether to check the feasibility of the problem with a max-flow, prior to // solving it. This uses about twice as much memory, but detects infeasible // problems (where the flow can't be satisfied) and makes Solve() return @@ -481,6 +472,12 @@ class GenericMinCostFlow : public MinCostFlowBase { void SetPriceScaling(bool value) { scale_prices_ = value; } private: + // Checks for feasibility, i.e., that all the supplies and demands can be + // matched without exceeding bottlenecks in the network. + // Note that CheckFeasibility is called by Solve() when SetCheckFeasibility() + // is set to true, which is the default. + bool CheckFeasibility(); + // Returns true if the given arc is admissible i.e. if its residual capacity // is strictly positive, and its reduced cost strictly negative, i.e., pushing // more flow into it will result in a reduction of the total cost. @@ -643,12 +640,6 @@ class GenericMinCostFlow : public MinCostFlowBase { // node. This is used to create the max-flow-based feasibility checker. std::vector initial_node_excess_; - // An array containing the best acceptable excesses for each of the - // nodes. These excesses are imposed by the result of the max-flow-based - // feasibility checker for the nodes with an initial supply != 0. For the - // other nodes, the excess is simply 0. - std::vector feasible_node_excess_; - // Statistics about this class. StatsGroup stats_; diff --git a/ortools/graph/min_cost_flow_test.cc b/ortools/graph/min_cost_flow_test.cc index 3391167d48..623166a3b6 100644 --- a/ortools/graph/min_cost_flow_test.cc +++ b/ortools/graph/min_cost_flow_test.cc @@ -34,6 +34,73 @@ #include "ortools/linear_solver/linear_solver.h" namespace operations_research { +namespace { + +TEST(SolveTest, CapacityTooLarge) { + using Graph = ::util::ReverseArcListGraph; + using Solver = + ::operations_research::GenericMinCostFlow; + + const int num_nodes = 6; + const int num_arcs = 10; + auto graph = std::make_unique(num_nodes, num_arcs); + auto solver = std::make_unique(graph.get()); + const std::vector tails = {1, 2, 3, 4, 5, 0, 1, 2, 3, 4}; + const std::vector heads = {0, 1, 2, 3, 4, 5, 5, 5, 5, 5}; + const std::vector capacities = { + 3184525836262886912, 3184525836262886912, 3184525836262886912, + 3184525836262886912, 3184525836262886912, 1025, + 3184525836262886914, 3184525836262886914, 3184525836262886914, + 3184525836262886914, + }; + const std::vector objectives(num_arcs, 0); + const std::vector supplies = {-3184525836262885888, 1, 1, 1, 1, + 3184525836262885884}; + for (int64_t e = 0; e < num_arcs; ++e) { + graph->AddArc(/*tail=*/tails[e], /*head=*/heads[e]); + solver->SetArcCapacity(/*arc=*/e, /*new_capacity=*/capacities[e]); + solver->SetArcUnitCost(/*arc=*/e, /*unit_cost=*/objectives[e]); + } + for (int64_t n = 0; n < num_nodes; ++n) { + solver->SetNodeSupply(/*node=*/n, /*supply=*/supplies[n]); + } + + // This one can actually be "corrected" by our simple heuristic. + EXPECT_TRUE(solver->Solve()); + EXPECT_EQ(solver->status(), Solver::OPTIMAL); +} + +TEST(SolveTest, CapacityTooLarge2) { + using Graph = ::util::ReverseArcListGraph; + using Solver = + ::operations_research::GenericMinCostFlow; + + const int num_nodes = 3; + const int num_arcs = 6; + auto graph = std::make_unique(num_nodes, num_arcs); + auto solver = std::make_unique(graph.get()); + + // We construct a double cycle so that the incoming/outgoing flow cannot be + // easily bounded. + const int64_t kint64max = std::numeric_limits::max(); + const std::vector tails = {0, 0, 1, 1, 2, 2}; + const std::vector heads = {1, 1, 2, 2, 0, 0}; + const std::vector capacities(num_arcs, kint64max - 10); + const std::vector objectives(num_arcs, 0); + const std::vector supplies = {-(kint64max - 10), kint64max - 10, 0}; + + for (int64_t e = 0; e < num_arcs; ++e) { + graph->AddArc(/*tail=*/tails[e], /*head=*/heads[e]); + solver->SetArcCapacity(/*arc=*/e, /*new_capacity=*/capacities[e]); + solver->SetArcUnitCost(/*arc=*/e, /*unit_cost=*/objectives[e]); + } + for (int64_t n = 0; n < num_nodes; ++n) { + solver->SetNodeSupply(/*node=*/n, /*supply=*/supplies[n]); + } + + EXPECT_FALSE(solver->Solve()); + EXPECT_EQ(solver->status(), Solver::BAD_CAPACITY_RANGE); +} template void GenericMinCostFlowTester( @@ -77,11 +144,6 @@ void GenericMinCostFlowTester( } } else if (expected_status == GenericMinCostFlow::INFEASIBLE) { EXPECT_FALSE(ok); - for (int node = 0; node < graph.num_nodes(); ++node) { - FlowQuantity delta = min_cost_flow.InitialSupply(node) - - min_cost_flow.FeasibleSupply(node); - EXPECT_EQ(0, delta) << "at node " << node; - } } } } @@ -169,66 +231,6 @@ TYPED_TEST(GenericMinCostFlowTest, Test3) { kExpectedFlowCost, kExpectedFlow, GenericMinCostFlow::OPTIMAL); } -TYPED_TEST(GenericMinCostFlowTest, Infeasible) { - const int kNumNodes = 6; - const int kNumArcs = 9; - const FlowQuantity kNodeSupply[kNumNodes] = {20, 0, 0, 0, 0, -20}; - const FlowQuantity kNodeInfeasibility[kNumNodes] = {4, 0, 0, 0, 0, -4}; - const NodeIndex kTail[kNumArcs] = {0, 0, 1, 1, 1, 2, 3, 4, 4}; - const NodeIndex kHead[kNumArcs] = {1, 2, 1, 3, 4, 3, 5, 3, 5}; - const CostValue kCost[kNumArcs] = {1, 6, 3, 5, 7, 3, 1, 6, 9}; - const FlowQuantity kCapacity[kNumArcs] = {10, 10, 10, 6, 6, 6, 10, 10, 10}; - const NodeIndex kInfeasibleSupplyNode[] = {0}; - const NodeIndex kInfeasibleDemandNode[] = {5}; - const NodeIndex kFeasibleSupply[] = {16}; - const NodeIndex kFeasibleDemand[] = {-16}; - - TypeParam graph(kNumNodes, kNumArcs); - for (ArcIndex arc = 0; arc < kNumArcs; ++arc) { - graph.AddArc(kTail[arc], kHead[arc]); - } - std::vector permutation; - Graphs::Build(&graph, &permutation); - EXPECT_TRUE(permutation.empty()); - - GenericMinCostFlow min_cost_flow(&graph); - for (ArcIndex arc = 0; arc < kNumArcs; ++arc) { - min_cost_flow.SetArcUnitCost(arc, kCost[arc]); - min_cost_flow.SetArcCapacity(arc, kCapacity[arc]); - } - for (ArcIndex arc = 0; arc < kNumNodes; ++arc) { - min_cost_flow.SetNodeSupply(arc, kNodeSupply[arc]); - } - std::vector infeasible_supply_node; - std::vector infeasible_demand_node; - bool feasible = min_cost_flow.CheckFeasibility(&infeasible_supply_node, - &infeasible_demand_node); - EXPECT_FALSE(feasible); - for (int i = 0; i < infeasible_supply_node.size(); ++i) { - const NodeIndex node = infeasible_supply_node[i]; - EXPECT_EQ(node, kInfeasibleSupplyNode[i]); - EXPECT_EQ(min_cost_flow.FeasibleSupply(node), kFeasibleSupply[i]); - } - for (int i = 0; i < infeasible_demand_node.size(); ++i) { - const NodeIndex node = infeasible_demand_node[i]; - EXPECT_EQ(node, kInfeasibleDemandNode[i]); - EXPECT_EQ(min_cost_flow.FeasibleSupply(node), kFeasibleDemand[i]); - } - bool ok = min_cost_flow.Solve(); - EXPECT_FALSE(ok); - EXPECT_EQ(GenericMinCostFlow::INFEASIBLE, min_cost_flow.status()); - for (NodeIndex node = 0; node < kNumNodes; ++node) { - FlowQuantity delta = - min_cost_flow.InitialSupply(node) - min_cost_flow.FeasibleSupply(node); - EXPECT_EQ(kNodeInfeasibility[node], delta); - } - EXPECT_EQ(min_cost_flow.GetOptimalCost(), 0); - min_cost_flow.MakeFeasible(); - ok = min_cost_flow.Solve(); - EXPECT_TRUE(ok); - EXPECT_EQ(GenericMinCostFlow::OPTIMAL, min_cost_flow.status()); -} - // Test on a 4x4 matrix. Example taken from // http://www.ee.oulu.fi/~mpa/matreng/eem1_2-1.htm TYPED_TEST(GenericMinCostFlowTest, Small4x4Matrix) { @@ -275,7 +277,7 @@ TYPED_TEST(GenericMinCostFlowTest, Small4x4Matrix) { TYPED_TEST(GenericMinCostFlowTest, TotalFlowCostOverflow) { const int kNumNodes = 2; const int kNumArcs = 1; - const FlowQuantity kNodeSupply[kNumNodes] = {1LL << 61, -1LL << 61}; + const FlowQuantity kNodeSupply[kNumNodes] = {1LL << 61, -(1LL << 61)}; const NodeIndex kTail[kNumArcs] = {0}; const NodeIndex kHead[kNumArcs] = {1}; const CostValue kCost[kNumArcs] = {10}; @@ -298,7 +300,7 @@ TEST(GenericMinCostFlowTest, OverflowPrevention1) { mcf.SetNodeSupply(1, -std::numeric_limits::max()); EXPECT_FALSE(mcf.Solve()); - EXPECT_EQ(mcf.status(), MinCostFlowBase::BAD_COST_RANGE); + EXPECT_EQ(mcf.status(), MinCostFlowBase::BAD_CAPACITY_RANGE); } TEST(GenericMinCostFlowTest, OverflowPrevention2) { @@ -976,4 +978,5 @@ static void BM_MinCostFlowOnMultiMatchingProblem(benchmark::State& state) { BENCHMARK(BM_MinCostFlowOnMultiMatchingProblem); +} // namespace } // namespace operations_research diff --git a/ortools/graph/minimum_vertex_cover.cc b/ortools/graph/minimum_vertex_cover.cc new file mode 100644 index 0000000000..c74d653474 --- /dev/null +++ b/ortools/graph/minimum_vertex_cover.cc @@ -0,0 +1,103 @@ +// Copyright 2010-2024 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/graph/minimum_vertex_cover.h" + +#include + +#include "absl/log/check.h" +#include "ortools/graph/ebert_graph.h" +#include "ortools/graph/max_flow.h" + +namespace operations_research { + +std::vector BipartiteMinimumVertexCover( + const std::vector>& left_to_right_arcs, int num_right) { + // This algorithm first uses the maximum flow to find a maximum matching. Then + // it uses the same method outlined in the proof of Konig's theorem to + // transform the maximum matching into a minimum vertex cover. + // + // More concretely, it uses a DFS starting with unmatched nodes and + // alternating matched/unmatched edges to find a minimum vertex cover. + SimpleMaxFlow max_flow; + const int num_left = left_to_right_arcs.size(); + std::vector arcs; + for (int i = 0; i < num_left; ++i) { + for (const int right_node : left_to_right_arcs[i]) { + DCHECK_GE(right_node, num_left); + DCHECK_LT(right_node, num_right + num_left); + arcs.push_back(max_flow.AddArcWithCapacity(i, right_node, 1)); + } + } + std::vector> adj_list = left_to_right_arcs; + adj_list.resize(num_left + num_right); + for (int i = 0; i < num_left; ++i) { + for (const int right_node : left_to_right_arcs[i]) { + adj_list[right_node].push_back(i); + } + } + const int sink = num_left + num_right; + const int source = num_left + num_right + 1; + for (int i = 0; i < num_left; ++i) { + max_flow.AddArcWithCapacity(source, i, 1); + } + for (int i = 0; i < num_right; ++i) { + max_flow.AddArcWithCapacity(i + num_left, sink, 1); + } + CHECK(max_flow.Solve(source, sink) == SimpleMaxFlow::OPTIMAL); + std::vector maximum_matching(num_left + num_right, -1); + for (const ArcIndex arc : arcs) { + if (max_flow.Flow(arc) > 0) { + maximum_matching[max_flow.Tail(arc)] = max_flow.Head(arc); + maximum_matching[max_flow.Head(arc)] = max_flow.Tail(arc); + } + } + // We do a DFS starting with unmatched nodes and alternating matched/unmatched + // edges. + std::vector in_alternating_path(num_left + num_right, false); + std::vector to_visit; + for (int i = 0; i < num_left; ++i) { + if (maximum_matching[i] == -1) { + to_visit.push_back(i); + } + } + while (!to_visit.empty()) { + const int current = to_visit.back(); + to_visit.pop_back(); + if (in_alternating_path[current]) { + continue; + } + in_alternating_path[current] = true; + for (const int j : adj_list[current]) { + if (current < num_left && maximum_matching[current] != j) { + to_visit.push_back(j); + } else if (current >= num_left && maximum_matching[current] == j) { + to_visit.push_back(j); + } + } + } + std::vector minimum_vertex_cover(num_left + num_right, false); + for (int i = 0; i < num_left; ++i) { + if (!in_alternating_path[i]) { + minimum_vertex_cover[i] = true; + } + } + for (int i = num_left; i < num_left + num_right; ++i) { + if (in_alternating_path[i]) { + minimum_vertex_cover[i] = true; + } + } + return minimum_vertex_cover; +} + +} // namespace operations_research diff --git a/ortools/graph/minimum_vertex_cover.h b/ortools/graph/minimum_vertex_cover.h new file mode 100644 index 0000000000..13b2c4f56b --- /dev/null +++ b/ortools/graph/minimum_vertex_cover.h @@ -0,0 +1,35 @@ +// Copyright 2010-2024 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_GRAPH_MINIMUM_VERTEX_COVER_H_ +#define OR_TOOLS_GRAPH_MINIMUM_VERTEX_COVER_H_ + +#include + +namespace operations_research { + +// This method computes a minimum vertex cover for the bipartite graph. +// +// If we define num_left=left_to_right_arcs.size(), the "left" nodes are +// integers in [0, num_left), and the "right" nodes are integers in [num_left, +// num_left + num_right). +// +// Returns a vector of size num_left+num_right, such that element #l is true if +// it is part of the minimum vertex cover and false if it is part of the maximum +// independent set (one is the complement of the other). +std::vector BipartiteMinimumVertexCover( + const std::vector>& left_to_right_arcs, int num_right); + +} // namespace operations_research + +#endif // OR_TOOLS_GRAPH_MINIMUM_VERTEX_COVER_H_ diff --git a/ortools/graph/minimum_vertex_cover_test.cc b/ortools/graph/minimum_vertex_cover_test.cc new file mode 100644 index 0000000000..6c63ae091f --- /dev/null +++ b/ortools/graph/minimum_vertex_cover_test.cc @@ -0,0 +1,78 @@ +// Copyright 2010-2024 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/graph/minimum_vertex_cover.h" + +#include + +#include "absl/algorithm/container.h" +#include "gtest/gtest.h" + +namespace operations_research { +namespace { + +// Creates the complete bipartite graph K(n, m). +std::vector> MakeCompleteBipartiteGraph(int num_left, + int num_right) { + std::vector adjacencies(num_right); + absl::c_iota(adjacencies, num_left); + return std::vector>(num_left, adjacencies); +} + +TEST(BipartiteMinimumVertexCoverTest, BasicBehavior) { + const int num_right = 4; + const std::vector> left_to_right = { + {5}, {4, 5, 6}, {5}, {5, 6, 7}}; + const auto cover = BipartiteMinimumVertexCover(left_to_right, num_right); + EXPECT_EQ(absl::c_count(cover, true), 3); + EXPECT_EQ(absl::c_count(cover, false), 5); +} + +TEST(BipartiteMinimumVertexCoverTest, StarGraph) { + const std::vector> left_to_right = + MakeCompleteBipartiteGraph(1, 4); + const auto cover = BipartiteMinimumVertexCover(left_to_right, 4); + EXPECT_EQ(absl::c_count(cover, true), 1); + EXPECT_EQ(absl::c_count(cover, false), 4); +} + +TEST(BipartiteMinimumVertexCoverTest, UtilityGraph) { + const std::vector> left_to_right = + MakeCompleteBipartiteGraph(3, 3); + const auto cover = BipartiteMinimumVertexCover(left_to_right, 3); + EXPECT_EQ(absl::c_count(cover, true), 3); + EXPECT_EQ(absl::c_count(cover, false), 3); +} + +TEST(BipartiteMinimumVertexCoverTest, DuplicateEdges) { + const int num_right = 4; + const std::vector> left_to_right = { + {5, 5}, {4, 4, 5, 6}, {5, 5, 5}, {5, 5, 5, 6, 6, 7}}; + EXPECT_EQ(absl::c_count(BipartiteMinimumVertexCover(left_to_right, num_right), + true), + 3); + EXPECT_EQ(absl::c_count(BipartiteMinimumVertexCover(left_to_right, num_right), + false), + 5); +} + +TEST(BipartiteMinimumVertexCoverTest, Empty) { + const int num_right = 4; + const std::vector> left_to_right = {{}, {}}; + EXPECT_EQ(absl::c_count(BipartiteMinimumVertexCover(left_to_right, num_right), + false), + 6); +} + +} // namespace +} // namespace operations_research diff --git a/ortools/graph/python/min_cost_flow.cc b/ortools/graph/python/min_cost_flow.cc index 98f62f511e..7d0f11028d 100644 --- a/ortools/graph/python/min_cost_flow.cc +++ b/ortools/graph/python/min_cost_flow.cc @@ -51,6 +51,7 @@ PYBIND11_MODULE(min_cost_flow, m) { pybind11::enum_(smcf, "Status") .value("BAD_COST_RANGE", MinCostFlowBase::Status::BAD_COST_RANGE) + .value("BAD_CAPACITY_RANGE", MinCostFlowBase::Status::BAD_CAPACITY_RANGE) .value("BAD_RESULT", MinCostFlowBase::Status::BAD_RESULT) .value("FEASIBLE", MinCostFlowBase::Status::FEASIBLE) .value("INFEASIBLE", MinCostFlowBase::Status::INFEASIBLE)