diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index 20f6b7836a..5c23c85401 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -814,6 +814,7 @@ cc_library( ":linear_constraint", ":model", ":sat_base", + ":synchronization", "//ortools/base", "//ortools/base:strong_vector", "//ortools/util:bitset", @@ -1031,7 +1032,7 @@ cc_library( ":linear_constraint", ":linear_programming_constraint", ":model", - ":presolve_util", + ":presolve_util", ":routing_cuts", ":sat_base", ":sat_parameters_cc_proto", diff --git a/ortools/sat/cp_model_checker.cc b/ortools/sat/cp_model_checker.cc index b308c84410..5efd598fce 100644 --- a/ortools/sat/cp_model_checker.cc +++ b/ortools/sat/cp_model_checker.cc @@ -1075,8 +1075,8 @@ namespace { class ConstraintChecker { public: - explicit ConstraintChecker(const std::vector& variable_values) - : variable_values_(variable_values) {} + explicit ConstraintChecker(absl::Span variable_values) + : variable_values_(variable_values.begin(), variable_values.end()) {} bool LiteralIsTrue(int l) const { if (l >= 0) return variable_values_[l] != 0; @@ -1524,13 +1524,13 @@ class ConstraintChecker { } private: - std::vector variable_values_; + const std::vector variable_values_; }; } // namespace bool SolutionIsFeasible(const CpModelProto& model, - const std::vector& variable_values, + absl::Span variable_values, const CpModelProto* mapping_proto, const std::vector* postsolve_mapping) { if (variable_values.size() != model.variables_size()) { diff --git a/ortools/sat/cp_model_checker.h b/ortools/sat/cp_model_checker.h index 1a0f6e7913..d97f01b54c 100644 --- a/ortools/sat/cp_model_checker.h +++ b/ortools/sat/cp_model_checker.h @@ -63,7 +63,7 @@ bool PossibleIntegerOverflow(const CpModelProto& model, // The last two arguments are optional and help debugging a failing constraint // due to presolve. bool SolutionIsFeasible(const CpModelProto& model, - const std::vector& variable_values, + absl::Span variable_values, const CpModelProto* mapping_proto = nullptr, const std::vector* postsolve_mapping = nullptr); diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 80c6e749d7..b4ab905877 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -1161,12 +1161,8 @@ void NeighborhoodGenerator::Synchronize() { // to a better "new objective" if the base solution wasn't the best one. // // This might not be a final solution, but it does work ok for now. - const IntegerValue best_objective_improvement = - IsRelaxationGenerator() - ? IntegerValue(CapSub(data.new_objective_bound.value(), - data.initial_best_objective_bound.value())) - : IntegerValue(CapSub(data.initial_best_objective.value(), - data.new_objective.value())); + const IntegerValue best_objective_improvement = IntegerValue(CapSub( + data.initial_best_objective.value(), data.new_objective.value())); if (best_objective_improvement > 0) { num_consecutive_non_improving_calls_ = 0; } else { @@ -1867,166 +1863,5 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate( return neighborhood; } -Neighborhood ConsecutiveConstraintsRelaxationNeighborhoodGenerator::Generate( - const CpSolverResponse& /*initial_solution*/, double difficulty, - absl::BitGenRef random) { - std::vector removable_constraints; - const int num_constraints = helper_.ModelProto().constraints_size(); - removable_constraints.reserve(num_constraints); - for (int c = 0; c < num_constraints; ++c) { - // Removing intervals is not easy because other constraint might require - // them, so for now, we don't remove them. - if (helper_.ModelProto().constraints(c).constraint_case() == - ConstraintProto::kInterval) { - continue; - } - removable_constraints.push_back(c); - } - - const int target_size = - std::round((1.0 - difficulty) * removable_constraints.size()); - - const int random_start_index = - absl::Uniform(random, 0, removable_constraints.size()); - std::vector removed_constraints; - removed_constraints.reserve(target_size); - int c = random_start_index; - while (removed_constraints.size() < target_size) { - removed_constraints.push_back(removable_constraints[c]); - ++c; - if (c == removable_constraints.size()) { - c = 0; - } - } - - return helper_.RemoveMarkedConstraints(removed_constraints); -} - -WeightedRandomRelaxationNeighborhoodGenerator:: - WeightedRandomRelaxationNeighborhoodGenerator( - NeighborhoodGeneratorHelper const* helper, const std::string& name) - : NeighborhoodGenerator(name, helper) { - std::vector removable_constraints; - const int num_constraints = helper_.ModelProto().constraints_size(); - constraint_weights_.reserve(num_constraints); - // TODO(user): Experiment with different starting weights. - for (int c = 0; c < num_constraints; ++c) { - switch (helper_.ModelProto().constraints(c).constraint_case()) { - case ConstraintProto::kCumulative: - case ConstraintProto::kAllDiff: - case ConstraintProto::kElement: - case ConstraintProto::kRoutes: - case ConstraintProto::kCircuit: - constraint_weights_.push_back(3.0); - num_removable_constraints_++; - break; - case ConstraintProto::kBoolOr: - case ConstraintProto::kBoolAnd: - case ConstraintProto::kBoolXor: - case ConstraintProto::kIntProd: - case ConstraintProto::kIntDiv: - case ConstraintProto::kIntMod: - case ConstraintProto::kLinMax: - case ConstraintProto::kNoOverlap: - case ConstraintProto::kNoOverlap2D: - constraint_weights_.push_back(2.0); - num_removable_constraints_++; - break; - case ConstraintProto::kLinear: - case ConstraintProto::kTable: - case ConstraintProto::kAutomaton: - case ConstraintProto::kInverse: - case ConstraintProto::kReservoir: - case ConstraintProto::kAtMostOne: - case ConstraintProto::kExactlyOne: - constraint_weights_.push_back(1.0); - num_removable_constraints_++; - break; - case ConstraintProto::CONSTRAINT_NOT_SET: - case ConstraintProto::kDummyConstraint: - case ConstraintProto::kInterval: - // Removing intervals is not easy because other constraint might require - // them, so for now, we don't remove them. - constraint_weights_.push_back(0.0); - break; - } - } -} - -void WeightedRandomRelaxationNeighborhoodGenerator:: - AdditionalProcessingOnSynchronize(const SolveData& solve_data) { - const IntegerValue best_objective_improvement = - solve_data.new_objective_bound - solve_data.initial_best_objective_bound; - - const std::vector& removed_constraints = - removed_constraints_[solve_data.neighborhood_id]; - - // Heuristic: We change the weights of the removed constraints if the - // neighborhood is solved (status is OPTIMAL or INFEASIBLE) or we observe an - // improvement in objective bounds. Otherwise we assume that the - // difficulty/time wasn't right for us to record feedbacks. - // - // If the objective bounds are improved, we bump up the weights. If the - // objective bounds are worse and the problem status is OPTIMAL, we bump down - // the weights. Otherwise if the new objective bounds are same as current - // bounds (which happens a lot on some instances), we do not update the - // weights as we do not have a clear signal whether the constraints removed - // were good choices or not. - // TODO(user): We can improve this heuristic with more experiments. - if (best_objective_improvement > 0) { - // Bump up the weights of all removed constraints. - for (int c : removed_constraints) { - if (constraint_weights_[c] <= 90.0) { - constraint_weights_[c] += 10.0; - } else { - constraint_weights_[c] = 100.0; - } - } - } else if (solve_data.status == CpSolverStatus::OPTIMAL && - best_objective_improvement < 0) { - // Bump down the weights of all removed constraints. - for (int c : removed_constraints) { - if (constraint_weights_[c] > 0.5) { - constraint_weights_[c] -= 0.5; - } - } - } - removed_constraints_.erase(solve_data.neighborhood_id); -} - -Neighborhood WeightedRandomRelaxationNeighborhoodGenerator::Generate( - const CpSolverResponse& /*initial_solution*/, double difficulty, - absl::BitGenRef random) { - const int target_size = - std::round((1.0 - difficulty) * num_removable_constraints_); - - std::vector removed_constraints; - - // Generate a random number between (0,1) = u[i] and use score[i] = - // u[i]^(1/w[i]) and then select top k items with largest scores. - // Reference: https://utopia.duth.gr/~pefraimi/research/data/2007EncOfAlg.pdf - std::vector> constraint_removal_scores; - std::uniform_real_distribution random_var(0.0, 1.0); - for (int c = 0; c < constraint_weights_.size(); ++c) { - if (constraint_weights_[c] <= 0) continue; - const double u = random_var(random); - const double score = std::pow(u, (1 / constraint_weights_[c])); - constraint_removal_scores.push_back({score, c}); - } - std::sort(constraint_removal_scores.rbegin(), - constraint_removal_scores.rend()); - for (int i = 0; i < target_size; ++i) { - removed_constraints.push_back(constraint_removal_scores[i].second); - } - - Neighborhood result = helper_.RemoveMarkedConstraints(removed_constraints); - - absl::MutexLock mutex_lock(&generator_mutex_); - result.id = next_available_id_; - next_available_id_++; - removed_constraints_.insert({result.id, removed_constraints}); - return result; -} - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index 9629e36e75..57d7921715 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -334,10 +334,6 @@ class NeighborhoodGenerator { // Returns true if the neighborhood generator can generate a neighborhood. virtual bool ReadyToGenerate() const; - // Returns true if the neighborhood generator generates relaxation of the - // given problem. - virtual bool IsRelaxationGenerator() const { return false; } - // Uses UCB1 algorithm to compute the score (Multi armed bandit problem). // Details are at // https://lilianweng.github.io/lil-log/2018/01/23/the-multi-armed-bandit-problem-and-its-solutions.html. @@ -350,10 +346,6 @@ class NeighborhoodGenerator { // Adds solve data about one "solved" neighborhood. struct SolveData { - // Neighborhood Id. Used to identify the neighborhood by a generator. - // Currently only used by WeightedRandomRelaxationNeighborhoodGenerator. - int64_t neighborhood_id = 0; - // The status of the sub-solve. CpSolverStatus status = CpSolverStatus::UNKNOWN; @@ -379,22 +371,14 @@ class NeighborhoodGenerator { IntegerValue base_objective = IntegerValue(0); IntegerValue new_objective = IntegerValue(0); - // Bounds data is only used by relaxation neighborhoods. - IntegerValue initial_best_objective_bound = IntegerValue(0); - IntegerValue new_objective_bound = IntegerValue(0); - // This is just used to construct a deterministic order for the updates. bool operator<(const SolveData& o) const { return std::tie(status, difficulty, deterministic_limit, deterministic_time, initial_best_objective, - base_objective, new_objective, - initial_best_objective_bound, new_objective_bound, - neighborhood_id) < + base_objective, new_objective) < std::tie(o.status, o.difficulty, o.deterministic_limit, o.deterministic_time, o.initial_best_objective, - o.base_objective, o.new_objective, - o.initial_best_objective_bound, o.new_objective_bound, - o.neighborhood_id); + o.base_objective, o.new_objective); } }; void AddSolveData(SolveData data) { @@ -690,60 +674,6 @@ class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator { SharedIncompleteSolutionManager* incomplete_solutions_; }; -// Generates a relaxation of the original model by removing a consecutive span -// of constraints starting at a random index. The number of constraints removed -// is in sync with the difficulty passed to the generator. -class ConsecutiveConstraintsRelaxationNeighborhoodGenerator - : public NeighborhoodGenerator { - public: - explicit ConsecutiveConstraintsRelaxationNeighborhoodGenerator( - NeighborhoodGeneratorHelper const* helper, const std::string& name) - : NeighborhoodGenerator(name, helper) {} - Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, absl::BitGenRef random) final; - - bool IsRelaxationGenerator() const override { return true; } - bool ReadyToGenerate() const override { return true; } -}; - -// Generates a relaxation of the original model by removing some constraints -// randomly with a given weight for each constraint that controls the -// probability of constraint getting removed. The number of constraints removed -// is in sync with the difficulty passed to the generator. Higher weighted -// constraints are more likely to get removed. -class WeightedRandomRelaxationNeighborhoodGenerator - : public NeighborhoodGenerator { - public: - WeightedRandomRelaxationNeighborhoodGenerator( - NeighborhoodGeneratorHelper const* helper, const std::string& name); - - // Generates the neighborhood as described above. Also stores the removed - // constraints indices for adjusting the weights. - Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, absl::BitGenRef random) final; - - bool IsRelaxationGenerator() const override { return true; } - bool ReadyToGenerate() const override { return true; } - - private: - // Adjusts the weights of the constraints removed to get the neighborhood - // based on the solve_data. - void AdditionalProcessingOnSynchronize(const SolveData& solve_data) override - ABSL_EXCLUSIVE_LOCKS_REQUIRED(generator_mutex_); - - // Higher weighted constraints are more likely to get removed. - std::vector constraint_weights_; - int num_removable_constraints_ = 0; - - // Indices of the removed constraints per generated neighborhood. - absl::flat_hash_map> removed_constraints_ - ABSL_GUARDED_BY(generator_mutex_); - - // TODO(user): Move this to parent class if other generators start using - // feedbacks. - int64_t next_available_id_ ABSL_GUARDED_BY(generator_mutex_) = 0; -}; - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_loader.cc b/ortools/sat/cp_model_loader.cc index 54dc2d2235..d489d34717 100644 --- a/ortools/sat/cp_model_loader.cc +++ b/ortools/sat/cp_model_loader.cc @@ -383,6 +383,7 @@ void ExtractEncoding(const CpModelProto& model_proto, Model* m) { // if some of the literal view used in the LP are created later, but that // should be fixable via calls to implied_bounds->NotifyNewIntegerView(). auto* implied_bounds = m->GetOrCreate(); + auto* detector = m->GetOrCreate(); // Detection of literal equivalent to (i_var >= bound). We also collect // all the half-refied part and we will sort the vector for detection of the @@ -468,6 +469,10 @@ void ExtractEncoding(const CpModelProto& model_proto, Model* m) { { const Domain inter = domain.IntersectionWith(domain_if_enforced); if (!inter.IsEmpty() && inter.Min() == inter.Max()) { + if (inter.Min() == 0) { + detector->ProcessConditionalZero(enforcement_literal, + mapping->Integer(var)); + } var_to_equalities[var].push_back( {&ct, enforcement_literal, inter.Min(), true}); if (domain.Contains(inter.Min())) { @@ -869,6 +874,9 @@ void LoadBoolOrConstraint(const ConstraintProto& ct, Model* m) { literals.push_back(mapping->Literal(ref).Negated()); } m->Add(ClauseConstraint(literals)); + if (literals.size() == 3) { + m->GetOrCreate()->ProcessTernaryClause(literals); + } } void LoadBoolAndConstraint(const ConstraintProto& ct, Model* m) { @@ -894,7 +902,11 @@ void LoadAtMostOneConstraint(const ConstraintProto& ct, Model* m) { void LoadExactlyOneConstraint(const ConstraintProto& ct, Model* m) { auto* mapping = m->GetOrCreate(); CHECK(!HasEnforcementLiteral(ct)) << "Not supported."; - m->Add(ExactlyOneConstraint(mapping->Literals(ct.exactly_one().literals()))); + const auto& literals = mapping->Literals(ct.exactly_one().literals()); + m->Add(ExactlyOneConstraint(literals)); + if (literals.size() == 3) { + m->GetOrCreate()->ProcessTernaryExactlyOne(literals); + } } void LoadBoolXorConstraint(const ConstraintProto& ct, Model* m) { @@ -975,6 +987,18 @@ void LoadEquivalenceNeqAC(const std::vector enforcement_literal, } } +bool IsPartOfProductEncoding(const ConstraintProto& ct) { + if (ct.enforcement_literal().size() != 1) return false; + if (ct.linear().vars().size() > 2) return false; + if (ct.linear().domain().size() != 2) return false; + if (ct.linear().domain(0) != 0) return false; + if (ct.linear().domain(1) != 0) return false; + for (const int64_t coeff : ct.linear().coeffs()) { + if (std::abs(coeff) != 1) return false; + } + return true; +} + } // namespace void LoadLinearConstraint(const ConstraintProto& ct, Model* m) { @@ -995,6 +1019,23 @@ void LoadLinearConstraint(const ConstraintProto& ct, Model* m) { return; } + if (IsPartOfProductEncoding(ct)) { + const Literal l = mapping->Literal(ct.enforcement_literal(0)); + auto* detector = m->GetOrCreate(); + if (ct.linear().vars().size() == 1) { + // TODO(user): Actually this should never be called since we process + // linear1 in ExtractEncoding(). + detector->ProcessConditionalZero(l, + mapping->Integer(ct.linear().vars(0))); + } else if (ct.linear().vars().size() == 2) { + const IntegerVariable x = mapping->Integer(ct.linear().vars(0)); + const IntegerVariable y = mapping->Integer(ct.linear().vars(1)); + detector->ProcessConditionalEquality( + l, x, + ct.linear().coeffs(0) == ct.linear().coeffs(1) ? NegationOf(y) : y); + } + } + auto* integer_trail = m->GetOrCreate(); const std::vector vars = mapping->Integers(ct.linear().vars()); diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index b67da793b1..d0d8b45de2 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -2483,6 +2483,15 @@ bool CpModelPresolver::PresolveDiophantine(ConstraintProto* ct) { // This tries to decompose the constraint into coeff * part1 + part2 and show // that the value that part2 take is not important, thus the constraint can // only be transformed on a constraint on the first part. +// +// TODO(user): Improve !! we miss simple case like x + 47 y + 50 z >= 50 +// for positive variables. We should remove x, and ideally we should rewrite +// this as y + 2z >= 2 if we can show that its relaxation is just better? +// We should at least see that it is the same as 47y + 50 z >= 48. +// +// TODO(user): One easy algo is to first remove all enforcement term (even +// non-Boolean one) before applying the algo here and then re-linearize the +// non-Boolean terms. void CpModelPresolver::TryToReduceCoefficientsOfLinearConstraint( int c, ConstraintProto* ct) { if (ct->constraint_case() != ConstraintProto::kLinear) return; @@ -3282,7 +3291,7 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint( // Note that for now the default is false, and also there are problem calling // GetOrCreateVarValueEncoding() after expansion because we might have removed // the variable used in the encoding. - const bool only_booleans = + const bool only_extract_booleans = !context_->params().presolve_extract_integer_enforcement() || context_->ModelIsExpanded(); @@ -3302,10 +3311,25 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint( const bool is_boolean = context_->CanBeUsedAsLiteral(ref); if (context_->IsFixed(ref) || coeff < threshold || - (only_booleans && !is_boolean)) { - // We keep this term. + (only_extract_booleans && !is_boolean)) { mutable_arg->set_vars(new_size, mutable_arg->vars(i)); - mutable_arg->set_coeffs(new_size, mutable_arg->coeffs(i)); + + // We keep this term but reduces its coeff. + // This is only for the case where only_extract_booleans == true. + if (coeff > threshold) { + context_->UpdateRuleStats("linear: coefficient strenghtening."); + if (lower_bounded) { + // coeff * (X - LB + LB) -> threshold * (X - LB) + coeff * LB + rhs_offset -= (coeff - threshold) * context_->MinOf(ref); + } else { + // coeff * (X - UB + UB) -> threshold * (X - UB) + coeff * UB + rhs_offset -= (coeff - threshold) * context_->MaxOf(ref); + } + mutable_arg->set_coeffs(new_size, + arg.coeffs(i) > 0 ? threshold : -threshold); + } else { + mutable_arg->set_coeffs(new_size, arg.coeffs(i)); + } ++new_size; continue; } diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index be1940170e..15072dd6cf 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -693,8 +693,7 @@ std::vector GetDiverseSetOfParameters( // worker for them. int target = base_params.num_workers(); if (!base_params.interleave_search() && - (base_params.use_rins_lns() || base_params.use_relaxation_lns() || - base_params.use_feasibility_pump())) { + (base_params.use_rins_lns() || base_params.use_feasibility_pump())) { target = std::max(1, base_params.num_workers() - 1); } if (!base_params.interleave_search() && result.size() > target) { diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 8e2c104152..4fe87bdac2 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -67,6 +67,7 @@ #include "ortools/sat/drat_checker.h" #include "ortools/sat/drat_proof_handler.h" #include "ortools/sat/feasibility_pump.h" +#include "ortools/sat/implied_bounds.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_expr.h" #include "ortools/sat/integer_search.h" @@ -598,16 +599,14 @@ void LoadDebugSolution(const CpModelProto& model_proto, Model* model) { const IntegerVariable objective_var = objective_def->objective_var; const int64_t objective_value = - ComputeInnerObjective(model_proto.objective(), response); + ComputeInnerObjective(model_proto.objective(), response.solution()); debug_solution[objective_var] = objective_value; debug_solution[NegationOf(objective_var)] = -objective_value; #endif // __PORTABLE_PLATFORM__ } -void FillSolutionInResponse(const CpModelProto& model_proto, const Model& model, - CpSolverResponse* response) { - response->clear_solution(); - +std::vector GetSolutionValues(const CpModelProto& model_proto, + const Model& model) { auto* mapping = model.Get(); auto* trail = model.Get(); @@ -636,9 +635,7 @@ void FillSolutionInResponse(const CpModelProto& model_proto, const Model& model, // TODO(user): Checks against initial model. CHECK(SolutionIsFeasible(model_proto, solution)); } - for (const int64_t value : solution) { - response->add_solution(value); - } + return solution; } namespace { @@ -1269,6 +1266,9 @@ void LoadBaseModel(const CpModelProto& model_proto, Model* model) { model->GetOrCreate() ->AddAllImplicationsBetweenAssociatedLiterals(); if (!sat_solver->FinishPropagation()) return unsat(); + + model->GetOrCreate()->ProcessImplicationGraph( + model->GetOrCreate()); } void LoadFeasibilityPump(const CpModelProto& model_proto, Model* model) { @@ -1522,10 +1522,9 @@ void LoadCpModel(const CpModelProto& model_proto, Model* model) { const std::string solution_info = model->Name(); const auto solution_observer = [&model_proto, model, solution_info, shared_response_manager]() { - CpSolverResponse response; - FillSolutionInResponse(model_proto, *model, &response); - response.set_solution_info(solution_info); - shared_response_manager->NewSolution(response, model); + const std::vector solution = + GetSolutionValues(model_proto, *model); + shared_response_manager->NewSolution(solution, solution_info, model); }; const auto& objective = *model->GetOrCreate(); @@ -1557,10 +1556,9 @@ void SolveLoadedCpModel(const CpModelProto& model_proto, Model* model) { const std::string& solution_info = model->Name(); const auto solution_observer = [&model_proto, &model, &solution_info, &shared_response_manager]() { - CpSolverResponse response; - FillSolutionInResponse(model_proto, *model, &response); - response.set_solution_info(solution_info); - shared_response_manager->NewSolution(response, model); + const std::vector solution = + GetSolutionValues(model_proto, *model); + shared_response_manager->NewSolution(solution, solution_info, model); }; // Reconfigure search heuristic if it was changed. @@ -1694,10 +1692,10 @@ void QuickSolveWithHint(const CpModelProto& model_proto, Model* model) { const std::string& solution_info = model->Name(); if (status == SatSolver::Status::FEASIBLE) { - CpSolverResponse response; - FillSolutionInResponse(model_proto, *model, &response); - response.set_solution_info(absl::StrCat(solution_info, " [hint]")); - shared_response_manager->NewSolution(response, model); + const std::vector solution = + GetSolutionValues(model_proto, *model); + shared_response_manager->NewSolution( + solution, absl::StrCat(solution_info, " [hint]"), model); if (!model_proto.has_objective()) { if (parameters->enumerate_all_solutions()) { @@ -1828,18 +1826,17 @@ void MinimizeL1DistanceWithHint(const CpModelProto& model_proto, Model* model) { const std::string& solution_info = model->Name(); if (status == SatSolver::Status::FEASIBLE) { - CpSolverResponse response; - FillSolutionInResponse(model_proto, local_model, &response); + const std::vector solution = + GetSolutionValues(model_proto, local_model); if (DEBUG_MODE) { - CpSolverResponse updated_response; - FillSolutionInResponse(updated_model_proto, local_model, - &updated_response); + const std::vector updated_solution = + GetSolutionValues(updated_model_proto, local_model); LOG(INFO) << "Found solution with repaired hint penalty = " << ComputeInnerObjective(updated_model_proto.objective(), - updated_response); + updated_solution); } - response.set_solution_info(absl::StrCat(solution_info, " [repaired]")); - shared_response_manager->NewSolution(response, &local_model); + shared_response_manager->NewSolution( + solution, absl::StrCat(solution_info, " [repaired]"), &local_model); } } @@ -2512,8 +2509,6 @@ class LnsSolver : public SubSolver { if (!neighborhood.is_generated) return; - data.neighborhood_id = neighborhood.id; - const int64_t num_calls = std::max(int64_t{1}, generator_->num_calls()); const double fully_solved_proportion = static_cast(generator_->num_fully_solved_calls()) / @@ -2522,7 +2517,7 @@ class LnsSolver : public SubSolver { if (!neighborhood.source_info.empty()) { absl::StrAppend(&source_info, "_", neighborhood.source_info); } - const std::string solution_info = absl::StrFormat( + const std::string lns_info = absl::StrFormat( "%s(d=%0.2f s=%i t=%0.2f p=%0.2f)", source_info, data.difficulty, task_id, data.deterministic_limit, fully_solved_proportion); @@ -2534,7 +2529,7 @@ class LnsSolver : public SubSolver { local_params.set_symmetry_level(0); local_params.set_solution_pool_size(0); - Model local_model(solution_info); + Model local_model(lns_info); *(local_model.GetOrCreate()) = local_params; TimeLimit* local_time_limit = local_model.GetOrCreate(); local_time_limit->ResetLimitFromParameters(local_params); @@ -2616,106 +2611,49 @@ class LnsSolver : public SubSolver { local_response_manager->SetResponseStatus(presolve_status); } CpSolverResponse local_response = local_response_manager->GetResponse(); + const std::string solution_info = local_response.solution_info(); + std::vector solution_values(local_response.solution().begin(), + local_response.solution().end()); + data.status = local_response.status(); // TODO(user): we actually do not need to postsolve if the solution is // not going to be used... if (local_params.cp_model_presolve() && - (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::FEASIBLE)) { - std::vector solution(local_response.solution().begin(), - local_response.solution().end()); - PostsolveResponseWrapper(local_params, - helper_->ModelProto().variables_size(), - mapping_proto, postsolve_mapping, &solution); - local_response.mutable_solution()->Assign(solution.begin(), - solution.end()); + (data.status == CpSolverStatus::OPTIMAL || + data.status == CpSolverStatus::FEASIBLE)) { + PostsolveResponseWrapper( + local_params, helper_->ModelProto().variables_size(), mapping_proto, + postsolve_mapping, &solution_values); + local_response.mutable_solution()->Assign(solution_values.begin(), + solution_values.end()); } - data.status = local_response.status(); data.deterministic_time = local_time_limit->GetElapsedDeterministicTime(); bool new_solution = false; bool display_lns_info = VLOG_IS_ON(2); - if (generator_->IsRelaxationGenerator()) { - bool has_feasible_solution = false; - if (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::FEASIBLE) { - has_feasible_solution = true; - } + const std::vector solution(local_response.solution().begin(), + local_response.solution().end()); - if (local_response.status() == CpSolverStatus::INFEASIBLE) { - shared_->response->NotifyThatImprovingProblemIsInfeasible( - local_response.solution_info()); - } - - if (shared_->model_proto->has_objective()) { - // TODO(user): This is not deterministic since it is updated without - // synchronization! So we shouldn't base the LNS score out of that. - const IntegerValue current_obj_lb = - shared_->response->GetInnerObjectiveLowerBound(); - - const IntegerValue local_obj_lb = - local_response_manager->GetInnerObjectiveLowerBound(); - - const double scaled_local_obj_bound = ScaleObjectiveValue( - lns_fragment.objective(), local_obj_lb.value()); - - // Update the bound. - const IntegerValue new_inner_obj_lb = IntegerValue( - std::ceil(UnscaleObjectiveValue(shared_->model_proto->objective(), - scaled_local_obj_bound) - - 1e-6)); - data.new_objective_bound = new_inner_obj_lb; - data.initial_best_objective_bound = current_obj_lb; - if (new_inner_obj_lb > current_obj_lb) { - shared_->response->UpdateInnerObjectiveBounds( - solution_info, new_inner_obj_lb, kMaxIntegerValue); - } - } - - // If we have a solution of the relaxed problem, we check if it is also - // a valid solution of the non-relaxed one. - if (has_feasible_solution) { - if (SolutionIsFeasible( - *shared_->model_proto, - std::vector(local_response.solution().begin(), - local_response.solution().end()))) { - shared_->response->NewSolution(local_response, - /*model=*/nullptr); - - // Mark the solution optimal if the relaxation status is optimal. - if (local_response.status() == CpSolverStatus::OPTIMAL) { - shared_->response->NotifyThatImprovingProblemIsInfeasible( - local_response.solution_info()); - shared_->time_limit->Stop(); - } - } - shared_->relaxation_solutions->NewRelaxationSolution(local_response); - } - } else { - const std::vector solution(local_response.solution().begin(), - local_response.solution().end()); - - if (!solution.empty()) { - // A solution that does not pass our validator indicates a bug. We - // abort and dump the problematic model to facilitate debugging. - // - // TODO(user): In a production environment, we should probably just - // ignore this fragment and continue. - const bool feasible = - SolutionIsFeasible(*shared_->model_proto, solution); - if (!feasible) { - 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"); - LOG(INFO) << "Dumping problematic LNS model to '" << name << "'."; - CHECK_OK(file::SetTextProto(name, debug_copy, file::Defaults())); - } - LOG(FATAL) << "Infeasible LNS solution! " << solution_info - << " solved with params " - << local_params.ShortDebugString(); + if (!solution.empty()) { + // A solution that does not pass our validator indicates a bug. We + // abort and dump the problematic model to facilitate debugging. + // + // TODO(user): In a production environment, we should probably just + // ignore this fragment and continue. + const bool feasible = + SolutionIsFeasible(*shared_->model_proto, solution_values); + if (!feasible) { + 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"); + LOG(INFO) << "Dumping problematic LNS model to '" << name << "'."; + CHECK_OK(file::SetTextProto(name, debug_copy, file::Defaults())); } + LOG(FATAL) << "Infeasible LNS solution! " << solution_info + << " solved with params " + << local_params.ShortDebugString(); } // Special case if we solved a part of the full problem! @@ -2729,41 +2667,42 @@ class LnsSolver : public SubSolver { // - It should be fine if all our generator are fully or not at // all included in the variable we are fixing. So we can relax the // test here. Try on z26.mps or spot5_1401.fzn. - if (local_response.status() == CpSolverStatus::OPTIMAL && - !shared_->model_proto->has_symmetry() && !solution.empty() && + if (data.status == CpSolverStatus::OPTIMAL && + !shared_->model_proto->has_symmetry() && !solution_values.empty() && neighborhood.is_simple && !neighborhood.variables_that_can_be_fixed_to_local_optimum .empty()) { display_lns_info = true; shared_->bounds->FixVariablesFromPartialSolution( - solution, + solution_values, neighborhood.variables_that_can_be_fixed_to_local_optimum); } // Finish to fill the SolveData now that the local solve is done. data.new_objective = data.base_objective; - if (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::FEASIBLE) { + if (data.status == CpSolverStatus::OPTIMAL || + data.status == CpSolverStatus::FEASIBLE) { data.new_objective = IntegerValue(ComputeInnerObjective( - shared_->model_proto->objective(), local_response)); + shared_->model_proto->objective(), solution_values)); } // Report any feasible solution we have. Optimization: We don't do that // if we just recovered the base solution. - if (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::FEASIBLE) { + if (data.status == CpSolverStatus::OPTIMAL || + data.status == CpSolverStatus::FEASIBLE) { const std::vector base_solution( base_response.solution().begin(), base_response.solution().end()); - if (solution != base_solution) { + if (solution_values != base_solution) { new_solution = true; - shared_->response->NewSolution(local_response, /*model=*/nullptr); + shared_->response->NewSolution(solution_values, solution_info, + /*model=*/nullptr); } } if (!neighborhood.is_reduced && - (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::INFEASIBLE)) { + (data.status == CpSolverStatus::OPTIMAL || + data.status == CpSolverStatus::INFEASIBLE)) { shared_->response->NotifyThatImprovingProblemIsInfeasible( - local_response.solution_info()); + solution_info); shared_->time_limit->Stop(); } } @@ -2777,11 +2716,11 @@ class LnsSolver : public SubSolver { const double base_obj = ScaleObjectiveValue( shared_->model_proto->objective(), ComputeInnerObjective(shared_->model_proto->objective(), - base_response)); + base_response.solution())); const double new_obj = ScaleObjectiveValue( shared_->model_proto->objective(), ComputeInnerObjective(shared_->model_proto->objective(), - local_response)); + solution_values)); absl::StrAppend(&s, " [new_sol:", base_obj, " -> ", new_obj, "]"); } if (neighborhood.is_simple) { @@ -2831,13 +2770,6 @@ void SolveCpModelParallel(const CpModelProto& model_proto, std::unique_ptr shared_relaxation_solutions; - if (params.use_relaxation_lns()) { - shared_relaxation_solutions = - std::make_unique( - /*num_solutions_to_keep=*/10); - global_model->Register( - shared_relaxation_solutions.get()); - } auto shared_lp_solutions = std::make_unique( /*num_solutions_to_keep=*/10); @@ -3041,19 +2973,6 @@ void SolveCpModelParallel(const CpModelProto& model_proto, absl::StrCat("rens_lns_", local_params.name())), local_params, helper, &shared)); } - - if (params.use_relaxation_lns()) { - subsolvers.push_back(std::make_unique( - std::make_unique< - ConsecutiveConstraintsRelaxationNeighborhoodGenerator>( - helper, absl::StrCat("rnd_rel_lns_", local_params.name())), - local_params, helper, &shared)); - - subsolvers.push_back(std::make_unique( - std::make_unique( - helper, absl::StrCat("wgt_rel_lns_", local_params.name())), - local_params, helper, &shared)); - } } // Add a synchronization point for the gap integral that is executed last. diff --git a/ortools/sat/cp_model_utils.cc b/ortools/sat/cp_model_utils.cc index 94bceeed4c..4d09ec82f7 100644 --- a/ortools/sat/cp_model_utils.cc +++ b/ortools/sat/cp_model_utils.cc @@ -517,14 +517,14 @@ std::vector UsedIntervals(const ConstraintProto& ct) { } int64_t ComputeInnerObjective(const CpObjectiveProto& objective, - const CpSolverResponse& response) { + absl::Span solution) { int64_t objective_value = 0; for (int i = 0; i < objective.vars_size(); ++i) { int64_t coeff = objective.coeffs(i); const int ref = objective.vars(i); const int var = PositiveRef(ref); if (!RefIsPositive(ref)) coeff = -coeff; - objective_value += coeff * response.solution()[var]; + objective_value += coeff * solution[var]; } return objective_value; } diff --git a/ortools/sat/cp_model_utils.h b/ortools/sat/cp_model_utils.h index c0f791a02a..a53e8919da 100644 --- a/ortools/sat/cp_model_utils.h +++ b/ortools/sat/cp_model_utils.h @@ -164,7 +164,7 @@ inline double UnscaleObjectiveValue(const CpObjectiveProto& proto, // This is the objective without offset and scaling. Call ScaleObjectiveValue() // to get the user facing objective. int64_t ComputeInnerObjective(const CpObjectiveProto& objective, - const CpSolverResponse& response); + absl::Span solution); // Returns true if a linear expression can be reduced to a single ref. bool ExpressionContainsSingleRef(const LinearExpressionProto& expr); diff --git a/ortools/sat/implied_bounds.cc b/ortools/sat/implied_bounds.cc index e4cd937deb..7094b1233e 100644 --- a/ortools/sat/implied_bounds.cc +++ b/ortools/sat/implied_bounds.cc @@ -16,6 +16,8 @@ #include #include +#include +#include #include #include #include @@ -41,10 +43,14 @@ namespace sat { // Just display some global statistics on destruction. ImpliedBounds::~ImpliedBounds() { - VLOG(1) << num_deductions_ << " enqueued deductions."; - VLOG(1) << bounds_.size() << " implied bounds stored."; - VLOG(1) << num_enqueued_in_var_to_bounds_ - << " implied bounds with view stored."; + if (!VLOG_IS_ON(1)) return; + if (shared_stats_ == nullptr) return; + std::vector> stats; + stats.push_back({"implied_bound/num_deductions", num_deductions_}); + stats.push_back({"implied_bound/num_stored", bounds_.size()}); + stats.push_back( + {"implied_bound/num_stored_with_view", num_enqueued_in_var_to_bounds_}); + shared_stats_->AddStats(stats); } void ImpliedBounds::Add(Literal literal, IntegerLiteral integer_literal) { @@ -124,7 +130,7 @@ void ImpliedBounds::Add(Literal literal, IntegerLiteral integer_literal) { level_zero_lower_bounds_[var] = deduction; new_level_zero_bounds_.Set(var); - VLOG(1) << "Deduction old: " + VLOG(2) << "Deduction old: " << IntegerLiteral::GreaterOrEqual( var, integer_trail_->LevelZeroLowerBound(var)) << " new: " << IntegerLiteral::GreaterOrEqual(var, deduction); @@ -339,7 +345,7 @@ std::vector TryToReconcileSize2Encodings( const std::vector& right_enc = integer_encoder->FullDomainEncoding(right.var); if (left_enc.size() != 2 || right_enc.size() != 2) { - VLOG(3) << "encodings are not fully propagated"; + VLOG(2) << "encodings are not fully propagated"; return terms; } @@ -426,7 +432,7 @@ std::vector TryToDecomposeProduct( } if (compatible_keys.size() > 1) { - VLOG(1) << "More than one exactly_one involved in the encoding of the two " + VLOG(3) << "More than one exactly_one involved in the encoding of the two " "variables"; } @@ -515,5 +521,248 @@ bool DetectLinearEncodingOfProducts(const AffineExpression& left, return true; } +ProductDetector::ProductDetector(Model* model) + : enabled_(model->GetOrCreate()->linearization_level() > 1), + sat_solver_(model->GetOrCreate()), + trail_(model->GetOrCreate()), + integer_trail_(model->GetOrCreate()), + integer_encoder_(model->GetOrCreate()), + shared_stats_(model->GetOrCreate()) {} + +ProductDetector::~ProductDetector() { + if (!VLOG_IS_ON(1)) return; + if (shared_stats_ == nullptr) return; + std::vector> stats; + stats.push_back( + {"product_detector/num_processed_binary", num_processed_binary_}); + stats.push_back( + {"product_detector/num_processed_exactly_one", num_processed_exo_}); + stats.push_back( + {"product_detector/num_processed_ternary", num_processed_ternary_}); + stats.push_back({"product_detector/num_trail_updates", num_trail_updates_}); + stats.push_back({"product_detector/num_products", num_products_}); + stats.push_back({"product_detector/num_conditional_equalities", + num_conditional_equalities_}); + stats.push_back( + {"product_detector/num_conditional_zeros", num_conditional_zeros_}); + stats.push_back({"product_detector/num_int_products", num_int_products_}); + shared_stats_->AddStats(stats); +} + +void ProductDetector::ProcessTernaryClause( + absl::Span ternary_clause) { + if (!enabled_) return; + if (ternary_clause.size() != 3) return; + ++num_processed_ternary_; + candidates_[GetKey(ternary_clause[0].Index(), ternary_clause[1].Index())] + .push_back(ternary_clause[2].Index()); + candidates_[GetKey(ternary_clause[0].Index(), ternary_clause[2].Index())] + .push_back(ternary_clause[1].Index()); + candidates_[GetKey(ternary_clause[1].Index(), ternary_clause[2].Index())] + .push_back(ternary_clause[0].Index()); + + // We mark the literal of the ternary clause as seen. + // Only a => b with a seen need to be looked at. + for (const Literal l : ternary_clause) { + if (l.Index() >= seen_.size()) seen_.resize(l.Index() + 1); + seen_[l.Index()] = true; + } +} + +void ProductDetector::ProcessTernaryExactlyOne( + absl::Span ternary_exo) { + if (!enabled_) return; + if (ternary_exo.size() != 3) return; + ++num_processed_exo_; + ProcessNewProduct(ternary_exo[0].Index(), ternary_exo[1].NegatedIndex(), + ternary_exo[2].NegatedIndex()); + ProcessNewProduct(ternary_exo[1].Index(), ternary_exo[0].NegatedIndex(), + ternary_exo[2].NegatedIndex()); + ProcessNewProduct(ternary_exo[2].Index(), ternary_exo[0].NegatedIndex(), + ternary_exo[1].NegatedIndex()); +} + +// TODO(user): As product are discovered, we could remove entries from our +// hash maps! +void ProductDetector::ProcessBinaryClause( + absl::Span binary_clause) { + if (!enabled_) return; + if (binary_clause.size() != 2) return; + ++num_processed_binary_; + const std::array key = + GetKey(binary_clause[0].NegatedIndex(), binary_clause[1].NegatedIndex()); + std::array ternary; + for (const LiteralIndex l : candidates_[key]) { + ternary[0] = key[0]; + ternary[1] = key[1]; + ternary[2] = l; + std::sort(ternary.begin(), ternary.end()); + const int l_index = ternary[0] == l ? 0 : ternary[1] == l ? 1 : 2; + std::bitset<3>& bs = detector_[ternary]; + if (bs[l_index]) continue; + bs[l_index] = true; + if (bs[0] && bs[1] && l_index != 2) { + ProcessNewProduct(ternary[2], Literal(ternary[0]).NegatedIndex(), + Literal(ternary[1]).NegatedIndex()); + } + if (bs[0] && bs[2] && l_index != 1) { + ProcessNewProduct(ternary[1], Literal(ternary[0]).NegatedIndex(), + Literal(ternary[2]).NegatedIndex()); + } + if (bs[1] && bs[2] && l_index != 0) { + ProcessNewProduct(ternary[0], Literal(ternary[1]).NegatedIndex(), + Literal(ternary[2]).NegatedIndex()); + } + } +} + +void ProductDetector::ProcessImplicationGraph(BinaryImplicationGraph* graph) { + if (!enabled_) return; + for (LiteralIndex a(0); a < seen_.size(); ++a) { + if (!seen_[a]) continue; + if (trail_->Assignment().LiteralIsAssigned(Literal(a))) continue; + const Literal not_a = Literal(a).Negated(); + for (const Literal b : graph->DirectImplications(Literal(a))) { + ProcessBinaryClause({not_a, b}); // a => b; + } + } +} + +void ProductDetector::ProcessTrailAtLevelOne() { + if (!enabled_) return; + if (trail_->CurrentDecisionLevel() != 1) return; + ++num_trail_updates_; + + const SatSolver::Decision decision = sat_solver_->Decisions()[0]; + if (decision.literal.Index() >= seen_.size() || + !seen_[decision.literal.Index()]) { + return; + } + const Literal not_a = decision.literal.Negated(); + const int current_index = trail_->Index(); + for (int i = decision.trail_index + 1; i < current_index; ++i) { + const Literal b = (*trail_)[i]; + ProcessBinaryClause({not_a, b}); + } +} + +LiteralIndex ProductDetector::GetProduct(Literal a, Literal b) const { + const auto it = products_.find(GetKey(a.Index(), b.Index())); + if (it == products_.end()) return kNoLiteralIndex; + return it->second; +} + +std::array ProductDetector::GetKey(LiteralIndex a, + LiteralIndex b) const { + std::array key{a, b}; + if (key[0] > key[1]) std::swap(key[0], key[1]); + return key; +} + +void ProductDetector::ProcessNewProduct(LiteralIndex p, LiteralIndex a, + LiteralIndex b) { + // If many literal correspond to the same product, we just keep one. + ++num_products_; + products_[GetKey(a, b)] = p; + + // This is used by ProductIsLinearizable(). + has_product_.insert( + GetKey(Literal(a).IsPositive() ? a : Literal(a).NegatedIndex(), + Literal(b).IsPositive() ? b : Literal(b).NegatedIndex())); +} + +bool ProductDetector::ProductIsLinearizable(IntegerVariable a, + IntegerVariable b) const { + if (a == b) return true; + if (a == NegationOf(b)) return true; + + // Otherwise, we need both a and b to be expressible as linear expression + // involving Booleans whose product is also expressible. + if (integer_trail_->InitialVariableDomain(a).Size() != 2) return false; + if (integer_trail_->InitialVariableDomain(b).Size() != 2) return false; + + const LiteralIndex la = + integer_encoder_->GetAssociatedLiteral(IntegerLiteral::GreaterOrEqual( + a, integer_trail_->LevelZeroUpperBound(a))); + if (la == kNoLiteralIndex) return false; + + const LiteralIndex lb = + integer_encoder_->GetAssociatedLiteral(IntegerLiteral::GreaterOrEqual( + b, integer_trail_->LevelZeroUpperBound(b))); + if (lb == kNoLiteralIndex) return false; + + // Any product involving la/not(la) * lb/not(lb) can be used. + return has_product_.contains( + GetKey(Literal(la).IsPositive() ? la : Literal(la).NegatedIndex(), + Literal(lb).IsPositive() ? lb : Literal(lb).NegatedIndex())); +} + +IntegerVariable ProductDetector::GetProduct(Literal a, + IntegerVariable b) const { + const auto it = int_products_.find({a.Index(), PositiveVariable(b)}); + if (it == int_products_.end()) return kNoIntegerVariable; + return VariableIsPositive(b) ? it->second : NegationOf(it->second); +} + +void ProductDetector::ProcessNewProduct(IntegerVariable p, Literal l, + IntegerVariable x) { + if (!VariableIsPositive(x)) { + x = NegationOf(x); + p = NegationOf(p); + } + + // We only store one product if there are many. + ++num_int_products_; + int_products_[{l.Index(), x}] = p; +} + +void ProductDetector::ProcessConditionalEquality(Literal l, IntegerVariable x, + IntegerVariable y) { + ++num_conditional_equalities_; + if (x == y) return; + + // We process both possibilities (product = x or product = y). + for (int i = 0; i < 2; ++i) { + if (!VariableIsPositive(x)) { + x = NegationOf(x); + y = NegationOf(y); + } + bool seen = false; + + // TODO(user): Linear scan can be bad if b => X = many other variables. + // Hopefully this will not be common. + std::vector& others = + conditional_equalities_[{l.Index(), x}]; + for (const IntegerVariable o : others) { + if (o == y) { + seen = true; + break; + } + } + + if (!seen) { + others.push_back(y); + if (conditional_zeros_.contains({l.NegatedIndex(), x})) { + ProcessNewProduct(/*p=*/x, l, y); + } + } + std::swap(x, y); + } +} + +void ProductDetector::ProcessConditionalZero(Literal l, IntegerVariable p) { + ++num_conditional_zeros_; + p = PositiveVariable(p); + auto [_, inserted] = conditional_zeros_.insert({l.Index(), p}); + if (inserted) { + const auto it = conditional_equalities_.find({l.NegatedIndex(), p}); + if (it != conditional_equalities_.end()) { + for (const IntegerVariable x : it->second) { + ProcessNewProduct(p, l.Negated(), x); + } + } + } +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/implied_bounds.h b/ortools/sat/implied_bounds.h index 4f31c3b905..4dc2d08baf 100644 --- a/ortools/sat/implied_bounds.h +++ b/ortools/sat/implied_bounds.h @@ -15,6 +15,8 @@ #define OR_TOOLS_SAT_IMPLIED_BOUNDS_H_ #include +#include +#include #include #include #include @@ -29,6 +31,7 @@ #include "ortools/sat/sat_base.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/sat_solver.h" +#include "ortools/sat/synchronization.h" #include "ortools/util/bitset.h" #include "ortools/util/strong_integers.h" @@ -85,7 +88,8 @@ class ImpliedBounds { : parameters_(*model->GetOrCreate()), sat_solver_(model->GetOrCreate()), integer_trail_(model->GetOrCreate()), - integer_encoder_(model->GetOrCreate()) {} + integer_encoder_(model->GetOrCreate()), + shared_stats_(model->GetOrCreate()) {} ~ImpliedBounds(); // Adds literal => integer_literal to the repository. @@ -157,6 +161,7 @@ class ImpliedBounds { SatSolver* sat_solver_; IntegerTrail* integer_trail_; IntegerEncoder* integer_encoder_; + SharedStatistics* shared_stats_; // TODO(user): Remove the need for this. std::vector tmp_integer_literals_; @@ -225,6 +230,116 @@ bool DetectLinearEncodingOfProducts(const AffineExpression& left, const AffineExpression& right, Model* model, LinearConstraintBuilder* builder); +// Class used to detect and hold all the information about a variable beeing the +// product of two others. This class is meant to be used by LP relaxation and +// cuts. +class ProductDetector { + public: + explicit ProductDetector(Model* model); + ~ProductDetector(); + + // Internally, a Boolean product is encoded in a linear fashion: + // p = a * b become: + // 1/ a and b => p, i.e. a clause (not(a), not(b), p). + // 2/ p => a and p => b, which is a clause (not(p), a) and (not(p), b). + // + // In particular if we have a+b+c==1 then we have a=b*c, b=a*c, and c=a*b !! + // + // For the detection to work, we must load all ternary clause first, then the + // implication. + void ProcessTernaryClause(absl::Span ternary_clause); + void ProcessTernaryExactlyOne(absl::Span ternary_exo); + void ProcessBinaryClause(absl::Span binary_clause); + + // Utility function to process a bunch of implication all at once. + void ProcessImplicationGraph(BinaryImplicationGraph* graph); + void ProcessTrailAtLevelOne(); + + // We also detect product of Boolean with IntegerVariable. + // After presolve, a product P = l * X should be encoded with: + // l => P = X + // not(l) => P = 0 + // + // TODO(user): Generalize to a * X + b = l * (Y + c) since these are also + // easy to linearize if we see l * Y. + void ProcessConditionalEquality(Literal l, IntegerVariable x, + IntegerVariable y); + void ProcessConditionalZero(Literal l, IntegerVariable p); + + // Query function mainly used for testing. + LiteralIndex GetProduct(Literal a, Literal b) const; + IntegerVariable GetProduct(Literal a, IntegerVariable b) const; + + // Query Functions. LinearizeProduct() should only be called if + // ProductIsLinearizable() is true. + bool ProductIsLinearizable(IntegerVariable a, IntegerVariable b) const; + + // TODO(user): Implement! + LinearExpression LinearizeProduct(IntegerVariable a, IntegerVariable b); + + // Returns an expression that is always lower or equal to the product a * b. + // This use the exact LinearizeProduct() if ProductIsLinearizable() otherwise + // it uses the simple McCormick lower bound. + // + // TODO(user): Implement! + LinearExpression ProductLowerBound(IntegerVariable a, IntegerVariable b); + + private: + std::array GetKey(LiteralIndex a, LiteralIndex b) const; + void ProcessNewProduct(LiteralIndex p, LiteralIndex a, LiteralIndex b); + void ProcessNewProduct(IntegerVariable p, Literal l, IntegerVariable x); + + // Fixed at creation time. + bool enabled_; + SatSolver* sat_solver_; + Trail* trail_; + IntegerTrail* integer_trail_; + IntegerEncoder* integer_encoder_; + SharedStatistics* shared_stats_; + + // No need to process implication a => b if a was never seen. + absl::StrongVector seen_; + + // For each clause of size 3 (l0, l1, l2) and a permutation of index (i, j, k) + // we bitset[i] to true if lj => not(lk) and lk => not(lj). + // + // The key is sorted. + absl::flat_hash_map, std::bitset<3>> detector_; + + // For each (l0, l1) we list all the l2 such that (l0, l1, l2) is a 3 clause. + absl::flat_hash_map, std::vector> + candidates_; + + // Products (a, b) -> p such that p == a * b. They key is sorted. + absl::flat_hash_map, LiteralIndex> products_; + + // Same keys has in products_ but canonicalized so we capture all 4 products + // a * b, (1 - a) * b, a * (1 - b) and (1 - a) * (1 - b) with one query. + absl::flat_hash_set> has_product_; + + // For bool * int detection. Note that we only use positive IntegerVariable + // in the key part. + absl::flat_hash_set> + conditional_zeros_; + absl::flat_hash_map, + std::vector> + conditional_equalities_; + + // Stores l * X = P. + absl::flat_hash_map, IntegerVariable> + int_products_; + + // Stats. + int64_t num_products_ = 0; + int64_t num_int_products_ = 0; + int64_t num_trail_updates_ = 0; + int64_t num_processed_binary_ = 0; + int64_t num_processed_ternary_ = 0; + int64_t num_processed_exo_ = 0; + int64_t num_conditional_zeros_ = 0; + int64_t num_conditional_equalities_ = 0; +}; + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index d17cfef0a2..39d2c27517 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -791,6 +791,7 @@ IntegerSearchHelper::IntegerSearchHelper(Model* model) integer_trail_(model->GetOrCreate()), encoder_(model->GetOrCreate()), implied_bounds_(model->GetOrCreate()), + product_detector_(model->GetOrCreate()), time_limit_(model->GetOrCreate()), pseudo_costs_(model->GetOrCreate()) { // This is needed for recording the pseudo-costs. @@ -907,6 +908,7 @@ bool IntegerSearchHelper::TakeDecision(Literal decision) { // This is "almost free", so we might as well do it. if (old_level == 0 && sat_solver_->CurrentDecisionLevel() == 1) { implied_bounds_->ProcessIntegerTrail(decision); + product_detector_->ProcessTrailAtLevelOne(); } // Update the pseudo costs. diff --git a/ortools/sat/integer_search.h b/ortools/sat/integer_search.h index bc0d9dee7b..73611dfc19 100644 --- a/ortools/sat/integer_search.h +++ b/ortools/sat/integer_search.h @@ -273,6 +273,7 @@ class IntegerSearchHelper { IntegerTrail* integer_trail_; IntegerEncoder* encoder_; ImpliedBounds* implied_bounds_; + ProductDetector* product_detector_; TimeLimit* time_limit_; PseudoCosts* pseudo_costs_; IntegerVariable objective_var_ = kNoIntegerVariable; diff --git a/ortools/sat/linear_propagation.h b/ortools/sat/linear_propagation.h index fa821d01ff..2fccdce83e 100644 --- a/ortools/sat/linear_propagation.h +++ b/ortools/sat/linear_propagation.h @@ -195,10 +195,11 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface { int rev_size; // The size of the non-fixed terms. IntegerValue rev_rhs; // The current rhs, updated on fixed terms. }; - #if !defined(_MSC_VER) + +#if !defined(_MSC_VER) static_assert(sizeof(ConstraintInfo) == 24, "ERROR_ConstraintInfo_is_not_well_compacted"); - #endif +#endif // !defined(_MSC_VER) absl::Span GetCoeffs(const ConstraintInfo& info); absl::Span GetVariables(const ConstraintInfo& info); diff --git a/ortools/sat/probing.cc b/ortools/sat/probing.cc index 57ec6c0493..661a9f2116 100644 --- a/ortools/sat/probing.cc +++ b/ortools/sat/probing.cc @@ -45,6 +45,7 @@ Prober::Prober(Model* model) assignment_(model->GetOrCreate()->Assignment()), integer_trail_(model->GetOrCreate()), implied_bounds_(model->GetOrCreate()), + product_detector_(model->GetOrCreate()), sat_solver_(model->GetOrCreate()), time_limit_(model->GetOrCreate()), implication_graph_(model->GetOrCreate()), @@ -83,6 +84,7 @@ bool Prober::ProbeOneVariableInternal(BooleanVariable b) { } implied_bounds_->ProcessIntegerTrail(decision); + product_detector_->ProcessTrailAtLevelOne(); integer_trail_->AppendNewBounds(&new_integer_bounds_); for (int i = saved_index + 1; i < trail_.Index(); ++i) { const Literal l = trail_[i]; diff --git a/ortools/sat/probing.h b/ortools/sat/probing.h index 64aa880230..19bef6bb0c 100644 --- a/ortools/sat/probing.h +++ b/ortools/sat/probing.h @@ -102,6 +102,7 @@ class Prober { const VariablesAssignment& assignment_; IntegerTrail* integer_trail_; ImpliedBounds* implied_bounds_; + ProductDetector* product_detector_; SatSolver* sat_solver_; TimeLimit* time_limit_; BinaryImplicationGraph* implication_graph_; diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index eaae897046..62a505a20f 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -1119,10 +1119,6 @@ message SatParameters { } optional FPRoundingMethod fp_rounding = 165 [default = PROPAGATION_ASSISTED]; - // Turns on a lns worker which solves relaxed version of the original problem - // by removing constraints from the problem in order to get better bounds. - optional bool use_relaxation_lns = 150 [default = false]; - // If true, registers more lns subsolvers with different parameters. optional bool diversify_lns_params = 137 [default = false]; diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index 74c0ee9bd1..c310cbf60f 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -64,21 +64,22 @@ namespace operations_research { namespace sat { void SharedRelaxationSolutionRepository::NewRelaxationSolution( - const CpSolverResponse& response) { + absl::Span solution_values, + IntegerValue inner_objective_value) { // Note that the Add() method already applies mutex lock. So we don't need it // here. - if (response.solution().empty()) return; + if (solution_values.empty()) return; // Add this solution to the pool. SharedSolutionRepository::Solution solution; - solution.variable_values.assign(response.solution().begin(), - response.solution().end()); + solution.variable_values.assign(solution_values.begin(), + solution_values.end()); // For now we use the negated lower bound as the "internal objective" to // prefer solution with an higher bound. // // Note: If the model doesn't have objective, the best_objective_bound is set // to default value 0. - solution.rank = -response.best_objective_bound(); + solution.rank = -inner_objective_value.value(); Add(solution); } @@ -507,28 +508,29 @@ void SharedResponseManager::FillObjectiveValuesInBestResponse() { best_response_.set_gap_integral(gap_integral_); } -void SharedResponseManager::NewSolution(const CpSolverResponse& response, - Model* model) { +void SharedResponseManager::NewSolution( + absl::Span solution_values, const std::string& solution_info, + Model* model) { absl::MutexLock mutex_lock(&mutex_); // Special case if the user asked to keep solutions in the pool. if (objective_or_null_ == nullptr && parameters_.enumerate_all_solutions() && parameters_.fill_additional_solutions_in_response()) { SharedSolutionRepository::Solution solution; - solution.variable_values.assign(response.solution().begin(), - response.solution().end()); + solution.variable_values.assign(solution_values.begin(), + solution_values.end()); solutions_.Add(solution); } if (objective_or_null_ != nullptr) { const int64_t objective_value = - ComputeInnerObjective(*objective_or_null_, response); + ComputeInnerObjective(*objective_or_null_, solution_values); // Add this solution to the pool, even if it is not improving. - if (!response.solution().empty()) { + if (!solution_values.empty()) { SharedSolutionRepository::Solution solution; - solution.variable_values.assign(response.solution().begin(), - response.solution().end()); + solution.variable_values.assign(solution_values.begin(), + solution_values.end()); solution.rank = objective_value; solutions_.Add(solution); } @@ -563,8 +565,9 @@ void SharedResponseManager::NewSolution(const CpSolverResponse& response, best_response_.set_status(CpSolverStatus::FEASIBLE); } - best_response_.set_solution_info(response.solution_info()); - *best_response_.mutable_solution() = response.solution(); + best_response_.set_solution_info(solution_info); + best_response_.mutable_solution()->Assign(solution_values.begin(), + solution_values.end()); // Mark model as OPTIMAL if the inner bound crossed. if (objective_or_null_ != nullptr && @@ -574,12 +577,13 @@ void SharedResponseManager::NewSolution(const CpSolverResponse& response, // Logging. ++num_solutions_; + // TODO(user): Remove this code and the need for model in this function. if (logger_->LoggingIsEnabled()) { - std::string solution_info = response.solution_info(); + std::string solution_message = solution_info; if (model != nullptr) { const int64_t num_bool = model->Get()->NumVariables(); const int64_t num_fixed = model->Get()->NumFixedVariables(); - absl::StrAppend(&solution_info, " fixed_bools:", num_fixed, "/", + absl::StrAppend(&solution_message, " fixed_bools:", num_fixed, "/", num_bool); } @@ -592,13 +596,14 @@ void SharedResponseManager::NewSolution(const CpSolverResponse& response, if (obj.scaling_factor() < 0) { std::swap(lb, ub); } - RegisterSolutionFound(solution_info); + RegisterSolutionFound(solution_message); SOLVER_LOG(logger_, ProgressMessage(absl::StrCat(num_solutions_), wall_timer_.Get(), best, lb, ub, - solution_info)); + solution_message)); } else { - SOLVER_LOG(logger_, SatProgressMessage(absl::StrCat(num_solutions_), - wall_timer_.Get(), solution_info)); + SOLVER_LOG(logger_, + SatProgressMessage(absl::StrCat(num_solutions_), + wall_timer_.Get(), solution_message)); } } diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 287d4a74e3..69420c5426 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -138,7 +138,8 @@ class SharedRelaxationSolutionRepository explicit SharedRelaxationSolutionRepository(int num_solutions_to_keep) : SharedSolutionRepository(num_solutions_to_keep) {} - void NewRelaxationSolution(const CpSolverResponse& response); + void NewRelaxationSolution(absl::Span solution_values, + IntegerValue inner_objective_value); }; class SharedLPSolutionRepository : public SharedSolutionRepository { @@ -291,12 +292,8 @@ class SharedResponseManager { // Reads the new solution from the response and update our state. For an // optimization problem, we only do something if the solution is strictly // improving. - // - // TODO(user): Only the following fields from response are accessed here, we - // might want a tighter API: - // - solution_info - // - solution - void NewSolution(const CpSolverResponse& response, Model* model); + void NewSolution(absl::Span solution_values, + const std::string& solution_info, Model* model = nullptr); // Changes the solution to reflect the fact that the "improving" problem is // infeasible. This means that if we have a solution, we have proven