From e8e53ae6eea804033c6322f0dbd31844ec87b602 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 18 Oct 2022 14:19:15 +0200 Subject: [PATCH] [CP-SAT] coefficient strenghtening --- ortools/sat/cp_model_presolve.cc | 57 ++++++++--- ortools/sat/linear_constraint_manager.cc | 122 +++++++++++++++-------- 2 files changed, 119 insertions(+), 60 deletions(-) diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index d0d8b45de2..d6cc40fbcb 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -2593,6 +2593,7 @@ void CpModelPresolver::TryToReduceCoefficientsOfLinearConstraint( } if (max_error == 0) break; // Last term. + if ((!use_ub || max_error <= PositiveRemainder(rhs_ub, IntegerValue(gcd))) && (!use_lb || @@ -2601,7 +2602,6 @@ void CpModelPresolver::TryToReduceCoefficientsOfLinearConstraint( // can be ignored. int64_t shift_lb = 0; int64_t shift_ub = 0; - int64_t local_max_variation = 0; context_->UpdateRuleStats("linear: simplify using partial gcd"); LinearConstraintProto* mutable_linear = ct->mutable_linear(); @@ -2612,7 +2612,6 @@ void CpModelPresolver::TryToReduceCoefficientsOfLinearConstraint( if (m < e.magnitude) continue; shift_lb += lbs[i] * m; shift_ub += ubs[i] * m; - local_max_variation += m * (ubs[i] - lbs[i]); mutable_linear->add_vars(vars[i]); mutable_linear->add_coeffs(coeffs[i]); } @@ -3215,12 +3214,14 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint( int64_t min_sum = 0; int64_t max_sum = 0; int64_t max_coeff_magnitude = 0; + int64_t min_coeff_magnitude = std::numeric_limits::max(); for (int i = 0; i < num_vars; ++i) { const int ref = arg.vars(i); const int64_t coeff = arg.coeffs(i); const int64_t term_a = coeff * context_->MinOf(ref); const int64_t term_b = coeff * context_->MaxOf(ref); max_coeff_magnitude = std::max(max_coeff_magnitude, std::abs(coeff)); + min_coeff_magnitude = std::min(min_coeff_magnitude, std::abs(coeff)); min_sum += std::min(term_a, term_b); max_sum += std::max(term_a, term_b); } @@ -3235,7 +3236,10 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint( const int64_t ub_threshold = domain[domain.size() - 2] - min_sum; const int64_t lb_threshold = max_sum - domain[1]; const Domain rhs_domain = ReadDomainFromProto(ct->linear()); - if (max_coeff_magnitude < std::max(ub_threshold, lb_threshold)) return; + if (max_coeff_magnitude + min_coeff_magnitude < + std::max(ub_threshold, lb_threshold)) { + return; + } // We need the constraint to be only bounded on one side in order to extract // enforcement literal. @@ -3258,6 +3262,9 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint( const bool upper_bounded = max_sum > rhs_domain.Max(); if (!lower_bounded && !upper_bounded) return; if (lower_bounded && upper_bounded) { + // Lets not split except if we extract enforcement. + if (max_coeff_magnitude < std::max(ub_threshold, lb_threshold)) return; + context_->UpdateRuleStats("linear: split boxed constraint"); ConstraintProto* new_ct1 = context_->working_model->add_constraints(); *new_ct1 = *ct; @@ -3286,6 +3293,15 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint( // remove coefficient, the threshold do not change! const int64_t threshold = lower_bounded ? ub_threshold : lb_threshold; + // All coeffs in [second_threshold, threshold) can be reduced to + // second_threshold. + // + // TODO(user): If 2 * min_coeff_magnitude >= bound, then the constraint can + // be completely rewriten to 2 * (enforcement_part) + sum var >= 2 which is + // what happen eventually when bound is even, but not if it is odd currently. + const int64_t second_threshold = std::max(CeilOfRatio(threshold, int64_t{2}), + threshold - min_coeff_magnitude); + // Do we only extract Booleans? // // Note that for now the default is false, and also there are problem calling @@ -3314,22 +3330,31 @@ void CpModelPresolver::ExtractEnforcementLiteralFromLinearConstraint( (only_extract_booleans && !is_boolean)) { mutable_arg->set_vars(new_size, mutable_arg->vars(i)); - // We keep this term but reduces its coeff. - // This is only for the case where only_extract_booleans == true. + int64_t new_magnitude = std::abs(arg.coeffs(i)); if (coeff > threshold) { + // We keep this term but reduces its coeff. + // This is only for the case where only_extract_booleans == true. + new_magnitude = threshold; context_->UpdateRuleStats("linear: coefficient strenghtening."); - if (lower_bounded) { - // coeff * (X - LB + LB) -> threshold * (X - LB) + coeff * LB - rhs_offset -= (coeff - threshold) * context_->MinOf(ref); - } else { - // coeff * (X - UB + UB) -> threshold * (X - UB) + coeff * UB - rhs_offset -= (coeff - threshold) * context_->MaxOf(ref); - } - mutable_arg->set_coeffs(new_size, - arg.coeffs(i) > 0 ? threshold : -threshold); - } else { - mutable_arg->set_coeffs(new_size, arg.coeffs(i)); + } else if (coeff > second_threshold && coeff < threshold) { + // This cover the special case where one big + on small is enough + // to satisfy the constraint, we can reduce the big. + new_magnitude = second_threshold; + context_->UpdateRuleStats( + "linear: advanced coefficient strenghtening."); } + if (coeff != new_magnitude) { + if (lower_bounded) { + // coeff * (X - LB + LB) -> new_magnitude * (X - LB) + coeff * LB + rhs_offset -= (coeff - new_magnitude) * context_->MinOf(ref); + } else { + // coeff * (X - UB + UB) -> new_magnitude * (X - UB) + coeff * UB + rhs_offset -= (coeff - new_magnitude) * context_->MaxOf(ref); + } + } + + mutable_arg->set_coeffs( + new_size, arg.coeffs(i) > 0 ? new_magnitude : -new_magnitude); ++new_size; continue; } diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 4565a74cfd..876ef9beb1 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -349,6 +349,7 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) { IntegerValue min_sum(0); IntegerValue max_sum(0); IntegerValue max_magnitude(0); + IntegerValue min_magnitude = kMaxIntegerValue; int new_size = 0; const int num_terms = ct->vars.size(); for (int i = 0; i < num_terms; ++i) { @@ -362,7 +363,9 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) { if (lb == ub) continue; ++new_size; - max_magnitude = std::max(max_magnitude, IntTypeAbs(coeff)); + const IntegerValue magnitude = IntTypeAbs(coeff); + max_magnitude = std::max(max_magnitude, magnitude); + min_magnitude = std::min(min_magnitude, magnitude); if (coeff > 0.0) { min_sum += coeff * lb; max_sum += coeff * ub; @@ -412,54 +415,85 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) { // TODO(user): Split constraint in two if it is boxed and there is possible // reduction? // - // TODO(user): Make sure there cannot be any overflow. They shouldn't, but - // I am not sure all the generated cuts are safe regarding min/max sum - // computation. We should check this. - if (ct->ub != kMaxIntegerValue && max_magnitude > max_sum - ct->ub) { - if (ct->lb != kMinIntegerValue) { - ++num_split_constraints_; - } else { - term_changed = true; - ++num_coeff_strenghtening_; - const int num_terms = ct->vars.size(); - const IntegerValue target = max_sum - ct->ub; - for (int i = 0; i < num_terms; ++i) { - const IntegerValue coeff = ct->coeffs[i]; - if (coeff > target) { - const IntegerVariable var = ct->vars[i]; - const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); - ct->coeffs[i] = target; - ct->ub -= (coeff - target) * ub; - } else if (coeff < -target) { - const IntegerVariable var = ct->vars[i]; - const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); - ct->coeffs[i] = -target; - ct->ub += (-target - coeff) * lb; + // TODO(user): We could cover more case of coefficient strenghtening. For + // example, if whe have 15 * X + 3 * Y >= 19, coeff of X can be reduced to 13. + if (ct->ub != kMaxIntegerValue) { + const IntegerValue threshold = max_sum - ct->ub; + const IntegerValue second_threshold = std::max( + CeilRatio(threshold, IntegerValue(2)), threshold - min_magnitude); + if (max_magnitude > second_threshold) { + if (ct->lb != kMinIntegerValue) { + ++num_split_constraints_; + } else { + term_changed = true; + ++num_coeff_strenghtening_; + const int num_terms = ct->vars.size(); + for (int i = 0; i < num_terms; ++i) { + // In all cases, we reason on a transformed constraint where the term + // is max_value - |coeff| * positive_X. If we change coeff, and + // retransform the constraint, we need to change the rhs by the + // constant term left. + const IntegerValue coeff = ct->coeffs[i]; + if (coeff > threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); + ct->coeffs[i] = threshold; + ct->ub -= (coeff - threshold) * ub; + } else if (coeff > second_threshold && coeff < threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); + ct->coeffs[i] = second_threshold; + ct->ub -= (coeff - second_threshold) * ub; + } else if (coeff < -threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); + ct->coeffs[i] = -threshold; + ct->ub -= (coeff + threshold) * lb; + } else if (coeff < -second_threshold && coeff > -threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); + ct->coeffs[i] = -second_threshold; + ct->ub -= (coeff + second_threshold) * lb; + } } } } } - if (ct->lb != kMinIntegerValue && max_magnitude > ct->lb - min_sum) { - if (ct->ub != kMaxIntegerValue) { - ++num_split_constraints_; - } else { - term_changed = true; - ++num_coeff_strenghtening_; - const int num_terms = ct->vars.size(); - const IntegerValue target = ct->lb - min_sum; - for (int i = 0; i < num_terms; ++i) { - const IntegerValue coeff = ct->coeffs[i]; - if (coeff > target) { - const IntegerVariable var = ct->vars[i]; - const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); - ct->coeffs[i] = target; - ct->lb -= (coeff - target) * lb; - } else if (coeff < -target) { - const IntegerVariable var = ct->vars[i]; - const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); - ct->coeffs[i] = -target; - ct->lb += (-target - coeff) * ub; + if (ct->lb != kMinIntegerValue) { + const IntegerValue threshold = ct->lb - min_sum; + const IntegerValue second_threshold = std::max( + CeilRatio(threshold, IntegerValue(2)), threshold - min_magnitude); + if (max_magnitude > second_threshold) { + if (ct->ub != kMaxIntegerValue) { + ++num_split_constraints_; + } else { + term_changed = true; + ++num_coeff_strenghtening_; + const int num_terms = ct->vars.size(); + for (int i = 0; i < num_terms; ++i) { + const IntegerValue coeff = ct->coeffs[i]; + if (coeff > threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); + ct->coeffs[i] = threshold; + ct->lb -= (coeff - threshold) * lb; + } else if (coeff > second_threshold && coeff < threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); + ct->coeffs[i] = second_threshold; + ct->lb -= (coeff - second_threshold) * lb; + } else if (coeff < -threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); + ct->coeffs[i] = -threshold; + ct->lb -= (coeff + threshold) * ub; + } else if (coeff < -second_threshold && coeff > -threshold) { + const IntegerVariable var = ct->vars[i]; + const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); + ct->coeffs[i] = -second_threshold; + ct->lb -= (coeff + second_threshold) * ub; + } } } }