always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

improve CP-SAT presolver, add Gomory (off by default) and MIR (naive, off by default) to the CP-SAT solver; improve CP-SAT logging

Browse Source
This commit is contained in:
Laurent Perron
2019-01-17 16:12:52 +01:00
parent 76a3d69498
commit 61a2b13fcf
15 changed files with 834 additions and 495 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
makefiles/Makefile.gen.mk
View File

@@ -1332,16 +1332,17 @@ objs/sat/cp_model_presolve.$O: ortools/sat/cp_model_presolve.cc \
ortools/base/commandlineflags.h ortools/base/logging.h \
ortools/base/integral_types.h ortools/base/macros.h ortools/base/timer.h \
ortools/base/basictypes.h ortools/util/running_stat.h \
ortools/base/hash.h ortools/base/map_util.h ortools/base/stl_util.h \
ortools/port/proto_utils.h ortools/sat/cp_model_checker.h \
ortools/sat/cp_model_loader.h ortools/base/int_type.h \
ortools/base/int_type_indexed_vector.h ortools/sat/cp_model_utils.h \
ortools/util/sorted_interval_list.h ortools/sat/integer.h \
ortools/graph/iterators.h ortools/sat/model.h ortools/base/typeid.h \
ortools/sat/sat_base.h ortools/util/bitset.h ortools/sat/sat_solver.h \
ortools/sat/clause.h ortools/sat/drat_proof_handler.h \
ortools/sat/drat_checker.h ortools/sat/drat_writer.h ortools/base/file.h \
ortools/base/status.h ortools/util/random_engine.h ortools/util/stats.h \
ortools/base/hash.h ortools/base/map_util.h ortools/base/mathutil.h \
ortools/base/stl_util.h ortools/port/proto_utils.h \
ortools/sat/cp_model_checker.h ortools/sat/cp_model_loader.h \
ortools/base/int_type.h ortools/base/int_type_indexed_vector.h \
ortools/sat/cp_model_utils.h ortools/util/sorted_interval_list.h \
ortools/sat/integer.h ortools/graph/iterators.h ortools/sat/model.h \
ortools/base/typeid.h ortools/sat/sat_base.h ortools/util/bitset.h \
ortools/sat/sat_solver.h ortools/sat/clause.h \
ortools/sat/drat_proof_handler.h ortools/sat/drat_checker.h \
ortools/sat/drat_writer.h ortools/base/file.h ortools/base/status.h \
ortools/util/random_engine.h ortools/util/stats.h \
ortools/sat/pb_constraint.h ortools/sat/restart.h \
ortools/sat/sat_decision.h ortools/util/integer_pq.h ortools/util/rev.h \
ortools/util/saturated_arithmetic.h ortools/sat/intervals.h \
@@ -3256,11 +3257,11 @@ objs/constraint_solver/routing_parameters.$O: \
objs/constraint_solver/routing_search.$O: \
ortools/constraint_solver/routing_search.cc ortools/base/small_map.h \
ortools/base/small_ordered_set.h ortools/base/stl_util.h \
ortools/base/integral_types.h ortools/base/macros.h \
ortools/constraint_solver/routing.h \
ortools/base/adjustable_priority_queue-inl.h \
ortools/base/adjustable_priority_queue.h ortools/base/basictypes.h \
ortools/base/integral_types.h ortools/base/logging.h \
ortools/base/macros.h ortools/base/commandlineflags.h \
ortools/base/logging.h ortools/base/commandlineflags.h \
ortools/base/hash.h ortools/base/int_type_indexed_vector.h \
ortools/base/int_type.h ortools/constraint_solver/constraint_solver.h \
ortools/base/map_util.h ortools/base/random.h ortools/base/sysinfo.h \

146
ortools/sat/cp_model_expand.cc
View File

