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

add initial implementation for zero half cut

Browse Source
This commit is contained in:
Laurent Perron
2020-09-03 09:08:10 +02:00
parent 0b8f2f2b5b
commit c52946b783
7 changed files with 479 additions and 10 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
ortools/sat/BUILD
View File

@@ -898,6 +898,7 @@ cc_library(
":linear_constraint",
":linear_constraint_manager",
":model",
":zero_half_cuts",
"//ortools/base",
"//ortools/base:int_type",
"//ortools/base:int_type_indexed_vector",
@@ -951,6 +952,19 @@ cc_library(
],
)
cc_library(
name = "zero_half_cuts",
srcs = ["zero_half_cuts.cc"],
hdrs = ["zero_half_cuts.h"],
deps = [
":integer",
":util",
"//ortools/base",
"//ortools/base:int_type",
"//ortools/lp_data:base",
],
)
cc_library(
name = "lp_utils",
srcs = ["lp_utils.cc"],

8
ortools/sat/cp_model_presolve.cc
View File

@@ -2927,8 +2927,8 @@ bool CpModelPresolver::PresolveCircuit(ConstraintProto* ct) {
// here (even if they are at false).
std::vector<std::vector<int>> incoming_arcs;
std::vector<std::vector<int>> outgoing_arcs;
const int num_arcs = proto.literals_size();
int num_nodes = 0;
const int num_arcs = proto.literals_size();
for (int i = 0; i < num_arcs; ++i) {
const int ref = proto.literals(i);
const int tail = proto.tails(i);
@@ -4479,12 +4479,16 @@ void CpModelPresolver::PresolveToFixPoint() {
if (constraints.size() != 1) continue;
const int c = *constraints.begin();
if (c < 0) continue;
if (in_queue[c]) continue;
// Note that to avoid bad complexity in problem like a TSP with just one
// big constraint. we mark all the singleton variables of a constraint
// even if this constraint is already in the queue.
if (gtl::ContainsKey(var_constraint_pair_already_called,
std::pair<int, int>(v, c))) {
continue;
}
var_constraint_pair_already_called.insert({v, c});
if (!in_queue[c]) {
in_queue[c] = true;
queue.push_back(c);

40
ortools/sat/linear_programming_constraint.cc
View File

@@ -37,6 +37,7 @@
#include "ortools/lp_data/lp_types.h"
#include "ortools/sat/implied_bounds.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/zero_half_cuts.h"
#include "ortools/util/saturated_arithmetic.h"
namespace operations_research {
@@ -600,7 +601,8 @@ bool LinearProgrammingConstraint::AddCutFromConstraints(
// divide by GCD or call PreventOverflow() here.
ConvertToLinearConstraint(tmp_dense_vector_, cut_ub, &cut_);
// This should be tight.
// Note that the base constraint we use are currently always tight.
// It is not a requirement though.
const double norm = ToDouble(ComputeInfinityNorm(cut_));
if (std::abs(ComputeActivity(cut_, expanded_lp_solution_) -
ToDouble(cut_.ub)) /
@@ -1136,6 +1138,41 @@ void LinearProgrammingConstraint::AddMirCuts() {
}
}
void LinearProgrammingConstraint::AddZeroHalfCuts() {
CHECK_EQ(trail_->CurrentDecisionLevel(), 0);
tmp_lp_values_.clear();
tmp_var_lbs_.clear();
tmp_var_ubs_.clear();
for (const IntegerVariable var : integer_variables_) {
tmp_lp_values_.push_back(expanded_lp_solution_[var]);
tmp_var_lbs_.push_back(integer_trail_->LowerBound(var));
tmp_var_ubs_.push_back(integer_trail_->UpperBound(var));
}
// TODO(user): See if it make sense to try to use implied bounds there.
zero_half_cut_helper_.ProcessVariables(tmp_lp_values_, tmp_var_lbs_,
tmp_var_ubs_);
for (glop::RowIndex row(0); row < integer_lp_.size(); ++row) {
// Even though we could use non-tight row, for now we prefer to use tight
// ones.
const auto status = simplex_.GetConstraintStatus(row);
if (status == glop::ConstraintStatus::BASIC) continue;
if (status == glop::ConstraintStatus::FREE) continue;
zero_half_cut_helper_.AddOneConstraint(
row, integer_lp_[row].terms, integer_lp_[row].lb, integer_lp_[row].ub);
}
for (const std::vector<std::pair<RowIndex, IntegerValue>>& multipliers :
zero_half_cut_helper_.InterestingCandidates(random_)) {
// TODO(user): Make sure that if the resulting linear coefficients are not
// too high, we do try a "divisor" of two and thus try a true zero-half cut
// instead of just using our best MIR heuristic (which might still be better
// though).
AddCutFromConstraints("ZERO_HALF", multipliers);
}
}
void LinearProgrammingConstraint::UpdateSimplexIterationLimit(
const int64 min_iter, const int64 max_iter) {
if (sat_parameters_.linearization_level() < 2) return;
@@ -1227,6 +1264,7 @@ bool LinearProgrammingConstraint::Propagate() {
implied_bounds_processor_.ClearCache();
if (sat_parameters_.add_mir_cuts()) AddMirCuts();
if (sat_parameters_.add_cg_cuts()) AddCGCuts();
if (sat_parameters_.add_zero_half_cuts()) AddZeroHalfCuts();
}
// Try to add cuts.

8
ortools/sat/linear_programming_constraint.h
View File

@@ -32,6 +32,7 @@
#include "ortools/sat/linear_constraint_manager.h"
#include "ortools/sat/model.h"
#include "ortools/sat/util.h"
#include "ortools/sat/zero_half_cuts.h"
#include "ortools/util/rev.h"
#include "ortools/util/time_limit.h"
@@ -196,13 +197,11 @@ class LinearProgrammingConstraint : public PropagatorInterface,
const std::vector<std::pair<glop::RowIndex, IntegerValue>>&
integer_multipliers);
// Computes and adds Chvatal-Gomory cuts.
// Computes and adds the corresponding type of 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();
void AddZeroHalfCuts();
// Updates the bounds of the LP variables from the CP bounds.
void UpdateBoundsOfLpVariables();
@@ -341,6 +340,7 @@ class LinearProgrammingConstraint : public PropagatorInterface,
glop::LpScalingHelper scaler_;
// Temporary data for cuts.
ZeroHalfCutHelper zero_half_cut_helper_;
IntegerRoundingCutHelper integer_rounding_cut_helper_;
LinearConstraint cut_;
gtl::ITIVector<glop::ColIndex, IntegerValue> tmp_dense_vector_;

10
ortools/sat/sat_parameters.proto
View File

@@ -21,7 +21,7 @@ option java_multiple_files = true;
// Contains the definitions for all the sat algorithm parameters and their
// default values.
//
// NEXT TAG: 169
// NEXT TAG: 170
message SatParameters {
// ==========================================================================
// Branching and polarity
@@ -588,13 +588,17 @@ message SatParameters {
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.
// Note that for now, this is not heavily tuned.
optional bool add_cg_cuts = 117 [default = true];
// Whether we generate MIR cuts at root node.
// Note that for now, this is not heavily tunned.
// Note that for now, this is not heavily tuned.
optional bool add_mir_cuts = 120 [default = true];
// Whether we generate Zero-Half cuts at root node.
// Note that for now, this is not heavily tuned.
optional bool add_zero_half_cuts = 169 [default = false];
// Cut generator for all diffs can add too many cuts for large all_diff
// constraints. This parameter restricts the large all_diff constraints to
// have a cut generator.

277
ortools/sat/zero_half_cuts.cc Normal file
View File

@@ -0,0 +1,277 @@
// Copyright 2010-2018 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/sat/zero_half_cuts.h"
namespace operations_research {
namespace sat {
void ZeroHalfCutHelper::Reset(int size) {
rows_.clear();
shifted_lp_values_.clear();
bound_parity_.clear();
col_to_rows_.clear();
col_to_rows_.resize(size);
tmp_marked_.resize(size);
}
void ZeroHalfCutHelper::ProcessVariables(
const std::vector<double>& lp_values,
const std::vector<IntegerValue>& lower_bounds,
const std::vector<IntegerValue>& upper_bounds) {
Reset(lp_values.size());
// Shift all variables to their closest bound.
lp_values_ = lp_values;
for (int i = 0; i < lp_values.size(); ++i) {
const double lb_dist = lp_values[i] - ToDouble(lower_bounds[i]);
const double ub_dist = ToDouble(upper_bounds[i]) - lp_values_[i];
if (lb_dist < ub_dist) {
shifted_lp_values_.push_back(lb_dist);
bound_parity_.push_back(lower_bounds[i].value() & 1);
} else {
shifted_lp_values_.push_back(ub_dist);
bound_parity_.push_back(upper_bounds[i].value() & 1);
}
}
}
void ZeroHalfCutHelper::AddBinaryRow(const CombinationOfRows& binary_row) {
// No point pushing an all zero row with a zero rhs.
if (binary_row.cols.empty() && !binary_row.rhs_parity) return;
for (const int col : binary_row.cols) {
col_to_rows_[col].push_back(rows_.size());
}
rows_.push_back(binary_row);
}
void ZeroHalfCutHelper::AddOneConstraint(
const glop::RowIndex row,
const std::vector<std::pair<glop::ColIndex, IntegerValue>>& terms,
IntegerValue lb, IntegerValue ub) {
if (terms.size() > kMaxInputConstraintSize) return;
double activity = 0.0;
IntegerValue magnitude(0);
CombinationOfRows binary_row;
int rhs_adjust = 0;
for (const auto& term : terms) {
const int col = term.first.value();
activity += ToDouble(term.second) * lp_values_[col];
magnitude = std::max(magnitude, IntTypeAbs(term.second));
// Only consider odd coefficient.
if ((term.second.value() & 1) == 0) continue;
// Ignore column in the binary matrix if its lp value is almost zero.
if (shifted_lp_values_[col] > 1e-2) {
binary_row.cols.push_back(col);
}
// Because we work on the shifted variable, the rhs needs to be updated.
rhs_adjust ^= bound_parity_[col];
}
// We ignore constraint with large coefficient, since there is little chance
// to cancel them and because of that the efficacity of a generated cut will
// be limited.
if (magnitude > kMaxInputConstraintMagnitude) return;
// TODO(user): experiment with the best value. probably only tight rows are
// best? and we could use the basis status rather than recomputing the
// activity for that.
//
// TODO(user): Avoid adding duplicates and just randomly pick one. Note
// that we should also remove duplicate in a generic way.
const double tighteness_threshold = 1e-2;
if (ToDouble(ub) - activity < tighteness_threshold) {
binary_row.multipliers = {{row, IntegerValue(1)}};
binary_row.slack = ToDouble(ub) - activity;
binary_row.rhs_parity = (ub.value() & 1) ^ rhs_adjust;
AddBinaryRow(binary_row);
}
if (activity - ToDouble(lb) < tighteness_threshold) {
binary_row.multipliers = {{row, IntegerValue(-1)}};
binary_row.slack = activity - ToDouble(lb);
binary_row.rhs_parity = (lb.value() & 1) ^ rhs_adjust;
AddBinaryRow(binary_row);
}
}
void ZeroHalfCutHelper::SymmetricDifference(
std::function<bool(int)> extra_condition, const std::vector<int>& a,
std::vector<int>* b) {
for (const int v : *b) tmp_marked_[v] = true;
for (const int v : a) {
if (tmp_marked_[v]) {
tmp_marked_[v] = false;
} else {
tmp_marked_[v] = true;
// TODO(user): optim by doing that at the end?
b->push_back(v);
}
}
// Remove position that are not marked, and clear tmp_marked_.
int new_size = 0;
for (const int v : *b) {
if (tmp_marked_[v]) {
if (extra_condition(v)) {
(*b)[new_size++] = v;
}
tmp_marked_[v] = false;
}
}
b->resize(new_size);
}
void ZeroHalfCutHelper::ProcessSingletonColumns() {
for (const int singleton_col : singleton_cols_) {
if (col_to_rows_[singleton_col].empty()) continue;
CHECK_EQ(col_to_rows_[singleton_col].size(), 1);
const int row = col_to_rows_[singleton_col][0];
int new_size = 0;
auto& mutable_cols = rows_[row].cols;
for (const int col : mutable_cols) {
if (col == singleton_col) continue;
mutable_cols[new_size++] = col;
}
CHECK_LT(new_size, mutable_cols.size());
mutable_cols.resize(new_size);
col_to_rows_[singleton_col].clear();
rows_[row].slack += shifted_lp_values_[singleton_col];
}
singleton_cols_.clear();
}
// This is basically one step of a Gaussian elimination with the given pivot.
void ZeroHalfCutHelper::EliminateVarUsingRow(int eliminated_col,
int eliminated_row) {
CHECK_LE(rows_[eliminated_row].slack, 1e-6);
CHECK(!rows_[eliminated_row].cols.empty());
// First update the row representation of the matrix.
tmp_marked_.resize(std::max(col_to_rows_.size(), rows_.size()));
DCHECK(std::all_of(tmp_marked_.begin(), tmp_marked_.end(),
[](bool b) { return !b; }));
int new_size = 0;
for (const int other_row : col_to_rows_[eliminated_col]) {
if (other_row == eliminated_row) continue;
col_to_rows_[eliminated_col][new_size++] = other_row;
SymmetricDifference([](int i) { return true; }, rows_[eliminated_row].cols,
&rows_[other_row].cols);
// Update slack & parity.
rows_[other_row].rhs_parity ^= rows_[eliminated_row].rhs_parity;
rows_[other_row].slack += rows_[eliminated_row].slack;
// Update the multipliers the same way.
{
auto& mutable_multipliers = rows_[other_row].multipliers;
mutable_multipliers.insert(mutable_multipliers.end(),
rows_[eliminated_row].multipliers.begin(),
rows_[eliminated_row].multipliers.end());
std::sort(mutable_multipliers.begin(), mutable_multipliers.end());
int new_size = 0;
for (const auto& entry : mutable_multipliers) {
if (new_size > 0 && entry == mutable_multipliers[new_size - 1]) {
// Cancel both.
--new_size;
} else {
mutable_multipliers[new_size++] = entry;
}
}
mutable_multipliers.resize(new_size);
}
}
col_to_rows_[eliminated_col].resize(new_size);
// Then update the col representation of the matrix.
//
// Note that we remove from the col-wise representation any rows with a large
// slack.
{
int new_size = 0;
for (const int other_col : rows_[eliminated_row].cols) {
if (other_col == eliminated_col) continue;
const int old_size = col_to_rows_[other_col].size();
rows_[eliminated_row].cols[new_size++] = other_col;
SymmetricDifference(
[this](int i) { return rows_[i].slack < kSlackThreshold; },
col_to_rows_[eliminated_col], &col_to_rows_[other_col]);
if (old_size != 1 && col_to_rows_[other_col].size() == 1) {
singleton_cols_.push_back(other_col);
}
}
rows_[eliminated_row].cols.resize(new_size);
}
// Clear col.
col_to_rows_[eliminated_col].clear();
rows_[eliminated_row].slack += shifted_lp_values_[eliminated_col];
}
std::vector<std::vector<std::pair<glop::RowIndex, IntegerValue>>>
ZeroHalfCutHelper::InterestingCandidates(ModelRandomGenerator* random) {
std::vector<std::vector<std::pair<glop::RowIndex, IntegerValue>>> result;
// Initialize singleton_cols_.
singleton_cols_.clear();
for (int col = 0; col < col_to_rows_.size(); ++col) {
if (col_to_rows_[col].size() == 1) singleton_cols_.push_back(col);
}
// Process rows by increasing size, but randomize if same size.
std::vector<int> to_process;
for (int row = 0; row < rows_.size(); ++row) to_process.push_back(row);
std::shuffle(to_process.begin(), to_process.end(), *random);
std::stable_sort(to_process.begin(), to_process.end(), [this](int a, int b) {
return rows_[a].cols.size() < rows_[b].cols.size();
});
for (const int row : to_process) {
ProcessSingletonColumns();
if (rows_[row].cols.empty()) continue;
if (rows_[row].slack > 1e-6) continue;
if (rows_[row].multipliers.size() > kMaxAggregationSize) continue;
// Heuristic: eliminate the variable with highest shifted lp value.
int eliminated_col = -1;
double max_lp_value = 0.0;
for (const int col : rows_[row].cols) {
if (shifted_lp_values_[col] > max_lp_value) {
max_lp_value = shifted_lp_values_[col];
eliminated_col = col;
}
}
if (eliminated_col == -1) continue;
EliminateVarUsingRow(eliminated_col, row);
}
// As an heuristic, we just try to add zero rows with an odd rhs and a low
// enough slack.
for (const auto& row : rows_) {
if (row.cols.empty() && row.rhs_parity && row.slack < kSlackThreshold) {
result.push_back(row.multipliers);
}
}
VLOG(1) << "#candidates: " << result.size() << " / " << rows_.size();
return result;
}
} // namespace sat
} // namespace operations_research

132
ortools/sat/zero_half_cuts.h Normal file
View File

@@ -0,0 +1,132 @@
// Copyright 2010-2018 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_SAT_ZERO_HALF_CUTS_H_
#define OR_TOOLS_SAT_ZERO_HALF_CUTS_H_
#include <vector>
#include "ortools/lp_data/lp_types.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/util.h"
namespace operations_research {
namespace sat {
// Heuristic to find a good sums of rows from the LP (with coeff -1, +1) that
// can lead to a violated zero-half cut (i.e. after integer rounding with a
// divisor 2).
//
// For this, all that matter is the parity of the coefficients and the rhs in
// the linear combination of the original problem constraint. So this class
// maintain a copy of the LP matrix modulo 2 on which simplification and
// heuristic are performed to find good cut candidates(s).
//
// Most of what is done here is described in the paper "Algorithms to Separate
// {0, 1/2}-Chvátal-Gomory Cuts", Arie M. C. A. Koster, Adrian Zymolka, Manuel
// Kutschka.
class ZeroHalfCutHelper {
public:
// Public API: ProcessVariables() must be called first and then constraints
// can be added one by one. Finally GetZeroHalfInterestingCuts() will return a
// set of good candidates.
//
// TODO(user): This is a first implementation, both the heuristic and the
// code performance can probably be improved uppon.
void ProcessVariables(const std::vector<double>& lp_values,
const std::vector<IntegerValue>& lower_bounds,
const std::vector<IntegerValue>& upper_bounds);
void AddOneConstraint(
glop::RowIndex,
const std::vector<std::pair<glop::ColIndex, IntegerValue>>& terms,
IntegerValue lb, IntegerValue ub);
std::vector<std::vector<std::pair<glop::RowIndex, IntegerValue>>>
InterestingCandidates(ModelRandomGenerator* random);
// Visible for testing.
void Reset(int size);
// Visible for testing.
//
// Boolean matrix. Each column correspond to one variable (col indices).
// Each row to a sum of the initial problem constraints. We store the
// coefficient modulo 2, so only the positions of the ones.
struct CombinationOfRows {
// How this row was formed from the initial problem constraints.
std::vector<std::pair<glop::RowIndex, IntegerValue>> multipliers;
// The index of the odd coefficient of this combination.
std::vector<int> cols;
// The parity of the rhs (1 for odd).
int rhs_parity;
// How tight this constraints is under the current LP solution.
double slack;
};
void AddBinaryRow(const CombinationOfRows& binary_row);
const CombinationOfRows& MatrixRow(int row) const { return rows_[row]; }
const std::vector<int>& MatrixCol(int col) const { return col_to_rows_[col]; }
// Visible for testing.
//
// Adds the given row to all other rows having an odd cofficient on the given
// column. This then eliminate the entry (col, row) that is now a singleton by
// incresing the slack of the given row.
void EliminateVarUsingRow(int col, int row);
// Visible for testing.
//
// Like std::set_symmetric_difference, but use a vector<bool> instead of sort.
// This assumes tmp_marked_ to be all false. We don't DCHECK it here for
// speed, but it DCHECKed on each EliminateVarUsingRow() call. In addition,
// the result is filtered using the extra_condition function.
void SymmetricDifference(std::function<bool(int)> extra_condition,
const std::vector<int>& a, std::vector<int>* b);
private:
void ProcessSingletonColumns();
// As we combine rows, when the activity of a combination get too far away
// from its bound, we just discard it. Note that the row will still be there
// but its index will not appear in the col-wise representation of the matrix.
const double kSlackThreshold = 0.5;
const int kMaxAggregationSize = 100;
// We don't consider long constraint or constraint with high magnitude, since
// the highest violation we can hope for is 1, and if the magnitude is large
// then the cut efficacity will not be great.
const int kMaxInputConstraintSize = 100;
const double kMaxInputConstraintMagnitude = 1e6;
// Variable information.
std::vector<double> lp_values_;
std::vector<double> shifted_lp_values_;
std::vector<int> bound_parity_;
// Binary matrix.
//
// Note that as we combine rows, we never move their indices. So after initial
// creation rows_ will always have the same size.
std::vector<CombinationOfRows> rows_;
std::vector<std::vector<int>> col_to_rows_;
std::vector<int> singleton_cols_;
// Temporary vector used by SymmetricDifference().
std::vector<bool> tmp_marked_;
};
} // namespace sat
} // namespace operations_research
#endif // OR_TOOLS_SAT_ZERO_HALF_CUTS_H_