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

add initial implementation of cover cuts

Browse Source
This commit is contained in:
Laurent Perron
2020-09-03 17:50:04 +02:00
parent a042fb2cae
commit 2bb2115cdc
9 changed files with 241 additions and 65 deletions
Download Patch File Download Diff File Expand all files Collapse all files

108
ortools/sat/cuts.cc
View File

@@ -1137,6 +1137,114 @@ void IntegerRoundingCutHelper::ComputeCut(
DivideByGCD(cut);
}
bool CoverCutHelper::TrySimpleKnapsack(
const LinearConstraint base_ct, const std::vector<double>& lp_values,
const std::vector<IntegerValue>& lower_bounds,
const std::vector<IntegerValue>& upper_bounds) {
const int base_size = lp_values.size();
// Fill terms with a rewrite of the base constraint where all coeffs &
// variables are positive by using either (X - LB) or (UB - X) as new
// variables.
terms_.clear();
IntegerValue rhs = base_ct.ub;
IntegerValue sum_of_diff(0);
for (int i = 0; i < base_size; ++i) {
const IntegerValue coeff = base_ct.coeffs[i];
const IntegerValue bound_diff = upper_bounds[i] - lower_bounds[i];
if (!AddProductTo(IntTypeAbs(coeff), bound_diff, &sum_of_diff)) {
return false;
}
const IntegerValue diff = IntTypeAbs(coeff) * bound_diff;
if (coeff > 0) {
if (!AddProductTo(-coeff, lower_bounds[i], &rhs)) return false;
terms_.push_back({i, ToDouble(upper_bounds[i]) - lp_values[i], diff});
} else {
if (!AddProductTo(-coeff, upper_bounds[i], &rhs)) return false;
terms_.push_back({i, lp_values[i] - ToDouble(lower_bounds[i]), diff});
}
}
// Try simple cover heuristic.
// Look for violated CUT if the form sum (UB - X) or (X - LB) >= 1.
double activity = 0.0;
int new_size = 0;
std::sort(terms_.begin(), terms_.end(), [](const Term& a, const Term& b) {
return a.dist_to_max_value < b.dist_to_max_value;
});
for (int i = 0; i < terms_.size(); ++i) {
const Term& term = terms_[i];
activity += term.dist_to_max_value;
if (activity > 0.99) {
new_size = i; // before this entry.
break;
}
rhs -= term.diff;
}
// If the rhs is now negative, we have a cut.
if (rhs >= 0) return false;
if (new_size == 0) return false;
// Transfrom to a minimal cover. We want to greedily remove the largest coeff
// first, so we have more chance for the "lifting" below which can increase
// the cut violation.
terms_.resize(new_size);
std::sort(terms_.begin(), terms_.end(),
[&base_ct](const Term& a, const Term& b) {
return IntTypeAbs(base_ct.coeffs[a.index]) >
IntTypeAbs(base_ct.coeffs[b.index]);
});
in_cut_.assign(base_ct.vars.size(), false);
// Remove terms greedily to get a minimal cover.
// Compute the cut at the same time.
cut_.ClearTerms();
cut_.lb = kMinIntegerValue;
cut_.ub = IntegerValue(-1);
IntegerValue max_coeff(0);
for (const Term term : terms_) {
if (term.diff + rhs < 0) {
rhs += term.diff;
continue;
}
in_cut_[term.index] = true;
max_coeff = std::max(max_coeff, IntTypeAbs(base_ct.coeffs[term.index]));
cut_.vars.push_back(base_ct.vars[term.index]);
if (base_ct.coeffs[term.index] > 0) {
cut_.coeffs.push_back(IntegerValue(1));
cut_.ub += upper_bounds[term.index];
} else {
cut_.coeffs.push_back(IntegerValue(-1));
cut_.ub -= lower_bounds[term.index];
}
}
// Basic Lifting. A move in activity of 1 in the cut correspond to a move of
// max_coeff in the original constraint.
num_lifting_ = 0;
for (int i = 0; i < base_size; ++i) {
if (in_cut_[i]) continue;
const IntegerValue magnitude = IntTypeAbs(base_ct.coeffs[i]);
if (magnitude < max_coeff) continue;
++num_lifting_;
const IntegerValue f = FloorRatio(magnitude, max_coeff);
if (base_ct.coeffs[i] > 0) { // Add f * (X - LB)
cut_.coeffs.push_back(IntegerValue(f));
cut_.vars.push_back(base_ct.vars[i]);
cut_.ub += lower_bounds[i] * f;
} else { // Add f * (UB - X)
cut_.coeffs.push_back(IntegerValue(-f));
cut_.vars.push_back(base_ct.vars[i]);
cut_.ub -= upper_bounds[i] * f;
}
}
return true;
}
CutGenerator CreatePositiveMultiplicationCutGenerator(IntegerVariable z,
IntegerVariable x,
IntegerVariable y,

30
ortools/sat/cuts.h
View File

@@ -226,6 +226,36 @@ class IntegerRoundingCutHelper {
std::vector<std::pair<IntegerVariable, IntegerValue>> tmp_terms_;
};
// Helper to find knapsack or flow cover cuts (not yet implemented).
class CoverCutHelper {
public:
// Try to find a cut with a knapsack heuristic.
// If this returns true, you can get the cut via cut().
bool TrySimpleKnapsack(const LinearConstraint base_ct,
const std::vector<double>& lp_values,
const std::vector<IntegerValue>& lower_bounds,
const std::vector<IntegerValue>& upper_bounds);
// If successful, info about the last generated cut.
LinearConstraint* mutable_cut() { return &cut_; }
const LinearConstraint& cut() const { return cut_; }
// Single line of text that we append to the cut log line.
const std::string Info() { return absl::StrCat("lift=", num_lifting_); }
private:
struct Term {
int index;
double dist_to_max_value;
IntegerValue diff;
};
std::vector<Term> terms_;
std::vector<bool> in_cut_;
LinearConstraint cut_;
int num_lifting_;
};
// If a variable is away from its upper bound by more than value 1.0, then it
// cannot be part of a cover that will violate the lp solution. This method
// returns a reduced constraint by removing such variables from the given

5
ortools/sat/linear_constraint.h
View File

@@ -50,6 +50,11 @@ struct LinearConstraint {
coeffs.push_back(coeff);
}
void ClearTerms() {
vars.clear();
coeffs.clear();
}
std::string DebugString() const {
std::string result;
if (lb.value() > kMinIntegerValue) {

75
ortools/sat/linear_programming_constraint.cc
View File

@@ -701,27 +701,53 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
}
}
// Get the cut using some integer rounding heuristic.
RoundingOptions options;
options.max_scaling = sat_parameters_.max_integer_rounding_scaling();
integer_rounding_cut_helper_.ComputeCut(options, tmp_lp_values_, tmp_var_lbs_,
tmp_var_ubs_,
&implied_bounds_processor_, &cut_);
bool at_least_one_added = false;
// Try cover appraoch to find cut.
{
if (cover_cut_helper_.TrySimpleKnapsack(cut_, tmp_lp_values_, tmp_var_lbs_,
tmp_var_ubs_)) {
at_least_one_added |= PostprocessAndAddCut(
absl::StrCat(name, "_K"), cover_cut_helper_.Info(), first_new_var,
first_slack, ib_slack_infos, cover_cut_helper_.mutable_cut());
}
}
// Try integer rounding heuristic to find cut.
{
RoundingOptions options;
options.max_scaling = sat_parameters_.max_integer_rounding_scaling();
integer_rounding_cut_helper_.ComputeCut(options, tmp_lp_values_,
tmp_var_lbs_, tmp_var_ubs_,
&implied_bounds_processor_, &cut_);
at_least_one_added |= PostprocessAndAddCut(
name,
absl::StrCat("num_lifted_booleans=",
integer_rounding_cut_helper_.NumLiftedBooleans()),
first_new_var, first_slack, ib_slack_infos, &cut_);
}
return at_least_one_added;
}
bool LinearProgrammingConstraint::PostprocessAndAddCut(
const std::string& name, const std::string& info,
IntegerVariable first_new_var, IntegerVariable first_slack,
const std::vector<ImpliedBoundsProcessor::SlackInfo>& ib_slack_infos,
LinearConstraint* cut) {
// Compute the activity. Warning: the cut no longer have the same size so we
// cannot use tmp_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_new_var) {
for (int i = 0; i < cut->vars.size(); ++i) {
if (cut->vars[i] < first_new_var) {
activity +=
ToDouble(cut_.coeffs[i]) * expanded_lp_solution_[cut_.vars[i]];
ToDouble(cut->coeffs[i]) * expanded_lp_solution_[cut->vars[i]];
}
}
const double kMinViolation = 1e-4;
const double violation = activity - ToDouble(cut_.ub);
const double violation = activity - ToDouble(cut->ub);
if (violation < kMinViolation) {
VLOG(3) << "Bad cut " << activity << " <= " << ToDouble(cut_.ub);
VLOG(3) << "Bad cut " << activity << " <= " << ToDouble(cut->ub);
return false;
}
@@ -729,26 +755,26 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
{
int num_slack = 0;
tmp_dense_vector_.assign(integer_variables_.size(), IntegerValue(0));
IntegerValue cut_ub = cut_.ub;
IntegerValue cut_ub = cut->ub;
bool overflow = false;
for (int i = 0; i < cut_.vars.size(); ++i) {
const IntegerVariable var = cut_.vars[i];
for (int i = 0; i < cut->vars.size(); ++i) {
const IntegerVariable var = cut->vars[i];
// Simple copy for non-slack variables.
if (var < first_new_var) {
const glop::ColIndex col =
gtl::FindOrDie(mirror_lp_variable_, PositiveVariable(var));
if (VariableIsPositive(var)) {
tmp_dense_vector_[col] += cut_.coeffs[i];
tmp_dense_vector_[col] += cut->coeffs[i];
} else {
tmp_dense_vector_[col] -= cut_.coeffs[i];
tmp_dense_vector_[col] -= cut->coeffs[i];
}
continue;
}
// Replace slack from bound substitution.
if (var < first_slack) {
const IntegerValue multiplier = cut_.coeffs[i];
const IntegerValue multiplier = cut->coeffs[i];
const int index = (var.value() - first_new_var.value()) / 2;
CHECK_LT(index, ib_slack_infos.size());
@@ -776,7 +802,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
++num_slack;
const int slack_index = (var.value() - first_slack.value()) / 2;
const glop::RowIndex row = tmp_slack_rows_[slack_index];
const IntegerValue multiplier = -cut_.coeffs[i];
const IntegerValue multiplier = -cut->coeffs[i];
if (!AddLinearExpressionMultiple(multiplier, integer_lp_[row].terms,
&tmp_dense_vector_)) {
overflow = true;
@@ -796,23 +822,22 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
}
VLOG(3) << " num_slack: " << num_slack;
ConvertToLinearConstraint(tmp_dense_vector_, cut_ub, &cut_);
ConvertToLinearConstraint(tmp_dense_vector_, cut_ub, cut);
}
// Display some stats used for investigation of cut generation.
const std::string extra_info = absl::StrCat(
"num_ib_substitutions=", ib_slack_infos.size(), " num_lifted_booleans=",
integer_rounding_cut_helper_.NumLiftedBooleans());
const std::string extra_info =
absl::StrCat(info, " num_ib_substitutions=", ib_slack_infos.size());
const double new_violation =
ComputeActivity(cut_, expanded_lp_solution_) - ToDouble(cut_.ub);
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_);
return constraint_manager_.AddCut(cut_, name, expanded_lp_solution_,
DivideByGCD(cut);
return constraint_manager_.AddCut(*cut, name, expanded_lp_solution_,
extra_info);
}

8
ortools/sat/linear_programming_constraint.h
View File

@@ -197,6 +197,13 @@ class LinearProgrammingConstraint : public PropagatorInterface,
const std::vector<std::pair<glop::RowIndex, IntegerValue>>&
integer_multipliers);
// Second half of AddCutFromConstraints().
bool PostprocessAndAddCut(
const std::string& name, const std::string& info,
IntegerVariable first_new_var, IntegerVariable first_slack,
const std::vector<ImpliedBoundsProcessor::SlackInfo>& ib_slack_infos,
LinearConstraint* cut);
// Computes and adds the corresponding type of cuts.
// This can currently only be called at the root node.
void AddCGCuts();
@@ -341,6 +348,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
// Temporary data for cuts.
ZeroHalfCutHelper zero_half_cut_helper_;
CoverCutHelper cover_cut_helper_;
IntegerRoundingCutHelper integer_rounding_cut_helper_;
LinearConstraint cut_;
gtl::ITIVector<glop::ColIndex, IntegerValue> tmp_dense_vector_;

31
ortools/sat/optimization.cc
View File

@@ -215,37 +215,6 @@ class FuMalikSymmetryBreaker {
} // namespace
void MinimizeCore(SatSolver* solver, std::vector<Literal>* core) {
std::vector<Literal> temp = *core;
std::reverse(temp.begin(), temp.end());
solver->Backtrack(0);
solver->SetAssumptionLevel(0);
// Note that this Solve() is really fast, since the solver should detect that
// the assumptions are unsat with unit propagation only. This is just a
// convenient way to remove assumptions that are propagated by the one before
// them.
const SatSolver::Status status =
solver->ResetAndSolveWithGivenAssumptions(temp);
if (status != SatSolver::ASSUMPTIONS_UNSAT) {
if (status != SatSolver::LIMIT_REACHED) {
CHECK_NE(status, SatSolver::FEASIBLE);
// This should almost never happen, but it is not impossible. The reason
// is that the solver may delete some learned clauses required by the unit
// propagation to show that the core is unsat.
LOG(WARNING) << "This should only happen rarely! otherwise, investigate. "
<< "Returned status is " << SatStatusString(status);
}
return;
}
temp = solver->GetLastIncompatibleDecisions();
if (temp.size() < core->size()) {
VLOG(1) << "minimization " << core->size() << " -> " << temp.size();
std::reverse(temp.begin(), temp.end());
*core = temp;
}
}
void MinimizeCoreWithPropagation(SatSolver* solver,
std::vector<Literal>* core) {
if (solver->IsModelUnsat()) return;

9
ortools/sat/optimization.h
View File

@@ -28,15 +28,6 @@
namespace operations_research {
namespace sat {
// Tries to minimize the given UNSAT core with a really simple heuristic.
// The idea is to remove literals that are consequences of others in the core.
// We already know that in the initial order, no literal is propagated by the
// one before it, so we just look for propagation in the reverse order.
//
// Important: The given SatSolver must be the one that just produced the given
// core.
void MinimizeCore(SatSolver* solver, std::vector<Literal>* core);
// Like MinimizeCore() with a slower but strictly better heuristic. This
// algorithm should produce a minimal core with respect to propagation. We put
// each literal of the initial core "last" at least once, so if such literal can

31
ortools/sat/sat_solver.cc
View File

@@ -2539,5 +2539,36 @@ std::string SatStatusString(SatSolver::Status status) {
return "UNKNOWN";
}
void MinimizeCore(SatSolver* solver, std::vector<Literal>* core) {
std::vector<Literal> temp = *core;
std::reverse(temp.begin(), temp.end());
solver->Backtrack(0);
solver->SetAssumptionLevel(0);
// Note that this Solve() is really fast, since the solver should detect that
// the assumptions are unsat with unit propagation only. This is just a
// convenient way to remove assumptions that are propagated by the one before
// them.
const SatSolver::Status status =
solver->ResetAndSolveWithGivenAssumptions(temp);
if (status != SatSolver::ASSUMPTIONS_UNSAT) {
if (status != SatSolver::LIMIT_REACHED) {
CHECK_NE(status, SatSolver::FEASIBLE);
// This should almost never happen, but it is not impossible. The reason
// is that the solver may delete some learned clauses required by the unit
// propagation to show that the core is unsat.
LOG(WARNING) << "This should only happen rarely! otherwise, investigate. "
<< "Returned status is " << SatStatusString(status);
}
return;
}
temp = solver->GetLastIncompatibleDecisions();
if (temp.size() < core->size()) {
VLOG(1) << "minimization " << core->size() << " -> " << temp.size();
std::reverse(temp.begin(), temp.end());
*core = temp;
}
}
} // namespace sat
} // namespace operations_research

9
ortools/sat/sat_solver.h
View File

@@ -823,6 +823,15 @@ class SatSolver {
DISALLOW_COPY_AND_ASSIGN(SatSolver);
};
// Tries to minimize the given UNSAT core with a really simple heuristic.
// The idea is to remove literals that are consequences of others in the core.
// We already know that in the initial order, no literal is propagated by the
// one before it, so we just look for propagation in the reverse order.
//
// Important: The given SatSolver must be the one that just produced the given
// core.
void MinimizeCore(SatSolver* solver, std::vector<Literal>* core);
// ============================================================================
// Model based functions.
//