diff --git a/ortools/sat/cp_model.proto b/ortools/sat/cp_model.proto index 77433fd288..d125d5bfd0 100644 --- a/ortools/sat/cp_model.proto +++ b/ortools/sat/cp_model.proto @@ -521,6 +521,7 @@ enum CpSolverStatus { // // TODO(user): support returning multiple solutions. Look at the Stubby // streaming API as we probably wants to get them as they are found. +// Next id: 23 message CpSolverResponse { // The status of the solve. CpSolverStatus status = 1; @@ -582,6 +583,7 @@ message CpSolverResponse { double wall_time = 15; double user_time = 16; double deterministic_time = 17; + double primal_integral = 22; // Additional information about how the solution was found. string solution_info = 20; diff --git a/ortools/sat/cp_model_loader.cc b/ortools/sat/cp_model_loader.cc index 0dd9b92fef..d25f2d0647 100644 --- a/ortools/sat/cp_model_loader.cc +++ b/ortools/sat/cp_model_loader.cc @@ -56,6 +56,77 @@ std::vector ValuesFromProto(const Values& values) { return std::vector(values.begin(), values.end()); } +bool IsTooLargeToEncode(int var, CpModelMapping* mapping, + IntegerTrail* integer_trail) { + const IntegerVariable sat_var = mapping->Integer(var); + const IntegerValue lb = integer_trail->LowerBound(sat_var); + const IntegerValue ub = integer_trail->UpperBound(sat_var); + return ub - lb > 1024; // Arbitrary limit value. +} + +bool IsSmallEnoughToAlwaysEncode(int var, CpModelMapping* mapping, + IntegerTrail* integer_trail) { + const IntegerVariable sat_var = mapping->Integer(var); + const IntegerValue lb = integer_trail->LowerBound(sat_var); + const IntegerValue ub = integer_trail->UpperBound(sat_var); + return ub - lb <= 64; // Arbitrary limit value. +} + +bool IsEqCst(const LinearConstraintProto& proto) { + if (proto.domain_size() == 2 && proto.domain(0) == proto.domain(1)) { + return true; + } + + return false; +} + +bool IsNeqCst(const LinearConstraintProto& proto, CpModelMapping* mapping, + IntegerTrail* integer_trail, int64* single_value) { + int64 sum_min = 0; + int64 sum_max = 0; + + for (int i = 0; i < proto.vars_size(); ++i) { + const int var = proto.vars(i); + const int64 coeff = proto.coeffs(i); + const IntegerVariable sat_var = mapping->Integer(var); + const int64 lb = integer_trail->LowerBound(sat_var).value(); + const int64 ub = integer_trail->UpperBound(sat_var).value(); + if (coeff >= 0) { + sum_min += coeff * lb; + sum_max += coeff * ub; + } else { + sum_min += coeff * ub; + sum_max += coeff * lb; + } + } + + if (proto.domain_size() == 4 && proto.domain(0) <= sum_min && + proto.domain(1) + 2 == proto.domain(2) && proto.domain(3) >= sum_max) { + if (single_value != nullptr) { + *single_value = proto.domain(1) + 1; + } + return true; + } + + if (proto.domain_size() == 2 && proto.domain(0) <= sum_min && + proto.domain(1) == sum_max - 1) { + if (single_value != nullptr) { + *single_value = sum_max; + } + return true; + } + + if (proto.domain_size() == 2 && proto.domain(0) == sum_min + 1 && + proto.domain(1) >= sum_max) { + if (single_value != nullptr) { + *single_value = sum_min; + } + return true; + } + + return false; +} + } // namespace void CpModelMapping::CreateVariables(const CpModelProto& model_proto, @@ -704,53 +775,47 @@ bool FullEncodingFixedPointComputer::ProcessLinear(ConstraintIndex ct_index) { const ConstraintProto& ct = model_proto_.constraints(ct_index.value()); if (parameters_.boolean_encoding_level() == 0) return true; - // Only act when the constraint is an equality. - if (ct.linear().domain(0) != ct.linear().domain(1)) return true; + // Only act when the constraint is an equality, or non-equality. + if (!IsEqCst(ct.linear()) && + !IsNeqCst(ct.linear(), mapping_, integer_trail_, nullptr)) { + return true; + } - // If some domain is too large, abort; + // If some domain is too large, abort. if (!constraint_is_registered_[ct_index]) { for (const int v : ct.linear().vars()) { - const IntegerVariable var = mapping_->Integer(v); - const IntegerValue lb = integer_trail_->LowerBound(var); - const IntegerValue ub = integer_trail_->UpperBound(var); - if (ub - lb > 1024) return true; // Arbitrary limit value. + if (IsTooLargeToEncode(v, mapping_, integer_trail_)) { + return true; + } } } const int num_vars = ct.linear().vars_size(); - if (HasEnforcementLiteral(ct)) { + if (HasEnforcementLiteral(ct) && num_vars == 1) { // Fully encode x in half-reified equality b => x == constant. - if (num_vars == 1) FullyEncode(ct.linear().vars(0)); - - // Skip any other form of half-reified linear constraint for now. + FullyEncode(ct.linear().vars(0)); + return true; + } else if (num_vars == 2) { + // Fully encode x and y in [b =>] ax + by == cst or [b =>] ax + by != cst. + for (const int v : ct.linear().vars()) { + if (!IsSmallEnoughToAlwaysEncode(v, mapping_, integer_trail_) && + !IsFullyEncoded(v)) { + // Register on remaining variables if not already done. + if (!constraint_is_registered_[ct_index]) { + for (const int var : ct.linear().vars()) { + if (!IsFullyEncoded(var)) Register(ct_index, var); + } + } + // Can continue looking at the constraint. + return false; + } + } + FullyEncode(ct.linear().vars(0)); + FullyEncode(ct.linear().vars(1)); return true; } - // If all variables but one are fully encoded, - // force the last one to be fully encoded. - int variable_not_fully_encoded; - int num_fully_encoded = 0; - for (const int var : ct.linear().vars()) { - if (IsFullyEncoded(var)) { - num_fully_encoded++; - } else { - variable_not_fully_encoded = var; - } - } - if (num_fully_encoded == num_vars - 1) { - FullyEncode(variable_not_fully_encoded); - return true; - } - if (num_fully_encoded == num_vars) return true; - - // Register on remaining variables if not already done. - if (!constraint_is_registered_[ct_index]) { - for (const int var : ct.linear().vars()) { - if (!IsFullyEncoded(var)) Register(ct_index, var); - } - } - - return false; + return true; } void MaybeFullyEncodeMoreVariables(const CpModelProto& model_proto, Model* m) { @@ -843,18 +908,44 @@ void LoadEquivalenceAC(const std::vector enforcement_literal, } } +// Boolean encoding of: +// enforcement_literal => coeff1 * var1 + coeff2 * var2 != rhs; +void LoadEquivalenceNeqAC(const std::vector enforcement_literal, + IntegerValue coeff1, IntegerVariable var1, + IntegerValue coeff2, IntegerVariable var2, + const IntegerValue rhs, Model* m) { + auto* encoder = m->GetOrCreate(); + CHECK(encoder->VariableIsFullyEncoded(var1)); + CHECK(encoder->VariableIsFullyEncoded(var2)); + absl::flat_hash_map term1_value_to_literal; + for (const auto value_literal : encoder->FullDomainEncoding(var1)) { + term1_value_to_literal[coeff1 * value_literal.value] = + value_literal.literal; + } + for (const auto value_literal : encoder->FullDomainEncoding(var2)) { + const IntegerValue target_value = rhs - value_literal.value * coeff2; + if (gtl::ContainsKey(term1_value_to_literal, target_value)) { + const Literal target_literal = term1_value_to_literal[target_value]; + m->Add(EnforcedClause( + enforcement_literal, + {value_literal.literal.Negated(), target_literal.Negated()})); + } + } +} + } // namespace void LoadLinearConstraint(const ConstraintProto& ct, Model* m) { auto* mapping = m->GetOrCreate(); + auto* integer_trail = m->GetOrCreate(); const std::vector vars = mapping->Integers(ct.linear().vars()); const std::vector coeffs = ValuesFromProto(ct.linear().coeffs()); const SatParameters& params = *m->GetOrCreate(); - if (params.boolean_encoding_level() > 0 && vars.size() == 2 && - ct.linear().domain_size() == 2 && - ct.linear().domain(0) == ct.linear().domain(1)) { + int64 single_value = 0; + if (params.boolean_encoding_level() > 0 && ct.linear().vars_size() == 2 && + IsEqCst(ct.linear())) { auto* encoder = m->GetOrCreate(); if (encoder->VariableIsFullyEncoded(vars[0]) && encoder->VariableIsFullyEncoded(vars[1])) { @@ -864,11 +955,22 @@ void LoadLinearConstraint(const ConstraintProto& ct, Model* m) { IntegerValue(ct.linear().domain(0)), m); } } + if (params.boolean_encoding_level() > 0 && ct.linear().vars_size() == 2 && + IsNeqCst(ct.linear(), mapping, integer_trail, &single_value)) { + auto* encoder = m->GetOrCreate(); + if (encoder->VariableIsFullyEncoded(vars[0]) && + encoder->VariableIsFullyEncoded(vars[1])) { + VLOG(2) << "Extract " << ct.DebugString(); + return LoadEquivalenceNeqAC(mapping->Literals(ct.enforcement_literal()), + IntegerValue(coeffs[0]), vars[0], + IntegerValue(coeffs[1]), vars[1], + IntegerValue(single_value), m); + } + } // Compute the min/max to relax the bounds if needed. IntegerValue min_sum(0); IntegerValue max_sum(0); - auto* integer_trail = m->GetOrCreate(); bool all_booleans = true; for (int i = 0; i < vars.size(); ++i) { if (all_booleans && !mapping->IsBoolean(ct.linear().vars(i))) { diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 9891972bc3..706f2f07ef 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -290,6 +290,7 @@ std::string CpSolverResponseStats(const CpSolverResponse& response) { absl::StrAppend(&result, "\nusertime: ", response.user_time()); absl::StrAppend(&result, "\ndeterministic_time: ", response.deterministic_time()); + absl::StrAppend(&result, "\nprimal_integral: ", response.primal_integral()); absl::StrAppend(&result, "\n"); return result; } @@ -1508,9 +1509,11 @@ void PostsolveResponse(const std::string& debug_info, postsolve_model.Add(operations_research::sat::NewSatParameters(params)); } + std::unique_ptr time_limit(TimeLimit::Infinite()); + SharedTimeLimit shared_time_limit(time_limit.get()); SharedResponseManager local_response_manager( /*log_updates=*/false, /*enumerate_all_solutions=*/false, - /*solution_limit=*/-1, &mapping_proto, wall_timer); + /*solution_limit=*/-1, &mapping_proto, wall_timer, &shared_time_limit); LoadCpModel(mapping_proto, &local_response_manager, &postsolve_model); SolveLoadedCpModel(mapping_proto, &local_response_manager, &postsolve_model); const CpSolverResponse postsolve_response = @@ -1969,7 +1972,8 @@ class LnsSolver : public SubSolver { // maybe we can just randomize them like for the base solution used. SharedResponseManager local_response_manager( /*log_updates=*/false, /*enumerate_all_solutions=*/false, - /*solution_limit=*/-1, &neighborhood.cp_model, shared_->wall_timer); + /*solution_limit=*/-1, &neighborhood.cp_model, shared_->wall_timer, + shared_->time_limit); LoadCpModel(neighborhood.cp_model, &local_response_manager, &local_model); QuickSolveWithHint(neighborhood.cp_model, &local_response_manager, &local_model); @@ -2371,7 +2375,7 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { SharedResponseManager shared_response_manager( log_search, params.enumerate_all_solutions(), params.stop_after_first_solution() ? 1 : -1, &new_cp_model_proto, - &wall_timer); + &wall_timer, &shared_time_limit); const auto& observers = model->GetOrCreate()->observers; if (!observers.empty()) { shared_response_manager.AddSolutionCallback( diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index 66dc86b7c6..2c114944b0 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -14,6 +14,7 @@ #include "ortools/sat/synchronization.h" #include "absl/container/flat_hash_set.h" +#include "ortools/base/integral_types.h" #include "ortools/base/stl_util.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/cp_model_search.h" @@ -21,6 +22,7 @@ #include "ortools/sat/integer.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_base.h" +#include "ortools/util/time_limit.h" namespace operations_research { namespace sat { @@ -80,16 +82,16 @@ void SharedSolutionRepository::Synchronize() { } // TODO(user): Experiments and play with the num_solutions_to_keep parameter. -SharedResponseManager::SharedResponseManager(bool log_updates, - bool enumerate_all_solutions, - int solution_limit, - const CpModelProto* proto, - const WallTimer* wall_timer) +SharedResponseManager::SharedResponseManager( + bool log_updates, bool enumerate_all_solutions, int solution_limit, + const CpModelProto* proto, const WallTimer* wall_timer, + const SharedTimeLimit* shared_time_limit) : log_updates_(log_updates), enumerate_all_solutions_(enumerate_all_solutions), solution_limit_(solution_limit), model_proto_(*proto), wall_timer_(*wall_timer), + shared_time_limit_(*shared_time_limit), solutions_(/*num_solutions_to_keep=*/10) {} namespace { @@ -113,6 +115,23 @@ void LogNewSatSolution(const std::string& event_or_solution_count, } // namespace +void SharedResponseManager::UpdatePrimalIntegral(int64 max_integral) { + if (!model_proto_.has_objective()) return; + const double current_time = shared_time_limit_.GetElapsedDeterministicTime(); + const double time_delta = current_time - last_primal_integral_time_stamp_; + last_primal_integral_time_stamp_ = current_time; + const CpObjectiveProto& obj = model_proto_.objective(); + double bounds_delta = 0.0; + if (inner_objective_upper_bound_ == kint64max || + inner_objective_lower_bound_ == kint64min) { + bounds_delta = ScaleObjectiveValue(obj, max_integral); + } else { + bounds_delta = ScaleObjectiveValue( + obj, inner_objective_upper_bound_ - inner_objective_lower_bound_); + } + primal_integral_ += time_delta * std::abs(bounds_delta); +} + void SharedResponseManager::UpdateInnerObjectiveBounds( const std::string& worker_info, IntegerValue lb, IntegerValue ub) { absl::MutexLock mutex_lock(&mutex_); @@ -127,13 +146,20 @@ void SharedResponseManager::UpdateInnerObjectiveBounds( return; } - bool change = false; + const bool change = + (lb > inner_objective_lower_bound_ || ub < inner_objective_upper_bound_); + if (change) { + IntegerValue max_integral(0); + if (!AddProductTo(IntegerValue(2), ub - lb, &max_integral)) { + // Overflow. + max_integral = kMaxIntegerValue; + } + UpdatePrimalIntegral(/*max_integral=*/std::max(int64{0}, max_integral.value())); + } if (lb > inner_objective_lower_bound_) { - change = true; inner_objective_lower_bound_ = lb.value(); } if (ub < inner_objective_upper_bound_) { - change = true; inner_objective_upper_bound_ = ub.value(); } if (inner_objective_lower_bound_ > inner_objective_upper_bound_) { @@ -178,6 +204,7 @@ void SharedResponseManager::NotifyThatImprovingProblemIsInfeasible( CHECK_EQ(num_solutions_, 0); best_response_.set_status(CpSolverStatus::INFEASIBLE); } + UpdatePrimalIntegral(/*max_integral=*/kint64max); if (log_updates_) LogNewSatSolution("Done", wall_timer_.Get(), worker_info); } @@ -196,6 +223,11 @@ IntegerValue SharedResponseManager::BestSolutionInnerObjectiveValue() { return IntegerValue(best_solution_objective_value_); } +double SharedResponseManager::PrimalIntegral() const { + absl::MutexLock mutex_lock(&mutex_); + return primal_integral_; +} + int SharedResponseManager::AddSolutionCallback( std::function callback) { absl::MutexLock mutex_lock(&mutex_); @@ -280,6 +312,14 @@ void SharedResponseManager::NewSolution(const CpSolverResponse& response, CHECK_NE(best_response_.status(), CpSolverStatus::OPTIMAL); best_solution_objective_value_ = objective_value; + IntegerValue max_integral(0); + if (!AddProductTo(IntegerValue(2), IntegerValue(std::abs(objective_value)), + &max_integral)) { + // Overflow. + max_integral = kMaxIntegerValue; + } + UpdatePrimalIntegral(/*max_integral=*/max_integral.value()); + // Update the new bound. inner_objective_upper_bound_ = objective_value - 1; } @@ -349,6 +389,7 @@ void SharedResponseManager::SetStatsFromModelInternal(Model* model) { best_response_.set_num_booleans(sat_solver->NumVariables()); best_response_.set_num_branches(sat_solver->num_branches()); best_response_.set_num_conflicts(sat_solver->num_failures()); + best_response_.set_primal_integral(primal_integral_); best_response_.set_num_binary_propagations(sat_solver->num_propagations()); best_response_.set_num_integer_propagations( integer_trail == nullptr ? 0 : integer_trail->num_enqueues()); diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 07cdb5909c..b27bf22710 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -26,6 +26,7 @@ #include "ortools/sat/sat_base.h" #include "ortools/util/bitset.h" #include "ortools/util/random_engine.h" +#include "ortools/util/time_limit.h" namespace operations_research { namespace sat { @@ -72,7 +73,7 @@ class SharedTimeLimit { time_limit_->AdvanceDeterministicTime(deterministic_duration); } - double GetElapsedDeterministicTime() { + double GetElapsedDeterministicTime() const { absl::MutexLock mutex_lock(&mutex_); return time_limit_->GetElapsedDeterministicTime(); } @@ -156,7 +157,8 @@ class SharedResponseManager { // logged. This class is responsible for our solver log progress. SharedResponseManager(bool log_updates, bool enumerate_all_solutions, int solution_limit, const CpModelProto* proto, - const WallTimer* wall_timer); + const WallTimer* wall_timer, + const SharedTimeLimit* shared_time_limit); // Returns the current solver response. That is the best known response at the // time of the call with the best feasible solution and objective bounds. @@ -186,6 +188,8 @@ class SharedResponseManager { // there is no solution. IntegerValue BestSolutionInnerObjectiveValue(); + double PrimalIntegral() const; + // Updates the inner objective bounds. void UpdateInnerObjectiveBounds(const std::string& worker_info, IntegerValue lb, IntegerValue ub); @@ -230,11 +234,17 @@ class SharedResponseManager { void FillObjectiveValuesInBestResponse() EXCLUSIVE_LOCKS_REQUIRED(mutex_); void SetStatsFromModelInternal(Model* model) EXCLUSIVE_LOCKS_REQUIRED(mutex_); + // Updates the primal integral using the old bounds on the objective. If the + // old bounds are not finite, it uses the 'max_integral' value instead of gap. + void UpdatePrimalIntegral(int64 max_integral) + EXCLUSIVE_LOCKS_REQUIRED(mutex_); + const bool log_updates_; const bool enumerate_all_solutions_; const int solution_limit_; const CpModelProto& model_proto_; const WallTimer& wall_timer_; + const SharedTimeLimit& shared_time_limit_; mutable absl::Mutex mutex_; @@ -245,6 +255,8 @@ class SharedResponseManager { int64 inner_objective_lower_bound_ GUARDED_BY(mutex_) = kint64min; int64 inner_objective_upper_bound_ GUARDED_BY(mutex_) = kint64max; int64 best_solution_objective_value_ GUARDED_BY(mutex_) = kint64max; + double primal_integral_ GUARDED_BY(mutex_) = 0.0; + double last_primal_integral_time_stamp_ GUARDED_BY(mutex_) = 0.0; int next_callback_id_ GUARDED_BY(mutex_) = 0; std::vector>>