diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 2ce249f1af..f5b49f9694 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -1287,6 +1287,11 @@ bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) { // list. // // This operation is similar to coefficient strengthening in the MIP world. +// +// TODO(user): If the constraint is boxed, split it in two if enforcement +// literals can be extracted this way. This would just be better for the CP +// propagation, and for the LP it will improve the relaxation at the cost of +// an extra row in the matrix, but this should still be better. void ExtractEnforcementLiteralFromLinearConstraint(ConstraintProto* ct, PresolveContext* context) { if (context->ModelIsUnsat()) return; diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index f9ef3337f0..dade242f69 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -1166,6 +1166,7 @@ bool IntegerTrail::EnqueueInternal( // Special case for level zero. if (integer_search_levels_.empty()) { + ++num_level_zero_enqueues_; vars_[i_lit.var].current_bound = i_lit.bound; integer_trail_[i_lit.var.value()].bound = i_lit.bound; diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index 8d477cdbfd..c2d1534825 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -692,6 +692,9 @@ class IntegerTrail : public SatPropagator { // doesn't change since we directly update the "fixed" bounds. int64 num_enqueues() const { return num_enqueues_; } + // Same as num_enqueues but only count the level zero changes. + int64 num_level_zero_enqueues() const { return num_level_zero_enqueues_; } + // All the registered bitsets will be set to one each time a LbVar is // modified. It is up to the client to clear it if it wants to be notified // with the newly modified variables. @@ -907,6 +910,7 @@ class IntegerTrail : public SatPropagator { std::vector boolean_trail_index_to_integer_one_; int64 num_enqueues_ = 0; + int64 num_level_zero_enqueues_ = 0; std::vector*> watchers_; std::vector reversible_classes_; diff --git a/ortools/sat/linear_constraint.h b/ortools/sat/linear_constraint.h index 7a9a601d80..d30e78fb19 100644 --- a/ortools/sat/linear_constraint.h +++ b/ortools/sat/linear_constraint.h @@ -58,10 +58,11 @@ struct LinearConstraint { for (int i = 0; i < vars.size(); ++i) { const IntegerValue coeff = VariableIsPositive(vars[i]) ? coeffs[i] : -coeffs[i]; - absl::StrAppend(&result, coeff.value(), "*X", vars[i].value() / 2, " "); + absl::StrAppend(&result, i > 0 ? " " : "", coeff.value(), "*X", + vars[i].value() / 2); } if (ub.value() < kMaxIntegerValue) { - absl::StrAppend(&result, "<= ", ub.value()); + absl::StrAppend(&result, " <= ", ub.value()); } return result; } diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index ff0a0754ef..a03c0dd91d 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -57,6 +57,15 @@ LinearConstraintManager::~LinearConstraintManager() { if (num_merged_constraints_ > 0) { VLOG(2) << "num_merged_constraints: " << num_merged_constraints_; } + if (num_shortened_constraints_ > 0) { + VLOG(2) << "num_shortened_constraints: " << num_shortened_constraints_; + } + if (num_splitted_constraints_ > 0) { + VLOG(2) << "num_splitted_constraints: " << num_splitted_constraints_; + } + if (num_coeff_strenghtening_ > 0) { + VLOG(2) << "num_coeff_strenghtening: " << num_coeff_strenghtening_; + } for (const auto entry : type_to_num_cuts_) { VLOG(1) << "Added " << entry.second << " cuts of type '" << entry.first << "'."; @@ -86,6 +95,7 @@ void LinearConstraintManager::RemoveMarkedConstraints() { // we regenerate identical cuts for some reason. // // TODO(user): Also compute GCD and divide by it. +// TODO(user): Call SimplifyConstraint() too. void LinearConstraintManager::Add(const LinearConstraint& ct) { LinearConstraint canonicalized = ct; const Terms terms = CanonicalizeConstraintAndGetTerms(&canonicalized); @@ -109,7 +119,7 @@ void LinearConstraintManager::Add(const LinearConstraint& ct) { ++num_merged_constraints_; } else { constraint_l2_norms_.push_back(ComputeL2Norm(canonicalized)); - equiv_constraints_[terms] = constraints_.size(); + equiv_constraints_[terms] = ConstraintIndex(constraints_.size()); constraint_is_in_lp_.push_back(false); constraint_is_cut_.push_back(false); constraint_inactive_count_.push_back(0); @@ -118,6 +128,21 @@ void LinearConstraintManager::Add(const LinearConstraint& ct) { } } +void LinearConstraintManager::RemoveConstraintFromEquivTable( + ConstraintIndex ct_index) { + const Terms terms = + CanonicalizeConstraintAndGetTerms(&constraints_[ct_index]); + equiv_constraints_.erase(terms); +} + +void LinearConstraintManager::NotifyConstraintChanged( + ConstraintIndex ct_index) { + const Terms terms = + CanonicalizeConstraintAndGetTerms(&constraints_[ct_index]); + equiv_constraints_[terms] = ct_index; + constraint_l2_norms_[ct_index] = ComputeL2Norm(constraints_[ct_index]); +} + // Same as Add(), but logs some information about the newly added constraint. // Cuts are also handled slightly differently than normal constraints. void LinearConstraintManager::AddCut( @@ -144,6 +169,140 @@ void LinearConstraintManager::AddCut( << " violation=" << violation << " eff=" << violation / l2_norm; } +void LinearConstraintManager::SimplifyConstraint(ConstraintIndex ct_index) { + LinearConstraint& ct = constraints_[ct_index]; + + IntegerValue min_sum(0); + IntegerValue max_sum(0); + IntegerValue max_magnitude(0); + int new_size = 0; + const int num_terms = ct.vars.size(); + for (int i = 0; i < num_terms; ++i) { + const IntegerVariable var = ct.vars[i]; + const IntegerValue coeff = ct.coeffs[i]; + const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); + const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); + + // For now we do not change ct, but just compute its new_size if we where + // to remove a fixed term. + if (lb == ub) continue; + ++new_size; + + max_magnitude = std::max(max_magnitude, IntTypeAbs(coeff)); + if (coeff > 0.0) { + min_sum += coeff * lb; + max_sum += coeff * ub; + } else { + min_sum += coeff * ub; + max_sum += coeff * lb; + } + } + + // Shorten the constraint if needed. + if (new_size < num_terms) { + ++num_shortened_constraints_; + RemoveConstraintFromEquivTable(ct_index); + + new_size = 0; + for (int i = 0; i < num_terms; ++i) { + const IntegerVariable var = ct.vars[i]; + const IntegerValue coeff = ct.coeffs[i]; + const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var); + const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var); + if (lb == ub) { + const IntegerValue rhs_adjust = lb * coeff; + if (ct.lb > kMinIntegerValue) ct.lb -= rhs_adjust; + if (ct.ub < kMaxIntegerValue) ct.ub -= rhs_adjust; + continue; + } + ct.vars[new_size] = var; + ct.coeffs[new_size] = coeff; + ++new_size; + } + ct.vars.resize(new_size); + ct.coeffs.resize(new_size); + + NotifyConstraintChanged(ct_index); + } + + // Relax the bound if needed, note that this doesn't require a change to + // the equiv map. + if (min_sum >= ct.lb) ct.lb = kMinIntegerValue; + if (max_sum <= ct.ub) ct.ub = kMaxIntegerValue; + + // Clear constraints that are always true. + // We rely on the deletion code to remove them eventually. + if (ct.lb == kMinIntegerValue && ct.ub == kMaxIntegerValue) { + RemoveConstraintFromEquivTable(ct_index); + ct.vars.clear(); + ct.coeffs.clear(); + constraint_l2_norms_[ct_index] = 0.0; + return; + } + + // 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 should'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_splitted_constraints_; + } else { + ++num_coeff_strenghtening_; + RemoveConstraintFromEquivTable(ct_index); + + 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; + } + } + + NotifyConstraintChanged(ct_index); + } + } + + if (ct.lb != kMinIntegerValue && max_magnitude > ct.lb - min_sum) { + if (ct.ub != kMaxIntegerValue) { + ++num_splitted_constraints_; + } else { + ++num_coeff_strenghtening_; + RemoveConstraintFromEquivTable(ct_index); + + 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; + } + } + + NotifyConstraintChanged(ct_index); + } + } +} + bool LinearConstraintManager::ChangeLp( const gtl::ITIVector& lp_solution) { const int old_num_constraints = lp_constraints_.size(); @@ -152,12 +311,19 @@ bool LinearConstraintManager::ChangeLp( std::vector new_constraints_efficacies; std::vector new_constraints_orthogonalities; + const bool simplify_constraints = + integer_trail_.num_level_zero_enqueues() > last_simplification_timestamp_; + last_simplification_timestamp_ = integer_trail_.num_level_zero_enqueues(); + // We keep any constraints that is already present, and otherwise, we add the // ones that are currently not satisfied by at least "tolerance". const double tolerance = 1e-6; for (ConstraintIndex i(0); i < constraints_.size(); ++i) { if (constraint_permanently_removed_[i]) continue; + // Inprocessing of the constraint. + if (simplify_constraints) SimplifyConstraint(i); + // TODO(user,user): Use constraint status from GLOP for constraints // already in LP. const double activity = ComputeActivity(constraints_[i], lp_solution); @@ -179,6 +345,7 @@ bool LinearConstraintManager::ChangeLp( new_constraints_orthogonalities.push_back(1.0); } } + // Remove constraints in a batch to avoid changing the LP too frequently. if (constraints_removal_list_.size() >= sat_parameters_.constraint_removal_batch_size()) { diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h index 9a27832dae..3b63ac14f9 100644 --- a/ortools/sat/linear_constraint_manager.h +++ b/ortools/sat/linear_constraint_manager.h @@ -39,7 +39,9 @@ namespace sat { // which constraint should go into the current LP. class LinearConstraintManager { public: - LinearConstraintManager() {} + explicit LinearConstraintManager(Model* model) + : sat_parameters_(*model->GetOrCreate()), + integer_trail_(*model->GetOrCreate()) {} ~LinearConstraintManager(); // Add a new constraint to the manager. Note that we canonicalize constraints @@ -75,19 +77,41 @@ class LinearConstraintManager { return lp_constraints_; } - void SetParameters(const SatParameters& params) { sat_parameters_ = params; } - int num_cuts() const { return num_cuts_; } + int64 num_shortened_constraints() const { return num_shortened_constraints_; } + int64 num_coeff_strenghtening() const { return num_coeff_strenghtening_; } + + // Visible for testing. + // + // Apply basic inprocessing simplification rules: + // - remove fixed variable + // - reduce large coefficient (i.e. coeff strenghtenning or big-M reduction). + // This uses level-zero bounds. + void SimplifyConstraint(ConstraintIndex index); private: // Removes the marked constraints from the LP. void RemoveMarkedConstraints(); - SatParameters sat_parameters_; + // Important: before modifying a constraint in equiv_constraints_, one + // must call RemoveConstraintFromEquivTable() and when done, call + // NotifyConstraintChanged(). + // + // TODO(user): That might be a bit slow and brittle. Redesign the equiv map + // so that we compare actual constraints (so we never merge different ones), + // and we don't need to keep duplicate of the constraints in memory. + void RemoveConstraintFromEquivTable(ConstraintIndex ct_index); + void NotifyConstraintChanged(ConstraintIndex ct_index); + + const SatParameters& sat_parameters_; + const IntegerTrail& integer_trail_; // Set at true by Add() and at false by ChangeLp(). bool some_lp_constraint_bounds_changed_ = false; + // Optimization to avoid calling SimplifyConstraint() when not needed. + int64 last_simplification_timestamp_ = -1; + // TODO(user): Merge all the constraint related info in a struct and store // a vector of struct instead. The global list of constraint. gtl::ITIVector constraints_; @@ -123,8 +147,12 @@ class LinearConstraintManager { return hash; } }; - absl::flat_hash_map equiv_constraints_; + absl::flat_hash_map equiv_constraints_; + int64 num_merged_constraints_ = 0; + int64 num_shortened_constraints_ = 0; + int64 num_splitted_constraints_ = 0; + int64 num_coeff_strenghtening_ = 0; int num_cuts_ = 0; std::map type_to_num_cuts_; diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index f00b4a7371..6646cd3c44 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -42,7 +42,8 @@ const double LinearProgrammingConstraint::kLpEpsilon = 1e-6; // TODO(user): make SatParameters singleton too, otherwise changing them after // a constraint was added will have no effect on this class. LinearProgrammingConstraint::LinearProgrammingConstraint(Model* model) - : sat_parameters_(*(model->GetOrCreate())), + : constraint_manager_(model), + sat_parameters_(*(model->GetOrCreate())), model_(model), time_limit_(model->GetOrCreate()), integer_trail_(model->GetOrCreate()), @@ -116,7 +117,6 @@ void LinearProgrammingConstraint::CreateLpFromConstraintManager() { // Fill integer_lp_. integer_lp_.clear(); infinity_norms_.clear(); - constraint_manager_.SetParameters(sat_parameters_); const auto& all_constraints = constraint_manager_.AllConstraints(); for (const auto index : constraint_manager_.LpConstraints()) { const LinearConstraint& ct = all_constraints[index]; @@ -1168,9 +1168,12 @@ bool LinearProgrammingConstraint::ExactLpReasonning() { gtl::ITIVector reduced_costs; IntegerValue rc_ub; - CHECK(ComputeNewLinearConstraint( - /*use_constraint_status=*/false, integer_multipliers, &reduced_costs, - &rc_ub)); + if (!ComputeNewLinearConstraint( + /*use_constraint_status=*/false, integer_multipliers, &reduced_costs, + &rc_ub)) { + VLOG(1) << "Issue while computing the exact LP reason. Aborting."; + return true; + } // The "objective constraint" behave like if the unscaled cp multiplier was // 1.0, so we will multiply it by this number and add it to reduced_costs. @@ -1215,9 +1218,12 @@ bool LinearProgrammingConstraint::FillExactDualRayReason() { gtl::ITIVector dense_new_constraint; IntegerValue new_constraint_ub; - CHECK(ComputeNewLinearConstraint( - /*use_constraint_status=*/false, integer_multipliers, - &dense_new_constraint, &new_constraint_ub)); + if (!ComputeNewLinearConstraint( + /*use_constraint_status=*/false, integer_multipliers, + &dense_new_constraint, &new_constraint_ub)) { + VLOG(1) << "Isse while computing the exact dual ray reason. Aborting."; + return false; + } AdjustNewLinearConstraint(&integer_multipliers, &dense_new_constraint, &new_constraint_ub);