From 7cd3ff1700a08cc21aa172f1842929920f2c6c0a Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 5 Oct 2021 17:11:30 +0200 Subject: [PATCH] [CP-SAT] better detection of encodings; use LinMax for AddAbsEquality --- ortools/sat/cp_model_postsolve.cc | 35 +++ ortools/sat/cp_model_presolve.cc | 424 ++++++++++++++++++++++-------- ortools/sat/cp_model_presolve.h | 6 +- ortools/sat/cp_model_utils.cc | 32 +++ ortools/sat/cp_model_utils.h | 19 ++ ortools/sat/presolve_context.cc | 7 + ortools/sat/presolve_context.h | 1 + ortools/sat/util.cc | 46 ++++ ortools/sat/util.h | 9 + 9 files changed, 475 insertions(+), 104 deletions(-) diff --git a/ortools/sat/cp_model_postsolve.cc b/ortools/sat/cp_model_postsolve.cc index 86c0f0f4d6..e0aad8b3a3 100644 --- a/ortools/sat/cp_model_postsolve.cc +++ b/ortools/sat/cp_model_postsolve.cc @@ -219,6 +219,38 @@ void PostsolveIntMax(const ConstraintProto& ct, std::vector* domains) { CHECK(!(*domains)[target_var].IsEmpty()); } +namespace { + +int64_t EvaluateLinearExpression(const LinearExpressionProto& expr, + const std::vector& domains) { + int64_t value = expr.offset(); + for (int i = 0; i < expr.vars_size(); ++i) { + const int ref = expr.vars(i); + const int64_t increment = + domains[PositiveRef(expr.vars(i))].FixedValue() * expr.coeffs(i); + value += RefIsPositive(ref) ? increment : -increment; + } + return value; +} + +} // namespace + +// Compute the max of each expression, and assign it to the target expr (which +// must be of the form +ref or -ref); +// We only support post-solving the case were the target is unassigned, +// but everything else is fixed. +void PostsolveLinMax(const ConstraintProto& ct, std::vector* domains) { + int64_t max_value = std::numeric_limits::min(); + for (const LinearExpressionProto& expr : ct.lin_max().exprs()) { + max_value = std::max(max_value, EvaluateLinearExpression(expr, *domains)); + } + const int target_ref = GetSingleRefFromExpression(ct.lin_max().target()); + const int target_var = PositiveRef(target_ref); + (*domains)[target_var] = (*domains)[target_var].IntersectionWith( + Domain(RefIsPositive(target_ref) ? max_value : -max_value)); + CHECK(!(*domains)[target_var].IsEmpty()); +} + // We only support 3 cases in the presolve currently. void PostsolveElement(const ConstraintProto& ct, std::vector* domains) { const int index_ref = ct.element().index(); @@ -358,6 +390,9 @@ void PostsolveResponse(const int64_t num_variables_in_original_model, case ConstraintProto::kIntMax: PostsolveIntMax(ct, &domains); break; + case ConstraintProto::kLinMax: + PostsolveLinMax(ct, &domains); + break; case ConstraintProto::kElement: PostsolveElement(ct, &domains); break; diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index cf5bc9d0bf..397d9f7a9b 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -761,8 +761,7 @@ bool CpModelPresolver::ConvertIntMax(ConstraintProto* ct) { return PresolveLinMax(ct); } -// TODO(user): Add all the missing presolve from PresolveIntMax(). -bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) { +bool CpModelPresolver::CanonicalizeLinMax(ConstraintProto* ct) { if (context_->ModelIsUnsat()) return false; // Canonicalize all involved expression. @@ -774,6 +773,12 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) { for (LinearExpressionProto& exp : *(ct->mutable_lin_max()->mutable_exprs())) { changed |= CanonicalizeLinearExpression(*ct, &exp); } + return changed; +} + +// TODO(user): Add all the missing presolve from PresolveIntMax(). +bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) { + if (context_->ModelIsUnsat()) return false; // Compute the infered min/max of the target. // Update target domain (if it is not a complex expression). @@ -812,6 +817,7 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) { // Filter the expressions which are smaller than target_min. const int64_t target_min = context_->MinOf(target); const int64_t target_max = context_->MaxOf(target); + bool changed = false; { int new_size = 0; for (int i = 0; i < ct->lin_max().exprs_size(); ++i) { @@ -1001,79 +1007,111 @@ bool CpModelPresolver::PresolveLinMax(ConstraintProto* ct) { return changed; } +// This presolve expect that the constraint only contains affine expressions. bool CpModelPresolver::PresolveIntAbs(ConstraintProto* ct) { CHECK_EQ(ct->enforcement_literal_size(), 0); if (context_->ModelIsUnsat()) return false; - const int target_ref = ct->int_max().target(); - const int var = PositiveRef(ct->int_max().vars(0)); + const LinearExpressionProto& target_expr = ct->lin_max().target(); + const LinearExpressionProto& expr = ct->lin_max().exprs(0); + DCHECK_EQ(expr.vars_size(), 1); - // Propagate from the variable domain to the target variable. - const Domain var_domain = context_->DomainOf(var); - const Domain new_target_domain = - var_domain.UnionWith(var_domain.Negation()) - .IntersectionWith({0, std::numeric_limits::max()}); - if (!context_->DomainOf(target_ref).IsIncludedIn(new_target_domain)) { - if (!context_->IntersectDomainWith(target_ref, new_target_domain)) { - return true; + // Propagate domain from the expression to the target. + { + const Domain expr_domain = context_->DomainSuperSetOf(expr); + const Domain new_target_domain = + expr_domain.UnionWith(expr_domain.Negation()) + .IntersectionWith({0, std::numeric_limits::max()}); + bool target_domain_modified = false; + if (!context_->IntersectDomainWith(target_expr, new_target_domain, + &target_domain_modified)) { + return false; + } + if (expr_domain.IsFixed()) { + context_->UpdateRuleStats("int_abs: fixed expression"); + return RemoveConstraint(ct); + } + if (target_domain_modified) { + context_->UpdateRuleStats("int_abs: propagate domain from x to abs(x)"); } - context_->UpdateRuleStats("int_abs: propagate domain x to abs(x)"); } // Propagate from target domain to variable. - const Domain target_domain = context_->DomainOf(target_ref); - const Domain new_var_domain = - target_domain.UnionWith(target_domain.Negation()); - if (!context_->DomainOf(var).IsIncludedIn(new_var_domain)) { - if (!context_->IntersectDomainWith(var, new_var_domain)) { + { + const Domain target_domain = + context_->DomainSuperSetOf(target_expr) + .IntersectionWith(Domain(0, std::numeric_limits::max())); + const Domain new_expr_domain = + target_domain.UnionWith(target_domain.Negation()); + bool expr_domain_modified = false; + if (!context_->IntersectDomainWith(expr, new_expr_domain, + &expr_domain_modified)) { return true; } - context_->UpdateRuleStats("int_abs: propagate domain abs(x) to x"); - } - - if (context_->MinOf(var) >= 0 && !context_->IsFixed(var)) { - context_->UpdateRuleStats("int_abs: converted to equality"); - ConstraintProto* new_ct = context_->working_model->add_constraints(); - new_ct->set_name(ct->name()); - auto* arg = new_ct->mutable_linear(); - arg->add_vars(target_ref); - arg->add_coeffs(1); - arg->add_vars(var); - arg->add_coeffs(-1); - arg->add_domain(0); - arg->add_domain(0); - context_->UpdateNewConstraintsVariableUsage(); - return RemoveConstraint(ct); - } - - if (context_->MaxOf(var) <= 0 && !context_->IsFixed(var)) { - context_->UpdateRuleStats("int_abs: converted to equality"); - ConstraintProto* new_ct = context_->working_model->add_constraints(); - new_ct->set_name(ct->name()); - auto* arg = new_ct->mutable_linear(); - arg->add_vars(target_ref); - arg->add_coeffs(1); - arg->add_vars(var); - arg->add_coeffs(1); - arg->add_domain(0); - arg->add_domain(0); - context_->UpdateNewConstraintsVariableUsage(); - return RemoveConstraint(ct); - } - - // Remove the abs constraint if the target is removable or fixed, as domains - // have been propagated. - if (context_->VariableIsUniqueAndRemovable(target_ref) || - context_->IsFixed(target_ref)) { - if (!context_->IsFixed(target_ref)) { - context_->MarkVariableAsRemoved(target_ref); - *context_->mapping_model->add_constraints() = *ct; + // This is the only reason why we don't support fully generic linear + // expression. + if (context_->IsFixed(target_expr)) { + context_->UpdateRuleStats("int_abs: fixed target"); + return RemoveConstraint(ct); } - context_->UpdateRuleStats("int_abs: remove constraint"); + if (expr_domain_modified) { + context_->UpdateRuleStats("int_abs: propagate domain from abs(x) to x"); + } + } + + // Convert to equality if the sign of expr is fixed. + if (context_->MinOf(expr) >= 0) { + context_->UpdateRuleStats("int_abs: converted to equality"); + ConstraintProto* new_ct = context_->working_model->add_constraints(); + new_ct->set_name(ct->name()); + auto* arg = new_ct->mutable_linear(); + arg->add_domain(0); + arg->add_domain(0); + AddLinearExpressionToLinearConstraint(target_expr, 1, arg); + AddLinearExpressionToLinearConstraint(expr, -1, arg); + if (!CanonicalizeLinear(new_ct)) return false; + context_->UpdateNewConstraintsVariableUsage(); return RemoveConstraint(ct); } - if (context_->StoreAbsRelation(target_ref, var)) { - context_->UpdateRuleStats("int_abs: store abs(x) == y"); + if (context_->MaxOf(expr) <= 0) { + context_->UpdateRuleStats("int_abs: converted to equality"); + ConstraintProto* new_ct = context_->working_model->add_constraints(); + new_ct->set_name(ct->name()); + auto* arg = new_ct->mutable_linear(); + arg->add_domain(0); + arg->add_domain(0); + AddLinearExpressionToLinearConstraint(target_expr, 1, arg); + AddLinearExpressionToLinearConstraint(expr, 1, arg); + if (!CanonicalizeLinear(new_ct)) return false; + context_->UpdateNewConstraintsVariableUsage(); + return RemoveConstraint(ct); + } + + // Remove the abs constraint if the target is removable and if domains have + // been propagated without loss. + // For now, we known that there is no loss if the target is a single ref. + // Since all the expression are affine, in this case we are fine. + if (ExpressionContainsSingleRef(target_expr) && + context_->VariableIsUniqueAndRemovable(target_expr.vars(0))) { + context_->MarkVariableAsRemoved(target_expr.vars(0)); + *context_->mapping_model->add_constraints() = *ct; + context_->UpdateRuleStats("int_abs: unused target"); + return RemoveConstraint(ct); + } + + // Store the x == abs(y) relation if expr and target_expr can be cast into a + // ref. + // TODO(user): Support general affine expression in for expr in the Store + // method call. + { + if (ExpressionContainsSingleRef(target_expr) && + ExpressionContainsSingleRef(expr)) { + const int target_ref = GetSingleRefFromExpression(target_expr); + const int expr_ref = GetSingleRefFromExpression(expr); + if (context_->StoreAbsRelation(target_ref, expr_ref)) { + context_->UpdateRuleStats("int_abs: store abs(x) == y"); + } + } } return false; @@ -1442,39 +1480,6 @@ void CpModelPresolver::DivideLinearByGcd(ConstraintProto* ct) { } } -void CpModelPresolver::PresolveLinearEqualityModuloTwo(ConstraintProto* ct) { - if (!ct->enforcement_literal().empty()) return; - if (ct->linear().domain().size() != 2) return; - if (ct->linear().domain(0) != ct->linear().domain(1)) return; - if (context_->ModelIsUnsat()) return; - - // Any equality must be true modulo n. - // The case modulo 2 is interesting if the non-zero terms are Booleans. - std::vector literals; - for (int i = 0; i < ct->linear().vars().size(); ++i) { - const int64_t coeff = ct->linear().coeffs(i); - const int ref = ct->linear().vars(i); - if (coeff % 2 == 0) continue; - if (!context_->CanBeUsedAsLiteral(ref)) return; - literals.push_back(PositiveRef(ref)); - if (literals.size() > 2) return; - } - if (literals.size() == 1) { - const int64_t rhs = std::abs(ct->linear().domain(0)); - context_->UpdateRuleStats("linear: only one odd Boolean in equality"); - if (!context_->IntersectDomainWith(literals[0], Domain(rhs % 2))) return; - } else if (literals.size() == 2) { - const int64_t rhs = std::abs(ct->linear().domain(0)); - context_->UpdateRuleStats("linear: only two odd Booleans in equality"); - if (rhs % 2) { - context_->StoreBooleanEqualityRelation(literals[0], - NegatedRef(literals[1])); - } else { - context_->StoreBooleanEqualityRelation(literals[0], literals[1]); - } - } -} - template bool CpModelPresolver::CanonicalizeLinearExpressionInternal( const ConstraintProto& ct, ProtoWithVarsAndCoeffs* proto, int64_t* offset) { @@ -1748,6 +1753,169 @@ bool CpModelPresolver::RemoveSingletonInLinear(ConstraintProto* ct) { return true; } +// If the gcd of all but one term (with index target_index) is not one, we can +// rewrite the last term using an affine representative. +bool CpModelPresolver::AddVarAffineRepresentativeFromLinearEquality( + int target_index, ConstraintProto* ct) { + int64_t gcd = 0; + const int num_variables = ct->linear().vars().size(); + for (int i = 0; i < num_variables; ++i) { + if (i == target_index) continue; + const int64_t magnitude = std::abs(ct->linear().coeffs(i)); + gcd = MathUtil::GCD64(gcd, magnitude); + if (gcd == 1) return false; + } + + // If we take the constraint % gcd, we have + // target_var * reduced_coeff = reduced_rhs + CHECK_GT(gcd, 1); + int target_var = ct->linear().vars(target_index); + int64_t reduced_coeff = ct->linear().coeffs(target_index); + if (!RefIsPositive(target_var)) { + target_var = NegatedRef(target_var); + reduced_coeff = -reduced_coeff; + } + reduced_coeff %= gcd; + if (reduced_coeff < 0) reduced_coeff += gcd; + int64_t reduced_rhs = ct->linear().domain(0) % gcd; + if (reduced_rhs < 0) reduced_rhs += gcd; + + // This should have been processed before by just dividing the whole + // constraint by the gcd. + if (reduced_coeff == 0) return false; + + // From target_var * reduced_coeff = reduced_rhs modulo gcd. + // We deduce that target_var = reduced_rhs * inverse + X * gcd; + CHECK_NE(reduced_coeff, 0); + const int64_t inverse = ModularInverse(reduced_coeff, gcd); + CHECK_NE(inverse, 0); + + // We abort if we have an overflow here. + const int64_t p = CapProd(inverse, reduced_rhs); + if (p == std::numeric_limits::max()) return false; + const int64_t offset = p % gcd; + + // We have target_var = gcd * X + offset ! + // Lets create a new integer variable and add the affine relation. + const Domain new_domain = context_->DomainOf(target_var) + .AdditionWith(Domain(-offset)) + .InverseMultiplicationBy(gcd); + if (new_domain.IsEmpty()) { + return context_->NotifyThatModelIsUnsat(); + } + if (new_domain.IsFixed()) { + if (!context_->IntersectDomainWith( + target_var, new_domain.ContinuousMultiplicationBy(gcd).AdditionWith( + Domain(offset)))) { + return false; + } + context_->UpdateRuleStats( + "linear: fixed variable due to affine representative."); + + // This will remove the now fixed variable and divide the rest by gcd. + return CanonicalizeLinear(ct); + } + + VLOG(2) << "In linear with " << num_variables << " terms, deduced " + << "(" << reduced_coeff << " * target ≡ " << reduced_rhs << " mod " + << gcd << ") => (target = " << gcd << " * X + " << offset + << "), target ∊ " << context_->DomainOf(target_var) << ", X ∊ " + << new_domain; + + // A potential problem with this and also the same code in + // TryToSimplifyDomain() is that it messes up the natural variable order + // chosen by the modeler. We try to correct that when mapping variables at the + // end of the presolve. + const int new_var = context_->NewIntVar(new_domain); + CHECK(context_->StoreAffineRelation(target_var, new_var, gcd, offset)); + context_->UpdateRuleStats("linear: canonicalize affine domain"); + context_->UpdateNewConstraintsVariableUsage(); + + // We use the new variable in the constraint. + // Note that we will divide everything by the gcd too. + return CanonicalizeLinear(ct); +} + +// Any equality must be true modulo n. +// +// If the gcd of all but one term is not one, we can rewrite the last term using +// an affine representative by considering the equality modulo that gcd. +// As an heuristic, we only test the smallest term or small primes 2, 3, and 5. +// +// We also handle the special case of having two non-zero literals modulo 2. +bool CpModelPresolver::PresolveLinearEqualityWithModulo(ConstraintProto* ct) { + if (context_->ModelIsUnsat()) return false; + if (ct->constraint_case() != ConstraintProto::kLinear) return false; + if (ct->linear().domain().size() != 2) return false; + if (ct->linear().domain(0) != ct->linear().domain(1)) return false; + if (!ct->enforcement_literal().empty()) return false; + + const int num_variables = ct->linear().vars().size(); + if (num_variables < 2) return false; + + std::vector mod2_indices; + std::vector mod3_indices; + std::vector mod5_indices; + + int64_t min_magnitude; + int num_smallest = 0; + int smallest_index; + for (int i = 0; i < num_variables; ++i) { + const int64_t magnitude = std::abs(ct->linear().coeffs(i)); + if (num_smallest == 0 || magnitude < min_magnitude) { + min_magnitude = magnitude; + num_smallest = 1; + smallest_index = i; + } else if (magnitude == min_magnitude) { + ++num_smallest; + } + + if (magnitude % 2 != 0) mod2_indices.push_back(i); + if (magnitude % 3 != 0) mod3_indices.push_back(i); + if (magnitude % 5 != 0) mod5_indices.push_back(i); + } + + if (mod2_indices.size() == 2) { + bool ok = true; + std::vector literals; + for (const int i : mod2_indices) { + const int ref = ct->linear().vars(i); + if (!context_->CanBeUsedAsLiteral(ref)) { + ok = false; + break; + } + literals.push_back(ref); + } + if (ok) { + const int64_t rhs = std::abs(ct->linear().domain(0)); + context_->UpdateRuleStats("linear: only two odd Booleans in equality"); + if (rhs % 2) { + context_->StoreBooleanEqualityRelation(literals[0], + NegatedRef(literals[1])); + } else { + context_->StoreBooleanEqualityRelation(literals[0], literals[1]); + } + } + } + + // TODO(user): More than one reduction might be possible, so we will need + // to call this again if we apply any of these reduction. + if (mod2_indices.size() == 1) { + return AddVarAffineRepresentativeFromLinearEquality(mod2_indices[0], ct); + } + if (mod3_indices.size() == 1) { + return AddVarAffineRepresentativeFromLinearEquality(mod3_indices[0], ct); + } + if (mod5_indices.size() == 1) { + return AddVarAffineRepresentativeFromLinearEquality(mod5_indices[0], ct); + } + if (num_smallest == 1) { + return AddVarAffineRepresentativeFromLinearEquality(smallest_index, ct); + } + + return false; +} + bool CpModelPresolver::PresolveSmallLinear(ConstraintProto* ct) { if (ct->constraint_case() != ConstraintProto::ConstraintCase::kLinear || context_->ModelIsUnsat()) { @@ -1892,6 +2060,9 @@ bool CpModelPresolver::PresolveSmallLinear(ConstraintProto* ct) { } else { // In this case, we can solve the diophantine equation, and write // both x and y in term of a new affine representative z. + // + // Note that PresolveLinearEqualityWithModularInverse() will have the + // same effect. context_->UpdateRuleStats("TODO linear: ax + by = cte"); } if (added) return RemoveConstraint(ct); @@ -5504,6 +5675,35 @@ void CpModelPresolver::TransformIntoMaxCliques() { } } +namespace { + +bool IsAffineIntAbs(ConstraintProto* ct) { + if (ct->constraint_case() != ConstraintProto::kLinMax || + ct->lin_max().exprs_size() != 2 || + ct->lin_max().target().vars_size() > 1 || + ct->lin_max().exprs(0).vars_size() != 1 || + ct->lin_max().exprs(1).vars_size() != 1) { + return false; + } + + const LinearArgumentProto& lin_max = ct->lin_max(); + if (lin_max.exprs(0).offset() != -lin_max.exprs(1).offset()) return false; + if (PositiveRef(lin_max.exprs(0).vars(0)) != + PositiveRef(lin_max.exprs(1).vars(0))) { + return false; + } + + const int64_t left_coeff = RefIsPositive(lin_max.exprs(0).vars(0)) + ? lin_max.exprs(0).coeffs(0) + : -lin_max.exprs(0).coeffs(0); + const int64_t right_coeff = RefIsPositive(lin_max.exprs(1).vars(0)) + ? lin_max.exprs(1).coeffs(0) + : -lin_max.exprs(1).coeffs(0); + return left_coeff == -right_coeff; +} + +} // namespace + bool CpModelPresolver::PresolveOneConstraint(int c) { if (context_->ModelIsUnsat()) return false; ConstraintProto* ct = context_->working_model->mutable_constraints(c); @@ -5531,16 +5731,18 @@ bool CpModelPresolver::PresolveOneConstraint(int c) { case ConstraintProto::ConstraintCase::kBoolXor: return PresolveBoolXor(ct); case ConstraintProto::ConstraintCase::kIntMax: - if (ct->int_max().vars_size() == 2 && - NegatedRef(ct->int_max().vars(0)) == ct->int_max().vars(1)) { - return PresolveIntAbs(ct); - } else { - return PresolveIntMax(ct); - } + return PresolveIntMax(ct); case ConstraintProto::ConstraintCase::kIntMin: return PresolveIntMin(ct); case ConstraintProto::ConstraintCase::kLinMax: - return PresolveLinMax(ct); + if (CanonicalizeLinMax(ct)) { + context_->UpdateConstraintVariableUsage(c); + } + if (IsAffineIntAbs(ct)) { + return PresolveIntAbs(ct); + } else { + return PresolveLinMax(ct); + } case ConstraintProto::ConstraintCase::kLinMin: return PresolveLinMin(ct); case ConstraintProto::ConstraintCase::kIntProd: @@ -5570,6 +5772,9 @@ bool CpModelPresolver::PresolveOneConstraint(int c) { if (PresolveSmallLinear(ct)) { context_->UpdateConstraintVariableUsage(c); } + if (PresolveLinearEqualityWithModulo(ct)) { + context_->UpdateConstraintVariableUsage(c); + } if (PropagateDomainsInLinear(c, ct)) { context_->UpdateConstraintVariableUsage(c); } @@ -5601,7 +5806,6 @@ bool CpModelPresolver::PresolveOneConstraint(int c) { PresolveSmallLinear(ct)) { context_->UpdateConstraintVariableUsage(c); } - PresolveLinearEqualityModuloTwo(ct); } // Note that it is important for this to be last, so that if a constraint @@ -7526,9 +7730,22 @@ bool CpModelPresolver::Presolve() { std::vector mapping(context_->working_model->variables_size(), -1); int num_free_variables = 0; for (int i = 0; i < context_->working_model->variables_size(); ++i) { + if (mapping[i] != -1) continue; // Already mapped. + if (context_->VariableIsNotUsedAnymore(i) && !context_->keep_all_feasible_solutions) { - if (!context_->VariableWasRemoved(i)) { + if (context_->VariableWasRemoved(i)) { + // Heuristic: If a variable is removed and has a representative that is + // not, we "move" the representative to the spot of that variable in the + // original order. This is to preserve any info encoded in the variable + // order by the modeler. + const int r = + PositiveRef(context_->GetAffineRelation(i).representative); + if (mapping[r] == -1 && !context_->VariableIsNotUsedAnymore(r)) { + mapping[r] = postsolve_mapping_->size(); + postsolve_mapping_->push_back(r); + } + } else { // Tricky. Variables that where not removed by a presolve rule should be // fixed first during postsolve, so that more complex postsolve rules // can use their values. One way to do that is to fix them here. @@ -7540,6 +7757,7 @@ bool CpModelPresolver::Presolve() { } continue; } + mapping[i] = postsolve_mapping_->size(); postsolve_mapping_->push_back(i); } diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h index 58a56086b4..6dc1d64638 100644 --- a/ortools/sat/cp_model_presolve.h +++ b/ortools/sat/cp_model_presolve.h @@ -136,6 +136,7 @@ class CpModelPresolver { int64_t* offset); bool CanonicalizeLinearExpression(const ConstraintProto& ct, LinearExpressionProto* exp); + bool CanonicalizeLinMax(ConstraintProto* ct); // For the linear constraints, we have more than one function. bool CanonicalizeLinear(ConstraintProto* ct); @@ -143,7 +144,10 @@ class CpModelPresolver { bool RemoveSingletonInLinear(ConstraintProto* ct); bool PresolveSmallLinear(ConstraintProto* ct); bool PresolveLinearOnBooleans(ConstraintProto* ct); - void PresolveLinearEqualityModuloTwo(ConstraintProto* ct); + bool AddVarAffineRepresentativeFromLinearEquality(int target_index, + ConstraintProto* ct); + bool PresolveLinearEqualityWithModulo(ConstraintProto* ct); + bool DetectAndProcessOneSidedLinearConstraint(int c, ConstraintProto* ct); // Scheduling helpers. diff --git a/ortools/sat/cp_model_utils.cc b/ortools/sat/cp_model_utils.cc index 6a0330ceb9..1eade1551f 100644 --- a/ortools/sat/cp_model_utils.cc +++ b/ortools/sat/cp_model_utils.cc @@ -575,5 +575,37 @@ int64_t ComputeInnerObjective(const CpObjectiveProto& objective, return objective_value; } +bool ExpressionContainsSingleRef(const LinearExpressionProto& expr) { + return expr.offset() == 0 && expr.vars_size() == 1 && + std::abs(expr.coeffs(0)) == 1; +} + +bool ExpressionIsAffine(const LinearExpressionProto& expr) { + return expr.vars_size() <= 1; +} + +// Returns the reference the expression can be reduced to. It will DCHECK that +// ExpressionContainsSingleRef(expr) is true. +int GetSingleRefFromExpression(const LinearExpressionProto& expr) { + DCHECK(ExpressionContainsSingleRef(expr)); + return expr.coeffs(0) == 1 ? expr.vars(0) : NegatedRef(expr.vars(0)); +} + +void AddLinearExpressionToLinearConstraint(const LinearExpressionProto& expr, + int64_t coefficient, + LinearConstraintProto* linear) { + for (int i = 0; i < expr.vars_size(); ++i) { + linear->add_vars(expr.vars(i)); + linear->add_coeffs(expr.coeffs(i) * coefficient); + } + DCHECK(!linear->domain().empty()); + const int64_t shift = coefficient * expr.offset(); + if (shift != 0) { + for (int64_t& d : *linear->mutable_domain()) { + d -= shift; + } + } +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_utils.h b/ortools/sat/cp_model_utils.h index c9f0dca1e7..3c191f6924 100644 --- a/ortools/sat/cp_model_utils.h +++ b/ortools/sat/cp_model_utils.h @@ -155,6 +155,25 @@ inline double UnscaleObjectiveValue(const CpObjectiveProto& proto, int64_t ComputeInnerObjective(const CpObjectiveProto& objective, const CpSolverResponse& response); +// Returns true if a linear expression can be reduced to a single ref. +bool ExpressionContainsSingleRef(const LinearExpressionProto& expr); + +// Checks if the expression is affine or constant. +bool ExpressionIsAffine(const LinearExpressionProto& expr); + +// Returns the reference the expression can be reduced to. It will DCHECK that +// ExpressionContainsSingleRef(expr) is true. +int GetSingleRefFromExpression(const LinearExpressionProto& expr); + +// Adds a linear expression proto to a linear constraint in place. +// +// Important: The domain must already be set, otherwise the offset will be lost. +// We also do not do any duplicate detection, so the constraint might need +// presolving afterwards. +void AddLinearExpressionToLinearConstraint(const LinearExpressionProto& expr, + int64_t coefficient, + LinearConstraintProto* linear); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index b965b72c48..9797a8dd67 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -147,6 +147,13 @@ int64_t PresolveContext::MaxOf(const LinearExpressionProto& expr) const { return result; } +bool PresolveContext::IsFixed(const LinearExpressionProto& expr) const { + for (int i = 0; i < expr.vars_size(); ++i) { + if (expr.coeffs(i) != 0 && !IsFixed(expr.vars(i))) return false; + } + return true; +} + Domain PresolveContext::DomainSuperSetOf( const LinearExpressionProto& expr) const { Domain result(expr.offset()); diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index 75ef957987..c5a94dce35 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -118,6 +118,7 @@ class PresolveContext { // should be such that this cannot happen (tested at validation). int64_t MinOf(const LinearExpressionProto& expr) const; int64_t MaxOf(const LinearExpressionProto& expr) const; + bool IsFixed(const LinearExpressionProto& expr) const; // Return a super-set of the domain of the linear expression. Domain DomainSuperSetOf(const LinearExpressionProto& expr) const; diff --git a/ortools/sat/util.cc b/ortools/sat/util.cc index cfa4140053..9a6db0a5ee 100644 --- a/ortools/sat/util.cc +++ b/ortools/sat/util.cc @@ -22,6 +22,52 @@ namespace operations_research { namespace sat { +namespace { +// This will be optimized into one division. I tested that in other places: +// +// Note that I am not 100% sure we need the indirection for the optimization +// to kick in though, but this seemed safer given our weird r[i ^ 1] inputs. +void QuotientAndRemainder(int64_t a, int64_t b, int64_t& q, int64_t& r) { + q = a / b; + r = a % b; +} +} // namespace + +// Using the extended Euclidian algo, we find a and b such that a x + b m = gcd. +// https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm +int64_t ModularInverse(int64_t x, int64_t m) { + DCHECK_GE(x, 0); + DCHECK_LT(x, m); + + int64_t r[2] = {m, x}; + int64_t t[2] = {0, 1}; + int64_t q; + + // We only keep the last two terms of the sequences with the "^1" trick: + // + // q = r[i-2] / r[i-1] + // r[i] = r[i-2] % r[i-1] + // t[i] = t[i-2] - t[i-1] * q + // + // We always have: + // - gcd(r[i], r[i - 1]) = gcd(r[i - 1], r[i - 2]) + // - x * t[i] + m * t[i - 1] = r[i] + int i = 0; + for (; r[i ^ 1] != 0; i ^= 1) { + QuotientAndRemainder(r[i], r[i ^ 1], q, r[i]); + t[i] -= t[i ^ 1] * q; + } + + // If the gcd is not one, there is no inverse, we returns 0. + if (r[i] != 1) return 0; + + // Correct the result so that it is in [0, m). Note that abs(t[i]) is known to + // be less than or equal to x / 2, and we have thorough unit-tests. + if (t[i] < 0) t[i] += m; + + return t[i]; +} + int MoveOneUnprocessedLiteralLast(const std::set& processed, int relevant_prefix_size, std::vector* literals) { diff --git a/ortools/sat/util.h b/ortools/sat/util.h index 404e11b5b3..a75692fb1b 100644 --- a/ortools/sat/util.h +++ b/ortools/sat/util.h @@ -32,6 +32,15 @@ namespace operations_research { namespace sat { +// Returns a in [0, m) such that a * x = 1 modulo m. +// If gcd(x, m) != 1, there is no inverse, and it returns 0. +// +// This DCHECK that x is in [0, m). +// This is integer overflow safe. +// +// Note(user): I didn't find this in a easily usable standard library. +int64_t ModularInverse(int64_t x, int64_t m); + // The model "singleton" random engine used in the solver. // // In test, we usually set use_absl_random() so that the sequence is changed at