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

simplify mccormick cut in CP-SAT; rearchitecture glop internals

Browse Source
This commit is contained in:
Laurent Perron
2019-06-21 13:39:09 +02:00
parent e5b9bad308
commit a66a929a4e
13 changed files with 243 additions and 219 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
ortools/glop/basis_representation.cc
View File

@@ -310,7 +310,7 @@ void BasisFactorization::LeftSolve(ScatteredRow* y) const {
if (use_middle_product_form_update_) {
lu_factorization_.LeftSolveUWithNonZeros(y);
rank_one_factorization_.LeftSolveWithNonZeros(y);
lu_factorization_.LeftSolveLWithNonZeros(y, nullptr);
lu_factorization_.LeftSolveLWithNonZeros(y);
y->SortNonZerosIfNeeded();
} else {
y->non_zeros.clear();
@@ -341,8 +341,8 @@ const DenseColumn& BasisFactorization::RightSolveForTau(
BumpDeterministicTimeForSolve(compact_matrix_.num_rows().value());
if (use_middle_product_form_update_) {
if (tau_computation_can_be_optimized_) {
// Once used, the intermediate result is overridden, so RightSolveForTau()
// can no longer use the optimized algorithm.
// Once used, the intermediate result is overwritten, so
// RightSolveForTau() can no longer use the optimized algorithm.
tau_computation_can_be_optimized_ = false;
lu_factorization_.RightSolveLWithPermutedInput(a.values, &tau_.values);
tau_.non_zeros.clear();
@@ -409,7 +409,7 @@ void BasisFactorization::LeftSolveForUnitRow(ColIndex j,
tau_.non_zeros.clear();
} else {
tau_computation_can_be_optimized_ = false;
lu_factorization_.LeftSolveLWithNonZeros(y, nullptr);
lu_factorization_.LeftSolveLWithNonZeros(y);
}
y->SortNonZerosIfNeeded();
}
@@ -423,7 +423,7 @@ void BasisFactorization::TemporaryLeftSolveForUnitRow(ColIndex j,
ClearAndResizeVectorWithNonZeros(RowToColIndex(compact_matrix_.num_rows()),
y);
lu_factorization_.LeftSolveUForUnitRow(j, y);
lu_factorization_.LeftSolveLWithNonZeros(y, nullptr);
lu_factorization_.LeftSolveLWithNonZeros(y);
y->SortNonZerosIfNeeded();
}

4
ortools/glop/lu_factorization.cc
View File

@@ -388,6 +388,10 @@ bool LuFactorization::LeftSolveLWithNonZeros(
}
}
void LuFactorization::LeftSolveLWithNonZeros(ScatteredRow* y) const {
LeftSolveLWithNonZeros(y, nullptr);
}
ColIndex LuFactorization::LeftSolveUForUnitRow(ColIndex col,
ScatteredRow* y) const {
SCOPED_TIME_STAT(&stats_);

1
ortools/glop/lu_factorization.h
View File

@@ -99,6 +99,7 @@ class LuFactorization {
// then false is returned.
bool LeftSolveLWithNonZeros(ScatteredRow* y,
ScatteredColumn* result_before_permutation) const;
void LeftSolveLWithNonZeros(ScatteredRow* y) const;
// Specialized version of RightSolveLWithNonZeros() that takes a SparseColumn
// or a ScatteredColumn as input. non_zeros will either be cleared or set to

4
ortools/glop/preprocessor.cc
View File

@@ -709,9 +709,9 @@ void ProportionalColumnPreprocessor::RecoverSolution(
const ColIndex representative = merged_columns_[col];
if (representative != kInvalidCol) {
if (IsFinite(distance_to_bound[representative])) {
// If the distance if finite, then each variable is set to its
// If the distance is finite, then each variable is set to its
// corresponding bound (the one from which the distance is computed) and
// is then changed by has much as possible until the distance is zero.
// is then changed by as much as possible until the distance is zero.
const Fractional bound_factor =
column_factors_[col] / column_factors_[representative];
const Fractional scaled_distance =

26
ortools/glop/reduced_costs.cc
View File

@@ -64,7 +64,7 @@ bool ReducedCosts::TestEnteringReducedCostPrecision(
}
const Fractional old_reduced_cost = reduced_costs_[entering_col];
const Fractional precise_reduced_cost =
objective_[entering_col] + objective_perturbation_[entering_col] -
objective_[entering_col] + cost_perturbations_[entering_col] -
PreciseScalarProduct(basic_objective_, direction);
// Update the reduced cost of the entering variable with the precise version.
@@ -128,7 +128,7 @@ Fractional ReducedCosts::ComputeMaximumDualResidual() const {
for (RowIndex row(0); row < num_rows; ++row) {
const ColIndex basic_col = basis_[row];
const Fractional residual =
objective_[basic_col] + objective_perturbation_[basic_col] -
objective_[basic_col] + cost_perturbations_[basic_col] -
matrix_.ColumnScalarProduct(basic_col, dual_values);
dual_residual_error = std::max(dual_residual_error, std::abs(residual));
}
@@ -246,7 +246,7 @@ void ReducedCosts::PerturbCosts() {
std::max(max_cost_magnitude, std::abs(objective_[col]));
}
objective_perturbation_.AssignToZero(matrix_.num_cols());
cost_perturbations_.AssignToZero(matrix_.num_cols());
for (ColIndex col(0); col < structural_size; ++col) {
const Fractional objective = objective_[col];
const Fractional magnitude =
@@ -264,10 +264,10 @@ void ReducedCosts::PerturbCosts() {
case VariableType::FIXED_VARIABLE:
break;
case VariableType::LOWER_BOUNDED:
objective_perturbation_[col] = magnitude;
cost_perturbations_[col] = magnitude;
break;
case VariableType::UPPER_BOUNDED:
objective_perturbation_[col] = -magnitude;
cost_perturbations_[col] = -magnitude;
break;
case VariableType::UPPER_AND_LOWER_BOUNDED:
// Here we don't necessarily maintain the dual-feasibility of a dual
@@ -276,9 +276,9 @@ void ReducedCosts::PerturbCosts() {
// done by MakeBoxedVariableDualFeasible() at the end of the dual
// phase-I algorithm.
if (objective > 0.0) {
objective_perturbation_[col] = magnitude;
cost_perturbations_[col] = magnitude;
} else if (objective < 0.0) {
objective_perturbation_[col] = -magnitude;
cost_perturbations_[col] = -magnitude;
}
break;
}
@@ -292,13 +292,13 @@ void ReducedCosts::ShiftCost(ColIndex col) {
dual_feasibility_tolerance_ *
(reduced_costs_[col] > 0.0 ? kToleranceFactor : -kToleranceFactor);
IF_STATS_ENABLED(stats_.cost_shift.Add(reduced_costs_[col] + small_step));
objective_perturbation_[col] -= reduced_costs_[col] + small_step;
cost_perturbations_[col] -= reduced_costs_[col] + small_step;
reduced_costs_[col] = -small_step;
}
void ReducedCosts::ClearAndRemoveCostShifts() {
SCOPED_TIME_STAT(&stats_);
objective_perturbation_.AssignToZero(matrix_.num_cols());
cost_perturbations_.AssignToZero(matrix_.num_cols());
recompute_basic_objective_ = true;
recompute_basic_objective_left_inverse_ = true;
recompute_reduced_costs_ = true;
@@ -339,12 +339,12 @@ void ReducedCosts::RecomputeReducedCostsAndPrimalEnteringCandidatesIfNeeded() {
void ReducedCosts::ComputeBasicObjective() {
SCOPED_TIME_STAT(&stats_);
const ColIndex num_cols_in_basis = RowToColIndex(matrix_.num_rows());
objective_perturbation_.resize(matrix_.num_cols(), 0.0);
cost_perturbations_.resize(matrix_.num_cols(), 0.0);
basic_objective_.resize(num_cols_in_basis, 0.0);
for (ColIndex col(0); col < num_cols_in_basis; ++col) {
const ColIndex basis_col = basis_[ColToRowIndex(col)];
basic_objective_[col] =
objective_[basis_col] + objective_perturbation_[basis_col];
objective_[basis_col] + cost_perturbations_[basis_col];
}
recompute_basic_objective_ = false;
recompute_basic_objective_left_inverse_ = true;
@@ -367,7 +367,7 @@ void ReducedCosts::ComputeReducedCosts() {
#endif
if (num_omp_threads == 1) {
for (ColIndex col(0); col < num_cols; ++col) {
reduced_costs_[col] = objective_[col] + objective_perturbation_[col] -
reduced_costs_[col] = objective_[col] + cost_perturbations_[col] -
matrix_.ColumnScalarProduct(
col, basic_objective_left_inverse_.values);
@@ -557,7 +557,7 @@ void ReducedCosts::UpdateBasicObjective(ColIndex entering_col,
RowIndex leaving_row) {
SCOPED_TIME_STAT(&stats_);
basic_objective_[RowToColIndex(leaving_row)] =
objective_[entering_col] + objective_perturbation_[entering_col];
objective_[entering_col] + cost_perturbations_[entering_col];
recompute_basic_objective_left_inverse_ = true;
}

9
ortools/glop/reduced_costs.h
View File

@@ -180,9 +180,7 @@ class ReducedCosts {
bool IsValidPrimalEnteringCandidate(ColIndex col) const;
// Visible for testing.
const DenseRow& GetCostPerturbations() const {
return objective_perturbation_;
}
const DenseRow& GetCostPerturbations() const { return cost_perturbations_; }
private:
// Statistics about this class.
@@ -256,10 +254,7 @@ class ReducedCosts {
// Perturbations to the objective function. This may be introduced to
// counter degenerecency. It will be removed at the end of the algorithm.
//
// TODO(user): rename this cost_perturbations_ to be more consistent with
// the literature.
DenseRow objective_perturbation_;
DenseRow cost_perturbations_;
// Reduced costs of the relevant columns of A.
DenseRow reduced_costs_;

4
ortools/lp_data/lp_data.h
View File

@@ -180,7 +180,7 @@ class LinearProgram {
// modifying the matrix does not change the result of any function in this
// class until UseTransposeMatrixAsReference() is called. This is because the
// transpose matrix is only used by GetTransposeSparseMatrix() and this
// function will recompute the whole tranpose from the matrix. In particular,
// function will recompute the whole transpose from the matrix. In particular,
// do not call GetTransposeSparseMatrix() while you modify the matrix returned
// by GetMutableTransposeSparseMatrix() otherwise all your changes will be
// lost.
@@ -570,7 +570,7 @@ class LinearProgram {
SparseMatrix matrix_;
// The transpose of matrix_. This will be lazily recomputed by
// GetTransposeSparseMatrix() if tranpose_matrix_is_consistent_ is false.
// GetTransposeSparseMatrix() if transpose_matrix_is_consistent_ is false.
mutable SparseMatrix transpose_matrix_;
// Constraint related quantities.

58
ortools/lp_data/sparse.cc
View File

@@ -859,16 +859,26 @@ void TriangularMatrix::TransposeLowerSolveInternal(DenseColumn* rhs) const {
}
}
// TODO(user): exploit all_diagonal_coefficients_are_one_ when true in
// all the hyper-sparse functions.
void TriangularMatrix::HyperSparseSolve(DenseColumn* rhs,
RowIndexVector* non_zero_rows) const {
if (all_diagonal_coefficients_are_one_) {
HyperSparseSolveInternal<true>(rhs, non_zero_rows);
} else {
HyperSparseSolveInternal<false>(rhs, non_zero_rows);
}
}
template <bool diagonal_of_ones>
void TriangularMatrix::HyperSparseSolveInternal(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const {
RETURN_IF_NULL(rhs);
int new_size = 0;
for (const RowIndex row : *non_zero_rows) {
if ((*rhs)[row] == 0.0) continue;
const ColIndex row_as_col = RowToColIndex(row);
const Fractional coeff = (*rhs)[row] / diagonal_coefficients_[row_as_col];
const Fractional coeff =
diagonal_of_ones ? (*rhs)[row]
: (*rhs)[row] / diagonal_coefficients_[row_as_col];
(*rhs)[row] = coeff;
for (const EntryIndex i : Column(row_as_col)) {
(*rhs)[EntryRow(i)] -= coeff * EntryCoefficient(i);
@@ -881,12 +891,24 @@ void TriangularMatrix::HyperSparseSolve(DenseColumn* rhs,
void TriangularMatrix::HyperSparseSolveWithReversedNonZeros(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const {
if (all_diagonal_coefficients_are_one_) {
HyperSparseSolveWithReversedNonZerosInternal<true>(rhs, non_zero_rows);
} else {
HyperSparseSolveWithReversedNonZerosInternal<false>(rhs, non_zero_rows);
}
}
template <bool diagonal_of_ones>
void TriangularMatrix::HyperSparseSolveWithReversedNonZerosInternal(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const {
RETURN_IF_NULL(rhs);
int new_start = non_zero_rows->size();
for (const RowIndex row : Reverse(*non_zero_rows)) {
if ((*rhs)[row] == 0.0) continue;
const ColIndex row_as_col = RowToColIndex(row);
const Fractional coeff = (*rhs)[row] / diagonal_coefficients_[row_as_col];
const Fractional coeff =
diagonal_of_ones ? (*rhs)[row]
: (*rhs)[row] / diagonal_coefficients_[row_as_col];
(*rhs)[row] = coeff;
for (const EntryIndex i : Column(row_as_col)) {
(*rhs)[EntryRow(i)] -= coeff * EntryCoefficient(i);
@@ -900,6 +922,16 @@ void TriangularMatrix::HyperSparseSolveWithReversedNonZeros(
void TriangularMatrix::TransposeHyperSparseSolve(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const {
if (all_diagonal_coefficients_are_one_) {
TransposeHyperSparseSolveInternal<true>(rhs, non_zero_rows);
} else {
TransposeHyperSparseSolveInternal<false>(rhs, non_zero_rows);
}
}
template <bool diagonal_of_ones>
void TriangularMatrix::TransposeHyperSparseSolveInternal(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const {
RETURN_IF_NULL(rhs);
int new_size = 0;
for (const RowIndex row : *non_zero_rows) {
@@ -908,7 +940,8 @@ void TriangularMatrix::TransposeHyperSparseSolve(
for (const EntryIndex i : Column(row_as_col)) {
sum -= EntryCoefficient(i) * (*rhs)[EntryRow(i)];
}
(*rhs)[row] = sum / diagonal_coefficients_[row_as_col];
(*rhs)[row] =
diagonal_of_ones ? sum : sum / diagonal_coefficients_[row_as_col];
if (sum != 0.0) {
(*non_zero_rows)[new_size] = row;
++new_size;
@@ -919,6 +952,18 @@ void TriangularMatrix::TransposeHyperSparseSolve(
void TriangularMatrix::TransposeHyperSparseSolveWithReversedNonZeros(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const {
if (all_diagonal_coefficients_are_one_) {
TransposeHyperSparseSolveWithReversedNonZerosInternal<true>(rhs,
non_zero_rows);
} else {
TransposeHyperSparseSolveWithReversedNonZerosInternal<false>(rhs,
non_zero_rows);
}
}
template <bool diagonal_of_ones>
void TriangularMatrix::TransposeHyperSparseSolveWithReversedNonZerosInternal(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const {
RETURN_IF_NULL(rhs);
int new_start = non_zero_rows->size();
for (const RowIndex row : Reverse(*non_zero_rows)) {
@@ -932,7 +977,8 @@ void TriangularMatrix::TransposeHyperSparseSolveWithReversedNonZeros(
for (; i >= i_end; --i) {
sum -= EntryCoefficient(i) * (*rhs)[EntryRow(i)];
}
(*rhs)[row] = sum / diagonal_coefficients_[row_as_col];
(*rhs)[row] =
diagonal_of_ones ? sum : sum / diagonal_coefficients_[row_as_col];
if (sum != 0.0) {
--new_start;
(*non_zero_rows)[new_start] = row;

12
ortools/lp_data/sparse.h
View File

@@ -734,6 +734,18 @@ class TriangularMatrix : private CompactSparseMatrix {
void TransposeLowerSolveInternal(DenseColumn* rhs) const;
template <bool diagonal_of_ones>
void TransposeUpperSolveInternal(DenseColumn* rhs) const;
template <bool diagonal_of_ones>
void HyperSparseSolveInternal(DenseColumn* rhs,
RowIndexVector* non_zero_rows) const;
template <bool diagonal_of_ones>
void HyperSparseSolveWithReversedNonZerosInternal(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const;
template <bool diagonal_of_ones>
void TransposeHyperSparseSolveInternal(DenseColumn* rhs,
RowIndexVector* non_zero_rows) const;
template <bool diagonal_of_ones>
void TransposeHyperSparseSolveWithReversedNonZerosInternal(
DenseColumn* rhs, RowIndexVector* non_zero_rows) const;
// Internal function used by the Add*() functions to finish adding
// a new column to a triangular matrix.

52
ortools/sat/cp_model_solver.cc
View File

@@ -379,11 +379,6 @@ std::string Summarize(const std::string& input) {
input.substr(input.size() - half, half));
}
bool IsPositive(IntegerVariable v, Model* m) {
IntegerTrail* const integer_trail = m->GetOrCreate<IntegerTrail>();
return integer_trail->LevelZeroLowerBound(v) >= 0;
}
} // namespace.
// =============================================================================
@@ -733,17 +728,46 @@ void TryToAddCutGenerators(const CpModelProto& model_proto,
if (HasEnforcementLiteral(ct)) return;
if (ct.int_prod().vars_size() != 2) return;
const int target = ct.int_prod().target();
const IntegerVariable v0 = mapping->Integer(ct.int_prod().vars(0));
const IntegerVariable v1 = mapping->Integer(ct.int_prod().vars(1));
// Constraint is z == x * y.
IntegerVariable z = mapping->Integer(ct.int_prod().target());
IntegerVariable x = mapping->Integer(ct.int_prod().vars(0));
IntegerVariable y = mapping->Integer(ct.int_prod().vars(1));
IntegerTrail* const integer_trail = m->GetOrCreate<IntegerTrail>();
IntegerValue x_lb = integer_trail->LowerBound(x);
IntegerValue x_ub = integer_trail->UpperBound(x);
IntegerValue y_lb = integer_trail->LowerBound(y);
IntegerValue y_ub = integer_trail->UpperBound(y);
if (x == y) {
// We currently only support variables with non-negative domains.
if (x_lb < 0 && x_ub > 0) return;
// Change the sigh of x if its domain is non-positive.
if (x_ub <= 0) {
x = NegationOf(x);
}
relaxation->cut_generators.push_back(CreateSquareCutGenerator(z, x, m));
} else {
// We currently only support variables with non-negative domains.
if (x_lb < 0 && x_ub > 0) return;
if (y_lb < 0 && y_ub > 0) return;
// Change signs to return to the case where all variables are a domain
// with non negative values only.
if (x_ub <= 0) {
x = NegationOf(x);
z = NegationOf(z);
}
if (y_ub <= 0) {
y = NegationOf(y);
z = NegationOf(z);
}
if (v0 == v1 && IsPositive(v0, m)) {
relaxation->cut_generators.push_back(
CreateSquareCutGenerator(mapping->Integer(target), v0, m));
} else if (IsPositive(v0, m) && IsPositive(v1, m)) {
relaxation->cut_generators.push_back(
CreatePositiveMultiplicationCutGenerator(mapping->Integer(target), v0,
v1, m));
CreatePositiveMultiplicationCutGenerator(z, x, y, m));
}
}
}

191
ortools/sat/cuts.cc
View File

@@ -784,140 +784,157 @@ void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
DivideByGCD(cut);
}
CutGenerator CreatePositiveMultiplicationCutGenerator(
IntegerVariable target_var, IntegerVariable v1, IntegerVariable v2,
Model* model) {
CutGenerator CreatePositiveMultiplicationCutGenerator(IntegerVariable z,
IntegerVariable x,
IntegerVariable y,
Model* model) {
CutGenerator result;
result.vars = {target_var, v1, v2};
result.vars = {z, x, y};
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
result.generate_cuts =
[target_var, v1, v2, model, integer_trail](
[z, x, y, integer_trail](
const gtl::ITIVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
const int64 v1_ub = integer_trail->LevelZeroUpperBound(v1).value();
const int64 v1_lb = integer_trail->LevelZeroLowerBound(v1).value();
const int64 v2_ub = integer_trail->LevelZeroUpperBound(v2).value();
const int64 v2_lb = integer_trail->LevelZeroLowerBound(v2).value();
const int64 x_lb = integer_trail->LevelZeroLowerBound(x).value();
const int64 x_ub = integer_trail->LevelZeroUpperBound(x).value();
const int64 y_lb = integer_trail->LevelZeroLowerBound(y).value();
const int64 y_ub = integer_trail->LevelZeroUpperBound(y).value();
const int64 kMaxSafeInteger = int64{9007199254740991}; // 2 ^ 53 - 1.
// TODO(user): Compute a better bound (int_max / 4 ?).
const int64 kMaxSafeInteger = (int64{1} << 53) - 1;
if (CapProd(v1_ub, v2_ub) >= kMaxSafeInteger) {
if (CapProd(x_ub, y_ub) >= kMaxSafeInteger) {
VLOG(3) << "Potential overflow in PositiveMultiplicationCutGenerator";
return;
}
const double target_value = lp_values[target_var];
const double v1_value = lp_values[v1];
const double v2_value = lp_values[v2];
const double x_value = lp_values[x];
const double y_value = lp_values[y];
const double z_value = lp_values[z];
if (target_value <= v1_value * v2_lb + v2_value * v1_lb -
v1_lb * v2_lb - kMinCutViolation) {
// cut: target >= v1 * v2_lb + v2 * v2_lb - v1_lb * v2_lb
LinearConstraint cut;
cut.vars.push_back(target_var);
cut.coeffs.push_back(IntegerValue(1));
if (v2_lb != 0) {
cut.vars.push_back(v1);
cut.coeffs.push_back(IntegerValue(-v2_lb));
// TODO: As the bounds change monotonically, these cuts dominate any
// previous one. try to keep a reference to the cut and replace it.
// Alternatively, add an API for a level-zero bound change callback.
// We implement the McCormick relaxation of bilinear constraints.
// Cut -z + x_coeff * x + y_coeff* y <= rhs
auto try_add_above_cut = [manager, z_value, x_value, y_value, x, y, z,
lp_values](int64 x_coeff, int64 y_coeff,
int64 rhs) {
if (-z_value + x_value * x_coeff + y_value * y_coeff >=
rhs + kMinCutViolation) {
// cut: -z + x * x_coeff + y * y_coeff <= rhs
LinearConstraint cut;
cut.vars.push_back(z);
cut.coeffs.push_back(IntegerValue(-1));
if (x_coeff != 0) {
cut.vars.push_back(x);
cut.coeffs.push_back(IntegerValue(x_coeff));
}
if (y_coeff != 0) {
cut.vars.push_back(y);
cut.coeffs.push_back(IntegerValue(y_coeff));
}
cut.lb = kMinIntegerValue;
cut.ub = IntegerValue(rhs);
manager->AddCut(cut, "PositiveProduct", lp_values);
}
if (v1_lb != 0) {
cut.vars.push_back(v2);
cut.coeffs.push_back(IntegerValue(-v1_lb));
};
// Cut -z + x_coeff * x + y_coeff* y >= rhs
auto try_add_below_cut = [manager, z_value, x_value, y_value, x, y, z,
lp_values](int64 x_coeff, int64 y_coeff,
int64 rhs) {
if (-z_value + x_value * x_coeff + y_value * y_coeff <=
rhs - kMinCutViolation) {
// cut: -z + x * x_coeff + y * y_coeff >= rhs
LinearConstraint cut;
cut.vars.push_back(z);
cut.coeffs.push_back(IntegerValue(-1));
if (x_coeff != 0) {
cut.vars.push_back(x);
cut.coeffs.push_back(IntegerValue(x_coeff));
}
if (y_coeff != 0) {
cut.vars.push_back(y);
cut.coeffs.push_back(IntegerValue(y_coeff));
}
cut.lb = IntegerValue(rhs);
cut.ub = kMaxIntegerValue;
manager->AddCut(cut, "PositiveProduct", lp_values);
}
cut.lb = IntegerValue(-v1_lb * v2_lb);
cut.ub = kMaxIntegerValue;
manager->AddCut(cut, "PositiveProduct", lp_values);
}
};
if (target_value >= v2_value * v1_ub + kMinCutViolation) {
// cut: target <= v2 * v1_ub.
LinearConstraint cut;
cut.vars.push_back(target_var);
cut.coeffs.push_back(IntegerValue(1));
cut.vars.push_back(v2);
cut.coeffs.push_back(IntegerValue(-v1_ub));
cut.lb = kMinIntegerValue;
cut.ub = IntegerValue(0);
manager->AddCut(cut, "PositiveProduct", lp_values);
}
// McCormick cuts.
try_add_above_cut(y_lb, x_lb, x_lb * y_lb);
try_add_above_cut(y_ub, x_ub, x_ub * y_ub);
if (target_value >= v1_value * v2_ub + kMinCutViolation) {
// cut: target <= v1 * v2_ub.
LinearConstraint cut;
cut.vars.push_back(target_var);
cut.coeffs.push_back(IntegerValue(1));
cut.vars.push_back(v1);
cut.coeffs.push_back(IntegerValue(-v2_ub));
cut.lb = kMinIntegerValue;
cut.ub = IntegerValue(0);
manager->AddCut(cut, "PositiveProduct", lp_values);
}
try_add_below_cut(y_ub, x_lb, x_lb * y_ub);
try_add_below_cut(y_lb, x_ub, x_ub * y_lb);
};
return result;
}
CutGenerator CreateSquareCutGenerator(IntegerVariable target_var,
IntegerVariable int_var, Model* model) {
CutGenerator CreateSquareCutGenerator(IntegerVariable y, IntegerVariable x,
Model* model) {
CutGenerator result;
result.vars = {target_var, int_var};
result.vars = {y, x};
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
result.generate_cuts =
[target_var, int_var, model, integer_trail](
[y, x, model, integer_trail](
const gtl::ITIVector<IntegerVariable, double>& lp_values,
LinearConstraintManager* manager) {
const int64 var_ub =
integer_trail->LevelZeroUpperBound(int_var).value();
const int64 var_lb =
integer_trail->LevelZeroLowerBound(int_var).value();
const int64 x_ub = integer_trail->LevelZeroUpperBound(x).value();
const int64 x_lb = integer_trail->LevelZeroLowerBound(x).value();
if (var_lb == var_ub) return;
if (x_lb == x_ub) return;
// Check for potential overflows.
if (var_ub > int64{1000000000}) return;
DCHECK_GE(var_lb, 0);
if (x_ub > int64{1000000000}) return;
DCHECK_GE(x_lb, 0);
const double target_value = lp_values[target_var];
const double var_value = lp_values[int_var];
const double y_value = lp_values[y];
const double x_value = lp_values[x];
// First cut: target should be below the line:
// (var_lb, val_lb ^ 2) to (var_ub, var_ub ^ 2).
// (x_lb, val_lb ^ 2) to (x_ub, x_ub ^ 2).
// The slope of that line is (ub^2 - lb^2) / (ub - lb) = ub + lb.
const int64 target_lb = var_lb * var_lb;
const int64 above_slope = var_ub + var_lb;
const double max_target =
target_lb + above_slope * (var_value - var_lb);
if (target_value >= max_target + kMinCutViolation) {
// cut: target <= (var_lb + var_ub) * int_var - var_lb * var_ub
const int64 y_lb = x_lb * x_lb;
const int64 above_slope = x_ub + x_lb;
const double max_y = y_lb + above_slope * (x_value - x_lb);
if (y_value >= max_y + kMinCutViolation) {
// cut: target <= (x_lb + x_ub) * x - x_lb * x_ub
LinearConstraint above_cut;
above_cut.vars.push_back(target_var);
above_cut.vars.push_back(y);
above_cut.coeffs.push_back(IntegerValue(1));
above_cut.vars.push_back(int_var);
above_cut.vars.push_back(x);
above_cut.coeffs.push_back(IntegerValue(-above_slope));
above_cut.lb = kMinIntegerValue;
above_cut.ub = IntegerValue(-var_lb * var_ub);
above_cut.ub = IntegerValue(-x_lb * x_ub);
manager->AddCut(above_cut, "SquareUpper", lp_values);
}
// Second cut: target should be above all the lines
// (value, value ^ 2) to (value + 1, (value + 1) ^ 2)
// The slope of that line is 2 * value + 1
const int64 var_floor = static_cast<int64>(std::floor(var_value));
const int64 below_slope = 2 * var_floor + 1;
const double min_target =
below_slope * var_value - var_floor - var_floor * var_floor;
if (min_target >= target_value + kMinCutViolation) {
// cut: target >= below_slope * (int_var - var_floor) +
// var_floor * var_floor
// : target >= below_slope * int_var - var_floor ^ 2 - var_floor
const int64 x_floor = static_cast<int64>(std::floor(x_value));
const int64 below_slope = 2 * x_floor + 1;
const double min_y =
below_slope * x_value - x_floor - x_floor * x_floor;
if (min_y >= y_value + kMinCutViolation) {
// cut: target >= below_slope * (x - x_floor) +
// x_floor * x_floor
// : target >= below_slope * x - x_floor ^ 2 - x_floor
LinearConstraint below_cut;
below_cut.vars.push_back(target_var);
below_cut.vars.push_back(y);
below_cut.coeffs.push_back(IntegerValue(1));
below_cut.vars.push_back(int_var);
below_cut.vars.push_back(x);
below_cut.coeffs.push_back(-IntegerValue(below_slope));
below_cut.lb = IntegerValue(-var_floor - var_floor * var_floor);
below_cut.lb = IntegerValue(-x_floor - x_floor * x_floor);
below_cut.ub = kMaxIntegerValue;
manager->AddCut(below_cut, "SquareLower", lp_values);
}

11
ortools/sat/cuts.h
View File

@@ -256,14 +256,15 @@ CutGenerator CreateKnapsackCoverCutGenerator(
const std::vector<IntegerVariable>& vars, Model* model);
// A cut generator for z = x * y (x and y >= 0)
CutGenerator CreatePositiveMultiplicationCutGenerator(
IntegerVariable target_var, IntegerVariable v1, IntegerVariable v2,
Model* model);
CutGenerator CreatePositiveMultiplicationCutGenerator(IntegerVariable z,
IntegerVariable x,
IntegerVariable y,
Model* model);
// A cut generator for y = x ^ 2. (x >= 0).
// It will dynamically add a linear inequality to push y closer to the parabola.
CutGenerator CreateSquareCutGenerator(IntegerVariable target_var,
IntegerVariable int_var, Model* model);
CutGenerator CreateSquareCutGenerator(IntegerVariable y, IntegerVariable x,
Model* model);
} // namespace sat
} // namespace operations_research

80
ortools/sat/linear_relaxation.cc
View File

@@ -249,6 +249,8 @@ void AppendEnforcedUpperBound(const Literal enforcing_lit,
// constraints. The highest linearization_level is, the more types of constraint
// we encode. This method should be called only for linearization_level > 0.
//
// Note: IntProd is linearized dynamically using the cut generators.
//
// TODO(user): In full generality, we could encode all the constraint as an LP.
// TODO(user,user): Add unit tests for this method.
void TryToLinearizeConstraint(const CpModelProto& model_proto,
@@ -320,84 +322,6 @@ void TryToLinearizeConstraint(const CpModelProto& model_proto,
NegationOf(mapping->Integers(ct.int_min().vars()));
AppendMaxRelaxation(negative_target, negative_vars, linearization_level,
model, relaxation);
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kIntProd) {
if (HasEnforcementLiteral(ct)) return;
if (ct.int_prod().vars_size() != 2) return;
IntegerVariable target = mapping->Integer(ct.int_prod().target());
IntegerVariable v0 = mapping->Integer(ct.int_prod().vars(0));
IntegerVariable v1 = mapping->Integer(ct.int_prod().vars(1));
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
// We just linearize x = y^2 by x >= y which is far from ideal but at
// least pushes x when y moves away from zero. Note that if y is negative,
// we should probably also add x >= -y, but then this do not happen in
// our test set.
if (v0 == v1) { // target = v0 ^ 2
if (integer_trail->UpperBound(v0) <= 0) {
v0 = NegationOf(v0);
}
const IntegerValue var_lb = integer_trail->LowerBound(v0);
const IntegerValue var_ub = integer_trail->UpperBound(v0);
CHECK_GE(var_lb, 0);
// target >= (2 * var_lb + 1) * v0 - var_lb - var_lb ^ 2.
LinearConstraintBuilder lc_below(model, -var_lb - var_lb * var_lb,
kMaxIntegerValue);
lc_below.AddTerm(target, IntegerValue(1));
lc_below.AddTerm(v0, 2 * var_lb + 1);
relaxation->linear_constraints.push_back(lc_below.Build());
// target <= (var_lb + var_ub) * v0 - var_lb * var_ub
LinearConstraintBuilder lc_above(model, kMinIntegerValue,
-var_lb * var_ub);
lc_above.AddTerm(target, IntegerValue(1));
lc_above.AddTerm(v0, -(var_lb + var_ub));
relaxation->linear_constraints.push_back(lc_above.Build());
} else { // target = v0 * v1
// Change signs if needed.
if (integer_trail->UpperBound(v0) <= 0) {
v0 = NegationOf(v0);
target = NegationOf(target);
}
if (integer_trail->UpperBound(v1) <= 0) {
v1 = NegationOf(v1);
target = NegationOf(target);
}
const IntegerValue v0_lb = integer_trail->LowerBound(v0);
const IntegerValue v0_ub = integer_trail->UpperBound(v0);
const IntegerValue v1_lb = integer_trail->LowerBound(v1);
const IntegerValue v1_ub = integer_trail->UpperBound(v1);
CHECK_GE(v0_lb, 0);
CHECK_GE(v1_lb, 0);
// target >= v1 * v0_lb + v0 * v1_lb - v0_lb * v1_lb
LinearConstraintBuilder lc_below(model, -v0_lb * v1_lb, kMaxIntegerValue);
lc_below.AddTerm(target, IntegerValue(1));
if (v1_lb > 0) {
lc_below.AddTerm(v0, v1_lb);
}
if (v0_lb > 0) {
lc_below.AddTerm(v1, v0_lb);
}
relaxation->linear_constraints.push_back(lc_below.Build());
// target <= v0 * v1_ub
LinearConstraintBuilder lc_v0(model, kMinIntegerValue, IntegerValue(0));
lc_v0.AddTerm(target, IntegerValue(1));
lc_v0.AddTerm(v0, -v1_ub);
relaxation->linear_constraints.push_back(lc_v0.Build());
// target <= v1 * v0_ub
LinearConstraintBuilder lc_v1(model, kMinIntegerValue, IntegerValue(0));
lc_v1.AddTerm(target, IntegerValue(1));
lc_v1.AddTerm(v1, -v0_ub);
relaxation->linear_constraints.push_back(lc_v1.Build());
}
} else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinear) {
AppendLinearConstraintRelaxation(ct, linearization_level, *model,
relaxation);