From d36c0d7ea51130550bba7b5c981d15583ba91453 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Wed, 19 Jul 2023 11:54:17 -0700 Subject: [PATCH] [CP-SAT] fix 1 overflow; minor speed up of the sat part; add no_overlap_2d dedicated lns --- ortools/sat/clause.cc | 7 ++-- ortools/sat/cp_model_lns.cc | 62 +++++++++++++++++++++++++++++++++ ortools/sat/cp_model_lns.h | 22 ++++++++++++ ortools/sat/cp_model_solver.cc | 58 ++++++++++++++++++++++++++++-- ortools/sat/presolve_context.cc | 20 +++++------ ortools/sat/sat_base.h | 10 ++++-- 6 files changed, 161 insertions(+), 18 deletions(-) diff --git a/ortools/sat/clause.cc b/ortools/sat/clause.cc index f68931cbc0a..d2021e535e3 100644 --- a/ortools/sat/clause.cc +++ b/ortools/sat/clause.cc @@ -721,7 +721,7 @@ bool BinaryImplicationGraph::PropagateOnTrue(Literal true_literal, // loop does slow down the code by a few percent. num_inspections_ += implications_[true_literal.Index()].size(); - for (Literal literal : implications_[true_literal.Index()]) { + for (const Literal literal : implications_[true_literal.Index()]) { if (assignment.LiteralIsTrue(literal)) { // Note(user): I tried to update the reason here if the literal was // enqueued after the true_literal on the trail. This property is @@ -739,7 +739,7 @@ bool BinaryImplicationGraph::PropagateOnTrue(Literal true_literal, } else { // Propagation. reasons_[trail->Index()] = true_literal.Negated(); - trail->Enqueue(literal, propagator_id_); + trail->FastEnqueue(literal); } } @@ -770,7 +770,7 @@ bool BinaryImplicationGraph::PropagateOnTrue(Literal true_literal, } else { // Propagation. reasons_[trail->Index()] = true_literal.Negated(); - trail->Enqueue(literal.Negated(), propagator_id_); + trail->FastEnqueue(literal.Negated()); } } } @@ -784,6 +784,7 @@ bool BinaryImplicationGraph::Propagate(Trail* trail) { propagation_trail_index_ = trail->Index(); return true; } + trail->SetCurrentPropagatorId(propagator_id_); while (propagation_trail_index_ < trail->Index()) { const Literal literal = (*trail)[propagation_trail_index_++]; if (!PropagateOnTrue(literal, trail)) return false; diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index c5b50c45578..ff3786bdc31 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -443,6 +443,32 @@ std::vector NeighborhoodGeneratorHelper::GetActiveIntervals( return active_intervals; } +std::vector> +NeighborhoodGeneratorHelper::GetActiveRectangles( + const CpSolverResponse& initial_solution) const { + const std::vector active_intervals = + GetActiveIntervals(initial_solution); + const absl::flat_hash_set active_intervals_set(active_intervals.begin(), + active_intervals.end()); + + std::vector> active_rectangles; + for (const int ct_index : TypeToConstraints(ConstraintProto::kNoOverlap2D)) { + const NoOverlap2DConstraintProto& ct = + model_proto_.constraints(ct_index).no_overlap_2d(); + for (int i = 0; i < ct.x_intervals_size(); ++i) { + const int x_i = ct.x_intervals(i); + const int y_i = ct.y_intervals(i); + if (!active_intervals_set.contains(x_i) && + !active_intervals_set.contains(y_i)) { + continue; + } + active_rectangles.push_back({x_i, y_i}); + } + } + + return active_rectangles; +} + std::vector> NeighborhoodGeneratorHelper::GetUniqueIntervalSets() const { std::vector> intervals_in_constraints; @@ -2012,6 +2038,42 @@ Neighborhood SchedulingResourceWindowsNeighborhoodGenerator::Generate( intervals, initial_solution, random, helper_); } +Neighborhood RandomRectanglesPackingNeighborhoodGenerator::Generate( + const CpSolverResponse& initial_solution, double difficulty, + absl::BitGenRef random) { + std::vector> rectangles_to_freeze = + helper_.GetActiveRectangles(initial_solution); + GetRandomSubset(1.0 - difficulty, &rectangles_to_freeze, random); + + // ConstraintToVar() is not initialised for intervals. We need to parse them + // manually. + auto insert_vars_from_intervals = + [this](int i, absl::flat_hash_set& vars_to_freeze) { + const ConstraintProto& ct = helper_.ModelProto().constraints(i); + for (const int lit : ct.enforcement_literal()) { + const int var = PositiveRef(lit); + vars_to_freeze.insert(var); + } + for (const int var : ct.interval().start().vars()) { + vars_to_freeze.insert(var); + } + for (const int var : ct.interval().end().vars()) { + vars_to_freeze.insert(var); + } + for (const int var : ct.interval().size().vars()) { + vars_to_freeze.insert(var); + } + }; + + absl::flat_hash_set variables_to_freeze; + for (const auto& [x, y] : rectangles_to_freeze) { + insert_vars_from_intervals(x, variables_to_freeze); + insert_vars_from_intervals(y, variables_to_freeze); + } + + return helper_.FixGivenVariables(initial_solution, variables_to_freeze); +} + Neighborhood RoutingRandomNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, absl::BitGenRef random) { diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index f3b4b11c80b..074c9cac502 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -208,6 +208,14 @@ class NeighborhoodGeneratorHelper : public SubSolver { std::vector GetActiveIntervals( const CpSolverResponse& initial_solution) const; + // Returns the list of active rectangles appearing in no_overlap_2d + // constraints according to the initial_solution and the parameter + // lns_focus_on_performed_intervals. If true, this method returns the list of + // performed rectangles in the solution. If false, it returns all rectangles + // of the model. + std::vector> GetActiveRectangles( + const CpSolverResponse& initial_solution) const; + // Returns the set of unique intervals list appearing in a no_overlap, // cumulative, or as a dimension of a no_overlap_2d constraint. std::vector> GetUniqueIntervalSets() const; @@ -672,6 +680,20 @@ class SchedulingResourceWindowsNeighborhoodGenerator absl::flat_hash_set intervals_to_relax_; }; +// Only make sense for problems with no_overlap_2d constraints. This select a +// random set of rectangles (i.e. a pair of intervals) of the problem according +// to the difficulty. Then, fix all variables in the selected intervals. +class RandomRectanglesPackingNeighborhoodGenerator + : public NeighborhoodGenerator { + public: + explicit RandomRectanglesPackingNeighborhoodGenerator( + NeighborhoodGeneratorHelper const* helper, const std::string& name) + : NeighborhoodGenerator(name, helper) {} + + Neighborhood Generate(const CpSolverResponse& initial_solution, + double difficulty, absl::BitGenRef random) final; +}; + // This routing based LNS generator will relax random arcs in all the paths of // the circuit or routes constraints. class RoutingRandomNeighborhoodGenerator : public NeighborhoodGenerator { diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index b4b73e9d6d3..f0d64dc2b2d 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -223,15 +223,20 @@ std::string CpModelStats(const CpModelProto& model_proto) { int no_overlap_2d_num_optional_rectangles = 0; int no_overlap_2d_num_linear_areas = 0; int no_overlap_2d_num_quadratic_areas = 0; + int no_overlap_2d_num_fixed_rectangles = 0; int cumulative_num_intervals = 0; int cumulative_num_optional_intervals = 0; int cumulative_num_variable_sizes = 0; int cumulative_num_variable_demands = 0; + int cumulative_num_fixed_intervals = 0; int no_overlap_num_intervals = 0; int no_overlap_num_optional_intervals = 0; int no_overlap_num_variable_sizes = 0; + int no_overlap_num_fixed_intervals = 0; + + int num_fixed_intervals = 0; for (const ConstraintProto& ct : model_proto.constraints()) { std::string name = ConstraintCaseName(ct.constraint_case()); @@ -278,6 +283,19 @@ std::string CpModelStats(const CpModelProto& model_proto) { return !model_proto.constraints(i).enforcement_literal().empty(); }; + auto interval_is_fixed = [&variable_is_fixed, + expression_is_fixed](const ConstraintProto& ct) { + if (ct.constraint_case() != ConstraintProto::ConstraintCase::kInterval) { + return false; + } + for (const int lit : ct.enforcement_literal()) { + if (!variable_is_fixed(lit)) return false; + } + return (expression_is_fixed(ct.interval().start()) && + expression_is_fixed(ct.interval().size()) && + expression_is_fixed(ct.interval().end())); + }; + // For pure Boolean constraints, we also display the total number of literal // involved as this gives a good idea of the problem size. if (ct.constraint_case() == ConstraintProto::ConstraintCase::kBoolOr) { @@ -295,6 +313,9 @@ std::string CpModelStats(const CpModelProto& model_proto) { } else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinMax) { name_to_num_expressions[name] += ct.lin_max().exprs().size(); + } else if (ct.constraint_case() == + ConstraintProto::ConstraintCase::kInterval) { + if (interval_is_fixed(ct)) num_fixed_intervals++; } else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kNoOverlap2D) { const int num_boxes = ct.no_overlap_2d().x_intervals_size(); @@ -313,6 +334,10 @@ std::string CpModelStats(const CpModelProto& model_proto) { } else if (num_fixed == 1) { no_overlap_2d_num_linear_areas++; } + if (interval_is_fixed(model_proto.constraints(x_interval)) && + interval_is_fixed(model_proto.constraints(y_interval))) { + no_overlap_2d_num_fixed_rectangles++; + } } } else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kNoOverlap) { @@ -326,6 +351,9 @@ std::string CpModelStats(const CpModelProto& model_proto) { if (!interval_has_fixed_size(interval)) { no_overlap_num_variable_sizes++; } + if (interval_is_fixed(model_proto.constraints(interval))) { + no_overlap_num_fixed_intervals++; + } } } else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kCumulative) { @@ -342,6 +370,9 @@ std::string CpModelStats(const CpModelProto& model_proto) { if (!expression_is_fixed(ct.cumulative().demands(i))) { cumulative_num_variable_demands++; } + if (interval_is_fixed(model_proto.constraints(interval))) { + cumulative_num_fixed_intervals++; + } } } @@ -513,7 +544,10 @@ std::string CpModelStats(const CpModelProto& model_proto) { " (#complex_domain: ", name_to_num_complex_domain[name], ")"); } - if (name == "kNoOverlap2D") { + if (name == "kInterval" && num_fixed_intervals > 0) { + absl::StrAppend(&constraints.back(), " (#fixed: ", num_fixed_intervals, + ")"); + } else if (name == "kNoOverlap2D") { absl::StrAppend(&constraints.back(), " (#rectangles: ", no_overlap_2d_num_rectangles); if (no_overlap_2d_num_optional_rectangles > 0) { @@ -528,6 +562,10 @@ std::string CpModelStats(const CpModelProto& model_proto) { absl::StrAppend(&constraints.back(), ", #quadratic_areas: ", no_overlap_2d_num_quadratic_areas); } + if (no_overlap_2d_num_fixed_rectangles > 0) { + absl::StrAppend(&constraints.back(), ", #fixed_rectangles: ", + no_overlap_2d_num_fixed_rectangles); + } absl::StrAppend(&constraints.back(), ")"); } else if (name == "kCumulative") { absl::StrAppend(&constraints.back(), @@ -544,6 +582,10 @@ std::string CpModelStats(const CpModelProto& model_proto) { absl::StrAppend(&constraints.back(), ", #variable_demands: ", cumulative_num_variable_demands); } + if (cumulative_num_fixed_intervals > 0) { + absl::StrAppend(&constraints.back(), + ", #fixed_intervals: ", cumulative_num_fixed_intervals); + } absl::StrAppend(&constraints.back(), ")"); } else if (name == "kNoOverlap") { absl::StrAppend(&constraints.back(), @@ -556,6 +598,10 @@ std::string CpModelStats(const CpModelProto& model_proto) { absl::StrAppend(&constraints.back(), ", #variable_sizes: ", no_overlap_num_variable_sizes); } + if (no_overlap_num_fixed_intervals > 0) { + absl::StrAppend(&constraints.back(), + ", #fixed_intervals: ", no_overlap_num_fixed_intervals); + } absl::StrAppend(&constraints.back(), ")"); } } @@ -3096,7 +3142,7 @@ class LnsSolver : public SubSolver { // TODO(user): export the delta too if needed. const std::string lns_name = absl::StrCat(absl::GetFlag(FLAGS_cp_model_dump_prefix), - lns_fragment.name(), ".pb.txt"); + lns_fragment.name(), "_", source_info, ".pb.txt"); LOG(INFO) << "Dumping LNS model to '" << lns_name << "'."; CHECK(WriteModelProtoToFile(lns_fragment, lns_name)); } @@ -3173,7 +3219,7 @@ class LnsSolver : public SubSolver { if (absl::GetFlag(FLAGS_cp_model_dump_problematic_lns)) { const std::string name = absl::StrCat(absl::GetFlag(FLAGS_cp_model_dump_prefix), - debug_copy.name(), ".pb.txt"); + debug_copy.name(), "_", source_info, ".pb.txt"); LOG(INFO) << "Dumping problematic LNS model to '" << name << "'."; CHECK(WriteModelProtoToFile(debug_copy, name)); } @@ -3639,6 +3685,12 @@ void SolveCpModelParallel(const CpModelProto& model_proto, "scheduling_resource_windows_lns"), params, helper, &shared)); } + if (!helper->TypeToConstraints(ConstraintProto::kNoOverlap2D).empty()) { + subsolvers.push_back(std::make_unique( + std::make_unique( + helper, "packing_rectangles_lns"), + params, helper, &shared)); + } } const int num_circuit = static_cast( diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index e152cff54be..9595f2116c1 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -1693,17 +1693,17 @@ bool PresolveContext::CanonicalizeObjective(bool simplify_domain) { objective_scaling_factor_ *= static_cast(gcd); // We update the offset accordingly. - const absl::int128 offset = - absl::int128(objective_integer_before_offset_) * - absl::int128(objective_integer_scaling_factor_) + - absl::int128(objective_integer_after_offset_); - - if (objective_domain_.IsFixed() && objective_domain_.FixedValue() * gcd == - -objective_integer_before_offset_) { - // We avoid a corner case where this would overflow but the objective is - // zero. In this case any factor work, so we just take 1 and avoid the - // overflow. + absl::int128 offset = absl::int128(objective_integer_before_offset_) * + absl::int128(objective_integer_scaling_factor_) + + absl::int128(objective_integer_after_offset_); + + if (objective_domain_.IsFixed()) { + // To avoid overflow in (fixed_value * gcd + before_offset) * factor + + // after_offset because the objective is constant (and should fit on an + // int64_t), we can rewrite it as fixed_value + offset. objective_integer_scaling_factor_ = 1; + offset += + absl::int128(gcd - 1) * absl::int128(objective_domain_.FixedValue()); } else { objective_integer_scaling_factor_ *= gcd; } diff --git a/ortools/sat/sat_base.h b/ortools/sat/sat_base.h index e28028025c1..e0113a385a6 100644 --- a/ortools/sat/sat_base.h +++ b/ortools/sat/sat_base.h @@ -260,14 +260,20 @@ class Trail { // Enqueues the assignment that make the given literal true on the trail. This // should only be called on unassigned variables. - void Enqueue(Literal true_literal, int propagator_id) { + void SetCurrentPropagatorId(int propagator_id) { + current_info_.type = propagator_id; + } + void FastEnqueue(Literal true_literal) { DCHECK(!assignment_.VariableIsAssigned(true_literal.Variable())); trail_[current_info_.trail_index] = true_literal; - current_info_.type = propagator_id; info_[true_literal.Variable()] = current_info_; assignment_.AssignFromTrueLiteral(true_literal); ++current_info_.trail_index; } + void Enqueue(Literal true_literal, int propagator_id) { + SetCurrentPropagatorId(propagator_id); + FastEnqueue(true_literal); + } // Specific Enqueue() version for the search decision. void EnqueueSearchDecision(Literal true_literal) {