diff --git a/ortools/glop/basis_representation.cc b/ortools/glop/basis_representation.cc index 0fbcd4d622..0f31b7991d 100644 --- a/ortools/glop/basis_representation.cc +++ b/ortools/glop/basis_representation.cc @@ -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(); } diff --git a/ortools/glop/lu_factorization.cc b/ortools/glop/lu_factorization.cc index fb0573fae6..a44e90f807 100644 --- a/ortools/glop/lu_factorization.cc +++ b/ortools/glop/lu_factorization.cc @@ -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_); diff --git a/ortools/glop/lu_factorization.h b/ortools/glop/lu_factorization.h index c16f026a59..a6176a49a3 100644 --- a/ortools/glop/lu_factorization.h +++ b/ortools/glop/lu_factorization.h @@ -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 diff --git a/ortools/glop/preprocessor.cc b/ortools/glop/preprocessor.cc index ce054b593e..240840bcac 100644 --- a/ortools/glop/preprocessor.cc +++ b/ortools/glop/preprocessor.cc @@ -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 = diff --git a/ortools/glop/reduced_costs.cc b/ortools/glop/reduced_costs.cc index 4fa2fef20d..5ad7b2a31f 100644 --- a/ortools/glop/reduced_costs.cc +++ b/ortools/glop/reduced_costs.cc @@ -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; } diff --git a/ortools/glop/reduced_costs.h b/ortools/glop/reduced_costs.h index cb5b0e2832..a392b0b32e 100644 --- a/ortools/glop/reduced_costs.h +++ b/ortools/glop/reduced_costs.h @@ -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_; diff --git a/ortools/lp_data/lp_data.h b/ortools/lp_data/lp_data.h index 8204732346..7200443ef1 100644 --- a/ortools/lp_data/lp_data.h +++ b/ortools/lp_data/lp_data.h @@ -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. diff --git a/ortools/lp_data/sparse.cc b/ortools/lp_data/sparse.cc index 372d4bec99..dfcdc9197f 100644 --- a/ortools/lp_data/sparse.cc +++ b/ortools/lp_data/sparse.cc @@ -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(rhs, non_zero_rows); + } else { + HyperSparseSolveInternal(rhs, non_zero_rows); + } +} + +template +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(rhs, non_zero_rows); + } else { + HyperSparseSolveWithReversedNonZerosInternal(rhs, non_zero_rows); + } +} + +template +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(rhs, non_zero_rows); + } else { + TransposeHyperSparseSolveInternal(rhs, non_zero_rows); + } +} + +template +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(rhs, + non_zero_rows); + } else { + TransposeHyperSparseSolveWithReversedNonZerosInternal(rhs, + non_zero_rows); + } +} + +template +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; diff --git a/ortools/lp_data/sparse.h b/ortools/lp_data/sparse.h index 0a4807c1b5..c63d2b4a73 100644 --- a/ortools/lp_data/sparse.h +++ b/ortools/lp_data/sparse.h @@ -734,6 +734,18 @@ class TriangularMatrix : private CompactSparseMatrix { void TransposeLowerSolveInternal(DenseColumn* rhs) const; template void TransposeUpperSolveInternal(DenseColumn* rhs) const; + template + void HyperSparseSolveInternal(DenseColumn* rhs, + RowIndexVector* non_zero_rows) const; + template + void HyperSparseSolveWithReversedNonZerosInternal( + DenseColumn* rhs, RowIndexVector* non_zero_rows) const; + template + void TransposeHyperSparseSolveInternal(DenseColumn* rhs, + RowIndexVector* non_zero_rows) const; + template + 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. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 8a19996ced..7ca1e50193 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -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(); - 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(); + 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)); } } } diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index c8d3217c94..b1f6933eb3 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -784,140 +784,157 @@ void IntegerRoundingCut(RoundingOptions options, std::vector 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* const integer_trail = model->GetOrCreate(); result.generate_cuts = - [target_var, v1, v2, model, integer_trail]( + [z, x, y, integer_trail]( const gtl::ITIVector& 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(); result.generate_cuts = - [target_var, int_var, model, integer_trail]( + [y, x, model, integer_trail]( const gtl::ITIVector& 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(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(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); } diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index a7593f496f..361a17a5c9 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -256,14 +256,15 @@ CutGenerator CreateKnapsackCoverCutGenerator( const std::vector& 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 diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc index 868ca81c03..0594d1eaa6 100644 --- a/ortools/sat/linear_relaxation.cc +++ b/ortools/sat/linear_relaxation.cc @@ -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(); - - // 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);