@@ -26,22 +26,24 @@ namespace operations_research {
namespace sat {
namespace {
struct ExpansionHelper {
CpModelProto expanded_proto;
struct ExpansionContext {
CpModelProto working_model;
absl::flat_hash_map<std::pair<int, int>, int> precedence_cache;
std::map<std::string, int> statistics;
std::map<std::string, int> stats_by_rule_name;
static const int kAlwaysTrue = kint32min;
void UpdateRuleStats(const std::string& rule) { stats_by_rule_name[rule]++; }
// a => b.
void AddImplication(int a, int b) {
ConstraintProto* const ct = expanded_proto.add_constraints();
ConstraintProto* const ct = working_model.add_constraints();
ct->add_enforcement_literal(a);
ct->mutable_bool_and()->add_literals(b);
}
// b => x in [lb, ub].
void AddImplyInDomain(int b, int x, int64 lb, int64 ub) {
ConstraintProto* const imply = expanded_proto.add_constraints();
ConstraintProto* const imply = working_model.add_constraints();
imply->add_enforcement_literal(b);
imply->mutable_linear()->add_vars(x);
imply->mutable_linear()->add_coeffs(1);
@@ -50,17 +52,17 @@ struct ExpansionHelper {
}
int AddIntVar(int64 lb, int64 ub) {
IntegerVariableProto* const var = expanded_proto.add_variables();
IntegerVariableProto* const var = working_model.add_variables();
var->add_domain(lb);
var->add_domain(ub);
return expanded_proto.variables_size() - 1;
return working_model.variables_size() - 1;
}
int AddBoolVar() { return AddIntVar(0, 1); }
void AddBoolOr(const std::vector<int>& literals) {
BoolArgumentProto* const bool_or =
expanded_proto.add_constraints()->mutable_bool_or();
working_model.add_constraints()->mutable_bool_or();
for (const int lit : literals) {
bool_or->add_literals(lit);
}
@@ -85,7 +87,7 @@ struct ExpansionHelper {
// x_lesseq_y <=> (x <= y && l_x is true && l_y is true).
void AddReifiedPrecedence(int x_lesseq_y, int x, int y, int l_x, int l_y) {
// x_lesseq_y => (x <= y) && l_x is true && l_y is true.
ConstraintProto* const lesseq = expanded_proto.add_constraints();
ConstraintProto* const lesseq = working_model.add_constraints();
lesseq->add_enforcement_literal(x_lesseq_y);
lesseq->mutable_linear()->add_vars(x);
lesseq->mutable_linear()->add_vars(y);
@@ -101,7 +103,7 @@ struct ExpansionHelper {
}
// Not(x_lesseq_y) && l_x && l_y => (x > y)
ConstraintProto* const greater = expanded_proto.add_constraints();
ConstraintProto* const greater = working_model.add_constraints();
greater->mutable_linear()->add_vars(x);
greater->mutable_linear()->add_vars(y);
greater->mutable_linear()->add_coeffs(-1);
@@ -131,21 +133,21 @@ struct ExpansionHelper {
}
};
void ExpandReservoir(ConstraintProto* ct, ExpansionHelper* helper) {
void ExpandReservoir(ConstraintProto* ct, ExpansionContext* context) {
const ReservoirConstraintProto& reservoir = ct->reservoir();
const int num_variables = reservoir.times_size();
CpModelProto& expanded_proto = helper->expanded_proto;
CpModelProto& working_model = context->working_model;
auto is_optional = [&expanded_proto, &reservoir](int index) {
auto is_optional = [&working_model, &reservoir](int index) {
if (reservoir.actives_size() == 0) return false;
const int literal = reservoir.actives(index);
const int ref = PositiveRef(literal);
const IntegerVariableProto& var_proto = expanded_proto.variables(ref);
const IntegerVariableProto& var_proto = working_model.variables(ref);
return var_proto.domain_size() != 2 ||
var_proto.domain(0) != var_proto.domain(1);
};
auto active = [&reservoir, &helper](int index) {
if (reservoir.actives_size() == 0) return helper->kAlwaysTrue;
auto active = [&reservoir, &context](int index) {
if (reservoir.actives_size() == 0) return context->kAlwaysTrue;
return reservoir.actives(index);
};
@@ -167,20 +169,20 @@ void ExpandReservoir(ConstraintProto* ct, ExpansionHelper* helper) {
const int time_j = reservoir.times(j);
const std::pair<int, int> p = std::make_pair(time_i, time_j);
const std::pair<int, int> rev_p = std::make_pair(time_j, time_i);
if (gtl::ContainsKey(helper->precedence_cache, p)) continue;
if (gtl::ContainsKey(context->precedence_cache, p)) continue;
const int i_lesseq_j = helper->AddBoolVar();
helper->precedence_cache[p] = i_lesseq_j;
const int j_lesseq_i = helper->AddBoolVar();
helper->precedence_cache[rev_p] = j_lesseq_i;
helper->AddReifiedPrecedence(i_lesseq_j, time_i, time_j, active(i),
active(j));
helper->AddReifiedPrecedence(j_lesseq_i, time_j, time_i, active(j),
active(i));
const int i_lesseq_j = context->AddBoolVar();
context->precedence_cache[p] = i_lesseq_j;
const int j_lesseq_i = context->AddBoolVar();
context->precedence_cache[rev_p] = j_lesseq_i;
context->AddReifiedPrecedence(i_lesseq_j, time_i, time_j, active(i),
active(j));
context->AddReifiedPrecedence(j_lesseq_i, time_j, time_i, active(j),
active(i));
// Consistency. This is redundant but should improves performance.
auto* const bool_or =
expanded_proto.add_constraints()->mutable_bool_or();
working_model.add_constraints()->mutable_bool_or();
bool_or->add_literals(i_lesseq_j);
bool_or->add_literals(j_lesseq_i);
if (is_optional(i)) {
@@ -199,12 +201,12 @@ void ExpandReservoir(ConstraintProto* ct, ExpansionHelper* helper) {
for (int i = 0; i < num_variables; ++i) {
const int time_i = reservoir.times(i);
// Accumulates demands of all predecessors.
ConstraintProto* const level = expanded_proto.add_constraints();
ConstraintProto* const level = working_model.add_constraints();
for (int j = 0; j < num_variables; ++j) {
if (i == j) continue;
const int time_j = reservoir.times(j);
level->mutable_linear()->add_vars(gtl::FindOrDieNoPrint(
helper->precedence_cache, std::make_pair(time_j, time_i)));
context->precedence_cache, std::make_pair(time_j, time_i)));
level->mutable_linear()->add_coeffs(reservoir.demands(j));
}
// Accounts for own demand.
@@ -221,7 +223,7 @@ void ExpandReservoir(ConstraintProto* ct, ExpansionHelper* helper) {
// If all demands have the same sign, we do not care about the order, just
// the sum.
int64 fixed_demand = 0;
auto* const sum = expanded_proto.add_constraints()->mutable_linear();
auto* const sum = working_model.add_constraints()->mutable_linear();
for (int i = 0; i < num_variables; ++i) {
const int64 demand = reservoir.demands(i);
if (demand == 0) continue;
@@ -240,11 +242,11 @@ void ExpandReservoir(ConstraintProto* ct, ExpansionHelper* helper) {
// We need to do it only if 0 is not in [min_level..max_level].
// Otherwise, the regular propagation will already check it.
if (reservoir.min_level() > 0 || reservoir.max_level() < 0) {
auto* const initial_ct = expanded_proto.add_constraints()->mutable_linear();
auto* const initial_ct = working_model.add_constraints()->mutable_linear();
for (int i = 0; i < num_variables; ++i) {
const int time_i = reservoir.times(i);
const int lesseq_0 = helper->AddBoolVar();
helper->AddReifiedLessOrEqualThanZero(lesseq_0, time_i, active(i));
const int lesseq_0 = context->AddBoolVar();
context->AddReifiedLessOrEqualThanZero(lesseq_0, time_i, active(i));
initial_ct->add_vars(lesseq_0);
initial_ct->add_coeffs(reservoir.demands(i));
}
@@ -253,15 +255,15 @@ void ExpandReservoir(ConstraintProto* ct, ExpansionHelper* helper) {
}
ct->Clear();
helper->statistics["kReservoir"]++;
context->UpdateRuleStats("reservoir: expanded");
}
void ExpandIntMod(ConstraintProto* ct, ExpansionHelper* helper) {
void ExpandIntMod(ConstraintProto* ct, ExpansionContext* context) {
const IntegerArgumentProto& int_mod = ct->int_mod();
const IntegerVariableProto& var_proto =
helper->expanded_proto.variables(int_mod.vars(0));
context->working_model.variables(int_mod.vars(0));
const IntegerVariableProto& mod_proto =
helper->expanded_proto.variables(int_mod.vars(1));
context->working_model.variables(int_mod.vars(1));
const int target_var = int_mod.target();
const int64 mod_lb = mod_proto.domain(0);
@@ -272,19 +274,19 @@ void ExpandIntMod(ConstraintProto* ct, ExpansionHelper* helper) {
const int64 var_ub = var_proto.domain(var_proto.domain_size() - 1);
// Compute domains of var / mod_proto.
const int div_var = helper->AddIntVar(var_lb / mod_ub, var_ub / mod_lb);
const int div_var = context->AddIntVar(var_lb / mod_ub, var_ub / mod_lb);
auto add_enforcement_literal_if_needed = [&]() {
if (ct->enforcement_literal_size() == 0) return;
const int literal = ct->enforcement_literal(0);
ConstraintProto* const last = helper->expanded_proto.mutable_constraints(
helper->expanded_proto.constraints_size() - 1);
ConstraintProto* const last = context->working_model.mutable_constraints(
context->working_model.constraints_size() - 1);
last->add_enforcement_literal(literal);
};
// div = var / mod.
IntegerArgumentProto* const div_proto =
helper->expanded_proto.add_constraints()->mutable_int_div();
context->working_model.add_constraints()->mutable_int_div();
div_proto->set_target(div_var);
div_proto->add_vars(int_mod.vars(0));
div_proto->add_vars(int_mod.vars(1));
@@ -294,7 +296,7 @@ void ExpandIntMod(ConstraintProto* ct, ExpansionHelper* helper) {
if (mod_lb == mod_ub) {
// var - div_var * mod = target.
LinearConstraintProto* const lin =
helper->expanded_proto.add_constraints()->mutable_linear();
context->working_model.add_constraints()->mutable_linear();
lin->add_vars(int_mod.vars(0));
lin->add_coeffs(1);
lin->add_vars(div_var);
@@ -308,9 +310,9 @@ void ExpandIntMod(ConstraintProto* ct, ExpansionHelper* helper) {
// Create prod_var = div_var * mod.
const int mod_var = int_mod.vars(1);
const int prod_var =
helper->AddIntVar(var_lb * mod_lb / mod_ub, var_ub * mod_ub / mod_lb);
context->AddIntVar(var_lb * mod_lb / mod_ub, var_ub * mod_ub / mod_lb);
IntegerArgumentProto* const int_prod =
helper->expanded_proto.add_constraints()->mutable_int_prod();
context->working_model.add_constraints()->mutable_int_prod();
int_prod->set_target(prod_var);
int_prod->add_vars(div_var);
int_prod->add_vars(mod_var);
@@ -318,7 +320,7 @@ void ExpandIntMod(ConstraintProto* ct, ExpansionHelper* helper) {
// var - prod_var = target.
LinearConstraintProto* const lin =
helper->expanded_proto.add_constraints()->mutable_linear();
context->working_model.add_constraints()->mutable_linear();
lin->add_vars(int_mod.vars(0));
lin->add_coeffs(1);
lin->add_vars(prod_var);
@@ -331,12 +333,12 @@ void ExpandIntMod(ConstraintProto* ct, ExpansionHelper* helper) {
}
ct->Clear();
helper->statistics["kIntMod"]++;
context->UpdateRuleStats("int_mod: expanded");
}
void ExpandIntProdWithBoolean(int bool_ref, int int_ref, int product_ref,
ExpansionHelper* helper) {
ConstraintProto* const one = helper->expanded_proto.add_constraints();
ExpansionContext* context) {
ConstraintProto* const one = context->working_model.add_constraints();
one->add_enforcement_literal(bool_ref);
one->mutable_linear()->add_vars(int_ref);
one->mutable_linear()->add_coeffs(1);
@@ -345,7 +347,7 @@ void ExpandIntProdWithBoolean(int bool_ref, int int_ref, int product_ref,
one->mutable_linear()->add_domain(0);
one->mutable_linear()->add_domain(0);
ConstraintProto* const zero = helper->expanded_proto.add_constraints();
ConstraintProto* const zero = context->working_model.add_constraints();
zero->add_enforcement_literal(NegatedRef(bool_ref));
zero->mutable_linear()->add_vars(product_ref);
zero->mutable_linear()->add_coeffs(1);
@@ -353,15 +355,15 @@ void ExpandIntProdWithBoolean(int bool_ref, int int_ref, int product_ref,
zero->mutable_linear()->add_domain(0);
}
void ExpandIntProd(ConstraintProto* ct, ExpansionHelper* helper) {
void ExpandIntProd(ConstraintProto* ct, ExpansionContext* context) {
const IntegerArgumentProto& int_prod = ct->int_prod();
if (int_prod.vars_size() != 2) return;
const int a = int_prod.vars(0);
const int b = int_prod.vars(1);
const IntegerVariableProto& a_proto =
helper->expanded_proto.variables(PositiveRef(a));
context->working_model.variables(PositiveRef(a));
const IntegerVariableProto& b_proto =
helper->expanded_proto.variables(PositiveRef(b));
context->working_model.variables(PositiveRef(b));
const int p = int_prod.target();
const bool a_is_boolean = RefIsPositive(a) && a_proto.domain_size() == 2 &&
a_proto.domain(0) == 0 && a_proto.domain(1) == 1;
@@ -371,49 +373,53 @@ void ExpandIntProd(ConstraintProto* ct, ExpansionHelper* helper) {
// We expand if exactly one of {a, b} is Boolean. If both are Boolean, it
// will be presolved into a better version.
if (a_is_boolean && !b_is_boolean) {
ExpandIntProdWithBoolean(a, b, p, helper);
ExpandIntProdWithBoolean(a, b, p, context);
ct->Clear();
helper->statistics["kIntProd"]++;
context->UpdateRuleStats("int_prod: expanded product with Boolean var");
} else if (b_is_boolean && !a_is_boolean) {
ExpandIntProdWithBoolean(b, a, p, helper);
ExpandIntProdWithBoolean(b, a, p, context);
ct->Clear();
helper->statistics["kIntProd"]++;
context->UpdateRuleStats("int_prod: expanded product with Boolean var");
}
}
} // namespace
CpModelProto ExpandCpModel(const CpModelProto& initial_model) {
ExpansionHelper helper;
helper.expanded_proto = initial_model;
const int num_constraints = helper.expanded_proto.constraints_size();
CpModelProto ExpandCpModel(const CpModelProto& initial_model, bool log) {
ExpansionContext context;
context.working_model = initial_model;
const int num_constraints = context.working_model.constraints_size();
for (int i = 0; i < num_constraints; ++i) {
ConstraintProto* const ct = helper.expanded_proto.mutable_constraints(i);
ConstraintProto* const ct = context.working_model.mutable_constraints(i);
switch (ct->constraint_case()) {
case ConstraintProto::ConstraintCase::kReservoir:
ExpandReservoir(ct, &helper);
ExpandReservoir(ct, &context);
break;
case ConstraintProto::ConstraintCase::kIntMod:
ExpandIntMod(ct, &helper);
ExpandIntMod(ct, &context);
break;
case ConstraintProto::ConstraintCase::kIntProd:
ExpandIntProd(ct, &helper);
ExpandIntProd(ct, &context);
break;
default:
break;
}
}
for (const auto& entry : helper.statistics) {
if (entry.second == 1) {
VLOG(1) << "Expanded 1 '" << entry.first << "' constraint.";
} else {
VLOG(1) << "Expanded " << entry.second << " '" << entry.first
<< "' constraints";
if (log) {
std::map<std::string, int> sorted_rules(context.stats_by_rule_name.begin(),
context.stats_by_rule_name.end());
for (const auto& entry : sorted_rules) {
if (entry.second == 1) {
LOG(INFO) << "- rule '" << entry.first << "' was applied 1 time.";
} else {
LOG(INFO) << "- rule '" << entry.first << "' was applied "
<< entry.second << " times.";
}
}
}
return helper.expanded_proto;
return context.working_model;
}
} // namespace sat
} // namespace operations_research

2
ortools/sat/cp_model_expand.h
View File

@@ -23,7 +23,7 @@ namespace sat {
// simpler constraints.
// This is different from PresolveCpModel() as there are no reduction or
// simplification of the model. Furthermore, this expansion is mandatory.
CpModelProto ExpandCpModel(const CpModelProto& initial_model);
CpModelProto ExpandCpModel(const CpModelProto& initial_model, bool log);
} // namespace sat
} // namespace operations_research

12
ortools/sat/cp_model_loader.cc
View File

@@ -55,16 +55,6 @@ std::vector<int64> ValuesFromProto(const Values& values) {
return std::vector<int64>(values.begin(), values.end());
}
// Returns the size of the given domain capped to int64max.
int64 DomainSize(const Domain& domain) {
int64 size = 0;
for (const ClosedInterval interval : domain) {
size += operations_research::CapAdd(
1, operations_research::CapSub(interval.end, interval.start));
}
return size;
}
} // namespace
void CpModelMapping::CreateVariables(const CpModelProto& model_proto,
@@ -389,7 +379,7 @@ void CpModelMapping::ExtractEncoding(const CpModelProto& model_proto,
//
// TODO(user): Also fully encode variable that are almost fully encoded?
const Domain domain = ReadDomainFromProto(model_proto.variables(i));
if (DomainSize(domain) == values.size()) {
if (domain.Size() == values.size()) {
++num_fully_encoded;
if (!encoder->VariableIsFullyEncoded(integers_[i])) {
encoder->FullyEncodeVariable(integers_[i]);

231
ortools/sat/cp_model_presolve.cc
View File

@@ -31,6 +31,7 @@
#include "ortools/base/integral_types.h"
#include "ortools/base/logging.h"
#include "ortools/base/map_util.h"
#include "ortools/base/mathutil.h"
#include "ortools/base/stl_util.h"
#include "ortools/port/proto_utils.h"
#include "ortools/sat/cp_model.pb.h"
@@ -305,6 +306,9 @@ struct PresolveContext {
var_to_constraints.resize(domains.size());
}
// Returns the number of mapped domains.
int NumDomainsStored() const { return domains.size(); }
// This regroup all the affine relations between variables. Note that the
// constraints used to detect such relations will not be removed from the
// model at detection time (thus allowing proper domain propagation). However,
@@ -1564,14 +1568,39 @@ bool PresolveTable(ConstraintProto* ct, PresolveContext* context) {
std::vector<std::vector<int64>> new_tuples;
new_tuples.reserve(num_tuples);
std::vector<absl::flat_hash_set<int64>> new_domains(num_vars);
std::vector<AffineRelation::Relation> affine_relations;
bool modified_variables = false;
for (int v = 0; v < num_vars; ++v) {
const int var = ct->table().vars(v);
if (!RefIsPositive(var)) {
affine_relations.push_back({var, 1, 0}); // Support opposite.
} else {
affine_relations.push_back(context->GetAffineRelation(var));
if (affine_relations[v].representative != var) {
modified_variables = true;
}
}
}
for (int i = 0; i < num_tuples; ++i) {
bool delete_row = false;
std::string tmp;
for (int j = 0; j < num_vars; ++j) {
const int ref = ct->table().vars(j);
const int64 v = ct->table().values(i * num_vars + j);
int64 v = ct->table().values(i * num_vars + j);
const AffineRelation::Relation& r = affine_relations[j];
if (r.representative != ref) {
const int64 inverse_value = (v - r.offset) / r.coeff;
if (inverse_value * r.coeff + r.offset != v) {
// Bad rounding.
delete_row = true;
break;
}
v = inverse_value;
}
tuple[j] = v;
if (!context->DomainOf(ref).Contains(v)) {
if (!context->DomainOf(r.representative).Contains(v)) {
delete_row = true;
break;
}
@@ -1587,14 +1616,26 @@ bool PresolveTable(ConstraintProto* ct, PresolveContext* context) {
gtl::STLSortAndRemoveDuplicates(&new_tuples);
// Update the list of tuples if needed.
if (new_tuples.size() < num_tuples) {
if (new_tuples.size() < num_tuples || modified_variables) {
ct->mutable_table()->clear_values();
for (const std::vector<int64>& t : new_tuples) {
for (const int64 v : t) {
ct->mutable_table()->add_values(v);
}
}
context->UpdateRuleStats("table: removed rows");
if (new_tuples.size() < num_tuples) {
context->UpdateRuleStats("table: removed rows");
}
}
if (modified_variables) {
for (int j = 0; j < num_vars; ++j) {
const AffineRelation::Relation& r = affine_relations[j];
if (r.representative != ct->table().vars(j)) {
ct->mutable_table()->set_vars(j, r.representative);
}
}
context->UpdateRuleStats("table: replace variable by canonical affine one");
}
// Filter the variable domains.
@@ -1655,7 +1696,7 @@ bool PresolveTable(ConstraintProto* ct, PresolveContext* context) {
}
context->UpdateRuleStats("table: negated");
}
return false;
return modified_variables;
}
bool PresolveAllDiff(ConstraintProto* ct, PresolveContext* context) {
@@ -2533,6 +2574,92 @@ void MergeNoOverlapConstraints(PresolveContext* context) {
}
}
// Extracts cliques from bool_and and small at_most_one constraints and
// transforms them into maximal cliques.
void TransformIntoMaxCliques(PresolveContext* context) {
if (context->is_unsat) return;
auto convert = [](int ref) {
if (RefIsPositive(ref)) return Literal(BooleanVariable(ref), true);
return Literal(BooleanVariable(NegatedRef(ref)), false);
};
const int num_constraints = context->working_model->constraints_size();
// Extract the bool_and and at_most_one constraints.
std::vector<std::vector<Literal>> cliques;
for (int c = 0; c < num_constraints; ++c) {
ConstraintProto* ct = context->working_model->mutable_constraints(c);
if (ct->constraint_case() == ConstraintProto::ConstraintCase::kAtMostOne) {
std::vector<Literal> clique;
// Process only small constraints.
for (const int ref : ct->at_most_one().literals()) {
clique.push_back(convert(ref));
}
cliques.push_back(clique);
if (RemoveConstraint(ct, context)) {
context->UpdateConstraintVariableUsage(c);
}
} else if (ct->constraint_case() ==
ConstraintProto::ConstraintCase::kBoolAnd) {
if (ct->enforcement_literal().size() != 1) continue;
const Literal enforcement = convert(ct->enforcement_literal(0));
for (const int ref : ct->bool_and().literals()) {
cliques.push_back({enforcement, convert(ref).Negated()});
}
if (RemoveConstraint(ct, context)) {
context->UpdateConstraintVariableUsage(c);
}
}
}
const int old_cliques = cliques.size();
// We reuse the max-clique code from sat.
Model local_model;
auto* graph = local_model.GetOrCreate<BinaryImplicationGraph>();
const int num_variables = context->working_model->variables().size();
graph->Resize(num_variables);
for (const std::vector<Literal>& clique : cliques) {
if (clique.size() <= 100) graph->AddAtMostOne(clique);
}
if (!graph->DetectEquivalences()) {
context->is_unsat = true;
return;
}
graph->TransformIntoMaxCliques(&cliques);
// Add the Boolean variable equivalence detected by DetectEquivalences().
// Those are needed because TransformIntoMaxCliques() will replace all
// variable by its representative.
for (int var = 0; var < num_variables; ++var) {
const Literal l = Literal(BooleanVariable(var), true);
if (graph->RepresentativeOf(l) != l) {
const Literal r = graph->RepresentativeOf(l);
context->AddBooleanEqualityRelation(
var, r.IsPositive() ? r.Variable().value()
: NegatedRef(r.Variable().value()));
}
}
int new_cliques = 0;
for (const std::vector<Literal>& clique : cliques) {
if (clique.empty()) continue;
new_cliques++;
ConstraintProto* ct = context->working_model->add_constraints();
for (const Literal literal : clique) {
if (literal.IsPositive()) {
ct->mutable_at_most_one()->add_literals(literal.Variable().value());
} else {
ct->mutable_at_most_one()->add_literals(
NegatedRef(literal.Variable().value()));
}
}
}
context->UpdateNewConstraintsVariableUsage();
VLOG(1) << "Merged " << old_cliques << " into " << new_cliques << " cliques";
}
bool PresolveOneConstraint(int c, PresolveContext* context) {
ConstraintProto* ct = context->working_model->mutable_constraints(c);
@@ -2801,6 +2928,83 @@ bool ProcessSetPPC(PresolveContext* context, TimeLimit* time_limit) {
return changed;
}
void TryToSimplifyDomains(PresolveContext* context) {
if (context->is_unsat) return;
for (int var = 0; var < context->NumDomainsStored(); ++var) {
if (context->IsFixed(var)) continue;
const AffineRelation::Relation r = context->GetAffineRelation(var);
if (r.representative != var) continue;
const Domain& domain = context->DomainOf(var);
if (domain.NumIntervals() == 1) continue;
if (domain.Size() == 2 && (domain.Min() != 0 || domain.Max() != 1)) {
const int index = context->working_model->variables_size();
IntegerVariableProto* const var_proto =
context->working_model->add_variables();
var_proto->add_domain(0);
var_proto->add_domain(1);
ConstraintProto* const ct = context->working_model->add_constraints();
LinearConstraintProto* const lin = ct->mutable_linear();
lin->add_vars(var);
lin->add_coeffs(1);
lin->add_vars(index);
lin->add_coeffs(domain.Min() - domain.Max());
lin->add_domain(domain.Min());
lin->add_domain(domain.Min());
context->InitializeNewDomains();
context->UpdateNewConstraintsVariableUsage();
context->AddAffineRelation(*ct, var, index, domain.Max() - domain.Min(),
domain.Min());
context->UpdateRuleStats("variables: canonicalize domain of size 2");
} else if (domain.NumIntervals() > 2 &&
domain.NumIntervals() == domain.Size()) { // Discrete domain.
const int64 var_min = domain.Min();
std::vector<int64> diffs(domain.NumIntervals());
int64 gcd = -1;
bool ok = true;
for (int index = 0; index < domain.NumIntervals(); ++index) {
const ClosedInterval& i = domain[index];
CHECK_EQ(i.start, i.end);
diffs[index] = i.start - var_min;
CHECK_GE(diffs[index], 0);
if (gcd == -1) {
gcd = diffs[index];
} else {
gcd = MathUtil::GCD64(gcd, diffs[index]);
}
if (gcd == 1) {
ok = false;
break;
}
}
if (!ok) {
continue;
}
std::vector<int64> scaled_values;
for (const int64 value : diffs) {
scaled_values.push_back(value / gcd);
}
const int index = context->working_model->variables_size();
IntegerVariableProto* const var_proto =
context->working_model->add_variables();
const Domain scaled_domain = Domain::FromValues(scaled_values);
FillDomainInProto(scaled_domain, var_proto);
ConstraintProto* const ct = context->working_model->add_constraints();
LinearConstraintProto* const lin = ct->mutable_linear();
lin->add_vars(var);
lin->add_coeffs(1);
lin->add_vars(index);
lin->add_coeffs(-gcd);
lin->add_domain(var_min);
lin->add_domain(var_min);
context->InitializeNewDomains();
context->UpdateNewConstraintsVariableUsage();
context->AddAffineRelation(*ct, var, index, gcd, var_min);
context->UpdateRuleStats("variables: canonicalize affine domain");
}
}
}
void PresolveToFixPoint(PresolveContext* context, TimeLimit* time_limit) {
if (context->is_unsat) return;
@@ -2842,6 +3046,20 @@ void PresolveToFixPoint(PresolveContext* context, TimeLimit* time_limit) {
}
}
// Look at variables to see if we can canonicalize the domain
const int num_constraints_before_simplify =
context->working_model->constraints_size();
TryToSimplifyDomains(context);
const int num_constraints_after_simplify =
context->working_model->constraints_size();
if (num_constraints_after_simplify > num_constraints_before_simplify) {
in_queue.resize(num_constraints_after_simplify, true);
for (int c = num_constraints_before_simplify;
c < num_constraints_after_simplify; ++c) {
queue.push_back(c);
}
}
// Re-add to the queue constraints that have unique variables. Note that to
// not enter an infinite loop, we call each (var, constraint) pair at most
// once.
@@ -2980,7 +3198,6 @@ void RemoveUnusedEquivalentVariables(PresolveContext* context) {
// Update the variable usage.
context->UpdateNewConstraintsVariableUsage();
}
} // namespace.
// =============================================================================
@@ -3059,6 +3276,8 @@ bool PresolveCpModel(const PresolveOptions& options,
PresolvePureSatPart(&context);
}
TransformIntoMaxCliques(&context);
// Process set packing, partitioning and covering constraint.
if (options.time_limit == nullptr || !options.time_limit->LimitReached()) {
ProcessSetPPC(&context, options.time_limit);

29
ortools/sat/cp_model_solver.cc
View File

@@ -357,6 +357,7 @@ std::string Summarize(const std::string& input) {
std::string CpModelStats(const CpModelProto& model_proto) {
std::map<std::string, int> num_constraints_by_name;
std::map<std::string, int> num_reif_constraints_by_name;
std::map<std::string, int> name_to_num_literals;
for (const ConstraintProto& ct : model_proto.constraints()) {
std::string name = ConstraintCaseName(ct.constraint_case());
@@ -372,6 +373,18 @@ std::string CpModelStats(const CpModelProto& model_proto) {
if (!ct.enforcement_literal().empty()) {
num_reif_constraints_by_name[name]++;
}
// For pure Boolean constraints, we also display the total number of literal
// involved as this gives a good idea of the problem size.
if (ct.constraint_case() == ConstraintProto::ConstraintCase::kBoolOr) {
name_to_num_literals[name] += ct.bool_or().literals().size();
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kBoolAnd) {
name_to_num_literals[name] += ct.bool_and().literals().size();
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kAtMostOne) {
name_to_num_literals[name] += ct.at_most_one().literals().size();
}
}
int num_constants = 0;
@@ -445,10 +458,16 @@ std::string CpModelStats(const CpModelProto& model_proto) {
std::vector<std::string> constraints;
constraints.reserve(num_constraints_by_name.size());
for (const auto entry : num_constraints_by_name) {
constraints.push_back(
absl::StrCat("#", entry.first, ": ", entry.second, " (",
num_reif_constraints_by_name[entry.first],
" with enforcement literal)"));
const std::string& name = entry.first;
constraints.push_back(absl::StrCat("#", name, ": ", entry.second));
if (gtl::ContainsKey(num_reif_constraints_by_name, name)) {
absl::StrAppend(&constraints.back(),
" (#enforced: ", num_reif_constraints_by_name[name], ")");
}
if (gtl::ContainsKey(name_to_num_literals, name)) {
absl::StrAppend(&constraints.back(),
" (#literals: ", name_to_num_literals[name], ")");
}
}
std::sort(constraints.begin(), constraints.end());
absl::StrAppend(&result, absl::StrJoin(constraints, "\n"));
@@ -2316,7 +2335,7 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) {
}
// Starts by expanding some constraints if needed.
CpModelProto new_model = ExpandCpModel(model_proto);
CpModelProto new_model = ExpandCpModel(model_proto, VLOG_IS_ON(1));
// Presolve?
std::function<void(CpSolverResponse * response)> postprocess_solution;

328
ortools/sat/cuts.cc
View File

@@ -553,74 +553,69 @@ CutGenerator CreateKnapsackCoverCutGenerator(
return result;
}
namespace {
// Returns 0 if there is none.
// Note that we normalize the fractionality by its coefficient.
IntegerValue MagnitudeOfMostFractionalVariable(
const std::vector<double>& lp_values, const LinearConstraint& cut) {
double best_score = 0;
IntegerValue best_magnitude(0);
int num_fractional_vars = 0;
const int size = lp_values.size();
for (int i = 0; i < size; ++i) {
const IntegerValue coeff = cut.coeffs[i];
const double value = lp_values[i];
const double score =
std::abs(value - std::round(value)) * ToDouble(IntTypeAbs(coeff));
if (score > best_score) {
best_score = score;
best_magnitude = IntTypeAbs(coeff);
}
if (std::abs(value - std::round(value)) > 0.01) {
++num_fractional_vars;
VLOG(3) << "value: " << value << " coeff: " << coeff
<< " score:" << score;
}
}
VLOG(2) << "num_fractional_vars: " << num_fractional_vars << "/" << size;
return best_magnitude;
}
} // namespace
std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
IntegerValue remainder, IntegerValue divisor, IntegerValue max_scaling) {
const IntegerValue target = CeilRatio(divisor, IntegerValue(2)) - 1;
const IntegerValue t = std::max(
IntegerValue(1),
remainder == 0 ? IntegerValue(1)
: std::min(max_scaling / 2, CeilRatio(target, remainder)));
const IntegerValue threshold = std::max(target, t * remainder);
return [t, threshold, divisor](IntegerValue coeff) {
const IntegerValue ratio = FloorRatio(t * coeff, divisor);
const IntegerValue remainder = t * coeff - ratio * divisor;
return 2 * ratio + (remainder > threshold ? 1 : 0);
};
}
bool use_letchford_lodi_version, IntegerValue rhs_remainder,
IntegerValue divisor, IntegerValue max_scaling) {
// Compute the larger t <= max_scaling such that
// t * rhs_remainder >= divisor / 2.
const IntegerValue t =
rhs_remainder == 0
? max_scaling
: std::min(max_scaling, CeilRatio(divisor / 2, rhs_remainder));
// Adjust after the multiplication by t.
rhs_remainder *= t;
max_scaling /= t;
// This is the only difference compared to a discretized MIR function.
if (use_letchford_lodi_version && max_scaling > 2) max_scaling = 2;
// TODO(user): if 2 * rhs_remainder < divisor, multiply by a factor t before
// rounding.
std::function<IntegerValue(IntegerValue)> GetMirFunction(
IntegerValue rhs_remainder, IntegerValue divisor,
IntegerValue max_scaling) {
CHECK_GE(max_scaling, 1);
if (divisor - rhs_remainder <= max_scaling) {
return [rhs_remainder, divisor](IntegerValue coeff) {
const IntegerValue ratio = FloorRatio(coeff, divisor);
const IntegerValue remainder = coeff - ratio * divisor;
return (divisor - rhs_remainder) * ratio +
std::max(IntegerValue(0), remainder - rhs_remainder);
const IntegerValue size = divisor - rhs_remainder;
if (max_scaling == 1) {
// TODO(user): Use everywhere a two step computation to avoid overflow?
// First divide by divisor, then multiply by t.
return [t, divisor](IntegerValue coeff) {
return FloorRatio(t * coeff, divisor);
};
} else if (size <= max_scaling) {
return [size, rhs_remainder, t, divisor](IntegerValue coeff) {
const IntegerValue ratio = FloorRatio(t * coeff, divisor);
const IntegerValue remainder = t * coeff - ratio * divisor;
const IntegerValue diff = remainder - rhs_remainder;
return size * ratio + std::max(IntegerValue(0), diff);
};
} else {
// TODO(user): This function is not maximal, improve?
return [rhs_remainder, divisor, max_scaling](IntegerValue coeff) {
const IntegerValue ratio = FloorRatio(coeff, divisor);
const IntegerValue remainder = coeff - ratio * divisor;
return max_scaling * ratio +
std::max(FloorRatio((remainder - rhs_remainder) * max_scaling,
divisor - rhs_remainder),
IntegerValue(0));
// We divide (size = divisor - rhs_remainder) into (max_scaling - 1) buckets
// and increase the function by 1 / max_scaling for each of them.
//
// Note that for different values of max_scaling, we get a family of
// functions that do not dominate each others. So potentially, a max scaling
// as low as 2 could lead to the better cut (this is exactly the Letchford &
// Lodi function).
///
// Another intersting fact, is that if we want to compute the maximum alpha
// for a constraint with 2 terms like:
// divisor * Y + (ratio * divisor + remainder) * X
// <= rhs_ratio * divisor + rhs_remainder
// so that we have the cut:
// Y + (ratio + alpha) * X <= rhs_ratio
// This is the same as computing the maximum alpha such that for all integer
// X > 0 we have CeilRatio(alpha * divisor * X, divisor)
// <= CeilRatio(remainder * X - rhs_remainder, divisor).
// We can prove that this alpha is of the form (n - 1) / n, and it will
// be reached by such function for a max_scaling of n.
//
// TODO(user): This function is not always maximal when
// size % (max_scaling - 1) == 0. Improve?
return [size, rhs_remainder, t, divisor, max_scaling](IntegerValue coeff) {
const IntegerValue ratio = FloorRatio(t * coeff, divisor);
const IntegerValue remainder = t * coeff - ratio * divisor;
const IntegerValue diff = remainder - rhs_remainder;
const IntegerValue bucket =
diff > 0 ? CeilRatio(diff * (max_scaling - 1), size)
: IntegerValue(0);
return max_scaling * ratio + bucket;
};
}
}
@@ -630,6 +625,7 @@ void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
std::vector<IntegerValue> upper_bounds,
LinearConstraint* cut) {
const int size = lp_values.size();
if (size == 0) return;
CHECK_EQ(lower_bounds.size(), size);
CHECK_EQ(upper_bounds.size(), size);
CHECK_EQ(cut->vars.size(), size);
@@ -653,67 +649,50 @@ void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
}
}
// Find the magnitude of the most fractional variable, note that we normalize
// the fractionality by its coefficient.
const IntegerValue divisor =
MagnitudeOfMostFractionalVariable(lp_values, *cut);
if (divisor == 0) {
VLOG(1) << "Issue, no fractional variables.";
*cut = LinearConstraint(IntegerValue(0), IntegerValue(0));
return;
}
// To simplify the code below, we make all coefficients positive.
std::vector<bool> change_sign_at_postprocessing(size, false);
for (int i = 0; i < size; ++i) {
if (cut->coeffs[i] > 0) continue;
change_sign_at_postprocessing[i] = true;
cut->coeffs[i] = -cut->coeffs[i];
lp_values[i] = -lp_values[i];
std::swap(lower_bounds[i], upper_bounds[i]);
lower_bounds[i] = -lower_bounds[i];
upper_bounds[i] = -upper_bounds[i];
}
// Shift each variable using its lower/upper bound so that no variable can
// change sign.
// change sign. We eventually do a change of variable to its negation so
// that all variable are non-negative.
bool overflow = false;
std::vector<IntegerValue> shifts(size, IntegerValue(0));
std::vector<bool> var_is_positive_or_zero(size, true);
std::vector<bool> change_sign_at_postprocessing(size, false);
IntegerValue max_initial_magnitude(1);
for (int i = 0; i < size && !overflow; ++i) {
const IntegerValue coeff = cut->coeffs[i];
if (coeff == 0) continue;
if (cut->coeffs[i] == 0) continue;
// Note that since we use ToDouble() this code works fine with lb/ub at
// min/max integer value.
const double value = lp_values[i];
const IntegerValue lb = lower_bounds[i];
const IntegerValue ub = upper_bounds[i];
if (std::abs(value - ToDouble(lb)) < std::abs(value - ToDouble(ub))) {
// We want coeff * (X - lb) so the new var is >= 0.
var_is_positive_or_zero[i] = true;
shifts[i] = lb;
} else {
// We want coeff * (X - ub) so the new var is <= 0.
var_is_positive_or_zero[i] = false;
shifts[i] = ub;
{
const double value = lp_values[i];
const IntegerValue lb = lower_bounds[i];
const IntegerValue ub = upper_bounds[i];
if (std::abs(value - ToDouble(lb)) > std::abs(value - ToDouble(ub))) {
// Change the variable sign.
change_sign_at_postprocessing[i] = true;
cut->coeffs[i] = -cut->coeffs[i];
lp_values[i] = -lp_values[i];
std::swap(lower_bounds[i], upper_bounds[i]);
lower_bounds[i] = -lower_bounds[i];
upper_bounds[i] = -upper_bounds[i];
}
}
// Always shift to lb.
// coeff * X = coeff * (X - shift) + coeff * shift.
if (!AddProductTo(-coeff, shifts[i], &cut->ub)) {
lp_values[i] -= ToDouble(lower_bounds[i]);
if (!AddProductTo(-cut->coeffs[i], lower_bounds[i], &cut->ub)) {
overflow = true;
break;
}
// Deal with fixed variable, no need to shift back in this case, we can
// just remove the term.
if (lb == ub) {
shifts[i] = IntegerValue(0);
if (lower_bounds[i] == upper_bounds[i]) {
cut->coeffs[i] = IntegerValue(0);
lp_values[i] = 0.0;
}
max_initial_magnitude =
std::max(max_initial_magnitude, IntTypeAbs(cut->coeffs[i]));
}
if (overflow) {
VLOG(1) << "Issue, overflow.";
@@ -721,66 +700,103 @@ void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
return;
}
// We will adjust coefficient that are close to an exact multiple of divisor
// to an exact multiple. This is meant to get rid of small errors that appears
// due to rounding error in our exact computation of the initial constraint
// given to this class.
//
// TODO(user): Tune the threshold or use a parameter. Maybe it should depend
// on the number of term in the constraint? But the basic idea is that we do
// not want to change the rhs_remainder (see below) by too much. So here we
// change it at most by: num_terms * 0.0002. Note that in practice we don't
// except a lot of terms to be close to a multiple of divisor.
const IntegerValue adjust_threshold = divisor / IntegerValue(5000);
// Our heuristic will try to generate a few different cuts, and we will keep
// the most violated one.
double best_scaled_violation = 0.01;
LinearConstraint best_cut(IntegerValue(0), IntegerValue(0));
for (int i = 0; i < size; ++i) {
const IntegerValue coeff = cut->coeffs[i];
const IntegerValue diff(
CapSub(upper_bounds[i].value(), lower_bounds[i].value()));
if (var_is_positive_or_zero[i]) {
// Adjust coeff of the form k * divisor - epsilon.
const IntegerValue remainder =
CeilRatio(coeff, divisor) * divisor - coeff;
if (CapProd(diff.value(), remainder.value()) > adjust_threshold) continue;
cut->ub += remainder * diff;
cut->coeffs[i] += remainder;
} else {
// Adjust coeff of the form k * divisor + epsilon.
const IntegerValue remainder =
coeff - FloorRatio(coeff, divisor) * divisor;
if (CapProd(diff.value(), remainder.value()) > adjust_threshold) continue;
cut->ub += remainder * diff;
cut->coeffs[i] -= remainder;
// Skip shifted variable almost at their lower bound.
if (lp_values[i] <= 1e-4) continue;
const IntegerValue divisor = IntTypeAbs(cut->coeffs[i]);
// Skip if we don't have the potential to generate a good enough cut.
const IntegerValue initial_rhs_remainder =
cut->ub - FloorRatio(cut->ub, divisor) * divisor;
if (ToDouble(initial_rhs_remainder) / ToDouble(max_initial_magnitude) <=
best_scaled_violation) {
continue;
}
}
// Create the super-additive function f().
const IntegerValue rhs_remainder =
cut->ub - FloorRatio(cut->ub, divisor) * divisor;
const auto f =
options.use_mir
? GetMirFunction(rhs_remainder, divisor, options.max_scaling)
: GetSuperAdditiveRoundingFunction(rhs_remainder, divisor,
options.max_scaling);
// TODO(user): We could avoid this copy.
LinearConstraint temp_cut = *cut;
// Apply f() to the cut.
cut->ub = f(cut->ub);
for (int i = 0; i < cut->coeffs.size(); ++i) {
const IntegerValue coeff = cut->coeffs[i];
if (coeff == 0) continue;
if (var_is_positive_or_zero[i]) {
cut->coeffs[i] = f(coeff);
} else {
cut->coeffs[i] = -f(-coeff);
// We will adjust coefficient that are just under an exact multiple of
// divisor to an exact multiple. This is meant to get rid of small errors
// that appears due to rounding error in our exact computation of the
// initial constraint given to this class.
//
// Each adjustement will cause the initial_rhs_remainder to increase, and we
// do not want to increase it above divisor. Our threshold below guarantees
// this. Note that the higher the rhs_remainder becomes, the more the
// function f() has a chance to reduce the violation, so it is not always a
// good idea to use all the slack we have between initial_rhs_remainder and
// divisor.
//
// TODO(user): If possible, it might be better to complement these
// variables. Even if the adjusted lp_values end up larger, if we loose less
// when taking f(), then we will have a better violation.
const IntegerValue adjust_threshold =
(divisor - initial_rhs_remainder - 1) / IntegerValue(size);
if (adjust_threshold > 0) {
for (int i = 0; i < size; ++i) {
const IntegerValue coeff = temp_cut.coeffs[i];
const IntegerValue diff(
CapSub(upper_bounds[i].value(), lower_bounds[i].value()));
// Adjust coeff of the form k * divisor - epsilon.
const IntegerValue remainder =
CeilRatio(coeff, divisor) * divisor - coeff;
if (CapProd(diff.value(), remainder.value()) > adjust_threshold) {
continue;
}
temp_cut.ub += remainder * diff;
temp_cut.coeffs[i] += remainder;
}
}
// Create the super-additive function f().
const IntegerValue rhs_remainder =
temp_cut.ub - FloorRatio(temp_cut.ub, divisor) * divisor;
if (rhs_remainder == 0) continue;
const auto f = GetSuperAdditiveRoundingFunction(
!options.use_mir, rhs_remainder, divisor, options.max_scaling);
// Apply f() to the cut and compute the cut violation.
temp_cut.ub = f(temp_cut.ub);
double violation = -ToDouble(temp_cut.ub);
double max_magnitude = 1.0;
for (int i = 0; i < temp_cut.coeffs.size(); ++i) {
const IntegerValue coeff = temp_cut.coeffs[i];
if (coeff == 0) continue;
const IntegerValue new_coeff = f(coeff);
temp_cut.coeffs[i] = new_coeff;
max_magnitude = std::max(max_magnitude, std::abs(ToDouble(new_coeff)));
violation += ToDouble(new_coeff) * lp_values[i];
}
violation /= max_magnitude;
if (violation > 0.0) {
VLOG(2) << "lp_value: " << lp_values[i] << " divisor: " << divisor
<< " cut_violation: " << violation;
}
if (violation > best_scaled_violation) {
best_scaled_violation = violation;
best_cut = temp_cut;
}
}
// Remove the bound shifts so the constraint is expressed in the original
// variables and do some basic post-processing.
for (int i = 0; i < size; ++i) {
*cut = best_cut;
for (int i = 0; i < cut->coeffs.size(); ++i) {
const IntegerValue coeff = cut->coeffs[i];
if (coeff == 0) continue;
cut->ub = IntegerValue(
CapAdd((cut->coeffs[i] * shifts[i]).value(), cut->ub.value()));
CapAdd((coeff * lower_bounds[i]).value(), cut->ub.value()));
if (change_sign_at_postprocessing[i]) {
cut->coeffs[i] = -cut->coeffs[i];
cut->coeffs[i] = -coeff;
}
}
RemoveZeroTerms(cut);

30
ortools/sat/cuts.h
View File

@@ -45,7 +45,7 @@ struct CutGenerator {
// Visible for testing. Returns a function f on integers such that:
// - f is non-decreasing.
// - f is super-additive: f(a) + f(b) <= f(a + b)
// - 1 <= f(divisor) <= 2 * max_scaling
// - 1 <= f(divisor) <= max_scaling
// - For all x, f(x * divisor) = x * f(divisor)
// - For all x, f(x * divisor + remainder) = x * f(divisor)
//
@@ -57,19 +57,23 @@ struct CutGenerator {
// the cut. Just taking f(x) = x / divisor result in the non-strengthened cut
// and using any function that stricly dominate this one is better.
//
// Note(user): I could have used the MIR (Mixed integer rounding function), but
// I prefered to use the function described in "Strenghtening Chvatal-Gomory
// cuts and Gomory fractional cuts", Adam N. Letchfrod, Andrea Lodi. This is
// because it gives better result for small value of max_scaling, and we do not
// want to have large integer coefficient.
// Algorithm:
// - We first scale by a factor t so that rhs_remainder >= divisor / 2.
// - Then, if use_letchford_lodi_version is true, we use the function described
// in "Strenghtening Chvatal-Gomory cuts and Gomory fractional cuts", Adam N.
// Letchfrod, Andrea Lodi.
// - Otherwise, we use a generalization of this which is a discretized version
// of the classical MIR rounding function that only take the value of the
// form "an_integer / max_scaling". As max_scaling goes to infinity, this
// converge to the real-valued MIR function.
//
// TODO(user): Still not clear to me if this can be improved or not.
// Note that for each value of max_scaling we will get a different function.
// And that there is no dominance relation between any of these functions. So
// it could be nice to try to generate a cut using different values of
// max_scaling.
std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue max_scaling);
// Same as GetSuperAdditiveRoundingFunction() but uses the classic MIR one.
std::function<IntegerValue(IntegerValue)> GetMirFunction(
IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue max_scaling);
bool use_letchford_lodi_version, IntegerValue rhs_remainder,
IntegerValue divisor, IntegerValue max_scaling);
// Given an upper bounded linear constraint, this function tries to transform it
// to a valid cut that violate the given LP solution using integer rounding.
@@ -102,7 +106,7 @@ std::function<IntegerValue(IntegerValue)> GetMirFunction(
// path, so we can spend more effort in finding good cuts.
struct RoundingOptions {
bool use_mir = false;
IntegerValue max_scaling = IntegerValue(1000);
IntegerValue max_scaling = IntegerValue(60);
};
void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
std::vector<IntegerValue> lower_bounds,

47
ortools/sat/linear_constraint_manager.cc
View File

@@ -134,32 +134,35 @@ void LinearConstraintManager::Add(const LinearConstraint& ct) {
void LinearConstraintManager::AddCut(
const LinearConstraint& ct, std::string type_name,
const gtl::ITIVector<IntegerVariable, double>& lp_solution) {
CHECK_NE(ct.vars.size(), 0);
if (ct.vars.empty()) return;
IntegerValue max_magnitude(0);
double activity = 0.0;
double l2_norm = 0.0;
const int size = ct.vars.size();
for (int i = 0; i < size; ++i) {
max_magnitude = std::max(max_magnitude, IntTypeAbs(ct.coeffs[i]));
const double coeff = ToDouble(ct.coeffs[i]);
activity += coeff * lp_solution[ct.vars[i]];
l2_norm += coeff * coeff;
}
l2_norm = sqrt(l2_norm);
double violation = 0.0;
violation = std::max(violation, activity - ToDouble(ct.ub));
violation = std::max(violation, ToDouble(ct.lb) - activity);
// Only add cut with sufficient efficacity.
if (violation / l2_norm < 1e-5) return;
Add(ct);
num_cuts_++;
type_to_num_cuts_[type_name]++;
if (VLOG_IS_ON(1)) {
IntegerValue max_magnitude(0);
double activity = 0.0;
double l2_norm = 0.0;
const int size = ct.vars.size();
for (int i = 0; i < size; ++i) {
max_magnitude = std::max(max_magnitude, IntTypeAbs(ct.coeffs[i]));
const double coeff = ToDouble(ct.coeffs[i]);
activity += coeff * lp_solution[ct.vars[i]];
l2_norm += coeff * coeff;
}
l2_norm = sqrt(l2_norm);
double violation = 0.0;
violation = std::max(violation, activity - ToDouble(ct.ub));
violation = std::max(violation, ToDouble(ct.lb) - activity);
LOG(INFO) << "Cut '" << type_name << "'"
<< " size=" << size << " max_magnitude=" << max_magnitude
<< " norm=" << l2_norm << " violation=" << violation
<< " eff=" << violation / l2_norm;
}
VLOG(1) << "Cut '" << type_name << "'"
<< " size=" << size << " max_magnitude=" << max_magnitude
<< " norm=" << l2_norm << " violation=" << violation
<< " eff=" << violation / l2_norm;
}
bool LinearConstraintManager::ChangeLp(

401
ortools/sat/linear_programming_constraint.cc
View File

@@ -375,8 +375,164 @@ bool AddLinearExpressionMultiple(
} // namespace
// TODO(user): By heuristically choosing the lp_multipliers we can easily
// adapt this code to generate MIR cuts. To investigate.
void LinearProgrammingConstraint::AddCutFromConstraints(
const std::string& name,
const std::vector<std::pair<RowIndex, IntegerValue>>& integer_multipliers) {
// This is initialized to a valid linear contraint (by taking linear
// combination of the LP rows) and will be transformed into a cut if
// possible.
//
// TODO(user): Ideally this linear combination should have only one
// fractional variable (basis_col). But because of imprecision, we get a
// bunch of fractional entry with small coefficient (relative to the one of
// basis_col). We try to handle that in IntegerRoundingCut(), but it might
// be better to add small multiple of the involved rows to get rid of them.
LinearConstraint cut;
{
gtl::ITIVector<ColIndex, IntegerValue> dense_cut;
IntegerValue cut_ub;
if (!ComputeNewLinearConstraint(
/*use_constraint_status=*/true, integer_multipliers, &dense_cut,
&cut_ub)) {
VLOG(1) << "Issue, overflow!";
return;
}
cut = ConvertToLinearConstraint(dense_cut, cut_ub);
}
// This should be tight!
if (std::abs(ComputeActivity(cut, expanded_lp_solution_) - ToDouble(cut.ub)) >
0.1) {
VLOG(1) << "Not tight " << ComputeActivity(cut, expanded_lp_solution_)
<< " " << ToDouble(cut.ub);
return;
}
// Fills data for IntegerRoundingCut().
//
// Note(user): we use the current bound here, so the reasonement will only
// produce locally valid cut if we call this at a non-root node. We could
// use the level zero bounds if we wanted to generate a globally valid cut
// at another level, but we will likely not genereate a constraint violating
// the current lp solution in that case.
std::vector<double> lp_values;
std::vector<IntegerValue> var_lbs;
std::vector<IntegerValue> var_ubs;
for (const IntegerVariable var : cut.vars) {
lp_values.push_back(expanded_lp_solution_[var]);
var_lbs.push_back(integer_trail_->LowerBound(var));
var_ubs.push_back(integer_trail_->UpperBound(var));
}
// Add slack.
// definition: integer_lp_[row] + slack_row == bound;
const IntegerVariable first_slack(expanded_lp_solution_.size());
for (const auto pair : integer_multipliers) {
const RowIndex row = pair.first;
const IntegerValue coeff = pair.second;
const auto status = simplex_.GetConstraintStatus(row);
if (status == glop::ConstraintStatus::FIXED_VALUE) continue;
lp_values.push_back(0.0);
cut.vars.push_back(first_slack + IntegerVariable(row.value()));
cut.coeffs.push_back(coeff);
const IntegerValue diff(CapSub(integer_lp_[row.value()].ub.value(),
integer_lp_[row.value()].lb.value()));
if (status == glop::ConstraintStatus::AT_UPPER_BOUND) {
var_lbs.push_back(IntegerValue(0));
var_ubs.push_back(diff);
} else {
CHECK_EQ(status, glop::ConstraintStatus::AT_LOWER_BOUND);
var_lbs.push_back(-diff);
var_ubs.push_back(IntegerValue(0));
}
}
// Get the cut using some integer rounding heuristic.
RoundingOptions options;
options.use_mir = sat_parameters_.use_mir_rounding();
options.max_scaling = sat_parameters_.max_integer_rounding_scaling();
IntegerRoundingCut(options, lp_values, var_lbs, var_ubs, &cut);
// Compute the activity. Warning: the cut no longer have the same size so we
// cannot use lp_values. Note that the substitution below shouldn't change
// the activity by definition.
double activity = 0.0;
for (int i = 0; i < cut.vars.size(); ++i) {
if (cut.vars[i] < first_slack) {
activity += ToDouble(cut.coeffs[i]) * expanded_lp_solution_[cut.vars[i]];
}
}
const double kMinViolation = 1e-4;
const double violation = activity - ToDouble(cut.ub);
if (violation < kMinViolation) {
VLOG(2) << "Bad cut " << activity << " <= " << ToDouble(cut.ub);
return;
}
// Substitute any slack left.
{
int num_slack = 0;
gtl::ITIVector<ColIndex, IntegerValue> dense_cut(integer_variables_.size(),
IntegerValue(0));
IntegerValue cut_ub = cut.ub;
bool overflow = false;
for (int i = 0; i < cut.vars.size(); ++i) {
if (cut.vars[i] < first_slack) {
CHECK(VariableIsPositive(cut.vars[i]));
const glop::ColIndex col =
gtl::FindOrDie(mirror_lp_variable_, cut.vars[i]);
dense_cut[col] = cut.coeffs[i];
} else {
++num_slack;
// Update the constraint.
const glop::RowIndex row(cut.vars[i].value() - first_slack.value());
const IntegerValue multiplier = -cut.coeffs[i];
if (!AddLinearExpressionMultiple(
multiplier, integer_lp_[row.value()].terms, &dense_cut)) {
overflow = true;
break;
}
// Update rhs.
const auto status = simplex_.GetConstraintStatus(row);
if (status == glop::ConstraintStatus::AT_LOWER_BOUND) {
if (!AddProductTo(multiplier, integer_lp_[row.value()].lb, &cut_ub)) {
overflow = true;
break;
}
} else {
CHECK_EQ(status, glop::ConstraintStatus::AT_UPPER_BOUND);
if (!AddProductTo(multiplier, integer_lp_[row.value()].ub, &cut_ub)) {
overflow = true;
break;
}
}
}
}
if (overflow) {
VLOG(1) << "Overflow in slack removal.";
return;
}
VLOG(2) << " num_slack: " << num_slack;
cut = ConvertToLinearConstraint(dense_cut, cut_ub);
}
const double new_violation =
ComputeActivity(cut, expanded_lp_solution_) - ToDouble(cut.ub);
if (std::abs(violation - new_violation) >= 1e-4) {
VLOG(1) << "Violation discrepancy after slack removal. "
<< " before = " << violation << " after = " << new_violation;
}
DivideByGCD(&cut);
constraint_manager_.AddCut(cut, name, expanded_lp_solution_);
}
void LinearProgrammingConstraint::AddCGCuts() {
CHECK_EQ(trail_->CurrentDecisionLevel(), 0);
const RowIndex num_rows = lp_data_.num_constraints();
@@ -423,169 +579,28 @@ void LinearProgrammingConstraint::AddCGCuts() {
}
if (num_non_zeros == 0) continue;
// This is initialized to a valid linear contraint (by taking linear
// combination of the LP rows) and will be transformed into a cut if
// possible.
//
// TODO(user): Ideally this linear combination should have only one
// fractional variable (basis_col). But because of imprecision, we get a
// bunch of fractional entry with small coefficient (relative to the one of
// basis_col). We try to handle that in IntegerRoundingCut(), but it might
// be better to add small multiple of the involved rows to get rid of them.
LinearConstraint cut;
Fractional scaling;
const std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers =
ScaleLpMultiplier(/*take_objective_into_account=*/false,
/*use_constraint_status=*/true, lp_multipliers,
&scaling);
AddCutFromConstraints("CG", integer_multipliers);
}
}
void LinearProgrammingConstraint::AddMirCuts() {
CHECK_EQ(trail_->CurrentDecisionLevel(), 0);
const RowIndex num_rows = lp_data_.num_constraints();
for (RowIndex row(0); row < num_rows; ++row) {
const auto status = simplex_.GetConstraintStatus(row);
if (status == glop::ConstraintStatus::BASIC) continue;
// TODO(user): Do not consider just one constraint, but take linear
// combination of a small number of constraints. There is a lot of
// literature on the possible heuristics here.
std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
{
Fractional scaling;
gtl::ITIVector<ColIndex, IntegerValue> dense_cut;
IntegerValue cut_ub;
if (!ComputeNewLinearConstraint(
/*take_objective_into_account=*/false,
/*use_constraint_status=*/true, lp_multipliers, &scaling,
&integer_multipliers, &dense_cut, &cut_ub)) {
VLOG(1) << "Issue, overflow!";
continue;
}
cut = ConvertToLinearConstraint(dense_cut, cut_ub);
}
VLOG(2) << " lp_multipliers num_non_zeros: " << num_non_zeros
<< " magnitude: " << magnitude
<< " num_after_scaling: " << integer_multipliers.size();
// This should be tight!
if (std::abs(ComputeActivity(cut, expanded_lp_solution_) -
ToDouble(cut.ub)) > 0.1) {
VLOG(1) << "Not tight " << ComputeActivity(cut, expanded_lp_solution_)
<< " " << ToDouble(cut.ub);
continue;
}
// Fills data for IntegerRoundingCut().
//
// Note(user): we use the current bound here, so the reasonement will only
// produce locally valid cut if we call this at a non-root node. We could
// use the level zero bounds if we wanted to generate a globally valid cut
// at another level, but we will likely not genereate a constraint violating
// the current lp solution in that case.
std::vector<double> lp_values;
std::vector<IntegerValue> var_lbs;
std::vector<IntegerValue> var_ubs;
for (const IntegerVariable var : cut.vars) {
lp_values.push_back(expanded_lp_solution_[var]);
var_lbs.push_back(integer_trail_->LowerBound(var));
var_ubs.push_back(integer_trail_->UpperBound(var));
}
// Add slack.
// definition: integer_lp_[row] + slack_row == bound;
const IntegerVariable first_slack(expanded_lp_solution_.size());
for (const auto pair : integer_multipliers) {
const RowIndex row = pair.first;
const IntegerValue coeff = pair.second;
const auto status = simplex_.GetConstraintStatus(row);
if (status == glop::ConstraintStatus::FIXED_VALUE) continue;
lp_values.push_back(0.0);
cut.vars.push_back(first_slack + IntegerVariable(row.value()));
cut.coeffs.push_back(coeff);
const IntegerValue diff(CapSub(integer_lp_[row.value()].ub.value(),
integer_lp_[row.value()].lb.value()));
if (status == glop::ConstraintStatus::AT_UPPER_BOUND) {
var_lbs.push_back(IntegerValue(0));
var_ubs.push_back(diff);
} else {
CHECK_EQ(status, glop::ConstraintStatus::AT_LOWER_BOUND);
var_lbs.push_back(-diff);
var_ubs.push_back(IntegerValue(0));
}
}
// Get the cut using some integer rounding heuristic.
RoundingOptions options;
options.use_mir = sat_parameters_.use_mir_rounding();
IntegerRoundingCut(options, lp_values, var_lbs, var_ubs, &cut);
// Compute the activity. Warning: the cut no longer have the same size so we
// cannot use lp_values. Note that the substitution below shouldn't change
// the activity by definition.
double activity = 0.0;
for (int i = 0; i < cut.vars.size(); ++i) {
if (cut.vars[i] < first_slack) {
activity +=
ToDouble(cut.coeffs[i]) * expanded_lp_solution_[cut.vars[i]];
}
}
const double kMinViolation = 1e-4;
const double violation = activity - ToDouble(cut.ub);
if (violation < kMinViolation) {
VLOG(2) << "Bad cut " << activity << " <= " << ToDouble(cut.ub);
continue;
}
// Substitute any slack left.
{
int num_slack = 0;
gtl::ITIVector<ColIndex, IntegerValue> dense_cut(
integer_variables_.size(), IntegerValue(0));
IntegerValue cut_ub = cut.ub;
bool overflow = false;
for (int i = 0; i < cut.vars.size(); ++i) {
if (cut.vars[i] < first_slack) {
CHECK(VariableIsPositive(cut.vars[i]));
const glop::ColIndex col =
gtl::FindOrDie(mirror_lp_variable_, cut.vars[i]);
dense_cut[col] = cut.coeffs[i];
} else {
++num_slack;
// Update the constraint.
const glop::RowIndex row(cut.vars[i].value() - first_slack.value());
const IntegerValue multiplier = -cut.coeffs[i];
if (!AddLinearExpressionMultiple(
multiplier, integer_lp_[row.value()].terms, &dense_cut)) {
overflow = true;
break;
}
// Update rhs.
const auto status = simplex_.GetConstraintStatus(row);
if (status == glop::ConstraintStatus::AT_LOWER_BOUND) {
if (!AddProductTo(multiplier, integer_lp_[row.value()].lb,
&cut_ub)) {
overflow = true;
break;
}
} else {
CHECK_EQ(status, glop::ConstraintStatus::AT_UPPER_BOUND);
if (!AddProductTo(multiplier, integer_lp_[row.value()].ub,
&cut_ub)) {
overflow = true;
break;
}
}
}
}
if (overflow) {
VLOG(1) << "Overflow in slack removal.";
continue;
}
VLOG(2) << " num_slack: " << num_slack;
cut = ConvertToLinearConstraint(dense_cut, cut_ub);
}
const double new_violation =
ComputeActivity(cut, expanded_lp_solution_) - ToDouble(cut.ub);
if (std::abs(violation - new_violation) < 1e-4) {
VLOG(1) << "Violation discrepancy after slack removal. "
<< " before = " << violation << " after = " << new_violation;
if (new_violation < kMinViolation) continue;
}
DivideByGCD(&cut);
constraint_manager_.AddCut(cut, "CG", expanded_lp_solution_);
integer_multipliers.push_back({row, IntegerValue(1)});
AddCutFromConstraints("MIR1", integer_multipliers);
}
}
@@ -635,16 +650,21 @@ bool LinearProgrammingConstraint::Propagate() {
if (constraint_manager_.ChangeLp(expanded_lp_solution_)) {
CreateLpFromConstraintManager();
if (!SolveLp()) return true;
} else {
} else if (constraint_manager_.num_cuts() <
sat_parameters_.max_num_cuts()) {
const int old_num_cuts = constraint_manager_.num_cuts();
if (sat_parameters_.add_cg_cuts() &&
trail_->CurrentDecisionLevel() == 0) {
AddCGCuts();
// The "generic" cuts are currently part of this class as they are using
// data from the current LP.
//
// TODO(user): Refactor so that they are just normal cut generators?
if (trail_->CurrentDecisionLevel() == 0) {
if (sat_parameters_.add_mir_cuts()) AddMirCuts();
if (sat_parameters_.add_cg_cuts()) AddCGCuts();
}
// Try to add cuts.
if (!cut_generators_.empty() &&
constraint_manager_.num_cuts() < sat_parameters_.max_num_cuts() &&
(trail_->CurrentDecisionLevel() == 0 ||
!sat_parameters_.only_add_cuts_at_level_zero())) {
for (const CutGenerator& generator : cut_generators_) {
@@ -868,12 +888,10 @@ void LinearProgrammingConstraint::SetImpliedLowerBoundReason(
}
// TODO(user): Provide a sparse interface.
bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
std::vector<std::pair<RowIndex, IntegerValue>>
LinearProgrammingConstraint::ScaleLpMultiplier(
bool take_objective_into_account, bool use_constraint_status,
const glop::DenseColumn& dense_lp_multipliers, Fractional* scaling,
std::vector<std::pair<RowIndex, IntegerValue>>* integer_multipliers,
gtl::ITIVector<ColIndex, IntegerValue>* dense_terms,
IntegerValue* upper_bound) const {
const glop::DenseColumn& dense_lp_multipliers, Fractional* scaling) const {
// Process the dense_lp_multipliers and compute their infinity norm.
std::vector<std::pair<RowIndex, Fractional>> lp_multipliers;
Fractional lp_multipliers_norm = take_objective_into_account ? 1.0 : 0.0;
@@ -947,12 +965,19 @@ bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
// TODO(user): Divide dual by gcd to limit overflow?
// TODO(user): To avoid overflow, we could lower scaling at the cost of
// loosing precision.
integer_multipliers->clear();
std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
for (const auto entry : cp_multipliers) {
const IntegerValue coeff(std::round(entry.second * (*scaling)));
if (coeff != 0) integer_multipliers->push_back({entry.first, coeff});
if (coeff != 0) integer_multipliers.push_back({entry.first, coeff});
}
return integer_multipliers;
}
bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
bool use_constraint_status,
const std::vector<std::pair<RowIndex, IntegerValue>>& integer_multipliers,
gtl::ITIVector<ColIndex, IntegerValue>* dense_terms,
IntegerValue* upper_bound) const {
// Initialize the new constraint.
*upper_bound = 0;
dense_terms->assign(integer_variables_.size(), IntegerValue(0));
@@ -960,7 +985,7 @@ bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
// Compute the new constraint by taking the linear combination given by
// integer_multipliers of the integer constraints in integer_lp_.
const ColIndex num_cols(integer_variables_.size());
for (const std::pair<RowIndex, IntegerValue> term : *integer_multipliers) {
for (const std::pair<RowIndex, IntegerValue> term : integer_multipliers) {
const RowIndex row = term.first;
const IntegerValue multiplier = term.second;
CHECK_LT(row, integer_lp_.size());
@@ -1022,13 +1047,16 @@ bool LinearProgrammingConstraint::ExactLpReasonning() {
}
Fractional scaling;
const std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers =
ScaleLpMultiplier(/*take_objective_into_account=*/true,
/*use_constraint_status=*/false, lp_multipliers,
&scaling);
gtl::ITIVector<ColIndex, IntegerValue> reduced_costs;
IntegerValue rc_ub;
std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
if (!ComputeNewLinearConstraint(
/*take_objective_into_account=*/true,
/*use_constraint_status=*/false, lp_multipliers, &scaling,
&integer_multipliers, &reduced_costs, &rc_ub)) {
/*use_constraint_status=*/false, integer_multipliers, &reduced_costs,
&rc_ub)) {
VLOG(2) << "Overflow during exact LP reasoning.";
return true;
}
@@ -1078,13 +1106,16 @@ bool LinearProgrammingConstraint::ExactLpReasonning() {
bool LinearProgrammingConstraint::FillExactDualRayReason() {
Fractional scaling;
const std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers =
ScaleLpMultiplier(/*take_objective_into_account=*/false,
/*use_constraint_status=*/false, simplex_.GetDualRay(),
&scaling);
gtl::ITIVector<ColIndex, IntegerValue> dense_new_constraint;
IntegerValue new_constraint_ub;
std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
if (!ComputeNewLinearConstraint(
/*take_objective_into_account=*/false,
/*use_constraint_status=*/false, simplex_.GetDualRay(), &scaling,
&integer_multipliers, &dense_new_constraint, &new_constraint_ub)) {
/*use_constraint_status=*/false, integer_multipliers,
&dense_new_constraint, &new_constraint_ub)) {
return false;
}
const LinearConstraint new_constraint =

33
ortools/sat/linear_programming_constraint.h
View File

@@ -155,10 +155,22 @@ class LinearProgrammingConstraint : public PropagatorInterface,
// Solve the LP, returns false if something went wrong in the LP solver.
bool SolveLp();
// Add a "MIR" cut obtained by first taking the linear combination of the
// row of the matrix according to "integer_multipliers" and then trying
// some integer rounding heuristic.
void AddCutFromConstraints(
const std::string& name,
const std::vector<std::pair<glop::RowIndex, IntegerValue>>&
integer_multipliers);
// Computes and adds Chvatal-Gomory cuts.
// This can currently only be called at the root node.
void AddCGCuts();
// Computes and adds MIR cuts.
// This can currently only be called at the root node.
void AddMirCuts();
// The factor to multiply a CP variable value to get the value in the LP side.
glop::Fractional CpToLpScalingFactor(glop::ColIndex col) const;
glop::Fractional LpToCpScalingFactor(glop::ColIndex col) const;
@@ -189,17 +201,22 @@ class LinearProgrammingConstraint : public PropagatorInterface,
// world and then make them integer (eventually multiplying them by a new
// scaling factor returned in *scaling).
//
// Then computes from this linear combination of the integer rows of the LP a
// new constraint of the form "sum terms <= upper_bound". Note that whatever
// lp_multipliers are given, the constraint will always be an exact valid
// constraint of the problem.
// Note that this will loose some precision, but our subsequent computation
// will still be exact as it will work for any set of multiplier.
std::vector<std::pair<glop::RowIndex, IntegerValue>> ScaleLpMultiplier(
bool take_objective_into_account, bool use_constraint_status,
const glop::DenseColumn& dense_lp_multipliers,
glop::Fractional* scaling) const;
// Computes from an integer linear combination of the integer rows of the LP a
// new constraint of the form "sum terms <= upper_bound". All computation are
// exact here.
//
// Returns false if we encountered any integer overflow.
bool ComputeNewLinearConstraint(
bool take_objective_into_account, // For the scaling.
bool use_constraint_status, const glop::DenseColumn& dense_lp_multipliers,
glop::Fractional* scaling,
std::vector<std::pair<glop::RowIndex, IntegerValue>>* integer_multipliers,
bool use_constraint_status,
const std::vector<std::pair<glop::RowIndex, IntegerValue>>&
integer_multipliers,
gtl::ITIVector<glop::ColIndex, IntegerValue>* dense_terms,
IntegerValue* upper_bound) const;

23
ortools/sat/sat_parameters.proto
View File

@@ -21,7 +21,7 @@ package operations_research.sat;
// Contains the definitions for all the sat algorithm parameters and their
// default values.
//
// NEXT TAG: 119
// NEXT TAG: 121
message SatParameters {
// ==========================================================================
// Branching and polarity
@@ -497,21 +497,36 @@ message SatParameters {
// For now, we don't have any manager for the LP cuts, so we just add any
// generated cuts until this limit is reached.
optional int32 max_num_cuts = 91 [default = 1000];
// For the cut that can be generated at any level, this control if we only
// try to generate them at the root node.
optional bool only_add_cuts_at_level_zero = 92 [default = false];
// Whether we generate knapsack cuts. Note that in our setting where all
// variables are integer and bounded on both side, such a cut could be applied
// to any constraint.
optional bool add_knapsack_cuts = 111 [default = false];
// Whether we generate and add Chvatal-Gomory cuts to the LP at root node.
// Note that for now, this is not heavily tunned.
//
// TODO(user): When we add more cuts, better to refactorize this in a list
// of enum that list enabled cuts, like "cuts:[KNAPSACK,GOMORY,CIRCUIT]".
optional bool add_cg_cuts = 117 [default = false];
// Whether we generate MIR cuts at root node.
// Note that for now, this is not heavily tunned.
optional bool add_mir_cuts = 120 [default = false];
// Whether we use the classical MIR rounding function when computing a cut
// using integer rounding. If false we use our own variant based on the strong
// fractional rouding of Letchford and Lodi.
optional bool use_mir_rounding = 118 [default = true];
// In the integer rounding procedure used for MIR and Gomory cut, the maximum
// "scaling" we use (must be positive). The lower this is, the lower the
// integer coefficients of the cut will be. Note that cut generated by lower
// values are not necessarily worse than cut generated by larger value. There
// is no strict dominance relationship.
optional int32 max_integer_rounding_scaling = 119 [default = 600];
// If true, we start by an empty LP, and only add constraints not satisfied
// by the current LP solution batch by batch. A constraint that is only added
// like this is known as a "lazy" constraint in the literature, except that we

10
ortools/sat/synchronization.cc
View File

@@ -218,6 +218,8 @@ void RegisterObjectiveBestBoundExport(
const CpObjectiveProto& obj = model_proto.objective();
const double new_best_bound = ScaleObjectiveValue(
obj, integer_trail->LevelZeroLowerBound(objective_var).value());
const double new_objective_value = ScaleObjectiveValue(
obj, integer_trail->LevelZeroUpperBound(objective_var).value());
const double current_best_bound = helper->get_external_best_bound();
const double current_objective_value =
helper->get_external_best_objective();
@@ -228,12 +230,16 @@ void RegisterObjectiveBestBoundExport(
(helper->scaling_factor < 0 &&
new_best_bound < current_best_bound)) {
if (log_progress) {
const double reported_objective_value =
std::isfinite(current_objective_value) ? current_objective_value
: new_objective_value;
if (new_best_bound > current_best_bound) { // minimization.
LogNewSolution("ObjLb", wall_timer->Get(), new_best_bound,
current_objective_value, worker_info->worker_name);
reported_objective_value,
worker_info->worker_name);
} else {
LogNewSolution("ObjUb", wall_timer->Get(),
current_objective_value, new_best_bound,
reported_objective_value, new_best_bound,
worker_info->worker_name);
}
}

9
ortools/util/sorted_interval_list.cc
View File

@@ -136,6 +136,15 @@ Domain Domain::FromIntervals(absl::Span<const ClosedInterval> intervals) {
bool Domain::IsEmpty() const { return intervals_.empty(); }
int64 Domain::Size() const {
int64 size = 0;
for (const ClosedInterval interval : intervals_) {
size += operations_research::CapAdd(
1, operations_research::CapSub(interval.end, interval.start));
}
return size;
}
int64 Domain::Min() const {
CHECK(!IsEmpty());
return intervals_.front().start;

3
ortools/util/sorted_interval_list.h
View File

@@ -87,6 +87,9 @@ class Domain {
// Returns true if this is the empty set.
bool IsEmpty() const;
// Returns the number of elements in the domain. It is capped at kint64max.
int64 Size() const;
// Returns the domain min/max value.
// This Checks that the domain is not empty.
int64 Min() const;