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

[CP-SAT] improve cut management

Browse Source
This commit is contained in:
Laurent Perron
2023-05-10 17:19:43 +02:00
parent 75c7eaa127
commit ff3bbcc0b3
5 changed files with 165 additions and 27 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
ortools/sat/cuts.cc
View File

@@ -1117,15 +1117,15 @@ int CoverCutHelper::GetCoverSizeForBooleans(int relevant_size) {
// Sorting can be slow, so we start by splitting the vector in 3 parts
// [can always be in cover, candidates, can never be in cover].
int part1 = 0;
const double threshold = 1.0 / static_cast<double>(relevant_size);
const double threshold = 1.0 - 1.0 / static_cast<double>(relevant_size);
for (int i = 0; i < relevant_size;) {
const double dist = base_ct_.terms[i].LpDistToMaxValue();
if (dist < threshold) {
const double lp_value = base_ct_.terms[i].lp_value;
if (lp_value >= threshold) {
// Move to part 1.
std::swap(base_ct_.terms[i], base_ct_.terms[part1]);
++i;
++part1;
} else if (dist < 0.9999) {
} else if (lp_value >= 0.001) {
// Keep in part 2.
++i;
} else {
@@ -1134,16 +1134,16 @@ int CoverCutHelper::GetCoverSizeForBooleans(int relevant_size) {
std::swap(base_ct_.terms[i], base_ct_.terms[relevant_size]);
}
}
// Sort by decreasing Lp value.
std::sort(base_ct_.terms.begin() + part1,
base_ct_.terms.begin() + relevant_size,
[](const CutTerm& a, const CutTerm& b) {
const double dist_a = a.LpDistToMaxValue();
const double dist_b = b.LpDistToMaxValue();
if (dist_a == dist_b) {
if (a.lp_value == b.lp_value) {
// Prefer low coefficients if the distance is the same.
return a.coeff < b.coeff;
}
return dist_a < dist_b;
return a.lp_value > b.lp_value;
});
double activity = 0.0;

2
ortools/sat/linear_constraint_manager.cc
View File

@@ -276,7 +276,7 @@ bool LinearConstraintManager::AddCut(
const double l2_norm = ComputeL2Norm(ct);
// Only add cut with sufficient efficacy.
if (violation / l2_norm < 1e-5) {
if (violation / l2_norm < 1e-4) {
VLOG(3) << "BAD Cut '" << type_name << "'"
<< " size=" << ct.vars.size()
<< " max_magnitude=" << ComputeInfinityNorm(ct)

107
ortools/sat/linear_programming_constraint.cc
View File

@@ -165,6 +165,35 @@ void ScatteredIntegerVector::ConvertToLinearConstraint(
result->ub = upper_bound;
}
void ScatteredIntegerVector::ConvertToCutData(
absl::int128 rhs, const std::vector<IntegerVariable>& integer_variables,
const std::vector<double>& lp_solution, IntegerTrail* integer_trail,
CutData* result) {
result->terms.clear();
result->rhs = rhs;
if (is_sparse_) {
std::sort(non_zeros_.begin(), non_zeros_.end());
for (const glop::ColIndex col : non_zeros_) {
const IntegerValue coeff = dense_vector_[col];
if (coeff == 0) continue;
const IntegerVariable var = integer_variables[col.value()];
CHECK(result->AppendOneTerm(var, coeff, lp_solution[col.value()],
integer_trail->LevelZeroLowerBound(var),
integer_trail->LevelZeroUpperBound(var)));
}
} else {
const int size = dense_vector_.size();
for (glop::ColIndex col(0); col < size; ++col) {
const IntegerValue coeff = dense_vector_[col];
if (coeff == 0) continue;
const IntegerVariable var = integer_variables[col.value()];
CHECK(result->AppendOneTerm(var, coeff, lp_solution[col.value()],
integer_trail->LevelZeroLowerBound(var),
integer_trail->LevelZeroUpperBound(var)));
}
}
}
std::vector<std::pair<glop::ColIndex, IntegerValue>>
ScatteredIntegerVector::GetTerms() {
std::vector<std::pair<glop::ColIndex, IntegerValue>> result;
@@ -356,12 +385,14 @@ bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
for (const LinearConstraintInternal& ct : integer_lp_) {
const ConstraintIndex row = lp_data_.CreateNewConstraint();
// Note that we use the bounds even if they are trivial.
//
// TODO(user): This seems good for things like sum bool <= 1 since setting
// the slack in [0, 1] can lead to bound flip in the simplex. However if the
// bound is large, maybe it make more sense to use +/- infinity.
lp_data_.SetConstraintBounds(row, ToDouble(ct.lb), ToDouble(ct.ub));
// TODO(user): Using trivial bound might be good for things like
// sum bool <= 1 since setting the slack in [0, 1] can lead to bound flip in
// the simplex. However if the bound is large, maybe it make more sense to
// use +/- infinity.
const double infinity = std::numeric_limits<double>::infinity();
lp_data_.SetConstraintBounds(
row, ct.lb_is_trivial ? -infinity : ToDouble(ct.lb),
ct.ub_is_trivial ? +infinity : ToDouble(ct.ub));
for (const auto& term : ct.terms) {
lp_data_.SetCoefficient(row, term.first, ToDouble(term.second));
}
@@ -831,7 +862,12 @@ bool LinearProgrammingConstraint::AnalyzeLp() {
return true;
}
// TODO(user): We could also strengthen the coefficients.
// Note that since we call this on the constraint with slack, we actually have
// linear expression == rhs, we can use this to propagate more!
//
// TODO(user): Also propagate on -cut ? in practice we already do that in many
// places were we try to generate the cut on -cut... But we coould do it sooner
// and more cleanly here.
bool LinearProgrammingConstraint::PreprocessCut(IntegerVariable first_slack,
CutData* cut) {
// Because of complement, all coeffs and all terms are positive after this.
@@ -841,15 +877,47 @@ bool LinearProgrammingConstraint::PreprocessCut(IntegerVariable first_slack,
return false;
}
// Limited DP to compute first few reachable values.
// Note that all coeff are positive.
reachable_.Reset();
for (const CutTerm& term : cut->terms) {
reachable_.Add(term.coeff.value());
}
// Extra propag since we know it is actually an equality constraint.
if (cut->rhs < absl::int128(reachable_.LastValue()) &&
!reachable_.MightBeReachable(static_cast<int64_t>(cut->rhs))) {
problem_proven_infeasible_by_cuts_ = true;
return false;
}
bool some_fixed_terms = false;
bool some_relevant_positions = false;
for (CutTerm& term : cut->terms) {
const absl::int128 magnitude128 = term.coeff.value();
const absl::int128 range =
absl::int128(term.bound_diff.value()) * magnitude128;
IntegerValue new_diff = term.bound_diff;
if (range > cut->rhs) {
// Update the term bound_diff.
term.bound_diff = static_cast<int64_t>(cut->rhs / magnitude128);
new_diff = static_cast<int64_t>(cut->rhs / magnitude128);
}
{
// Extra propag since we know it is actually an equality constraint.
absl::int128 rest128 =
cut->rhs - absl::int128(new_diff.value()) * magnitude128;
while (rest128 < absl::int128(reachable_.LastValue()) &&
!reachable_.MightBeReachable(static_cast<int64_t>(rest128))) {
++total_num_eq_propagations_;
CHECK_GT(new_diff, 0);
--new_diff;
rest128 += magnitude128;
}
}
if (new_diff < term.bound_diff) {
term.bound_diff = new_diff;
const IntegerVariable var = term.expr_vars[0];
if (var < first_slack) {
// Normal variable.
@@ -940,19 +1008,16 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
return false;
}
// Important: because we use integer_multipliers below, we cannot just
// divide by GCD or call PreventOverflow() here.
// Important: because we use integer_multipliers below to create the slack,
// we cannot try to adjust this linear combination easily.
//
// TODO(user): the conversion col_index -> IntegerVariable is slow and could
// in principle be removed. Easy for cuts, but not so much for
// implied_bounds_processor_. Note that in theory this could allow us to
// use Literal directly without the need to have an IntegerVariable for them.
tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, cut_ub,
&cut_);
// TODO(user): Remove one conversion step.
CHECK(base_ct_.FillFromLinearConstraint(cut_, expanded_lp_solution_,
integer_trail_));
tmp_scattered_vector_.ConvertToCutData(cut_ub.value(), integer_variables_,
lp_solution_, integer_trail_,
&base_ct_);
// If there are no integer (all Booleans), no need to try implied bounds
// heurititics. By setting this to nullptr, we are a bit faster.
@@ -1006,6 +1071,9 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
// This also make all coefficients positive.
if (!PreprocessCut(first_slack, &base_ct_)) return false;
// We cannot cut sum Bool <= 1.
if (base_ct_.rhs == 1) return false;
// Constraint with slack should be tight.
if (DEBUG_MODE) {
double activity = 0.0;
@@ -1033,6 +1101,7 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
}
// Try cover approach to find cut.
// TODO(user): Share common computation between two kinds.
{
DCHECK(base_ct_.AllCoefficientsArePositive());
@@ -1206,7 +1275,7 @@ void LinearProgrammingConstraint::AddCGCuts() {
tmp_integer_multipliers_ = ScaleLpMultiplier(
/*take_objective_into_account=*/false,
/*ignore_trivial_constraints=*/old_gomory, tmp_lp_multipliers_,
&scaling, int64_t{1} << 52);
&scaling);
if (scaling != 0) {
AddCutFromConstraints("CG", tmp_integer_multipliers_);
}
@@ -2633,6 +2702,8 @@ std::string LinearProgrammingConstraint::Statistics() const {
FormatCounter(total_num_simplex_iterations_), "\n");
absl::StrAppend(&result, " total num cut propagation: ",
FormatCounter(total_num_cut_propagations_), "\n");
absl::StrAppend(&result, " total num eq cut propagation: ",
FormatCounter(total_num_eq_propagations_), "\n");
absl::StrAppend(&result, " total num cut overflow: ",
FormatCounter(num_cut_overflows_), "\n");
absl::StrAppend(&result, " total num adjust: ", FormatCounter(num_adjusts_),

7
ortools/sat/linear_programming_constraint.h
View File

@@ -97,6 +97,11 @@ class ScatteredIntegerVector {
const std::vector<IntegerVariable>& integer_variables,
IntegerValue upper_bound, LinearConstraint* result);
void ConvertToCutData(absl::int128 rhs,
const std::vector<IntegerVariable>& integer_variables,
const std::vector<double>& lp_solution,
IntegerTrail* integer_trail, CutData* result);
// Similar to ConvertToLinearConstraint().
std::vector<std::pair<glop::ColIndex, IntegerValue>> GetTerms();
@@ -575,7 +580,9 @@ class LinearProgrammingConstraint : public PropagatorInterface,
// As we form candidate form cuts, sometimes we can propagate level zero
// bounds with them.
FirstFewValues<10> reachable_;
int64_t total_num_cut_propagations_ = 0;
int64_t total_num_eq_propagations_ = 0;
// Some stats on the LP statuses encountered.
int64_t num_solves_ = 0;

60
ortools/sat/util.h
View File

@@ -34,6 +34,7 @@
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_parameters.pb.h"
#include "ortools/util/random_engine.h"
#include "ortools/util/saturated_arithmetic.h"
#include "ortools/util/sorted_interval_list.h"
#include "ortools/util/time_limit.h"
@@ -236,6 +237,65 @@ class MaxBoundedSubsetSum {
std::vector<int64_t> filtered_values_;
};
// Simple DP to keep the set of the first n reachable value (n > 1).
//
// TODO(user): Maybe modulo some prime number we can keep more info.
// TODO(user): Another common case is a bunch of really small values and larger
// ones, so we could bound the sum of the small values and keep the first few
// reachable by the big ones. This is similar to some presolve transformations.
template <int n>
class FirstFewValues {
public:
FirstFewValues() { Reset(); }
void Reset() {
reachable_.fill(std::numeric_limits<int64_t>::max());
reachable_[0] = 0;
new_reachable_[0] = 0;
}
// We assume the given positive value can be added as many time as wanted.
//
// TODO(user): Implement Add() with an upper bound on the multiplicity.
void Add(const int64_t positive_value) {
DCHECK_GT(positive_value, 0);
if (positive_value >= reachable_.back()) return;
// We copy from reachable_[i] to new_reachable_[j].
// The position zero is already copied.
int i = 1;
int j = 1;
for (int base = 0; j < n && base < n; ++base) {
const int64_t candidate = CapAdd(new_reachable_[base], positive_value);
while (j < n && i < n && reachable_[i] < candidate) {
new_reachable_[j++] = reachable_[i++];
}
if (j < n) {
// Eliminate duplicates.
while (i < n && reachable_[i] == candidate) i++;
// insert candidate in its final place.
new_reachable_[j++] = candidate;
}
}
std::swap(reachable_, new_reachable_);
}
// Returns true iff sum might be expressible as a weighted sum of the added
// value. Any sum >= LastValue() is always considered potentially reachable.
bool MightBeReachable(int64_t sum) const {
if (sum >= reachable_.back()) return true;
return std::binary_search(reachable_.begin(), reachable_.end(), sum);
}
const std::array<int64_t, n>& reachable() const { return reachable_; }
int64_t LastValue() const { return reachable_.back(); }
private:
std::array<int64_t, n> reachable_;
std::array<int64_t, n> new_reachable_;
};
// Use Dynamic programming to solve a single knapsack. This is used by the
// presolver to simplify variables appearing in a single linear constraint.
//