From d0081e54093d74a1962346f2dcaa8705e7b6c57a Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Fri, 27 Sep 2019 17:01:54 +0200 Subject: [PATCH] more work on CP-SAT internals --- ortools/sat/cp_model_expand.cc | 47 ++++----- ortools/sat/cp_model_loader.cc | 8 ++ ortools/sat/cp_model_presolve.cc | 115 ++++++++++++++++++++++- ortools/sat/cp_model_solver.cc | 3 +- ortools/sat/integer.cc | 12 ++- ortools/sat/integer.h | 3 + ortools/sat/linear_constraint_manager.cc | 3 +- ortools/sat/linear_relaxation.cc | 8 +- ortools/sat/presolve_context.cc | 110 ++++++++++++++-------- ortools/sat/presolve_context.h | 2 +- ortools/sat/presolve_util.cc | 85 +++++++++++++++++ ortools/sat/presolve_util.h | 6 ++ ortools/sat/sat_parameters.proto | 11 ++- ortools/sat/synchronization.h | 6 +- 14 files changed, 344 insertions(+), 75 deletions(-) diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index 6d0d4f2401..3dec12d65d 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -446,12 +446,9 @@ void ExpandElement(ConstraintProto* ct, PresolveContext* context) { } const int64 value = var_domain.Min(); - - if (!gtl::ContainsKey(constant_var_values_usage, value)) { + if (constant_var_values_usage[value]++ == 0) { constant_var_values.push_back(value); } - - constant_var_values_usage[value]++; } } @@ -489,12 +486,13 @@ void ExpandElement(ConstraintProto* ct, PresolveContext* context) { for (const ClosedInterval& interval : target_domain) { for (int64 v = interval.start; v <= interval.end; ++v) { - if (constant_var_values_usage[v] > 1) { + const int usage = gtl::FindOrDie(constant_var_values_usage, v); + if (usage > 1) { const int lit = context->GetOrCreateVarValueEncoding(target_ref, v); - CHECK(gtl::ContainsKey(constant_var_values_usage, v)); - supports[v] = + BoolArgumentProto* const support = context->working_model->add_constraints()->mutable_bool_or(); - supports[v]->add_literals(NegatedRef(lit)); + supports[v] = support; + support->add_literals(NegatedRef(lit)); } } } @@ -518,7 +516,7 @@ void ExpandElement(ConstraintProto* ct, PresolveContext* context) { context->AddImplyInDomain(index_lit, var, Domain(v)); } else if (target_domain.Size() == 1) { // TODO(user): If we know all variables are different, then this - // becomes an equivalence. + // becomes an equivalence. context->AddImplyInDomain(index_lit, var, target_domain); } else if (var_domain.Size() == 1) { if (all_constants) { @@ -529,8 +527,7 @@ void ExpandElement(ConstraintProto* ct, PresolveContext* context) { const int target_lit = context->GetOrCreateVarValueEncoding(target_ref, value); context->AddImplication(index_lit, target_lit); - BoolArgumentProto* const support = gtl::FindOrDie(supports, value); - support->add_literals(index_lit); + gtl::FindOrDie(supports, value)->add_literals(index_lit); } else { // Try to reuse the literal of the index. context->InsertVarValueEncoding(index_lit, target_ref, value); @@ -555,26 +552,30 @@ void ExpandElement(ConstraintProto* ct, PresolveContext* context) { const int64 var_min = target_domain.Min(); // Scan all values to find the one with the most literals attached. - int64 value = kint64max; + int64 most_frequent_value = kint64max; int usage = -1; for (const auto it : constant_var_values_usage) { - if (it.second > usage || (it.second == usage && it.first < value)) { + if (it.second > usage || + (it.second == usage && it.first < most_frequent_value)) { usage = it.second; - value = it.first; + most_frequent_value = it.first; } } - if (value != var_min) { - VLOG(3) << "expand element: choose " << value << " with usage " << usage - << " over " << var_min << " among " << size << " values."; - } - // Add a linear constraint. This helps the linear relaxation. // // We try to minimize the size of the linear constraint (if the gain is // meaningful compared to using the min that has the advantage that all // coefficients are positive). - const int64 base = usage > 2 && usage > size / 10 ? value : var_min; + // TODO(user): Benchmark if using base is always beneficial. + // TODO(user): Try not to create this if max_usage == 1. + const int64 base = + usage > 2 && usage > size / 10 ? most_frequent_value : var_min; + if (base != var_min) { + VLOG(3) << "expand element: choose " << base << " with usage " << usage + << " over " << var_min << " among " << size << " values."; + } + LinearConstraintProto* const linear = context->working_model->add_constraints()->mutable_linear(); int64 rhs = -base; @@ -582,11 +583,11 @@ void ExpandElement(ConstraintProto* ct, PresolveContext* context) { linear->add_coeffs(-1); for (const ClosedInterval& interval : index_domain) { for (int64 v = interval.start; v <= interval.end; ++v) { - const int var = element.vars(v); + const int ref = element.vars(v); const int index_lit = context->GetOrCreateVarValueEncoding(index_ref, v); - const int64 delta = context->DomainOf(var).Min() - base; - if (index_lit >= 0) { + const int64 delta = context->DomainOf(ref).Min() - base; + if (RefIsPositive(index_lit)) { linear->add_vars(index_lit); linear->add_coeffs(delta); } else { diff --git a/ortools/sat/cp_model_loader.cc b/ortools/sat/cp_model_loader.cc index 04d99c969b..36b589f835 100644 --- a/ortools/sat/cp_model_loader.cc +++ b/ortools/sat/cp_model_loader.cc @@ -1050,6 +1050,14 @@ void LoadEquivalenceNeqAC(const std::vector enforcement_literal, } // namespace void LoadLinearConstraint(const ConstraintProto& ct, Model* m) { + if (ct.linear().vars().empty()) { + const Domain rhs = ReadDomainFromProto(ct.linear()); + if (rhs.Contains(0)) return; + VLOG(1) << "Trivially UNSAT constraint: " << ct.DebugString(); + m->GetOrCreate()->NotifyThatModelIsUnsat(); + return; + } + auto* mapping = m->GetOrCreate(); auto* integer_trail = m->GetOrCreate(); auto* encoder = m->GetOrCreate(); diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 5d886f505c..d16fb51a35 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -1010,6 +1010,8 @@ bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) { } FillDomainInProto(rhs, ct->mutable_linear()); + const bool was_affine = context_->affine_constraints.contains(ct); + // Propagate the variable bounds. if (ct->enforcement_literal().size() <= 1) { bool new_bounds = false; @@ -1017,11 +1019,12 @@ bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) { Domain right_domain(0, 0); term_domains[num_vars] = rhs.Negation(); for (int i = num_vars - 1; i >= 0; --i) { + const int64 var_coeff = ct->linear().coeffs(i); right_domain = right_domain.AdditionWith(term_domains[i + 1]).RelaxIfTooComplex(); new_domain = left_domains[i] .AdditionWith(right_domain) - .InverseMultiplicationBy(-ct->linear().coeffs(i)); + .InverseMultiplicationBy(-var_coeff); if (ct->enforcement_literal().size() == 1) { // We cannot push the new domain, but we can add some deduction. @@ -1042,6 +1045,94 @@ bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) { if (domain_modified) { new_bounds = true; } + + // Can we perform some substitution? + // + // TODO(user): there is no guarante we will not miss some since we might + // not reprocess a constraint once other have been deleted. + + // Skip affine constraint. Note however that we are before the affine + // detection for this constraint. If this variable only appear in linear + // constraints, the final result should be the same though. It should be + // even better since we don't need to track any affine relations because + // we remove the variable usage right away. + if (was_affine) continue; + + // TODO(user): We actually do not need a strict equality when + // keep_all_feasible_solutions is false, but that simplifies things as the + // SubstituteVariable() function cannot fail this way. + if (rhs.Min() != rhs.Max()) continue; + + // Only consider "implied free" variables. Note that the coefficient of + // magnitude 1 is important otherwise we can't easily remove the + // constraint since the fact that the sum of the other terms must be a + // multiple of coeff will not be enforced anymore. + const int var = ct->linear().vars(i); + if (context_->DomainOf(var) != new_domain) continue; + if (std::abs(var_coeff) != 1) continue; + + // Skip if the variable is in the objective. + if (context_->var_to_constraints[var].contains(-1)) continue; + + // Only consider low degree columns. + if (context_->var_to_constraints[var].size() < 2) continue; + if (context_->var_to_constraints[var].size() > + options_.parameters.presolve_substitution_level()) { + continue; + } + + // Check pre-conditions on all the constraints in which this variable + // appear. Basically they must all be linear and not already used for + // affine relation. + std::vector others; + for (const int c : context_->var_to_constraints[var]) { + if (context_->working_model->mutable_constraints(c) == ct) continue; + others.push_back(c); + } + bool abort = false; + for (const int c : others) { + if (context_->working_model->constraints(c).constraint_case() != + ConstraintProto::ConstraintCase::kLinear) { + abort = true; + break; + } + if (context_->affine_constraints.contains( + &context_->working_model->constraints(c))) { + abort = true; + break; + } + for (const int ref : + context_->working_model->constraints(c).enforcement_literal()) { + if (PositiveRef(ref) == var) { + abort = true; + break; + } + } + if (abort) break; + } + if (abort) continue; + + // Do the actual substitution. + for (const int c : others) { + // TODO(user): In some corner cases, this might create integer overflow + // issues. The danger is limited since the range of the linear + // expression used in the definition do not exceed the domain of the + // variable we substitute. + SubstituteVariable(var, var_coeff, *ct, + context_->working_model->mutable_constraints(c)); + + // TODO(user): We should re-enqueue these constraints for presolve. + context_->UpdateConstraintVariableUsage(c); + } + + context_->UpdateRuleStats( + absl::StrCat("linear: variable substitution ", others.size())); + + // The variable now only appear in its definition and we can remove it + // because it was implied free. + CHECK_EQ(context_->var_to_constraints[var].size(), 1); + *context_->mapping_model->add_constraints() = *ct; + return RemoveConstraint(ct); } if (new_bounds) { context_->UpdateRuleStats("linear: reduced variable domains"); @@ -1052,7 +1143,6 @@ bool CpModelPresolver::PresolveLinear(ConstraintProto* ct) { // // TODO(user): it might be better to first add only the affine relation with // a coefficient of magnitude 1, and later the one with larger coeffs. - const bool was_affine = gtl::ContainsKey(context_->affine_constraints, ct); if (!was_affine && !HasEnforcementLiteral(*ct)) { const LinearConstraintProto& arg = ct->linear(); const int64 rhs_min = rhs.Min(); @@ -3653,6 +3743,22 @@ void CpModelPresolver::TryToSimplifyDomains() { context_->UpdateNewConstraintsVariableUsage(); } +int ScoreConstraint(PresolveContext* context, int c) { + const int num_constraints = context->working_model->constraints_size(); + const ConstraintProto& ct = context->working_model->constraints(c); + if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinear) { + if (ct.enforcement_literal_size() == 0 && ct.linear().vars_size() == 2 && + ct.linear().domain_size() == 2 && + ct.linear().domain(0) == ct.linear().domain(1)) { + return c; + } else if (ct.enforcement_literal_size() == 1 && + ct.linear().vars_size() == 1) { + return num_constraints + c; + } + } + return num_constraints * 2 + c; +} + void CpModelPresolver::PresolveToFixPoint() { if (context_->ModelIsUnsat()) return; @@ -3666,6 +3772,9 @@ void CpModelPresolver::PresolveToFixPoint() { std::vector in_queue(context_->working_model->constraints_size(), true); std::deque queue(context_->working_model->constraints_size()); std::iota(queue.begin(), queue.end(), 0); + std::sort(queue.begin(), queue.end(), [this](int a, int b) { + return ScoreConstraint(context_, a) < ScoreConstraint(context_, b); + }); while (!queue.empty() && !context_->ModelIsUnsat()) { if (time_limit != nullptr && time_limit->LimitReached()) break; while (!queue.empty() && !context_->ModelIsUnsat()) { @@ -3865,6 +3974,8 @@ void CpModelPresolver::RemoveUnusedEquivalentVariables() { arg->add_domain(r.offset); } + VLOG(1) << "num_affine_relations kept = " << num_affine_relations; + // Update the variable usage. context_->UpdateNewConstraintsVariableUsage(); } diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 58deda23ba..0a7302c3be 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -1858,7 +1858,8 @@ class FullProblemSolver : public SubSolver { bool solving_first_chunk_ = true; absl::Mutex mutex_; - double deterministic_time_since_last_synchronize_ GUARDED_BY(mutex_) = 0.0; + double deterministic_time_since_last_synchronize_ GUARDED_BY(mutex_) = + 0.0; bool previous_task_is_completed_ GUARDED_BY(mutex_) = true; }; diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index 1566a503e5..ea7cabe5e5 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -43,11 +43,21 @@ void IntegerEncoder::FullyEncodeVariable(IntegerVariable var) { // TODO(user): Maybe we can optimize the literal creation order and their // polarity as our default SAT heuristics initially depends on this. + // + // TODO(user): Currently, in some corner cases, + // GetOrCreateLiteralAssociatedToEquality() might trigger some propagation + // that update the domain of var, so we need to cache the values to not read + // garbage. Note that it is okay to call the function on values no longer + // reachable, as this will just do nothing. + tmp_values_.clear(); for (const ClosedInterval interval : (*domains_)[var]) { for (IntegerValue v(interval.start); v <= interval.end; ++v) { - GetOrCreateLiteralAssociatedToEquality(var, v); + tmp_values_.push_back(v); } } + for (const IntegerValue v : tmp_values_) { + GetOrCreateLiteralAssociatedToEquality(var, v); + } // Mark var and Negation(var) as fully encoded. CHECK_LT(GetPositiveOnlyIndex(var), is_fully_encoded_.size()); diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index 1f75b94cf0..dada8e4841 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -473,6 +473,9 @@ class IntegerEncoder { // This will be lazily created when needed. LiteralIndex literal_index_true_ = kNoLiteralIndex; + // Temporary memory used by FullyEncodeVariable(). + std::vector tmp_values_; + DISALLOW_COPY_AND_ASSIGN(IntegerEncoder); }; diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 3a088318df..006145dc11 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -192,8 +192,9 @@ void LinearConstraintManager::AddCut( // Add the constraint. We only mark the constraint as a cut if it is not an // update of an already existing one. + const int64 prev_size = constraint_infos_.size(); const ConstraintIndex ct_index = Add(std::move(ct)); - if (ct_index + 1 == constraint_infos_.size()) { + if (prev_size + 1 == constraint_infos_.size()) { num_cuts_++; type_to_num_cuts_[type_name]++; constraint_infos_[ct_index].is_cut = true; diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc index cea70f3358..a74c4f3200 100644 --- a/ortools/sat/linear_relaxation.cc +++ b/ortools/sat/linear_relaxation.cc @@ -33,11 +33,7 @@ bool AppendFullEncodingRelaxation(IntegerVariable var, const Model& model, if (!encoder->VariableIsFullyEncoded(var)) return false; const auto& encoding = encoder->FullDomainEncoding(var); - // Compute min_value - IntegerValue var_min(kint64max); - for (const auto value_literal : encoding) { - var_min = std::min(var_min, value_literal.value); - } + const IntegerValue var_min = model.Get()->LowerBound(var); LinearConstraintBuilder at_least_one(&model, IntegerValue(1), kMaxIntegerValue); @@ -50,7 +46,7 @@ bool AppendFullEncodingRelaxation(IntegerVariable var, const Model& model, for (const auto value_literal : encoding) { const Literal lit = value_literal.literal; const IntegerValue delta = value_literal.value - var_min; - CHECK_GE(delta, IntegerValue(0)); + DCHECK_GE(delta, IntegerValue(0)); at_most_one.push_back(lit); if (!at_least_one.AddLiteralTerm(lit, IntegerValue(1))) return false; if (delta != IntegerValue(0)) { diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index 744ebb9112..4cbf296e67 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -351,70 +351,106 @@ void PresolveContext::InitializeNewDomains() { void PresolveContext::InsertVarValueEncoding(int literal, int ref, int64 value) { const int var = PositiveRef(ref); - const int64 s_value = RefIsPositive(ref) ? value : -value; - std::pair key{var, s_value}; - if (encoding.contains(key)) { - const int previous_literal = encoding[key]; + const int64 var_value = RefIsPositive(ref) ? value : -value; + const std::pair, int> key = + std::make_pair(std::make_pair(var, var_value), literal); + const auto& insert = encoding.insert(key); + if (insert.second) { + if (DomainOf(var).Size() == 2) { + // Encode the other literal. + const int64 var_min = MinOf(var); + const int64 var_max = MaxOf(var); + const int64 other_value = value == var_min ? var_max : var_min; + const std::pair other_key{var, other_value}; + auto other_it = encoding.find(other_key); + if (other_it != encoding.end()) { + // Other value in the domain was already encoded. + const int previous_other_literal = other_it->second; + if (previous_other_literal != NegatedRef(literal)) { + AddImplication(NegatedRef(literal), previous_other_literal); + AddImplication(previous_other_literal, NegatedRef(literal)); + } + } else { + encoding[other_key] = NegatedRef(literal); + // Add affine relation. + if (var_min != 0 || var_max != 1) { + ConstraintProto* const ct = working_model->add_constraints(); + LinearConstraintProto* const lin = ct->mutable_linear(); + lin->add_vars(var); + lin->add_coeffs(1); + lin->add_vars(literal); + lin->add_coeffs(var_min - var_max); + lin->add_domain(var_min); + lin->add_domain(var_min); + StoreAffineRelation(*ct, var, literal, var_max - var_min, var_min); + } + } + } else { + AddImplyInDomain(literal, var, Domain(var_value)); + AddImplyInDomain(NegatedRef(literal), var, + Domain(var_value).Complement()); + } + + } else { + const int previous_literal = insert.first->second; if (literal != previous_literal) { AddImplication(literal, previous_literal); AddImplication(previous_literal, literal); } - } else { - AddImplyInDomain(literal, var, Domain(s_value)); - AddImplyInDomain(NegatedRef(literal), var, Domain(s_value).Complement()); - encoding[key] = literal; } } int PresolveContext::GetOrCreateVarValueEncoding(int ref, int64 value) { // TODO(user,user): use affine relation here. const int var = PositiveRef(ref); - const int64 s_value = RefIsPositive(ref) ? value : -value; - if (!domains[var].Contains(s_value)) { + const int64 var_value = RefIsPositive(ref) ? value : -value; + + // Returns the false literal if the value is not in the domain. + if (!domains[var].Contains(var_value)) { return GetOrCreateConstantVar(0); } - std::pair key{var, s_value}; + // Returns the associated literal if already present. + const std::pair key{var, var_value}; auto it = encoding.find(key); - if (it != encoding.end()) return it->second; + if (it != encoding.end()) { + return it->second; + } + // Special case for fixed domains. if (domains[var].Size() == 1) { const int true_literal = GetOrCreateConstantVar(1); encoding[key] = true_literal; return true_literal; } + // Special case for domains of size 2. + const int64 var_min = MinOf(var); + const int64 var_max = MaxOf(var); if (domains[var].Size() == 2) { - const int64 var_min = MinOf(var); - const int64 var_max = MaxOf(var); - - if (var_min == 0 && var_max == 1) { - encoding[std::make_pair(var, 0)] = NegatedRef(var); - encoding[std::make_pair(var, 1)] = var; - } else { - const int literal = NewBoolVar(); - encoding[std::make_pair(var, var_min)] = NegatedRef(literal); - encoding[std::make_pair(var, var_max)] = literal; - - ConstraintProto* const ct = working_model->add_constraints(); - LinearConstraintProto* const lin = ct->mutable_linear(); - lin->add_vars(var); - lin->add_coeffs(1); - lin->add_vars(literal); - lin->add_coeffs(var_min - var_max); - lin->add_domain(var_min); - lin->add_domain(var_min); - StoreAffineRelation(*ct, var, literal, var_max - var_min, var_min); + // Checks if the other value is already encoded. + const int64 other_value = var_value == var_min ? var_max : var_min; + const std::pair other_key{var, other_value}; + auto other_it = encoding.find(other_key); + if (other_it != encoding.end()) { + // Fill in other value. + encoding[key] = NegatedRef(other_it->second); + return NegatedRef(other_it->second); } - return gtl::FindOrDieNoPrint(encoding, key); + if (var_min == 0 && var_max == 1) { + encoding[{var, 1}] = var; + encoding[{var, 0}] = NegatedRef(var); + return value == 1 ? var : NegatedRef(var); + } else { + const int literal = NewBoolVar(); + InsertVarValueEncoding(literal, var, var_max); + return var_value == var_max ? literal : NegatedRef(literal); + } } const int literal = NewBoolVar(); - AddImplyInDomain(literal, var, Domain(s_value)); - AddImplyInDomain(NegatedRef(literal), var, Domain(s_value).Complement()); - encoding[key] = literal; - + InsertVarValueEncoding(literal, var, var_value); return literal; } diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index 1e38efeda1..6989275ac9 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -121,7 +121,7 @@ struct PresolveContext { // Clears the "rules" statistics. void ClearStats(); - // Insert the given literal to encode ref == value. + // Inserts the given literal to encode ref == value. // If an encoding already exists, it adds the two implications between // the previous encoding and the new encoding. void InsertVarValueEncoding(int literal, int ref, int64 value); diff --git a/ortools/sat/presolve_util.cc b/ortools/sat/presolve_util.cc index 0d8a55bd85..db28ebcf78 100644 --- a/ortools/sat/presolve_util.cc +++ b/ortools/sat/presolve_util.cc @@ -14,6 +14,7 @@ #include "ortools/sat/presolve_util.h" #include "ortools/base/map_util.h" +#include "ortools/sat/cp_model_utils.h" namespace operations_research { namespace sat { @@ -95,5 +96,89 @@ std::vector> DomainDeductions::ProcessClause( return result; } +void SubstituteVariable(int var, int64 var_coeff_in_definition, + const ConstraintProto& definition, + ConstraintProto* ct) { + CHECK(RefIsPositive(var)); + CHECK_EQ(std::abs(var_coeff_in_definition), 1); + + // Copy all the terms (except the one refering to var). + std::vector> terms; + bool found = false; + int64 var_coeff = 0; + const int size = ct->linear().vars().size(); + for (int i = 0; i < size; ++i) { + int ref = ct->linear().vars(i); + int64 coeff = ct->linear().coeffs(i); + if (!RefIsPositive(ref)) { + ref = NegatedRef(ref); + coeff = -coeff; + } + + if (ref == var) { + CHECK(!found); + found = true; + var_coeff = coeff; + continue; + } else { + terms.push_back({ref, coeff}); + } + } + CHECK(found); + + if (var_coeff_in_definition < 0) var_coeff *= -1; + + // Add all the terms in the definition of var. + const int definition_size = definition.linear().vars().size(); + for (int i = 0; i < definition_size; ++i) { + int ref = definition.linear().vars(i); + int64 coeff = definition.linear().coeffs(i); + if (!RefIsPositive(ref)) { + ref = NegatedRef(ref); + coeff = -coeff; + } + + if (ref == var) { + CHECK_EQ(coeff, var_coeff_in_definition); + } else { + terms.push_back({ref, -coeff * var_coeff}); + } + } + + // The substitution is correct only if we don't loose information here. + // But for a constant definition rhs that is always the case. + bool exact = false; + Domain offset = ReadDomainFromProto(definition.linear()); + offset = offset.MultiplicationBy(-var_coeff, &exact); + CHECK(exact); + + const Domain rhs = ReadDomainFromProto(ct->linear()); + FillDomainInProto(rhs.AdditionWith(offset), ct->mutable_linear()); + + // Sort and merge terms refering to the same variable. + ct->mutable_linear()->clear_vars(); + ct->mutable_linear()->clear_coeffs(); + std::sort(terms.begin(), terms.end()); + int current_var = 0; + int64 current_coeff = 0; + for (const auto entry : terms) { + CHECK(RefIsPositive(entry.first)); + if (entry.first == current_var) { + current_coeff += entry.second; + } else { + if (current_coeff != 0) { + ct->mutable_linear()->add_vars(current_var); + ct->mutable_linear()->add_coeffs(current_coeff); + } + current_var = entry.first; + current_coeff = entry.second; + } + } + if (current_coeff != 0) { + ct->mutable_linear()->add_vars(current_var); + ct->mutable_linear()->add_coeffs(current_coeff); + } +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/presolve_util.h b/ortools/sat/presolve_util.h index ed1eed668a..70ce082cb2 100644 --- a/ortools/sat/presolve_util.h +++ b/ortools/sat/presolve_util.h @@ -22,6 +22,7 @@ #include "ortools/base/int_type_indexed_vector.h" #include "ortools/base/integral_types.h" #include "ortools/base/logging.h" +#include "ortools/sat/cp_model.pb.h" #include "ortools/util/bitset.h" #include "ortools/util/sorted_interval_list.h" @@ -81,6 +82,11 @@ class DomainDeductions { absl::flat_hash_map, Domain> deductions_; }; +// Replaces the variable var in ct using the definition constraint. +// Currently the coefficient in the definition must be 1 or -1. +void SubstituteVariable(int var, int64 var_coeff_in_definition, + const ConstraintProto& definition, ConstraintProto* ct); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index de441ceb0c..d896cc7deb 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -21,7 +21,7 @@ option java_multiple_files = true; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 147 +// NEXT TAG: 148 message SatParameters { // ========================================================================== // Branching and polarity @@ -420,6 +420,15 @@ message SatParameters { optional double merge_no_overlap_work_limit = 145 [default = 1e12]; optional double merge_at_most_one_work_limit = 146 [default = 1e8]; + // How much substitution (also called free variable aggregation in MIP + // litterature) should we perform at presolve. This currently only concerns + // variable appearing only in linear constraints. + // + // TODO(user): Implement a proper heuristic, maybe in term of non-zero + // created. Right now, we basically substitute all variables than can be + // easily substitued and that appear in at most this number of constraints. + optional int32 presolve_substitution_level = 147 [default = 0]; + // ========================================================================== // Max-sat parameters // ========================================================================== diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 898bfcd63a..b785022ac0 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -227,8 +227,10 @@ class SharedResponseManager { SharedSolutionRepository* MutableSolutionsRepository() { return &solutions_; } private: - void FillObjectiveValuesInBestResponse() EXCLUSIVE_LOCKS_REQUIRED(mutex_); - void SetStatsFromModelInternal(Model* model) EXCLUSIVE_LOCKS_REQUIRED(mutex_); + void FillObjectiveValuesInBestResponse() + EXCLUSIVE_LOCKS_REQUIRED(mutex_); + void SetStatsFromModelInternal(Model* model) + EXCLUSIVE_LOCKS_REQUIRED(mutex_); // Updates the primal integral using the old bounds on the objective. If the // old bounds are not finite, it uses the 'max_integral' value instead of gap.