diff --git a/ortools/sat/BUILD b/ortools/sat/BUILD index c355f6d354..be4bd607af 100644 --- a/ortools/sat/BUILD +++ b/ortools/sat/BUILD @@ -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"], diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 9aece62d3c..62bc49b866 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -2927,8 +2927,8 @@ bool CpModelPresolver::PresolveCircuit(ConstraintProto* ct) { // here (even if they are at false). std::vector> incoming_arcs; std::vector> 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(v, c))) { continue; } var_constraint_pair_already_called.insert({v, c}); + if (!in_queue[c]) { in_queue[c] = true; queue.push_back(c); diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 10922947f8..577586c11b 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -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>& 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. diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 0b4078faa8..afa2df5a40 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -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>& 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 tmp_dense_vector_; diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 26865cf5b1..3bd824be6b 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -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. diff --git a/ortools/sat/zero_half_cuts.cc b/ortools/sat/zero_half_cuts.cc new file mode 100644 index 0000000000..b1998556e2 --- /dev/null +++ b/ortools/sat/zero_half_cuts.cc @@ -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& lp_values, + const std::vector& lower_bounds, + const std::vector& 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>& 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 extra_condition, const std::vector& a, + std::vector* 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>> +ZeroHalfCutHelper::InterestingCandidates(ModelRandomGenerator* random) { + std::vector>> 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 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 diff --git a/ortools/sat/zero_half_cuts.h b/ortools/sat/zero_half_cuts.h new file mode 100644 index 0000000000..9c5ce1db13 --- /dev/null +++ b/ortools/sat/zero_half_cuts.h @@ -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 + +#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& lp_values, + const std::vector& lower_bounds, + const std::vector& upper_bounds); + void AddOneConstraint( + glop::RowIndex, + const std::vector>& terms, + IntegerValue lb, IntegerValue ub); + std::vector>> + 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> multipliers; + + // The index of the odd coefficient of this combination. + std::vector 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& 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 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 extra_condition, + const std::vector& a, std::vector* 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 lp_values_; + std::vector shifted_lp_values_; + std::vector 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 rows_; + std::vector> col_to_rows_; + std::vector singleton_cols_; + + // Temporary vector used by SymmetricDifference(). + std::vector tmp_marked_; +}; + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_ZERO_HALF_CUTS_H_