From 198a4d3f120bf1a5676e90e5e1ebd7fae0be43bb Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Mon, 16 Dec 2024 14:20:09 +0100 Subject: [PATCH] [CP-SAT] fir rare bug in dependency graph and implications; speed up no_overlap_2d and cumulative relaxation --- ortools/sat/BUILD.bazel | 3 +- ortools/sat/circuit.cc | 7 +- ortools/sat/circuit.h | 6 +- ortools/sat/clause.cc | 10 ++ ortools/sat/clause_test.cc | 10 ++ ortools/sat/cp_model_lns.cc | 152 +++++++++---------- ortools/sat/cp_model_lns.h | 7 +- ortools/sat/cp_model_postsolve.cc | 2 +- ortools/sat/cp_model_presolve.cc | 31 +++- ortools/sat/cp_model_presolve_random_test.cc | 1 + ortools/sat/cp_model_search.cc | 2 +- ortools/sat/cp_model_search.h | 2 +- ortools/sat/cp_model_solver.cc | 29 ++-- ortools/sat/cp_model_solver_fuzz.cc | 5 +- ortools/sat/cp_model_utils.h | 21 ++- ortools/sat/diffn.cc | 62 +++----- ortools/sat/diffn_util.cc | 3 + ortools/sat/integer_expr.cc | 38 ++--- ortools/sat/integer_expr.h | 10 +- ortools/sat/integer_expr_test.cc | 26 ++++ ortools/sat/lp_utils.cc | 2 +- ortools/sat/presolve_context.h | 4 + ortools/sat/routing_cuts.cc | 2 +- ortools/sat/routing_cuts.h | 2 +- ortools/sat/var_domination_test.cc | 6 +- 25 files changed, 262 insertions(+), 181 deletions(-) diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index c7d1383407..4bf59d5d7d 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -125,6 +125,7 @@ cc_library( "//ortools/base:file", "//ortools/base:hash", "//ortools/base:stl_util", + "//ortools/util:bitset", "//ortools/util:saturated_arithmetic", "//ortools/util:sorted_interval_list", "@com_google_absl//absl/container:flat_hash_map", @@ -3109,7 +3110,7 @@ cc_test( ":integer_base", ":util", "//ortools/base", - "//ortools/base:gmock", + "//ortools/base:gmock_main", "//ortools/graph:connected_components", "//ortools/graph:strongly_connected_components", "//ortools/util:saturated_arithmetic", diff --git a/ortools/sat/circuit.cc b/ortools/sat/circuit.cc index 5208f52c39..92ba0a8f4c 100644 --- a/ortools/sat/circuit.cc +++ b/ortools/sat/circuit.cc @@ -358,10 +358,9 @@ bool CircuitPropagator::Propagate() { return true; } -NoCyclePropagator::NoCyclePropagator(int num_nodes, - const std::vector& tails, - const std::vector& heads, - const std::vector& literals, +NoCyclePropagator::NoCyclePropagator(int num_nodes, absl::Span tails, + absl::Span heads, + absl::Span literals, Model* model) : num_nodes_(num_nodes), trail_(model->GetOrCreate()), diff --git a/ortools/sat/circuit.h b/ortools/sat/circuit.h index 74f49c9961..3c74e70a67 100644 --- a/ortools/sat/circuit.h +++ b/ortools/sat/circuit.h @@ -119,9 +119,9 @@ class CircuitPropagator : PropagatorInterface, ReversibleInterface { // Enforce the fact that there is no cycle in the given directed graph. class NoCyclePropagator : PropagatorInterface, ReversibleInterface { public: - NoCyclePropagator(int num_nodes, const std::vector& tails, - const std::vector& heads, - const std::vector& literals, Model* model); + NoCyclePropagator(int num_nodes, absl::Span tails, + absl::Span heads, + absl::Span literals, Model* model); void SetLevel(int level) final; bool Propagate() final; diff --git a/ortools/sat/clause.cc b/ortools/sat/clause.cc index 305821ff6a..dc1d174095 100644 --- a/ortools/sat/clause.cc +++ b/ortools/sat/clause.cc @@ -1446,7 +1446,17 @@ bool BinaryImplicationGraph::DetectEquivalences(bool log_info) { // that this might result in more implications when we expand small at most // one. at_most_ones_.clear(); + int saved_trail_index = propagation_trail_index_; CleanUpAndAddAtMostOnes(/*base_index=*/0); + // This might have run the propagation on a few variables without taking + // into account the AMOs. Propagate again. + // + // TODO(user): Maybe a better alternative is to not propagate when we fix + // variables inside CleanUpAndAddAtMostOnes(). + if (propagation_trail_index_ != saved_trail_index) { + propagation_trail_index_ = saved_trail_index; + Propagate(trail_); + } num_implications_ = 0; for (LiteralIndex i(0); i < size; ++i) { diff --git a/ortools/sat/clause_test.cc b/ortools/sat/clause_test.cc index c950e6bb92..b7d1ff0895 100644 --- a/ortools/sat/clause_test.cc +++ b/ortools/sat/clause_test.cc @@ -108,6 +108,16 @@ TEST(BinaryImplicationGraphTest, BasicUnsatSccTest) { EXPECT_FALSE(graph->DetectEquivalences()); } +TEST(BinaryImplicationGraphTest, IssueFoundOnMipLibTest) { + Model model; + auto* graph = model.GetOrCreate(); + model.GetOrCreate()->SetNumVariables(10); + EXPECT_TRUE(graph->AddAtMostOne(Literals({+1, +2, -3, -4}))); + EXPECT_TRUE(graph->AddAtMostOne(Literals({+1, +2, +3, +4, +5}))); + EXPECT_TRUE(graph->AddAtMostOne(Literals({+1, -2, -3, +4, +5}))); + EXPECT_TRUE(graph->DetectEquivalences()); +} + TEST(BinaryImplicationGraphTest, DetectEquivalences) { // We take a bunch of random permutations, equivalence classes will be cycles. // We make sure the representative of x and not(x) are always negation of diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 9c22c5950f..3482ab2f8d 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -558,7 +558,7 @@ struct TimePartition { // Selects all intervals in a random time window to meet the difficulty // requirement. TimePartition PartitionIndicesAroundRandomTimeWindow( - const std::vector& intervals, const CpModelProto& model_proto, + absl::Span intervals, const CpModelProto& model_proto, const CpSolverResponse& initial_solution, double difficulty, absl::BitGenRef random) { std::vector start_end_indices; @@ -1036,15 +1036,10 @@ std::vector> NeighborhoodGeneratorHelper::GetRoutingPaths( Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( const CpSolverResponse& base_solution, - const absl::flat_hash_set& variables_to_fix) const { - int initial_num_variables = 0; - { - absl::ReaderMutexLock domain_lock(&domain_mutex_); - - initial_num_variables = - model_proto_with_only_variables_->variables().size(); - } - Neighborhood neighborhood(initial_num_variables); + const Bitset64& variables_to_fix) const { + const int num_variables = variables_to_fix.size(); + Neighborhood neighborhood(num_variables); + neighborhood.delta.mutable_variables()->Reserve(num_variables); // TODO(user): Maybe relax all variables in the objective when the number // is small or negligible compared to the number of variables. @@ -1054,12 +1049,9 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( : -1; // Fill in neighborhood.delta all variable domains. + int num_fixed = 0; { absl::ReaderMutexLock domain_lock(&domain_mutex_); - - const int num_variables = - model_proto_with_only_variables_->variables().size(); - neighborhood.delta.mutable_variables()->Reserve(num_variables); for (int i = 0; i < num_variables; ++i) { const IntegerVariableProto& current_var = model_proto_with_only_variables_->variables(i); @@ -1068,17 +1060,20 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( // We only copy the name in debug mode. if (DEBUG_MODE) new_var->set_name(current_var.name()); - const Domain domain = ReadDomainFromProto(current_var); - const int64_t base_value = base_solution.solution(i); + if (variables_to_fix[i] && i != unique_objective_variable) { + ++num_fixed; - if (variables_to_fix.contains(i) && i != unique_objective_variable) { - if (domain.Contains(base_value)) { + // Note the use of DomainInProtoContains() instead of + // ReadDomainFromProto() as the later is slower and allocate memory. + const int64_t base_value = base_solution.solution(i); + if (DomainInProtoContains(current_var, base_value)) { new_var->add_domain(base_value); new_var->add_domain(base_value); } else { // If under the updated domain, the base solution is no longer valid, // We should probably regenerate this neighborhood. But for now we // just do a best effort and take the closest value. + const Domain domain = ReadDomainFromProto(current_var); int64_t closest_value = domain.Min(); int64_t closest_dist = std::abs(closest_value - base_value); for (const ClosedInterval interval : domain) { @@ -1093,7 +1088,7 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( FillDomainInProto(Domain(closest_value, closest_value), new_var); } } else { - FillDomainInProto(domain, new_var); + *new_var->mutable_domain() = current_var.domain(); } } } @@ -1138,10 +1133,15 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( neighborhood.variables_that_can_be_fixed_to_local_optimum.clear(); } + const int num_relaxed = num_variables - num_fixed; + neighborhood.delta.mutable_solution_hint()->mutable_vars()->Reserve( + num_relaxed); + neighborhood.delta.mutable_solution_hint()->mutable_values()->Reserve( + num_relaxed); AddSolutionHinting(base_solution, &neighborhood.delta); neighborhood.is_generated = true; - neighborhood.is_reduced = !variables_to_fix.empty(); + neighborhood.is_reduced = num_fixed > 0; neighborhood.is_simple = true; // TODO(user): force better objective? Note that this is already done when the @@ -1172,25 +1172,26 @@ void NeighborhoodGeneratorHelper::AddSolutionHinting( Neighborhood NeighborhoodGeneratorHelper::RelaxGivenVariables( const CpSolverResponse& initial_solution, absl::Span relaxed_variables) const { - std::vector relaxed_variables_set(model_proto_.variables_size(), false); - for (const int var : relaxed_variables) relaxed_variables_set[var] = true; - absl::flat_hash_set fixed_variables; + Bitset64 fixed_variables(NumVariables()); { absl::ReaderMutexLock graph_lock(&graph_mutex_); for (const int i : active_variables_) { - if (!relaxed_variables_set[i]) { - fixed_variables.insert(i); - } + fixed_variables.Set(i); } } + for (const int var : relaxed_variables) fixed_variables.Clear(var); return FixGivenVariables(initial_solution, fixed_variables); } Neighborhood NeighborhoodGeneratorHelper::FixAllVariables( const CpSolverResponse& initial_solution) const { - const std::vector& all_variables = ActiveVariables(); - const absl::flat_hash_set fixed_variables(all_variables.begin(), - all_variables.end()); + Bitset64 fixed_variables(NumVariables()); + { + absl::ReaderMutexLock graph_lock(&graph_mutex_); + for (const int i : active_variables_) { + fixed_variables.Set(i); + } + } return FixGivenVariables(initial_solution, fixed_variables); } @@ -1319,8 +1320,10 @@ Neighborhood RelaxRandomVariablesGenerator::Generate( absl::BitGenRef random) { std::vector fixed_variables = helper_.ActiveVariables(); GetRandomSubset(1.0 - data.difficulty, &fixed_variables, random); - return helper_.FixGivenVariables( - initial_solution, {fixed_variables.begin(), fixed_variables.end()}); + + Bitset64 to_fix(helper_.NumVariables()); + for (const int var : fixed_variables) to_fix.Set(var); + return helper_.FixGivenVariables(initial_solution, to_fix); } Neighborhood RelaxRandomConstraintsGenerator::Generate( @@ -2302,14 +2305,13 @@ Neighborhood RandomRectanglesPackingNeighborhoodGenerator::Generate( helper_.GetActiveRectangles(initial_solution); GetRandomSubset(1.0 - data.difficulty, &rectangles_to_freeze, random); - absl::flat_hash_set variables_to_freeze; + Bitset64 variables_to_freeze(helper_.NumVariables()); for (const ActiveRectangle& rectangle : rectangles_to_freeze) { - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval, - variables_to_freeze); - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval, - variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), rectangle.x_interval, + variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), rectangle.y_interval, + variables_to_freeze); } - return helper_.FixGivenVariables(initial_solution, variables_to_freeze); } @@ -2351,15 +2353,15 @@ Neighborhood RectanglesPackingRelaxOneNeighborhoodGenerator::Generate( // // Note that we only consider two rectangles as potential neighbors if they // are part of the same no_overlap_2d constraint. - absl::flat_hash_set variables_to_freeze; + Bitset64 variables_to_freeze(helper_.NumVariables()); std::vector> distances; distances.reserve(all_active_rectangles.size()); for (int i = 0; i < all_active_rectangles.size(); ++i) { const ActiveRectangle& rectangle = all_active_rectangles[i]; - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval, - variables_to_freeze); - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval, - variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), rectangle.x_interval, + variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), rectangle.y_interval, + variables_to_freeze); const Rectangle rect = get_rectangle(rectangle); const bool same_no_overlap_as_center_rect = absl::c_any_of( @@ -2398,14 +2400,10 @@ Neighborhood RectanglesPackingRelaxOneNeighborhoodGenerator::Generate( for (const int b : boxes_to_relax) { const ActiveRectangle& rectangle = all_active_rectangles[b]; - absl::flat_hash_set variables_to_relax; - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval, - variables_to_relax); - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval, - variables_to_relax); - for (const int v : variables_to_relax) { - variables_to_freeze.erase(v); - } + RemoveVariablesFromInterval(helper_.ModelProto(), rectangle.x_interval, + variables_to_freeze); + RemoveVariablesFromInterval(helper_.ModelProto(), rectangle.y_interval, + variables_to_freeze); } Neighborhood neighborhood = helper_.FixGivenVariables(initial_solution, variables_to_freeze); @@ -2489,17 +2487,17 @@ Neighborhood RectanglesPackingRelaxTwoNeighborhoodsGenerator::Generate( // TODO(user): This computes the distance between the center of the // rectangles. We could use the real distance between the closest points, but // not sure it is worth the extra complexity. - absl::flat_hash_set variables_to_freeze; + Bitset64 variables_to_freeze(helper_.NumVariables()); std::vector> distances1; std::vector> distances2; distances1.reserve(all_active_rectangles.size()); distances2.reserve(all_active_rectangles.size()); for (int i = 0; i < all_active_rectangles.size(); ++i) { const ActiveRectangle& rectangle = all_active_rectangles[i]; - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval, - variables_to_freeze); - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval, - variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), rectangle.x_interval, + variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), rectangle.y_interval, + variables_to_freeze); const Rectangle rect = get_rectangle(rectangle); const bool same_no_overlap_as_rect1 = @@ -2525,22 +2523,18 @@ Neighborhood RectanglesPackingRelaxTwoNeighborhoodsGenerator::Generate( [](const auto& a, const auto& b) { return a.second < b.second; }); std::sort(distances2.begin(), distances2.end(), [](const auto& a, const auto& b) { return a.second < b.second; }); - absl::flat_hash_set variables_to_relax; for (auto& samples : {distances1, distances2}) { const int num_potential_samples = samples.size(); for (int i = 0; i < std::min(num_potential_samples, num_to_sample_each); ++i) { const int rectangle_idx = samples[i].first; const ActiveRectangle& rectangle = all_active_rectangles[rectangle_idx]; - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.x_interval, - variables_to_relax); - InsertVariablesFromConstraint(helper_.ModelProto(), rectangle.y_interval, - variables_to_relax); + RemoveVariablesFromInterval(helper_.ModelProto(), rectangle.x_interval, + variables_to_freeze); + RemoveVariablesFromInterval(helper_.ModelProto(), rectangle.y_interval, + variables_to_freeze); } } - for (const int v : variables_to_relax) { - variables_to_freeze.erase(v); - } return helper_.FixGivenVariables(initial_solution, variables_to_freeze); } @@ -2583,15 +2577,15 @@ Neighborhood SlicePackingNeighborhoodGenerator::Generate( indices_to_fix[index] = false; } - absl::flat_hash_set variables_to_freeze; + Bitset64 variables_to_freeze(helper_.NumVariables()); for (int index = 0; index < active_rectangles.size(); ++index) { if (indices_to_fix[index]) { - InsertVariablesFromConstraint(helper_.ModelProto(), - active_rectangles[index].x_interval, - variables_to_freeze); - InsertVariablesFromConstraint(helper_.ModelProto(), - active_rectangles[index].y_interval, - variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), + active_rectangles[index].x_interval, + variables_to_freeze); + InsertVariablesFromInterval(helper_.ModelProto(), + active_rectangles[index].y_interval, + variables_to_freeze); } } @@ -2613,8 +2607,10 @@ Neighborhood RoutingRandomNeighborhoodGenerator::Generate( all_path_variables.end()); std::sort(fixed_variables.begin(), fixed_variables.end()); GetRandomSubset(1.0 - data.difficulty, &fixed_variables, random); - return helper_.FixGivenVariables( - initial_solution, {fixed_variables.begin(), fixed_variables.end()}); + + Bitset64 to_fix(helper_.NumVariables()); + for (const int var : fixed_variables) to_fix.Set(var); + return helper_.FixGivenVariables(initial_solution, to_fix); } Neighborhood RoutingPathNeighborhoodGenerator::Generate( @@ -2656,11 +2652,11 @@ Neighborhood RoutingPathNeighborhoodGenerator::Generate( } // Compute the set of variables to fix. - absl::flat_hash_set fixed_variables; + Bitset64 to_fix(helper_.NumVariables()); for (const int var : all_path_variables) { - if (!relaxed_variables.contains(var)) fixed_variables.insert(var); + if (!relaxed_variables.contains(var)) to_fix.Set(var); } - return helper_.FixGivenVariables(initial_solution, fixed_variables); + return helper_.FixGivenVariables(initial_solution, to_fix); } Neighborhood RoutingFullPathNeighborhoodGenerator::Generate( @@ -2722,11 +2718,11 @@ Neighborhood RoutingFullPathNeighborhoodGenerator::Generate( } // Compute the set of variables to fix. - absl::flat_hash_set fixed_variables; + Bitset64 to_fix(helper_.NumVariables()); for (const int var : all_path_variables) { - if (!relaxed_variables.contains(var)) fixed_variables.insert(var); + if (!relaxed_variables.contains(var)) to_fix.Set(var); } - return helper_.FixGivenVariables(initial_solution, fixed_variables); + return helper_.FixGivenVariables(initial_solution, to_fix); } bool RelaxationInducedNeighborhoodGenerator::ReadyToGenerate() const { diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index f9051366f7..dc36e83415 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -115,11 +115,12 @@ class NeighborhoodGeneratorHelper : public SubSolver { } void Synchronize() override; + int NumVariables() const { return model_proto_.variables().size(); } + // Returns the LNS fragment where the given variables are fixed to the value // they take in the given solution. - Neighborhood FixGivenVariables( - const CpSolverResponse& base_solution, - const absl::flat_hash_set& variables_to_fix) const; + Neighborhood FixGivenVariables(const CpSolverResponse& base_solution, + const Bitset64& variables_to_fix) const; // Returns the LNS fragment which will relax all inactive variables and all // variables in relaxed_variables. diff --git a/ortools/sat/cp_model_postsolve.cc b/ortools/sat/cp_model_postsolve.cc index 4a522f554d..662151ec29 100644 --- a/ortools/sat/cp_model_postsolve.cc +++ b/ortools/sat/cp_model_postsolve.cc @@ -215,7 +215,7 @@ namespace { // Note that if a domain is not fixed, we just take its Min() value. int64_t EvaluateLinearExpression(const LinearExpressionProto& expr, - const std::vector& domains) { + absl::Span domains) { int64_t value = expr.offset(); for (int i = 0; i < expr.vars_size(); ++i) { const int ref = expr.vars(i); diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 59a42516f5..7f0f37a505 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -11834,7 +11834,21 @@ void CpModelPresolver::ProcessVariableOnlyUsedInEncoding(int var) { // Nothing else to do in this case, the constraint will be reduced to // a pure Boolean constraint later. context_->UpdateRuleStats("variables: reduced domain to two values"); - return (void)context_->IntersectDomainWithAndUpdateHint( + // If the hint enforces `ct`, then the hint of `var` must be in the + // implied domain. In any case its new value must not increase the + // objective value (the objective domain is non-constraining, but this + // only guarantees that `var` can freely *decrease* the objective). The + // code below ensures this (`value2` is the 'cheapest' value the implied + // domain, and `value1` the cheapest value in the variable's domain). + bool enforcing_hint = true; + for (const int enforcement_lit : ct.enforcement_literal()) { + if (context_->LiteralSolutionHintIs(enforcement_lit, false)) { + enforcing_hint = false; + break; + } + } + context_->UpdateVarSolutionHint(var, enforcing_hint ? value2 : value1); + return (void)context_->IntersectDomainWith( var, Domain::FromValues({value1, value2})); } } @@ -12938,6 +12952,13 @@ bool ModelCopy::CopyLinear(const ConstraintProto& ct) { DCHECK(!non_fixed_variables_.empty()); + if (non_fixed_variables_.size() == 1 && ct.enforcement_literal().empty()) { + context_->UpdateRuleStats("linear1: x in domain"); + return context_->IntersectDomainWith( + non_fixed_variables_[0], + new_rhs.InverseMultiplicationBy(non_fixed_coefficients_[0])); + } + ConstraintProto* new_ct = context_->working_model->add_constraints(); FinishEnforcementCopy(new_ct); LinearConstraintProto* linear = new_ct->mutable_linear(); @@ -14365,7 +14386,15 @@ bool CpModelPresolver::MaybeRemoveFixedVariables( } mutable_vars->Truncate(new_size); mutable_coeffs->Truncate(new_size); + + // Tricky: The objective domain is on the linear objective expression + // without the offset, so we need to shift it if there is one. objective->set_offset(objective->offset() + offset); + if (!objective->domain().empty()) { + Domain obj_domain = + ReadDomainFromProto(*objective).AdditionWith(Domain(-offset)); + FillDomainInProto(obj_domain, objective); + } context_->ReadObjectiveFromProto(); if (!context_->CanonicalizeObjective()) return false; diff --git a/ortools/sat/cp_model_presolve_random_test.cc b/ortools/sat/cp_model_presolve_random_test.cc index dca4a12157..a16503b7bd 100644 --- a/ortools/sat/cp_model_presolve_random_test.cc +++ b/ortools/sat/cp_model_presolve_random_test.cc @@ -239,6 +239,7 @@ TEST_P(RandomPreprocessorTest, TestHintSurvivePresolve) { // We just check that the hint is correct. SatParameters params; params.set_debug_crash_on_bad_hint(true); + params.set_debug_crash_if_presolve_breaks_hint(true); params.set_stop_after_first_solution(true); const CpSolverResponse with_hint = SolveWithParameters(model_proto, params); diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index ba580e810b..ee6ec2e685 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -357,7 +357,7 @@ std::function ConstructHeuristicSearchStrategy( std::function ConstructIntegerCompletionSearchStrategy( - const std::vector& variable_mapping, + absl::Span variable_mapping, IntegerVariable objective_var, Model* model) { const auto& params = *model->GetOrCreate(); if (!params.instantiate_all_variables()) { diff --git a/ortools/sat/cp_model_search.h b/ortools/sat/cp_model_search.h index 3892e97ac9..0518ebcd68 100644 --- a/ortools/sat/cp_model_search.h +++ b/ortools/sat/cp_model_search.h @@ -78,7 +78,7 @@ std::function ConstructHeuristicSearchStrategy( // Constructs an integer completion search strategy. std::function ConstructIntegerCompletionSearchStrategy( - const std::vector& variable_mapping, + absl::Span variable_mapping, IntegerVariable objective_var, Model* model); // Constructs a search strategy that follows the hints from the model. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 4aa90acc7a..60242a452b 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -810,12 +810,12 @@ void RestrictObjectiveUsingHint(CpModelProto* model_proto) { } } -// Returns false if there is a complete solution hint that is infeasible, or -// true otherwise. -bool TestSolutionHintForFeasibility(const CpModelProto& model_proto, - SolverLogger* logger = nullptr, - SharedResponseManager* manager = nullptr) { - if (!model_proto.has_solution_hint()) return true; +// Returns true iff there is a hint, and (ignoring fixed variables) if it is +// complete and feasible. +bool SolutionHintIsCompleteAndFeasible( + const CpModelProto& model_proto, SolverLogger* logger = nullptr, + SharedResponseManager* manager = nullptr) { + if (!model_proto.has_solution_hint()) return false; int num_active_variables = 0; int num_hinted_variables = 0; @@ -837,7 +837,7 @@ bool TestSolutionHintForFeasibility(const CpModelProto& model_proto, logger, "The solution hint is incomplete: ", num_hinted_variables, " out of ", num_active_variables, " non fixed variables hinted."); } - return true; + return false; } std::vector solution(model_proto.variables_size(), 0); @@ -1287,8 +1287,8 @@ class LnsSolver : public SubSolver { absl::MutexLock l(&next_arena_size_mutex_); buffer_size = next_arena_size_; } - std::vector arena_buffer(buffer_size); - google::protobuf::Arena arena(arena_buffer.data(), arena_buffer.size()); + auto arena_buffer = std::make_unique(buffer_size); + google::protobuf::Arena arena(arena_buffer.get(), buffer_size); CpModelProto& lns_fragment = *google::protobuf::Arena::Create(&arena); CpModelProto& mapping_proto = @@ -1359,7 +1359,7 @@ class LnsSolver : public SubSolver { bool hint_feasible_before_presolve = false; if (lns_parameters_.debug_crash_if_presolve_breaks_hint()) { hint_feasible_before_presolve = - TestSolutionHintForFeasibility(lns_fragment, /*logger=*/nullptr); + SolutionHintIsCompleteAndFeasible(lns_fragment, /*logger=*/nullptr); } // If we use a hint, we will restrict the objective to be <= to the one @@ -1398,7 +1398,8 @@ class LnsSolver : public SubSolver { if (lns_parameters_.debug_crash_if_presolve_breaks_hint() && hint_feasible_before_presolve && - !TestSolutionHintForFeasibility(lns_fragment, /*logger=*/nullptr)) { + !SolutionHintIsCompleteAndFeasible(lns_fragment, + /*logger=*/nullptr)) { LOG(FATAL) << "Presolve broke a feasible LNS hint. The model name is '" << lns_fragment.name() << "' (use the --cp_model_dump_submodels flag to dump it)."; @@ -2390,7 +2391,7 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { bool hint_feasible_before_presolve = false; if (!context->ModelIsUnsat()) { hint_feasible_before_presolve = - TestSolutionHintForFeasibility(model_proto, logger); + SolutionHintIsCompleteAndFeasible(model_proto, logger); } // If the objective was a floating point one, do some postprocessing on the @@ -2726,11 +2727,11 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { // different. bool hint_feasible_after_presolve; if (!params.enumerate_all_solutions()) { - hint_feasible_after_presolve = TestSolutionHintForFeasibility( + hint_feasible_after_presolve = SolutionHintIsCompleteAndFeasible( *new_cp_model_proto, logger, shared_response_manager); } else { hint_feasible_after_presolve = - TestSolutionHintForFeasibility(*new_cp_model_proto, logger, nullptr); + SolutionHintIsCompleteAndFeasible(*new_cp_model_proto, logger, nullptr); } if (hint_feasible_before_presolve && !hint_feasible_after_presolve && diff --git a/ortools/sat/cp_model_solver_fuzz.cc b/ortools/sat/cp_model_solver_fuzz.cc index 35d69966d0..f2cad86d0e 100644 --- a/ortools/sat/cp_model_solver_fuzz.cc +++ b/ortools/sat/cp_model_solver_fuzz.cc @@ -32,8 +32,9 @@ std::string GetTestDataDir() { void Solve(const CpModelProto& proto) { const CpSolverResponse response = - operations_research::sat::SolveWithParameters(proto, - "max_time_in_seconds: 4.0"); + operations_research::sat::SolveWithParameters( + proto, + "max_time_in_seconds: 4.0,debug_crash_if_presolve_breaks_hint:true"); const CpSolverResponse response_no_presolve = operations_research::sat::SolveWithParameters( diff --git a/ortools/sat/cp_model_utils.h b/ortools/sat/cp_model_utils.h index 59d901cba6..898d9455b0 100644 --- a/ortools/sat/cp_model_utils.h +++ b/ortools/sat/cp_model_utils.h @@ -34,6 +34,7 @@ #include "ortools/base/hash.h" #include "ortools/base/options.h" #include "ortools/sat/cp_model.pb.h" +#include "ortools/util/bitset.h" #include "ortools/util/sorted_interval_list.h" namespace operations_research { @@ -99,12 +100,22 @@ std::vector UsedVariables(const ConstraintProto& ct); // Returns the sorted list of interval used by a constraint. std::vector UsedIntervals(const ConstraintProto& ct); -// Insert variables in a constraint into a set. -template -void InsertVariablesFromConstraint(const CpModelProto& model_proto, int index, - Set& output) { +// Insert/Remove variables from an interval constraint into a bitset. +inline void InsertVariablesFromInterval(const CpModelProto& model_proto, + int index, Bitset64& output) { const ConstraintProto& ct = model_proto.constraints(index); - for (const int var : UsedVariables(ct)) output.insert(var); + for (const int ref : ct.enforcement_literal()) output.Set(PositiveRef(ref)); + for (const int var : ct.interval().start().vars()) output.Set(var); + for (const int var : ct.interval().size().vars()) output.Set(var); + for (const int var : ct.interval().end().vars()) output.Set(var); +} +inline void RemoveVariablesFromInterval(const CpModelProto& model_proto, + int index, Bitset64& output) { + const ConstraintProto& ct = model_proto.constraints(index); + for (const int ref : ct.enforcement_literal()) output.Clear(PositiveRef(ref)); + for (const int var : ct.interval().start().vars()) output.Clear(var); + for (const int var : ct.interval().size().vars()) output.Clear(var); + for (const int var : ct.interval().end().vars()) output.Clear(var); } // Returns true if a proto.domain() contain the given value. diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index ef74f4ecc4..a81853c2b2 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -66,51 +66,30 @@ IntegerVariable CreateVariableWithTightDomain( return integer_trail->AddIntegerVariable(min, max); } -// TODO(user): Use the faster variable only version if all expressions reduce -// to a single variable? -IntegerVariable CreateVariableEqualToMinOf( +IntegerVariable CreateVariableAtOrAboveMinOf( absl::Span exprs, Model* model) { - std::vector converted; - converted.reserve(exprs.size()); - for (const AffineExpression& affine : exprs) { - LinearExpression e; - e.offset = affine.constant; - if (affine.var != kNoIntegerVariable) { - e.vars.push_back(affine.var); - e.coeffs.push_back(affine.coeff); - } - converted.push_back(e); - } - const IntegerVariable var = CreateVariableWithTightDomain(exprs, model); - LinearExpression target; - target.vars.push_back(var); - target.coeffs.push_back(IntegerValue(1)); - AddIsEqualToMinOf(target, std::move(converted), model); + auto* constraint = new MinPropagator({exprs.begin(), exprs.end()}, var, + model->GetOrCreate()); + constraint->RegisterWith(model->GetOrCreate()); + model->TakeOwnership(constraint); + return var; } -IntegerVariable CreateVariableEqualToMaxOf( +IntegerVariable CreateVariableAtOrBelowMaxOf( absl::Span exprs, Model* model) { - std::vector converted; - converted.reserve(exprs.size()); - for (const AffineExpression& affine : exprs) { - // We take the negation of affine. - LinearExpression e; - e.offset = -affine.constant; - if (affine.var != kNoIntegerVariable) { - e.vars = {NegationOf(affine.var)}; - e.coeffs = {affine.coeff}; - } - converted.push_back(std::move(e)); - } + std::vector negated_exprs(exprs.begin(), exprs.end()); + for (AffineExpression& affine : negated_exprs) affine = affine.Negated(); + + const IntegerVariable var = + CreateVariableWithTightDomain(negated_exprs, model); + auto* constraint = new MinPropagator(std::move(negated_exprs), var, + model->GetOrCreate()); + constraint->RegisterWith(model->GetOrCreate()); + model->TakeOwnership(constraint); - const IntegerVariable var = CreateVariableWithTightDomain(exprs, model); - LinearExpression target; - target.vars = {NegationOf(var)}; - target.coeffs = {IntegerValue(1)}; - AddIsEqualToMinOf(target, std::move(converted), model); - return var; + return NegationOf(var); } // Add a cumulative relaxation. That is, on one dimension, it does not enforce @@ -118,11 +97,14 @@ IntegerVariable CreateVariableEqualToMaxOf( void AddDiffnCumulativeRelationOnX(SchedulingConstraintHelper* x, SchedulingConstraintHelper* y, Model* model) { + // Note that we only need one side! + // We want something <= max_end - min_start + // // TODO(user): Use conditional affine min/max !! const IntegerVariable min_start_var = - CreateVariableEqualToMinOf(y->Starts(), model); + CreateVariableAtOrAboveMinOf(y->Starts(), model); const IntegerVariable max_end_var = - CreateVariableEqualToMaxOf(y->Ends(), model); + CreateVariableAtOrBelowMaxOf(y->Ends(), model); // (max_end - min_start) >= capacity. auto* integer_trail = model->GetOrCreate(); diff --git a/ortools/sat/diffn_util.cc b/ortools/sat/diffn_util.cc index 02ba333c50..704e6faac1 100644 --- a/ortools/sat/diffn_util.cc +++ b/ortools/sat/diffn_util.cc @@ -1972,6 +1972,8 @@ std::vector> FindPartialRectangleIntersections( } const int max_y_index = cur_index + 1; + gtl::STLClearObject(&y_events); + std::sort(x_events.begin(), x_events.end()); IntegerValue prev_x = 0; prev_event = Event::kEnd; @@ -2008,6 +2010,7 @@ std::vector> FindPartialRectangleIntersections( } } + gtl::STLClearObject(&x_events); return FindPartialRectangleIntersectionsImpl( absl::MakeSpan(sorted_rectangles32), max_y_index); } diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index a6cbc8ab85..6c50250762 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -550,10 +550,12 @@ bool LevelZeroEquality::Propagate() { return true; } -MinPropagator::MinPropagator(const std::vector& vars, - IntegerVariable min_var, +MinPropagator::MinPropagator(std::vector vars, + AffineExpression min_var, IntegerTrail* integer_trail) - : vars_(vars), min_var_(min_var), integer_trail_(integer_trail) {} + : vars_(std::move(vars)), + min_var_(min_var), + integer_trail_(integer_trail) {} bool MinPropagator::Propagate() { if (vars_.empty()) return true; @@ -579,11 +581,11 @@ bool MinPropagator::Propagate() { // Propagation a) if (min > integer_trail_->LowerBound(min_var_)) { integer_reason_.clear(); - for (const IntegerVariable var : vars_) { - integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, min)); + for (const AffineExpression& var : vars_) { + integer_reason_.push_back(var.GreaterOrEqual(min)); } - if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(min_var_, min), - {}, integer_reason_)) { + if (!integer_trail_->SafeEnqueue(min_var_.GreaterOrEqual(min), + integer_reason_)) { return false; } } @@ -598,15 +600,13 @@ bool MinPropagator::Propagate() { // The reason is that all the other interval start after current_min_ub. // And that min_ub has its current value. integer_reason_.push_back(min_ub_literal); - for (const IntegerVariable var : vars_) { + for (const AffineExpression& var : vars_) { if (var == vars_[last_possible_min_interval]) continue; - integer_reason_.push_back( - IntegerLiteral::GreaterOrEqual(var, current_min_ub + 1)); + integer_reason_.push_back(var.GreaterOrEqual(current_min_ub + 1)); } - if (!integer_trail_->Enqueue( - IntegerLiteral::LowerOrEqual(vars_[last_possible_min_interval], - current_min_ub), - {}, integer_reason_)) { + if (!integer_trail_->SafeEnqueue( + vars_[last_possible_min_interval].LowerOrEqual(current_min_ub), + integer_reason_)) { return false; } } @@ -623,9 +623,11 @@ bool MinPropagator::Propagate() { // Almost the same as propagation b). integer_reason_.push_back(min_ub_literal); - for (const IntegerVariable var : vars_) { - integer_reason_.push_back( - IntegerLiteral::GreaterOrEqual(var, current_min_ub + 1)); + for (const AffineExpression& var : vars_) { + IntegerLiteral lit = var.GreaterOrEqual(current_min_ub + 1); + if (lit != IntegerLiteral::TrueLiteral()) { + integer_reason_.push_back(lit); + } } return integer_trail_->ReportConflict(integer_reason_); } @@ -635,7 +637,7 @@ bool MinPropagator::Propagate() { void MinPropagator::RegisterWith(GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); - for (const IntegerVariable& var : vars_) { + for (const AffineExpression& var : vars_) { watcher->WatchLowerBound(var, id); } watcher->WatchUpperBound(min_var_, id); diff --git a/ortools/sat/integer_expr.h b/ortools/sat/integer_expr.h index 2de0b04b8d..bb9096519a 100644 --- a/ortools/sat/integer_expr.h +++ b/ortools/sat/integer_expr.h @@ -220,8 +220,8 @@ class LevelZeroEquality : PropagatorInterface { // TODO(user): Implement a more efficient algorithm when the need arise. class MinPropagator : public PropagatorInterface { public: - MinPropagator(const std::vector& vars, - IntegerVariable min_var, IntegerTrail* integer_trail); + MinPropagator(std::vector vars, AffineExpression min_var, + IntegerTrail* integer_trail); // This type is neither copyable nor movable. MinPropagator(const MinPropagator&) = delete; @@ -231,8 +231,8 @@ class MinPropagator : public PropagatorInterface { void RegisterWith(GenericLiteralWatcher* watcher); private: - const std::vector vars_; - const IntegerVariable min_var_; + const std::vector vars_; + const AffineExpression min_var_; IntegerTrail* integer_trail_; std::vector integer_reason_; @@ -729,6 +729,8 @@ inline void AddIsEqualToMinOf(const LinearExpression& min_expr, builder.AddTerm(min_var, -1); LoadLinearConstraint(builder.Build(), model); } + + // Add for all i, min <= exprs[i]. for (const LinearExpression& expr : exprs) { LinearConstraintBuilder builder(0, kMaxIntegerValue); builder.AddLinearExpression(expr, 1); diff --git a/ortools/sat/integer_expr_test.cc b/ortools/sat/integer_expr_test.cc index faf12759d2..566206d017 100644 --- a/ortools/sat/integer_expr_test.cc +++ b/ortools/sat/integer_expr_test.cc @@ -434,6 +434,32 @@ TEST(LinMinMaxTest, LevelZeroPropagation) { EXPECT_BOUNDS_EQ(vars[2], 5, 8); } +TEST(AffineMinTest, LevelZeroPropagation) { + Model model; + std::vector vars{model.Add(NewIntegerVariable(4, 9)), + model.Add(NewIntegerVariable(2, 7)), + model.Add(NewIntegerVariable(3, 8))}; + const AffineExpression min = model.Add(NewIntegerVariable(-100, 100)); + + auto* constraint = + new MinPropagator(vars, min, model.GetOrCreate()); + constraint->RegisterWith(model.GetOrCreate()); + model.TakeOwnership(constraint); + + auto* integer_trail = model.GetOrCreate(); + EXPECT_EQ(SatSolver::FEASIBLE, model.GetOrCreate()->Solve()); + EXPECT_EQ(integer_trail->LowerBound(min), 2); + + EXPECT_TRUE(integer_trail->Enqueue(min.LowerOrEqual(2))); + EXPECT_EQ(SatSolver::FEASIBLE, model.GetOrCreate()->Solve()); + EXPECT_EQ(integer_trail->UpperBound(min), 2); + + // Vars 1 is the only candidate. + EXPECT_EQ(integer_trail->UpperBound(vars[0]), 9); + EXPECT_EQ(integer_trail->UpperBound(vars[1]), 2); + EXPECT_EQ(integer_trail->UpperBound(vars[2]), 8); +} + TEST(LinMinTest, OnlyOnePossibleCandidate) { Model model; std::vector vars{model.Add(NewIntegerVariable(4, 7)), diff --git a/ortools/sat/lp_utils.cc b/ortools/sat/lp_utils.cc index 4578764ac1..031e924c8a 100644 --- a/ortools/sat/lp_utils.cc +++ b/ortools/sat/lp_utils.cc @@ -67,7 +67,7 @@ void ScaleConstraint(absl::Span var_scaling, } } -void ApplyVarScaling(const std::vector& var_scaling, +void ApplyVarScaling(absl::Span var_scaling, MPModelProto* mp_model) { const int num_variables = mp_model->variable_size(); for (int i = 0; i < num_variables; ++i) { diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index 270372f78b..1929e71128 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -542,6 +542,10 @@ class PresolveContext { const auto it = objective_map_.find(var); return it == objective_map_.end() ? 0 : it->second; } + + // Returns true if the variables in the objective with a positive (resp. + // negative) coefficient can freely decrease (resp. increase) within their + // domain (if we ignore the other constraints). bool ObjectiveDomainIsConstraining() const { return objective_domain_is_constraining_; } diff --git a/ortools/sat/routing_cuts.cc b/ortools/sat/routing_cuts.cc index aa6fa9559f..0d162ad2bf 100644 --- a/ortools/sat/routing_cuts.cc +++ b/ortools/sat/routing_cuts.cc @@ -403,7 +403,7 @@ bool OutgoingCutHelper::TryBlossomSubsetCut( } // namespace void GenerateInterestingSubsets(int num_nodes, - const std::vector>& arcs, + absl::Span> arcs, int stop_at_num_components, std::vector* subset_data, std::vector>* subsets) { diff --git a/ortools/sat/routing_cuts.h b/ortools/sat/routing_cuts.h index 6dad882654..5039e3ae97 100644 --- a/ortools/sat/routing_cuts.h +++ b/ortools/sat/routing_cuts.h @@ -49,7 +49,7 @@ namespace sat { // Note that this is mainly a "symmetric" case algo, but it does still work for // the asymmetric case. void GenerateInterestingSubsets(int num_nodes, - const std::vector>& arcs, + absl::Span> arcs, int stop_at_num_components, std::vector* subset_data, std::vector>* subsets); diff --git a/ortools/sat/var_domination_test.cc b/ortools/sat/var_domination_test.cc index f354a3aa31..44590b87c1 100644 --- a/ortools/sat/var_domination_test.cc +++ b/ortools/sat/var_domination_test.cc @@ -975,11 +975,13 @@ TEST(DualBoundReductionTest, AddImplication) { // not(a) => not(b) and not(a) => not(c) should be added. ASSERT_EQ(context.working_model->constraints_size(), 3); const ConstraintProto expected_constraint_proto1 = - ParseTestProto(R"pb(enforcement_literal: -1 bool_and { literals: -2 })pb"); + ParseTestProto(R"pb(enforcement_literal: -1 + bool_and { literals: -2 })pb"); EXPECT_THAT(context.working_model->constraints(1), EqualsProto(expected_constraint_proto1)); const ConstraintProto expected_constraint_proto2 = - ParseTestProto(R"pb(enforcement_literal: -1 bool_and { literals: -3 })pb"); + ParseTestProto(R"pb(enforcement_literal: -1 + bool_and { literals: -3 })pb"); EXPECT_THAT(context.working_model->constraints(2), EqualsProto(expected_constraint_proto2)); EXPECT_EQ(context.DomainOf(0).ToString(), "[0]");