diff --git a/ortools/set_cover/BUILD.bazel b/ortools/set_cover/BUILD.bazel index 286720493c..81af52cf08 100644 --- a/ortools/set_cover/BUILD.bazel +++ b/ortools/set_cover/BUILD.bazel @@ -45,6 +45,8 @@ cc_library( deps = [ "//ortools/base:intops", "//ortools/base:strong_vector", + "//ortools/base:timer", + "@abseil-cpp//absl/time", ], ) @@ -54,12 +56,14 @@ cc_library( hdrs = ["set_cover_lagrangian.h"], deps = [ ":base_types", + ":set_cover_heuristics", ":set_cover_invariant", ":set_cover_model", "//ortools/algorithms:adjustable_k_ary_heap", "//ortools/base:threadpool", "@abseil-cpp//absl/log:check", "@abseil-cpp//absl/synchronization", + "@abseil-cpp//absl/types:span", ], ) @@ -115,6 +119,8 @@ cc_library( "@abseil-cpp//absl/log:check", "@abseil-cpp//absl/numeric:bits", "@abseil-cpp//absl/random:distributions", + "@abseil-cpp//absl/strings:string_view", + "@abseil-cpp//absl/time", "@abseil-cpp//absl/types:span", ], ) @@ -125,12 +131,14 @@ cc_library( hdrs = ["set_cover_mip.h"], deps = [ ":base_types", + ":set_cover_heuristics", ":set_cover_invariant", ":set_cover_model", "//ortools/linear_solver", "//ortools/lp_data:base", "@abseil-cpp//absl/log", "@abseil-cpp//absl/log:check", + "@abseil-cpp//absl/strings:string_view", "@abseil-cpp//absl/types:span", ], ) @@ -189,8 +197,10 @@ cc_binary( ":set_cover_invariant", ":set_cover_model", ":set_cover_reader", + "//ortools/base", "//ortools/base:timer", + "@abseil-cpp//absl/base:core_headers", "@abseil-cpp//absl/flags:flag", "@abseil-cpp//absl/log", "@abseil-cpp//absl/log:check", diff --git a/ortools/set_cover/base_types.h b/ortools/set_cover/base_types.h index 222d0725fb..678cb84626 100644 --- a/ortools/set_cover/base_types.h +++ b/ortools/set_cover/base_types.h @@ -16,8 +16,10 @@ #include +#include "absl/time/time.h" #include "ortools/base/strong_int.h" #include "ortools/base/strong_vector.h" +#include "ortools/base/timer.h" namespace operations_research { @@ -57,14 +59,37 @@ using SparseRow = util_intops::StrongVector; using ElementToIntVector = util_intops::StrongVector; using SubsetToIntVector = util_intops::StrongVector; -// Views of the sparse vectors. These need not be aligned as it's their contents -// that need to be aligned. +// Views of the sparse vectors. using SparseColumnView = util_intops::StrongVector; using SparseRowView = util_intops::StrongVector; using SubsetBoolVector = util_intops::StrongVector; using ElementBoolVector = util_intops::StrongVector; +// Simple stopwatch class that enables recording the time spent in various +// functions in the code. +// It uses RAII to automatically record the time spent in the constructor and +// destructor, independently of the path taken by the code. +class StopWatch { + public: + explicit StopWatch(absl::Duration* duration) : duration_(duration), timer_() { + timer_.Start(); + } + ~StopWatch() { + timer_.Stop(); + *duration_ = timer_.GetDuration(); + } + // Returns the elapsed time in seconds at a given moment. To be use to + // implement time limits in the derived classes. + double elapsed_time_in_seconds() const { return timer_.Get(); } + + absl::Duration GetElapsedDuration() const { return timer_.GetDuration(); } + + private: + absl::Duration* duration_; + WallTimer timer_; +}; + } // namespace operations_research #endif // OR_TOOLS_SET_COVER_BASE_TYPES_H_ diff --git a/ortools/set_cover/python/set_cover.cc b/ortools/set_cover/python/set_cover.cc index bde95db86e..d91bc33c5d 100644 --- a/ortools/set_cover/python/set_cover.cc +++ b/ortools/set_cover/python/set_cover.cc @@ -130,10 +130,13 @@ PYBIND11_MODULE(set_cover, m) { .def_readwrite("median", &SetCoverModel::Stats::median) .def_readwrite("mean", &SetCoverModel::Stats::mean) .def_readwrite("stddev", &SetCoverModel::Stats::stddev) - .def("debug_string", &SetCoverModel::Stats::DebugString); + .def_property_readonly("to_string", &SetCoverModel::Stats::ToString) + .def_property_readonly("to_verbose_string", + &SetCoverModel::Stats::ToVerboseString); py::class_(m, "SetCoverModel") .def(py::init<>()) + .def_property_readonly("name", &SetCoverModel::name) .def_property_readonly("num_elements", &SetCoverModel::num_elements) .def_property_readonly("num_subsets", &SetCoverModel::num_subsets) .def_property_readonly("num_nonzeros", &SetCoverModel::num_nonzeros) @@ -204,6 +207,7 @@ PYBIND11_MODULE(set_cover, m) { }); return subsets; }) + .def("set_name", &SetCoverModel::SetName) .def("add_empty_subset", &SetCoverModel::AddEmptySubset, arg("cost")) .def( "add_element_to_last_subset", @@ -227,9 +231,9 @@ PYBIND11_MODULE(set_cover, m) { .def("sort_elements_in_subsets", &SetCoverModel::SortElementsInSubsets) .def("compute_feasibility", &SetCoverModel::ComputeFeasibility) .def( - "reserve_num_subsets", + "resize_num_subsets", [](SetCoverModel& model, BaseInt num_subsets) { - model.ReserveNumSubsets(num_subsets); + model.ResizeNumSubsets(num_subsets); }, arg("num_subsets")) .def( @@ -379,7 +383,8 @@ PYBIND11_MODULE(set_cover, m) { const std::vector& in_focus) -> bool { return heuristic.NextSolution( BoolVectorToSubsetBoolVector(in_focus)); - }); + }) + .def("name", &TrivialSolutionGenerator::name); py::class_(m, "RandomSolutionGenerator") .def(py::init()) diff --git a/ortools/set_cover/set_cover_heuristics.cc b/ortools/set_cover/set_cover_heuristics.cc index 13bd5139ad..739dfae68e 100644 --- a/ortools/set_cover/set_cover_heuristics.cc +++ b/ortools/set_cover/set_cover_heuristics.cc @@ -46,6 +46,7 @@ using CL = SetCoverInvariant::ConsistencyLevel; bool TrivialSolutionGenerator::NextSolution( absl::Span focus) { + StopWatch stop_watch(&run_time_); const SubsetIndex num_subsets(model()->num_subsets()); SubsetBoolVector choices(num_subsets, false); for (const SubsetIndex subset : focus) { @@ -60,6 +61,7 @@ bool TrivialSolutionGenerator::NextSolution( bool RandomSolutionGenerator::NextSolution( absl::Span focus) { + StopWatch stop_watch(&run_time_); inv()->ClearTrace(); std::vector shuffled(focus.begin(), focus.end()); std::shuffle(shuffled.begin(), shuffled.end(), absl::BitGen()); @@ -78,6 +80,7 @@ bool RandomSolutionGenerator::NextSolution( bool GreedySolutionGenerator::NextSolution( absl::Span focus) { + StopWatch stop_watch(&run_time_); DCHECK(inv()->CheckConsistency(CL::kCostAndCoverage)); inv()->Recompute(CL::kFreeAndUncovered); inv()->ClearTrace(); @@ -358,6 +361,7 @@ double Determinant(Cost c1, BaseInt n1, Cost c2, BaseInt n2) { bool ElementDegreeSolutionGenerator::NextSolution( const SubsetBoolVector& in_focus) { + StopWatch stop_watch(&run_time_); DVLOG(1) << "Entering ElementDegreeSolutionGenerator::NextSolution"; inv()->Recompute(CL::kFreeAndUncovered); // Create the list of all the indices in the problem. @@ -415,6 +419,7 @@ bool ElementDegreeSolutionGenerator::NextSolution( bool LazyElementDegreeSolutionGenerator::NextSolution( const SubsetBoolVector& in_focus) { + StopWatch stop_watch(&run_time_); inv()->CompressTrace(); DVLOG(1) << "Entering LazyElementDegreeSolutionGenerator::NextSolution"; DCHECK(inv()->CheckConsistency(CL::kCostAndCoverage)); @@ -464,6 +469,7 @@ bool LazyElementDegreeSolutionGenerator::NextSolution( // SteepestSearch. bool SteepestSearch::NextSolution(const SubsetBoolVector& in_focus) { + StopWatch stop_watch(&run_time_); const int64_t num_iterations = max_iterations(); DCHECK(inv()->CheckConsistency(CL::kCostAndCoverage)); inv()->Recompute(CL::kFreeAndUncovered); @@ -519,6 +525,7 @@ bool SteepestSearch::NextSolution(const SubsetBoolVector& in_focus) { // LazySteepestSearch. bool LazySteepestSearch::NextSolution(const SubsetBoolVector& in_focus) { + StopWatch stop_watch(&run_time_); const int64_t num_iterations = max_iterations(); DCHECK(inv()->CheckConsistency(CL::kCostAndCoverage)); DVLOG(1) << "Entering LazySteepestSearch::NextSolution, num_iterations = " @@ -614,6 +621,7 @@ void GuidedTabuSearch::UpdatePenalties(absl::Span focus) { } bool GuidedTabuSearch::NextSolution(absl::Span focus) { + StopWatch stop_watch(&run_time_); const int64_t num_iterations = max_iterations(); DCHECK(inv()->CheckConsistency(CL::kFreeAndUncovered)); DVLOG(1) << "Entering GuidedTabuSearch::NextSolution, num_iterations = " @@ -723,6 +731,7 @@ Cost GuidedLocalSearch::ComputeDelta(SubsetIndex subset) const { } bool GuidedLocalSearch::NextSolution(absl::Span focus) { + StopWatch stop_watch(&run_time_); const int64_t num_iterations = max_iterations(); inv()->Recompute(CL::kRedundancy); Cost best_cost = inv()->cost(); diff --git a/ortools/set_cover/set_cover_heuristics.h b/ortools/set_cover/set_cover_heuristics.h index 105f78d72b..03743c53b3 100644 --- a/ortools/set_cover/set_cover_heuristics.h +++ b/ortools/set_cover/set_cover_heuristics.h @@ -18,8 +18,11 @@ #include #include +#include #include +#include "absl/strings/string_view.h" +#include "absl/time/time.h" #include "absl/types/span.h" #include "ortools/algorithms/adjustable_k_ary_heap.h" #include "ortools/set_cover/base_types.h" @@ -44,7 +47,6 @@ namespace operations_research { // The term 'focus' hereafter means a subset of the S_j's designated by their // indices. Focus make it possible to run the algorithms on the corresponding // subproblems. -// // Base class for all set cover solution generators. This is almost an // interface. @@ -53,13 +55,27 @@ class SetCoverSolutionGenerator { // By default, the maximum number of iterations is set to infinity, and the // maximum time in seconds is set to infinity as well (and the time limit is // not yet implemented). - explicit SetCoverSolutionGenerator(SetCoverInvariant* inv) - : inv_(inv), - max_iterations_(std::numeric_limits::max()), - max_time_in_seconds_(std::numeric_limits::infinity()) {} + SetCoverSolutionGenerator(SetCoverInvariant* inv, + const absl::string_view name) + : run_time_(absl::ZeroDuration()), + inv_(inv), + name_(name), + time_limit_in_seconds_(std::numeric_limits::infinity()), + max_iterations_(std::numeric_limits::max()) {} virtual ~SetCoverSolutionGenerator() = default; + void SetName(const absl::string_view name) { name_ = name; } + + SetCoverInvariant* inv() const { return inv_; } + + // Resets the limits to their default values. + virtual SetCoverSolutionGenerator& ResetLimits() { + time_limit_in_seconds_ = std::numeric_limits::infinity(); + max_iterations_ = std::numeric_limits::max(); + return *this; + } + // Sets the maximum number of iterations. SetCoverSolutionGenerator& SetMaxIterations(int64_t max_iterations) { max_iterations_ = max_iterations; @@ -69,7 +85,29 @@ class SetCoverSolutionGenerator { // Returns the maximum number of iterations. int64_t max_iterations() const { return max_iterations_; } - // TODO(user): add a max_time_in_seconds_ getter and setter. + // Sets the time limit in seconds. + SetCoverSolutionGenerator& SetTimeLimitInSeconds(double seconds) { + time_limit_in_seconds_ = seconds; + return *this; + } + + absl::Duration run_time() const { return run_time_; } + + // Returns the total elapsed runtime in seconds. + double run_time_in_seconds() const { + return absl::ToDoubleSeconds(run_time_); + } + + // Returns the total elapsed runtime in microseconds. + double run_time_in_microseconds() const { + return absl::ToInt64Microseconds(run_time_); + } + + // Returns the name of the heuristic. + std::string name() const { return name_; } + + // Returns the current cost of the solution in the invariant. + Cost cost() const { return inv_->cost(); } // Virtual methods that must be implemented by derived classes. @@ -85,10 +123,52 @@ class SetCoverSolutionGenerator { protected: // Accessors. - SetCoverInvariant* inv() const { return inv_; } SetCoverModel* model() const { return inv_->model(); } BaseInt num_subsets() const { return model()->num_subsets(); } + // The time limit in seconds. + double time_limit_in_seconds() const { return time_limit_in_seconds_; } + + // run_time_ is an abstract duration for the time spent in NextSolution(). + absl::Duration run_time_; + + private: + // The data structure that will maintain the invariant for the model. + SetCoverInvariant* inv_; + + // The name of the heuristic. + std::string name_; + + // The time limit in seconds. + double time_limit_in_seconds_; + + // The maximum number of iterations. + int64_t max_iterations_; +}; + +// Now we define two classes that are used to implement the solution +// generators. The first one uses a vector of subset indices as focus, while +// the second one uses a vector of Booleans. This makes the declaration of the +// solution generators shorter, more readable and improves type correctness. + +// The class of solution generators that use a vector of subset indices as +// focus, with a transformation from a vector of Booleans to a vector of +// subset indices if needed. +class SubsetListBasedSolutionGenerator : public SetCoverSolutionGenerator { + public: + explicit SubsetListBasedSolutionGenerator(SetCoverInvariant* inv, + absl::string_view name) + : SetCoverSolutionGenerator(inv, name) {} + + bool NextSolution(absl::Span _) override { return false; } + + bool NextSolution() final { return NextSolution(model()->all_subsets()); } + + bool NextSolution(const SubsetBoolVector& in_focus) final { + return NextSolution(MakeSubsetIndexSpan(in_focus)); + } + + private: // Converts a vector of Booleans to a vector of subset indices. // TODO(user): this should not be, but a better iterator system should be // implemented. @@ -104,7 +184,28 @@ class SetCoverSolutionGenerator { } return absl::MakeConstSpan(result); } +}; +// The class of solution generators that use a vector of Booleans as focus, +// with a transformation from a vector of subset indices to a vector of +// Booleans if needed. +class BoolVectorBasedSolutionGenerator : public SetCoverSolutionGenerator { + public: + explicit BoolVectorBasedSolutionGenerator(SetCoverInvariant* inv, + absl::string_view name) + : SetCoverSolutionGenerator(inv, name) {} + + bool NextSolution(const SubsetBoolVector& _) override { return false; } + + bool NextSolution(absl::Span focus) final { + return NextSolution(MakeBoolVector(focus, num_subsets())); + } + + bool NextSolution() final { + return NextSolution(SubsetBoolVector(num_subsets(), true)); + } + + private: // Converts a vector of subset indices to a vector of Booleans. // TODO(user): this should not be, but a better iterator system should be // implemented. @@ -116,16 +217,6 @@ class SetCoverSolutionGenerator { } return result; } - - private: - // The data structure that will maintain the invariant for the model. - SetCoverInvariant* inv_; - - // The maximum number of iterations. - int64_t max_iterations_; - - // The maximum time in seconds. - double max_time_in_seconds_; }; // An obvious idea is to take all the S_j's (or equivalently to set all the @@ -133,19 +224,14 @@ class SetCoverSolutionGenerator { // using local search. // The consistency level is maintained up to kFreeAndUncovered. -class TrivialSolutionGenerator : public SetCoverSolutionGenerator { +class TrivialSolutionGenerator : public SubsetListBasedSolutionGenerator { public: - explicit TrivialSolutionGenerator(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv) {} + explicit TrivialSolutionGenerator(SetCoverInvariant* inv, + absl::string_view name = "TrivialGenerator") + : SubsetListBasedSolutionGenerator(inv, name) {} - bool NextSolution() final { return NextSolution(model()->all_subsets()); } - - // The implementation works on a list of subsets. + using SubsetListBasedSolutionGenerator::NextSolution; bool NextSolution(absl::Span focus) final; - - bool NextSolution(const SubsetBoolVector& in_focus) final { - return NextSolution(MakeSubsetIndexSpan(in_focus)); - } }; // A slightly more complicated but better way to compute a first solution is @@ -155,19 +241,14 @@ class TrivialSolutionGenerator : public SetCoverSolutionGenerator { // the generator towards the columns with the least marginal costs. // The consistency level is maintained up to kFreeAndUncovered. -class RandomSolutionGenerator : public SetCoverSolutionGenerator { +class RandomSolutionGenerator : public SubsetListBasedSolutionGenerator { public: - explicit RandomSolutionGenerator(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv) {} + explicit RandomSolutionGenerator(SetCoverInvariant* inv, + absl::string_view name = "RandomGenerator") + : SubsetListBasedSolutionGenerator(inv, name) {} - bool NextSolution() final { return NextSolution(model()->all_subsets()); } - - // The implementation works on a list of subsets. + using SubsetListBasedSolutionGenerator::NextSolution; bool NextSolution(absl::Span focus) final; - - bool NextSolution(const SubsetBoolVector& in_focus) final { - return NextSolution(MakeSubsetIndexSpan(in_focus)); - } }; // The first solution is obtained using the Chvatal heuristic, that @@ -191,19 +272,14 @@ class RandomSolutionGenerator : public SetCoverSolutionGenerator { // https://doi.org/10.1145/1871437.1871501. // The consistency level is maintained up to kFreeAndUncovered. -class GreedySolutionGenerator : public SetCoverSolutionGenerator { +class GreedySolutionGenerator : public SubsetListBasedSolutionGenerator { public: - explicit GreedySolutionGenerator(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv) {} + explicit GreedySolutionGenerator(SetCoverInvariant* inv, + absl::string_view name = "GreedyGenerator") + : SubsetListBasedSolutionGenerator(inv, name) {} - bool NextSolution() final { return NextSolution(model()->all_subsets()); } - - // The implementation works on a list of subsets. + using SubsetListBasedSolutionGenerator::NextSolution; bool NextSolution(absl::Span focus) final; - - bool NextSolution(const SubsetBoolVector& in_focus) final { - return NextSolution(MakeSubsetIndexSpan(in_focus)); - } }; // Solution generator based on the degree of elements. @@ -211,50 +287,38 @@ class GreedySolutionGenerator : public SetCoverSolutionGenerator { // The generator consists in iteratively choosing a non-covered element with // the smallest degree, and selecting a subset that covers it with the least // ratio cost / number of uncovered elements. The number of uncovered elements -// are updated for each impacted subset. The newly-covered elements degrees are -// also updated and set to zero. +// are updated for each impacted subset. The newly-covered elements degrees +// are also updated and set to zero. // The consistency level is maintained up to kFreeAndUncovered. -class ElementDegreeSolutionGenerator : public SetCoverSolutionGenerator { +class ElementDegreeSolutionGenerator : public BoolVectorBasedSolutionGenerator { public: - explicit ElementDegreeSolutionGenerator(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv) {} + explicit ElementDegreeSolutionGenerator( + SetCoverInvariant* inv, absl::string_view name = "ElementDegreeGenerator") + : BoolVectorBasedSolutionGenerator(inv, name) {} - bool NextSolution() final { - return NextSolution(SubsetBoolVector(num_subsets(), true)); - } - - bool NextSolution(absl::Span focus) final { - return NextSolution(MakeBoolVector(focus, num_subsets())); - } - - // The implementation works on a vector of Booleans. + using BoolVectorBasedSolutionGenerator::NextSolution; bool NextSolution(const SubsetBoolVector& in_focus) final; }; // Solution generator based on the degree of elements. -// The heuristic is the same as ElementDegreeSolutionGenerator, but the number -// of uncovered elements for a subset is computed on-demand. In empirical tests, -// this is faster than ElementDegreeSolutionGenerator because a very small -// percentage of need to be computed, and even fewer among them need to be -// computed again later on. +// The heuristic is the same as ElementDegreeSolutionGenerator, but the number +// of uncovered elements for a subset is computed on-demand. In empirical +// tests, this is faster than ElementDegreeSolutionGenerator because a very +// small percentage of need to be computed, and even fewer among them need to +// be computed again later on. // Because the number of uncovered elements is computed on-demand, the // consistency level only needs to be set to kCostAndCoverage. -class LazyElementDegreeSolutionGenerator : public SetCoverSolutionGenerator { +class LazyElementDegreeSolutionGenerator + : public BoolVectorBasedSolutionGenerator { public: - explicit LazyElementDegreeSolutionGenerator(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv) {} + explicit LazyElementDegreeSolutionGenerator( + SetCoverInvariant* inv, + absl::string_view name = "LazyElementDegreeGenerator") + : BoolVectorBasedSolutionGenerator(inv, name) {} - bool NextSolution() final { - return NextSolution(SubsetBoolVector(num_subsets(), true)); - } - - bool NextSolution(absl::Span focus) final { - return NextSolution(MakeBoolVector(focus, num_subsets())); - } - - // The implementation works on a vector of Booleans. + using BoolVectorBasedSolutionGenerator::NextSolution; bool NextSolution(const SubsetBoolVector& in_focus) final; }; @@ -266,20 +330,13 @@ class LazyElementDegreeSolutionGenerator : public SetCoverSolutionGenerator { // direction, taking the S_j with the largest total cost. // The consistency level is maintained up to kFreeAndUncovered. -class SteepestSearch : public SetCoverSolutionGenerator { +class SteepestSearch : public BoolVectorBasedSolutionGenerator { public: - explicit SteepestSearch(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv) {} + explicit SteepestSearch(SetCoverInvariant* inv, + absl::string_view name = "SteepestSearch") + : BoolVectorBasedSolutionGenerator(inv, name) {} - bool NextSolution() final { - return NextSolution(SubsetBoolVector(num_subsets(), true)); - } - - bool NextSolution(absl::Span focus) final { - return NextSolution(MakeBoolVector(focus, num_subsets())); - } - - // The implementation works on a vector of Booleans. + using BoolVectorBasedSolutionGenerator::NextSolution; bool NextSolution(const SubsetBoolVector& in_focus) final; }; @@ -288,20 +345,13 @@ class SteepestSearch : public SetCoverSolutionGenerator { // priorities are computed when needed. It is faster to compute because // there are relatively few subsets in the solution, because the cardinality // of the solution is bounded by the number of elements. -class LazySteepestSearch : public SetCoverSolutionGenerator { +class LazySteepestSearch : public BoolVectorBasedSolutionGenerator { public: - explicit LazySteepestSearch(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv) {} + explicit LazySteepestSearch(SetCoverInvariant* inv, + absl::string_view name = "LazySteepestSearch") + : BoolVectorBasedSolutionGenerator(inv, name) {} - bool NextSolution() final { - return NextSolution(SubsetBoolVector(num_subsets(), true)); - } - - bool NextSolution(absl::Span focus) final { - return NextSolution(MakeBoolVector(focus, num_subsets())); - } - - // The implementation works on a vector of Booleans. + using BoolVectorBasedSolutionGenerator::NextSolution; bool NextSolution(const SubsetBoolVector& in_focus) final; }; @@ -377,10 +427,11 @@ class TabuList { // 2 (1): 4–32. doi:10.1287/ijoc.2.1.4. // The consistency level is maintained up to kFreeAndUncovered. -class GuidedTabuSearch : public SetCoverSolutionGenerator { +class GuidedTabuSearch : public SubsetListBasedSolutionGenerator { public: - explicit GuidedTabuSearch(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv), + explicit GuidedTabuSearch(SetCoverInvariant* inv, + absl::string_view name = "GuidedTabuSearch") + : SubsetListBasedSolutionGenerator(inv, name), lagrangian_factor_(kDefaultLagrangianFactor), penalty_factor_(kDefaultPenaltyFactor), epsilon_(kDefaultEpsilon), @@ -393,15 +444,9 @@ class GuidedTabuSearch : public SetCoverSolutionGenerator { // Initializes the Guided Tabu Search algorithm. void Initialize(); - bool NextSolution() final { return NextSolution(model()->all_subsets()); } - - // The implementation works on a list of subsets. + using SubsetListBasedSolutionGenerator::NextSolution; bool NextSolution(absl::Span focus) final; - bool NextSolution(const SubsetBoolVector& in_focus) final { - return NextSolution(MakeSubsetIndexSpan(in_focus)); - } - // TODO(user): re-introduce this is the code. It was used to favor // subsets with the same marginal costs but that would cover more // elements. But first, see if it makes sense to compute it. @@ -470,10 +515,11 @@ class GuidedTabuSearch : public SetCoverSolutionGenerator { // Colchester, UK, July, 1997. // The consistency level is maintained up to kRedundancy. -class GuidedLocalSearch : public SetCoverSolutionGenerator { +class GuidedLocalSearch : public SubsetListBasedSolutionGenerator { public: - explicit GuidedLocalSearch(SetCoverInvariant* inv) - : SetCoverSolutionGenerator(inv), + explicit GuidedLocalSearch(SetCoverInvariant* inv, + absl::string_view name = "GuidedLocalSearch") + : SubsetListBasedSolutionGenerator(inv, name), epsilon_(kDefaultEpsilon), alpha_(kDefaultAlpha) { Initialize(); @@ -482,15 +528,9 @@ class GuidedLocalSearch : public SetCoverSolutionGenerator { // Initializes the Guided Local Search algorithm. void Initialize(); - bool NextSolution() final { return NextSolution(model()->all_subsets()); } - - // The implementation works on a list of subsets. + using SubsetListBasedSolutionGenerator::NextSolution; bool NextSolution(absl::Span focus) final; - bool NextSolution(const SubsetBoolVector& in_focus) final { - return NextSolution(MakeSubsetIndexSpan(in_focus)); - } - private: // Setters and getters for the Guided Local Search algorithm parameters. void SetEpsilon(double r) { epsilon_ = r; } diff --git a/ortools/set_cover/set_cover_invariant.cc b/ortools/set_cover/set_cover_invariant.cc index 36f7653dec..3deb63975d 100644 --- a/ortools/set_cover/set_cover_invariant.cc +++ b/ortools/set_cover/set_cover_invariant.cc @@ -42,6 +42,9 @@ void SetCoverInvariant::Initialize() { void SetCoverInvariant::Clear() { cost_ = 0.0; + lower_bound_ = 0.0; + is_cost_consistent_ = true; + const BaseInt num_subsets = model_->num_subsets(); const BaseInt num_elements = model_->num_elements(); diff --git a/ortools/set_cover/set_cover_invariant.h b/ortools/set_cover/set_cover_invariant.h index 37e0b8e338..4f0f0c425a 100644 --- a/ortools/set_cover/set_cover_invariant.h +++ b/ortools/set_cover/set_cover_invariant.h @@ -100,9 +100,32 @@ class SetCoverInvariant { const SetCoverModel* const_model() const { return model_; } + // Reports the lower bound of the problem. + // If is_cost_consistent is true, the lower bound is ignored and the cost is + // reported instead. + void ReportLowerBound(Cost lower_bound, bool is_cost_consistent) { + lower_bound_ = lower_bound; + is_cost_consistent_ = is_cost_consistent; + } + // Returns the cost of current solution. Cost cost() const { return cost_; } + // Returns the cost of the current solution, or the lower bound if the + // solution is a relaxation. + Cost CostOrLowerBound() const { + // For the time being we assume that we either have a feasible solution or + // a relaxation. + return is_cost_consistent_ ? cost_ : lower_bound_; + } + + // Returns true if the cost of the current solution is consistent with the + // solution. + bool is_cost_consistent() const { return is_cost_consistent_; } + + // Returns the lower bound of the problem. + Cost LowerBound() const { return lower_bound_; } + // Returns the number of uncovered elements. BaseInt num_uncovered_elements() const { return num_uncovered_elements_; } @@ -138,6 +161,19 @@ class SetCoverInvariant { // Clears the trace. void ClearTrace() { trace_.clear(); } + // Returns the cardinality of the solution in the invariant. + BaseInt ComputeCardinality() const { + BaseInt cardinality = 0; + for (BaseInt i = 0; i < trace().size(); ++i) { + const SubsetIndex subset = trace()[i].subset(); + DCHECK_EQ(trace()[i].decision(), is_selected_[subset]); + if (is_selected_[subset]) { + ++cardinality; + } + } + return cardinality; + } + // Clear the removability information. void ClearRemovabilityInformation() { newly_removable_subsets_.clear(); @@ -237,6 +273,16 @@ class SetCoverInvariant { // Current cost. Cost cost_; + // True it the cost is consistent with the solution. Used to decide whether to + // report the cost or the lower bound. + // This is initialized to true, and set to false when a lower bound is + // reported without a feasible solution. + bool is_cost_consistent_; + + // The lower bound of the problem, when has_relaxation_ is true. + // By default it is 0.0. + Cost lower_bound_; + // The number of uncovered (or free) elements in the current solution. BaseInt num_uncovered_elements_; diff --git a/ortools/set_cover/set_cover_lagrangian.cc b/ortools/set_cover/set_cover_lagrangian.cc index b01c36cffc..cd8393abaf 100644 --- a/ortools/set_cover/set_cover_lagrangian.cc +++ b/ortools/set_cover/set_cover_lagrangian.cc @@ -46,19 +46,19 @@ namespace operations_research { // // Concerning the fundamental ideas behind this approach, one may refer to [2]. ElementCostVector SetCoverLagrangian::InitializeLagrangeMultipliers() const { - ElementCostVector multipliers(model_.num_elements(), + ElementCostVector multipliers(model()->num_elements(), std::numeric_limits::infinity()); - SubsetCostVector marginal_costs(model_.num_subsets()); + SubsetCostVector marginal_costs(model()->num_subsets()); // TODO(user): Parallelize this. - for (const SubsetIndex subset : model_.SubsetRange()) { + for (const SubsetIndex subset : model()->SubsetRange()) { marginal_costs[subset] = - model_.subset_costs()[subset] / model_.columns()[subset].size(); + model()->subset_costs()[subset] / model()->columns()[subset].size(); } // TODO(user): Parallelize this. - for (const ElementIndex element : model_.ElementRange()) { + for (const ElementIndex element : model()->ElementRange()) { // Minimum marginal cost to cover this element. Cost min_marginal_cost = std::numeric_limits::infinity(); - const SparseRowView& rows = model_.rows(); + const SparseRowView& rows = model()->rows(); // TODO(user): use std::min_element on rows[element] with a custom // comparator that gets marginal_costs[subset]. Check performance. for (const SubsetIndex subset : rows[element]) { @@ -108,8 +108,8 @@ BaseInt BlockSize(BaseInt size, int num_threads) { // Computes the reduced costs for all subsets in parallel using ThreadPool. SubsetCostVector SetCoverLagrangian::ParallelComputeReducedCosts( const SubsetCostVector& costs, const ElementCostVector& multipliers) const { - const SubsetIndex num_subsets(model_.num_subsets()); - const SparseColumnView& columns = model_.columns(); + const SubsetIndex num_subsets(model()->num_subsets()); + const SparseColumnView& columns = model()->columns(); SubsetCostVector reduced_costs(num_subsets); // TODO(user): compute a close-to-optimal k-subset partitioning of the columns // based on their sizes. [***] @@ -138,7 +138,7 @@ SubsetCostVector SetCoverLagrangian::ParallelComputeReducedCosts( // c_j(u) = c_j - sum_{i \in I_j} a_{ij}.u_i SubsetCostVector SetCoverLagrangian::ComputeReducedCosts( const SubsetCostVector& costs, const ElementCostVector& multipliers) const { - const SparseColumnView& columns = model_.columns(); + const SparseColumnView& columns = model()->columns(); SubsetCostVector reduced_costs(costs.size()); FillReducedCostsSlice(SubsetIndex(0), SubsetIndex(reduced_costs.size()), costs, multipliers, columns, &reduced_costs); @@ -172,23 +172,23 @@ void FillSubgradientSlice(SubsetIndex slice_start, SubsetIndex slice_end, ElementCostVector SetCoverLagrangian::ComputeSubgradient( const SubsetCostVector& reduced_costs) const { // NOTE(user): Should the initialization be done with coverage[element]? - ElementCostVector subgradient(model_.num_elements(), 1.0); + ElementCostVector subgradient(model()->num_elements(), 1.0); FillSubgradientSlice(SubsetIndex(0), SubsetIndex(reduced_costs.size()), - model_.columns(), reduced_costs, &subgradient); + model()->columns(), reduced_costs, &subgradient); return subgradient; } ElementCostVector SetCoverLagrangian::ParallelComputeSubgradient( const SubsetCostVector& reduced_costs) const { - const SubsetIndex num_subsets(model_.num_subsets()); - const SparseColumnView& columns = model_.columns(); - ElementCostVector subgradient(model_.num_elements(), 1.0); + const SubsetIndex num_subsets(model()->num_subsets()); + const SparseColumnView& columns = model()->columns(); + ElementCostVector subgradient(model()->num_elements(), 1.0); // The subgradient has one component per element, each thread processes // several subsets. Hence, have one vector per thread to avoid data races. // TODO(user): it may be better to split the elements among the threads, // although this might be less well-balanced. std::vector subgradients( - num_threads_, ElementCostVector(model_.num_elements())); + num_threads_, ElementCostVector(model()->num_elements())); absl::BlockingCounter num_threads_running(num_threads_); const SubsetIndex block_size(BlockSize(num_subsets.value(), num_threads_)); SubsetIndex slice_start(0); @@ -206,7 +206,7 @@ ElementCostVector SetCoverLagrangian::ParallelComputeSubgradient( } num_threads_running.Wait(); for (int thread_index = 0; thread_index < num_threads_; ++thread_index) { - for (const ElementIndex element : model_.ElementRange()) { + for (const ElementIndex element : model()->ElementRange()) { subgradient[element] += subgradients[thread_index][element]; } } @@ -265,8 +265,8 @@ Cost SetCoverLagrangian::ParallelComputeLagrangianValue( } std::vector lagrangian_values(num_threads_, 0.0); absl::BlockingCounter num_threads_running(num_threads_); - const SubsetIndex block_size(BlockSize(model_.num_subsets(), num_threads_)); - const SubsetIndex num_subsets(model_.num_subsets()); + const SubsetIndex block_size(BlockSize(model()->num_subsets(), num_threads_)); + const SubsetIndex num_subsets(model()->num_subsets()); SubsetIndex slice_start(0); for (int thread_index = 0; thread_index < num_threads_; ++thread_index) { const SubsetIndex slice_end = @@ -318,7 +318,7 @@ void SetCoverLagrangian::UpdateMultipliers( // First compute lambda_k * (UB - L(u^k)). const Cost factor = step_size * (upper_bound - lagrangian_value) / subgradient_square_norm; - for (const ElementIndex element : model_.ElementRange()) { + for (const ElementIndex element : model()->ElementRange()) { // Avoid multipliers to go negative and to go over the roof. 1e6 chosen // arbitrarily. [***] (*multipliers)[element] = std::clamp( @@ -342,7 +342,7 @@ void SetCoverLagrangian::ParallelUpdateMultipliers( // First compute lambda_k * (UB - L(u^k)). const Cost factor = step_size * (upper_bound - lagrangian_value) / subgradient_square_norm; - for (const ElementIndex element : model_.ElementRange()) { + for (const ElementIndex element : model()->ElementRange()) { // Avoid multipliers to go negative and to go through the roof. 1e6 chosen const Cost kRoof = 1e6; // Arbitrary value, from [1]. (*multipliers)[element] = std::clamp( @@ -355,7 +355,7 @@ Cost SetCoverLagrangian::ComputeGap( const ElementCostVector& multipliers) const { Cost gap = 0.0; // TODO(user): Parallelize this, if need be. - for (const SubsetIndex subset : model_.SubsetRange()) { + for (const SubsetIndex subset : model()->SubsetRange()) { if (solution[subset] && reduced_costs[subset] > 0.0) { gap += reduced_costs[subset]; } else if (!solution[subset] && reduced_costs[subset] < 0.0) { @@ -363,8 +363,8 @@ Cost SetCoverLagrangian::ComputeGap( gap -= reduced_costs[subset]; } } - const ElementToIntVector& coverage = inv_->coverage(); - for (const ElementIndex element : model_.ElementRange()) { + const ElementToIntVector& coverage = inv()->coverage(); + for (const ElementIndex element : model()->ElementRange()) { gap += (coverage[element] - 1) * multipliers[element]; } return gap; @@ -373,12 +373,12 @@ Cost SetCoverLagrangian::ComputeGap( SubsetCostVector SetCoverLagrangian::ComputeDelta( const SubsetCostVector& reduced_costs, const ElementCostVector& multipliers) const { - SubsetCostVector delta(model_.num_subsets()); - const ElementToIntVector& coverage = inv_->coverage(); + SubsetCostVector delta(model()->num_subsets()); + const ElementToIntVector& coverage = inv()->coverage(); // This is definition (9) in [1]. - const SparseColumnView& columns = model_.columns(); + const SparseColumnView& columns = model()->columns(); // TODO(user): Parallelize this. - for (const SubsetIndex subset : model_.SubsetRange()) { + for (const SubsetIndex subset : model()->SubsetRange()) { delta[subset] = std::max(reduced_costs[subset], 0.0); for (const ElementIndex element : columns[subset]) { const double size = coverage[element]; @@ -491,6 +491,7 @@ class TopKHeap { std::tuple SetCoverLagrangian::ComputeLowerBound(const SubsetCostVector& costs, Cost upper_bound) { + StopWatch stop_watch(&run_time_); Cost lower_bound = 0.0; ElementCostVector multipliers = InitializeLagrangeMultipliers(); double step_size = 0.1; // [***] arbitrary, from [1]. @@ -515,6 +516,7 @@ SetCoverLagrangian::ComputeLowerBound(const SubsetCostVector& costs, // break; // } } + inv()->ReportLowerBound(lower_bound, /*is_cost_consistent=*/false); return std::make_tuple(lower_bound, reduced_costs, multipliers); } diff --git a/ortools/set_cover/set_cover_lagrangian.h b/ortools/set_cover/set_cover_lagrangian.h index 56b32dec14..c8893c711b 100644 --- a/ortools/set_cover/set_cover_lagrangian.h +++ b/ortools/set_cover/set_cover_lagrangian.h @@ -18,8 +18,10 @@ #include #include +#include "absl/types/span.h" #include "ortools/base/threadpool.h" #include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_heuristics.h" #include "ortools/set_cover/set_cover_invariant.h" #include "ortools/set_cover/set_cover_model.h" @@ -42,24 +44,28 @@ namespace operations_research { // Algorithms.” Mathematical Programming, 91 (3): 447–78. // https://link.springer.com/article/10.1007/s101070100262 -class SetCoverLagrangian { +// The implementation benefits from the features of SetCoverSolutionGenerator. +// Notice however that NextSolution() is not implemented, and the class is +// intended to be used only for ComputeLowerBound(). FOR THE TIME BEING. + +class SetCoverLagrangian : public SubsetListBasedSolutionGenerator { public: - explicit SetCoverLagrangian(SetCoverInvariant* inv, int num_threads = 1) - : inv_(inv), - model_(*inv->model()), - num_threads_(num_threads), - thread_pool_(new ThreadPool(num_threads)) { - thread_pool_->StartWorkers(); + explicit SetCoverLagrangian(SetCoverInvariant* inv, + const absl::string_view name = "Lagrangian") + : SubsetListBasedSolutionGenerator(inv, name), + num_threads_(1), + thread_pool_(nullptr) {} + + SetCoverLagrangian& UseNumThreads(int num_threads) { + num_threads_ = num_threads; + thread_pool_ = std::make_unique(num_threads); + return *this; } - // Returns true if a solution was found. - // TODO(user): Add time-outs and exit with a partial solution. This seems - // unlikely, though. - bool NextSolution(); + using SubsetListBasedSolutionGenerator::NextSolution; - // Computes the next partial solution considering only the subsets whose - // indices are in focus. - bool NextSolution(const std::vector& focus); + // This is a dummy implementation of NextSolution() that is not used. + bool NextSolution(absl::Span _) final { return false; } // Initializes the multipliers vector (u) based on the cost per subset. ElementCostVector InitializeLagrangeMultipliers() const; @@ -134,12 +140,6 @@ class SetCoverLagrangian { const SubsetCostVector& costs, Cost upper_bound); private: - // The invariant on which the algorithm will run. - SetCoverInvariant* inv_; - - // The model on which the invariant is defined. - const SetCoverModel& model_; - // The number of threads to use for parallelization. int num_threads_; diff --git a/ortools/set_cover/set_cover_mip.cc b/ortools/set_cover/set_cover_mip.cc index 1c7d346187..d0eff02b59 100644 --- a/ortools/set_cover/set_cover_mip.cc +++ b/ortools/set_cover/set_cover_mip.cc @@ -43,47 +43,40 @@ ElementToIntVector Subtract(const ElementToIntVector& a, template using StrictVector = glop::StrictITIVector; -bool SetCoverMip::NextSolution(bool use_integers, - double time_limit_in_seconds) { - return NextSolution(inv_->model()->all_subsets(), use_integers, - time_limit_in_seconds); -} - -bool SetCoverMip::NextSolution(absl::Span focus, - bool use_integers, - double time_limit_in_seconds) { - SetCoverModel* model = inv_->model(); - const SubsetIndex num_subsets(model->num_subsets()); - const ElementIndex num_elements(model->num_elements()); - SubsetBoolVector choices = inv_->is_selected(); +bool SetCoverMip::NextSolution(absl::Span focus) { + inv()->ReportLowerBound(0.0, true); + StopWatch stop_watch(&run_time_); + const SubsetIndex num_subsets(model()->num_subsets()); + const ElementIndex num_elements(model()->num_elements()); + SubsetBoolVector choices = inv()->is_selected(); MPSolver::OptimizationProblemType problem_type; switch (mip_solver_) { case SetCoverMipSolver::SCIP: problem_type = MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING; break; case SetCoverMipSolver::GUROBI: - if (use_integers) { + if (use_integers_) { problem_type = MPSolver::GUROBI_MIXED_INTEGER_PROGRAMMING; } else { problem_type = MPSolver::GUROBI_LINEAR_PROGRAMMING; } break; case SetCoverMipSolver::SAT: - if (!use_integers) { - LOG(INFO) << "Defaulting to integer variables with SAT"; - use_integers = true; + if (!use_integers_) { + VLOG(1) << "Defaulting to integer variables with SAT"; + use_integers_ = true; } problem_type = MPSolver::SAT_INTEGER_PROGRAMMING; break; case SetCoverMipSolver::GLOP: - LOG(INFO) << "Defaulting to linear relaxation with GLOP"; - use_integers = false; + VLOG(1) << "Defaulting to linear relaxation with GLOP"; + use_integers_ = false; problem_type = MPSolver::GLOP_LINEAR_PROGRAMMING; break; case SetCoverMipSolver::PDLP: - if (use_integers) { - LOG(INFO) << "Defaulting to linear relaxation with PDLP"; - use_integers = false; + if (use_integers_) { + VLOG(1) << "Defaulting to linear relaxation with PDLP"; + use_integers_ = false; } problem_type = MPSolver::PDLP_LINEAR_PROGRAMMING; break; @@ -104,13 +97,13 @@ bool SetCoverMip::NextSolution(absl::Span focus, StrictVector constraints(num_elements, nullptr); StrictVector vars(num_subsets, nullptr); ElementToIntVector coverage_outside_focus = - Subtract(inv_->coverage(), inv_->ComputeCoverageInFocus(focus)); + Subtract(inv()->coverage(), inv()->ComputeCoverageInFocus(focus)); for (const SubsetIndex subset : focus) { - vars[subset] = solver.MakeVar(0, 1, use_integers, ""); - objective->SetCoefficient(vars[subset], model->subset_costs()[subset]); - for (const ElementIndex element : model->columns()[subset]) { - // The model should only contain elements that are not forcibly covered by - // subsets outside the focus. + vars[subset] = solver.MakeVar(0, 1, use_integers_, ""); + objective->SetCoefficient(vars[subset], model()->subset_costs()[subset]); + for (const ElementIndex element : model()->columns()[subset]) { + // The model should only contain elements that are not forcibly covered + // by subsets outside the focus. if (coverage_outside_focus[element] != 0) continue; if (constraints[element] == nullptr) { constexpr double kInfinity = std::numeric_limits::infinity(); @@ -120,7 +113,7 @@ bool SetCoverMip::NextSolution(absl::Span focus, } } // set_time_limit takes milliseconds as a unit. - solver.set_time_limit(static_cast(time_limit_in_seconds * 1000)); + solver.set_time_limit(static_cast(time_limit_in_seconds() * 1000)); // Call the solver. const MPSolver::ResultStatus solve_status = solver.Solve(); @@ -139,15 +132,17 @@ bool SetCoverMip::NextSolution(absl::Span focus, LOG(ERROR) << "Solving resulted in an error."; return false; } - if (use_integers) { + if (use_integers_) { using CL = SetCoverInvariant::ConsistencyLevel; for (const SubsetIndex subset : focus) { if (vars[subset]->solution_value() > 0.9) { - inv_->Select(subset, CL::kCostAndCoverage); + inv()->Select(subset, CL::kCostAndCoverage); } } } else { - lower_bound_ = solver.Objective().Value(); + // Report the objective value as a lower bound, and mention that the cost is + // not consistent with the solution. + inv()->ReportLowerBound(solver.Objective().Value(), false); } return true; } diff --git a/ortools/set_cover/set_cover_mip.h b/ortools/set_cover/set_cover_mip.h index f06c7613ed..53a8de1341 100644 --- a/ortools/set_cover/set_cover_mip.h +++ b/ortools/set_cover/set_cover_mip.h @@ -14,8 +14,10 @@ #ifndef OR_TOOLS_SET_COVER_SET_COVER_MIP_H_ #define OR_TOOLS_SET_COVER_SET_COVER_MIP_H_ +#include "absl/strings/string_view.h" #include "absl/types/span.h" #include "ortools/set_cover/base_types.h" +#include "ortools/set_cover/set_cover_heuristics.h" #include "ortools/set_cover/set_cover_invariant.h" namespace operations_research { @@ -27,43 +29,37 @@ enum class SetCoverMipSolver : int { PDLP = 4 }; -class SetCoverMip { +class SetCoverMip : public SubsetListBasedSolutionGenerator { public: // Simpler constructor that uses SCIP by default. - explicit SetCoverMip(SetCoverInvariant* inv) - : inv_(inv), mip_solver_(SetCoverMipSolver::SCIP) {} + explicit SetCoverMip(SetCoverInvariant* inv, + const absl::string_view name = "Mip") + : SubsetListBasedSolutionGenerator(inv, name), + mip_solver_(SetCoverMipSolver::SCIP), + use_integers_(true) {} - // The constructor takes a SetCoverInvariant that will store the resulting - // variable choices, and a MIP Solver. - SetCoverMip(SetCoverInvariant* inv, SetCoverMipSolver mip_solver) - : inv_(inv), mip_solver_(mip_solver) {} + SetCoverMip& UseMipSolver(SetCoverMipSolver mip_solver) { + mip_solver_ = mip_solver; + return *this; + } - // Returns true if a solution was found. - // If use_integers is false, lower_bound_ is populated with a linear - // lower bound. - // time_limit_in_seconds is a (rather soft) time limit for the execution time. - // TODO(user): Add time-outs and exit with a partial solution. This seems - // unlikely, though. - bool NextSolution(bool use_integers, double time_limit_in_seconds); + SetCoverMip& UseIntegers(bool use_integers) { + use_integers_ = use_integers; + return *this; + } + + using SubsetListBasedSolutionGenerator::NextSolution; // Computes the next partial solution considering only the subsets whose // indices are in focus. - bool NextSolution(absl::Span focus, bool use_integers, - double time_limit_in_seconds); - - // Returns the lower bound of the linear relaxation of the problem. - double lower_bound() const { return lower_bound_; } + bool NextSolution(absl::Span focus) final; private: - // The invariant used to maintain the state of the problem. - SetCoverInvariant* inv_; - // The MIP solver flavor used by the instance. SetCoverMipSolver mip_solver_; - // The lower bound of the problem, when use_integers is false. The MIP with - // continuous variables becomes a computationally simpler linear program. - double lower_bound_; + // Whether to use integer variables in the MIP. + bool use_integers_; }; } // namespace operations_research diff --git a/ortools/set_cover/set_cover_model.cc b/ortools/set_cover/set_cover_model.cc index c86f5cfe55..cd81873fc6 100644 --- a/ortools/set_cover/set_cover_model.cc +++ b/ortools/set_cover/set_cover_model.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -30,6 +31,8 @@ #include "absl/random/distributions.h" #include "absl/random/random.h" #include "absl/strings/str_format.h" +#include "absl/strings/str_join.h" +#include "absl/strings/string_view.h" #include "absl/types/span.h" #include "ortools/algorithms/radix_sort.h" #include "ortools/set_cover/base_types.h" @@ -88,12 +91,13 @@ SetCoverModel SetCoverModel::GenerateRandomModelFrom( const SetCoverModel& seed_model, BaseInt num_elements, BaseInt num_subsets, double row_scale, double column_scale, double cost_scale) { SetCoverModel model; + StopWatch stop_watch(&(model.generation_duration_)); DCHECK_GT(row_scale, 0.0); DCHECK_GT(column_scale, 0.0); DCHECK_GT(cost_scale, 0.0); model.num_elements_ = num_elements; model.num_nonzeros_ = 0; - model.ReserveNumSubsets(num_subsets); + model.ResizeNumSubsets(num_subsets); model.UpdateAllSubsetsList(); absl::BitGen bitgen; @@ -304,30 +308,23 @@ void SetCoverModel::AddElementToSubset(ElementIndex element, AddElementToSubset(element.value(), subset.value()); } -// Reserves num_subsets columns in the model. -void SetCoverModel::ReserveNumSubsets(BaseInt num_subsets) { +void SetCoverModel::ResizeNumSubsets(BaseInt num_subsets) { num_subsets_ = std::max(num_subsets_, num_subsets); columns_.resize(num_subsets_, SparseColumn()); subset_costs_.resize(num_subsets_, 0.0); UpdateAllSubsetsList(); } -void SetCoverModel::ReserveNumSubsets(SubsetIndex num_subsets) { - ReserveNumSubsets(num_subsets.value()); +void SetCoverModel::ResizeNumSubsets(SubsetIndex num_subsets) { + ResizeNumSubsets(num_subsets.value()); } -// Reserves num_elements rows in the column indexed by subset. void SetCoverModel::ReserveNumElementsInSubset(BaseInt num_elements, BaseInt subset) { - ReserveNumSubsets(subset); + ResizeNumSubsets(subset); columns_[SubsetIndex(subset)].reserve(ColumnEntryIndex(num_elements)); } -void SetCoverModel::ReserveNumElementsInSubset(ElementIndex num_elements, - SubsetIndex subset) { - ReserveNumElementsInSubset(num_elements.value(), subset.value()); -} - void SetCoverModel::SortElementsInSubsets() { for (const SubsetIndex subset : SubsetRange()) { // std::sort(columns_[subset].begin(), columns_[subset].end()); @@ -338,6 +335,7 @@ void SetCoverModel::SortElementsInSubsets() { } void SetCoverModel::CreateSparseRowView() { + StopWatch stop_watch(&create_sparse_row_view_duration_); if (row_view_is_valid_) { VLOG(1) << "CreateSparseRowView: already valid"; return; @@ -371,7 +369,108 @@ void SetCoverModel::CreateSparseRowView() { VLOG(1) << "CreateSparseRowView finished"; } -bool SetCoverModel::ComputeFeasibility() const { +std::vector SetCoverModel::ComputeSliceIndices( + int num_partitions) { + if (num_partitions <= 1 || columns_.empty()) { + return {SubsetIndex(columns_.size())}; + } + + std::vector partial_sum_nnz(columns_.size()); + partial_sum_nnz[0] = columns_[SubsetIndex(0)].size(); + for (BaseInt col = 1; col < columns_.size(); ++col) { + partial_sum_nnz[col] = + partial_sum_nnz[col - 1] + columns_[SubsetIndex(col)].size(); + } + const BaseInt total_nnz = partial_sum_nnz.back(); + const BaseInt target_nnz = (total_nnz + num_partitions - 1) / num_partitions; + + std::vector partition_points(num_partitions); + BaseInt threshold = target_nnz; + BaseInt pos = 0; + for (const SubsetIndex col : SubsetRange()) { + if (partial_sum_nnz[col.value()] >= threshold) { + partition_points[pos] = col; + ++pos; + threshold += target_nnz; + } + } + partition_points[num_partitions - 1] = SubsetIndex(columns_.size()); + return partition_points; +} + +SparseRowView SetCoverModel::ComputeSparseRowViewSlice(SubsetIndex begin_subset, + SubsetIndex end_subset) { + SparseRowView rows; + rows.reserve(num_elements_); + ElementToIntVector row_sizes(num_elements_, 0); + for (SubsetIndex subset = begin_subset; subset < end_subset; ++subset) { + BaseInt* data = reinterpret_cast(columns_[subset].data()); + RadixSort(absl::MakeSpan(data, columns_[subset].size())); + + ElementIndex preceding_element(-1); + for (const ElementIndex element : columns_[subset]) { + CHECK_GT(element, preceding_element) + << "Repetition in column " + << subset; // Fail if there is a repetition. + ++row_sizes[element]; + preceding_element = element; + } + } + for (const ElementIndex element : ElementRange()) { + rows[element].reserve(RowEntryIndex(row_sizes[element])); + } + for (SubsetIndex subset = begin_subset; subset < end_subset; ++subset) { + for (const ElementIndex element : columns_[subset]) { + rows[element].emplace_back(subset); + } + } + return rows; +} + +std::vector SetCoverModel::CutSparseRowViewInSlices( + const std::vector& partition_points) { + std::vector row_views; + row_views.reserve(partition_points.size()); + SubsetIndex begin_subset(0); + // This should be done in parallel. It is a bottleneck. + for (SubsetIndex end_subset : partition_points) { + row_views.push_back(ComputeSparseRowViewSlice(begin_subset, end_subset)); + begin_subset = end_subset; + } + return row_views; +} + +SparseRowView SetCoverModel::ReduceSparseRowViewSlices( + const std::vector& slices) { + SparseRowView result_rows; + // This is not a ReduceTree. This will be done later through parallelism. + result_rows.reserve(num_elements_); + ElementToIntVector row_sizes(num_elements_, 0); + for (const SparseRowView& slice : slices) { + // This should done as a reduce tree, in parallel. + for (const ElementIndex element : ElementRange()) { + result_rows[element].insert(result_rows[element].end(), + slice[element].begin(), slice[element].end()); + } + } + return result_rows; +} + +SparseRowView SetCoverModel::ComputeSparseRowViewUsingSlices() { + StopWatch stop_watch(&compute_sparse_row_view_using_slices_duration_); + SparseRowView rows; + VLOG(1) << "CreateSparseRowViewUsingSlices started"; + const std::vector partition_points = + ComputeSliceIndices(num_subsets()); + const std::vector slices = + CutSparseRowViewInSlices(partition_points); + rows = ReduceSparseRowViewSlices(slices); + VLOG(1) << "CreateSparseRowViewUsingSlices finished"; + return rows; +} + +bool SetCoverModel::ComputeFeasibility() { + StopWatch stop_watch(&feasibility_duration_); CHECK_GT(num_elements(), 0); CHECK_GT(num_subsets(), 0); CHECK_EQ(columns_.size(), num_subsets()); @@ -438,7 +537,7 @@ SetCoverProto SetCoverModel::ExportModelAsProto() const { void SetCoverModel::ImportModelFromProto(const SetCoverProto& message) { columns_.clear(); subset_costs_.clear(); - ReserveNumSubsets(message.subset_size()); + ResizeNumSubsets(message.subset_size()); SubsetIndex subset_index(0); for (const SetCoverProto::Subset& subset_proto : message.subset()) { subset_costs_[subset_index] = subset_proto.cost(); @@ -456,6 +555,20 @@ void SetCoverModel::ImportModelFromProto(const SetCoverProto& message) { CreateSparseRowView(); } +std::string SetCoverModel::ToVerboseString(absl::string_view sep) const { + return absl::StrJoin( + std::make_tuple("num_elements", num_elements(), "num_subsets", + num_subsets(), "num_nonzeros", num_nonzeros(), + "fill_rate", FillRate()), + sep); +} + +std::string SetCoverModel::ToString(absl::string_view sep) const { + return absl::StrJoin(std::make_tuple(num_elements(), num_subsets(), + num_nonzeros(), FillRate()), + sep); +} + namespace { // Returns the standard deviation of the vector v, excluding those values that // are zero. @@ -580,6 +693,18 @@ SetCoverModel::Stats SetCoverModel::ComputeColumnStats() const { return ComputeStats(std::move(column_sizes)); } +std::string SetCoverModel::Stats::ToString(absl::string_view sep) const { + return absl::StrJoin(std::make_tuple(min, max, median, mean, stddev, iqr), + sep); +} + +std::string SetCoverModel::Stats::ToVerboseString(absl::string_view sep) const { + return absl::StrJoin( + std::make_tuple("min", min, "max", max, "median", median, "mean", mean, + "stddev", stddev, "iqr", iqr), + sep); +} + std::vector SetCoverModel::ComputeRowDeciles() const { std::vector row_sizes(num_elements(), 0); for (const SparseColumn& column : columns_) { @@ -631,5 +756,4 @@ SetCoverModel::Stats SetCoverModel::ComputeRowDeltaSizeStats() const { } return acc.ComputeStats(); } - } // namespace operations_research diff --git a/ortools/set_cover/set_cover_model.h b/ortools/set_cover/set_cover_model.h index b49410b58f..54ae32946c 100644 --- a/ortools/set_cover/set_cover_model.h +++ b/ortools/set_cover/set_cover_model.h @@ -19,7 +19,7 @@ #include #include "absl/log/check.h" -#include "absl/strings/str_cat.h" +#include "absl/strings/string_view.h" #include "ortools/base/strong_int.h" #include "ortools/base/strong_vector.h" #include "ortools/set_cover/base_types.h" @@ -56,8 +56,9 @@ namespace operations_research { class SetCoverModel { public: // Constructs an empty weighted set-covering problem. - SetCoverModel() - : num_elements_(0), + explicit SetCoverModel(const absl::string_view name = "SetCoverModel") + : name_(name), + num_elements_(0), num_subsets_(0), num_nonzeros_(0), row_view_is_valid_(false), @@ -69,6 +70,10 @@ class SetCoverModel { rows_(), all_subsets_() {} + std::string name() const { return name_; } + + void SetName(const absl::string_view name) { name_ = name; } + // Constructs a weighted set-covering problem from a seed model, with // num_elements elements and num_subsets subsets. // - The distributions of the degrees of the elements and the cardinalities of @@ -103,15 +108,43 @@ class SetCoverModel { // Current number of subsets in the model. In matrix terms, this is the // number of columns. - BaseInt num_subsets() const { return num_subsets_; } + BaseInt num_subsets() const { + CHECK_NE(this, nullptr); + return num_subsets_; + } // Current number of nonzeros in the matrix. int64_t num_nonzeros() const { return num_nonzeros_; } + // Returns the fill rate of the matrix. double FillRate() const { return 1.0 * num_nonzeros() / (1.0 * num_elements() * num_subsets()); } + // Computes the number of singleton columns in the model, i.e. subsets + // covering only one element. + BaseInt ComputeNumSingletonColumns() const { + BaseInt num_singleton_columns = 0; + for (const SparseColumn& column : columns_) { + if (column.size() == 1) { + ++num_singleton_columns; + } + } + return num_singleton_columns; + } + + // Computes the number of singleton rows in the model, i.e. elements in the + // model that can be covered by one subset only. + BaseInt ComputeNumSingletonRows() const { + BaseInt num_singleton_rows = 0; + for (const SparseRow& row : rows_) { + if (row.size() == 1) { + ++num_singleton_rows; + } + } + return num_singleton_rows; + } + // Vector of costs for each subset. const SubsetCostVector& subset_costs() const { return subset_costs_; } @@ -157,6 +190,17 @@ class SetCoverModel { // Returns the list of indices for all the subsets in the model. std::vector all_subsets() const { return all_subsets_; } + // Accessors to the durations of the different stages of the model creation. + absl::Duration generation_duration() const { return generation_duration_; } + + absl::Duration create_sparse_row_view_duration() const { + return create_sparse_row_view_duration_; + } + + absl::Duration compute_sparse_row_view_using_slices_duration() const { + return compute_sparse_row_view_using_slices_duration_; + } + // Adds an empty subset with a cost to the problem. In matrix terms, this // adds a column to the matrix. void AddEmptySubset(Cost cost); @@ -184,18 +228,50 @@ class SetCoverModel { // the elements in each subset. void CreateSparseRowView(); - // Returns true if the problem is feasible, i.e. if the subsets cover all - // the elements. - bool ComputeFeasibility() const; + // Same as CreateSparseRowView, but uses a slicing algorithm, more prone to + // parallelism. + SparseRowView ComputeSparseRowViewUsingSlices(); - // Reserves num_subsets columns in the model. - void ReserveNumSubsets(BaseInt num_subsets); - void ReserveNumSubsets(SubsetIndex num_subsets); + // The following functions are used by CreateSparseRowViewUsingSlices. + // They are exposed for testing purposes. + + // Returns a vector of subset indices that partition columns into + // `num_partitions` partitions of roughly equal size in number of non-zeros. + // The resulting vector contains the index of the first subset in the next + // partition. Therefore the last element is the total number of subsets. + // Also the vector is sorted. + std::vector ComputeSliceIndices(int num_partitions); + + // Returns a view of the rows of the problem with subset in the range + // [begin_subset, end_subset). + SparseRowView ComputeSparseRowViewSlice(SubsetIndex begin_subset, + SubsetIndex end_subset); + + // Returns a vector of row views, each corresponding to a partition of the + // problem. The partitions are defined by the given partition points. + std::vector CutSparseRowViewInSlices( + const std::vector& partition_points); + + // Returns the union of the rows of the given row views. + // The returned view is valid only as long as the given row views are valid. + // The indices in the rows are sorted. + SparseRowView ReduceSparseRowViewSlices( + const std::vector& row_slices); + + // Returns true if the problem is feasible, i.e. if the subsets cover all + // the elements. Could be const, but it updates the feasibility_duration_ + // field. + bool ComputeFeasibility(); + + // Resizes the model to have at least num_subsets columns. The new columns + // correspond to empty subsets. + // This function will never remove a column, i.e. if num_subsets is smaller + // than the current instance size the function does nothing. + void ResizeNumSubsets(BaseInt num_subsets); + void ResizeNumSubsets(SubsetIndex num_subsets); // Reserves num_elements rows in the column indexed by subset. void ReserveNumElementsInSubset(BaseInt num_elements, BaseInt subset); - void ReserveNumElementsInSubset(ElementIndex num_elements, - SubsetIndex subset); // Returns the model as a SetCoverProto. Note that the elements of each subset // are sorted locally before being exported to the proto. This is done to @@ -207,6 +283,13 @@ class SetCoverModel { // Imports the model from a SetCoverProto. void ImportModelFromProto(const SetCoverProto& message); + // Returns a verbose string representation of the model, using the given + // separator. + std::string ToVerboseString(absl::string_view sep) const; + + // Returns a string representation of the model, using the given separator. + std::string ToString(absl::string_view sep) const; + // A struct enabling to show basic statistics on rows and columns. // The meaning of the fields is obvious. struct Stats { @@ -217,11 +300,16 @@ class SetCoverModel { double stddev; double iqr; // Interquartile range. - std::string DebugString() const { - return absl::StrCat("min = ", min, ", max = ", max, ", mean = ", mean, - ", median = ", median, ", stddev = ", stddev, ", ", - "iqr = ", iqr); - } + // Returns a string representation of the stats. TODO(user): remove this as + // it is obsolete. + std::string DebugString() const { return ToVerboseString(", "); } + + // Returns a string representation of the stats, using the given separator. + std::string ToString(absl::string_view sep) const; + + // Returns a verbose string representation of the stats, using the given + // separator. + std::string ToVerboseString(absl::string_view sep) const; }; // Computes basic statistics on costs and returns a Stats structure. @@ -252,6 +340,9 @@ class SetCoverModel { // columns.size() - 1 void UpdateAllSubsetsList(); + // The name of the model, "SetCoverModel" as default. + std::string name_; + // Number of elements. BaseInt num_elements_; @@ -277,6 +368,20 @@ class SetCoverModel { // True when is_unicost_ is up-to-date. bool is_unicost_valid_; + // Stores the run time for CreateSparseRowView. Interesting to compare with + // the time spent to actually generate a solution to the model. + absl::Duration create_sparse_row_view_duration_; + + // Stores the run time for ComputeSparseRowViewUsingSlices. + absl::Duration compute_sparse_row_view_using_slices_duration_; + + // Stores the run time for GenerateRandomModelFrom. + absl::Duration generation_duration_; + + // Stores the run time for ComputeFeasibility. A good estimate of the time + // needed to traverse the complete model. + absl::Duration feasibility_duration_; + // Vector of columns. Each column corresponds to a subset and contains the // elements of the given subset. // This takes NNZ (number of non-zeros) BaseInts, or |E| * |S| * fill_rate. @@ -375,7 +480,6 @@ class IntersectingSubsetsIterator { // already seen by the iterator. SubsetBoolVector subset_seen_; }; - } // namespace operations_research #endif // OR_TOOLS_SET_COVER_SET_COVER_MODEL_H_ diff --git a/ortools/set_cover/set_cover_reader.cc b/ortools/set_cover/set_cover_reader.cc index 284ee5a503..6ac534c6fd 100644 --- a/ortools/set_cover/set_cover_reader.cc +++ b/ortools/set_cover/set_cover_reader.cc @@ -126,7 +126,7 @@ SetCoverModel ReadOrlibScp(absl::string_view filename) { SetCoverReader reader(file); const ElementIndex num_rows(reader.ParseNextInteger()); const SubsetIndex num_cols(reader.ParseNextInteger()); - model.ReserveNumSubsets(num_cols); + model.ResizeNumSubsets(num_cols); for (SubsetIndex subset : SubsetRange(num_cols)) { const double cost(reader.ParseNextDouble()); model.SetSubsetCost(subset, cost); @@ -157,7 +157,7 @@ SetCoverModel ReadOrlibRail(absl::string_view filename) { SetCoverReader reader(file); const ElementIndex num_rows(reader.ParseNextInteger()); const SubsetIndex num_cols(reader.ParseNextInteger()); - model.ReserveNumSubsets(num_cols); + model.ResizeNumSubsets(num_cols); for (SubsetIndex subset : SubsetRange(num_cols)) { LOG_EVERY_N_SEC(INFO, 5) << absl::StrFormat("Reading subset %d (%.1f%%)", subset.value(), diff --git a/ortools/set_cover/set_cover_solve.cc b/ortools/set_cover/set_cover_solve.cc index 75356ea98e..eccccea78a 100644 --- a/ortools/set_cover/set_cover_solve.cc +++ b/ortools/set_cover/set_cover_solve.cc @@ -11,13 +11,21 @@ // See the License for the specific language governing permissions and // limitations under the License. +#include +#include +#include +#include #include +#include #include +#include "absl/base/macros.h" #include "absl/flags/flag.h" #include "absl/log/check.h" #include "absl/log/log.h" #include "absl/strings/match.h" +#include "absl/strings/str_cat.h" +#include "absl/strings/str_format.h" #include "absl/strings/str_join.h" #include "absl/strings/string_view.h" #include "absl/time/time.h" @@ -29,20 +37,39 @@ #include "ortools/set_cover/set_cover_model.h" #include "ortools/set_cover/set_cover_reader.h" +// Example usages: +// +// Solve all the problems in the benchmarks directory and produce LaTeX output. +// Run the classic algorithms on problem with up to 100,000 elements. Display +// summaries (with geomean ratios) for each group of problems. +/* Copy-paste to a terminal: + set_cover_solve --benchmarks --benchmarks_dir ~/benchmarks \ + --max_elements_for_classic 100000 --solve --latex --summarize +*/ +// Generate a new model from the rail4284 problem, with 100,000 elements and +// 1,000,000,000 subsets, with row_scale = 1.1, column_scale = 1.1, and +// cost_scale = 10.0: +/* Copy-paste to a terminal: + set_cover_solve --input ~/scp-orlib/rail4284.txt --input_fmt rail \ + --output ~/rail4284_1B.txt --output_fmt orlibrail \ + --num_elements_wanted 100000 --num_subsets_wanted 100000000 \ + --cost_scale 10.0 --row_scale 1.1 --column_scale 1.1 --generate +*/ +// Display statistics about trail4284_1B.txt: +/* Copy-paste to a terminal: + set_cover_solve --input ~/rail4284_1B.txt --input_fmt orlib --stats +*/ +// + 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."); + "REQUIRED: Input file format. Either proto, proto_bin, rail, " + "orlib or fimi."); 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(bool, generate, false, "Generate a new model from the input model."); ABSL_FLAG(int, num_elements_wanted, 0, @@ -54,84 +81,102 @@ 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."); +ABSL_FLAG(bool, unicost, false, "Set all costs to 1.0."); + +ABSL_FLAG(bool, latex, false, "Output in LaTeX format. CSV otherwise."); ABSL_FLAG(bool, solve, false, "Solve the model."); ABSL_FLAG(bool, stats, false, "Log stats about the model."); +ABSL_FLAG(bool, summarize, false, + "Display the comparison of the solution generators."); +ABSL_FLAG(bool, tlns, false, "Run thrifty LNS."); +ABSL_FLAG(int, max_elements_for_classic, 5000, + "Do not use classic on larger problems."); + +ABSL_FLAG(bool, benchmarks, false, "Run benchmarks."); +ABSL_FLAG(std::string, benchmarks_dir, "", "Benchmarks directory."); + +// TODO(user): Add a flags to: +// - Choose problems by name or by size: filter_name, max_elements, max_subsets. +// - Exclude problems by name: exclude_name. +// - Choose which solution generators to run. +// - Parameterize the number of threads. num_threads. 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(); +int64_t RunTimeInMicroseconds(const SetCoverSolutionGenerator& gen) { + return absl::ToInt64Microseconds(gen.run_time()); +} - 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; - } +int64_t RunTimeInNanoseconds(const SetCoverSolutionGenerator& gen) { + return absl::ToInt64Nanoseconds(gen.run_time()); +} + +void LogStats(const SetCoverModel& model) { + const std::string sep = absl::GetFlag(FLAGS_latex) ? " & " : ", "; + std::string header = absl::StrCat(model.name(), sep); + // Lines start with a comma to make it easy to copy-paste the output to a + // spreadsheet as CSV. + if (!absl::GetFlag(FLAGS_latex)) { + header = absl::StrCat(sep, header); } - 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; + LOG(INFO) << header << model.ToVerboseString(sep); + LOG(INFO) << header << "cost" << sep + << model.ComputeCostStats().ToVerboseString(sep); + LOG(INFO) << header << "row stats" << sep + << model.ComputeRowStats().ToVerboseString(sep); + LOG(INFO) << header << "row size deciles" << sep + << absl::StrJoin(model.ComputeRowDeciles(), sep); + LOG(INFO) << header << "row delta byte size stats" << sep + << model.ComputeRowDeltaSizeStats().ToVerboseString(sep); + LOG(INFO) << header << "column stats" << sep + << model.ComputeColumnStats().ToVerboseString(sep); + LOG(INFO) << header << "column size deciles" << sep + << absl::StrJoin(model.ComputeColumnDeciles(), sep); + LOG(INFO) << header << "column delta byte size stats" << sep + << model.ComputeColumnDeltaSizeStats().ToVerboseString(sep); + LOG(INFO) << header << "num_singleton_rows" << sep + << model.ComputeNumSingletonRows() << sep << "num_singleton_columns" + << sep << model.ComputeNumSingletonColumns(); } -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(const absl::string_view problem_name, + absl::string_view alg_name, const SetCoverInvariant& inv, + int64_t run_time) { + LOG(INFO) << ", " << problem_name << ", " << alg_name << ", cost, " + << inv.CostOrLowerBound() << ", solution_cardinality, " + << inv.ComputeCardinality() << ", " << run_time << "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); +void LogCostAndTiming(const SetCoverSolutionGenerator& generator) { + const SetCoverInvariant& inv = *generator.inv(); + const SetCoverModel& model = *inv.model(); + LogCostAndTiming(model.name(), generator.name(), inv, + RunTimeInMicroseconds(generator)); } +// TODO(user): Move this set_cover_reader enum class FileFormat { EMPTY, - ORLIB_SCP, - ORLIB_RAIL, - FIMI_DAT, + ORLIB, + RAIL, + FIMI, PROTO, PROTO_BIN, TXT, }; +// TODO(user): Move this set_cover_reader 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, "orlib")) { + return FileFormat::ORLIB; + } else if (absl::EqualsIgnoreCase(format_name, "rail")) { + return FileFormat::RAIL; + } else if (absl::EqualsIgnoreCase(format_name, "fimi")) { + return FileFormat::FIMI; } else if (absl::EqualsIgnoreCase(format_name, "proto")) { return FileFormat::PROTO; } else if (absl::EqualsIgnoreCase(format_name, "proto_bin")) { @@ -143,107 +188,465 @@ FileFormat ParseFileFormat(const std::string& format_name) { } } -SetCoverModel ReadModel(absl::string_view 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); +// TODO(user): Move this set_cover_reader +SetCoverModel ReadModel(absl::string_view filename, FileFormat format) { + switch (format) { + case FileFormat::ORLIB: + return ReadOrlibScp(filename); + case FileFormat::RAIL: + return ReadOrlibRail(filename); + case FileFormat::FIMI: + return ReadFimiDat(filename); case FileFormat::PROTO: - return ReadSetCoverProto(input_file, /*binary=*/false); + return ReadSetCoverProto(filename, /*binary=*/false); case FileFormat::PROTO_BIN: - return ReadSetCoverProto(input_file, /*binary=*/true); + return ReadSetCoverProto(filename, /*binary=*/true); default: - LOG(FATAL) << "Unsupported input format: " - << static_cast(input_format); + LOG(FATAL) << "Unsupported input format: " << static_cast(format); } } -SubsetBoolVector ReadSolution(absl::string_view input_file, - FileFormat input_format) { - switch (input_format) { +// TODO(user): Move this set_cover_reader +SubsetBoolVector ReadSolution(absl::string_view filename, FileFormat format) { + switch (format) { case FileFormat::TXT: - return ReadSetCoverSolutionText(input_file); + return ReadSetCoverSolutionText(filename); case FileFormat::PROTO: - return ReadSetCoverSolutionProto(input_file, /*binary=*/false); + return ReadSetCoverSolutionProto(filename, /*binary=*/false); case FileFormat::PROTO_BIN: - return ReadSetCoverSolutionProto(input_file, /*binary=*/true); + return ReadSetCoverSolutionProto(filename, /*binary=*/true); default: - LOG(FATAL) << "Unsupported input format: " - << static_cast(input_format); + LOG(FATAL) << "Unsupported input format: " << static_cast(format); } } -void WriteModel(const SetCoverModel& model, const std::string& output_file, - FileFormat output_format) { - LOG(INFO) << "Writing model to " << output_file; - switch (output_format) { - case FileFormat::ORLIB_SCP: - WriteOrlibScp(model, output_file); +// TODO(user): Move this set_cover_reader +void WriteModel(const SetCoverModel& model, const std::string& filename, + FileFormat format) { + LOG(INFO) << "Writing model to " << filename; + switch (format) { + case FileFormat::ORLIB: + WriteOrlibScp(model, filename); break; - case FileFormat::ORLIB_RAIL: - WriteOrlibRail(model, output_file); + case FileFormat::RAIL: + WriteOrlibRail(model, filename); break; case FileFormat::PROTO: - WriteSetCoverProto(model, output_file, /*binary=*/false); + WriteSetCoverProto(model, filename, /*binary=*/false); break; case FileFormat::PROTO_BIN: - WriteSetCoverProto(model, output_file, /*binary=*/true); + WriteSetCoverProto(model, filename, /*binary=*/true); break; default: - LOG(FATAL) << "Unsupported output format: " - << static_cast(output_format); + LOG(FATAL) << "Unsupported output format: " << static_cast(format); } } +// TODO(user): Move this set_cover_reader void WriteSolution(const SetCoverModel& model, const SubsetBoolVector& solution, - absl::string_view output_file, FileFormat output_format) { - switch (output_format) { + absl::string_view filename, FileFormat format) { + switch (format) { case FileFormat::TXT: - WriteSetCoverSolutionText(model, solution, output_file); + WriteSetCoverSolutionText(model, solution, filename); break; case FileFormat::PROTO: - WriteSetCoverSolutionProto(model, solution, output_file, + WriteSetCoverSolutionProto(model, solution, filename, /*binary=*/false); break; case FileFormat::PROTO_BIN: - WriteSetCoverSolutionProto(model, solution, output_file, + WriteSetCoverSolutionProto(model, solution, filename, /*binary=*/true); break; default: - LOG(FATAL) << "Unsupported output format: " - << static_cast(output_format); + LOG(FATAL) << "Unsupported output format: " << static_cast(format); } } -SetCoverInvariant RunLazyElementDegree(std::string name, SetCoverModel* model) { +SetCoverInvariant RunLazyElementDegree(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); + LogCostAndTiming(element_degree); return inv; } +SetCoverInvariant RunGreedy(SetCoverModel* model) { + SetCoverInvariant inv(model); + GreedySolutionGenerator classic(&inv); + CHECK(classic.NextSolution()); + DCHECK(inv.CheckConsistency(CL::kCostAndCoverage)); + LogCostAndTiming(classic); + return inv; +} + +// Formatting and reporting functions with LaTeX and CSV support. +namespace { +std::string Separator() { return absl::GetFlag(FLAGS_latex) ? " & " : ", "; } + +std::string Eol() { return absl::GetFlag(FLAGS_latex) ? " \\\\\n" : "\n"; } + +// Computes the ratio of the cost of the new solution generator to the cost of +// the reference solution generator. +double CostRatio(SetCoverSolutionGenerator& ref_gen, + SetCoverSolutionGenerator& new_gen) { + return new_gen.cost() / ref_gen.cost(); +} + +// Computes the speedup of the new solution generator compared to the reference +// solution generator, where the speedup is defined as the ratio of the +// cumulated time of the reference solution generator to the cumulated time of +// the new solution generator. +double Speedup(SetCoverSolutionGenerator& ref_gen, + SetCoverSolutionGenerator& new_gen) { + // Avoid division by zero by considering the case where the new generator took + // less than 1 nanosecond (!) to run. + return 1.0 * RunTimeInNanoseconds(ref_gen) / + std::max(1, RunTimeInNanoseconds(new_gen)); +} + +// Same as above, but the cumulated time is the sum of the initialization and +// search times for two pairs of solution generators. +double SpeedupOnCumulatedTimes(SetCoverSolutionGenerator& ref_init, + SetCoverSolutionGenerator& ref_search, + SetCoverSolutionGenerator& new_init, + SetCoverSolutionGenerator& new_search) { + return 1.0 * + (RunTimeInNanoseconds(ref_init) + RunTimeInNanoseconds(ref_search)) / + std::max(1, RunTimeInNanoseconds(new_init) + + RunTimeInNanoseconds(new_search)); +} + +// In the case of LaTeX, the stats are printed in the format: +// & 123.456 (123) +/- 12.34 (12) & [123, 456] corresponding to +// & mean (median) +/- stddev (iqr) & [min, max] +// In the case of CSV, the stats are printed as a VerboseString. +std::string FormatStats(const SetCoverModel::Stats& stats) { + return absl::GetFlag(FLAGS_latex) + ? absl::StrFormat( + " & $%#.2f (%.0f) \\pm %#.2f (%.0f)$ & $[%.0f, %.0f]$", + stats.mean, stats.median, stats.stddev, stats.iqr, stats.min, + stats.max) + : stats.ToVerboseString(", "); +} + +// Returns the string str in LaTeX bold if condition is true and --latex is +// true. +std::string BoldIf(bool condition, absl::string_view str) { + return condition && absl::GetFlag(FLAGS_latex) + ? absl::StrCat("\\textbf{", str, "}") + : std::string(str); +} + +// Formats time in microseconds for LaTeX. For CSV, it is formatted as +// seconds by adding the suffix "e-6, s". +std::string FormatTime(int64_t time_us) { + return absl::GetFlag(FLAGS_latex) ? absl::StrFormat("%d", time_us) + : absl::StrFormat("%de-6, s", time_us); +} + +// Formats the cost and time, with cost in bold if the condition is true. +std::string FormatCostAndTimeIf(bool condition, double cost, int64_t time_us) { + const std::string cost_display = + BoldIf(condition, absl::StrFormat("%.9g", cost)); + return absl::StrCat(cost_display, Separator(), FormatTime(time_us)); +} + +// Formats the speedup with one decimal place. +// TODO(user): Bolden the speedup if it is better than 1.0. +std::string FormattedSpeedup(SetCoverSolutionGenerator& ref_gen, + SetCoverSolutionGenerator& new_gen) { + return absl::StrFormat("%.1f", Speedup(ref_gen, new_gen)); +} + +// Reports the second of the comparison of two solution generators, with only +// the cost and time of the new solution generator. +std::string ReportSecondPart(SetCoverSolutionGenerator& ref_gen, + SetCoverSolutionGenerator& new_gen) { + const double ref_cost = ref_gen.cost(); + const double new_cost = new_gen.cost(); + return absl::StrCat(Separator(), + FormatCostAndTimeIf(new_cost <= ref_cost, new_cost, + RunTimeInMicroseconds(new_gen)), + Separator(), FormattedSpeedup(ref_gen, new_gen)); +} + +// Reports the cost and time of both solution generators. +std::string ReportBothParts(SetCoverSolutionGenerator& ref_gen, + SetCoverSolutionGenerator& new_gen) { + const double ref_cost = ref_gen.cost(); + const double new_cost = new_gen.cost(); + return absl::StrCat(Separator(), + FormatCostAndTimeIf(ref_cost <= new_cost, ref_cost, + RunTimeInMicroseconds(ref_gen)), + ReportSecondPart(ref_gen, new_gen)); +} + +// Formats the model and stats in LaTeX or CSV format. +std::string FormatModelAndStats(const SetCoverModel& model) { + if (absl::GetFlag(FLAGS_latex)) { + return absl::StrCat(model.name(), Separator(), model.ToString(Separator()), + FormatStats(model.ComputeColumnStats()), + FormatStats(model.ComputeRowStats())); + } else { // CSV + const std::string header = + absl::StrCat(Separator(), model.name(), Separator()); + return absl::StrCat( + header, model.ToString(Separator()), Eol(), header, "column stats", + Separator(), FormatStats(model.ComputeColumnStats()), Eol(), header, + "row stats", Separator(), FormatStats(model.ComputeRowStats())); + } +} + +// Sets the name of the model to the filename, with a * suffix if the model is +// unicost. Removes the .txt suffix from the filename if any. +void SetModelName(absl::string_view filename, SetCoverModel* model) { + std::string problem_name = std::string(filename); + constexpr absl::string_view kTxtSuffix = ".txt"; + // Remove the .txt suffix from the problem name if any. + if (absl::EndsWith(problem_name, kTxtSuffix)) { + problem_name.resize(problem_name.size() - kTxtSuffix.size()); + } + if (absl::GetFlag(FLAGS_unicost) || model->is_unicost()) { + absl::StrAppend(&problem_name, "*"); + } + model->SetName(problem_name); +} + +// Accumulates data to compute the geometric mean of a sequence of values. +class GeometricMean { + public: + GeometricMean() : sum_(0.0), count_(0) {} + void Add(double value) { + sum_ += std::log(value); + ++count_; + } + double Get() const { return std::exp(sum_ / count_); } + + private: + double sum_; + int count_; +}; + +std::vector BuildVector(const char* const files[], int size) { + return std::vector(files, files + size); +} +} // namespace + +// List all the files from the literature. +static const char* const kRailFiles[] = { + "rail507.txt", "rail516.txt", "rail582.txt", "rail2536.txt", + "rail2586.txt", "rail4284.txt", "rail4872.txt", +}; + +static const char* const kScp4To6Files[] = { + "scp41.txt", "scp42.txt", "scp43.txt", "scp44.txt", "scp45.txt", + "scp46.txt", "scp47.txt", "scp48.txt", "scp49.txt", "scp410.txt", + "scp51.txt", "scp52.txt", "scp53.txt", "scp54.txt", "scp55.txt", + "scp56.txt", "scp57.txt", "scp58.txt", "scp59.txt", "scp510.txt", + "scp61.txt", "scp62.txt", "scp63.txt", "scp64.txt", "scp65.txt", +}; + +static const char* const kScpAToEFiles[] = { + "scpa1.txt", "scpa2.txt", "scpa3.txt", "scpa4.txt", "scpa5.txt", + "scpb1.txt", "scpb2.txt", "scpb3.txt", "scpb4.txt", "scpb5.txt", + "scpc1.txt", "scpc2.txt", "scpc3.txt", "scpc4.txt", "scpc5.txt", + "scpd1.txt", "scpd2.txt", "scpd3.txt", "scpd4.txt", "scpd5.txt", + "scpe1.txt", "scpe2.txt", "scpe3.txt", "scpe4.txt", "scpe5.txt", +}; + +static const char* const kScpNrFiles[] = { + "scpnre1.txt", "scpnre2.txt", "scpnre3.txt", "scpnre4.txt", "scpnre5.txt", + "scpnrf1.txt", "scpnrf2.txt", "scpnrf3.txt", "scpnrf4.txt", "scpnrf5.txt", + "scpnrg1.txt", "scpnrg2.txt", "scpnrg3.txt", "scpnrg4.txt", "scpnrg5.txt", + "scpnrh1.txt", "scpnrh2.txt", "scpnrh3.txt", "scpnrh4.txt", "scpnrh5.txt", +}; + +static const char* const kScpClrFiles[] = { + "scpclr10.txt", + "scpclr11.txt", + "scpclr12.txt", + "scpclr13.txt", +}; + +static const char* const kScpCycFiles[] = { + "scpcyc06.txt", "scpcyc07.txt", "scpcyc08.txt", + "scpcyc09.txt", "scpcyc10.txt", "scpcyc11.txt", +}; + +static const char* const kWedelinFiles[] = { + "a320_coc.txt", "a320.txt", "alitalia.txt", + "b727.txt", "sasd9imp2.txt", "sasjump.txt", +}; + +static const char* const kBalasFiles[] = { + "aa03.txt", "aa04.txt", "aa05.txt", "aa06.txt", "aa11.txt", "aa12.txt", + "aa13.txt", "aa14.txt", "aa15.txt", "aa16.txt", "aa17.txt", "aa18.txt", + "aa19.txt", "aa20.txt", "bus1.txt", "bus2.txt", +}; + +static const char* const kFimiFiles[] = { + "accidents.dat", "chess.dat", "connect.dat", "kosarak.dat", + "mushroom.dat", // "pumsb.dat", "pumsb_star.dat", + "retail.dat", "webdocs.dat", +}; + +using BenchmarksTableRow = + std::tuple, FileFormat>; + +std::vector BenchmarksTable() { +// This creates a vector of tuples, where each tuple contains the directory +// name, the vector of files and the file format. +// It is assumed that the scp* files are in BENCHMARKS_DIR/scp-orlib, the rail +// files are in BENCHMARKS_DIR/scp-rail, etc., with BENCHMARKS_DIR being the +// directory specified by the --benchmarks_dir flag. +// Use a macro to be able to compute the size of the array at compile time. +#define BUILD_VECTOR(files) BuildVector(files, ABSL_ARRAYSIZE(files)) + return std::vector{ + {"scp-orlib", BUILD_VECTOR(kScp4To6Files), FileFormat::ORLIB}, + {"scp-orlib", BUILD_VECTOR(kScpAToEFiles), FileFormat::ORLIB}, + {"scp-orlib", BUILD_VECTOR(kScpNrFiles), FileFormat::ORLIB}, + {"scp-orlib", BUILD_VECTOR(kScpClrFiles), FileFormat::ORLIB}, + {"scp-orlib", BUILD_VECTOR(kScpCycFiles), FileFormat::ORLIB}, + {"scp-rail", BUILD_VECTOR(kRailFiles), FileFormat::RAIL}, + {"scp-wedelin", BUILD_VECTOR(kWedelinFiles), FileFormat::ORLIB}, + {"scp-balas", BUILD_VECTOR(kBalasFiles), FileFormat::ORLIB}, + {"scp-fimi", BUILD_VECTOR(kFimiFiles), FileFormat::FIMI}, + }; +#undef BUILD_VECTOR +} + +void Benchmarks() { + QCHECK(!absl::GetFlag(FLAGS_benchmarks_dir).empty()) + << "Benchmarks directory must be specified."; + const std::vector kBenchmarks = BenchmarksTable(); + const std::string kSep = Separator(); + const std::string kEol = Eol(); + if (absl::GetFlag(FLAGS_stats)) { + std::cout << absl::StrJoin(std::make_tuple("Problem", "|S|", "|U|", "nnz", + "Fill", "Col size", "Row size"), + kSep) + << kEol; + } + for (const auto& [dir, files, input_format] : kBenchmarks) { + GeometricMean first_cost_ratio_geomean; + GeometricMean first_speedup_geomean; + GeometricMean improvement_cost_ratio_geomean; + GeometricMean improvement_speedup_geomean; + GeometricMean combined_cost_ratio_geomean; + GeometricMean combined_speedup_geomean; + GeometricMean tlns_cost_ratio_geomean; + for (const std::string& filename : files) { + const std::string filespec = absl::StrCat( + absl::GetFlag(FLAGS_benchmarks_dir), "/", dir, "/", filename); + + LOG(INFO) << "Reading " << filespec; + SetCoverModel model = ReadModel(filespec, input_format); + if (absl::GetFlag(FLAGS_unicost)) { + for (const SubsetIndex subset : model.SubsetRange()) { + model.SetSubsetCost(subset, 1.0); + } + } + SetModelName(filename, &model); + if (absl::GetFlag(FLAGS_stats)) { + std::cout << FormatModelAndStats(model) << kEol; + } + if (!absl::GetFlag(FLAGS_solve)) continue; + LOG(INFO) << "Solving " << model.name(); + std::string output = + absl::StrJoin(std::make_tuple(model.name(), model.num_subsets(), + model.num_elements()), + kSep); + model.CreateSparseRowView(); + + SetCoverInvariant classic_inv(&model); + GreedySolutionGenerator classic(&classic_inv); + SteepestSearch steepest(&classic_inv); + const bool run_classic_solvers = + model.num_elements() <= absl::GetFlag(FLAGS_max_elements_for_classic); + if (run_classic_solvers) { + CHECK(classic.NextSolution()); + CHECK(steepest.NextSolution()); + DCHECK(classic_inv.CheckConsistency(CL::kCostAndCoverage)); + } + + SetCoverInvariant inv(&model); + LazyElementDegreeSolutionGenerator lazy_element_degree(&inv); + CHECK(lazy_element_degree.NextSolution()); + DCHECK(inv.CheckConsistency(CL::kCostAndCoverage)); + + absl::StrAppend(&output, ReportBothParts(classic, lazy_element_degree)); + LazySteepestSearch lazy_steepest(&inv); + lazy_steepest.NextSolution(); + absl::StrAppend(&output, ReportBothParts(steepest, lazy_steepest)); + + if (run_classic_solvers) { + first_cost_ratio_geomean.Add(CostRatio(classic, lazy_element_degree)); + first_speedup_geomean.Add(Speedup(classic, lazy_element_degree)); + improvement_cost_ratio_geomean.Add(CostRatio(steepest, lazy_steepest)); + improvement_speedup_geomean.Add(Speedup(steepest, lazy_steepest)); + combined_cost_ratio_geomean.Add(CostRatio(classic, lazy_steepest)); + combined_speedup_geomean.Add(SpeedupOnCumulatedTimes( + classic, steepest, lazy_element_degree, lazy_steepest)); + } + + Cost best_cost = inv.cost(); + SubsetBoolVector best_assignment = inv.is_selected(); + if (absl::GetFlag(FLAGS_tlns)) { + WallTimer timer; + LazySteepestSearch steepest(&inv); + timer.Start(); + for (int i = 0; i < 500; ++i) { + std::vector range = + ClearRandomSubsets(0.1 * inv.trace().size(), &inv); + CHECK(lazy_element_degree.NextSolution()); + CHECK(steepest.NextSolution()); + if (inv.cost() < best_cost) { + best_cost = inv.cost(); + best_assignment = inv.is_selected(); + } + } + timer.Stop(); + absl::StrAppend( + &output, kSep, + FormatCostAndTimeIf( + best_cost <= lazy_element_degree.cost() && + best_cost <= classic.cost(), + best_cost, absl::ToInt64Microseconds(timer.GetDuration()))); + tlns_cost_ratio_geomean.Add(best_cost / classic.cost()); + } + std::cout << output << kEol; + } + + if (absl::GetFlag(FLAGS_summarize)) { + std::cout << "First speedup geometric mean: " + << first_speedup_geomean.Get() << kEol + << "First cost ratio geometric mean: " + << first_cost_ratio_geomean.Get() << kEol + << "Improvement speedup geometric mean: " + << improvement_speedup_geomean.Get() << kEol + << "Improvement cost ratio geometric mean: " + << improvement_cost_ratio_geomean.Get() << kEol + << "Combined speedup geometric mean: " + << combined_speedup_geomean.Get() << kEol; + if (absl::GetFlag(FLAGS_tlns)) { + std::cout << "TLNS cost ratio geometric mean: " + << tlns_cost_ratio_geomean.Get() << kEol; + } + } + } +} + 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)); - if (input.empty()) { - LOG(FATAL) << "No input file specified."; - } - if (!input.empty() && input_format == FileFormat::EMPTY) { - LOG(FATAL) << "Input format cannot be empty."; - } - if (!output.empty() && output_format == FileFormat::EMPTY) { - LOG(FATAL) << "Output format cannot be empty."; - } + QCHECK(!input.empty()) << "No input file specified."; + QCHECK(input.empty() || input_format != FileFormat::EMPTY) + << "Input format cannot be empty."; + QCHECK(output.empty() || output_format != FileFormat::EMPTY) + << "Output format cannot be empty."; SetCoverModel model = ReadModel(input, input_format); if (absl::GetFlag(FLAGS_generate)) { model.CreateSparseRowView(); @@ -253,25 +656,29 @@ void Run() { absl::GetFlag(FLAGS_column_scale), absl::GetFlag(FLAGS_cost_scale)); } if (!output.empty()) { - if (output_format == FileFormat::ORLIB_SCP) { + if (output_format == FileFormat::ORLIB) { model.CreateSparseRowView(); } WriteModel(model, output, output_format); } auto problem = output.empty() ? input : output; if (absl::GetFlag(FLAGS_stats)) { - LogStats(problem, &model); + LogStats(model); } if (absl::GetFlag(FLAGS_solve)) { LOG(INFO) << "Solving " << problem; model.CreateSparseRowView(); - SetCoverInvariant inv = RunLazyElementDegree(problem, &model); + SetCoverInvariant inv = RunLazyElementDegree(&model); } } } // namespace operations_research int main(int argc, char** argv) { InitGoogle(argv[0], &argc, &argv, true); - operations_research::Run(); + if (absl::GetFlag(FLAGS_benchmarks)) { + operations_research::Benchmarks(); + } else { + operations_research::Run(); + } return 0; } diff --git a/ortools/set_cover/set_cover_test.cc b/ortools/set_cover/set_cover_test.cc index 6a7e7ca1ec..7a8dfdd7dd 100644 --- a/ortools/set_cover/set_cover_test.cc +++ b/ortools/set_cover/set_cover_test.cc @@ -474,7 +474,8 @@ TEST(SetCoverTest, KnightsCoverRandomClearMip) { for (int i = 0; i < 1'000; ++i) { auto focus = ClearRandomSubsets(0.1 * inv.trace().size(), &inv); SetCoverMip mip(&inv); - mip.NextSolution(focus, true, 1); + mip.UseIntegers(true).SetTimeLimitInSeconds(1.0); + mip.NextSolution(focus); EXPECT_TRUE(inv.CheckConsistency(CL::kCostAndCoverage)); if (inv.cost() < best_cost) { best_cost = inv.cost(); @@ -502,7 +503,7 @@ TEST(SetCoverTest, KnightsCoverMip) { SetCoverModel model = knights.model(); SetCoverInvariant inv(&model); SetCoverMip mip(&inv); - mip.NextSolution(true, .5); + mip.UseIntegers(true).SetTimeLimitInSeconds(0.5).NextSolution(); LOG(INFO) << "Mip cost: " << inv.cost(); knights.DisplaySolution(inv.is_selected()); if (BoardSize == 50) {