diff --git a/ortools/algorithms/BUILD.bazel b/ortools/algorithms/BUILD.bazel index 7662f76dc5..a39b86759b 100644 --- a/ortools/algorithms/BUILD.bazel +++ b/ortools/algorithms/BUILD.bazel @@ -101,6 +101,7 @@ cc_library( "@com_google_absl//absl/base", "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", + "@com_google_absl//absl/numeric:bits", "@com_google_absl//absl/types:span", ], ) @@ -294,6 +295,7 @@ cc_library( "@com_google_absl//absl/random", "@com_google_absl//absl/random:distributions", "@com_google_absl//absl/strings", + "@com_google_absl//absl/strings:str_format", ], ) @@ -305,6 +307,7 @@ cc_library( ":set_cover_cc_proto", ":set_cover_model", "//ortools/base", + "//ortools/base:mathutil", "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", "@com_google_absl//absl/types:span", @@ -320,9 +323,10 @@ cc_library( ":set_cover_invariant", ":set_cover_model", "//ortools/base", - "@com_google_absl//absl/base:nullability", + "@com_google_absl//absl/base", "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", + "@com_google_absl//absl/numeric:bits", "@com_google_absl//absl/random:distributions", "@com_google_absl//absl/types:span", ], @@ -355,10 +359,30 @@ cc_library( "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", "@com_google_absl//absl/strings", + "@com_google_absl//absl/strings:str_format", "@com_google_absl//absl/strings:string_view", ], ) +cc_binary( + name = "set_cover_solve", + srcs = ["set_cover_solve.cc"], + deps = [ + ":set_cover_heuristics", + ":set_cover_invariant", + ":set_cover_model", + ":set_cover_reader", + #"//file/placer", + "//ortools/base", + "//ortools/base:timer", + "@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/time", + ], +) + cc_test( name = "set_cover_test", size = "medium", diff --git a/ortools/algorithms/n_choose_k.cc b/ortools/algorithms/n_choose_k.cc index 8799a43c9f..dcf350cabb 100644 --- a/ortools/algorithms/n_choose_k.cc +++ b/ortools/algorithms/n_choose_k.cc @@ -146,11 +146,10 @@ absl::StatusOr NChooseK(int64_t n, int64_t k) { if (k < 0) { return absl::InvalidArgumentError(absl::StrFormat("k is negative (%d)", k)); } - if (k > n) { - return absl::InvalidArgumentError( - absl::StrFormat("k=%d is greater than n=%d", k, n)); + if (k > n / 2) { + if (k > n) return 0; // No way to choose more than n elements from n. + k = n - k; } - if (k > n / 2) k = n - k; if (k == 0) return 1; if (n < std::numeric_limits::max() && !NChooseKIntermediateComputationOverflowsInt(n, k)) { diff --git a/ortools/algorithms/n_choose_k.h b/ortools/algorithms/n_choose_k.h index 937faffdce..8863bbbae1 100644 --- a/ortools/algorithms/n_choose_k.h +++ b/ortools/algorithms/n_choose_k.h @@ -23,10 +23,11 @@ namespace operations_research { // i.e., the binomial coefficient (n, k). // This is like std::exp(MathUtil::LogCombinations(n, k)), but faster, with // perfect accuracy, and returning an error iff the result would overflow an -// int64_t or if an argument is invalid (i.e., n < 0, k < 0, or k > n). +// int64_t or if an argument is invalid (i.e., n < 0, k < 0). Returns 0 for k > +// n. // // NOTE(user): If you need a variation of this, ask the authors: it's very easy -// to add. E.g., other int types, other behaviors (e.g., return 0 if k > n, or +// to add. E.g., other int types, other behaviors (e.g., return // std::numeric_limits::max() on overflow, etc). absl::StatusOr NChooseK(int64_t n, int64_t k); } // namespace operations_research diff --git a/ortools/algorithms/n_choose_k_test.cc b/ortools/algorithms/n_choose_k_test.cc index a57a1887c5..6ae6ff1bcf 100644 --- a/ortools/algorithms/n_choose_k_test.cc +++ b/ortools/algorithms/n_choose_k_test.cc @@ -28,6 +28,7 @@ #include "benchmark/benchmark.h" #include "gtest/gtest.h" #include "ortools/base/dump_vars.h" +//#include "ortools/base/fuzztest.h" #include "ortools/base/gmock.h" #include "ortools/base/mathutil.h" #include "ortools/util/flat_matrix.h" @@ -50,11 +51,7 @@ TEST(NChooseKTest, TrivialErrorCases) { HasSubstr("n is negative"))); EXPECT_THAT(NChooseK(x, -1), StatusIs(absl::StatusCode::kInvalidArgument, HasSubstr("k is negative"))); - if (x != kint64max) { - EXPECT_THAT(NChooseK(x, x + 1), - StatusIs(absl::StatusCode::kInvalidArgument, - HasSubstr("greater than n"))); - } + if (x != kint64max) EXPECT_THAT(NChooseK(x, x + 1), IsOkAndHolds(0)); ASSERT_FALSE(HasFailure()) << DUMP_VARS(t, x); } } @@ -310,13 +307,15 @@ void BM_NChooseK(benchmark::State& state) { } state.SetItemsProcessed(state.iterations() * kNumInputs); } -BENCHMARK(BM_NChooseK<30, operations_research::NChooseK>); // int32_t domain. -BENCHMARK( - BM_NChooseK<60, operations_research::NChooseK>); // int{32,64} domain. -BENCHMARK( - BM_NChooseK<100, operations_research::NChooseK>); // int{32,64,128} domain. -BENCHMARK( - BM_NChooseK<100, MathUtil::LogCombinations>); // int{32,64,128} domain. +// int32_t domain. +BENCHMARK(BM_NChooseK<30, operations_research::NChooseK>); + +// int{32,64} domain. +BENCHMARK(BM_NChooseK<60, operations_research::NChooseK>); + +// int{32,64,128} domain. +BENCHMARK(BM_NChooseK<100, operations_research::NChooseK>); +BENCHMARK(BM_NChooseK<100, MathUtil::LogCombinations>); } // namespace } // namespace operations_research diff --git a/ortools/algorithms/python/set_cover.cc b/ortools/algorithms/python/set_cover.cc index 141380358e..01ae1e0679 100644 --- a/ortools/algorithms/python/set_cover.cc +++ b/ortools/algorithms/python/set_cover.cc @@ -368,7 +368,7 @@ PYBIND11_MODULE(set_cover, m) { }) .def("next_solution", [](TrivialSolutionGenerator& heuristic, - const std::vector& focus) -> bool { + absl::Span focus) -> bool { return heuristic.NextSolution(VectorIntToVectorSubsetIndex(focus)); }); @@ -380,7 +380,7 @@ PYBIND11_MODULE(set_cover, m) { }) .def("next_solution", [](RandomSolutionGenerator& heuristic, - const std::vector& focus) -> bool { + absl::Span focus) -> bool { return heuristic.NextSolution(VectorIntToVectorSubsetIndex(focus)); }); @@ -392,13 +392,13 @@ PYBIND11_MODULE(set_cover, m) { }) .def("next_solution", [](GreedySolutionGenerator& heuristic, - const std::vector& focus) -> bool { + absl::Span focus) -> bool { return heuristic.NextSolution(VectorIntToVectorSubsetIndex(focus)); }) .def("next_solution", [](GreedySolutionGenerator& heuristic, - const std::vector& focus, - const std::vector& costs) -> bool { + absl::Span focus, + absl::Span costs) -> bool { return heuristic.NextSolution( VectorIntToVectorSubsetIndex(focus), VectorDoubleToSubsetCostVector(costs)); @@ -413,7 +413,7 @@ PYBIND11_MODULE(set_cover, m) { }) .def("next_solution", [](ElementDegreeSolutionGenerator& heuristic, - const std::vector& focus) -> bool { + absl::Span focus) -> bool { return heuristic.NextSolution(VectorIntToVectorSubsetIndex(focus)); }) .def("next_solution", @@ -439,8 +439,8 @@ PYBIND11_MODULE(set_cover, m) { }) .def("next_solution", [](LazyElementDegreeSolutionGenerator& heuristic, - const std::vector& focus, - const std::vector& costs) -> bool { + absl::Span focus, + absl::Span costs) -> bool { return heuristic.NextSolution( VectorIntToVectorSubsetIndex(focus), VectorDoubleToSubsetCostVector(costs)); @@ -474,7 +474,7 @@ PYBIND11_MODULE(set_cover, m) { return heuristic.NextSolution(num_iterations); }) .def("next_solution", - [](GuidedLocalSearch& heuristic, const std::vector& focus, + [](GuidedLocalSearch& heuristic, absl::Span focus, int num_iterations) -> bool { return heuristic.NextSolution(VectorIntToVectorSubsetIndex(focus), num_iterations); @@ -552,7 +552,7 @@ PYBIND11_MODULE(set_cover, m) { return {cleared.begin(), cleared.end()}; }); m.def("clear_most_covered_elements", - [](const std::vector& focus, BaseInt num_subsets, + [](absl::Span focus, BaseInt num_subsets, SetCoverInvariant* inv) -> std::vector { const std::vector cleared = ClearMostCoveredElements( VectorIntToVectorSubsetIndex(focus), num_subsets, inv); diff --git a/ortools/algorithms/samples/knapsack.cc b/ortools/algorithms/samples/knapsack.cc index c6fa5456bb..fc32b5df67 100644 --- a/ortools/algorithms/samples/knapsack.cc +++ b/ortools/algorithms/samples/knapsack.cc @@ -14,6 +14,7 @@ // [START program] // [START import] #include +#include #include #include #include diff --git a/ortools/algorithms/set_cover_heuristics.cc b/ortools/algorithms/set_cover_heuristics.cc index 5f5b666943..f403cb040b 100644 --- a/ortools/algorithms/set_cover_heuristics.cc +++ b/ortools/algorithms/set_cover_heuristics.cc @@ -14,12 +14,16 @@ #include "ortools/algorithms/set_cover_heuristics.h" #include +#include +#include #include #include #include #include +#include "absl/base/casts.h" #include "absl/log/check.h" +#include "absl/numeric/bits.h" #include "absl/random/distributions.h" #include "absl/random/random.h" #include "absl/types/span.h" @@ -218,26 +222,130 @@ class ComputationUsefulnessStats { SubsetToIntVector num_free_elements_; }; +namespace { +// Clearly not the fastest radix sort, but its complexity is the right one. +// Furthermore: +// - it is as memory-safe as std::vectors can be (no pointers), +// - no multiplication is performed, +// - it is stable +// - it handles the cases of signed and unsigned integers automatically, +// - bounds on the keys are optional, or they can be computed automatically, +// - based on those bounds, the number of passes is automatically computed, +// - a payload is associated to each key, and it is sorted in the same way +// as the keys. This payload can be a vector of integers or a vector of +// pointers to larger objects. +// TODO(user): Make it an independent library. +// - add support for decreasing counting sort, +// - make payloads optional, +// - support floats and doubles, +// - improve performance. +// - use vectorized code. +namespace internal { +uint32_t RawBits(uint32_t x) { return x; } // NOLINT +uint32_t RawBits(int x) { return absl::bit_cast(x); } // NOLINT +uint32_t RawBits(float x) { return absl::bit_cast(x); } // NOLINT +uint64_t RawBits(uint64_t x) { return x; } // NOLINT +uint64_t RawBits(int64_t x) { return absl::bit_cast(x); } // NOLINT +uint64_t RawBits(double x) { return absl::bit_cast(x); } // NOLINT + +inline uint32_t Bucket(uint32_t x, uint32_t shift, uint32_t radix) { + DCHECK_EQ(0, radix & (radix - 1)); // Must be a power of two. + // NOMUTANTS -- a way to compute the remainder of a division when radix is a + // power of two. + return (RawBits(x) >> shift) & (radix - 1); +} + +template +int NumBitsToRepresent(T value) { + DCHECK_LE(absl::countl_zero(RawBits(value)), sizeof(T) * CHAR_BIT); + return sizeof(T) * CHAR_BIT - absl::countl_zero(RawBits(value)); +} + +template +void UpdateCounters(uint32_t radix, int shift, std::vector& keys, + std::vector& counts) { + std::fill(counts.begin(), counts.end(), 0); + DCHECK_EQ(counts[0], 0); + DCHECK_EQ(0, radix & (radix - 1)); // Must be a power of two. + const auto num_keys = keys.size(); + for (int64_t i = 0; i < num_keys; ++i) { + ++counts[Bucket(keys[i], shift, radix)]; + } + // Now the counts will contain the sum of the sizes below and including each + // bucket. + for (uint64_t i = 1; i < radix; ++i) { + counts[i] += counts[i - 1]; + } +} + +template +void IncreasingCountingSort(uint32_t radix, int shift, std::vector& keys, + std::vector& payloads, + std::vector& scratch_keys, + std::vector& scratch_payloads, + std::vector& counts) { + DCHECK_EQ(0, radix & (radix - 1)); // Must be a power of two. + UpdateCounters(radix, shift, keys, counts); + const auto num_keys = keys.size(); + // In this order for stability. + for (int64_t i = num_keys - 1; i >= 0; --i) { + Counter& c = counts[Bucket(keys[i], shift, radix)]; + scratch_keys[c - 1] = keys[i]; + scratch_payloads[c - 1] = payloads[i]; + --c; + } + std::swap(keys, scratch_keys); + std::swap(payloads, scratch_payloads); +} +} // namespace internal + +template +void RadixSort(int radix_log, std::vector& keys, + std::vector& payloads, Key min_key, Key max_key) { + // range_log is the number of bits necessary to represent the max_key + // We could as well use max_key - min_key, but it is more expensive to + // compute. + const int range_log = internal::NumBitsToRepresent(max_key); + DCHECK_EQ(internal::NumBitsToRepresent(0), 0); + DCHECK_LE(internal::NumBitsToRepresent(std::numeric_limits::max()), + sizeof(Key) * CHAR_BIT); + const int radix = 1 << radix_log; // By definition. + std::vector counters(radix, 0); + std::vector scratch_keys(keys.size()); + std::vector scratch_payloads(payloads.size()); + for (int shift = 0; shift < range_log; shift += radix_log) { + DCHECK_LE(1 << shift, max_key); + internal::IncreasingCountingSort(radix, shift, keys, payloads, scratch_keys, + scratch_payloads, counters); + } +} +} // namespace + std::vector GetUncoveredElementsSortedByDegree( const SetCoverInvariant* const inv) { const BaseInt num_elements = inv->model()->num_elements(); - std::vector degree_sorted_elements; + std::vector degree_sorted_elements; // payloads degree_sorted_elements.reserve(num_elements); + std::vector keys; + keys.reserve(num_elements); + const SparseRowView& rows = inv->model()->rows(); + BaseInt max_degree = 0; for (ElementIndex element : inv->model()->ElementRange()) { // Already covered elements should not be considered. if (inv->coverage()[element] != 0) continue; degree_sorted_elements.push_back(element); + const BaseInt size = rows[element].size(); + max_degree = std::max(max_degree, size); + keys.push_back(size); } - const SparseRowView& rows = inv->model()->rows(); - // Sort indices by degree i.e. the size of the row corresponding to an - // element. - // NOMUTANTS -- The code still works if the sort is removed. - std::sort(degree_sorted_elements.begin(), degree_sorted_elements.end(), - [&rows](const ElementIndex a, const ElementIndex b) { - if (rows[a].size() < rows[b].size()) return true; - if (rows[a].size() == rows[b].size()) return a < b; - return false; - }); + RadixSort(11, keys, degree_sorted_elements, 1, max_degree); +#ifndef NDEBUG + BaseInt prev_key = -1; + for (const auto key : keys) { + DCHECK_LE(prev_key, key); + prev_key = key; + } +#endif return degree_sorted_elements; } @@ -358,6 +466,7 @@ bool LazyElementDegreeSolutionGenerator::NextSolution( std::vector degree_sorted_elements = GetUncoveredElementsSortedByDegree(inv_); const SparseRowView& rows = inv_->model()->rows(); + const SparseColumnView& columns = inv_->model()->columns(); ComputationUsefulnessStats stats(inv_, false); for (const ElementIndex element : degree_sorted_elements) { // No need to cover an element that is already covered. @@ -367,6 +476,12 @@ bool LazyElementDegreeSolutionGenerator::NextSolution( BaseInt best_subset_num_free_elts = 0; for (const SubsetIndex subset : rows[element]) { if (!in_focus[subset]) continue; + const Cost filtering_det = + Determinant(costs[subset], columns[subset].size(), best_subset_cost, + best_subset_num_free_elts); + // If the ratio with the initial number elements is greater, we skip this + // subset. + if (filtering_det > 0) continue; const BaseInt num_free_elements = inv_->ComputeNumFreeElements(subset); stats.Update(subset, num_free_elements); const Cost det = Determinant(costs[subset], num_free_elements, diff --git a/ortools/algorithms/set_cover_heuristics.h b/ortools/algorithms/set_cover_heuristics.h index 06fa732a06..f52b715490 100644 --- a/ortools/algorithms/set_cover_heuristics.h +++ b/ortools/algorithms/set_cover_heuristics.h @@ -16,7 +16,6 @@ #include -#include "absl/base/nullability.h" #include "absl/types/span.h" #include "ortools/algorithms/adjustable_k_ary_heap.h" #include "ortools/algorithms/set_cover_invariant.h" diff --git a/ortools/algorithms/set_cover_invariant.cc b/ortools/algorithms/set_cover_invariant.cc index d997bf5c67..a87b640525 100644 --- a/ortools/algorithms/set_cover_invariant.cc +++ b/ortools/algorithms/set_cover_invariant.cc @@ -22,6 +22,7 @@ #include "absl/types/span.h" #include "ortools/algorithms/set_cover_model.h" #include "ortools/base/logging.h" +#include "ortools/base/mathutil.h" namespace operations_research { @@ -71,7 +72,7 @@ bool SetCoverInvariant::CheckConsistency(ConsistencyLevel consistency) const { return true; } auto [cst, cvrg] = ComputeCostAndCoverage(is_selected_); - CHECK_EQ(cost_, cst); + CHECK(MathUtil::AlmostEquals(cost_, cst)); for (const ElementIndex element : model_->ElementRange()) { CHECK_EQ(cvrg[element], coverage_[element]); } @@ -448,6 +449,6 @@ void SetCoverInvariant::ImportSolutionFromProto( is_selected_[SubsetIndex(s)] = true; } Cost cost = message.cost(); - CHECK_EQ(cost, cost_); + CHECK(MathUtil::AlmostEquals(cost, cost_)); } } // namespace operations_research diff --git a/ortools/algorithms/set_cover_model.cc b/ortools/algorithms/set_cover_model.cc index 9c4c47d0a4..e0a79caed8 100644 --- a/ortools/algorithms/set_cover_model.cc +++ b/ortools/algorithms/set_cover_model.cc @@ -28,6 +28,7 @@ #include "absl/random/discrete_distribution.h" #include "absl/random/distributions.h" #include "absl/random/random.h" +#include "absl/strings/str_format.h" #include "ortools/algorithms/set_cover.pb.h" #include "ortools/base/logging.h" @@ -119,11 +120,17 @@ SetCoverModel SetCoverModel::GenerateRandomModelFrom( // Number of elements in the generated model, using the above vector. BaseInt num_elements_covered(0); + // Maximum number of tries to generate a random element that is not yet in + // the subset, or a random subset that does not contain the element. + const int kMaxTries = 10; + // Loop-local vector indicating whether the currently generated subset // contains an element. - ElementBoolVector subset_contains_element(num_elements, false); - + ElementBoolVector subset_already_contains_element(num_elements, false); for (SubsetIndex subset(0); subset.value() < num_subsets; ++subset) { + LOG_EVERY_N_SEC(INFO, 5) + << absl::StrFormat("Generating subset %d (%.1f%%)", subset.value(), + 100.0 * subset.value() / num_subsets); const BaseInt cardinality = DiscreteAffine(bitgen, column_dist, min_column_size, column_scale); model.columns_[subset].reserve(cardinality); @@ -136,22 +143,67 @@ SetCoverModel SetCoverModel::GenerateRandomModelFrom( element = ElementIndex(degree_dist(bitgen)); CHECK_LT(element.value(), num_elements); ++num_tries; - if (num_tries > 10) { - return SetCoverModel(); - } - } while (subset_contains_element[element]); + } while (num_tries < kMaxTries && + subset_already_contains_element[element]); ++model.num_nonzeros_; model.columns_[subset].push_back(element); - subset_contains_element[element] = true; + subset_already_contains_element[element] = true; if (!contains_element[element]) { contains_element[element] = true; ++num_elements_covered; } } for (const ElementIndex element : model.columns_[subset]) { - subset_contains_element[element] = false; + subset_already_contains_element[element] = false; } } + LOG(INFO) << "Finished genreating the model with " << num_elements_covered + << " elements covered."; + + // It can happen -- rarely in practice -- that some of the elements cannot be + // covered. Let's add them to randomly chosen subsets. + if (num_elements_covered != num_elements) { + LOG(INFO) << "Generated model with " << num_elements - num_elements_covered + << " elements that cannot be covered. Adding them to random " + "subsets."; + SubsetBoolVector element_already_in_subset(num_subsets, false); + for (ElementIndex element(0); element.value() < num_elements; ++element) { + LOG_EVERY_N_SEC(INFO, 5) << absl::StrFormat( + "Generating subsets for element %d (%.1f%%)", element.value(), + 100.0 * element.value() / num_elements); + if (!contains_element[element]) { + const BaseInt degree = + DiscreteAffine(bitgen, row_dist, min_row_size, row_scale); + std::vector generated_subsets; + generated_subsets.reserve(degree); + for (BaseInt i = 0; i < degree; ++i) { + int num_tries = 0; + SubsetIndex subset_index; + // The method is the same as above. + do { + subset_index = SubsetIndex(DiscreteAffine( + bitgen, column_dist, min_column_size, column_scale)); + ++num_tries; + } while (num_tries < kMaxTries && + element_already_in_subset[subset_index]); + ++model.num_nonzeros_; + model.columns_[subset_index].push_back(element); + element_already_in_subset[subset_index] = true; + generated_subsets.push_back(subset_index); + } + for (const SubsetIndex subset_index : generated_subsets) { + element_already_in_subset[subset_index] = false; + } + contains_element[element] = true; + ++num_elements_covered; + } + } + LOG(INFO) << "Finished generating subsets for elements that were not " + "covered in the original model."; + } + LOG(INFO) << "Finished generating the model. There are " + << num_elements - num_elements_covered << " uncovered elements."; + CHECK_EQ(num_elements_covered, num_elements); // TODO(user): if necessary, use a better distribution for the costs. @@ -335,6 +387,9 @@ SetCoverProto SetCoverModel::ExportModelAsProto() const { CHECK(elements_in_subsets_are_sorted_); SetCoverProto message; for (const SubsetIndex subset : SubsetRange()) { + LOG_EVERY_N_SEC(INFO, 5) + << absl::StrFormat("Exporting subset %d (%.1f%%)", subset.value(), + 100.0 * subset.value() / num_subsets()); SetCoverProto::Subset* subset_proto = message.add_subset(); subset_proto->set_cost(subset_costs_[subset]); SparseColumn column = columns_[subset]; @@ -343,6 +398,7 @@ SetCoverProto SetCoverModel::ExportModelAsProto() const { subset_proto->add_element(element.value()); } } + LOG(INFO) << "Finished exporting the model."; return message; } @@ -411,7 +467,7 @@ class StatsAccumulator { } SetCoverModel::Stats ComputeStats() const { - const BaseInt n = count_; + const int64_t n = count_; // Since the code is used on a known number of values, we can compute the // standard deviation exactly, even if the values are not all stored. const double stddev = @@ -445,13 +501,12 @@ SetCoverModel::Stats ComputeStats(std::vector sizes) { template std::vector ComputeDeciles(std::vector values) { const int kNumDeciles = 10; - std::vector deciles; - deciles.reserve(kNumDeciles); - const float step = values.size() / kNumDeciles; - for (int i = 1; i <= kNumDeciles; ++i) { - const size_t point = std::max(0, i * step - 1); + std::vector deciles(kNumDeciles, T{0}); + const double step = values.size() / kNumDeciles; + for (int i = 1; i < kNumDeciles; ++i) { + const size_t point = std::clamp(i * step, 0, values.size() - 1); std::nth_element(values.begin(), values.begin() + point, values.end()); - deciles.push_back(values[point]); + deciles[i] = values[point]; } return deciles; } @@ -463,7 +518,7 @@ SetCoverModel::Stats SetCoverModel::ComputeCostStats() { } SetCoverModel::Stats SetCoverModel::ComputeRowStats() { - std::vector row_sizes(num_elements(), 0); + std::vector row_sizes(num_elements(), 0); for (const SparseColumn& column : columns_) { for (const ElementIndex element : column) { ++row_sizes[element.value()]; @@ -473,15 +528,15 @@ SetCoverModel::Stats SetCoverModel::ComputeRowStats() { } SetCoverModel::Stats SetCoverModel::ComputeColumnStats() { - std::vector column_sizes(columns_.size()); + std::vector column_sizes(columns_.size()); for (const SubsetIndex subset : SubsetRange()) { column_sizes[subset.value()] = columns_[subset].size(); } return ComputeStats(std::move(column_sizes)); } -std::vector SetCoverModel::ComputeRowDeciles() const { - std::vector row_sizes(num_elements(), 0); +std::vector SetCoverModel::ComputeRowDeciles() const { + std::vector row_sizes(num_elements(), 0); for (const SparseColumn& column : columns_) { for (const ElementIndex element : column) { ++row_sizes[element.value()]; @@ -490,8 +545,8 @@ std::vector SetCoverModel::ComputeRowDeciles() const { return ComputeDeciles(std::move(row_sizes)); } -std::vector SetCoverModel::ComputeColumnDeciles() const { - std::vector column_sizes(columns_.size()); +std::vector SetCoverModel::ComputeColumnDeciles() const { + std::vector column_sizes(columns_.size()); for (const SubsetIndex subset : SubsetRange()) { column_sizes[subset.value()] = columns_[subset].size(); } @@ -502,16 +557,16 @@ namespace { // Returns the number of bytes needed to store x with a base-128 encoding. BaseInt Base128SizeInBytes(BaseInt x) { const uint64_t u = x == 0 ? 1 : static_cast(x); - return (64 - absl::countl_zero(u) + 6) / 7; + return (std::numeric_limits::digits - absl::countl_zero(u) + 6) / 7; } } // namespace SetCoverModel::Stats SetCoverModel::ComputeColumnDeltaSizeStats() const { StatsAccumulator acc; for (const SparseColumn& column : columns_) { - BaseInt previous = 0; + int64_t previous = 0; for (const ElementIndex element : column) { - const BaseInt delta = element.value() - previous; + const int64_t delta = element.value() - previous; previous = element.value(); acc.Register(Base128SizeInBytes(delta)); } @@ -522,9 +577,9 @@ SetCoverModel::Stats SetCoverModel::ComputeColumnDeltaSizeStats() const { SetCoverModel::Stats SetCoverModel::ComputeRowDeltaSizeStats() const { StatsAccumulator acc; for (const SparseRow& row : rows_) { - BaseInt previous = 0; + int64_t previous = 0; for (const SubsetIndex subset : row) { - const BaseInt delta = subset.value() - previous; + const int64_t delta = subset.value() - previous; previous = subset.value(); acc.Register(Base128SizeInBytes(delta)); } diff --git a/ortools/algorithms/set_cover_model.h b/ortools/algorithms/set_cover_model.h index 56e1950b4d..6d07f96c1f 100644 --- a/ortools/algorithms/set_cover_model.h +++ b/ortools/algorithms/set_cover_model.h @@ -261,10 +261,10 @@ class SetCoverModel { Stats ComputeColumnStats(); // Computes deciles on rows and returns a vector of deciles. - std::vector ComputeRowDeciles() const; + std::vector ComputeRowDeciles() const; // Computes deciles on columns and returns a vector of deciles. - std::vector ComputeColumnDeciles() const; + std::vector ComputeColumnDeciles() const; // Computes basic statistics on the deltas of the row and column elements and // returns a Stats structure. The deltas are computed as the difference diff --git a/ortools/algorithms/set_cover_reader.cc b/ortools/algorithms/set_cover_reader.cc index 73e2cfa945..eb309b04db 100644 --- a/ortools/algorithms/set_cover_reader.cc +++ b/ortools/algorithms/set_cover_reader.cc @@ -22,6 +22,7 @@ #include "absl/log/check.h" #include "absl/strings/numbers.h" #include "absl/strings/str_cat.h" +#include "absl/strings/str_format.h" #include "absl/strings/str_split.h" #include "absl/strings/string_view.h" #include "ortools/algorithms/set_cover.pb.h" @@ -111,6 +112,9 @@ SetCoverModel ReadOrlibScp(absl::string_view filename) { model.SetSubsetCost(subset.value(), cost); } for (ElementIndex element : ElementRange(num_rows)) { + LOG_EVERY_N_SEC(INFO, 5) + << absl::StrFormat("Reading element %d (%.1f%%)", element.value(), + 100.0 * element.value() / model.num_elements()); const RowEntryIndex row_size(reader.ParseNextInteger()); for (RowEntryIndex entry(0); entry < row_size; ++entry) { // Correct the 1-indexing. @@ -118,6 +122,7 @@ SetCoverModel ReadOrlibScp(absl::string_view filename) { model.AddElementToSubset(element.value(), subset); } } + LOG(INFO) << "Finished reading the model."; file->Close(file::Defaults()).IgnoreError(); model.CreateSparseRowView(); return model; @@ -131,7 +136,10 @@ SetCoverModel ReadOrlibRail(absl::string_view filename) { const ElementIndex num_rows(reader.ParseNextInteger()); const BaseInt num_cols(reader.ParseNextInteger()); model.ReserveNumSubsets(num_cols); - for (int subset(0); subset < num_cols; ++subset) { + for (BaseInt subset(0); subset < num_cols; ++subset) { + LOG_EVERY_N_SEC(INFO, 5) + << absl::StrFormat("Reading subset %d (%.1f%%)", subset, + 100.0 * subset / model.num_subsets()); const double cost(reader.ParseNextDouble()); model.SetSubsetCost(subset, cost); const ColumnEntryIndex column_size(reader.ParseNextInteger()); @@ -142,6 +150,7 @@ SetCoverModel ReadOrlibRail(absl::string_view filename) { model.AddElementToSubset(element.value(), subset); } } + LOG(INFO) << "Finished reading the model."; file->Close(file::Defaults()).IgnoreError(); model.CreateSparseRowView(); return model; @@ -149,8 +158,11 @@ SetCoverModel ReadOrlibRail(absl::string_view filename) { SetCoverModel ReadFimiDat(absl::string_view filename) { SetCoverModel model; + BaseInt subset(0); for (const std::string& line : FileLines(filename)) { - SubsetIndex subset(0); + LOG_EVERY_N_SEC(INFO, 5) + << absl::StrFormat("Reading subset %d (%.1f%%)", subset, + 100.0 * subset / model.num_subsets()); std::vector elements = absl::StrSplit(line, ' '); if (elements.back().empty() || elements.back()[0] == '\0') { elements.pop_back(); @@ -165,6 +177,7 @@ SetCoverModel ReadFimiDat(absl::string_view filename) { } ++subset; } + LOG(INFO) << "Finished reading the model."; model.CreateSparseRowView(); return model; } @@ -189,19 +202,27 @@ class LineWriter { LineWriter(File* file, int max_cols) : num_cols_(0), max_cols_(max_cols), line_(), file_(file) {} ~LineWriter() { Close(); } - void Write(BaseInt value) { - const std::string text = absl::StrCat(value, " "); + + void Write(absl::string_view text) { const int text_size = text.size(); - if (num_cols_ + text_size > max_cols_) { + if (!text.empty() && text_size + num_cols_ > max_cols_) { CHECK_OK(file::WriteString(file_, absl::StrCat(line_, "\n"), file::Defaults())); line_.clear(); - } else { - absl::StrAppend(&line_, text); - num_cols_ += text_size; + num_cols_ = 0; } + absl::StrAppend(&line_, text); + num_cols_ += text_size; + } + + void Write(BaseInt value) { Write(absl::StrCat(value, " ")); } + + void Write(double value) { Write(absl::StrFormat("%.17g ", value)); } + + void Close() { + CHECK_OK( + file::WriteString(file_, absl::StrCat(line_, "\n"), file::Defaults())); } - void Close() { CHECK_OK(file::WriteString(file_, line_, file::Defaults())); } private: int num_cols_; @@ -217,19 +238,25 @@ void WriteOrlibScp(const SetCoverModel& model, absl::string_view filename) { CHECK_OK(file::WriteString( file, absl::StrCat(model.num_elements(), " ", model.num_subsets(), "\n"), file::Defaults())); - LineWriter cost_writer(file, kMaxCols); - for (const SubsetIndex subset : model.SubsetRange()) { - cost_writer.Write(model.subset_costs()[subset]); - } - for (const ElementIndex element : model.ElementRange()) { - CHECK_OK(file::WriteString(file, - absl::StrCat(model.rows()[element].size(), "\n"), - file::Defaults())); - LineWriter row_writer(file, kMaxCols); - for (const SubsetIndex subset : model.rows()[element]) { - row_writer.Write(subset.value() + 1); + { // RAII for the file writer. + LineWriter cost_writer(file, kMaxCols); + for (const SubsetIndex subset : model.SubsetRange()) { + cost_writer.Write(model.subset_costs()[subset]); + } + for (const ElementIndex element : model.ElementRange()) { + LOG_EVERY_N_SEC(INFO, 5) + << absl::StrFormat("Writing element %d (%.1f%%)", element.value(), + 100.0 * element.value() / model.num_elements()); + CHECK_OK(file::WriteString( + file, absl::StrCat(model.rows()[element].size(), "\n"), + file::Defaults())); + LineWriter row_writer(file, kMaxCols); + for (const SubsetIndex subset : model.rows()[element]) { + row_writer.Write(subset.value() + 1); + } } } + LOG(INFO) << "Finished writing the model."; file->Close(file::Defaults()).IgnoreError(); } @@ -241,6 +268,9 @@ void WriteOrlibRail(const SetCoverModel& model, absl::string_view filename) { file, absl::StrCat(model.num_elements(), " ", model.num_subsets(), "\n"), file::Defaults())); for (const SubsetIndex subset : model.SubsetRange()) { + LOG_EVERY_N_SEC(INFO, 5) + << absl::StrFormat("Writing subset %d (%.1f%%)", subset.value(), + 100.0 * subset.value() / model.num_subsets()); CHECK_OK( file::WriteString(file, absl::StrCat(model.subset_costs()[subset], " ", @@ -251,6 +281,7 @@ void WriteOrlibRail(const SetCoverModel& model, absl::string_view filename) { writer.Write(element.value() + 1); } } + LOG(INFO) << "Finished writing the model."; file->Close(file::Defaults()).IgnoreError(); } diff --git a/ortools/algorithms/set_cover_solve.cc b/ortools/algorithms/set_cover_solve.cc new file mode 100644 index 0000000000..9cc010e982 --- /dev/null +++ b/ortools/algorithms/set_cover_solve.cc @@ -0,0 +1,257 @@ +// 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 +#include + +#include "absl/flags/flag.h" +#include "absl/log/check.h" +#include "absl/strings/match.h" +#include "absl/strings/str_join.h" +#include "absl/time/time.h" +#include "ortools/algorithms/set_cover_heuristics.h" +#include "ortools/algorithms/set_cover_invariant.h" +#include "ortools/algorithms/set_cover_model.h" +#include "ortools/algorithms/set_cover_reader.h" +#include "ortools/base/init_google.h" +#include "ortools/base/logging.h" +#include "ortools/base/timer.h" + +ABSL_FLAG(std::string, input, "", "REQUIRED: Input file name."); +ABSL_FLAG(std::string, input_fmt, "", + "REQUIRED: Input file format. Either proto, proto_bin, OrlibRail, " + "OrlibScp or FimiDat."); + +ABSL_FLAG(std::string, hint_sol, "", "Input file name for solution."); +ABSL_FLAG(std::string, hint_fmt, "", "Input file format for solution."); + +ABSL_FLAG(std::string, output, "", + "If non-empty, write the returned solution to the given file."); +ABSL_FLAG(std::string, output_fmt, "", + "If out is non-empty, use the given format for the output."); +ABSL_FLAG(std::string, output_model, "", + "If non-empty, write the set cover model to the given file. "); +ABSL_FLAG(std::string, output_model_fmt, "", + "If output_model is non-empty, use the given format for the output " + "model file. Can be proto, proto_bin, OrlibRail, OrlibScp."); + +ABSL_FLAG(bool, generate, false, "Generate a new model from the input model."); +ABSL_FLAG(int, num_elements_wanted, 0, + "Number of elements wanted in the new generated model."); +ABSL_FLAG(int, num_subsets_wanted, 0, + "Number of subsets wanted in the new generated model."); +ABSL_FLAG(float, row_scale, 1.0, "Row scale for the new generated model."); +ABSL_FLAG(float, column_scale, 1.0, + "Column scale for the new generated model."); +ABSL_FLAG(float, cost_scale, 1.0, "Cost scale for the new generated model."); + +ABSL_FLAG(std::string, generation, "", "First-solution generation method."); +ABSL_FLAG(std::string, improvement, "", "Solution improvement method."); +ABSL_FLAG(int, num_threads, 1, + "Number of threads to use by the underlying solver."); + +namespace operations_research { +using CL = SetCoverInvariant::ConsistencyLevel; + +void LogStats(std::string name, SetCoverModel* model) { + LOG(INFO) << ", " << name << ", num_elements, " << model->num_elements() + << ", num_subsets, " << model->num_subsets(); + LOG(INFO) << ", " << name << ", num_nonzeros, " << model->num_nonzeros() + << ", fill rate, " << model->FillRate(); + LOG(INFO) << ", " << name << ", cost, " + << model->ComputeCostStats().DebugString(); + + LOG(INFO) << ", " << name << ", num_rows, " << model->num_elements() + << ", rows sizes, " << model->ComputeRowStats().DebugString(); + LOG(INFO) << ", " << name << ", row size deciles, " + << absl::StrJoin(model->ComputeRowDeciles(), ", "); + LOG(INFO) << ", " << name << ", row delta byte size stats, " + << model->ComputeRowDeltaSizeStats().DebugString(); + LOG(INFO) << ", " << name << ", num_columns, " << model->num_subsets() + << ", columns sizes, " << model->ComputeColumnStats().DebugString(); + LOG(INFO) << ", " << name << ", column size deciles, " + << absl::StrJoin(model->ComputeColumnDeciles(), ", "); + LOG(INFO) << ", " << name << ", column delta byte size stats, " + << model->ComputeColumnDeltaSizeStats().DebugString(); + BaseInt num_singleton_rows = 0; + for (const ElementIndex element : model->ElementRange()) { + if (model->rows()[element].size() == 1) { + ++num_singleton_rows; + } + } + BaseInt num_singleton_columns = 0; + for (const SubsetIndex subset : model->SubsetRange()) { + if (model->columns()[subset].size() == 1) { + ++num_singleton_columns; + } + } + LOG(INFO) << ", " << name << ", num_singleton_rows, " << num_singleton_rows + << ", num_singleton_columns, " << num_singleton_columns; +} + +void LogCostAndTiming(std::string name, std::string algo, double cost, + BaseInt solution_cardinality, const WallTimer& timer) { + LOG(INFO) << ", " << name << ", " << algo << ", cost, " << cost + << ", solution_cardinality, " << solution_cardinality << ", " + << absl::ToInt64Microseconds(timer.GetDuration()) << "e-6, s"; +} + +void LogCostAndTiming(std::string name, std::string algo, + const SetCoverInvariant& inv, const WallTimer& timer) { + LogCostAndTiming(name, algo, inv.cost(), inv.trace().size(), timer); +} + +enum class FileFormat { + EMPTY, + ORLIB_SCP, + ORLIB_RAIL, + FIMI_DAT, + PROTO, + PROTO_BIN, + TXT, +}; + +FileFormat ParseFileFormat(const std::string& format_name) { + if (format_name.empty()) { + return FileFormat::EMPTY; + } else if (absl::EqualsIgnoreCase(format_name, "OrlibScp")) { + return FileFormat::ORLIB_SCP; + } else if (absl::EqualsIgnoreCase(format_name, "OrlibRail")) { + return FileFormat::ORLIB_RAIL; + } else if (absl::EqualsIgnoreCase(format_name, "FimiDat")) { + return FileFormat::FIMI_DAT; + } else if (absl::EqualsIgnoreCase(format_name, "proto")) { + return FileFormat::PROTO; + } else if (absl::EqualsIgnoreCase(format_name, "proto_bin")) { + return FileFormat::PROTO_BIN; + } else if (absl::EqualsIgnoreCase(format_name, "txt")) { + return FileFormat::TXT; + } else { + LOG(FATAL) << "Unsupported input format: " << format_name; + } +} + +SetCoverModel ReadModel(const std::string& input_file, + FileFormat input_format) { + switch (input_format) { + case FileFormat::ORLIB_SCP: + return ReadOrlibScp(input_file); + case FileFormat::ORLIB_RAIL: + return ReadOrlibRail(input_file); + case FileFormat::FIMI_DAT: + return ReadFimiDat(input_file); + case FileFormat::PROTO: + return ReadSetCoverProto(input_file, /*binary=*/false); + case FileFormat::PROTO_BIN: + return ReadSetCoverProto(input_file, /*binary=*/true); + default: + LOG(FATAL) << "Unsupported input format: " + << static_cast(input_format); + } +} + +SubsetBoolVector ReadSolution(const std::string& input_file, + FileFormat input_format) { + switch (input_format) { + case FileFormat::TXT: + return ReadSetCoverSolutionText(input_file); + case FileFormat::PROTO: + return ReadSetCoverSolutionProto(input_file, /*binary=*/false); + case FileFormat::PROTO_BIN: + return ReadSetCoverSolutionProto(input_file, /*binary=*/true); + default: + LOG(FATAL) << "Unsupported input format: " + << static_cast(input_format); + } +} + +void WriteModel(const SetCoverModel& model, const std::string& output_file, + FileFormat output_format) { + switch (output_format) { + case FileFormat::ORLIB_SCP: + WriteOrlibScp(model, output_file); + break; + case FileFormat::ORLIB_RAIL: + WriteOrlibRail(model, output_file); + break; + case FileFormat::PROTO: + WriteSetCoverProto(model, output_file, /*binary=*/false); + break; + case FileFormat::PROTO_BIN: + WriteSetCoverProto(model, output_file, /*binary=*/true); + break; + default: + LOG(FATAL) << "Unsupported output format: " + << static_cast(output_format); + } +} + +void WriteSolution(const SetCoverModel& model, const SubsetBoolVector& solution, + const std::string& output_file, FileFormat output_format) { + switch (output_format) { + case FileFormat::TXT: + WriteSetCoverSolutionText(model, solution, output_file); + break; + case FileFormat::PROTO: + WriteSetCoverSolutionProto(model, solution, output_file, + /*binary=*/false); + break; + case FileFormat::PROTO_BIN: + WriteSetCoverSolutionProto(model, solution, output_file, + /*binary=*/true); + break; + default: + LOG(FATAL) << "Unsupported output format: " + << static_cast(output_format); + } +} + +SetCoverInvariant RunLazyElementDegree(std::string name, SetCoverModel* model) { + SetCoverInvariant inv(model); + LazyElementDegreeSolutionGenerator element_degree(&inv); + WallTimer timer; + timer.Start(); + CHECK(element_degree.NextSolution()); + DCHECK(inv.CheckConsistency(CL::kCostAndCoverage)); + LogCostAndTiming(name, "LazyElementDegreeSolutionGenerator", inv.cost(), + inv.trace().size(), timer); + return inv; +} + +void Run() { + const auto& input = absl::GetFlag(FLAGS_input); + const auto& input_format = ParseFileFormat(absl::GetFlag(FLAGS_input_fmt)); + const auto& output = absl::GetFlag(FLAGS_output); + const auto& output_format = ParseFileFormat(absl::GetFlag(FLAGS_output_fmt)); + SetCoverModel model = ReadModel(input, input_format); + model.CreateSparseRowView(); + if (absl::GetFlag(FLAGS_generate)) { + model = SetCoverModel::GenerateRandomModelFrom( + model, absl::GetFlag(FLAGS_num_elements_wanted), + absl::GetFlag(FLAGS_num_subsets_wanted), absl::GetFlag(FLAGS_row_scale), + absl::GetFlag(FLAGS_column_scale), absl::GetFlag(FLAGS_cost_scale)); + } + if (!output.empty()) { + CHECK(output_format != FileFormat::EMPTY); + WriteModel(model, output, output_format); + } + LogStats(input, &model); + SetCoverInvariant inv = RunLazyElementDegree(input, &model); +} +} // namespace operations_research + +int main(int argc, char** argv) { + InitGoogle(argv[0], &argc, &argv, true); + operations_research::Run(); + return 0; +} diff --git a/ortools/base/mathutil.h b/ortools/base/mathutil.h index 0d24888d32..ab08ba8592 100644 --- a/ortools/base/mathutil.h +++ b/ortools/base/mathutil.h @@ -278,6 +278,61 @@ class MathUtil { // Note that if k > 15, this uses Stirling's approximation of log(n!). // The relative error is about 1/(1260*k^5) (which is 7.6e-10 when k=16). static double LogCombinations(int n, int k); + + // Tests whether two values are close enough to each other to be considered + // equal. This function is intended to be used mainly as a replacement for + // equality tests of floating point values in CHECK()s, and as a replacement + // for equality comparison against golden files. It is the same as == for + // integer types. The purpose of AlmostEquals() is to avoid false positive + // error reports in unit tests and regression tests due to minute differences + // in floating point arithmetic (for example, due to a different compiler). + // + // We cannot use simple equality to compare floating point values + // because floating point expressions often accumulate inaccuracies, and + // new compilers may introduce further variations in the values. + // + // Two values x and y are considered "almost equals" if: + // (a) Both values are very close to zero: x and y are in the range + // [-standard_error, standard_error] + // Normal calculations producing these values are likely to be dealing + // with meaningless residual values. + // -or- + // (b) The difference between the values is small: + // abs(x - y) <= standard_error + // -or- + // (c) The values are finite and the relative difference between them is + // small: + // abs (x - y) <= standard_error * max(abs(x), abs(y)) + // -or- + // (d) The values are both positive infinity or both negative infinity. + // + // Cases (b) and (c) are the same as MathUtils::NearByFractionOrMargin(x, y). + // + // standard_error is the corresponding MathLimits::kStdError constant. + // It is equivalent to 5 bits of mantissa error. See util/math/mathlimits.cc. + // + // Caveat: + // AlmostEquals() is not appropriate for checking long sequences of + // operations where errors may cascade (like extended sums or matrix + // computations), or where significant cancellation may occur + // (e.g., the expression (x+y)-x, where x is much larger than y). + // Both cases may produce errors in excess of standard_error. + // In any case, you should not test the results of calculations which have + // not been vetted for possible cancellation errors and the like. + template + static bool AlmostEquals(const T x, const T y) { + static_assert(std::is_arithmetic_v); + if constexpr (std::numeric_limits::is_integer) { + return x == y; + } else { + if (x == y) // Covers +inf and -inf, and is a shortcut for finite values. + return true; + + if (std::abs(x) <= 1e-6 && std::abs(y) <= 1e-6) return true; + if (std::abs(x - y) <= 1e-6) return true; + return std::abs(x - y) / std::max(std::abs(x), std::abs(y)) <= 1e-6; + } + } }; } // namespace operations_research diff --git a/ortools/constraint_solver/samples/nurses_cp.cc b/ortools/constraint_solver/samples/nurses_cp.cc index 7af24d3dd2..c270ff6dfe 100644 --- a/ortools/constraint_solver/samples/nurses_cp.cc +++ b/ortools/constraint_solver/samples/nurses_cp.cc @@ -12,6 +12,7 @@ // limitations under the License. #include +#include #include #include // std::iota #include diff --git a/ortools/glpk/glpk_formatters.cc b/ortools/glpk/glpk_formatters.cc index a024c92956..9f8630b19e 100644 --- a/ortools/glpk/glpk_formatters.cc +++ b/ortools/glpk/glpk_formatters.cc @@ -13,6 +13,7 @@ #include "ortools/glpk/glpk_formatters.h" +#include #include #include #include diff --git a/ortools/glpk/glpk_formatters.h b/ortools/glpk/glpk_formatters.h index 67290ae76a..231efac8fd 100644 --- a/ortools/glpk/glpk_formatters.h +++ b/ortools/glpk/glpk_formatters.h @@ -15,6 +15,7 @@ #ifndef OR_TOOLS_GLPK_GLPK_FORMATTERS_H_ #define OR_TOOLS_GLPK_GLPK_FORMATTERS_H_ +#include #include #include diff --git a/ortools/graph/BUILD.bazel b/ortools/graph/BUILD.bazel index 79918ea430..c06aa21bcc 100644 --- a/ortools/graph/BUILD.bazel +++ b/ortools/graph/BUILD.bazel @@ -114,6 +114,7 @@ cc_test( ":minimum_vertex_cover", "//ortools/base:gmock_main", "@com_google_absl//absl/algorithm:container", + "@com_google_benchmark//:benchmark", ], ) @@ -495,6 +496,7 @@ cc_library( "//ortools/base", "//ortools/util:stats", "//ortools/util:zvector", + "@com_google_absl//absl/memory", "@com_google_absl//absl/strings", ], ) @@ -579,6 +581,7 @@ cc_binary( ":min_cost_flow", "//ortools/base", "//ortools/base:file", + "//ortools/base:path", "//ortools/base:status_macros", "//ortools/base:timer", "//ortools/util:filelineiter", diff --git a/ortools/graph/generic_max_flow.h b/ortools/graph/generic_max_flow.h index 8976696f13..194226151c 100644 --- a/ortools/graph/generic_max_flow.h +++ b/ortools/graph/generic_max_flow.h @@ -123,11 +123,14 @@ #ifndef OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_ #define OR_TOOLS_GRAPH_GENERIC_MAX_FLOW_H_ +#include #include +#include #include #include #include +#include "absl/memory/memory.h" #include "absl/strings/string_view.h" #include "ortools/base/logging.h" #include "ortools/graph/ebert_graph.h" @@ -427,7 +430,7 @@ class GenericMaxFlow : public MaxFlowStatusClass { const Graph* graph_; // An array representing the excess for each node in graph_. - std::vector node_excess_; + std::unique_ptr 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 @@ -440,7 +443,7 @@ class GenericMaxFlow : public MaxFlowStatusClass { // 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_; + std::unique_ptr 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 @@ -461,7 +464,7 @@ class GenericMaxFlow : public MaxFlowStatusClass { ZVector residual_arc_capacity_; // An array representing the first admissible arc for each node in graph_. - std::vector first_admissible_arc_; + std::unique_ptr 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(). @@ -566,9 +569,12 @@ GenericMaxFlow::GenericMaxFlow(const Graph* graph, NodeIndex 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); + 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); bfs_queue_.reserve(max_num_nodes); } const ArcIndex max_num_arcs = Graphs::ArcReservation(*graph_); @@ -771,20 +777,21 @@ 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.) - node_excess_.assign(max_num_nodes, 0); + std::fill(node_excess_.get(), node_excess_.get() + 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); + 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(). - first_admissible_arc_.assign(max_num_nodes, Graph::kNilArc); + std::fill(first_admissible_arc_.get(), + first_admissible_arc_.get() + max_num_nodes, Graph::kNilArc); } // Note(user): Calling this function will break the property on the node @@ -956,7 +963,7 @@ void GenericMaxFlow::GlobalUpdate() { 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(); + FlowQuantity* const node_excess = node_excess_.get(); // We do a BFS in the reverse residual graph, starting from the sink. // Because all the arcs from the source are saturated (except in @@ -1222,9 +1229,9 @@ void GenericMaxFlow::Discharge(const NodeIndex node) { 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(); + FlowQuantity* const node_excess = node_excess_.get(); + NodeHeight* const node_potentials = node_potential_.get(); + ArcIndex* const first_admissible_arc = first_admissible_arc_.get(); while (true) { DCHECK(IsActive(node)); diff --git a/ortools/graph/minimum_vertex_cover_test.cc b/ortools/graph/minimum_vertex_cover_test.cc index 6c63ae091f..821b0e1677 100644 --- a/ortools/graph/minimum_vertex_cover_test.cc +++ b/ortools/graph/minimum_vertex_cover_test.cc @@ -16,6 +16,7 @@ #include #include "absl/algorithm/container.h" +#include "benchmark/benchmark.h" #include "gtest/gtest.h" namespace operations_research { @@ -74,5 +75,21 @@ TEST(BipartiteMinimumVertexCoverTest, Empty) { 6); } +void BM_CompleteBipartite(benchmark::State& state) { + const int num_left = state.range(0); + const int num_right = state.range(1); + std::vector> left_to_right = + MakeCompleteBipartiteGraph(num_left, num_right); + for (auto _ : state) { + BipartiteMinimumVertexCover(left_to_right, num_right); + } +} +BENCHMARK(BM_CompleteBipartite) + ->ArgPair(1, 128) + ->ArgPair(128, 1) + ->ArgPair(32, 32) + ->ArgPair(8, 64) + ->ArgPair(64, 8); + } // namespace } // namespace operations_research diff --git a/ortools/graph/samples/assignment_min_flow.cc b/ortools/graph/samples/assignment_min_flow.cc index 4586be1bfc..c9056d1498 100644 --- a/ortools/graph/samples/assignment_min_flow.cc +++ b/ortools/graph/samples/assignment_min_flow.cc @@ -13,6 +13,7 @@ // [START program] // [START import] +#include #include #include diff --git a/ortools/graph/samples/balance_min_flow.cc b/ortools/graph/samples/balance_min_flow.cc index f34408d26e..ebce98ea4b 100644 --- a/ortools/graph/samples/balance_min_flow.cc +++ b/ortools/graph/samples/balance_min_flow.cc @@ -13,6 +13,7 @@ // [START program] // [START import] +#include #include #include diff --git a/ortools/graph/samples/simple_max_flow_program.cc b/ortools/graph/samples/simple_max_flow_program.cc index b3e5f44728..39eaf7bdfc 100644 --- a/ortools/graph/samples/simple_max_flow_program.cc +++ b/ortools/graph/samples/simple_max_flow_program.cc @@ -14,6 +14,7 @@ // [START program] // From Taha 'Introduction to Operations Research', example 6.4-2.""" // [START import] +#include #include #include diff --git a/ortools/graph/solve_flow_model.cc b/ortools/graph/solve_flow_model.cc index 7d3118bd86..1165d00d84 100644 --- a/ortools/graph/solve_flow_model.cc +++ b/ortools/graph/solve_flow_model.cc @@ -37,6 +37,7 @@ #include "ortools/base/helpers.h" #include "ortools/base/init_google.h" #include "ortools/base/logging.h" +#include "ortools/base/path.h" #include "ortools/base/timer.h" #include "ortools/graph/flow_problem.pb.h" #include "ortools/graph/graph.h" @@ -303,7 +304,7 @@ int main(int argc, char** argv) { "file_name, parsing_time, loading_time, solving_time, optimal_cost\n"); for (int i = 0; i < file_list.size(); ++i) { const std::string file_name = file_list[i]; - absl::PrintF("%s,", file_name); + absl::PrintF("%s,", file::Basename(file_name)); fflush(stdout); // Parse the input as a proto. diff --git a/ortools/linear_solver/model_exporter.cc b/ortools/linear_solver/model_exporter.cc index b8ecab58ab..dbcf974b8c 100644 --- a/ortools/linear_solver/model_exporter.cc +++ b/ortools/linear_solver/model_exporter.cc @@ -167,8 +167,8 @@ class MPModelProtoExporter { // a name, and a value. If the line is not empty, then only the pair // (name, value) is appended. The number of columns, limited to 2 by the MPS // format is also taken care of. - void AppendMpsTermWithContext(const std::string& head_name, - const std::string& name, double value, + void AppendMpsTermWithContext(absl::string_view head_name, + absl::string_view name, double value, std::string* output); // Appends a new-line if two columns are already present on the MPS line. @@ -693,9 +693,10 @@ void MPModelProtoExporter::AppendMpsLineHeaderWithNewLine( absl::StrAppend(output, "\n"); } -void MPModelProtoExporter::AppendMpsTermWithContext( - const std::string& head_name, const std::string& name, double value, - std::string* output) { +void MPModelProtoExporter::AppendMpsTermWithContext(absl::string_view head_name, + absl::string_view name, + double value, + std::string* output) { if (current_mps_column_ == 0) { AppendMpsLineHeader("", head_name, output); } diff --git a/ortools/linear_solver/samples/stigler_diet.cc b/ortools/linear_solver/samples/stigler_diet.cc index 75608631ad..bef17b49bc 100644 --- a/ortools/linear_solver/samples/stigler_diet.cc +++ b/ortools/linear_solver/samples/stigler_diet.cc @@ -15,6 +15,7 @@ // The Stigler diet problem. // [START import] #include +#include #include #include #include // std::pair