From 8b4a24a0a3c48c9cc8484fba216c31edb268ee32 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Fri, 7 Feb 2025 22:24:22 +0100 Subject: [PATCH] [CP-SAT] polish new cumulative cuts code; more experimental code on routing cuts; cleaning, fixes --- ortools/sat/BUILD.bazel | 3 + ortools/sat/cp_model_search.cc | 49 +- ortools/sat/cp_model_solver.cc | 85 +- ortools/sat/cuts.h | 4 +- ortools/sat/integer_expr.cc | 5 + ortools/sat/linear_programming_constraint.cc | 2 +- ortools/sat/parameters_validation.cc | 2 + ortools/sat/precedences.cc | 53 ++ ortools/sat/precedences.h | 11 + ortools/sat/routing_cuts.cc | 849 ++++++++++++------- ortools/sat/routing_cuts.h | 32 +- ortools/sat/routing_cuts_test.cc | 24 +- ortools/sat/sat_parameters.proto | 22 +- ortools/sat/scheduling_cuts.cc | 210 +++-- ortools/sat/scheduling_helpers.cc | 24 +- ortools/sat/util.cc | 105 +++ ortools/sat/util.h | 48 ++ ortools/sat/util_test.cc | 77 ++ 18 files changed, 1131 insertions(+), 474 deletions(-) diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index a61e041b5e..e436d225da 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -430,6 +430,7 @@ cc_library( "//ortools/base", "//ortools/base:stl_util", "//ortools/util:strong_integers", + "@com_google_absl//absl/algorithm:container", "@com_google_absl//absl/container:flat_hash_map", "@com_google_absl//absl/container:flat_hash_set", "@com_google_absl//absl/log:check", @@ -1849,6 +1850,7 @@ cc_library( ":synchronization", ":util", "//ortools/base", + "//ortools/base:mathutil", "//ortools/base:stl_util", "//ortools/base:strong_vector", "//ortools/graph", @@ -2609,6 +2611,7 @@ cc_library( "@com_google_absl//absl/container:btree", "@com_google_absl//absl/container:flat_hash_map", "@com_google_absl//absl/container:flat_hash_set", + "@com_google_absl//absl/functional:any_invocable", "@com_google_absl//absl/log", "@com_google_absl//absl/log:check", "@com_google_absl//absl/meta:type_traits", diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index 5e16b41795..f8b210d70f 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -22,6 +22,7 @@ #include #include +#include "absl/algorithm/container.h" #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" #include "absl/log/check.h" @@ -741,20 +742,46 @@ absl::flat_hash_map GetNamedParameters( // Base parameters for LNS worker. { - SatParameters new_params = base_params; - new_params.set_stop_after_first_solution(false); - new_params.set_cp_model_presolve(true); + SatParameters lns_params = base_params; + lns_params.set_stop_after_first_solution(false); + lns_params.set_cp_model_presolve(true); // We disable costly presolve/inprocessing. - new_params.set_use_sat_inprocessing(false); - new_params.set_cp_model_probing_level(0); - new_params.set_symmetry_level(0); - new_params.set_find_big_linear_overlap(false); + lns_params.set_use_sat_inprocessing(false); + lns_params.set_cp_model_probing_level(0); + lns_params.set_symmetry_level(0); + lns_params.set_find_big_linear_overlap(false); - new_params.set_log_search_progress(false); - new_params.set_debug_crash_on_bad_hint(false); // Can happen in lns. - new_params.set_solution_pool_size(1); // Keep the best solution found. - strategies["lns"] = new_params; + lns_params.set_log_search_progress(false); + lns_params.set_debug_crash_on_bad_hint(false); // Can happen in lns. + lns_params.set_solution_pool_size(1); // Keep the best solution found. + strategies["lns"] = lns_params; + + // Note that we only do this for the derived parameters. The strategy "lns" + // will be handled along with the other ones. + auto it = absl::c_find_if( + base_params.subsolver_params(), + [](const SatParameters& params) { return params.name() == "lns"; }); + if (it != base_params.subsolver_params().end()) { + lns_params.MergeFrom(*it); + } + + SatParameters lns_params_base = lns_params; + lns_params_base.set_linearization_level(0); + lns_params_base.set_search_branching(SatParameters::AUTOMATIC_SEARCH); + strategies["lns_base"] = lns_params_base; + + SatParameters lns_params_stalling = lns_params; + lns_params_stalling.set_search_branching(SatParameters::PORTFOLIO_SEARCH); + lns_params_stalling.set_search_random_variable_pool_size(5); + strategies["lns_stalling"] = lns_params_stalling; + + // For routing, the LP relaxation seems pretty important, so we prefer an + // high linearization level to solve LNS subproblems. + SatParameters lns_params_routing = lns_params; + lns_params_routing.set_linearization_level(2); + lns_params_routing.set_search_branching(SatParameters::AUTOMATIC_SEARCH); + strategies["lns_routing"] = lns_params_routing; } // Add user defined ones. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 5c6a0bb572..67d3ac6723 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -46,6 +47,7 @@ #include "absl/strings/str_format.h" #include "absl/strings/str_join.h" #include "absl/strings/string_view.h" +#include "absl/strings/strip.h" #include "absl/synchronization/mutex.h" #include "absl/types/span.h" #include "google/protobuf/arena.h" @@ -1293,14 +1295,14 @@ class FeasibilityPumpSolver : public SubSolver { class LnsSolver : public SubSolver { public: LnsSolver(std::unique_ptr generator, - const SatParameters& lns_parameters, - NeighborhoodGeneratorHelper* helper, SharedClasses* shared, - int preferred_linearization_level = 0) + const SatParameters& lns_parameters_base, + const SatParameters& lns_parameters_stalling, + NeighborhoodGeneratorHelper* helper, SharedClasses* shared) : SubSolver(generator->name(), INCOMPLETE), - preferred_linearization_level_(preferred_linearization_level), generator_(std::move(generator)), helper_(helper), - lns_parameters_(lns_parameters), + lns_parameters_base_(lns_parameters_base), + lns_parameters_stalling_(lns_parameters_stalling), shared_(shared) {} ~LnsSolver() override { @@ -1328,7 +1330,7 @@ class LnsSolver : public SubSolver { // change the LNS behavior. const int32_t low = static_cast(task_id); const int32_t high = static_cast(task_id >> 32); - std::seed_seq seed{low, high, lns_parameters_.random_seed()}; + std::seed_seq seed{low, high, lns_parameters_base_.random_seed()}; random_engine_t random(seed); NeighborhoodGenerator::SolveData data; @@ -1372,26 +1374,22 @@ class LnsSolver : public SubSolver { if (!neighborhood.is_generated) return; - SatParameters local_params(lns_parameters_); - local_params.set_max_deterministic_time(data.deterministic_limit); + SatParameters local_params; - // TODO(user): Tune these. // TODO(user): This could be a good candidate for bandits. const int64_t stall = generator_->num_consecutive_non_improving_calls(); - std::string search_info; const int search_index = stall < 10 ? 0 : task_id % 2; switch (search_index) { case 0: - local_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH); - local_params.set_linearization_level(preferred_linearization_level_); - search_info = absl::StrCat("auto_l", preferred_linearization_level_); + local_params = lns_parameters_base_; break; default: - local_params.set_search_branching(SatParameters::PORTFOLIO_SEARCH); - local_params.set_search_random_variable_pool_size(5); - search_info = "folio_rnd"; + local_params = lns_parameters_stalling_; break; } + const std::string_view search_info = + absl::StripPrefix(std::string_view(local_params.name()), "lns_"); + local_params.set_max_deterministic_time(data.deterministic_limit); std::string source_info = neighborhood.source_info.empty() ? name() : neighborhood.source_info; @@ -1486,7 +1484,7 @@ class LnsSolver : public SubSolver { } } bool hint_feasible_before_presolve = false; - if (lns_parameters_.debug_crash_if_presolve_breaks_hint()) { + if (lns_parameters_base_.debug_crash_if_presolve_breaks_hint()) { hint_feasible_before_presolve = SolutionHintIsCompleteAndFeasible(lns_fragment, /*logger=*/nullptr); } @@ -1525,7 +1523,7 @@ class LnsSolver : public SubSolver { context.reset(nullptr); neighborhood.delta.Clear(); - if (lns_parameters_.debug_crash_if_presolve_breaks_hint() && + if (lns_parameters_base_.debug_crash_if_presolve_breaks_hint() && hint_feasible_before_presolve && !SolutionHintIsCompleteAndFeasible(lns_fragment, /*logger=*/nullptr)) { @@ -1720,10 +1718,10 @@ class LnsSolver : public SubSolver { } private: - int preferred_linearization_level_ = 0; std::unique_ptr generator_; NeighborhoodGeneratorHelper* helper_; - const SatParameters lns_parameters_; + const SatParameters lns_parameters_base_; + const SatParameters lns_parameters_stalling_; SharedClasses* shared_; // This is a optimization to allocate the arena for the LNS fragment already // at roughly the right size. We will update it with the last size of the @@ -1780,7 +1778,9 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { })); const auto name_to_params = GetNamedParameters(params); - const SatParameters& lns_params = name_to_params.at("lns"); + const SatParameters& lns_params_base = name_to_params.at("lns_base"); + const SatParameters& lns_params_stalling = name_to_params.at("lns_stalling"); + const SatParameters& lns_params_routing = name_to_params.at("lns_routing"); // Add the NeighborhoodGeneratorHelper as a special subsolver so that its // Synchronize() is called before any LNS neighborhood solvers. @@ -1854,7 +1854,7 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { std::make_unique( helper, shared->response, shared->lp_solutions.get(), shared->incomplete_solutions.get(), name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } // Add incomplete subsolvers that require an objective. @@ -1870,37 +1870,37 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("rnd_cst_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("graph_var_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("graph_arc_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("graph_cst_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("graph_dec_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (params.use_lb_relax_lns() && params.num_workers() >= params.lb_relax_num_workers_threshold() && @@ -1908,7 +1908,7 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName(), shared->time_limit, shared), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } const bool has_no_overlap_or_cumulative = @@ -1921,13 +1921,13 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("scheduling_time_window_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } const std::vector> intervals_in_constraints = helper->GetUniqueIntervalSets(); @@ -1936,7 +1936,7 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, intervals_in_constraints, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } } @@ -1948,31 +1948,31 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("packing_square_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("packing_swap_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("packing_precedences_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } if (name_filter.Keep("packing_slice_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } } @@ -1982,13 +1982,10 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared)); + lns_params_base, lns_params_stalling, helper, shared)); } } - // For routing, the LP relaxation seems pretty important, so we prefer an - // high linearization level to solve LNS subproblems. - const int routing_lin_level = 2; const int num_circuit = static_cast( helper->TypeToConstraints(ConstraintProto::kCircuit).size()); const int num_routes = static_cast( @@ -1998,13 +1995,13 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared, routing_lin_level)); + lns_params_routing, lns_params_stalling, helper, shared)); } if (name_filter.Keep("routing_path_lns")) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared, routing_lin_level)); + lns_params_routing, lns_params_stalling, helper, shared)); } } if (num_routes > 0 || num_circuit > 1) { @@ -2012,7 +2009,7 @@ void SolveCpModelParallel(SharedClasses* shared, Model* global_model) { reentrant_interleaved_subsolvers.push_back(std::make_unique( std::make_unique( helper, name_filter.LastName()), - lns_params, helper, shared, routing_lin_level)); + lns_params_routing, lns_params_stalling, helper, shared)); } } } diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index ab5563494b..37060cb6d7 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -29,11 +29,11 @@ #include "absl/container/btree_set.h" #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" +#include "absl/functional/any_invocable.h" #include "absl/numeric/int128.h" #include "absl/strings/str_cat.h" #include "absl/types/span.h" #include "ortools/base/strong_vector.h" -#include "ortools/lp_data/lp_types.h" #include "ortools/sat/implied_bounds.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_base.h" @@ -56,7 +56,7 @@ namespace sat { struct CutGenerator { bool only_run_at_level_zero = false; std::vector vars; - std::function generate_cuts; + absl::AnyInvocable generate_cuts; }; // To simplify cut generation code, we use a more complex data structure than diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index 011ec4468e..c217d0fb8c 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -828,6 +828,7 @@ bool LinMinPropagator::Propagate() { void LinMinPropagator::RegisterWith(GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); + bool has_var_also_in_exprs = false; for (const LinearExpression& expr : exprs_) { for (int i = 0; i < expr.vars.size(); ++i) { const IntegerVariable& var = expr.vars[i]; @@ -837,10 +838,14 @@ void LinMinPropagator::RegisterWith(GenericLiteralWatcher* watcher) { } else { watcher->WatchUpperBound(var, id); } + has_var_also_in_exprs |= (var == min_var_); } } watcher->WatchUpperBound(min_var_, id); watcher->RegisterReversibleInt(id, &rev_unique_candidate_); + if (has_var_also_in_exprs) { + watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); + } } ProductPropagator::ProductPropagator(AffineExpression a, AffineExpression b, diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index af082d5100..33521bebc6 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -2192,7 +2192,7 @@ bool LinearProgrammingConstraint::Propagate() { // Try to add cuts. if (level == 0 || !parameters_.only_add_cuts_at_level_zero()) { - for (const CutGenerator& generator : cut_generators_) { + for (CutGenerator& generator : cut_generators_) { if (level > 0 && generator.only_run_at_level_zero) continue; if (!generator.generate_cuts(&constraint_manager_)) { return false; diff --git a/ortools/sat/parameters_validation.cc b/ortools/sat/parameters_validation.cc index c9ed2efb56..4d64b1aaad 100644 --- a/ortools/sat/parameters_validation.cc +++ b/ortools/sat/parameters_validation.cc @@ -92,6 +92,8 @@ std::string ValidateParameters(const SatParameters& params) { TEST_IS_FINITE(strategy_change_increase_ratio); TEST_IS_FINITE(symmetry_detection_deterministic_time_limit); TEST_IS_FINITE(variable_activity_decay); + TEST_IN_RANGE(routing_cut_subset_size_for_tight_binary_relation_bound, 0, 32); + TEST_IS_FINITE(routing_cut_dp_effort); TEST_IS_FINITE(lns_initial_difficulty); TEST_IS_FINITE(lns_initial_deterministic_limit); diff --git a/ortools/sat/precedences.cc b/ortools/sat/precedences.cc index 5151fba7d0..7d30519f07 100644 --- a/ortools/sat/precedences.cc +++ b/ortools/sat/precedences.cc @@ -31,6 +31,7 @@ #include "absl/meta/type_traits.h" #include "absl/types/span.h" #include "ortools/base/logging.h" +#include "ortools/base/mathutil.h" #include "ortools/base/stl_util.h" #include "ortools/base/strong_vector.h" #include "ortools/graph/graph.h" @@ -1387,5 +1388,57 @@ int GreaterThanAtLeastOneOfDetector::AddGreaterThanAtLeastOneOfConstraints( return num_added_constraints; } +bool BinaryRelationRepository::PropagateLocalBounds( + const IntegerTrail& integer_trail, Literal lit, + const absl::flat_hash_map& input, + absl::flat_hash_map* output) const { + output->clear(); + if (lit.Index() >= lit_to_relations_.size()) return true; + + auto get_lower_bound = [&](IntegerVariable var) { + const auto it = input.find(var); + if (it != input.end()) return it->second; + return integer_trail.LevelZeroLowerBound(var); + }; + auto get_upper_bound = [&](IntegerVariable var) { + return -get_lower_bound(NegationOf(var)); + }; + auto update_lower_bound_by_var = [&](IntegerVariable var, IntegerValue lb) { + if (lb <= integer_trail.LevelZeroLowerBound(var)) return; + const auto [it, inserted] = output->insert({var, lb}); + if (!inserted) { + it->second = std::max(it->second, lb); + } + }; + auto update_upper_bound_by_var = [&](IntegerVariable var, IntegerValue ub) { + update_lower_bound_by_var(NegationOf(var), -ub); + }; + auto update_var_bounds = [&](const LinearTerm& a, const LinearTerm& b, + IntegerValue lhs, IntegerValue rhs) { + if (a.coeff == 0) return; + + // lb(b.y) <= b.y <= ub(b.y) and lhs <= a.x + b.y <= rhs imply + // ceil((lhs - ub(b.y)) / a) <= x <= floor((rhs - lb(b.y)) / a) + lhs = lhs - b.coeff * get_upper_bound(b.var); + rhs = rhs - b.coeff * get_lower_bound(b.var); + update_lower_bound_by_var(a.var, MathUtil::CeilOfRatio(lhs, a.coeff)); + update_upper_bound_by_var(a.var, MathUtil::FloorOfRatio(rhs, a.coeff)); + }; + for (const int relation_index : lit_to_relations_[lit]) { + auto r = relations_[relation_index]; + r.a.MakeCoeffPositive(); + r.b.MakeCoeffPositive(); + update_var_bounds(r.a, r.b, r.lhs, r.rhs); + update_var_bounds(r.b, r.a, r.lhs, r.rhs); + } + + // Check feasibility. + // TODO(user): we might do that earlier? + for (const auto [var, lb] : *output) { + if (lb > integer_trail.LevelZeroUpperBound(var)) return false; + } + return true; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/precedences.h b/ortools/sat/precedences.h index 2640ab6803..f79cc665f2 100644 --- a/ortools/sat/precedences.h +++ b/ortools/sat/precedences.h @@ -509,6 +509,17 @@ class BinaryRelationRepository { // relations have been added. void Build(); + // Assuming level-zero bounds + any (var >= value) in the input map, + // fills "output" with a "propagated" set of bounds assuming lit is true (by + // using the relations enforced by lit). Note that we will only fill bounds > + // level-zero ones in output. + // + // Returns false if the new bounds are infeasible at level zero. + bool PropagateLocalBounds( + const IntegerTrail& integer_trail, Literal lit, + const absl::flat_hash_map& input, + absl::flat_hash_map* output) const; + private: bool is_built_ = false; std::vector relations_; diff --git a/ortools/sat/routing_cuts.cc b/ortools/sat/routing_cuts.cc index 2071b8250b..ce38eda1c7 100644 --- a/ortools/sat/routing_cuts.cc +++ b/ortools/sat/routing_cuts.cc @@ -18,7 +18,10 @@ #include #include #include +#include +#include #include +#include #include #include @@ -33,6 +36,7 @@ #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" #include "ortools/base/strong_vector.h" +#include "ortools/graph/graph.h" #include "ortools/graph/max_flow.h" #include "ortools/sat/cuts.h" #include "ortools/sat/integer.h" @@ -42,6 +46,7 @@ #include "ortools/sat/model.h" #include "ortools/sat/precedences.h" #include "ortools/sat/sat_base.h" +#include "ortools/sat/util.h" #include "ortools/util/strong_integers.h" namespace operations_research { @@ -49,129 +54,38 @@ namespace sat { namespace { -// Helper to compute lower bounds based on binary relations enforced by arc -// literals in a RoutesConstraint. -class LowerBoundsHelper { - public: - // See MinOutgoingFlowHelper. - LowerBoundsHelper( - std::vector>& - node_var_lower_bounds, - std::vector>& - next_node_var_lower_bounds, - const IntegerTrail& integer_trail) - : integer_trail_(integer_trail), - node_var_lower_bounds_(node_var_lower_bounds), - next_node_var_lower_bounds_(next_node_var_lower_bounds) {} - - // Updates the bounds of `a.var` if the `arc_index` arc were selected, based - // on the given relation `a + b \in [lhs, rhs]` (supposedly enforced by the - // arc literal), and on the bounds of `b.var` (if the tail node appears at the - // "current" position -- see MinOutgoingFlowHelper). - void UpdateNextVarBoundsByArc(int tail, int arc_index, LinearTerm a, - LinearTerm b, IntegerValue lhs, - IntegerValue rhs) { - if (a.coeff == 0) return; - a.MakeCoeffPositive(); - b.MakeCoeffPositive(); - // lb(b.y) <= b.y <= ub(b.y) and lhs <= a.x + b.y <= rhs imply - // ceil((lhs - ub(b.y)) / a) <= x <= floor((rhs - lb(b.y)) / a) - lhs = lhs + b.coeff * GetCurrentVarLowerBound(tail, NegationOf(b.var)); - rhs = rhs - b.coeff * GetCurrentVarLowerBound(tail, b.var); - UpdateLowerBoundByArc(a.var, arc_index, - MathUtil::CeilOfRatio(lhs, a.coeff)); - UpdateUpperBoundByArc(a.var, arc_index, - MathUtil::FloorOfRatio(rhs, a.coeff)); - } - - // Updates the bounds of the variables if `head` appears at the "next" - // position in a feasible path, and if it can be reached by - // `num_reachable_incoming_arcs` arcs. `head` must be the head of the arcs - // used in all the previous calls to UpdateNextVarBoundsByArc(). - // - // Returns whether the new bounds of the variables are feasible. If not, - // `head` cannot appear at the "next" position in a feasible path. - bool UpdateNextVarBounds(int head, int num_reachable_incoming_arcs) { - if (num_reachable_incoming_arcs == 0) return false; - for (const auto& [var, lower_bound_by_arc_index] : - lower_bound_by_var_and_arc_index_) { - // If each arc which can reach `head` enforces some lower bound for `var`, - // then the lower bound of `var` can be increased to the minimum of these - // arc-specific lower bounds (since at least one of these arcs must be - // selected to reach `head`). - if (lower_bound_by_arc_index.size() == num_reachable_incoming_arcs) { - IntegerValue lb = lower_bound_by_arc_index.begin()->second; - for (const auto& [arc_index, lower_bound] : lower_bound_by_arc_index) { - lb = std::min(lb, lower_bound); - } - if (!UpdateNextVarLowerBound(head, var, lb)) return false; - } - } - return true; - } - - private: - // Returns the lower bound of `var` if `node` appears at the current position - // in a feasible path. - IntegerValue GetCurrentVarLowerBound(int node, IntegerVariable var) const { - auto it = node_var_lower_bounds_[node].find(var); - if (it != node_var_lower_bounds_[node].end()) return it->second; - return integer_trail_.LevelZeroLowerBound(var); - } - - // Returns the lower bound of `var` if `node` appears at the next position in - // a feasible path. - IntegerValue GetNextVarLowerBound(int node, IntegerVariable var) const { - auto it = next_node_var_lower_bounds_[node].find(var); - if (it != next_node_var_lower_bounds_[node].end()) return it->second; - return integer_trail_.LevelZeroLowerBound(var); - } - - // Sets the lower bound of `var` if `node` appears at the next position in a - // feasible path to `lb`, if `lb` is greater than the current value. Returns - // whether the new bounds of `var` are feasible. - bool UpdateNextVarLowerBound(int node, IntegerVariable var, IntegerValue lb) { - if (lb <= integer_trail_.LevelZeroLowerBound(var)) return true; - auto it = next_node_var_lower_bounds_[node].find(var); - if (it != next_node_var_lower_bounds_[node].end()) { - lb = std::max(lb, it->second); - it->second = lb; +// If we have many sets of (var >= values relations) and we know at least one +// set is true, then we can derive valid global lower bounds on the variable +// that appear in all such sets. +// +// This represents "one step" of computing such bounds by adding the sets one at +// the time. When we process a set 'new_lbs': +// - For the first set, then 'new_lbs' is globally valid. +// - For a new set, we need to erase variables not appearing in new_lbs and +// take the min for the ones appearing in both. +// +// TODO(user): The operation is symmetric, so we could swap both hash_map if +// new_lbs is smaller than current_lbs as a small optimization. +void ComputeMinLowerBoundOfSharedVariables( + const absl::flat_hash_map& new_lbs, + absl::flat_hash_map* current_lbs) { + if (new_lbs.empty()) current_lbs->clear(); + for (auto it = current_lbs->begin(); it != current_lbs->end();) { + const IntegerVariable var = it->first; + auto other_it = new_lbs.find(var); + if (other_it == new_lbs.end()) { + // We have no bounds about var in the new_fact, we need to erase it. + // + // Note: it = current_lbs->erase(it) does not work for flat_hash_map, this + // is the recommended way. + current_lbs->erase(it++); } else { - next_node_var_lower_bounds_[node][var] = lb; + // Update. + it->second = std::min(it->second, other_it->second); + ++it; } - const IntegerValue ub = -GetNextVarLowerBound(node, NegationOf(var)); - return lb <= ub; } - - void UpdateLowerBoundByArc(IntegerVariable var, int arc_index, - IntegerValue lb) { - auto& lower_bound_by_arc_index = lower_bound_by_var_and_arc_index_[var]; - auto it = lower_bound_by_arc_index.find(arc_index); - if (it != lower_bound_by_arc_index.end()) { - it->second = std::max(it->second, lb); - } else { - lower_bound_by_arc_index[arc_index] = lb; - } - }; - - void UpdateUpperBoundByArc(IntegerVariable var, int arc_index, - IntegerValue ub) { - UpdateLowerBoundByArc(NegationOf(var), arc_index, -ub); - }; - - const IntegerTrail& integer_trail_; - // See MinOutgoingFlowHelper. - std::vector>& - node_var_lower_bounds_; - std::vector>& - next_node_var_lower_bounds_; - // For each variable, a map from arc index to the lower bound of that variable - // if this arc were selected, deduced from the binary relations enforced by - // the arc literal. - // TODO(user): see if we can use a more efficient data structure. - absl::flat_hash_map> - lower_bound_by_var_and_arc_index_; -}; +} } // namespace @@ -186,7 +100,9 @@ MinOutgoingFlowHelper::MinOutgoingFlowHelper( trail_(*model->GetOrCreate()), integer_trail_(*model->GetOrCreate()), in_subset_(num_nodes, false), + index_in_subset_(num_nodes, -1), incoming_arc_indices_(num_nodes), + outgoing_arc_indices_(num_nodes), reachable_(num_nodes, false), next_reachable_(num_nodes, false), node_var_lower_bounds_(num_nodes), @@ -220,28 +136,35 @@ int MinOutgoingFlowHelper::ComputeMinOutgoingFlow( // Maximum number of nodes of a feasible path inside the subset. int longest_path_length = 1; + absl::flat_hash_map tmp_lbs; while (longest_path_length < subset.size()) { bool at_least_one_next_node_reachable = false; for (const int head : subset) { - LowerBoundsHelper lower_bounds_helper( - node_var_lower_bounds_, next_node_var_lower_bounds_, integer_trail_); - int num_reachable_incoming_arcs = 0; + bool is_reachable = false; for (const int incoming_arc_index : incoming_arc_indices_[head]) { const int tail = tails_[incoming_arc_index]; const Literal lit = literals_[incoming_arc_index]; if (!reachable_[tail]) continue; - ++num_reachable_incoming_arcs; - for (const int relation_index : - binary_relation_repository_.relation_indices(lit)) { - const auto& r = binary_relation_repository_.relation(relation_index); - lower_bounds_helper.UpdateNextVarBoundsByArc(tail, incoming_arc_index, - r.a, r.b, r.lhs, r.rhs); - lower_bounds_helper.UpdateNextVarBoundsByArc(tail, incoming_arc_index, - r.b, r.a, r.lhs, r.rhs); + + // If this arc cannot be taken skip. + tmp_lbs.clear(); + if (!binary_relation_repository_.PropagateLocalBounds( + integer_trail_, lit, node_var_lower_bounds_[tail], &tmp_lbs)) { + continue; + } + + if (!is_reachable) { + // This is the first arc that reach this node. + is_reachable = true; + next_node_var_lower_bounds_[head] = tmp_lbs; + } else { + // Combine information with previous one. + ComputeMinLowerBoundOfSharedVariables( + tmp_lbs, &next_node_var_lower_bounds_[head]); } } - next_reachable_[head] = lower_bounds_helper.UpdateNextVarBounds( - head, num_reachable_incoming_arcs); + + next_reachable_[head] = is_reachable; if (next_reachable_[head]) at_least_one_next_node_reachable = true; } if (!at_least_one_next_node_reachable) { @@ -267,58 +190,256 @@ int MinOutgoingFlowHelper::ComputeMinOutgoingFlow( node_var_lower_bounds_[n].clear(); next_node_var_lower_bounds_[n].clear(); } - if (max_longest_paths * longest_path_length < subset.size()) { - // If `longest_path_length` is 1, then `reachable_` has not been modified by - // the above while loop, and thus `max_longest_paths` is the subset size. - // This branch can thus not be taken in this case. + return GetMinOutgoingFlow(subset.size(), longest_path_length, + max_longest_paths); +} + +int MinOutgoingFlowHelper::GetMinOutgoingFlow(int subset_size, + int longest_path_length, + int max_longest_paths) { + if (max_longest_paths * longest_path_length < subset_size) { + // If `longest_path_length` is 1, there should be one such path per node, + // and thus `max_longest_paths` should be the subset size. This branch can + // thus not be taken in this case. DCHECK_GT(longest_path_length, 1); // TODO(user): Consider using information about the number of shorter // paths to derive an even better bound. return max_longest_paths + - MathUtil::CeilOfRatio(static_cast(subset.size()) - - max_longest_paths * longest_path_length, - longest_path_length - 1); + MathUtil::CeilOfRatio( + subset_size - max_longest_paths * longest_path_length, + longest_path_length - 1); } - return MathUtil::CeilOfRatio(static_cast(subset.size()), - longest_path_length); + return MathUtil::CeilOfRatio(subset_size, longest_path_length); +} + +namespace { +struct Path { + uint32_t node_set; // Bit i is set iif node subset[i] is in the path. + int last_node; // The last node in the path. + + bool operator==(const Path& p) const { + return node_set == p.node_set && last_node == p.last_node; + } + + template + friend H AbslHashValue(H h, const Path& p) { + return H::combine(std::move(h), p.node_set, p.last_node); + } +}; + +struct PathVariableBounds { + absl::flat_hash_set incoming_arc_indices; + absl::flat_hash_map> + lower_bound_by_var_and_arc_index; +}; +} // namespace + +int MinOutgoingFlowHelper::ComputeTightMinOutgoingFlow( + absl::Span subset) { + DCHECK_GE(subset.size(), 1); + DCHECK_LE(subset.size(), 32); + DCHECK(absl::c_all_of(index_in_subset_, [](int i) { return i == -1; })); + DCHECK(absl::c_all_of(outgoing_arc_indices_, + [](const auto& v) { return v.empty(); })); + + std::vector longest_path_length_by_end_node(subset.size(), 1); + for (int i = 0; i < subset.size(); ++i) { + index_in_subset_[subset[i]] = i; + } + for (int i = 0; i < tails_.size(); ++i) { + if (index_in_subset_[tails_[i]] != -1 && + index_in_subset_[heads_[i]] != -1 && heads_[i] != tails_[i]) { + outgoing_arc_indices_[tails_[i]].push_back(i); + } + } + + absl::flat_hash_map path_var_bounds; + std::vector paths; + std::vector next_paths; + for (int i = 0; i < subset.size(); ++i) { + paths.push_back( + {.node_set = static_cast(1 << i), .last_node = subset[i]}); + } + int longest_path_length = 1; + for (int path_length = 1; path_length <= subset.size(); ++path_length) { + for (const Path& path : paths) { + // Merge the bounds by variable and arc incoming to the last node of the + // path into bounds by variable, if possible, and check whether they are + // feasible or not. + const auto& var_bounds = path_var_bounds[path]; + absl::flat_hash_map lower_bound_by_var; + for (const auto& [var, lower_bound_by_arc_index] : + var_bounds.lower_bound_by_var_and_arc_index) { + // If each arc which can reach the last node of the path enforces some + // lower bound for `var`, then the lower bound of `var` can be increased + // to the minimum of these arc-specific lower bounds (since at least one + // of these arcs must be selected to reach this node). + if (lower_bound_by_arc_index.size() != + var_bounds.incoming_arc_indices.size()) { + continue; + } + IntegerValue lb = lower_bound_by_arc_index.begin()->second; + for (const auto& [_, lower_bound] : lower_bound_by_arc_index) { + lb = std::min(lb, lower_bound); + } + if (lb > integer_trail_.LevelZeroLowerBound(var)) { + lower_bound_by_var[var] = lb; + } + } + path_var_bounds.erase(path); + auto get_lower_bound = [&](IntegerVariable var) { + const auto it = lower_bound_by_var.find(var); + if (it != lower_bound_by_var.end()) return it->second; + return integer_trail_.LevelZeroLowerBound(var); + }; + auto get_upper_bound = [&](IntegerVariable var) { + return -get_lower_bound(NegationOf(var)); + }; + bool feasible_path = true; + for (const auto& [var, lb] : lower_bound_by_var) { + if (get_upper_bound(var) < lb) { + feasible_path = false; + break; + } + } + if (!feasible_path) continue; + // We have found a feasible path, update the path length statistics... + longest_path_length = path_length; + longest_path_length_by_end_node[index_in_subset_[path.last_node]] = + path_length; + + // ... and try to extend the path with each arc going out of its last + // node. + for (const int outgoing_arc_index : + outgoing_arc_indices_[path.last_node]) { + const int head = heads_[outgoing_arc_index]; + const int head_index_in_subset = index_in_subset_[head]; + DCHECK_NE(head_index_in_subset, -1); + if (path.node_set & (1 << head_index_in_subset)) continue; + const Path new_path = { + .node_set = path.node_set | (1 << head_index_in_subset), + .last_node = head}; + if (!path_var_bounds.contains(new_path)) { + next_paths.push_back(new_path); + } + auto& new_var_bounds = path_var_bounds[new_path]; + new_var_bounds.incoming_arc_indices.insert(outgoing_arc_index); + auto update_lower_bound_by_var_and_arc_index = + [&](IntegerVariable var, int arc_index, IntegerValue lb) { + auto& lower_bound_by_arc_index = + new_var_bounds.lower_bound_by_var_and_arc_index[var]; + auto it = lower_bound_by_arc_index.find(arc_index); + if (it != lower_bound_by_arc_index.end()) { + it->second = std::max(it->second, lb); + } else { + lower_bound_by_arc_index[arc_index] = lb; + } + }; + auto update_upper_bound_by_var_and_arc_index = + [&](IntegerVariable var, int arc_index, IntegerValue ub) { + update_lower_bound_by_var_and_arc_index(NegationOf(var), + arc_index, -ub); + }; + auto update_var_bounds = [&](int arc_index, LinearTerm a, LinearTerm b, + IntegerValue lhs, IntegerValue rhs) { + if (a.coeff == 0) return; + a.MakeCoeffPositive(); + b.MakeCoeffPositive(); + // lb(b.y) <= b.y <= ub(b.y) and lhs <= a.x + b.y <= rhs imply + // ceil((lhs - ub(b.y)) / a) <= x <= floor((rhs - lb(b.y)) / a) + lhs = lhs - b.coeff * get_upper_bound(b.var); + rhs = rhs - b.coeff * get_lower_bound(b.var); + update_lower_bound_by_var_and_arc_index( + a.var, arc_index, MathUtil::CeilOfRatio(lhs, a.coeff)); + update_upper_bound_by_var_and_arc_index( + a.var, arc_index, MathUtil::FloorOfRatio(rhs, a.coeff)); + }; + const Literal lit = literals_[outgoing_arc_index]; + for (const int relation_index : + binary_relation_repository_.relation_indices(lit)) { + const auto& r = binary_relation_repository_.relation(relation_index); + update_var_bounds(outgoing_arc_index, r.a, r.b, r.lhs, r.rhs); + update_var_bounds(outgoing_arc_index, r.b, r.a, r.lhs, r.rhs); + } + } + } + std::swap(paths, next_paths); + next_paths.clear(); + } + + int max_longest_paths = 0; + for (int i = 0; i < subset.size(); ++i) { + if (longest_path_length_by_end_node[i] == longest_path_length) { + ++max_longest_paths; + } + } + // Reset the temporary data structures for the next call. + for (const int n : subset) { + index_in_subset_[n] = -1; + outgoing_arc_indices_[n].clear(); + } + return GetMinOutgoingFlow(subset.size(), longest_path_length, + max_longest_paths); } namespace { class OutgoingCutHelper { public: - OutgoingCutHelper(int num_nodes, int64_t capacity, + OutgoingCutHelper(int num_nodes, bool is_route_constraint, int64_t capacity, absl::Span demands, - const std::vector& tails, - const std::vector& heads, - const std::vector& literals, - const std::vector& literal_lp_values, - const std::vector& relevant_arcs, - LinearConstraintManager* manager, Model* model) + absl::Span tails, absl::Span heads, + absl::Span literals, Model* model) : num_nodes_(num_nodes), + is_route_constraint_(is_route_constraint), capacity_(capacity), - demands_(demands), - tails_(tails), - heads_(heads), - literals_(literals), - literal_lp_values_(literal_lp_values), - relevant_arcs_(relevant_arcs), - manager_(manager), + demands_(demands.begin(), demands.end()), + tails_(tails.begin(), tails.end()), + heads_(heads.begin(), heads.end()), + literals_(literals.begin(), literals.end()), + literal_lp_values_(literals.size()), + params_(*model->GetOrCreate()), + trail_(*model->GetOrCreate()), random_(model->GetOrCreate()), encoder_(model->GetOrCreate()), - routing_cut_subset_size_for_binary_relation_bound_( - model->GetOrCreate() - ->routing_cut_subset_size_for_binary_relation_bound()), in_subset_(num_nodes, false), - min_outgoing_flow_helper_(num_nodes, tails, heads, literals, model) { + min_outgoing_flow_helper_(num_nodes, tails_, heads_, literals_, model) { // Compute the total demands in order to know the minimum incoming/outgoing // flow. for (const int64_t demand : demands) total_demand_ += demand; complement_of_subset_.reserve(num_nodes_); } + int num_nodes() const { return num_nodes_; } + + bool is_route_constraint() const { return is_route_constraint_; } + + // Returns the arcs computed by InitializeForNewLpSolution(). + absl::Span relevant_arcs() const { + return relevant_arcs_; + } + + // Returns the arcs computed in InitializeForNewLpSolution(), sorted by + // decreasing lp_value. + absl::Span> ordered_arcs() const { + return ordered_arcs_; + } + + // This should be called each time a new LP solution is available, before + // using the other methods. + void InitializeForNewLpSolution(LinearConstraintManager* manager); + + // Merge the relevant arcs (tail, head) and (head, tail) into a single (tail, + // head) arc, with the sum of their lp_values. + absl::Span SymmetrizedRelevantArcs() { + symmetrized_relevant_arcs_ = relevant_arcs_; + SymmetrizeArcs(&symmetrized_relevant_arcs_); + return symmetrized_relevant_arcs_; + } + // Try to add an outgoing cut from the given subset. - bool TrySubsetCut(std::string name, absl::Span subset); + bool TrySubsetCut(std::string name, absl::Span subset, + LinearConstraintManager* manager); // If we look at the symmetrized version (tail <-> head = tail->head + // head->tail) and we split all the edges between a subset of nodes S and the @@ -333,39 +454,109 @@ class OutgoingCutHelper { // to consider. bool TryBlossomSubsetCut(std::string name, absl::Span symmetrized_edges, - absl::Span subset); + absl::Span subset, + LinearConstraintManager* manager); private: + // Removes the arcs with a literal fixed at false at level zero. This is + // especially useful to remove fixed self loop. + void FilterFalseArcsAtLevelZero(); + // Add a cut of the form Sum_{outgoing arcs from S} lp >= rhs_lower_bound. // // Note that we used to also add the same cut for the incoming arcs, but // because of flow conservation on these problems, the outgoing flow is always // the same as the incoming flow, so adding this extra cut doesn't seem // relevant. - bool AddOutgoingCut(std::string name, int subset_size, - const std::vector& in_subset, + bool AddOutgoingCut(LinearConstraintManager* manager, std::string name, + int subset_size, const std::vector& in_subset, int64_t rhs_lower_bound, int ignore_arcs_with_head); const int num_nodes_; + const bool is_route_constraint_; const int64_t capacity_; - const absl::Span demands_; - const std::vector& tails_; - const std::vector& heads_; - const std::vector& literals_; - const std::vector& literal_lp_values_; - const std::vector& relevant_arcs_; - LinearConstraintManager* manager_; + std::vector demands_; + std::vector tails_; + std::vector heads_; + std::vector literals_; + std::vector literal_lp_values_; + std::vector relevant_arcs_; + std::vector symmetrized_relevant_arcs_; + std::vector> ordered_arcs_; + + const SatParameters& params_; + const Trail& trail_; ModelRandomGenerator* random_; IntegerEncoder* encoder_; - const int routing_cut_subset_size_for_binary_relation_bound_; int64_t total_demand_ = 0; std::vector in_subset_; std::vector complement_of_subset_; + + MaxBoundedSubsetSum max_bounded_subset_sum_; + MaxBoundedSubsetSumExact max_bounded_subset_sum_exact_; MinOutgoingFlowHelper min_outgoing_flow_helper_; }; -bool OutgoingCutHelper::AddOutgoingCut(std::string name, int subset_size, +void OutgoingCutHelper::FilterFalseArcsAtLevelZero() { + if (trail_.CurrentDecisionLevel() != 0) return; + + int new_size = 0; + const int size = static_cast(tails_.size()); + const VariablesAssignment& assignment = trail_.Assignment(); + for (int i = 0; i < size; ++i) { + if (assignment.LiteralIsFalse(literals_[i])) continue; + tails_[new_size] = tails_[i]; + heads_[new_size] = heads_[i]; + literals_[new_size] = literals_[i]; + ++new_size; + } + if (new_size < size) { + tails_.resize(new_size); + heads_.resize(new_size); + literals_.resize(new_size); + literal_lp_values_.resize(new_size); + } +} + +void OutgoingCutHelper::InitializeForNewLpSolution( + LinearConstraintManager* manager) { + FilterFalseArcsAtLevelZero(); + + // We will collect only the arcs with a positive lp_values to speed up some + // computation below. + relevant_arcs_.clear(); + + // Sort the arcs by non-increasing lp_values. + const auto& lp_values = manager->LpValues(); + std::vector> relevant_arc_by_decreasing_lp_values; + for (int i = 0; i < literals_.size(); ++i) { + double lp_value; + const IntegerVariable direct_view = encoder_->GetLiteralView(literals_[i]); + if (direct_view != kNoIntegerVariable) { + lp_value = lp_values[direct_view]; + } else { + lp_value = + 1.0 - lp_values[encoder_->GetLiteralView(literals_[i].Negated())]; + } + literal_lp_values_[i] = lp_value; + + if (lp_value < 1e-6) continue; + relevant_arcs_.push_back({tails_[i], heads_[i], lp_value}); + relevant_arc_by_decreasing_lp_values.push_back({lp_value, i}); + } + std::sort(relevant_arc_by_decreasing_lp_values.begin(), + relevant_arc_by_decreasing_lp_values.end(), + std::greater>()); + + ordered_arcs_.clear(); + for (const auto& [score, arc] : relevant_arc_by_decreasing_lp_values) { + ordered_arcs_.push_back({tails_[arc], heads_[arc]}); + } +} + +bool OutgoingCutHelper::AddOutgoingCut(LinearConstraintManager* manager, + std::string name, int subset_size, const std::vector& in_subset, int64_t rhs_lower_bound, int ignore_arcs_with_head) { @@ -444,24 +635,21 @@ bool OutgoingCutHelper::AddOutgoingCut(std::string name, int subset_size, } } - return manager_->AddCut(outgoing.Build(), name); + return manager->AddCut(outgoing.Build(), name); } bool OutgoingCutHelper::TrySubsetCut(std::string name, - absl::Span subset) { + absl::Span subset, + LinearConstraintManager* manager) { DCHECK_GE(subset.size(), 1); DCHECK_LT(subset.size(), num_nodes_); - // These fields will be left untouched if demands.empty(). + // Initialize "in_subset" and contain_depot. bool contain_depot = false; - int64_t subset_demand = 0; - - // Initialize "in_subset" and the subset demands. for (const int n : subset) { in_subset_[n] = true; - if (!demands_.empty()) { - if (n == 0) contain_depot = true; - subset_demand += demands_[n]; + if (n == 0 && is_route_constraint_) { + contain_depot = true; } } @@ -478,63 +666,149 @@ bool OutgoingCutHelper::TrySubsetCut(std::string name, // TODO(user): This is still not as good as the "capacity" bounds below in // some cases. Fix! we should be able to use the same relation to infer the // capacity bounds somehow. - if ((contain_depot ? num_nodes_ - subset.size() : subset.size()) <= - routing_cut_subset_size_for_binary_relation_bound_) { - int automatic_bound; + const int subset_or_complement_size = + contain_depot ? num_nodes_ - subset.size() : subset.size(); + if (subset_or_complement_size <= + params_.routing_cut_subset_size_for_binary_relation_bound() && + subset_or_complement_size >= + params_.routing_cut_subset_size_for_tight_binary_relation_bound()) { + int bound; if (contain_depot) { complement_of_subset_.clear(); for (int i = 0; i < num_nodes_; ++i) { if (!in_subset_[i]) complement_of_subset_.push_back(i); } - automatic_bound = min_outgoing_flow_helper_.ComputeMinOutgoingFlow( + bound = min_outgoing_flow_helper_.ComputeMinOutgoingFlow( complement_of_subset_); } else { - automatic_bound = - min_outgoing_flow_helper_.ComputeMinOutgoingFlow(subset); + bound = min_outgoing_flow_helper_.ComputeMinOutgoingFlow(subset); } - if (automatic_bound > min_outgoing_flow) { + if (bound > min_outgoing_flow) { absl::StrAppend(&name, "Automatic"); - min_outgoing_flow = automatic_bound; + min_outgoing_flow = bound; + } + } + if (subset_or_complement_size < + params_.routing_cut_subset_size_for_tight_binary_relation_bound()) { + int bound; + if (contain_depot) { + complement_of_subset_.clear(); + for (int i = 0; i < num_nodes_; ++i) { + if (!in_subset_[i]) complement_of_subset_.push_back(i); + } + bound = min_outgoing_flow_helper_.ComputeTightMinOutgoingFlow( + complement_of_subset_); + } else { + bound = min_outgoing_flow_helper_.ComputeTightMinOutgoingFlow(subset); + } + if (bound > min_outgoing_flow) { + absl::StrAppend(&name, "AutomaticTight"); + min_outgoing_flow = bound; } } // Bounds coming from the demands_/capacity_ fields (if set). std::vector to_ignore_candidates; if (!demands_.empty()) { - const int64_t capacity_bound = - contain_depot - ? MathUtil::CeilOfRatio(total_demand_ - subset_demand, capacity_) - : MathUtil::CeilOfRatio(subset_demand, capacity_); - if (capacity_bound > min_outgoing_flow) { - min_outgoing_flow = capacity_bound; - absl::StrAppend(&name, "Capacity"); + // If subset contains depot, we actually look at the subset complement to + // derive a bound on the outgoing flow. If we cannot reach the capacity + // given the demands in the subset, we can derive tighter bounds. + int64_t has_excessive_demands = false; + int64_t has_negative_demands = false; + int64_t sum_of_elements = 0; + std::vector elements; + const auto process_demand = [&](int64_t d) { + if (d < 0) has_negative_demands = true; + if (d > capacity_) has_excessive_demands = true; + sum_of_elements += d; + elements.push_back(d); + }; + if (contain_depot) { + for (int n = 0; n < num_nodes_; ++n) { + if (!in_subset_[n]) { + process_demand(demands_[n]); + } + } + } else { + for (const int n : subset) { + process_demand(demands_[n]); + } } - // It is possible to make this cut stronger, using similar reasoning to the - // Multistar CVRP cuts: if there is a node n (other than the depot) outside - // of `subset`, with a demand that causes subset_demand + n_demand to - // strictly increase the `capacity_bound` above, then we still need - // `capacity_bound` on the outgoing arcs of `subset` other than those - // towards n (noted A in the diagram below): + // Lets wait for these to disappear before adding cuts. + if (has_excessive_demands) return false; + + // Try to tighten the capacity using DP. Note that there is no point doing + // anything if one route can serve all demands since then the bound is + // already tight. // - // ^ ^ - // | | - // ---------------- - // <--| subset | n |--> - // ---------------- - // | | - // v v - // - // \------A------/\--B--/ - // - // By hypothesis, outgoing_flow(A) + outgoing_flow(B) > capacity_bound - // and, since n is not the depot, outgoing_flow(B) <= 1. Hence - // outgoing_flow(A) >= capacity_bound. - if (!contain_depot && capacity_bound >= min_outgoing_flow) { + // TODO(user): Compute a bound in the presence of negative demands? + bool exact_was_used = false; + int64_t tightened_capacity = capacity_; + if (!has_negative_demands && sum_of_elements > capacity_) { + max_bounded_subset_sum_.Reset(capacity_); + for (const int64_t e : elements) { + max_bounded_subset_sum_.Add(e); + } + tightened_capacity = max_bounded_subset_sum_.CurrentMax(); + + // If the complexity looks ok, try a more expensive DP than the quick one + // above. + if (max_bounded_subset_sum_exact_.ComplexityEstimate( + elements.size(), capacity_) < params_.routing_cut_dp_effort()) { + const int64_t exact = + max_bounded_subset_sum_exact_.MaxSubsetSum(elements, capacity_); + CHECK_LE(exact, tightened_capacity); + if (exact < tightened_capacity) { + tightened_capacity = exact; + exact_was_used = true; + } + } + } + + const int64_t flow_lower_bound = + MathUtil::CeilOfRatio(sum_of_elements, tightened_capacity); + if (flow_lower_bound > min_outgoing_flow) { + min_outgoing_flow = flow_lower_bound; + absl::StrAppend(&name, exact_was_used ? "Tightened" : "Capacity"); + } + + if (!contain_depot && flow_lower_bound >= min_outgoing_flow) { + // We compute the biggest extra item that could fit in 'flow_lower_bound' + // bins. If the first (flow_lower_bound - 1) bins are tight, i.e. all + // their tightened_capacity is filled, then the last bin will have + // 'last_bin_fillin' stuff, which will leave 'space_left' to fit an extra + // 'item. + const int64_t last_bin_fillin = + sum_of_elements - (flow_lower_bound - 1) * tightened_capacity; + const int64_t space_left = capacity_ - last_bin_fillin; + DCHECK_GE(space_left, 0); + DCHECK_LT(space_left, capacity_); + + // It is possible to make this cut stronger, using similar reasoning to + // the Multistar CVRP cuts: if there is a node n (other than the depot) + // outside of `subset`, with a demand that is greater than space_left, + // then the outgoing flow of (subset + n) is >= flow_lower_bound + 1. + // Using this, we can show that we still need `flow_lower_bound` on the + // outgoing arcs of `subset` other than those towards n (noted A in the + // diagram below): + // + // ^ ^ + // | | + // ---------------- + // <--| subset | n |--> + // ---------------- + // | | + // v v + // + // \------A------/\--B--/ + // + // By hypothesis, outgoing_flow(A) + outgoing_flow(B) > flow_lower_bound + // and, since n is not the depot, outgoing_flow(B) <= 1. Hence + // outgoing_flow(A) >= flow_lower_bound. for (int n = 1; n < num_nodes_; ++n) { if (in_subset_[n]) continue; - if (MathUtil::CeilOfRatio(subset_demand + demands_[n], capacity_) > - capacity_bound) { + if (demands_[n] > space_left) { to_ignore_candidates.push_back(n); } } @@ -599,7 +873,7 @@ bool OutgoingCutHelper::TrySubsetCut(std::string name, // Add a cut if the current outgoing flow is not enough. bool result = false; if (outgoing_flow + 1e-2 < min_outgoing_flow) { - result = AddOutgoingCut(name, subset.size(), in_subset_, + result = AddOutgoingCut(manager, name, subset.size(), in_subset_, /*rhs_lower_bound=*/min_outgoing_flow, ignore_arcs_with_head); } @@ -612,7 +886,7 @@ bool OutgoingCutHelper::TrySubsetCut(std::string name, bool OutgoingCutHelper::TryBlossomSubsetCut( std::string name, absl::Span symmetrized_edges, - absl::Span subset) { + absl::Span subset, LinearConstraintManager* manager) { DCHECK_GE(subset.size(), 1); DCHECK_LT(subset.size(), num_nodes_); @@ -736,7 +1010,7 @@ bool OutgoingCutHelper::TryBlossomSubsetCut( builder.AddConstant(IntegerValue(num_actual_inverted)); if (num_actual_inverted % 2 == 0) return false; - return manager_->AddCut(builder.Build(), name); + return manager->AddCut(builder.Build(), name); } } // namespace @@ -917,60 +1191,26 @@ void SymmetrizeArcs(std::vector* arcs) { // // Note that this is mainly a "symmetric" case algo, but it does still work for // the asymmetric case. -void SeparateSubtourInequalities( - int num_nodes, const std::vector& tails, const std::vector& heads, - const std::vector& literals, absl::Span demands, - int64_t capacity, LinearConstraintManager* manager, Model* model) { +void SeparateSubtourInequalities(OutgoingCutHelper& helper, + LinearConstraintManager* manager) { + const int num_nodes = helper.num_nodes(); if (num_nodes <= 2) return; - // We will collect only the arcs with a positive lp_values to speed up some - // computation below. - std::vector relevant_arcs; + helper.InitializeForNewLpSolution(manager); - // Sort the arcs by non-increasing lp_values. - const auto& lp_values = manager->LpValues(); - std::vector literal_lp_values(literals.size()); - std::vector> arc_by_decreasing_lp_values; - auto* encoder = model->GetOrCreate(); - for (int i = 0; i < literals.size(); ++i) { - double lp_value; - const IntegerVariable direct_view = encoder->GetLiteralView(literals[i]); - if (direct_view != kNoIntegerVariable) { - lp_value = lp_values[direct_view]; - } else { - lp_value = - 1.0 - lp_values[encoder->GetLiteralView(literals[i].Negated())]; - } - literal_lp_values[i] = lp_value; - - if (lp_value < 1e-6) continue; - relevant_arcs.push_back({tails[i], heads[i], lp_value}); - arc_by_decreasing_lp_values.push_back({lp_value, i}); - } - std::sort(arc_by_decreasing_lp_values.begin(), - arc_by_decreasing_lp_values.end(), - std::greater>()); - - std::vector> ordered_arcs; - for (const auto& [score, arc] : arc_by_decreasing_lp_values) { - ordered_arcs.push_back({tails[arc], heads[arc]}); - } std::vector subset_data; std::vector> subsets; - GenerateInterestingSubsets(num_nodes, ordered_arcs, + GenerateInterestingSubsets(num_nodes, helper.ordered_arcs(), /*stop_at_num_components=*/2, &subset_data, &subsets); const int depot = 0; - if (!demands.empty()) { + if (helper.is_route_constraint()) { // Add the depot so that we have a trivial bound on the number of // vehicle. subsets.push_back(absl::MakeSpan(&depot, 1)); } - OutgoingCutHelper helper(num_nodes, capacity, demands, tails, heads, literals, - literal_lp_values, relevant_arcs, manager, model); - // Hack/optim: we exploit the tree structure of the subsets to not add a cut // for a larger subset if we added a cut from one included in it. // @@ -984,7 +1224,7 @@ void SeparateSubtourInequalities( if (subset.size() <= 1) continue; const int start = static_cast(subset.data() - subset_data.data()); if (start <= last_added_start) continue; - if (helper.TrySubsetCut("Circuit", subset)) { + if (helper.TrySubsetCut("Circuit", subset, manager)) { ++num_added; last_added_start = start; } @@ -1005,8 +1245,10 @@ void SeparateSubtourInequalities( // the course of the algorithm. This could also be interesting. But it is // hard to tell with our current benchmark setup. if (num_added != 0) return; - SymmetrizeArcs(&relevant_arcs); - std::vector parent = ComputeGomoryHuTree(num_nodes, relevant_arcs); + absl::Span symmetrized_relevant_arcs = + helper.SymmetrizedRelevantArcs(); + std::vector parent = + ComputeGomoryHuTree(num_nodes, symmetrized_relevant_arcs); // Try all interesting subset from the Gomory-Hu tree. ExtractAllSubsetsFromForest(parent, &subset_data, &subsets); @@ -1016,7 +1258,7 @@ void SeparateSubtourInequalities( if (subset.size() == num_nodes) continue; const int start = static_cast(subset.data() - subset_data.data()); if (start <= last_added_start) continue; - if (helper.TrySubsetCut("CircuitExact", subset)) { + if (helper.TrySubsetCut("CircuitExact", subset, manager)) { ++num_added; last_added_start = start; } @@ -1025,13 +1267,11 @@ void SeparateSubtourInequalities( // Exact separation of symmetric Blossom cut. We use the algorithm in the // paper: "A Faster Exact Separation Algorithm for Blossom Inequalities", Adam // N. Letchford, Gerhard Reinelt, Dirk Oliver Theis, 2004. - // - // Note that the "relevant_arcs" were symmetrized above. if (num_added != 0) return; if (num_nodes <= 2) return; std::vector for_blossom; - for_blossom.reserve(relevant_arcs.size()); - for (ArcWithLpValue arc : relevant_arcs) { + for_blossom.reserve(symmetrized_relevant_arcs.size()); + for (ArcWithLpValue arc : symmetrized_relevant_arcs) { if (arc.lp_value > 0.5) arc.lp_value = 1.0 - arc.lp_value; if (arc.lp_value < 1e-6) continue; for_blossom.push_back(arc); @@ -1044,7 +1284,8 @@ void SeparateSubtourInequalities( if (subset.size() == num_nodes) continue; const int start = static_cast(subset.data() - subset_data.data()); if (start <= last_added_start) continue; - if (helper.TryBlossomSubsetCut("CircuitBlossom", relevant_arcs, subset)) { + if (helper.TryBlossomSubsetCut("CircuitBlossom", symmetrized_relevant_arcs, + subset, manager)) { ++num_added; last_added_start = start; } @@ -1070,62 +1311,42 @@ std::vector GetAssociatedVariables( return result; } -// This is especially useful to remove fixed self loop. -void FilterFalseArcsAtLevelZero(std::vector& tails, - std::vector& heads, - std::vector& literals, Model* model) { - const Trail& trail = *model->GetOrCreate(); - if (trail.CurrentDecisionLevel() != 0) return; - - int new_size = 0; - const int size = static_cast(tails.size()); - const VariablesAssignment& assignment = trail.Assignment(); - for (int i = 0; i < size; ++i) { - if (assignment.LiteralIsFalse(literals[i])) continue; - tails[new_size] = tails[i]; - heads[new_size] = heads[i]; - literals[new_size] = literals[i]; - ++new_size; - } - if (new_size < size) { - tails.resize(new_size); - heads.resize(new_size); - literals.resize(new_size); - } -} - } // namespace // We use a basic algorithm to detect components that are not connected to the // rest of the graph in the LP solution, and add cuts to force some arcs to // enter and leave this component from outside. CutGenerator CreateStronglyConnectedGraphCutGenerator( - int num_nodes, std::vector tails, std::vector heads, - std::vector literals, Model* model) { + int num_nodes, absl::Span tails, absl::Span heads, + absl::Span literals, Model* model) { + auto helper = std::make_unique( + num_nodes, /*is_route_constraint=*/false, /*capacity=*/0, + /*demands=*/absl::Span(), tails, heads, literals, model); CutGenerator result; result.vars = GetAssociatedVariables(literals, model); - result.generate_cuts = [=](LinearConstraintManager* manager) mutable { - FilterFalseArcsAtLevelZero(tails, heads, literals, model); - SeparateSubtourInequalities(num_nodes, tails, heads, literals, - /*demands=*/{}, /*capacity=*/0, manager, model); - return true; - }; + result.generate_cuts = + [helper = std::move(helper)](LinearConstraintManager* manager) { + SeparateSubtourInequalities(*helper, manager); + return true; + }; return result; } -CutGenerator CreateCVRPCutGenerator(int num_nodes, std::vector tails, - std::vector heads, - std::vector literals, - std::vector demands, +CutGenerator CreateCVRPCutGenerator(int num_nodes, absl::Span tails, + absl::Span heads, + absl::Span literals, + absl::Span demands, int64_t capacity, Model* model) { + auto helper = std::make_unique( + num_nodes, /*is_route_constraint=*/true, capacity, demands, tails, heads, + literals, model); CutGenerator result; result.vars = GetAssociatedVariables(literals, model); - result.generate_cuts = [=](LinearConstraintManager* manager) mutable { - FilterFalseArcsAtLevelZero(tails, heads, literals, model); - SeparateSubtourInequalities(num_nodes, tails, heads, literals, demands, - capacity, manager, model); - return true; - }; + result.generate_cuts = + [helper = std::move(helper)](LinearConstraintManager* manager) { + SeparateSubtourInequalities(*helper, manager); + return true; + }; return result; } @@ -1159,7 +1380,7 @@ void SeparateFlowInequalities( for (int i = 0; i < arc_capacities.size(); ++i) { const double lp_value = arc_capacities[i].LpValue(lp_values); if (!arc_capacities[i].IsConstant()) { - gcd = MathUtil::GCD64(gcd, std::abs(arc_capacities[i].coeff.value())); + gcd = std::gcd(gcd, std::abs(arc_capacities[i].coeff.value())); } if (lp_value < 1e-6 && arc_capacities[i].constant == 0) continue; relevant_arcs.push_back( diff --git a/ortools/sat/routing_cuts.h b/ortools/sat/routing_cuts.h index aadfb007cc..915abd3217 100644 --- a/ortools/sat/routing_cuts.h +++ b/ortools/sat/routing_cuts.h @@ -45,9 +45,23 @@ class MinOutgoingFlowHelper { // Returns the minimum flow going out of `subset`, based on a conservative // estimate of the maximum number of nodes of a feasible path inside this // subset. `subset` must not be empty and must not contain the depot (node 0). + // Paths are approximated by their length and their last node, and can thus + // contain cycles. The complexity is O(subset.size() ^ 3). int ComputeMinOutgoingFlow(absl::Span subset); + // Same as above, but uses less conservative estimates (paths are approximated + // by their set of nodes and their last node -- hence they can't contain + // cycles). The complexity is O(2 ^ subset.size()). + int ComputeTightMinOutgoingFlow(absl::Span subset); + private: + // Returns the minimum flow going out of a subset of size `subset_size`, + // assuming that the longest feasible path inside this subset has + // `longest_path_length` nodes and that there are at most `max_longest_paths` + // such paths. + int GetMinOutgoingFlow(int subset_size, int longest_path_length, + int max_longest_paths); + const std::vector& tails_; const std::vector& heads_; const std::vector& literals_; @@ -66,9 +80,11 @@ class MinOutgoingFlowHelper { // position i+1. std::vector in_subset_; - // For each node n, the indices (in tails_, heads_) of the m->n arcs inside - // the subset (self arcs excepted). + std::vector index_in_subset_; + // For each node n, the indices (in tails_, heads_) of the m->n and n->m arcs + // inside the subset (self arcs excepted). std::vector> incoming_arc_indices_; + std::vector> outgoing_arc_indices_; // For each node n, whether it can appear at the current and next position in // a feasible path. std::vector reachable_; @@ -171,17 +187,17 @@ std::vector ComputeGomoryHuTree( // connected. Note that we already assume basic constraint to be in the lp, so // we do not add any cuts for components of size 1. CutGenerator CreateStronglyConnectedGraphCutGenerator( - int num_nodes, std::vector tails, std::vector heads, - std::vector literals, Model* model); + int num_nodes, absl::Span tails, absl::Span heads, + absl::Span literals, Model* model); // Almost the same as CreateStronglyConnectedGraphCutGenerator() but for each // components, computes the demand needed to serves it, and depending on whether // it contains the depot (node zero) or not, compute the minimum number of // vehicle that needs to cross the component border. -CutGenerator CreateCVRPCutGenerator(int num_nodes, std::vector tails, - std::vector heads, - std::vector literals, - std::vector demands, +CutGenerator CreateCVRPCutGenerator(int num_nodes, absl::Span tails, + absl::Span heads, + absl::Span literals, + absl::Span demands, int64_t capacity, Model* model); // Try to find a subset where the current LP capacity of the outgoing or diff --git a/ortools/sat/routing_cuts_test.cc b/ortools/sat/routing_cuts_test.cc index edcdecf938..2862ecb59d 100644 --- a/ortools/sat/routing_cuts_test.cc +++ b/ortools/sat/routing_cuts_test.cc @@ -45,6 +45,22 @@ namespace { using ::testing::ElementsAre; +TEST(MinOutgoingFlowHelperTest, TwoNodesWithoutConstraints) { + Model model; + const std::vector tails = {0, 1}; + const std::vector heads = {1, 0}; + const std::vector literals = { + Literal(model.Add(NewBooleanVariable()), true), + Literal(model.Add(NewBooleanVariable()), true)}; + MinOutgoingFlowHelper helper(2, tails, heads, literals, &model); + + const int min_flow = helper.ComputeMinOutgoingFlow({0, 1}); + const int tight_min_flow = helper.ComputeTightMinOutgoingFlow({0, 1}); + + EXPECT_EQ(min_flow, 1); + EXPECT_EQ(tight_min_flow, 1); +} + TEST(MinOutgoingFlowHelperTest, CapacityConstraints) { Model model; const int num_nodes = 5; @@ -84,6 +100,7 @@ TEST(MinOutgoingFlowHelperTest, CapacityConstraints) { MinOutgoingFlowHelper helper(num_nodes, tails, heads, literals, &model); const int min_flow = helper.ComputeMinOutgoingFlow({1, 2, 3, 4}); + const int tight_min_flow = helper.ComputeTightMinOutgoingFlow({1, 2, 3, 4}); // Due to the capacity constraints, a feasible path can have at most 3 nodes, // hence at least two paths are needed. The lower bound of the vehicle load @@ -96,6 +113,7 @@ TEST(MinOutgoingFlowHelperTest, CapacityConstraints) { // 3 | 0 13 24 - // 4 | 0 14 24 - EXPECT_EQ(min_flow, 2); + EXPECT_EQ(tight_min_flow, 2); } TEST(MinOutgoingFlowHelperTest, TimeWindows) { @@ -136,6 +154,7 @@ TEST(MinOutgoingFlowHelperTest, TimeWindows) { MinOutgoingFlowHelper helper(num_nodes, tails, heads, literals, &model); const int min_flow = helper.ComputeMinOutgoingFlow({1, 2, 3, 4}); + const int tight_min_flow = helper.ComputeTightMinOutgoingFlow({1, 2, 3, 4}); // Due to the time window constraints, a feasible path can have at most 3 // nodes, hence at least two paths are needed. The earliest departure time @@ -148,6 +167,7 @@ TEST(MinOutgoingFlowHelperTest, TimeWindows) { // 3 | 18 18 - - // 4 | 28 28 28 - EXPECT_EQ(min_flow, 2); + EXPECT_EQ(tight_min_flow, 2); } // Test on a simple tree: @@ -473,7 +493,7 @@ TEST(CreateFlowCutGeneratorTest, BasicExample) { *min_incoming_flow = std::max(IntegerValue(0), demand); *min_outgoing_flow = std::max(IntegerValue(0), -demand); }; - const CutGenerator generator = CreateFlowCutGenerator( + CutGenerator generator = CreateFlowCutGenerator( num_nodes, tails, heads, capacities, get_flows, &model); LinearConstraintManager manager(&model); @@ -517,7 +537,7 @@ TEST(CreateFlowCutGeneratorTest, WithMinusOneArcs) { *min_incoming_flow = std::max(IntegerValue(0), demand); *min_outgoing_flow = std::max(IntegerValue(0), -demand); }; - const CutGenerator generator = CreateFlowCutGenerator( + CutGenerator generator = CreateFlowCutGenerator( num_nodes, tails, heads, capacities, get_flows, &model); LinearConstraintManager manager(&model); diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 33b586e2bf..25f0d5a60f 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -23,7 +23,7 @@ option java_multiple_files = true; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 313 +// NEXT TAG: 315 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -914,12 +914,26 @@ message SatParameters { // If the size of a subset of nodes of a RoutesConstraint is less than this // value, use linear constraints of size 1 and 2 (such as capacity and time // window constraints) enforced by the arc literals to compute cuts for this - // subset. The algorithm for these cuts has a O(n^3) complexity, where n is - // the subset size. Hence the value of this parameter should not be too large - // (e.g. 10 or 20). + // subset (unless the subset size is less than + // routing_cut_subset_size_for_tight_binary_relation_bound, in which case the + // corresponding algorithm is used instead). The algorithm for these cuts has + // a O(n^3) complexity, where n is the subset size. Hence the value of this + // parameter should not be too large (e.g. 10 or 20). optional int32 routing_cut_subset_size_for_binary_relation_bound = 312 [default = 0]; + // Similar to above, but with a different algorithm producing better cuts, at + // the price of a higher O(2^n) complexity, where n is the subset size. Hence + // the value of this parameter should be small (e.g. less than 10). + optional int32 routing_cut_subset_size_for_tight_binary_relation_bound = 313 + [default = 0]; + + // The amount of "effort" to spend in dynamic programming for computing + // routing cuts. This is in term of basic operations needed by the algorithm + // in the worst case, so a value like 1e8 should take less than a second to + // compute. + optional double routing_cut_dp_effort = 314 [default = 1e7]; + // The search branching will be used to decide how to branch on unfixed nodes. enum SearchBranching { // Try to fix all literals using the underlying SAT solver's heuristics, diff --git a/ortools/sat/scheduling_cuts.cc b/ortools/sat/scheduling_cuts.cc index 42064fccf0..99bef978e1 100644 --- a/ortools/sat/scheduling_cuts.cc +++ b/ortools/sat/scheduling_cuts.cc @@ -1266,49 +1266,126 @@ void GenerateShortCompletionTimeCutsWithExactBound( top_n_cuts.TransferToManager(manager); } +namespace { + +// Returns a copy of the event with the start time increased to time. +// Energy (min and decomposed) are updated accordingly. +CtEvent TrimEventAfter(IntegerValue time, const CtEvent& old_event) { + DCHECK_GT(time, old_event.x_start_min); + CtEvent event = old_event; // Copy. + event.lifted = true; + // Build the vector of energies as the vector of sizes. + event.energy_min = ComputeEnergyMinInWindow( + event.x_start_min, event.x_start_max, event.x_end_min, event.x_end_max, + event.x_size_min, event.y_size_min, event.decomposed_energy, time, + event.x_end_max); + event.x_size_min = event.x_size_min + event.x_start_min - time; + event.x_start_min = time; + if (event.energy_min > event.x_size_min * event.y_size_min) { + event.use_energy = true; + } + DCHECK_GE(event.energy_min, event.x_size_min * event.y_size_min); + return event; +} + +// Collects all possible demand values for the event, and adds them to the +// subset sum of the reachable capacity of the cumulative constraint. +void AddEventDemandsToCapacitySubsetSum( + const CtEvent& event, const VariablesAssignment& assignment, + IntegerValue capacity_max, std::vector& tmp_possible_demands, + MaxBoundedSubsetSum& dp) { + if (dp.CurrentMax() != capacity_max) { + if (event.y_size_is_fixed) { + dp.Add(event.y_size_min.value()); + } else if (!event.decomposed_energy.empty()) { + tmp_possible_demands.clear(); + for (const auto& [literal, size, demand] : event.decomposed_energy) { + if (assignment.LiteralIsFalse(literal)) continue; + if (assignment.LiteralIsTrue(literal)) { + tmp_possible_demands.assign({demand.value()}); + break; + } + tmp_possible_demands.push_back(demand.value()); + } + dp.AddChoices(tmp_possible_demands); + } else { // event.y_size_min is not fixed, we abort. + // TODO(user): We could still add all events in the range if the range + // is small. + dp.Add(capacity_max.value()); + } + } +} + +// In the no_overlap case, we have: +// area = event.x_size_min ^ 2 +// In the simple cumulative case, we split split the task (demand, size) into +// demand tasks in the no_overlap case. +// area = event.y_size_min * event.x_size_min * event.x_size_min +// In the cumulative case, the minimum energy of a task can be greater than +// the product of its size and demand (energy_min > side * demand). +// We use energy_min * size. +IntegerValue SmithRuleIndividualContribution(const CtEvent& event) { + if (event.use_energy) { + return event.energy_min * event.x_size_min; + } else { + return event.x_size_min * event.x_size_min * event.y_size_min; + } +} + +} // namespace + // We generate the cut from the Smith's rule from: // M. Queyranne, Structure of a simple scheduling polyhedron, // Mathematical Programming 58 (1993), 263–285 // // The original cut is: -// sum(end_min_i * duration_min_i) >= -// (sum(duration_min_i^2) + sum(duration_min_i)^2) / 2 +// sum(end_min_i * size_min_i) >= +// (sum(size_min_i^2) + sum(size_min_i)^2) / 2 // We strengthen this cuts by noticing that if all tasks starts after S, // then replacing end_min_i by (end_min_i - S) is still valid. // // A second difference is that we lift intervals that starts before a given // value, but are forced to cross it. // -// In the case of a cumulative constraint, we compute the cuts like if it -// was a disjunctive cut with all the tasks being `demand` tasks with -// duration `size_min`, all these spread out on `reachable_capacity` -// no_overlap machines, counted from `current_start_min`. +// In the case of a cumulative constraint with a capacity of C, we compute a +// valid equation by splitting the task (size_min si, demand_min di) into di +// tasks of size si and demand 1, that we spread across C no_overlap +// constraint. When doing so, the lhs of the equation is the same, the first +// term of the rhs is also unchanged. A lower bound of the second term of the +// rhs is reached when the split is exact (each no_overlap sees a long demand of +// sum(si * di / C). Thus the second term is greater or equal to +// (sum (si * di) ^ 2) / (2 * C) // -// Thus, for each task (size si, end ei, demand di) and capacity c: -// sum (si * di * (ei - current_start_min)) >= -// sum(di * si ^ 2) c * (sum (si * di) / c) ^ 2 -// ---------------- + --------------------------- -// 2 2 +// Sometimes, the min energy of the task i is greater than si * di. +// Let's introduce ai the minimum energy of the task and rewrite the previous +// equation. In that new setting, we can rewrite the cumulative transformation +// by splitting each tasks into at least di tasks of size at least si and demand +// 1. // -// sum (si * di * ei) - sum (si * di) * current_start_min >= -// sum(di * si ^ 2) sum (si * di) ^ 2 -// ---------------- + ----------------- -// 2 2 * c +// In that setting, the lhs is rewritten as sum(ai * ei) and the second term of +// the rhs is rewritten as sum(ai) ^ 2 / (2 * C). // -// By introducing ai the min energy of the task, we reinforce the -// constraint into: +// The question is how to rewrite the term `sum(di * si * si). The minimum +// contribution is when the task has size si and demand ai / si. (note that si +// is the minimum size of the task, and di its minimum demand). We can +// replace the small rectangle area term by ai * si. // -// sum (ai * ei) >= -// sum(di * ai) sum (ai) ^ 2 -// ------------ + ------------ + sum (ai) * current_start_min -// 2 2 * c - +// sum (ai * ei) - sum (ai) * current_start_min >= +// sum(si * ai) / 2 + (sum (ai) ^ 2) / (2 * C) +// +// The separation procedure is implemented using two loops: +// - first, we loop on potential start times in increasing order. +// - second loop, we add tasks that must contribute after this start time +// ordered by increasing end time in the LP relaxation. void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, std::vector events, IntegerValue capacity_max, bool skip_low_sizes, Model* model, LinearConstraintManager* manager) { TopNCuts top_n_cuts(5); + const VariablesAssignment& assignment = + model->GetOrCreate()->Assignment(); + std::vector tmp_possible_demands; // Sort by start min to bucketize by start_min. std::sort(events.begin(), events.end(), @@ -1316,6 +1393,8 @@ void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, return std::tie(e1.x_start_min, e1.y_size_min, e1.x_lp_end) < std::tie(e2.x_start_min, e2.y_size_min, e2.x_lp_end); }); + + // First loop: we loop on potential start times. for (int start = 0; start + 1 < events.size(); ++start) { // Skip to the next start_min value. if (start > 0 && @@ -1325,8 +1404,6 @@ void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, const IntegerValue sequence_start_min = events[start].x_start_min; std::vector residual_tasks(events.begin() + start, events.end()); - const VariablesAssignment& assignment = - model->GetOrCreate()->Assignment(); // We look at events that start before sequence_start_min, but are forced // to cross this time point. In that case, we replace this event by a @@ -1336,83 +1413,56 @@ void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, for (int before = 0; before < start; ++before) { if (events[before].x_start_min + events[before].x_size_min > sequence_start_min) { - // Build the vector of energies as the vector of sizes. - CtEvent event = events[before]; // Copy. - event.lifted = true; - event.energy_min = ComputeEnergyMinInWindow( - event.x_start_min, event.x_start_max, event.x_end_min, - event.x_end_max, event.x_size_min, event.y_size_min, - event.decomposed_energy, sequence_start_min, event.x_end_max); - event.x_size_min = - event.x_size_min + event.x_start_min - sequence_start_min; - event.x_start_min = sequence_start_min; - if (event.energy_min > event.x_size_min * event.y_size_min) { - event.use_energy = true; - } - DCHECK_GE(event.energy_min, event.x_size_min * event.y_size_min); + CtEvent event = TrimEventAfter(sequence_start_min, events[before]); if (event.energy_min <= 0) continue; - residual_tasks.push_back(event); + residual_tasks.push_back(std::move(event)); } } - std::sort(residual_tasks.begin(), residual_tasks.end(), - [](const CtEvent& e1, const CtEvent& e2) { - return e1.x_lp_end < e2.x_lp_end; - }); - // Best cut so far for this loop. int best_end = -1; double best_efficacy = 0.01; IntegerValue best_min_contrib = 0; IntegerValue best_capacity = 0; - bool loop_is_valid = true; - // Area of the large rectangle. + // Used in the first term of the rhs of the equation. + IntegerValue sum_event_contributions = 0; + // Used in the second term of the rhs of the equation. IntegerValue sum_energy = 0; // For normalization. IntegerValue sum_square_energy = 0; - // Area of the small rectangles. - IntegerValue weighted_sum_square_duration = 0; double lp_contrib = 0.0; IntegerValue current_start_min(kMaxIntegerValue); MaxBoundedSubsetSum dp(capacity_max.value()); - std::vector possible_demands; + + // We will add tasks one by one, sorted by end time, and evaluate the + // potential cut at each step. + std::sort(residual_tasks.begin(), residual_tasks.end(), + [](const CtEvent& e1, const CtEvent& e2) { + return e1.x_lp_end < e2.x_lp_end; + }); + + // Second loop: we add tasks one by one. for (int i = 0; i < residual_tasks.size(); ++i) { const CtEvent& event = residual_tasks[i]; DCHECK_GE(event.x_start_min, sequence_start_min); - const IntegerValue tmp_contrib = - CapProdI(event.x_size_min, event.energy_min); - // Make sure we do not overflow, and we do not generate bogus cuts if we - // do. - loop_is_valid = false; + // Make sure we do not overflow. if (!AddTo(event.energy_min, &sum_energy)) break; - if (!AddTo(tmp_contrib, &weighted_sum_square_duration)) break; + if (!AddTo(SmithRuleIndividualContribution(event), + &sum_event_contributions)) { + break; + } if (!AddSquareTo(event.energy_min, &sum_square_energy)) break; lp_contrib += event.x_lp_end * ToDouble(event.energy_min); current_start_min = std::min(current_start_min, event.x_start_min); - if (dp.CurrentMax() != capacity_max) { - if (event.y_size_is_fixed) { - dp.Add(event.y_size_min.value()); - } else if (!event.decomposed_energy.empty()) { - possible_demands.clear(); - for (const auto& [literal, size, demand] : event.decomposed_energy) { - if (assignment.LiteralIsFalse(literal)) continue; - if (assignment.LiteralIsTrue(literal)) { - possible_demands.assign({demand.value()}); - break; - } - possible_demands.push_back(demand.value()); - } - dp.AddChoices(possible_demands); - } else { - dp.Add(capacity_max.value()); - } - } + // Maintain the reachable capacity with a bounded complexity subset sum. + AddEventDemandsToCapacitySubsetSum(event, assignment, capacity_max, + tmp_possible_demands, dp); // This is competing with the brute force approach. Skip cases covered // by the other code. @@ -1420,6 +1470,7 @@ void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, const IntegerValue reachable_capacity = dp.CurrentMax(); + // Do we have a violated cut ? const IntegerValue large_rectangle_contrib = CapProdI(sum_energy, sum_energy); if (AtMinOrMaxInt64I(large_rectangle_contrib)) break; @@ -1427,14 +1478,12 @@ void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, CeilRatio(large_rectangle_contrib, reachable_capacity); IntegerValue min_contrib = - CapAddI(weighted_sum_square_duration, mean_large_rectangle_contrib); + CapAddI(sum_event_contributions, mean_large_rectangle_contrib); if (AtMinOrMaxInt64I(min_contrib)) break; - - // The above is the double of the area. min_contrib = CeilRatio(min_contrib, 2); // shift contribution by current_start_min. - if (!AddTo(CapProdI(sum_energy, current_start_min), &min_contrib)) break; + if (!AddProductTo(sum_energy, current_start_min, &min_contrib)) break; // The efficacy of the cut is the normalized violation of the above // equation. We will normalize by the sqrt of the sum of squared energies. @@ -1453,13 +1502,10 @@ void GenerateCompletionTimeCutsWithEnergy(absl::string_view cut_name, best_min_contrib = min_contrib; best_capacity = reachable_capacity; } - - // We got there, the loop is valid. - loop_is_valid = true; } - if (!loop_is_valid) continue; - + // We have inserted all tasks. Have we found a violated cut ? + // If so, add the most violated one to the top_n cut container. if (best_end != -1) { LinearConstraintBuilder cut(model, best_min_contrib, kMaxIntegerValue); bool is_lifted = false; diff --git a/ortools/sat/scheduling_helpers.cc b/ortools/sat/scheduling_helpers.cc index 80aa4bd1fb..b1307e2e62 100644 --- a/ortools/sat/scheduling_helpers.cc +++ b/ortools/sat/scheduling_helpers.cc @@ -371,16 +371,28 @@ bool SchedulingConstraintHelper::PropagatePrecedence(int a, int b) { if (before.coeff != 1) return true; if (after.var == kNoIntegerVariable) return true; if (before.var == kNoIntegerVariable) return true; + if (before.var == after.var) { + if (before.constant <= after.constant) { + return true; + } else { + sat_solver_->NotifyThatModelIsUnsat(); + return false; + } + } const IntegerValue offset = before.constant - after.constant; if (precedence_relations_->Add(before.var, after.var, offset)) { VLOG(2) << "new relation " << TaskDebugString(a) << " <= " << TaskDebugString(b); - - // TODO(user): Adding new constraint during propagation might not be the - // best idea as it can create some complication. - AddWeightedSumLowerOrEqual({}, {before.var, after.var}, - {int64_t{1}, int64_t{-1}}, -offset.value(), - model_); + if (before.var == NegationOf(after.var)) { + AddWeightedSumLowerOrEqual({}, {before.var}, {int64_t{2}}, + -offset.value(), model_); + } else { + // TODO(user): Adding new constraint during propagation might not be the + // best idea as it can create some complication. + AddWeightedSumLowerOrEqual({}, {before.var, after.var}, + {int64_t{1}, int64_t{-1}}, -offset.value(), + model_); + } if (sat_solver_->ModelIsUnsat()) return false; } return true; diff --git a/ortools/sat/util.cc b/ortools/sat/util.cc index 50ca686304..9d047a4f6b 100644 --- a/ortools/sat/util.cc +++ b/ortools/sat/util.cc @@ -903,5 +903,110 @@ std::vector> AtMostOneDecomposition( return decomposer.decomposition(); } +absl::Span SortedSubsetSums::Compute( + absl::Span elements, int64_t maximum_sum, + bool abort_if_maximum_reached) { + sums_.clear(); + sums_.push_back(0); + for (const int64_t e : elements) { + DCHECK_GE(e, 0); + if (e == 0 || e > maximum_sum) continue; + + // Optimization: If all the sums in [0, maximum_sum] are already reachable + // we can abort early since no new reachable sum wil be discovered. + if (sums_.size() == maximum_sum + 1) return sums_; + + // Early abort when asked if we already reached maximum_sum. + if (abort_if_maximum_reached && sums_.back() == maximum_sum) return sums_; + + // We merge sort sums_ and sums_ + e into new_sums_. + int i = 0; + int j = 0; + new_sums_.clear(); + const int size = sums_.size(); + const int64_t* const data = sums_.data(); + int64_t last_pushed = -1; + while (i < size) { + DCHECK_LE(j, i); // Since we add a positive e. + const int64_t a = data[i]; + const int64_t b = data[j] + e; + int64_t to_push; + if (a <= b) { + ++i; + if (a == b) ++j; + to_push = a; + } else { + to_push = b; + ++j; + } + if (to_push == last_pushed) continue; + if (to_push > maximum_sum) { + j = size; // so we don't keep pushing below. + break; + } + last_pushed = to_push; + new_sums_.push_back(to_push); + } + + // We are sure of this since we will break only when to_push > maximum_sum + // and we are guarantee to have pushed all sums below "to_push" before, that + // includes all the initial sums in sums_. + DCHECK_EQ(i, size); + + for (; j < size; ++j) { + const int64_t to_push = data[j] + e; + if (to_push == last_pushed) continue; + if (to_push > maximum_sum) break; + last_pushed = to_push; + new_sums_.push_back(to_push); + } + std::swap(sums_, new_sums_); + } + return sums_; +} + +double MaxBoundedSubsetSumExact::ComplexityEstimate(int num_elements, + int64_t bin_size) { + double estimate = std::numeric_limits::infinity(); + if (num_elements < 100) { + estimate = 2 * pow(2.0, num_elements / 2); + } + return std::min(estimate, static_cast(bin_size) * + static_cast(num_elements)); +} + +int64_t MaxBoundedSubsetSumExact::MaxSubsetSum( + absl::Span elements, int64_t bin_size) { + // Get rid of a few easy to decide corner cases. + if (elements.empty()) return 0; + if (elements.size() == 1) { + if (elements[0] > bin_size) return 0; + return elements[0]; + } + + // We split the elements in two equal-sized parts. + const int middle = elements.size() / 2; + const auto span_a = sums_a_.Compute(elements.subspan(0, middle), bin_size, + /*abort_if_maximum_reached=*/true); + if (span_a.back() == bin_size) return bin_size; + + const auto span_b = sums_b_.Compute(elements.subspan(middle), bin_size, + /*abort_if_maximum_reached=*/true); + if (span_b.back() == bin_size) return bin_size; + + // For all possible sum a, we compute the largest sum b that fits. + // We do that in linear time thanks to the sorted partial sums. + int64_t result = 0; + CHECK(!span_a.empty()); + CHECK(!span_b.empty()); + int j = span_b.size() - 1; + for (int i = 0; i < span_a.size(); ++i) { + while (j >= 0 && span_a[i] + span_b[j] > bin_size) --j; + result = std::max(result, span_a[i] + span_b[j]); + if (result == bin_size) return bin_size; + } + return result; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/util.h b/ortools/sat/util.h index b455fc1d7e..7f5b1b83cf 100644 --- a/ortools/sat/util.h +++ b/ortools/sat/util.h @@ -502,6 +502,54 @@ class FirstFewValues { std::unique_ptr new_reachable_; }; +// Yet another variant of FirstFewValues or MaxBoundedSubsetSum. +class SortedSubsetSums { + public: + // Computes all the possible subset sums in [0, maximum_sum]. + // Returns them sorted. All elements must be non-negative. + // + // If abort_if_maximum_reached is true, we might not return all possible + // subset sums as we stop the exploration as soon as a subset sum is equal to + // maximum_sum. When this happen, we guarantee that the last element returned + // will be maximum_sum though. + // + // Worst case complexity is in O(2^num_elements) if maximum_sum is large or + // O(maximum_sum * num_elements) if that is lower. + // + // TODO(user): We could optimize even further the case of a small maximum_sum. + absl::Span Compute(absl::Span elements, + int64_t maximum_sum, + bool abort_if_maximum_reached = false); + + // Returns the possible subset sums sorted. + absl::Span SortedSums() const { return sums_; } + + private: + std::vector sums_; + std::vector new_sums_; +}; + +// Similar to MaxBoundedSubsetSum() above but use a different algo. +class MaxBoundedSubsetSumExact { + public: + // If we pack the given elements into a bin of size 'bin_size', returns + // largest possible sum that can be reached. + // + // This implementation allow to solve this in O(2^(num_elements/2)) allowing + // to go easily to 30 or 40 elements. If bin_size is small, complexity is more + // like O(num_element * bin_size). + int64_t MaxSubsetSum(absl::Span elements, int64_t bin_size); + + // Returns an estimate of how many elementary operations + // MaxSubsetSum() is going to take. + double ComplexityEstimate(int num_elements, int64_t bin_size); + + private: + // We use a class just to reuse the memory and not allocate it on each query. + SortedSubsetSums sums_a_; + SortedSubsetSums sums_b_; +}; + // Use Dynamic programming to solve a single knapsack. This is used by the // presolver to simplify variables appearing in a single linear constraint. // diff --git a/ortools/sat/util_test.cc b/ortools/sat/util_test.cc index 91f5b4bc0f..fe4aabd91b 100644 --- a/ortools/sat/util_test.cc +++ b/ortools/sat/util_test.cc @@ -36,6 +36,7 @@ #include "ortools/base/gmock.h" #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" +#include "ortools/base/stl_util.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_solver.h" #include "ortools/sat/cp_model_utils.h" @@ -972,6 +973,82 @@ TEST(WeightedPickTest, SimpleTest) { } } +TEST(SortedSubsetSums, RandomTest) { + const int maximum = 100; + const int num_elements = 20; + absl::BitGen random; + std::vector elements(num_elements); + for (int i = 0; i < num_elements; ++i) { + elements[i] = absl::Uniform(random, 0, maximum); + } + + std::vector brute_force_result; + for (int mask = 0; mask < (1 << num_elements); ++mask) { + int64_t sum = 0; + for (int i = 0; i < num_elements; ++i) { + if ((mask >> i) & 1) sum += elements[i]; + } + if (sum <= maximum) { + brute_force_result.push_back(sum); + } + } + gtl::STLSortAndRemoveDuplicates(&brute_force_result); + + SortedSubsetSums tested; + tested.Compute(elements, maximum); + EXPECT_EQ(tested.SortedSums(), absl::MakeSpan(brute_force_result)); +} + +TEST(SortedSubsetSums, CornerCase) { + SortedSubsetSums helper; + EXPECT_THAT(helper.Compute({0, 5}, 4), ElementsAre(0)); +} + +TEST(MaxBoundedSubsetSumExactTest, RandomTest) { + const int bin_size = 100; + const int num_elements = 20; + absl::BitGen random; + std::vector elements(num_elements); + const int average_size = 2 * bin_size / num_elements; + int64_t sum_of_all = 0; + for (int i = 0; i < num_elements; ++i) { + elements[i] = absl::Uniform(random, average_size - 3, average_size + 3); + sum_of_all += elements[i]; + } + + // Lets compute the maximum by brute force. + int64_t brute_force_result = 0; + for (int mask = 0; mask < (1 << num_elements); ++mask) { + int64_t sum = 0; + for (int i = 0; i < num_elements; ++i) { + if ((mask >> i) & 1) sum += elements[i]; + } + if (sum > bin_size) continue; + brute_force_result = std::max(brute_force_result, sum); + } + + MaxBoundedSubsetSumExact helper; + EXPECT_EQ(brute_force_result, helper.MaxSubsetSum(elements, bin_size)); +} + +TEST(MaxBoundedSubsetSumExactTest, UsedToFail) { + const int64_t capacity = 109; + const std::vector elements = { + 9, 2, 1, 9, 7, 10, 3, 4, 5, 5, 5, 8, 1, 3, 1, 7, 9, 5, 2, 7, 9, 6, 1, + 5, 2, 2, 2, 3, 1, 2, 3, 5, 9, 8, 4, 1, 5, 3, 5, 4, 6, 7, 9, 1, 1, 3}; + + int64_t sum_of_elements = 0; + for (const int64_t e : elements) sum_of_elements += e; + + MaxBoundedSubsetSumExact helper; + EXPECT_EQ(helper.MaxSubsetSum(elements, capacity), capacity); +} + +TEST(MaxBoundedSubsetSumExactTest, CornerCase) { + MaxBoundedSubsetSumExact helper; + EXPECT_EQ(helper.MaxSubsetSum({0, 5, 6}, 4), 0); +} + } // namespace } // namespace sat } // namespace operations_research