From d77467d11ab9f08893171cac71a4ed7f0e0645fb Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Sat, 28 Dec 2024 11:22:06 +0100 Subject: [PATCH] revamp graph code; remove deprecated graph implementations --- examples/cpp/dimacs_assignment.cc | 58 ++--- ortools/graph/BUILD.bazel | 17 +- ortools/graph/assignment.cc | 13 +- ortools/graph/assignment.h | 15 +- ortools/graph/assignment_test.cc | 40 ++- ortools/graph/ebert_graph.h | 60 ++++- ortools/graph/ebert_graph_test.cc | 6 - ortools/graph/generic_max_flow.h | 224 +++++++++-------- ortools/graph/generic_max_flow_test.cc | 231 +++++++++++------- ortools/graph/graph.h | 33 ++- ortools/graph/graph_io.h | 12 + ortools/graph/linear_assignment.h | 118 ++++----- ortools/graph/linear_assignment_test.cc | 32 +-- ortools/graph/min_cost_flow.h | 12 +- ortools/graph/min_cost_flow_test.cc | 39 +-- ortools/graph/python/linear_sum_assignment.cc | 2 +- .../assignment_linear_sum_assignment.cc | 2 + 17 files changed, 551 insertions(+), 363 deletions(-) diff --git a/examples/cpp/dimacs_assignment.cc b/examples/cpp/dimacs_assignment.cc index 4a3f37e9b1..165071a1a0 100644 --- a/examples/cpp/dimacs_assignment.cc +++ b/examples/cpp/dimacs_assignment.cc @@ -13,20 +13,21 @@ #include #include +#include #include #include #include #include "absl/container/flat_hash_map.h" +#include "absl/flags/flag.h" #include "absl/strings/str_format.h" #include "examples/cpp/parse_dimacs_assignment.h" #include "examples/cpp/print_dimacs_assignment.h" #include "ortools/algorithms/hungarian.h" -#include "ortools/base/commandlineflags.h" #include "ortools/base/init_google.h" #include "ortools/base/logging.h" #include "ortools/base/timer.h" -#include "ortools/graph/ebert_graph.h" +#include "ortools/graph/graph.h" #include "ortools/graph/linear_assignment.h" ABSL_FLAG(bool, assignment_compare_hungarian, false, @@ -35,16 +36,18 @@ ABSL_FLAG(std::string, assignment_problem_output_file, "", "Print the problem to this file in DIMACS format (after layout " "is optimized, if applicable)."); ABSL_FLAG(bool, assignment_reverse_arcs, false, - "Ignored if --assignment_static_graph=true. Use StarGraph " - "if true, ForwardStarGraph if false."); + "Ignored if --assignment_static_graph=true. Use ReverseArcListGraph " + "if true, ListGraph if false."); ABSL_FLAG(bool, assignment_static_graph, true, - "Use the ForwardStarStaticGraph representation, " - "otherwise ForwardStarGraph or StarGraph according " + "Use the StaticGraph representation, " + "otherwise ListGraph or ReverseArcListGraph according " "to --assignment_reverse_arcs."); namespace operations_research { -typedef ForwardStarStaticGraph GraphType; +using NodeIndex = int32_t; +using ArcIndex = int32_t; +using CostValue = int64_t; template CostValue BuildAndSolveHungarianInstance( @@ -57,9 +60,7 @@ CostValue BuildAndSolveHungarianInstance( // First we have to find the biggest cost magnitude so we can // initialize the arc costs that aren't really there. CostValue largest_cost_magnitude = 0; - for (typename GraphType::ArcIterator arc_it(graph); arc_it.Ok(); - arc_it.Next()) { - ArcIndex arc = arc_it.Index(); + for (const auto arc : graph.AllForwardArcs()) { CostValue cost_magnitude = std::abs(assignment.ArcCost(arc)); largest_cost_magnitude = std::max(largest_cost_magnitude, cost_magnitude); } @@ -77,15 +78,9 @@ CostValue BuildAndSolveHungarianInstance( // hungarian algorithm). We opt for the alternative of iterating // over hte arcs via adjacency lists, which gives us the arc tails // implicitly. - for (typename GraphType::NodeIterator node_it(graph); node_it.Ok(); - node_it.Next()) { - NodeIndex node = node_it.Index(); - NodeIndex tail = (node - GraphType::kFirstNode); - for (typename GraphType::OutgoingArcIterator arc_it(graph, node); - arc_it.Ok(); arc_it.Next()) { - ArcIndex arc = arc_it.Index(); - NodeIndex head = - (graph.Head(arc) - assignment.NumLeftNodes() - GraphType::kFirstNode); + for (const auto tail : graph.AllNodes()) { + for (const auto arc : graph.OutgoingArcs(tail)) { + NodeIndex head = graph.Head(arc) - assignment.NumLeftNodes(); double cost = static_cast(assignment.ArcCost(arc)); hungarian_cost[tail][head] = cost; } @@ -132,17 +127,8 @@ int SolveDimacsAssignment(int argc, char* argv[]) { LOG(FATAL) << error_message; } if (!absl::GetFlag(FLAGS_assignment_problem_output_file).empty()) { - // The following tail array management stuff is done in a generic - // way so we can plug in different types of graphs for which the - // TailArrayManager template can be instantiated, even though we - // know the type of the graph explicitly. In this way, the type of - // the graph can be switched just by changing the graph type in - // this file and making no other changes to the code. - TailArrayManager tail_array_manager(graph); PrintDimacsAssignmentProblem( - *assignment, tail_array_manager, - absl::GetFlag(FLAGS_assignment_problem_output_file)); - tail_array_manager.ReleaseTailArrayIfForwardGraph(); + *assignment, absl::GetFlag(FLAGS_assignment_problem_output_file)); } CostValue hungarian_cost = 0.0; bool hungarian_solved = false; @@ -175,10 +161,9 @@ int SolveDimacsAssignment(int argc, char* argv[]) { static const char* const kUsageTemplate = "usage: %s "; -using ::operations_research::ForwardStarGraph; -using ::operations_research::ForwardStarStaticGraph; +using ::operations_research::ArcIndex; +using ::operations_research::NodeIndex; using ::operations_research::SolveDimacsAssignment; -using ::operations_research::StarGraph; int main(int argc, char* argv[]) { std::string usage; @@ -194,10 +179,13 @@ int main(int argc, char* argv[]) { } if (absl::GetFlag(FLAGS_assignment_static_graph)) { - return SolveDimacsAssignment(argc, argv); + return SolveDimacsAssignment<::util::StaticGraph>( + argc, argv); } else if (absl::GetFlag(FLAGS_assignment_reverse_arcs)) { - return SolveDimacsAssignment(argc, argv); + return SolveDimacsAssignment< + ::util::ReverseArcListGraph>(argc, argv); } else { - return SolveDimacsAssignment(argc, argv); + return SolveDimacsAssignment<::util::ListGraph>(argc, + argv); } } diff --git a/ortools/graph/BUILD.bazel b/ortools/graph/BUILD.bazel index c06aa21bcc..ee0f445798 100644 --- a/ortools/graph/BUILD.bazel +++ b/ortools/graph/BUILD.bazel @@ -351,12 +351,12 @@ cc_library( name = "ebert_graph", hdrs = ["ebert_graph.h"], deps = [ + ":iterators", "//ortools/base", - "//ortools/base:types", "//ortools/util:permutation", "//ortools/util:zvector", + "@com_google_absl//absl/base:core_headers", "@com_google_absl//absl/strings", - "@com_google_googletest//:gtest_prod", ], ) @@ -399,7 +399,7 @@ cc_library( cc_test( name = "shortest_paths_test", - timeout = "eternal", + size = "enormous", srcs = ["shortest_paths_test.cc"], tags = ["noasan"], # Times out occasionally in ASAN mode. deps = [ @@ -510,6 +510,7 @@ cc_test( ":generic_max_flow", ":graph", ":graphs", + "@com_google_absl//absl/random:bit_gen_ref", "//ortools/base", "//ortools/base:gmock_main", "//ortools/linear_solver", @@ -533,6 +534,7 @@ cc_library( }), deps = [ ":ebert_graph", + ":generic_max_flow", ":graph", ":graphs", ":max_flow", @@ -600,9 +602,8 @@ cc_library( srcs = ["assignment.cc"], hdrs = ["assignment.h"], deps = [ - ":ebert_graph", + ":graph", ":linear_assignment", - "@com_google_absl//absl/flags:flag", ], ) @@ -612,8 +613,10 @@ cc_test( srcs = ["assignment_test.cc"], deps = [ ":assignment", - ":ebert_graph", "//ortools/base:gmock_main", + "@com_google_absl//absl/log:check", + "@com_google_absl//absl/random:distributions", + "@com_google_benchmark//:benchmark", ], ) @@ -626,12 +629,10 @@ cc_library( deps = [ ":ebert_graph", "//ortools/base", - "//ortools/base:types", "//ortools/util:permutation", "//ortools/util:zvector", "@com_google_absl//absl/flags:flag", "@com_google_absl//absl/strings:str_format", - "@com_google_googletest//:gtest_prod", ], ) diff --git a/ortools/graph/assignment.cc b/ortools/graph/assignment.cc index f685a882ae..7ec3e97f52 100644 --- a/ortools/graph/assignment.cc +++ b/ortools/graph/assignment.cc @@ -14,14 +14,18 @@ #include "ortools/graph/assignment.h" #include +#include #include -#include "absl/flags/flag.h" -#include "ortools/graph/ebert_graph.h" +#include "ortools/graph/graph.h" #include "ortools/graph/linear_assignment.h" namespace operations_research { +using ArcIndex = int32_t; +using NodeIndex = int32_t; +using CostValue = int64_t; + SimpleLinearSumAssignment::SimpleLinearSumAssignment() : num_nodes_(0) {} ArcIndex SimpleLinearSumAssignment::AddArcWithCost(NodeIndex left_node, @@ -66,8 +70,9 @@ SimpleLinearSumAssignment::Status SimpleLinearSumAssignment::Solve() { } const ArcIndex num_arcs = arc_cost_.size(); - ForwardStarGraph graph(2 * num_nodes_, num_arcs); - LinearSumAssignment assignment(graph, num_nodes_); + ::util::ListGraph<> graph(2 * num_nodes_, num_arcs); + LinearSumAssignment<::util::ListGraph<>, CostValue> assignment(graph, + num_nodes_); for (ArcIndex arc = 0; arc < num_arcs; ++arc) { graph.AddArc(arc_tail_[arc], num_nodes_ + arc_head_[arc]); assignment.SetArcCost(arc, arc_cost_[arc]); diff --git a/ortools/graph/assignment.h b/ortools/graph/assignment.h index 34599949b2..5ca0b90bbd 100644 --- a/ortools/graph/assignment.h +++ b/ortools/graph/assignment.h @@ -48,14 +48,17 @@ #ifndef OR_TOOLS_GRAPH_ASSIGNMENT_H_ #define OR_TOOLS_GRAPH_ASSIGNMENT_H_ +#include #include -#include "ortools/graph/ebert_graph.h" - namespace operations_research { class SimpleLinearSumAssignment { public: + typedef int32_t NodeIndex; + typedef int32_t ArcIndex; + typedef int64_t CostValue; + // The constructor takes no size. // New node indices will be created lazily by AddArcWithCost(). SimpleLinearSumAssignment(); @@ -67,6 +70,14 @@ class SimpleLinearSumAssignment { delete; #endif + // Reserves space for the given number of arcs, to avoid reallocation in + // `AddArcWithCost(). + void ReserveArcs(ArcIndex num_arcs) { + arc_tail_.reserve(num_arcs); + arc_head_.reserve(num_arcs); + arc_cost_.reserve(num_arcs); + } + // Adds an arc from a left node to a right node with a given cost. // * Node indices must be non-negative (>= 0). For a perfect // matching to exist on n nodes, the values taken by "left_node" diff --git a/ortools/graph/assignment_test.cc b/ortools/graph/assignment_test.cc index 4c169f9d45..6675c45db8 100644 --- a/ortools/graph/assignment_test.cc +++ b/ortools/graph/assignment_test.cc @@ -14,9 +14,14 @@ #include "ortools/graph/assignment.h" #include +#include +#include +#include +#include "absl/log/check.h" +#include "absl/random/distributions.h" +#include "benchmark/benchmark.h" #include "gtest/gtest.h" -#include "ortools/graph/ebert_graph.h" namespace operations_research { @@ -69,6 +74,7 @@ TEST(SimpleLinearSumAssignmentTest, InfeasibleProblem) { } TEST(SimpleLinearSumAssignmentTest, Overflow) { + using CostValue = SimpleLinearSumAssignment::CostValue; SimpleLinearSumAssignment assignment; assignment.AddArcWithCost(0, 0, std::numeric_limits::max()); assignment.AddArcWithCost(0, 1, std::numeric_limits::max()); @@ -78,4 +84,36 @@ TEST(SimpleLinearSumAssignmentTest, Overflow) { EXPECT_EQ(0, assignment.OptimalCost()); } +void BM_SimpleLinearSumAssignment(benchmark::State& state) { + using CostValue = SimpleLinearSumAssignment::CostValue; + constexpr CostValue kCostLimit = 1000000; + + const int num_left = state.range(0); + const int degree = state.range(1); + + std::mt19937 rng(12345); + + // Each left node is connected to `degree` right nodes. + std::vector> arcs; + arcs.reserve(num_left * degree); + for (int left = 0; left < num_left; ++left) { + for (int i = 0; i < degree; ++i) { + const int right = (left + i) % num_left; + const CostValue cost = + absl::Uniform(rng, 0, absl::Uniform(rng, 0, kCostLimit)); + arcs.emplace_back(left, right, cost); + } + } + + for (auto _ : state) { + SimpleLinearSumAssignment assignment; + assignment.ReserveArcs(num_left * degree); + for (const auto& [left, right, cost] : arcs) { + assignment.AddArcWithCost(left, right, cost); + } + CHECK_EQ(assignment.Solve(), SimpleLinearSumAssignment::OPTIMAL); + } +} +BENCHMARK(BM_SimpleLinearSumAssignment)->ArgPair(100, 2)->ArgPair(100, 20); + } // namespace operations_research diff --git a/ortools/graph/ebert_graph.h b/ortools/graph/ebert_graph.h index 812b876718..2589ec0fcb 100644 --- a/ortools/graph/ebert_graph.h +++ b/ortools/graph/ebert_graph.h @@ -168,15 +168,19 @@ // the arc tail array adds another m * sizeof(NodeIndexType). #include +#include #include +#include #include #include #include #include #include +#include "absl/base/attributes.h" #include "absl/strings/str_cat.h" #include "ortools/base/logging.h" +#include "ortools/graph/iterators.h" #include "ortools/util/permutation.h" #include "ortools/util/zvector.h" @@ -202,12 +206,28 @@ typedef int64_t FlowQuantity; typedef int64_t CostValue; typedef EbertGraph StarGraph; typedef ForwardEbertGraph ForwardStarGraph; -typedef ForwardStaticGraph ForwardStarStaticGraph; typedef ZVector NodeIndexArray; typedef ZVector ArcIndexArray; typedef ZVector QuantityArray; typedef ZVector CostArray; +// Adapt our old iteration style to support range-based for loops. Add typedefs +// required by std::iterator_traits. +#define DEFINE_STL_ITERATOR_FUNCTIONS(iterator_class_name) \ + using iterator_category = std::input_iterator_tag; \ + using difference_type = ptrdiff_t; \ + using pointer = const ArcIndexType*; \ + using value_type = ArcIndexType; \ + using reference = value_type; \ + bool operator!=(const iterator_class_name& other) const { \ + return this->arc_ != other.arc_; \ + } \ + bool operator==(const iterator_class_name& other) const { \ + return this->arc_ == other.arc_; \ + } \ + ArcIndexType operator*() const { return this->arc_; } \ + void operator++() { this->Next(); } + template class StarGraphBase { public: @@ -400,6 +420,8 @@ class StarGraphBase { // Returns the index of the arc currently pointed to by the iterator. ArcIndexType Index() const { return arc_; } + DEFINE_STL_ITERATOR_FUNCTIONS(OutgoingArcIterator); + private: // Returns true if the invariant for the iterator is verified. // To be used in a DCHECK. @@ -530,7 +552,7 @@ class PermutationIndexComparisonByArcHead { }; template -class ForwardStaticGraph +class ABSL_DEPRECATED("Use `::util::StaticGraph<>` instead.") ForwardStaticGraph : public StarGraphBase > { typedef StarGraphBase, // and other type shortcuts; see the bottom of this file. template -class EbertGraph +class ABSL_DEPRECATED("Use `::util::ListGraph<>` instead.") EbertGraph : public EbertGraphBase > { typedef EbertGraphBase* permutation) { permutation->clear(); } + ArcIndexType OppositeArc(ArcIndex arc) const { return Opposite(arc); } + util::BeginEndWrapper OutgoingArcs( + NodeIndexType node) const { + return util::BeginEndWrapper( + typename Base::OutgoingArcIterator(*this, node), + typename Base::OutgoingArcIterator(*this, node, kNilArc)); + } + util::BeginEndWrapper + OutgoingOrOppositeIncomingArcs(NodeIndexType node) const { + return util::BeginEndWrapper( + OutgoingOrOppositeIncomingArcIterator(*this, node), + OutgoingOrOppositeIncomingArcIterator(*this, node, kNilArc)); + } + util::BeginEndWrapper + OutgoingOrOppositeIncomingArcsStartingFrom(NodeIndexType node, + ArcIndex arc) const { + return util::BeginEndWrapper( + OutgoingOrOppositeIncomingArcIterator(*this, node, arc), + OutgoingOrOppositeIncomingArcIterator(*this, node, kNilArc)); + } + // Utility function to check that an arc index is within the bounds. // It is exported so that users of the EbertGraph class can use it. // To be used in a DCHECK. @@ -1566,7 +1616,7 @@ class EbertGraph // A forward-star-only graph representation for greater efficiency in // those algorithms that don't need reverse arcs. template -class ForwardEbertGraph +class ABSL_DEPRECATED("Use `::util::ListGraph<>` instead.") ForwardEbertGraph : public EbertGraphBase > { typedef EbertGraphBase static void BM_RandomAnnotatedArcs(benchmark::State& state) { @@ -1229,11 +1227,9 @@ static void BM_RandomAnnotatedArcs(benchmark::State& state) { BENCHMARK_TEMPLATE2(BM_RandomAnnotatedArcs, StarGraph, false); BENCHMARK_TEMPLATE2(BM_RandomAnnotatedArcs, ForwardStarGraph, false); -BENCHMARK_TEMPLATE2(BM_RandomAnnotatedArcs, ForwardStarStaticGraph, false); BENCHMARK_TEMPLATE2(BM_RandomAnnotatedArcs, StarGraph, true); BENCHMARK_TEMPLATE2(BM_RandomAnnotatedArcs, ForwardStarGraph, true); -BENCHMARK_TEMPLATE2(BM_RandomAnnotatedArcs, ForwardStarStaticGraph, true); template static void BM_AddRandomArcsAndDoNotRetrieveGraph(benchmark::State& state) { @@ -1256,7 +1252,5 @@ static void BM_AddRandomArcsAndDoNotRetrieveGraph(benchmark::State& state) { BENCHMARK_TEMPLATE(BM_AddRandomArcsAndDoNotRetrieveGraph, StarGraph); BENCHMARK_TEMPLATE(BM_AddRandomArcsAndDoNotRetrieveGraph, ForwardStarGraph); -BENCHMARK_TEMPLATE(BM_AddRandomArcsAndDoNotRetrieveGraph, - ForwardStarStaticGraph); } // namespace operations_research diff --git a/ortools/graph/generic_max_flow.h b/ortools/graph/generic_max_flow.h index 194226151c..c1ad4bba79 100644 --- a/ortools/graph/generic_max_flow.h +++ b/ortools/graph/generic_max_flow.h @@ -211,18 +211,33 @@ class MaxFlowStatusClass { }; }; -// 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. +// Generic MaxFlow that works on graphs with the notion of "reverse" arcs. +// That is, the graph is directed, and each arc 'tail -> head' must be +// associated to an unique reverse arc going in the opposite direction +// 'head -> tail'. We must also have reverse[reverse[arc]] = arc. +// +// This works with all the reverse arc graphs from 'util/graph/graph.h' and +// uses the API defined there. +// +// We actually support two kind of graphs with "reverse" arcs depending on the +// value of Graph::kHasNegativeReverseArcs: +// - If it is true, the arcs from the "user graph" are called "direct" arcs and +// are assumed to be indexed in [0, num_arcs). These are the only ones that +// can have a positive capacity. All the "reverse/opposite" of these arcs +// will always have negative indices in [-num_arcs, 0) and a capacity of +// zero. +// - If it is false, then we only have "direct" arcs in [0, num_arcs), the +// reverse of each arc must lives in the same space, and both an arc and its +// reverse can have a positive initial capacity. This can lead to twice fewer +// arcs and a faster algo if the "user graph" as a lot of reverse arcs +// already. 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; + // TODO(b/385094969): This should be a template parameter. + typedef FlowQuantity FlowQuantityT; // The height of a node never excess 2 times the number of node, so we // use the same type as a Node index. @@ -257,6 +272,9 @@ class GenericMaxFlow : public MaxFlowStatusClass { NodeIndex GetSinkNodeIndex() const { return sink_; } // Sets the capacity for arc to new_capacity. + // + // Note that this will be ignored for self-arc, so do not be surprised if you + // get zero when reading the capacity of a self-arc back. void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity); // Returns true if a maximum flow was solved. @@ -265,25 +283,27 @@ class GenericMaxFlow : public MaxFlowStatusClass { // 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_. + // Returns the current flow on the given arc. + // + // Note that on a couple (arc, opposite_arc) the flow only goes in one + // direction (where it is positive) and the other direction will have the + // negation of that flow. FlowQuantity Flow(ArcIndex arc) const { - if (IsArcDirect(arc)) { - return residual_arc_capacity_[Opposite(arc)]; - } else { + if constexpr (Graph::kHasNegativeReverseArcs) { + if (IsArcDirect(arc)) return residual_arc_capacity_[Opposite(arc)]; return -residual_arc_capacity_[arc]; } + return initial_capacity_[arc] - residual_arc_capacity_[arc]; } - // Returns the capacity of arc using the equations given in the comment on - // residual_arc_capacity_. + // Returns the initial capacity of an arc. FlowQuantity Capacity(ArcIndex arc) const { - if (IsArcDirect(arc)) { + if constexpr (Graph::kHasNegativeReverseArcs) { + if (!IsArcDirect(arc)) return 0; return residual_arc_capacity_[arc] + residual_arc_capacity_[Opposite(arc)]; - } else { - return 0; } + return initial_capacity_[arc]; } // Returns the nodes reachable from the source in the residual graph, the @@ -328,12 +348,6 @@ class GenericMaxFlow : public MaxFlowStatusClass { 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. @@ -463,6 +477,11 @@ class GenericMaxFlow : public MaxFlowStatusClass { // reduces the amount of memory for this information by a factor 2. ZVector residual_arc_capacity_; + // The initial capacity as set by SetArcCapacity(), unused if + // `Graph::kNegativeReverseArcs`, as we can always recover the initial + // capacity from the residual capacities of an arc and its reverse. + std::unique_ptr initial_capacity_; + // An array representing the first admissible arc for each node in graph_. std::unique_ptr first_admissible_arc_; @@ -569,17 +588,24 @@ GenericMaxFlow::GenericMaxFlow(const Graph* graph, NodeIndex source, DCHECK(graph->IsNodeValid(sink)); const NodeIndex max_num_nodes = Graphs::NodeReservation(*graph_); if (max_num_nodes > 0) { + // We will initialize them in InitializePreflow(), so no need for memset. + // + // TODO(user): Unfortunately std/absl::make_unique_for_overwrite is not + // yet available in open-source. node_excess_ = std::make_unique(max_num_nodes); node_potential_ = std::make_unique(max_num_nodes); - first_admissible_arc_ = - absl::make_unique(max_num_nodes); - std::fill(first_admissible_arc_.get(), - first_admissible_arc_.get() + max_num_nodes, Graph::kNilArc); + first_admissible_arc_ = std::make_unique(max_num_nodes); 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); + if constexpr (Graph::kHasNegativeReverseArcs) { + residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1); + } else { + // We will need to store the initial capacity in this case. + initial_capacity_ = std::make_unique(max_num_arcs); + residual_arc_capacity_.Reserve(0, max_num_arcs - 1); + } residual_arc_capacity_.SetAll(0); } } @@ -590,31 +616,20 @@ void GenericMaxFlow::SetArcCapacity(ArcIndex arc, 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. - } + + // This serves no purpose from a max-flow point of view, so it is safer to + // just leave the capacity of all self-arc at zero. + if (Head(arc) == Tail(arc)) return; + 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]); + residual_arc_capacity_[arc] = new_capacity; + + // Since the class might have already be used, we need to make sure we clear + // any previous state on this arc and its reverse. + if constexpr (Graph::kHasNegativeReverseArcs) { + residual_arc_capacity_[Opposite(arc)] = 0; } 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); + initial_capacity_[arc] = new_capacity; } } @@ -690,9 +705,7 @@ bool GenericMaxFlow::AugmentingPathExists() const { 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(); + for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) { if (residual_arc_capacity_[arc] > 0) { const NodeIndex head = graph_->Head(arc); if (!is_reached[head]) { @@ -708,9 +721,7 @@ bool GenericMaxFlow::AugmentingPathExists() const { template bool GenericMaxFlow::CheckRelabelPrecondition(NodeIndex node) const { DCHECK(IsActive(node)); - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok(); - it.Next()) { - const ArcIndex arc = it.Index(); + for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) { DCHECK(!IsAdmissible(node, arc, node_potential_.data())) << DebugString("CheckRelabelPrecondition:", arc); } @@ -724,13 +735,13 @@ std::string GenericMaxFlow::DebugString(absl::string_view context, 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, " + "Residual capacity = %d, " + "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]); + context, arc, tail, head, residual_arc_capacity_[arc], + residual_arc_capacity_[Opposite(arc)], node_potential_[tail], + node_potential_[head], node_excess_[tail], node_excess_[head]); } template @@ -770,7 +781,6 @@ void GenericMaxFlow::InitializePreflow() { // 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. @@ -778,8 +788,19 @@ void GenericMaxFlow::InitializePreflow() { // 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.) std::fill(node_excess_.get(), node_excess_.get() + max_num_nodes, 0); - for (ArcIndex arc = 0; arc < num_arcs; ++arc) { - SetCapacityAndClearFlow(arc, Capacity(arc)); + + // Restart from a clear state with no flow, and an initial arc capacity. + const ArcIndex num_arcs = graph_->num_arcs(); + if constexpr (Graph::kHasNegativeReverseArcs) { + for (ArcIndex arc = 0; arc < num_arcs; ++arc) { + const ArcIndex opposite_arc = Opposite(arc); + residual_arc_capacity_[arc] += residual_arc_capacity_[opposite_arc]; + residual_arc_capacity_[opposite_arc] = 0; + } + } else { + for (ArcIndex arc = 0; arc < num_arcs; ++arc) { + residual_arc_capacity_[arc] = initial_capacity_[arc]; + } } // All the initial heights are zero except for the source whose height is @@ -787,11 +808,14 @@ void GenericMaxFlow::InitializePreflow() { std::fill(node_potential_.get(), node_potential_.get() + 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(). - std::fill(first_admissible_arc_.get(), - first_admissible_arc_.get() + max_num_nodes, Graph::kNilArc); + // Initially we set first_admissible_arc_ to the first arc in the iteration. + for (NodeIndex node = 0; node < num_nodes; ++node) { + first_admissible_arc_[node] = Graph::kNilArc; + for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) { + first_admissible_arc_[node] = arc; + break; + } + } } // Note(user): Calling this function will break the property on the node @@ -832,10 +856,8 @@ void GenericMaxFlow::PushFlowExcessBackToSource() { // 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) { + for (const ArcIndex arc : graph_->OutgoingArcs(source_)) { + if (Flow(arc) > 0) { arc_stack.push_back(arc); } } @@ -843,6 +865,9 @@ void GenericMaxFlow::PushFlowExcessBackToSource() { // Start the dfs on the subgraph formed by the direct arcs with positive flow. while (!arc_stack.empty()) { + DCHECK_GT(Flow(arc_stack.back()), 0) + << arc_stack.size() - 1 << " arc " << arc_stack.back() << " " + << Tail(arc_stack.back()) << "->" << Head(arc_stack.back()); const NodeIndex node = Head(arc_stack.back()); // If the node is visited, it means we have explored all its arcs and we @@ -866,9 +891,7 @@ void GenericMaxFlow::PushFlowExcessBackToSource() { (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(); + for (const ArcIndex arc : graph_->OutgoingArcs(node)) { const FlowQuantity flow = Flow(arc); const NodeIndex head = Head(arc); if (flow > 0 && !stored[head]) { @@ -887,12 +910,13 @@ void GenericMaxFlow::PushFlowExcessBackToSource() { // 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; + FlowQuantity flow_on_cycle = 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); + if (const FlowQuantity arc_flow = Flow(arc_on_cycle); + arc_flow <= flow_on_cycle) { + flow_on_cycle = arc_flow; first_saturated_index = i; } } @@ -902,10 +926,10 @@ void GenericMaxFlow::PushFlowExcessBackToSource() { // Cancel the flow on the cycle, and set visited[node] = false for // the node that will be backtracked over. - PushFlow(-max_flow, node, arc); + PushFlow(-flow_on_cycle, 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); + PushFlow(-flow_on_cycle, Tail(arc_on_cycle), arc_on_cycle); if (i >= first_saturated_index) { DCHECK(visited[Head(arc_on_cycle)]); visited[Head(arc_on_cycle)] = false; @@ -934,17 +958,17 @@ void GenericMaxFlow::PushFlowExcessBackToSource() { DCHECK(arc_stack.empty()); DCHECK(index_branch.empty()); - // Return the flow to the sink. Note that the sink_ and the source_ are not + // Return the flow to the source. 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); + for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) { + const FlowQuantity flow = Flow(arc); + if (flow < 0) { + DCHECK_GT(residual_arc_capacity_[arc], 0); + const FlowQuantity to_push = std::min(node_excess_[node], -flow); + PushFlow(to_push, node, arc); if (node_excess_[node] == 0) break; } } @@ -977,9 +1001,7 @@ void GenericMaxFlow::GlobalUpdate() { 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(); + for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) { const NodeIndex head = Head(arc); // Skip the arc if the height of head was already set to the correct @@ -1066,8 +1088,7 @@ bool GenericMaxFlow::SaturateOutgoingArcsFromSource() { 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(); + for (const ArcIndex arc : graph_->OutgoingArcs(source_)) { const FlowQuantity flow = residual_arc_capacity_[arc]; // This is a special IsAdmissible() condition for the source. @@ -1235,10 +1256,9 @@ void GenericMaxFlow::Discharge(const NodeIndex node) { while (true) { DCHECK(IsActive(node)); - for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node, - first_admissible_arc[node]); - it.Ok(); it.Next()) { - const ArcIndex arc = it.Index(); + for (const ArcIndex arc : + graph_->OutgoingOrOppositeIncomingArcsStartingFrom( + node, first_admissible_arc_[node])) { if (IsAdmissible(node, arc, node_potentials)) { DCHECK(IsActive(node)); const NodeIndex head = Head(arc); @@ -1272,9 +1292,7 @@ void GenericMaxFlow::Relabel(NodeIndex node) { // 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(); + for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) { if (residual_arc_capacity_[arc] > 0) { // Update min_height only for arcs with available capacity. NodeHeight head_height = node_potential_[Head(arc)]; @@ -1334,9 +1352,7 @@ void GenericMaxFlow::ComputeReachableNodes( 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(); + for (const ArcIndex arc : graph_->OutgoingOrOppositeIncomingArcs(node)) { const NodeIndex head = Head(arc); if (node_in_bfs_queue_[head]) continue; if (residual_arc_capacity_[reverse ? Opposite(arc) : arc] == 0) continue; diff --git a/ortools/graph/generic_max_flow_test.cc b/ortools/graph/generic_max_flow_test.cc index 8f3fb5982c..56470f0413 100644 --- a/ortools/graph/generic_max_flow_test.cc +++ b/ortools/graph/generic_max_flow_test.cc @@ -14,13 +14,16 @@ #include "ortools/graph/generic_max_flow.h" #include +#include #include #include #include +#include #include #include #include +#include "absl/random/bit_gen_ref.h" #include "absl/random/random.h" #include "absl/strings/str_format.h" #include "absl/types/span.h" @@ -56,19 +59,22 @@ typename GenericMaxFlow::Status MaxFlowTester( } 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]); + const int image = arc < permutation.size() ? permutation[arc] : arc; + + max_flow.SetArcCapacity(image, capacity[arc]); + EXPECT_EQ(max_flow.Capacity(image), 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; + for (int arc = 0; arc < num_arcs; ++arc) { + const int image = arc < permutation.size() ? permutation[arc] : arc; + + EXPECT_EQ(expected_flow[arc], max_flow.Flow(image)) << " arc = " << arc; } } @@ -94,7 +100,7 @@ class GenericMaxFlowTest : public ::testing::Test {}; typedef ::testing::Types, util::ReverseArcStaticGraph<>, - util::ReverseArcMixedGraph<> > + util::ReverseArcMixedGraph<>> GraphTypes; TYPED_TEST_SUITE(GenericMaxFlowTest, GraphTypes); @@ -284,8 +290,11 @@ void GenerateCompleteGraph(const typename Graph::NodeIndex num_tails, AddSourceAndSink(num_tails, num_heads, graph); } +// Generate a bipartite graph where each right node is connected to `degree` +// random nodes on the left. template -void GeneratePartialRandomGraph(const typename Graph::NodeIndex num_tails, +void GeneratePartialRandomGraph(absl::BitGenRef random, + const typename Graph::NodeIndex num_tails, const typename Graph::NodeIndex num_heads, const typename Graph::NodeIndex degree, Graph* graph) { @@ -293,11 +302,10 @@ void GeneratePartialRandomGraph(const typename Graph::NodeIndex num_tails, 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); + absl::Uniform(random, 0, num_heads); graph->AddArc(tail, head + num_tails); } } @@ -305,13 +313,13 @@ void GeneratePartialRandomGraph(const typename Graph::NodeIndex num_tails, } template -void GenerateRandomArcValuations(const Graph& graph, const int64_t max_range, +void GenerateRandomArcValuations(absl::BitGenRef random, 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); + (*arc_valuation)[arc] = absl::Uniform(random, 0, max_range); } } @@ -330,12 +338,14 @@ FlowQuantity SolveMaxFlow(GenericMaxFlow* max_flow) { 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; + const typename Graph::ArcIndex opposite_arc = graph->OppositeArc(arc); + EXPECT_EQ(max_flow->Flow(arc), -max_flow->Flow(opposite_arc)); + if (max_flow->Flow(arc) > 0) { + EXPECT_LE(max_flow->Flow(arc), max_flow->Capacity(arc)); + } else { + EXPECT_LE(0, max_flow->Flow(opposite_arc)); + EXPECT_LE(max_flow->Flow(opposite_arc), max_flow->Capacity(opposite_arc)); + } } return max_flow->GetOptimalFlow(); } @@ -365,13 +375,11 @@ FlowQuantity SolveMaxFlowWithLP(GenericMaxFlow* max_flow) { 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(); + for (const typename Graph::ArcIndex arc : graph->OutgoingArcs(source_index)) { objective->SetCoefficient(var[arc], -1.0); } solver.Solve(); - return static_cast(-objective->Value() + .5); + return static_cast(std::round(-objective->Value())); } template @@ -380,39 +388,56 @@ struct MaxFlowSolver { }; template -void FullRandomAssignment(typename MaxFlowSolver::Solver f, - typename Graph::NodeIndex num_tails, - typename Graph::NodeIndex num_heads, - FlowQuantity expected_flow1, - FlowQuantity expected_flow2) { +void FullAssignment(std::optional unused, + typename MaxFlowSolver::Solver f, + typename Graph::NodeIndex num_tails, + typename Graph::NodeIndex num_heads) { 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( + 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); + + // In a complete graph we should always reach the maximum flow. + const FlowQuantity flow = f(max_flow.get()); + EXPECT_EQ(std::min(num_tails, num_heads), flow); } template -void PartialRandomAssignment(typename MaxFlowSolver::Solver f, +void PartialRandomAssignment(std::optional expected_flow, + 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; + typename Graph::NodeIndex num_heads) { + absl::BitGen absl_random; + std::mt19937 mt_random(0); + absl::BitGenRef random = expected_flow != std::nullopt + ? absl::BitGenRef(mt_random) + : absl::BitGenRef(absl_random); + + const typename Graph::NodeIndex kDegree = 3; Graph graph; - GeneratePartialRandomGraph(num_tails, num_heads, kDegree, &graph); - Graphs::Build(&graph); - CHECK_EQ(graph.num_arcs(), num_tails * kDegree + num_tails + num_heads); + GeneratePartialRandomGraph(random, num_tails, num_heads, kDegree, &graph); std::vector arc_capacity(graph.num_arcs(), 1); - std::unique_ptr > max_flow(new GenericMaxFlow( + + std::vector permutation; + graph.Build(&permutation); + arc_capacity.resize(graph.num_arcs(), 0); + 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); + + const FlowQuantity flow = f(max_flow.get()); + if (expected_flow != std::nullopt) { + EXPECT_EQ(*expected_flow, flow); + } else { + // Use the LP as reference value. + const FlowQuantity expected_flow1 = SolveMaxFlowWithLP(max_flow.get()); + EXPECT_EQ(expected_flow1, flow); + } } template @@ -426,104 +451,127 @@ void ChangeCapacities(absl::Span arc_capacity, } template -void PartialRandomFlow(typename MaxFlowSolver::Solver f, +void PartialRandomFlow(std::optional expected_flow, + typename MaxFlowSolver::Solver f, typename Graph::NodeIndex num_tails, - typename Graph::NodeIndex num_heads, - FlowQuantity expected_flow1, - FlowQuantity expected_flow2) { + typename Graph::NodeIndex num_heads) { + absl::BitGen absl_random; + std::mt19937 mt_random(0); + absl::BitGenRef random = expected_flow != std::nullopt + ? absl::BitGenRef(mt_random) + : absl::BitGenRef(absl_random); + const typename Graph::NodeIndex kDegree = 10; const FlowQuantity kCapacityRange = 10000; const FlowQuantity kCapacityDelta = 1000; Graph graph; - GeneratePartialRandomGraph(num_tails, num_heads, kDegree, &graph); + GeneratePartialRandomGraph(random, num_tails, num_heads, kDegree, &graph); std::vector arc_capacity(graph.num_arcs()); - GenerateRandomArcValuations(graph, kCapacityRange, &arc_capacity); + GenerateRandomArcValuations(random, graph, kCapacityRange, &arc_capacity); std::vector permutation; Graphs::Build(&graph, &permutation); + arc_capacity.resize(graph.num_arcs(), 0); // In case Build() adds more arcs. util::Permute(permutation, &arc_capacity); - std::unique_ptr > max_flow(new GenericMaxFlow( + std::unique_ptr> max_flow(new GenericMaxFlow( &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); SetUpNetworkData(arc_capacity, max_flow.get()); + + if (expected_flow != std::nullopt) { + FlowQuantity flow = f(max_flow.get()); // Just solve once. + EXPECT_EQ(flow, *expected_flow); + return; + } + + const FlowQuantity expected_flow1 = SolveMaxFlowWithLP(max_flow.get()); FlowQuantity flow = f(max_flow.get()); EXPECT_EQ(expected_flow1, flow); + ChangeCapacities(arc_capacity, kCapacityDelta, max_flow.get()); + + const FlowQuantity expected_flow2 = SolveMaxFlowWithLP(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, +void FullRandomFlow(std::optional expected_flow, + typename MaxFlowSolver::Solver f, typename Graph::NodeIndex num_tails, - typename Graph::NodeIndex num_heads, - FlowQuantity expected_flow1, FlowQuantity expected_flow2) { + typename Graph::NodeIndex num_heads) { + absl::BitGen absl_random; + std::mt19937 mt_random(0); + absl::BitGenRef random = expected_flow != std::nullopt + ? absl::BitGenRef(mt_random) + : absl::BitGenRef(absl_random); + 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); + GenerateRandomArcValuations(random, graph, kCapacityRange, &arc_capacity); std::vector permutation; Graphs::Build(&graph, &permutation); + arc_capacity.resize(graph.num_arcs(), 0); // In case Build() adds more arcs. util::Permute(permutation, &arc_capacity); - std::unique_ptr > max_flow(new GenericMaxFlow( + std::unique_ptr> max_flow(new GenericMaxFlow( &graph, graph.num_nodes() - 2, graph.num_nodes() - 1)); SetUpNetworkData(arc_capacity, max_flow.get()); + + if (expected_flow != std::nullopt) { + const FlowQuantity flow = f(max_flow.get()); // Just solve once. + EXPECT_EQ(flow, expected_flow); + return; + } + + const FlowQuantity expected_flow1 = SolveMaxFlowWithLP(max_flow.get()); FlowQuantity flow = f(max_flow.get()); EXPECT_EQ(expected_flow1, flow); + ChangeCapacities(arc_capacity, kCapacityDelta, max_flow.get()); + const FlowQuantity expected_flow2 = SolveMaxFlowWithLP(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 LP_AND_FLOW_TEST(test_name, size) \ + TEST(MaxFlowStaticGraphTest, test_name##size) { \ + test_name>(std::nullopt, SolveMaxFlow, size, \ + size); \ + } \ + TEST(MaxFlowListGraphTest, test_name##size) { \ + test_name>(std::nullopt, SolveMaxFlow, size, \ + size); \ + } \ + TEST(MaxFlowStarGraphTest, test_name##size) { \ + test_name(std::nullopt, SolveMaxFlow, size, size); \ } -#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 +// These are absl::BitGen random test, so they will always work on different +// graphs. +LP_AND_FLOW_TEST(FullAssignment, 300); +LP_AND_FLOW_TEST(PartialRandomAssignment, 100); +LP_AND_FLOW_TEST(PartialRandomAssignment, 1000); +LP_AND_FLOW_TEST(PartialRandomFlow, 400); +LP_AND_FLOW_TEST(FullRandomFlow, 100); template static void BM_FullRandomAssignment(benchmark::State& state) { const int kSize = 3000; for (auto _ : state) { - FullRandomAssignment(SolveMaxFlow, kSize, kSize, kSize, kSize); + FullAssignment(std::nullopt, SolveMaxFlow, kSize, kSize); } state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); } @@ -532,7 +580,7 @@ template static void BM_PartialRandomAssignment(benchmark::State& state) { const int kSize = 10100; for (auto _ : state) { - PartialRandomAssignment(SolveMaxFlow, kSize, kSize, kSize, kSize); + PartialRandomAssignment(9512, SolveMaxFlow, kSize, kSize); } state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); } @@ -541,7 +589,7 @@ template static void BM_PartialRandomFlow(benchmark::State& state) { const int kSize = 800; for (auto _ : state) { - PartialRandomFlow(SolveMaxFlow, kSize, kSize, 3884850, 3112123); + PartialRandomFlow(3939172, SolveMaxFlow, kSize, kSize); } state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); } @@ -550,34 +598,33 @@ template static void BM_FullRandomFlow(benchmark::State& state) { const int kSize = 800; for (auto _ : state) { - FullRandomFlow(SolveMaxFlow, kSize, kSize, 4000549, 3239512); + FullRandomFlow(3952652, SolveMaxFlow, kSize, kSize); } state.SetItemsProcessed(static_cast(state.max_iterations) * kSize); } +// Note that these benchmark include the graph creation and generation... 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. diff --git a/ortools/graph/graph.h b/ortools/graph/graph.h index 943c5b29ae..4ea2a94948 100644 --- a/ortools/graph/graph.h +++ b/ortools/graph/graph.h @@ -167,6 +167,7 @@ #include "absl/base/port.h" #include "absl/debugging/leak_check.h" +#include "absl/log/check.h" #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/macros.h" @@ -194,6 +195,7 @@ class BaseGraph { // that you define a typedef ... Graph; for readability. typedef NodeIndexType NodeIndex; typedef ArcIndexType ArcIndex; + static constexpr bool kHasNegativeReverseArcs = HasReverseArcs; BaseGraph() : num_nodes_(0), @@ -2321,6 +2323,11 @@ class CompleteBipartiteGraph num_arcs_ = left_nodes * right_nodes; } + // Returns the arc index for the arc from `left` to `right`, where `left` is + // in `[0, left_nodes)` and `right` is in + // `[left_nodes, left_nodes + right_nodes)`. + ArcIndexType GetArc(NodeIndexType left_node, NodeIndexType right_node) const; + NodeIndexType Head(ArcIndexType arc) const; NodeIndexType Tail(ArcIndexType arc) const; ArcIndexType OutDegree(NodeIndexType node) const; @@ -2333,9 +2340,11 @@ class CompleteBipartiteGraph class OutgoingArcIterator { public: OutgoingArcIterator(const CompleteBipartiteGraph& graph, NodeIndexType node) - : index_(graph.right_nodes_ * node), - limit_(node >= graph.left_nodes_ ? index_ - : graph.right_nodes_ * (node + 1)) {} + : index_(static_cast(graph.right_nodes_) * node), + limit_(node >= graph.left_nodes_ + ? index_ + : static_cast(graph.right_nodes_) * + (node + 1)) {} bool Ok() const { return index_ < limit_; } ArcIndexType Index() const { return index_; } @@ -2351,6 +2360,16 @@ class CompleteBipartiteGraph const NodeIndexType right_nodes_; }; +template +ArcIndexType CompleteBipartiteGraph::GetArc( + NodeIndexType left_node, NodeIndexType right_node) const { + DCHECK_LT(left_node, left_nodes_); + DCHECK_GE(right_node, left_nodes_); + DCHECK_LT(right_node, num_nodes_); + return left_node * static_cast(right_nodes_) + + (right_node - left_nodes_); +} + template NodeIndexType CompleteBipartiteGraph::Head( ArcIndexType arc) const { @@ -2376,8 +2395,9 @@ IntegerRange CompleteBipartiteGraph::OutgoingArcs( NodeIndexType node) const { if (node < left_nodes_) { - return IntegerRange(right_nodes_ * node, - right_nodes_ * (node + 1)); + return IntegerRange( + static_cast(right_nodes_) * node, + static_cast(right_nodes_) * (node + 1)); } else { return IntegerRange(0, 0); } @@ -2388,7 +2408,8 @@ IntegerRange CompleteBipartiteGraph::OutgoingArcsStartingFrom( NodeIndexType node, ArcIndexType from) const { if (node < left_nodes_) { - return IntegerRange(from, right_nodes_ * (node + 1)); + return IntegerRange( + from, static_cast(right_nodes_) * (node + 1)); } else { return IntegerRange(0, 0); } diff --git a/ortools/graph/graph_io.h b/ortools/graph/graph_io.h index 39c3b2a759..81753667c2 100644 --- a/ortools/graph/graph_io.h +++ b/ortools/graph/graph_io.h @@ -48,6 +48,9 @@ enum GraphToStringFormat { // Ditto, but the adjacency lists are sorted. PRINT_GRAPH_ADJACENCY_LISTS_SORTED, + + // Dot format, can be visualized with Graphviz. + PRINT_GRAPH_DOT, }; template std::string GraphToString(const Graph& graph, GraphToStringFormat format); @@ -74,6 +77,15 @@ absl::Status WriteGraphToFile(const Graph& graph, const std::string& filename, template std::string GraphToString(const Graph& graph, GraphToStringFormat format) { std::string out; + if (format == PRINT_GRAPH_DOT) { + absl::StrAppend(&out, "digraph {\n"); + for (const auto arc : graph.AllForwardArcs()) { + absl::StrAppend(&out, " ", static_cast(graph.Tail(arc)), "->", + static_cast(graph.Head(arc)), ";\n"); + } + absl::StrAppend(&out, "}\n"); + return out; + } std::vector adj; for (const typename Graph::NodeIndex node : graph.AllNodes()) { if (format == PRINT_GRAPH_ARCS) { diff --git a/ortools/graph/linear_assignment.h b/ortools/graph/linear_assignment.h index 635fd275c9..4905e778ae 100644 --- a/ortools/graph/linear_assignment.h +++ b/ortools/graph/linear_assignment.h @@ -221,21 +221,22 @@ ABSL_DECLARE_FLAG(bool, assignment_stack_order); namespace operations_research { // This class does not take ownership of its underlying graph. -template +template class LinearSumAssignment { public: typedef typename GraphType::NodeIndex NodeIndex; typedef typename GraphType::ArcIndex ArcIndex; + typedef CostValue CostValueT; // Constructor for the case in which we will build the graph // incrementally as we discover arc costs, as might be done with any // of the dynamic graph representations such as StarGraph or ForwardStarGraph. LinearSumAssignment(const GraphType& graph, NodeIndex num_left_nodes); - // Constructor for the case in which the underlying graph cannot be - // built until after all the arc costs are known, as is the case - // with ForwardStarStaticGraph. In this case, the graph is passed to - // us later via the SetGraph() method, below. + // Constructor for the case in which the underlying graph cannot be built + // until after all the arc costs are known, as is the case with `StaticGraph`. + // In this case, the graph is passed to us later via the SetGraph() method, + // below. LinearSumAssignment(NodeIndex num_left_nodes, ArcIndex num_arcs); // This type is neither copyable nor movable. @@ -244,9 +245,9 @@ class LinearSumAssignment { ~LinearSumAssignment() {} - // Sets the graph used by the LinearSumAssignment instance, for use - // when the graph layout can be determined only after arc costs are - // set. This happens, for example, when we use a ForwardStarStaticGraph. + // Sets the graph used by the `LinearSumAssignment` instance, for use when the + // graph layout can be determined only after arc costs are set. This happens, + // for example, when we use a `StaticGraph`. void SetGraph(const GraphType* graph) { DCHECK(graph_ == nullptr); graph_ = graph; @@ -558,7 +559,7 @@ class LinearSumAssignment { // Minimum value of epsilon. When a flow is epsilon-optimal for // epsilon == kMinEpsilon, the flow is optimal. - static const CostValue kMinEpsilon; + static constexpr CostValue kMinEpsilon = 1; // Current value of epsilon, the cost scaling parameter. CostValue epsilon_; @@ -959,11 +960,8 @@ class LinearSumAssignment { // Implementation of out-of-line LinearSumAssignment template member // functions. -template -const CostValue LinearSumAssignment::kMinEpsilon = 1; - -template -LinearSumAssignment::LinearSumAssignment( +template +LinearSumAssignment::LinearSumAssignment( const GraphType& graph, const NodeIndex num_left_nodes) : graph_(&graph), num_left_nodes_(num_left_nodes), @@ -985,8 +983,8 @@ LinearSumAssignment::LinearSumAssignment( : static_cast( new ActiveNodeQueue())) {} -template -LinearSumAssignment::LinearSumAssignment( +template +LinearSumAssignment::LinearSumAssignment( const NodeIndex num_left_nodes, const ArcIndex num_arcs) : graph_(nullptr), num_left_nodes_(num_left_nodes), @@ -1008,8 +1006,9 @@ LinearSumAssignment::LinearSumAssignment( : static_cast( new ActiveNodeQueue())) {} -template -void LinearSumAssignment::SetArcCost(ArcIndex arc, CostValue cost) { +template +void LinearSumAssignment::SetArcCost(ArcIndex arc, + CostValue cost) { if (graph_ != nullptr) { DCHECK_GE(arc, 0); DCHECK_LT(arc, graph_->num_arcs()); @@ -1023,7 +1022,7 @@ void LinearSumAssignment::SetArcCost(ArcIndex arc, CostValue cost) { scaled_arc_cost_[arc] = cost; } -template +template class CostValueCycleHandler : public PermutationCycleHandler { public: explicit CostValueCycleHandler(std::vector* cost) @@ -1082,21 +1081,22 @@ class ArcIndexOrderingByTailNode { }; // Passes ownership of the cycle handler to the caller. -template +template PermutationCycleHandler* -LinearSumAssignment::ArcAnnotationCycleHandler() { - return new CostValueCycleHandler( +LinearSumAssignment::ArcAnnotationCycleHandler() { + return new CostValueCycleHandler( &scaled_arc_cost_); } -template -void LinearSumAssignment::OptimizeGraphLayout(GraphType* graph) { +template +void LinearSumAssignment::OptimizeGraphLayout( + GraphType* graph) { // The graph argument is only to give us a non-const-qualified // handle on the graph we already have. Any different graph is // nonsense. DCHECK_EQ(graph_, graph); const ArcIndexOrderingByTailNode compare(*graph_); - CostValueCycleHandler cycle_handler( + CostValueCycleHandler cycle_handler( &scaled_arc_cost_); TailArrayManager tail_array_manager(graph); tail_array_manager.BuildTailArrayFromAdjacencyListsIfForwardGraph(); @@ -1104,14 +1104,14 @@ void LinearSumAssignment::OptimizeGraphLayout(GraphType* graph) { tail_array_manager.ReleaseTailArrayIfForwardGraph(); } -template -CostValue LinearSumAssignment::NewEpsilon( +template +CostValue LinearSumAssignment::NewEpsilon( const CostValue current_epsilon) const { return std::max(current_epsilon / alpha_, kMinEpsilon); } -template -bool LinearSumAssignment::UpdateEpsilon() { +template +bool LinearSumAssignment::UpdateEpsilon() { CostValue new_epsilon = NewEpsilon(epsilon_); slack_relabeling_price_ = PriceChangeBound(epsilon_, new_epsilon, nullptr); epsilon_ = new_epsilon; @@ -1125,8 +1125,8 @@ bool LinearSumAssignment::UpdateEpsilon() { } // For production code that checks whether a left-side node is active. -template -inline bool LinearSumAssignment::IsActive( +template +inline bool LinearSumAssignment::IsActive( NodeIndex left_node) const { DCHECK_LT(left_node, num_left_nodes_); return matched_arc_[left_node] == GraphType::kNilArc; @@ -1135,8 +1135,8 @@ inline bool LinearSumAssignment::IsActive( // Only for debugging. Separate from the production IsActive() method // so that method can assert that its argument is a left-side node, // while for debugging we need to be able to test any node. -template -inline bool LinearSumAssignment::IsActiveForDebugging( +template +inline bool LinearSumAssignment::IsActiveForDebugging( NodeIndex node) const { if (node < num_left_nodes_) { return IsActive(node); @@ -1145,8 +1145,9 @@ inline bool LinearSumAssignment::IsActiveForDebugging( } } -template -void LinearSumAssignment::InitializeActiveNodeContainer() { +template +void LinearSumAssignment::InitializeActiveNodeContainer() { DCHECK(active_nodes_->Empty()); for (BipartiteLeftNodeIterator node_it(*graph_, num_left_nodes_); node_it.Ok(); node_it.Next()) { @@ -1167,8 +1168,8 @@ void LinearSumAssignment::InitializeActiveNodeContainer() { // nodes we unmatch here. If a matching arc is priced out, we will not // unmatch its endpoints since that element of the matching is // guaranteed not to change. -template -void LinearSumAssignment::SaturateNegativeArcs() { +template +void LinearSumAssignment::SaturateNegativeArcs() { total_excess_ = 0; for (BipartiteLeftNodeIterator node_it(*graph_, num_left_nodes_); node_it.Ok(); node_it.Next()) { @@ -1188,8 +1189,8 @@ void LinearSumAssignment::SaturateNegativeArcs() { } // Returns true for success, false for infeasible. -template -bool LinearSumAssignment::DoublePush(NodeIndex source) { +template +bool LinearSumAssignment::DoublePush(NodeIndex source) { DCHECK_GT(num_left_nodes_, source); DCHECK(IsActive(source)) << "Node " << source << "must be active (unmatched)!"; @@ -1226,8 +1227,8 @@ bool LinearSumAssignment::DoublePush(NodeIndex source) { return new_price >= price_lower_bound_; } -template -bool LinearSumAssignment::Refine() { +template +bool LinearSumAssignment::Refine() { SaturateNegativeArcs(); InitializeActiveNodeContainer(); while (total_excess_ > 0) { @@ -1267,9 +1268,10 @@ bool LinearSumAssignment::Refine() { // // This function is large enough that our suggestion that the compiler // inline it might be pointless. -template -inline typename LinearSumAssignment::ImplicitPriceSummary -LinearSumAssignment::BestArcAndGap(NodeIndex left_node) const { +template +inline typename LinearSumAssignment::ImplicitPriceSummary +LinearSumAssignment::BestArcAndGap( + NodeIndex left_node) const { DCHECK(IsActive(left_node)) << "Node " << left_node << " must be active (unmatched)!"; DCHECK_GT(epsilon_, 0); @@ -1307,8 +1309,8 @@ LinearSumAssignment::BestArcAndGap(NodeIndex left_node) const { // // Requires the precondition, explicitly computed in FinalizeSetup(), // that every left-side node has at least one incident arc. -template -inline CostValue LinearSumAssignment::ImplicitPrice( +template +inline CostValue LinearSumAssignment::ImplicitPrice( NodeIndex left_node) const { DCHECK_GT(num_left_nodes_, left_node); DCHECK_GT(epsilon_, 0); @@ -1343,8 +1345,8 @@ inline CostValue LinearSumAssignment::ImplicitPrice( } // Only for debugging. -template -bool LinearSumAssignment::AllMatched() const { +template +bool LinearSumAssignment::AllMatched() const { for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) { if (IsActiveForDebugging(node)) { return false; @@ -1354,8 +1356,8 @@ bool LinearSumAssignment::AllMatched() const { } // Only for debugging. -template -bool LinearSumAssignment::EpsilonOptimal() const { +template +bool LinearSumAssignment::EpsilonOptimal() const { for (BipartiteLeftNodeIterator node_it(*graph_, num_left_nodes_); node_it.Ok(); node_it.Next()) { const NodeIndex left_node = node_it.Index(); @@ -1389,8 +1391,8 @@ bool LinearSumAssignment::EpsilonOptimal() const { return true; } -template -bool LinearSumAssignment::FinalizeSetup() { +template +bool LinearSumAssignment::FinalizeSetup() { incidence_precondition_satisfied_ = true; // epsilon_ must be greater than kMinEpsilon so that in the case // where the largest arc cost is zero, we still do a Refine() @@ -1442,15 +1444,15 @@ bool LinearSumAssignment::FinalizeSetup() { return in_range; } -template -void LinearSumAssignment::ReportAndAccumulateStats() { +template +void LinearSumAssignment::ReportAndAccumulateStats() { total_stats_.Add(iteration_stats_); VLOG(3) << "Iteration stats: " << iteration_stats_.StatsString(); iteration_stats_.Clear(); } -template -bool LinearSumAssignment::ComputeAssignment() { +template +bool LinearSumAssignment::ComputeAssignment() { CHECK(graph_ != nullptr); bool ok = graph_->num_nodes() == 2 * num_left_nodes_; if (!ok) return false; @@ -1474,8 +1476,8 @@ bool LinearSumAssignment::ComputeAssignment() { return ok; } -template -CostValue LinearSumAssignment::GetCost() const { +template +CostValue LinearSumAssignment::GetCost() const { // It is illegal to call this method unless we successfully computed // an optimum assignment. DCHECK(success_); diff --git a/ortools/graph/linear_assignment_test.cc b/ortools/graph/linear_assignment_test.cc index 5212ceefd0..148de43f1d 100644 --- a/ortools/graph/linear_assignment_test.cc +++ b/ortools/graph/linear_assignment_test.cc @@ -107,7 +107,7 @@ typedef ::testing::Types< ForwardEbertGraph, ForwardStaticGraph, EbertGraph, ForwardEbertGraph, ForwardStaticGraph, StarGraph, ForwardStarGraph, - ForwardStarStaticGraph, util::ListGraph<>, util::ReverseArcListGraph<>> + util::ListGraph<>, util::ReverseArcListGraph<>> GraphTypesForAssignmentTestingWithGraphBuilder; TYPED_TEST_SUITE(LinearSumAssignmentTestWithGraphBuilder, @@ -479,23 +479,23 @@ class MacholWien // The test parameter specifies the problem size and the stack/queue active node // list flag. TEST_P(MacholWien, SolveHardProblem) { - typedef ForwardStarStaticGraph GraphType; + using Graph = ::util::CompleteBipartiteGraph<>; NodeIndex n = ::testing::get<0>(GetParam()); absl::SetFlag(&FLAGS_assignment_stack_order, ::testing::get<1>(GetParam())); - AssignmentProblemSetup setup(n, n * n, false); + Graph graph(n, n); + LinearSumAssignment assignment(graph, n); for (NodeIndex i = 0; i < n; ++i) { for (NodeIndex j = 0; j < n; ++j) { - CreateArcWithCost(i, n + j, i * j, setup.builder, - setup.assignment); + const ArcIndex arc = graph.GetArc(i, n + j); + assignment.SetArcCost(arc, i * j); } } - setup.Finalize(); - EXPECT_TRUE(setup.assignment->ComputeAssignment()); - for (LinearSumAssignment::BipartiteLeftNodeIterator node_it( - *setup.assignment); + EXPECT_TRUE(assignment.ComputeAssignment()); + for (LinearSumAssignment::BipartiteLeftNodeIterator node_it( + assignment); node_it.Ok(); node_it.Next()) { const NodeIndex left_node = node_it.Index(); - const NodeIndex right_node = setup.assignment->GetMate(left_node); + const NodeIndex right_node = assignment.GetMate(left_node); EXPECT_EQ(2 * n - 1, left_node + right_node); } } @@ -596,13 +596,9 @@ void BM_ConstructRandomAssignmentProblem(benchmark::State& state) { BENCHMARK_TEMPLATE2(BM_ConstructRandomAssignmentProblem, StarGraph, false); BENCHMARK_TEMPLATE2(BM_ConstructRandomAssignmentProblem, ForwardStarGraph, false); -BENCHMARK_TEMPLATE2(BM_ConstructRandomAssignmentProblem, ForwardStarStaticGraph, - false); BENCHMARK_TEMPLATE2(BM_ConstructRandomAssignmentProblem, StarGraph, true); BENCHMARK_TEMPLATE2(BM_ConstructRandomAssignmentProblem, ForwardStarGraph, true); -BENCHMARK_TEMPLATE2(BM_ConstructRandomAssignmentProblem, ForwardStarStaticGraph, - true); template void BM_ConstructRandomAssignmentProblemWithNewGraphApi( @@ -646,12 +642,8 @@ void BM_SolveRandomAssignmentProblem(benchmark::State& state) { BENCHMARK_TEMPLATE2(BM_SolveRandomAssignmentProblem, StarGraph, false); BENCHMARK_TEMPLATE2(BM_SolveRandomAssignmentProblem, ForwardStarGraph, false); -BENCHMARK_TEMPLATE2(BM_SolveRandomAssignmentProblem, ForwardStarStaticGraph, - false); BENCHMARK_TEMPLATE2(BM_SolveRandomAssignmentProblem, StarGraph, true); BENCHMARK_TEMPLATE2(BM_SolveRandomAssignmentProblem, ForwardStarGraph, true); -BENCHMARK_TEMPLATE2(BM_SolveRandomAssignmentProblem, ForwardStarStaticGraph, - true); template void BM_SolveRandomAssignmentProblemWithNewGraphApi(benchmark::State& state) { @@ -698,14 +690,10 @@ BENCHMARK_TEMPLATE2(BM_ConstructAndSolveRandomAssignmentProblem, StarGraph, false); BENCHMARK_TEMPLATE2(BM_ConstructAndSolveRandomAssignmentProblem, ForwardStarGraph, false); -BENCHMARK_TEMPLATE2(BM_ConstructAndSolveRandomAssignmentProblem, - ForwardStarStaticGraph, false); BENCHMARK_TEMPLATE2(BM_ConstructAndSolveRandomAssignmentProblem, StarGraph, true); BENCHMARK_TEMPLATE2(BM_ConstructAndSolveRandomAssignmentProblem, ForwardStarGraph, true); -BENCHMARK_TEMPLATE2(BM_ConstructAndSolveRandomAssignmentProblem, - ForwardStarStaticGraph, true); template void BM_ConstructAndSolveRandomAssignmentProblemWithNewGraphApi( diff --git a/ortools/graph/min_cost_flow.h b/ortools/graph/min_cost_flow.h index 3ac55df0c8..d12ad6a357 100644 --- a/ortools/graph/min_cost_flow.h +++ b/ortools/graph/min_cost_flow.h @@ -257,6 +257,11 @@ class MinCostFlowBase { // GenericMinCostFlow<> interface. class SimpleMinCostFlow : public MinCostFlowBase { public: + typedef int32_t NodeIndex; + typedef int32_t ArcIndex; + typedef int64_t FlowQuantity; + typedef int64_t CostValue; + // By default, the constructor takes no size. New node indices are created // lazily by AddArcWithCapacityAndUnitCost() or SetNodeSupply() such that the // set of valid nodes will always be [0, NumNodes()). @@ -389,14 +394,15 @@ class SimpleMinCostFlow : public MinCostFlowBase { // Note that the latter two are different than FlowQuantity and CostValue, which // are used for global, aggregated values and may need to be larger. // -// TODO(user): Avoid using the globally defined type CostValue and FlowQuantity. // Also uses the Arc*Type where there is no risk of overflow in more places. -template +template class GenericMinCostFlow : public MinCostFlowBase { public: typedef typename Graph::NodeIndex NodeIndex; typedef typename Graph::ArcIndex ArcIndex; + typedef int64_t CostValue; + typedef int64_t FlowQuantity; typedef typename Graph::IncomingArcIterator IncomingArcIterator; typedef typename Graph::OutgoingArcIterator OutgoingArcIterator; typedef typename Graph::OutgoingOrOppositeIncomingArcIterator diff --git a/ortools/graph/min_cost_flow_test.cc b/ortools/graph/min_cost_flow_test.cc index 623166a3b6..3d8d301b59 100644 --- a/ortools/graph/min_cost_flow_test.cc +++ b/ortools/graph/min_cost_flow_test.cc @@ -864,10 +864,10 @@ FLOW_ONLY_TEST_SG(FullRandomFlow, 3000, 5588622, 7140712); // 230s // blaze run -c opt --linkopt=-static [--run_under=perflab] -- // ortools/graph/min_cost_flow_test --benchmarks=all // --benchmark_min_iters=1 --heap_check= --benchmark_memory_usage -static void BM_MinCostFlowOnMultiMatchingProblem(benchmark::State& state) { +template +void BM_MinCostFlowOnMultiMatchingProblem(benchmark::State& state) { std::mt19937 my_random(12345); - const int kNumChannels = 20000; - const int kNumUsers = 20000; // Average probability of a user-channel pair being matched. const double kDensity = 1.0 / 200; const int kMaxChannelsPerUser = 5 * static_cast(kDensity * kNumChannels); @@ -909,12 +909,6 @@ static void BM_MinCostFlowOnMultiMatchingProblem(benchmark::State& state) { CHECK_LE(total_demand, kNumUsers * kMaxChannelsPerUser); for (auto _ : state) { - // We don't have more than 65536 nodes, so we use 16-bit integers to spare - // memory (and potentially speed up the code). Arc indices must be 32 bits - // though, since we have much more. - typedef util::ReverseArcStaticGraph< - /*NodeIndexType*/ uint16_t, /*ArcIndexType*/ int32_t> - Graph; Graph graph(/*num_nodes=*/kNumUsers + kNumChannels + 1, /*arc_capacity=*/kNumChannels * kNumUsers + kNumUsers); // We model this problem as a graph (on which we'll do a min-cost flow): @@ -935,17 +929,15 @@ static void BM_MinCostFlowOnMultiMatchingProblem(benchmark::State& state) { for (int j = 0; j < kNumUsers; ++j) { graph.AddArc(/*tail=*/kNumChannels + 1 + j, /*head=*/0); } - std::vector permutation; - graph.Build(&permutation); + std::vector permutation; + Graphs::Build(&graph, &permutation); // To spare memory, we added arcs in the right order, so that no permutation // is needed. See graph.h. CHECK(permutation.empty()); // To spare memory, the types are chosen as small as possible. - GenericMinCostFlow - min_cost_flow(&graph); + GenericMinCostFlow min_cost_flow( + &graph); // We also disable the feasibility check which takes a huge amount of // memory. @@ -954,7 +946,7 @@ static void BM_MinCostFlowOnMultiMatchingProblem(benchmark::State& state) { min_cost_flow.SetNodeSupply(/*node=*/0, /*supply=*/-total_demand); // Now, set the arcs capacity and cost, in the same order as we created them // above. - Graph::ArcIndex arc_index = 0; + typename Graph::ArcIndex arc_index = 0; for (int i = 0; i < kNumChannels; ++i) { min_cost_flow.SetNodeSupply( /*node=*/i + 1, /*supply=*/num_users_per_channel[i]); @@ -976,7 +968,20 @@ static void BM_MinCostFlowOnMultiMatchingProblem(benchmark::State& state) { } } -BENCHMARK(BM_MinCostFlowOnMultiMatchingProblem); +// We don't have more than 65536 nodes, so we use 16-bit integers to spare +// memory (and potentially speed up the code). Arc indices must be 32 bits +// though, since we have much more. +BENCHMARK(BM_MinCostFlowOnMultiMatchingProblem< + util::ReverseArcStaticGraph, int16_t, int32_t, + /*kNumChannels=*/20000, /*kNumUsers=*/20000>); +// We also benchmark with default parameter types and StarGraph for reference. +// We use fewer channels and users to avoid running out of memory. +BENCHMARK(BM_MinCostFlowOnMultiMatchingProblem< + ::util::ReverseArcListGraph<>, int64_t, int64_t, + /*kNumChannels=*/5000, /*kNumUsers=*/5000>); +BENCHMARK(BM_MinCostFlowOnMultiMatchingProblem); } // namespace } // namespace operations_research diff --git a/ortools/graph/python/linear_sum_assignment.cc b/ortools/graph/python/linear_sum_assignment.cc index b4b83340e1..568a86dfd8 100644 --- a/ortools/graph/python/linear_sum_assignment.cc +++ b/ortools/graph/python/linear_sum_assignment.cc @@ -16,8 +16,8 @@ #include "pybind11/pybind11.h" #include "pybind11/stl.h" -using ::operations_research::NodeIndex; using ::operations_research::SimpleLinearSumAssignment; +using NodeIndex = ::operations_research::SimpleLinearSumAssignment::NodeIndex; using ::pybind11::arg; PYBIND11_MODULE(linear_sum_assignment, m) { diff --git a/ortools/graph/samples/assignment_linear_sum_assignment.cc b/ortools/graph/samples/assignment_linear_sum_assignment.cc index 4c7fa9937c..0ef63fc9a8 100644 --- a/ortools/graph/samples/assignment_linear_sum_assignment.cc +++ b/ortools/graph/samples/assignment_linear_sum_assignment.cc @@ -19,6 +19,8 @@ #include #include #include + +#include "ortools/base/logging.h" // [END import] namespace operations_research {