diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index ca94bbc6de..7895d6016d 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -43,7 +43,7 @@ void ExpandReservoir(ConstraintProto* ct, PresolveContext* context) { auto is_optional = [&context, &reservoir](int index) { if (reservoir.actives_size() == 0) return false; const int literal = reservoir.actives(index); - return !context->DomainOf(PositiveRef(literal)).IsFixed(); + return !context->IsFixed(literal); }; const int true_literal = context->GetOrCreateConstantVar(1); auto active = [&reservoir, true_literal](int index) { diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 77ed1073cf..c709fc3daa 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -82,22 +82,48 @@ LinearConstraintManager::~LinearConstraintManager() { } } -// TODO(user,user): Also update the revised simplex starting basis for the -// next solve. -void LinearConstraintManager::RemoveMarkedConstraints() { - int new_index = 0; - for (int i = 0; i < lp_constraints_.size(); ++i) { - const ConstraintIndex constraint = lp_constraints_[i]; - if (constraints_removal_list_.contains(constraint)) { - constraint_infos_[constraint].is_in_lp = false; - continue; +bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints( + glop::BasisState* solution_state) { + if (solution_state->IsEmpty()) return false; // Mainly to simplify tests. + const glop::RowIndex num_rows(lp_constraints_.size()); + const glop::ColIndex num_cols = + solution_state->statuses.size() - RowToColIndex(num_rows); + + int new_size = 0; + for (int i = 0; i < num_rows; ++i) { + const ConstraintIndex constraint_index = lp_constraints_[i]; + + // Constraints that are not tight in the current solution have a basic + // status. We remove the ones that have been inactive in the last recent + // solves. + // + // TODO(user): More advanced heuristics might perform better, I didn't do + // a lot of tuning experiments yet. + const glop::VariableStatus row_status = + solution_state->statuses[num_cols + glop::ColIndex(i)]; + if (row_status == glop::VariableStatus::BASIC) { + constraint_infos_[constraint_index].inactive_count++; + if (constraint_infos_[constraint_index].inactive_count > + sat_parameters_.max_consecutive_inactive_count()) { + constraint_infos_[constraint_index].is_in_lp = false; + continue; // Remove it. + } + } else { + // Only count consecutive inactivities. + constraint_infos_[constraint_index].inactive_count = 0; } - lp_constraints_[new_index] = constraint; - new_index++; + + lp_constraints_[new_size] = constraint_index; + solution_state->statuses[num_cols + glop::ColIndex(new_size)] = row_status; + new_size++; } - lp_constraints_.resize(new_index); - VLOG(2) << "Removed " << constraints_removal_list_.size() << " constraints."; - constraints_removal_list_.clear(); + const int num_removed_constraints = lp_constraints_.size() - new_size; + lp_constraints_.resize(new_size); + solution_state->statuses.resize(num_cols + glop::ColIndex(new_size)); + if (num_removed_constraints > 0) { + VLOG(2) << "Removed " << num_removed_constraints << " constraints"; + } + return num_removed_constraints > 0; } // Because sometimes we split a == constraint in two (>= and <=), it makes sense @@ -343,7 +369,8 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) { } bool LinearConstraintManager::ChangeLp( - const gtl::ITIVector& lp_solution) { + const gtl::ITIVector& lp_solution, + glop::BasisState* solution_state) { VLOG(3) << "Enter ChangeLP, scan " << constraint_infos_.size() << " constraints"; std::vector new_constraints; @@ -375,25 +402,16 @@ bool LinearConstraintManager::ChangeLp( equiv_constraints_[constraint_infos_[i].hash] = i; } - // TODO(user,user): Use constraint status from GLOP for constraints - // already in LP. + if (constraint_infos_[i].is_in_lp) continue; + const double activity = ComputeActivity(constraint_infos_[i].constraint, lp_solution); - const double lb_violation = ToDouble(constraint_infos_[i].constraint.lb) - activity; const double ub_violation = activity - ToDouble(constraint_infos_[i].constraint.ub); const double violation = std::max(lb_violation, ub_violation); - if (constraint_infos_[i].is_in_lp && violation < tolerance) { - constraint_infos_[i].inactive_count++; - if (constraint_infos_[i].is_cut && - constraint_infos_[i].inactive_count > - sat_parameters_.max_inactive_count()) { - // Mark cut for removal. - constraints_removal_list_.insert(i); - } - } else if (!constraint_infos_[i].is_in_lp && violation >= tolerance) { + if (violation >= tolerance) { constraint_infos_[i].inactive_count = 0; new_constraints.push_back(i); new_constraints_efficacies.push_back(violation / @@ -413,10 +431,11 @@ bool LinearConstraintManager::ChangeLp( } } - // Remove constraints in a batch to avoid changing the LP too frequently. - if (constraints_removal_list_.size() >= - sat_parameters_.constraint_removal_batch_size()) { - RemoveMarkedConstraints(); + // Remove constraints from the current LP that have been inactive for a while. + // We do that after we computed new_constraints so we do not need to iterate + // over the just deleted constraints. + if (MaybeRemoveSomeInactiveConstraints(solution_state)) { + current_lp_is_changed_ = true; } // Note that the algo below is in O(limit * new_constraint). In order to @@ -430,7 +449,8 @@ bool LinearConstraintManager::ChangeLp( // TODO(user): This blowup factor could be adaptative w.r.t. the constraint // limit. const int kBlowupFactor = 4; - int constraint_limit = std::min(50, static_cast(new_constraints.size())); + int constraint_limit = std::min(sat_parameters_.new_constraints_batch_size(), + static_cast(new_constraints.size())); if (lp_constraints_.empty()) { constraint_limit = std::min(1000, static_cast(new_constraints.size())); } @@ -448,6 +468,7 @@ bool LinearConstraintManager::ChangeLp( new_constraints.resize(kBlowupFactor * constraint_limit); } + int num_added = 0; int num_skipped_checks = 0; const int kCheckFrequency = 100; ConstraintIndex last_added_candidate = kInvalidConstraintIndex; @@ -507,13 +528,20 @@ bool LinearConstraintManager::ChangeLp( // Note that it is important for LP incremental solving that the old // constraints stays at the same position in this list (and thus in the // returned GetLp()). + ++num_added; current_lp_is_changed_ = true; lp_constraints_.push_back(best_candidate); last_added_candidate = best_candidate; } } - VLOG(3) << " - Exit ChangeLP"; + if (num_added > 0) { + // We update the solution sate to match the new LP size. + VLOG(2) << "Added " << num_added << " constraints."; + solution_state->statuses.resize(solution_state->statuses.size() + num_added, + glop::VariableStatus::BASIC); + } + // The LP changed only if we added new constraints or if some constraints // already inside changed (simplification or tighter bounds). if (current_lp_is_changed_) { diff --git a/ortools/sat/linear_constraint_manager.h b/ortools/sat/linear_constraint_manager.h index 0724116eee..785581d817 100644 --- a/ortools/sat/linear_constraint_manager.h +++ b/ortools/sat/linear_constraint_manager.h @@ -19,6 +19,7 @@ #include "absl/container/flat_hash_map.h" #include "absl/container/flat_hash_set.h" +#include "ortools/glop/revised_simplex.h" #include "ortools/sat/linear_constraint.h" #include "ortools/sat/model.h" #include "ortools/sat/sat_parameters.pb.h" @@ -81,8 +82,13 @@ class LinearConstraintManager { // should usually be the optimal solution of the LP returned by GetLp() before // this call, but is just used as an heuristic. // + // The current solution state is used for detecting inactive constraints. It + // is also updated correctly on constraint deletion/addition so that the + // simplex can be fully iterative on restart by loading this modified state. + // // Returns true iff LpConstraints() will return a different LP than before. - bool ChangeLp(const gtl::ITIVector& lp_solution); + bool ChangeLp(const gtl::ITIVector& lp_solution, + glop::BasisState* solution_state); // This can be called initially to add all the current constraint to the LP // returned by GetLp(). @@ -105,8 +111,16 @@ class LinearConstraintManager { int64 num_coeff_strenghtening() const { return num_coeff_strenghtening_; } private: - // Removes the marked constraints from the LP. - void RemoveMarkedConstraints(); + // Heuristic that decide which constraints we should remove from the current + // LP. Note that such constraints can be added back later by the heuristic + // responsible for adding new constraints from the pool. + // + // Returns true iff one or more constraints where removed. + // + // If the solutions_state is empty, then this function does nothing and + // returns false (this is used for tests). Otherwise, the solutions_state is + // assumed to correspond to the current LP and to be of the correct size. + bool MaybeRemoveSomeInactiveConstraints(glop::BasisState* solution_state); // Apply basic inprocessing simplification rules: // - remove fixed variable @@ -133,10 +147,6 @@ class LinearConstraintManager { gtl::ITIVector constraint_infos_; - // Temporary list of constraints marked for removal. Note that we remove - // constraints in batch to avoid changing LP too frequently. - absl::flat_hash_set constraints_removal_list_; - // The subset of constraints currently in the lp. std::vector lp_constraints_; diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 70ec6e60fd..756086124e 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -886,15 +886,10 @@ bool LinearProgrammingConstraint::Propagate() { // this beeing called again on the next IncrementalPropagate() call, but that // might not always happen at level zero. if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) { - // First add any new lazy constraints or cuts that where previsouly - // generated and are now cutting the current solution. - if (constraint_manager_.ChangeLp(expanded_lp_solution_)) { - CreateLpFromConstraintManager(); - if (!SolveLp()) return true; - } else if (constraint_manager_.num_cuts() < - sat_parameters_.max_num_cuts()) { - const int old_num_cuts = constraint_manager_.num_cuts(); - + // We wait for the first batch of problem constraints to be added before we + // begin to generate cuts. + if (!integer_lp_.empty() && + constraint_manager_.num_cuts() < sat_parameters_.max_num_cuts()) { // The "generic" cuts are currently part of this class as they are using // data from the current LP. // @@ -912,22 +907,23 @@ bool LinearProgrammingConstraint::Propagate() { generator.generate_cuts(expanded_lp_solution_, &constraint_manager_); } } + } - if (constraint_manager_.num_cuts() > old_num_cuts && - constraint_manager_.ChangeLp(expanded_lp_solution_)) { - CreateLpFromConstraintManager(); - const double old_obj = simplex_.GetObjectiveValue(); - if (!SolveLp()) return true; - if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) { - VLOG(1) << "Cuts relaxation improvement " << old_obj << " -> " - << simplex_.GetObjectiveValue() - << " diff: " << simplex_.GetObjectiveValue() - old_obj - << " level: " << trail_->CurrentDecisionLevel(); - } - } else { - if (trail_->CurrentDecisionLevel() == 0) { - lp_at_level_zero_is_final_ = true; - } + glop::BasisState state = simplex_.GetState(); + if (constraint_manager_.ChangeLp(expanded_lp_solution_, &state)) { + simplex_.LoadStateForNextSolve(state); + CreateLpFromConstraintManager(); + const double old_obj = simplex_.GetObjectiveValue(); + if (!SolveLp()) return true; + if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) { + VLOG(1) << "Relaxation improvement " << old_obj << " -> " + << simplex_.GetObjectiveValue() + << " diff: " << simplex_.GetObjectiveValue() - old_obj + << " level: " << trail_->CurrentDecisionLevel(); + } + } else { + if (trail_->CurrentDecisionLevel() == 0) { + lp_at_level_zero_is_final_ = true; } } } @@ -1408,41 +1404,41 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint( for (const auto entry : integer_lp_[row].terms) { const ColIndex col = entry.first; const IntegerValue coeff = entry.second; + const IntegerValue abs_coef = IntTypeAbs(coeff); CHECK_NE(coeff, 0); + const IntegerVariable var = integer_variables_[col.value()]; + const IntegerValue lb = integer_trail_->LowerBound(var); + const IntegerValue ub = integer_trail_->UpperBound(var); + // Moving a variable away from zero seems to improve the bound even // if it reduces the number of non-zero. Note that this is because of // this that positive_diff and negative_diff are not the same. const IntegerValue current = (*dense_terms)[col]; if (current == 0) { const IntegerValue overflow_limit( - FloorRatio(kMaxWantedCoeff, IntTypeAbs(coeff))); + FloorRatio(kMaxWantedCoeff, abs_coef)); positive_limit = std::min(positive_limit, overflow_limit); negative_limit = std::min(negative_limit, overflow_limit); - const IntegerVariable var = integer_variables_[col.value()]; if (coeff > 0) { - positive_diff -= coeff * integer_trail_->LowerBound(var); - negative_diff -= coeff * integer_trail_->UpperBound(var); + positive_diff -= coeff * lb; + negative_diff -= coeff * ub; } else { - positive_diff -= coeff * integer_trail_->UpperBound(var); - negative_diff -= coeff * integer_trail_->LowerBound(var); + positive_diff -= coeff * ub; + negative_diff -= coeff * lb; } continue; } - // We don't want to change the sign of current or to have an overflow. - IntegerValue before_sign_change( - FloorRatio(IntTypeAbs(current), IntTypeAbs(coeff))); - - // If the variable is fixed, we don't actually care about changing the - // sign. - const IntegerVariable var = integer_variables_[col.value()]; - if (integer_trail_->LowerBound(var) == integer_trail_->UpperBound(var)) { - before_sign_change = kMaxWantedCoeff; + // We don't want to change the sign of current (except if the variable is + // fixed) or to have an overflow. + IntegerValue before_sign_change = kMaxWantedCoeff; + if (lb != ub) { + before_sign_change = FloorRatio(IntTypeAbs(current), abs_coef); } const IntegerValue overflow_limit( - FloorRatio(kMaxWantedCoeff - IntTypeAbs(current), IntTypeAbs(coeff))); + FloorRatio(kMaxWantedCoeff - IntTypeAbs(current), abs_coef)); if (current > 0 == coeff > 0) { // Same sign. negative_limit = std::min(negative_limit, before_sign_change); positive_limit = std::min(positive_limit, overflow_limit); @@ -1452,12 +1448,11 @@ void LinearProgrammingConstraint::AdjustNewLinearConstraint( } // This is how diff change. - const IntegerValue implied = current > 0 - ? integer_trail_->LowerBound(var) - : integer_trail_->UpperBound(var); - - positive_diff -= coeff * implied; - negative_diff -= coeff * implied; + const IntegerValue implied = current > 0 ? lb : ub; + if (implied != 0) { + positive_diff -= coeff * implied; + negative_diff -= coeff * implied; + } } // Only add a multiple of this row if it tighten the final constraint. @@ -1553,6 +1548,11 @@ bool LinearProgrammingConstraint::ExactLpReasonning() { IntegerSumLE* cp_constraint = new IntegerSumLE({}, new_constraint.vars, new_constraint.coeffs, new_constraint.ub, model_); + if (trail_->CurrentDecisionLevel() == 0) { + // Since we will never ask the reason for a constraint at level 0, we just + // keep the last one. + optimal_constraints_.clear(); + } optimal_constraints_.emplace_back(cp_constraint); rev_optimal_constraints_size_ = optimal_constraints_.size(); return cp_constraint->Propagate(); diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index a2e50d525b..be86ca61a6 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -590,13 +590,14 @@ message SatParameters { // feature. optional double min_orthogonality_for_lp_constraints = 115 [default = 0.0]; - // If a constraint/cut in LP is not active for the 'max_inactive_count', - // remove it from the LP. - optional int64 max_inactive_count = 121 [default = 1000]; + // If a constraint/cut in LP is not active for that many consecutive OPTIMAL + // solves, remove it from the LP. Note that it might be added again later if + // it become violated by the current LP solution. + optional int32 max_consecutive_inactive_count = 121 [default = 100]; - // Remove constraints only if at least 'constraint_removal_batch_size' - // constraints are marked for removal. - optional int64 constraint_removal_batch_size = 122 [default = 100]; + // Add that many lazy contraints (or cuts) at once in the LP. Note that at the + // beginning of the solve, we do add more than this. + optional int32 new_constraints_batch_size = 122 [default = 50]; // The search branching will be used to decide how to branch on unfixed nodes. enum SearchBranching {