add new info in cp-sat solving

This commit is contained in:
Laurent Perron
2019-08-06 15:56:48 -07:00
parent c36d1d7dc1
commit 6c8cedcf22
5 changed files with 214 additions and 53 deletions

View File

@@ -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;

View File

@@ -56,6 +56,77 @@ std::vector<int64> ValuesFromProto(const Values& values) {
return std::vector<int64>(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<Literal> enforcement_literal,
}
}
// Boolean encoding of:
// enforcement_literal => coeff1 * var1 + coeff2 * var2 != rhs;
void LoadEquivalenceNeqAC(const std::vector<Literal> enforcement_literal,
IntegerValue coeff1, IntegerVariable var1,
IntegerValue coeff2, IntegerVariable var2,
const IntegerValue rhs, Model* m) {
auto* encoder = m->GetOrCreate<IntegerEncoder>();
CHECK(encoder->VariableIsFullyEncoded(var1));
CHECK(encoder->VariableIsFullyEncoded(var2));
absl::flat_hash_map<IntegerValue, Literal> 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<CpModelMapping>();
auto* integer_trail = m->GetOrCreate<IntegerTrail>();
const std::vector<IntegerVariable> vars =
mapping->Integers(ct.linear().vars());
const std::vector<int64> coeffs = ValuesFromProto(ct.linear().coeffs());
const SatParameters& params = *m->GetOrCreate<SatParameters>();
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<IntegerEncoder>();
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<IntegerEncoder>();
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<IntegerTrail>();
bool all_booleans = true;
for (int i = 0; i < vars.size(); ++i) {
if (all_booleans && !mapping->IsBoolean(ct.linear().vars(i))) {

View File

@@ -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<TimeLimit> 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<SolutionObservers>()->observers;
if (!observers.empty()) {
shared_response_manager.AddSolutionCallback(

View File

@@ -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<void(const CpSolverResponse&)> 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());

View File

@@ -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<std::pair<int, std::function<void(const CpSolverResponse&)>>>