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

continue working on cuts; fix bug in CP-SAT Presolve (IntProd)

Browse Source
This commit is contained in:
Laurent Perron
2019-01-11 15:31:01 +01:00
parent 179d31baf3
commit c4cf58701f
8 changed files with 85 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
ortools/sat/cp_model_presolve.cc
View File

@@ -688,19 +688,31 @@ bool PresolveIntProd(ConstraintProto* ct, PresolveContext* context) {
if (ct->int_prod().vars_size() == 2) {
int a = ct->int_prod().vars(0);
int b = ct->int_prod().vars(1);
const int p = ct->int_prod().target();
const int product = ct->int_prod().target();
if (context->IsFixed(b)) std::swap(a, b);
if (context->IsFixed(a)) {
ConstraintProto* const lin = context->working_model->add_constraints();
lin->mutable_linear()->add_vars(b);
lin->mutable_linear()->add_coeffs(context->MinOf(a));
lin->mutable_linear()->add_vars(p);
lin->mutable_linear()->add_coeffs(-1);
lin->mutable_linear()->add_domain(0);
lin->mutable_linear()->add_domain(0);
if (b != product) {
ConstraintProto* const lin = context->working_model->add_constraints();
lin->mutable_linear()->add_vars(b);
lin->mutable_linear()->add_coeffs(context->MinOf(a));
lin->mutable_linear()->add_vars(product);
lin->mutable_linear()->add_coeffs(-1);
lin->mutable_linear()->add_domain(0);
lin->mutable_linear()->add_domain(0);
context->UpdateRuleStats("int_prod: linearize product by constant.");
context->UpdateRuleStats("int_prod: linearize product by constant.");
return RemoveConstraint(ct, context);
} else if (context->MinOf(a) != 1) {
context->IntersectDomainWith(product, Domain(0, 0));
context->UpdateRuleStats("int_prod: fix variable to zero.");
} else {
context->UpdateRuleStats("int_prod: remove identity.");
return RemoveConstraint(ct, context);
}
} else if (a == b && a == product) { // x = x * x, only true for {0, 1}.
context->IntersectDomainWith(product, Domain(0, 1));
context->UpdateRuleStats("int_prod: fix variable to zero or one.");
return RemoveConstraint(ct, context);
}
}

2
ortools/sat/cp_model_solver.cc
View File

@@ -2079,7 +2079,7 @@ CpSolverResponse SolveCpModelParallel(
}
VLOG(1) << "Starting parallel search with " << num_search_workers
<< " workers: " << absl::StrJoin(worker_names, ", ");
<< " workers: [" << absl::StrJoin(worker_names, ", ") << "]";
if (!model_proto.has_objective()) {
{

36
ortools/sat/cuts.cc
View File

@@ -590,7 +590,7 @@ std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
const IntegerValue t = std::max(
IntegerValue(1),
remainder == 0 ? IntegerValue(1)
: std::min(max_scaling, CeilRatio(target, remainder)));
: 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);
@@ -599,7 +599,33 @@ std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
};
}
void IntegerRoundingCut(std::vector<double> lp_values,
// 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);
};
} 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));
};
}
}
void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
std::vector<IntegerValue> lower_bounds,
std::vector<IntegerValue> upper_bounds,
LinearConstraint* cut) {
@@ -730,9 +756,11 @@ void IntegerRoundingCut(std::vector<double> lp_values,
// Create the super-additive function f().
const IntegerValue rhs_remainder =
cut->ub - FloorRatio(cut->ub, divisor) * divisor;
const IntegerValue kMaxScaling(4);
const auto f =
GetSuperAdditiveRoundingFunction(rhs_remainder, divisor, kMaxScaling);
options.use_mir
? GetMirFunction(rhs_remainder, divisor, options.max_scaling)
: GetSuperAdditiveRoundingFunction(rhs_remainder, divisor,
options.max_scaling);
// Apply f() to the cut.
cut->ub = f(cut->ub);

12
ortools/sat/cuts.h
View File

@@ -65,7 +65,11 @@ struct CutGenerator {
//
// TODO(user): Still not clear to me if this can be improved or not.
std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
IntegerValue remainder, IntegerValue divisor, IntegerValue max_scaling);
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);
// 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.
@@ -96,7 +100,11 @@ std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
// more effort tunning them. In particular, one can try many heuristics and keep
// the best looking cut (or more than one). This is not on the critical code
// path, so we can spend more effort in finding good cuts.
void IntegerRoundingCut(std::vector<double> lp_values,
struct RoundingOptions {
bool use_mir = false;
IntegerValue max_scaling = IntegerValue(1000);
};
void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
std::vector<IntegerValue> lower_bounds,
std::vector<IntegerValue> upper_bounds,
LinearConstraint* cut);

14
ortools/sat/linear_constraint_manager.cc
View File

@@ -142,17 +142,23 @@ void LinearConstraintManager::AddCut(
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]));
activity += ToDouble(ct.coeffs[i]) * lp_solution[ct.vars[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
<< " violation: " << violation;
LOG(INFO) << "Cut '" << type_name << "'"
<< " size=" << size << " max_magnitude=" << max_magnitude
<< " norm=" << l2_norm << " violation=" << violation
<< " eff=" << violation / l2_norm;
}
}

4
ortools/sat/linear_programming_constraint.cc
View File

@@ -502,7 +502,9 @@ void LinearProgrammingConstraint::AddCGCuts() {
}
// Get the cut using some integer rounding heuristic.
IntegerRoundingCut(lp_values, var_lbs, var_ubs, &cut);
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

4
ortools/sat/lp_utils.cc
View File

@@ -201,7 +201,7 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
// on an int64, if the scaled bounds are too large, the constraint is either
// always true or always false.
Fractional lb = mp_constraint.lower_bound();
lb -= std::max(1.0, std::abs(lb)) * 1e-6;
lb -= std::max(1.0, std::abs(lb)) * 1e-8;
const Fractional scaled_lb = std::ceil(lb * scaling_factor);
if (lb == -kInfinity || scaled_lb <= kint64min) {
arg->add_domain(kint64min);
@@ -211,7 +211,7 @@ bool ConvertMPModelProtoToCpModelProto(const MPModelProto& mp_model,
.value());
}
Fractional ub = mp_constraint.upper_bound();
ub += std::max(1.0, std::abs(ub)) * 1e-6;
ub += std::max(1.0, std::abs(ub)) * 1e-8;
const Fractional scaled_ub = std::floor(ub * scaling_factor);
if (ub == kInfinity || scaled_ub >= kint64max) {
arg->add_domain(kint64max);

7
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: 118
// NEXT TAG: 119
message SatParameters {
// ==========================================================================
// Branching and polarity
@@ -507,6 +507,11 @@ message SatParameters {
// of enum that list enabled cuts, like "cuts:[KNAPSACK,GOMORY,CIRCUIT]".
optional bool add_cg_cuts = 117 [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];
// 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