diff --git a/ortools/sat/clause.cc b/ortools/sat/clause.cc index 1b57101588..273a0f8965 100644 --- a/ortools/sat/clause.cc +++ b/ortools/sat/clause.cc @@ -194,7 +194,7 @@ SatClause* LiteralWatchers::ReasonClause(int trail_index) const { return reasons_[trail_index]; } -bool LiteralWatchers::AddClause(const std::vector& literals, +bool LiteralWatchers::AddClause(absl::Span literals, Trail* trail) { SatClause* clause = SatClause::Create(literals); clauses_.push_back(clause); @@ -999,7 +999,7 @@ std::vector BinaryImplicationGraph::ExpandAtMostOne( // ----- SatClause ----- // static -SatClause* SatClause::Create(const std::vector& literals) { +SatClause* SatClause::Create(absl::Span literals) { CHECK_GE(literals.size(), 2); SatClause* clause = reinterpret_cast( ::operator new(sizeof(SatClause) + literals.size() * sizeof(Literal))); diff --git a/ortools/sat/clause.h b/ortools/sat/clause.h index ebffd0e6df..7af3976a23 100644 --- a/ortools/sat/clause.h +++ b/ortools/sat/clause.h @@ -52,7 +52,7 @@ class SatClause { // treated separatly and never constructed. In practice, we do use // BinaryImplicationGraph for the clause of size 2, so this is mainly used for // size at least 3. - static SatClause* Create(const std::vector& literals); + static SatClause* Create(absl::Span literals); // Non-sized delete because this is a tail-padded class. void operator delete(void* p) { @@ -168,7 +168,7 @@ class LiteralWatchers : public SatPropagator { SatClause* ReasonClause(int trail_index) const; // Adds a new clause and perform initial propagation for this clause only. - bool AddClause(const std::vector& literals, Trail* trail); + bool AddClause(absl::Span literals, Trail* trail); // Same as AddClause() for a removable clause. This is only called on learned // conflict, so this should never have all its literal at false (CHECKED). diff --git a/ortools/sat/cp_constraints.cc b/ortools/sat/cp_constraints.cc index e221038c53..d23dd01360 100644 --- a/ortools/sat/cp_constraints.cc +++ b/ortools/sat/cp_constraints.cc @@ -77,17 +77,25 @@ void BooleanXorPropagator::RegisterWith(GenericLiteralWatcher* watcher) { } GreaterThanAtLeastOneOfPropagator::GreaterThanAtLeastOneOfPropagator( - IntegerVariable target_var, const std::vector& vars, - const std::vector& offsets, - const std::vector& selectors, Model* model) + IntegerVariable target_var, const absl::Span vars, + const absl::Span offsets, + const absl::Span selectors, + const absl::Span enforcements, Model* model) : target_var_(target_var), - vars_(vars), - offsets_(offsets), - selectors_(selectors), + vars_(vars.begin(), vars.end()), + offsets_(offsets.begin(), offsets.end()), + selectors_(selectors.begin(), selectors.end()), + enforcements_(enforcements.begin(), enforcements.end()), trail_(model->GetOrCreate()), integer_trail_(model->GetOrCreate()) {} bool GreaterThanAtLeastOneOfPropagator::Propagate() { + // TODO(user): In case of a conflict, we could push one of them to false if + // it is the only one. + for (const Literal l : enforcements_) { + if (!trail_->Assignment().LiteralIsTrue(l)) return true; + } + // Compute the min of the lower-bound for the still possible variables. // TODO(user): This could be optimized by keeping more info from the last // Propagate() calls. @@ -110,6 +118,9 @@ bool GreaterThanAtLeastOneOfPropagator::Propagate() { literal_reason_.clear(); integer_reason_.clear(); + for (const Literal l : enforcements_) { + literal_reason_.push_back(l.Negated()); + } for (int i = 0; i < vars_.size(); ++i) { if (trail_->Assignment().LiteralIsFalse(selectors_[i])) { literal_reason_.push_back(selectors_[i]); @@ -127,6 +138,7 @@ void GreaterThanAtLeastOneOfPropagator::RegisterWith( GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); for (const Literal l : selectors_) watcher->WatchLiteral(l.Negated(), id); + for (const Literal l : enforcements_) watcher->WatchLiteral(l, id); for (const IntegerVariable v : vars_) watcher->WatchLowerBound(v, id); } diff --git a/ortools/sat/cp_constraints.h b/ortools/sat/cp_constraints.h index 10472bd644..5dabae0cc3 100644 --- a/ortools/sat/cp_constraints.h +++ b/ortools/sat/cp_constraints.h @@ -65,13 +65,15 @@ class BooleanXorPropagator : public PropagatorInterface { // alternatives. // // This constraint take care of this case when no selectors[i] is chosen yet. +// +// This constraint support duplicate selectors. class GreaterThanAtLeastOneOfPropagator : public PropagatorInterface { public: - GreaterThanAtLeastOneOfPropagator(IntegerVariable target_var, - const std::vector& vars, - const std::vector& offsets, - const std::vector& selectors, - Model* model); + GreaterThanAtLeastOneOfPropagator( + IntegerVariable target_var, const absl::Span vars, + const absl::Span offsets, + const absl::Span selectors, + const absl::Span enforcements, Model* model); bool Propagate() final; void RegisterWith(GenericLiteralWatcher* watcher); @@ -81,6 +83,7 @@ class GreaterThanAtLeastOneOfPropagator : public PropagatorInterface { const std::vector vars_; const std::vector offsets_; const std::vector selectors_; + const std::vector enforcements_; Trail* trail_; IntegerTrail* integer_trail_; @@ -118,13 +121,27 @@ inline std::function LiteralXorIs( } inline std::function GreaterThanAtLeastOneOf( - IntegerVariable target_var, const std::vector& vars, - const std::vector& offsets, - const std::vector& selectors) { + IntegerVariable target_var, const absl::Span vars, + const absl::Span offsets, + const absl::Span selectors) { return [=](Model* model) { GreaterThanAtLeastOneOfPropagator* constraint = new GreaterThanAtLeastOneOfPropagator(target_var, vars, offsets, - selectors, model); + selectors, {}, model); + constraint->RegisterWith(model->GetOrCreate()); + model->TakeOwnership(constraint); + }; +} + +inline std::function GreaterThanAtLeastOneOf( + IntegerVariable target_var, const absl::Span vars, + const absl::Span offsets, + const absl::Span selectors, + const absl::Span enforcements) { + return [=](Model* model) { + GreaterThanAtLeastOneOfPropagator* constraint = + new GreaterThanAtLeastOneOfPropagator(target_var, vars, offsets, + selectors, enforcements, model); constraint->RegisterWith(model->GetOrCreate()); model->TakeOwnership(constraint); }; diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 8ad55fe450..7c194643c2 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -25,8 +25,13 @@ namespace operations_research { namespace sat { NeighborhoodGeneratorHelper::NeighborhoodGeneratorHelper( - CpModelProto const* model, bool focus_on_decision_variables) - : model_proto_(*model) { + CpModelProto const* model_proto, bool focus_on_decision_variables) + : model_proto_(*model_proto) { + UpdateHelperData(focus_on_decision_variables); +} + +void NeighborhoodGeneratorHelper::UpdateHelperData( + bool focus_on_decision_variables) { var_to_constraint_.resize(model_proto_.variables_size()); constraint_to_var_.resize(model_proto_.constraints_size()); for (int ct_index = 0; ct_index < model_proto_.constraints_size(); @@ -114,8 +119,9 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( } // TODO(user): force better objective? Note that this is already done when the - // hint above is sucessfully loaded (i.e. if it passes the presolve correctly) - // since the solver will try to find better solution than the current one. + // hint above is successfully loaded (i.e. if it passes the presolve + // correctly) since the solver will try to find better solution than the + // current one. return neighborhood; } diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index a8e0b0e622..1b506ce1d8 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -43,9 +43,12 @@ struct Neighborhood { // Thread-safe. class NeighborhoodGeneratorHelper { public: - NeighborhoodGeneratorHelper(CpModelProto const* model, + NeighborhoodGeneratorHelper(CpModelProto const* model_proto, bool focus_on_decision_variables); + // Updates private variables using the current 'model_proto_'. + void UpdateHelperData(bool focus_on_decision_variables); + // Returns the LNS fragment where the given variables are fixed to the value // they take in the given solution. Neighborhood FixGivenVariables( diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index f7e7666781..7bc3fc3ef3 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -984,12 +984,12 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto, relaxation.linear_constraints.resize(new_size); } - VLOG(2) << "num_full_encoding_relaxations: " << num_full_encoding_relaxations; - VLOG(2) << "num_partial_encoding_relaxations: " + VLOG(3) << "num_full_encoding_relaxations: " << num_full_encoding_relaxations; + VLOG(3) << "num_partial_encoding_relaxations: " << num_partial_encoding_relaxations; - VLOG(2) << relaxation.linear_constraints.size() + VLOG(3) << relaxation.linear_constraints.size() << " constraints in the LP relaxation."; - VLOG(2) << relaxation.cut_generators.size() << " cuts generators."; + VLOG(3) << relaxation.cut_generators.size() << " cuts generators."; // The bipartite graph of LP constraints might be disconnected: // make a partition of the variables into connected components. @@ -1134,7 +1134,7 @@ IntegerVariable AddLPConstraints(const CpModelProto& model_proto, // constraints have been added. for (auto* lp_constraint : lp_constraints) { lp_constraint->RegisterWith(m); - VLOG(2) << "LP constraint: " << lp_constraint->DimensionString() << "."; + VLOG(3) << "LP constraint: " << lp_constraint->DimensionString() << "."; } VLOG(2) << top_level_cp_terms.size() @@ -1287,19 +1287,19 @@ CpSolverResponse SolveCpModelInternal( Trail* trail = model->GetOrCreate(); const int old_num_fixed = trail->Index(); if (trail->Index() > old_num_fixed) { - VLOG(2) << "Constraint fixed " << trail->Index() - old_num_fixed + VLOG(3) << "Constraint fixed " << trail->Index() - old_num_fixed << " Boolean variable(s): " << ProtobufDebugString(ct); } } if (model->GetOrCreate()->IsModelUnsat()) { - VLOG(2) << "UNSAT during extraction (after adding '" + VLOG(3) << "UNSAT during extraction (after adding '" << ConstraintCaseName(ct.constraint_case()) << "'). " << ProtobufDebugString(ct); break; } } if (num_ignored_constraints > 0) { - VLOG(2) << num_ignored_constraints << " constraints were skipped."; + VLOG(3) << num_ignored_constraints << " constraints were skipped."; } if (!unsupported_types.empty()) { VLOG(1) << "There is unsuported constraints types in this model: "; @@ -1531,32 +1531,55 @@ CpSolverResponse SolveCpModelInternal( status = SolveProblemWithPortfolioSearch(decision_policies, {no_restart}, model); if (status == SatSolver::Status::FEASIBLE) { - solution_info = "hint"; - solution_observer(*model); - CHECK(SolutionIsFeasible(model_proto, - std::vector(response.solution().begin(), - response.solution().end()))); - if (!model_proto.has_objective()) { - if (parameters.enumerate_all_solutions()) { - model->Add( - ExcludeCurrentSolutionWithoutIgnoredVariableAndBacktrack()); - } else { - return response; + bool hint_is_valid = true; + if (model_proto.has_objective() && parameters.optimize_with_core()) { + // We need to fix the lower bound of the objective variable. + // If linearization_level = 0, the objective_var is not linked with + // the model. We recompute its value from scratch. + int64 objective_value = 0; + for (int i = 0; i < model_proto.objective().vars_size(); ++i) { + objective_value += model_proto.objective().coeffs(i) * + model->Get(LowerBound(mapping->Integer( + model_proto.objective().vars(i)))); } - } else { IntegerTrail* integer_trail = model->GetOrCreate(); - const IntegerValue current_internal_objective = - integer_trail->LowerBound(objective_var); - // Restrict the objective. - model->GetOrCreate()->Backtrack(0); if (!integer_trail->Enqueue( - IntegerLiteral::LowerOrEqual(objective_var, - current_internal_objective - 1), + IntegerLiteral::GreaterOrEqual(objective_var, + IntegerValue(objective_value)), {}, {})) { - response.set_best_objective_bound(response.objective_value()); - response.set_status(CpSolverStatus::OPTIMAL); - fill_response_statistics(); - return response; + hint_is_valid = false; + model->GetOrCreate()->Backtrack(0); + } + } + if (hint_is_valid) { + solution_info = "hint"; + solution_observer(*model); + num_solutions++; + CHECK(SolutionIsFeasible( + model_proto, std::vector(response.solution().begin(), + response.solution().end()))); + if (!model_proto.has_objective()) { + if (parameters.enumerate_all_solutions()) { + model->Add( + ExcludeCurrentSolutionWithoutIgnoredVariableAndBacktrack()); + } else { + return response; + } + } else { + IntegerTrail* integer_trail = model->GetOrCreate(); + IntegerValue current_internal_objective = + integer_trail->LowerBound(objective_var); + // Restrict the objective. + model->GetOrCreate()->Backtrack(0); + if (!integer_trail->Enqueue( + IntegerLiteral::LowerOrEqual(objective_var, + current_internal_objective - 1), + {}, {})) { + response.set_best_objective_bound(response.objective_value()); + response.set_status(CpSolverStatus::OPTIMAL); + fill_response_statistics(); + return response; + } } } } @@ -1584,8 +1607,8 @@ CpSolverResponse SolveCpModelInternal( } else { // Optimization problem. const CpObjectiveProto& obj = model_proto.objective(); - VLOG(2) << obj.vars_size() << " terms in the proto objective."; - VLOG(2) << "Initial num_bool: " << model->Get()->NumVariables(); + VLOG(3) << obj.vars_size() << " terms in the proto objective."; + VLOG(3) << "Initial num_bool: " << model->Get()->NumVariables(); if (parameters.optimize_with_core()) { std::vector linear_vars; @@ -1594,13 +1617,16 @@ CpSolverResponse SolveCpModelInternal( &linear_vars, &linear_coeffs); if (parameters.optimize_with_max_hs()) { status = MinimizeWithHittingSetAndLazyEncoding( - VLOG_IS_ON(2), objective_var, linear_vars, linear_coeffs, + VLOG_IS_ON(3), objective_var, linear_vars, linear_coeffs, next_decision, solution_observer, model); } else { status = MinimizeWithCoreAndLazyEncoding( - VLOG_IS_ON(2), objective_var, linear_vars, linear_coeffs, + VLOG_IS_ON(3), objective_var, linear_vars, linear_coeffs, next_decision, solution_observer, model); } + if (num_solutions > 0 && status == SatSolver::INFEASIBLE) { + status = SatSolver::FEASIBLE; + } } else { if (parameters.binary_search_num_conflicts() >= 0) { RestrictObjectiveDomainWithBinarySearch(objective_var, next_decision, @@ -1922,8 +1948,8 @@ CpSolverResponse SolveCpModelWithLNS( const bool focus_on_decision_variables = parameters->lns_focus_on_decision_variables(); - const NeighborhoodGeneratorHelper helper(&mutable_model_proto, - focus_on_decision_variables); + NeighborhoodGeneratorHelper helper(&mutable_model_proto, + focus_on_decision_variables); // For now we will just alternate between our possible neighborhoods. // @@ -1983,18 +2009,7 @@ CpSolverResponse SolveCpModelWithLNS( } } - // If we didn't see any progress recently, bump the time limit. - // TODO(user): Tune the logic and expose the parameters. - if (num_no_progress > 100) { - deterministic_time *= 1.1; - num_no_progress = 0; - } - return limit->LimitReached() || - response.objective_value() == response.best_objective_bound(); - }, - [&](int64 seed) { // Update the bounds on mutable model proto. - bool model_is_unsat = false; if (shared_bounds_manager != nullptr) { std::vector model_variables; std::vector new_lower_bounds; @@ -2007,20 +2022,35 @@ CpSolverResponse SolveCpModelWithLNS( const int var = model_variables[i]; const int64 new_lb = new_lower_bounds[i]; const int64 new_ub = new_upper_bounds[i]; - if (VLOG_IS_ON(2)) { + if (VLOG_IS_ON(3)) { const auto& domain = mutable_model_proto.variables(var).domain(); const int64 old_lb = domain.Get(0); const int64 old_ub = domain.Get(domain.size() - 1); - VLOG(2) << "Variable: " << var << " old domain: [" << old_lb + VLOG(3) << "Variable: " << var << " old domain: [" << old_lb << ", " << old_ub << "] new domain: [" << new_lb << ", " << new_ub << "]"; } if (!UpdateDomain(new_lb, new_ub, mutable_model_proto.mutable_variables(var))) { - model_is_unsat = true; + // Model is unsat. + return true; } } + if (!model_variables.empty()) { + helper.UpdateHelperData(focus_on_decision_variables); + } } + + // If we didn't see any progress recently, bump the time limit. + // TODO(user): Tune the logic and expose the parameters. + if (num_no_progress > 100) { + deterministic_time *= 1.1; + num_no_progress = 0; + } + return limit->LimitReached() || + response.objective_value() == response.best_objective_bound(); + }, + [&](int64 seed) { const int selected_generator = seed % generators.size(); AdaptiveParameterValue& difficulty = difficulties[selected_generator]; const double saved_difficulty = difficulty.value(); @@ -2052,24 +2082,21 @@ CpSolverResponse SolveCpModelWithLNS( // Presolve and solve the LNS fragment. CpSolverResponse local_response; - if (model_is_unsat) { - local_response.set_status(CpSolverStatus::INFEASIBLE); - } else { - CpModelProto mapping_proto; - std::vector postsolve_mapping; - PresolveOptions options; - options.log_info = VLOG_IS_ON(2); - options.parameters = *local_model.GetOrCreate(); - options.time_limit = local_model.GetOrCreate(); - PresolveCpModel(options, &local_problem, &mapping_proto, - &postsolve_mapping); - local_response = SolveCpModelInternal( - local_problem, /*is_real_solve=*/true, - [](const CpSolverResponse& response) {}, - /*shared_bounds_manager=*/nullptr, wall_timer, &local_model); - PostsolveResponse(model_proto, mapping_proto, postsolve_mapping, - wall_timer, &local_response); - } + + CpModelProto mapping_proto; + std::vector postsolve_mapping; + PresolveOptions options; + options.log_info = VLOG_IS_ON(3); + options.parameters = *local_model.GetOrCreate(); + options.time_limit = local_model.GetOrCreate(); + PresolveCpModel(options, &local_problem, &mapping_proto, + &postsolve_mapping); + local_response = SolveCpModelInternal( + local_problem, /*is_real_solve=*/true, + [](const CpSolverResponse& response) {}, + /*shared_bounds_manager=*/nullptr, wall_timer, &local_model); + PostsolveResponse(model_proto, mapping_proto, postsolve_mapping, + wall_timer, &local_response); const bool neighborhood_is_reduced = neighborhood.is_reduced; return [neighborhood_is_reduced, &num_no_progress, &model_proto, diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index 079d26f1d7..f9ef3337f0 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -205,6 +205,7 @@ Literal IntegerEncoder::GetOrCreateAssociatedLiteral(IntegerLiteral i_lit) { ++num_created_variables_; const Literal literal(sat_solver_->NewBooleanVariable(), true); AssociateToIntegerLiteral(literal, new_lit); + CHECK(!sat_solver_->Assignment().LiteralIsAssigned(literal)); return literal; } @@ -218,9 +219,25 @@ Literal IntegerEncoder::GetOrCreateLiteralAssociatedToEquality( } } + // Check for trivial true/false literal to avoid creating variable for no + // reasons. + const Domain& domain = (*domains_)[var]; + if (!domain.Contains(value.value())) { + AssociateToIntegerEqualValue(GetFalseLiteral(), var, value); + return GetFalseLiteral(); + } + if (value == domain.Min() && value == domain.Max()) { + AssociateToIntegerEqualValue(GetTrueLiteral(), var, value); + return GetTrueLiteral(); + } + ++num_created_variables_; const Literal literal(sat_solver_->NewBooleanVariable(), true); AssociateToIntegerEqualValue(literal, var, value); + + // TODO(user): on some problem the check below fail. We should probably + // make sure that we don't create extra fixed Boolean variable for no reason. + // CHECK(!sat_solver_->Assignment().LiteralIsAssigned(literal)); return literal; } diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index 21fd88d5b5..4775fb6f79 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -34,9 +34,11 @@ LiteralIndex AtMinValue(IntegerVariable var, Model* model) { const IntegerValue lb = integer_trail->LowerBound(var); DCHECK_LE(lb, integer_trail->UpperBound(var)); if (lb == integer_trail->UpperBound(var)) return kNoLiteralIndex; - return integer_encoder - ->GetOrCreateAssociatedLiteral(IntegerLiteral::LowerOrEqual(var, lb)) - .Index(); + + const Literal result = integer_encoder->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(var, lb)); + CHECK(!model->GetOrCreate()->Assignment().LiteralIsAssigned(result)); + return result.Index(); } LiteralIndex GreaterOrEqualToMiddleValue(IntegerVariable var, Model* model) { @@ -245,8 +247,10 @@ std::function SatSolverHeuristic(Model* model) { SatDecisionPolicy* decision_policy = model->GetOrCreate(); return [sat_solver, trail, decision_policy] { const bool all_assigned = trail->Index() == sat_solver->NumVariables(); - return all_assigned ? kNoLiteralIndex - : decision_policy->NextBranch().Index(); + if (all_assigned) return kNoLiteralIndex; + const Literal result = decision_policy->NextBranch(); + CHECK(!sat_solver->Assignment().LiteralIsAssigned(result)); + return result.Index(); }; } diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index e5dc17d04d..2b21483c69 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -1728,7 +1728,7 @@ LiteralIndex LinearProgrammingConstraint::LPReducedCostAverageDecision() { } if (selected_index == -1) return kNoLiteralIndex; - const IntegerVariable var = this->integer_variables_[selected_index]; + const IntegerVariable var = integer_variables_[selected_index]; // If ceil(value) is current upper bound, try var == upper bound first. // Guarding with >= prevents numerical problems. @@ -1738,9 +1738,10 @@ LiteralIndex LinearProgrammingConstraint::LPReducedCostAverageDecision() { const IntegerValue value_ceil( std::ceil(this->GetSolutionValue(var) - kCpEpsilon)); if (value_ceil >= ub) { - return integer_encoder_ - ->GetOrCreateAssociatedLiteral(IntegerLiteral::GreaterOrEqual(var, ub)) - .Index(); + const Literal result = integer_encoder_->GetOrCreateAssociatedLiteral( + IntegerLiteral::GreaterOrEqual(var, ub)); + CHECK(!trail_->Assignment().LiteralIsAssigned(result)); + return result.Index(); } // If floor(value) is current lower bound, try var == lower bound first. @@ -1749,24 +1750,26 @@ LiteralIndex LinearProgrammingConstraint::LPReducedCostAverageDecision() { const IntegerValue value_floor( std::floor(this->GetSolutionValue(var) + kCpEpsilon)); if (value_floor <= lb) { - return integer_encoder_ - ->GetOrCreateAssociatedLiteral(IntegerLiteral::LowerOrEqual(var, lb)) - .Index(); + const Literal result = integer_encoder_->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(var, lb)); + CHECK(!trail_->Assignment().LiteralIsAssigned(result)) + << " " << lb << " " << ub; + return result.Index(); } // Here lb < value_floor <= value_ceil < ub. // Try the most promising split between var <= floor or var >= ceil. if (sum_cost_down_[selected_index] / num_cost_down_[selected_index] < sum_cost_up_[selected_index] / num_cost_up_[selected_index]) { - return integer_encoder_ - ->GetOrCreateAssociatedLiteral( - IntegerLiteral::LowerOrEqual(var, value_floor)) - .Index(); + const Literal result = integer_encoder_->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(var, value_floor)); + CHECK(!trail_->Assignment().LiteralIsAssigned(result)); + return result.Index(); } else { - return integer_encoder_ - ->GetOrCreateAssociatedLiteral( - IntegerLiteral::GreaterOrEqual(var, value_ceil)) - .Index(); + const Literal result = integer_encoder_->GetOrCreateAssociatedLiteral( + IntegerLiteral::GreaterOrEqual(var, value_ceil)); + CHECK(!trail_->Assignment().LiteralIsAssigned(result)); + return result.Index(); } } diff --git a/ortools/sat/precedences.cc b/ortools/sat/precedences.cc index 5b51f6deb3..735e2cbc7b 100644 --- a/ortools/sat/precedences.cc +++ b/ortools/sat/precedences.cc @@ -19,6 +19,7 @@ #include "ortools/base/cleanup.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" +#include "ortools/sat/clause.h" #include "ortools/sat/cp_constraints.h" namespace operations_research { @@ -720,11 +721,92 @@ bool PrecedencesPropagator::BellmanFordTarjan(Trail* trail) { return true; } -void PrecedencesPropagator::AddGreaterThanAtLeastOneOfConstraints( - Model* model) { - VLOG(1) << "Detecting GreaterThanAtLeastOneOf() constraints..."; +int PrecedencesPropagator::AddGreaterThanAtLeastOneOfConstraintsFromClause( + const absl::Span clause, Model* model) { + CHECK_EQ(model->GetOrCreate()->CurrentDecisionLevel(), 0); + if (clause.size() < 2) return 0; + + // Collect all arcs impacted by this clause. + std::vector infos; + for (const Literal l : clause) { + if (l.Index() >= literal_to_new_impacted_arcs_.size()) continue; + for (const ArcIndex arc_index : literal_to_new_impacted_arcs_[l.Index()]) { + const ArcInfo& arc = arcs_[arc_index]; + if (arc.presence_literals.size() != 1) continue; + + // TODO(user): Support variable offset. + if (arc.offset_var != kNoIntegerVariable) continue; + infos.push_back(arc); + } + } + if (infos.size() <= 1) return 0; + + // Stable sort by head_var so that for a same head_var, the entry are sorted + // by Literal as they appear in clause. + std::stable_sort(infos.begin(), infos.end(), + [](const ArcInfo& a, const ArcInfo& b) { + return a.head_var < b.head_var; + }); + + // We process ArcInfo with the same head_var toghether. + int num_added_constraints = 0; auto* solver = model->GetOrCreate(); + for (int i = 0; i < infos.size();) { + const int start = i; + const IntegerVariable head_var = infos[start].head_var; + for (i++; i < infos.size() && infos[i].head_var == head_var; ++i) { + } + const absl::Span arcs(&infos[start], i - start); + + // Skip single arcs since it will already be fully propagated. + if (arcs.size() < 2) continue; + + // Heuristic. Look for full or almost full clauses. We could add + // GreaterThanAtLeastOneOf() with more enforcement literals. TODO(user): + // experiments. + if (arcs.size() + 1 < clause.size()) continue; + + std::vector vars; + std::vector offsets; + std::vector selectors; + std::vector enforcements; + + int j = 0; + for (const Literal l : clause) { + bool added = false; + for (; j < arcs.size() && l == arcs[j].presence_literals.front(); ++j) { + added = true; + vars.push_back(arcs[j].tail_var); + offsets.push_back(arcs[j].offset); + + // Note that duplicate selector are supported. + // + // TODO(user): If we support variable offset, we should regroup the arcs + // into one (tail + offset <= head) though, instead of having too + // identical entries. + selectors.push_back(l); + } + if (!added) { + enforcements.push_back(l.Negated()); + } + } + + // No point adding a constraint if there is not at least two different + // literals in selectors. + if (enforcements.size() + 1 == clause.size()) continue; + + ++num_added_constraints; + model->Add(GreaterThanAtLeastOneOf(head_var, vars, offsets, selectors, + enforcements)); + if (!solver->FinishPropagation()) return num_added_constraints; + } + return num_added_constraints; +} + +int PrecedencesPropagator:: + AddGreaterThanAtLeastOneOfConstraintsWithClauseAutoDetection(Model* model) { auto* time_limit = model->GetOrCreate(); + auto* solver = model->GetOrCreate(); // Fill the set of incoming conditional arcs for each variables. gtl::ITIVector> incoming_arcs_; @@ -745,13 +827,13 @@ void PrecedencesPropagator::AddGreaterThanAtLeastOneOfConstraints( int num_added_constraints = 0; for (IntegerVariable target(0); target < incoming_arcs_.size(); ++target) { if (incoming_arcs_[target].size() <= 1) continue; - if (time_limit->LimitReached()) return; + if (time_limit->LimitReached()) return num_added_constraints; // Detect set of incoming arcs for which at least one must be present. // TODO(user): Find more than one disjoint set of incoming arcs. // TODO(user): call MinimizeCoreWithPropagation() on the clause. solver->Backtrack(0); - if (solver->IsModelUnsat()) return; + if (solver->IsModelUnsat()) return num_added_constraints; std::vector clause; for (const ArcIndex arc_index : incoming_arcs_[target]) { const Literal literal = arcs_[arc_index].presence_literals.front(); @@ -759,7 +841,7 @@ void PrecedencesPropagator::AddGreaterThanAtLeastOneOfConstraints( const int old_level = solver->CurrentDecisionLevel(); solver->EnqueueDecisionAndBacktrackOnConflict(literal.Negated()); - if (solver->IsModelUnsat()) return; + if (solver->IsModelUnsat()) return num_added_constraints; const int new_level = solver->CurrentDecisionLevel(); if (new_level <= old_level) { clause = solver->GetLastIncompatibleDecisions(); @@ -791,12 +873,48 @@ void PrecedencesPropagator::AddGreaterThanAtLeastOneOfConstraints( selectors.push_back(Literal(arcs_[a].presence_literals.front())); } model->Add(GreaterThanAtLeastOneOf(target, vars, offsets, selectors)); - if (!solver->FinishPropagation()) return; + if (!solver->FinishPropagation()) return num_added_constraints; } } + return num_added_constraints; +} + +int PrecedencesPropagator::AddGreaterThanAtLeastOneOfConstraints(Model* model) { + VLOG(1) << "Detecting GreaterThanAtLeastOneOf() constraints..."; + auto* time_limit = model->GetOrCreate(); + auto* solver = model->GetOrCreate(); + auto* clauses = model->GetOrCreate(); + int num_added_constraints = 0; + + // We have two possible approaches. For now, we prefer the first one except if + // there is too many clauses in the problem. + // + // TODO(user): Do more extensive experiment. Remove the second approach as + // it is more time consuming? or identify when it make sense. Note that the + // first approach also allows to use "incomplete" at least one between arcs. + if (clauses->AllClausesInCreationOrder().size() < 1e6) { + // TODO(user): This does not take into account clause of size 2 since they + // are stored in the BinaryImplicationGraph instead. Some ideas specific + // to size 2: + // - There can be a lot of such clauses, but it might be nice to consider + // them. we need to experiments. + // - The automatic clause detection might be a better approach and it + // could be combined with probing. + for (const SatClause* clause : clauses->AllClausesInCreationOrder()) { + if (time_limit->LimitReached()) return num_added_constraints; + if (solver->IsModelUnsat()) return num_added_constraints; + num_added_constraints += AddGreaterThanAtLeastOneOfConstraintsFromClause( + clause->AsSpan(), model); + } + } else { + num_added_constraints += + AddGreaterThanAtLeastOneOfConstraintsWithClauseAutoDetection(model); + } + VLOG(1) << "Added " << num_added_constraints << " GreaterThanAtLeastOneOf() constraints."; + return num_added_constraints; } } // namespace sat diff --git a/ortools/sat/precedences.h b/ortools/sat/precedences.h index 049d932ad3..de51e094e6 100644 --- a/ortools/sat/precedences.h +++ b/ortools/sat/precedences.h @@ -118,18 +118,30 @@ class PrecedencesPropagator : public SatPropagator, PropagatorInterface { // Advanced usage. To be called once all the constraints have been added to // the model. This will loop over all "node" in this class, and if one of its // optional incoming arcs must be chosen, it will add a corresponding - // GreaterThanAtLeastOneOfConstraint(). Note that this might be a bit slow as - // it relies on the propagation engine to detect clauses between incoming arcs - // presence literals. + // GreaterThanAtLeastOneOfConstraint(). Returns the number of added + // constraint. // // TODO(user): This can be quite slow, add some kind of deterministic limit // so that we can use it all the time. - void AddGreaterThanAtLeastOneOfConstraints(Model* model); + int AddGreaterThanAtLeastOneOfConstraints(Model* model); private: DEFINE_INT_TYPE(ArcIndex, int); DEFINE_INT_TYPE(OptionalArcIndex, int); + // Given an existing clause, sees if it can be used to add "greater than at + // least one of" type of constraints. Returns the number of such constraint + // added. + int AddGreaterThanAtLeastOneOfConstraintsFromClause( + const absl::Span clause, Model* model); + + // Another approach for AddGreaterThanAtLeastOneOfConstraints(), this one + // might be a bit slow as it relies on the propagation engine to detect + // clauses between incoming arcs presence literals. + // Returns the number of added constraints. + int AddGreaterThanAtLeastOneOfConstraintsWithClauseAutoDetection( + Model* model); + // Information about an individual arc. struct ArcInfo { IntegerVariable tail_var; diff --git a/ortools/sat/sat_solver.cc b/ortools/sat/sat_solver.cc index 88d2b9e6c7..bbf1699682 100644 --- a/ortools/sat/sat_solver.cc +++ b/ortools/sat/sat_solver.cc @@ -187,7 +187,7 @@ bool SatSolver::AddProblemClause(absl::Span literals) { &tmp_pb_constraint_); } -bool SatSolver::AddProblemClauseInternal(const std::vector& literals) { +bool SatSolver::AddProblemClauseInternal(absl::Span literals) { SCOPED_TIME_STAT(&stats_); CHECK_EQ(CurrentDecisionLevel(), 0); @@ -321,7 +321,14 @@ bool SatSolver::AddLinearConstraint(bool use_lower_bound, return SetModelUnsat(); } } - if (!Propagate()) return SetModelUnsat(); + + // Tricky: The PropagationIsDone() condition shouldn't change anything for a + // pure SAT problem, however in the CP-SAT context, calling Propagate() can + // tigger computation (like the LP) even if no domain changed since the last + // call. We do not want to do that. + if (!PropagationIsDone() && !Propagate()) { + return SetModelUnsat(); + } return true; } diff --git a/ortools/sat/sat_solver.h b/ortools/sat/sat_solver.h index 90109ce736..a1b688c5aa 100644 --- a/ortools/sat/sat_solver.h +++ b/ortools/sat/sat_solver.h @@ -512,7 +512,7 @@ class SatSolver { // Add a problem clause. The clause is assumed to be "cleaned", that is no // duplicate variables (not strictly required) and not empty. - bool AddProblemClauseInternal(const std::vector& literals); + bool AddProblemClauseInternal(absl::Span literals); // This is used by all the Add*LinearConstraint() functions. It detects // infeasible/trivial constraints or clause constraints and takes the proper