From 9931205f7cb6e46f7db9d59ebd93053638edf550 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Fri, 9 Sep 2022 16:48:39 +0200 Subject: [PATCH] [CP-SAT] prepare for reservoir with variable demand; internal tweakes --- ortools/sat/cp_model.cc | 4 +- ortools/sat/cp_model.proto | 13 ++- ortools/sat/cp_model_checker.cc | 26 +++-- ortools/sat/cp_model_expand.cc | 158 ++++++++++++++++++++++++---- ortools/sat/cp_model_expand.h | 7 ++ ortools/sat/cp_model_lns.cc | 1 + ortools/sat/cp_model_loader.cc | 25 +++++ ortools/sat/cp_model_loader.h | 1 + ortools/sat/cp_model_presolve.cc | 40 +++++-- ortools/sat/cp_model_utils.cc | 7 ++ ortools/sat/csharp/Constraints.cs | 4 +- ortools/sat/cumulative.cc | 6 +- ortools/sat/disjunctive.cc | 8 +- ortools/sat/disjunctive.h | 1 - ortools/sat/integer_expr.cc | 72 ++++++------- ortools/sat/integer_expr.h | 6 -- ortools/sat/precedences.h | 80 -------------- ortools/sat/python/cp_model.py | 9 +- ortools/sat/sat_parameters.proto | 11 +- ortools/sat/scheduling_cuts.cc | 1 + ortools/sat/timetable.cc | 86 +++++++++------ ortools/sat/timetable.h | 6 +- ortools/sat/util.cc | 68 ++++++++++++ ortools/sat/util.h | 4 + ortools/util/parse_proto.cc | 91 ++++++++++++++++ ortools/util/parse_proto.h | 36 +++++++ ortools/util/saturated_arithmetic.h | 5 + 27 files changed, 558 insertions(+), 218 deletions(-) create mode 100644 ortools/util/parse_proto.cc create mode 100644 ortools/util/parse_proto.h diff --git a/ortools/sat/cp_model.cc b/ortools/sat/cp_model.cc index f41926645e..911974ac05 100644 --- a/ortools/sat/cp_model.cc +++ b/ortools/sat/cp_model.cc @@ -540,7 +540,7 @@ ReservoirConstraint::ReservoirConstraint(ConstraintProto* proto, void ReservoirConstraint::AddEvent(LinearExpr time, int64_t level_change) { *proto_->mutable_reservoir()->add_time_exprs() = builder_->LinearExprToProto(time); - proto_->mutable_reservoir()->add_level_changes(level_change); + proto_->mutable_reservoir()->add_level_changes()->set_offset(level_change); proto_->mutable_reservoir()->add_active_literals( builder_->IndexFromConstant(1)); } @@ -550,7 +550,7 @@ void ReservoirConstraint::AddOptionalEvent(LinearExpr time, BoolVar is_active) { *proto_->mutable_reservoir()->add_time_exprs() = builder_->LinearExprToProto(time); - proto_->mutable_reservoir()->add_level_changes(level_change); + proto_->mutable_reservoir()->add_level_changes()->set_offset(level_change); proto_->mutable_reservoir()->add_active_literals(is_active.index_); } diff --git a/ortools/sat/cp_model.proto b/ortools/sat/cp_model.proto index c87dac13c3..3e867addbb 100644 --- a/ortools/sat/cp_model.proto +++ b/ortools/sat/cp_model.proto @@ -163,11 +163,12 @@ message CumulativeConstraintProto { // Maintain a reservoir level within bounds. The water level starts at 0, and at // any time, it must be within [min_level, max_level]. // -// If the variable actives[i] is true, and if the expression times[i] is -// assigned a value t, then the current level changes by level_changes[i] (which -// is constant) at the time t. Therefore, at any time t: +// If the variable active_literals[i] is true, and if the expression +// time_exprs[i] is assigned a value t, then the current level changes by +// level_changes[i] at the time t. Therefore, at any time t: // -// sum(level_changes[i] * actives[i] if times[i] <= t) in [min_level, max_level] +// sum(level_changes[i] * active_literals[i] if time_exprs[i] <= t) +// in [min_level, max_level] // // Note that min level must be <= 0, and the max level must be >= 0. Please use // fixed level_changes to simulate initial state. @@ -179,8 +180,10 @@ message ReservoirConstraintProto { int64 min_level = 1; int64 max_level = 2; repeated LinearExpressionProto time_exprs = 3; // affine expressions. - repeated int64 level_changes = 4; // constants, can be negative. + // Currently, we only support constant level changes. + repeated LinearExpressionProto level_changes = 6; // affine expressions. repeated int32 active_literals = 5; + reserved 4; } // The circuit constraint is defined on a graph where the arc presence are diff --git a/ortools/sat/cp_model_checker.cc b/ortools/sat/cp_model_checker.cc index dc8c315fca..b308c84410 100644 --- a/ortools/sat/cp_model_checker.cc +++ b/ortools/sat/cp_model_checker.cc @@ -266,6 +266,15 @@ std::string ValidateAffineExpression(const CpModelProto& model, return ValidateLinearExpression(model, expr); } +std::string ValidateConstantAffineExpression( + const CpModelProto& model, const LinearExpressionProto& expr) { + if (!expr.vars().empty()) { + return absl::StrCat("expression must be constant: ", + ProtobufShortDebugString(expr)); + } + return ValidateLinearExpression(model, expr); +} + std::string ValidateLinearConstraint(const CpModelProto& model, const ConstraintProto& ct) { if (!DomainInProtoIsValid(ct.linear())) { @@ -666,6 +675,9 @@ std::string ValidateReservoirConstraint(const CpModelProto& model, for (const LinearExpressionProto& expr : ct.reservoir().time_exprs()) { RETURN_IF_NOT_EMPTY(ValidateAffineExpression(model, expr)); } + for (const LinearExpressionProto& expr : ct.reservoir().level_changes()) { + RETURN_IF_NOT_EMPTY(ValidateConstantAffineExpression(model, expr)); + } if (ct.reservoir().min_level() > 0) { return absl::StrCat( "The min level of a reservoir must be <= 0. Please use fixed events to " @@ -680,13 +692,11 @@ std::string ValidateReservoirConstraint(const CpModelProto& model, } int64_t sum_abs = 0; - for (const int64_t demand : ct.reservoir().level_changes()) { + for (const LinearExpressionProto& demand : ct.reservoir().level_changes()) { // We test for min int64_t before the abs(). - if (demand == std::numeric_limits::min()) { - return "Possible integer overflow in constraint: " + - ProtobufDebugString(ct); - } - sum_abs = CapAdd(sum_abs, std::abs(demand)); + const int64_t demand_min = MinOfExpression(model, demand); + const int64_t demand_max = MaxOfExpression(model, demand); + sum_abs = CapAdd(sum_abs, std::max(CapAbs(demand_min), CapAbs(demand_max))); if (sum_abs == std::numeric_limits::max()) { return "Possible integer overflow in constraint: " + ProtobufDebugString(ct); @@ -1496,7 +1506,9 @@ class ConstraintChecker { const int64_t time = LinearExpressionValue(ct.reservoir().time_exprs(i)); if (!has_active_variables || Value(ct.reservoir().active_literals(i)) == 1) { - deltas[time] += ct.reservoir().level_changes(i); + const int64_t level = + LinearExpressionValue(ct.reservoir().level_changes(i)); + deltas[time] += level; } } int64_t current_level = 0; diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index 302d49ac07..8628f54914 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -111,7 +111,8 @@ void ExpandReservoir(ConstraintProto* ct, PresolveContext* context) { int num_positives = 0; int num_negatives = 0; - for (const int64_t demand : reservoir.level_changes()) { + for (const LinearExpressionProto& demand_expr : reservoir.level_changes()) { + const int64_t demand = context->FixedValue(demand_expr); if (demand > 0) { num_positives++; } else if (demand < 0) { @@ -168,7 +169,7 @@ void ExpandReservoir(ConstraintProto* ct, PresolveContext* context) { const auto prec_it = precedence_cache.find({j, i}); CHECK(prec_it != precedence_cache.end()); const int prec_lit = prec_it->second; - const int64_t demand = reservoir.level_changes(j); + const int64_t demand = context->FixedValue(reservoir.level_changes(j)); if (RefIsPositive(prec_lit)) { level->mutable_linear()->add_vars(prec_lit); level->mutable_linear()->add_coeffs(demand); @@ -180,7 +181,7 @@ void ExpandReservoir(ConstraintProto* ct, PresolveContext* context) { } // Accounts for own demand in the domain of the sum. - const int64_t demand_i = reservoir.level_changes(i); + const int64_t demand_i = context->FixedValue(reservoir.level_changes(i)); level->mutable_linear()->add_domain( CapAdd(CapSub(reservoir.min_level(), demand_i), offset)); level->mutable_linear()->add_domain( @@ -193,7 +194,7 @@ void ExpandReservoir(ConstraintProto* ct, PresolveContext* context) { context->working_model->add_constraints()->mutable_linear(); for (int i = 0; i < num_events; ++i) { sum->add_vars(is_active_literal(i)); - sum->add_coeffs(reservoir.level_changes(i)); + sum->add_coeffs(context->FixedValue(reservoir.level_changes(i))); } sum->add_domain(reservoir.min_level()); sum->add_domain(reservoir.max_level()); @@ -1663,11 +1664,13 @@ void ExpandSomeLinearOfSizeTwo(ConstraintProto* ct, PresolveContext* context) { const int64_t coeff1 = arg.coeffs(0); const int64_t coeff2 = arg.coeffs(1); - const Domain reachable_rhs_superset = - context->DomainOf(var1).MultiplicationBy(coeff1).AdditionWith( - context->DomainOf(var2).MultiplicationBy(coeff2)); - + context->DomainOf(var1) + .MultiplicationBy(coeff1) + .RelaxIfTooComplex() + .AdditionWith(context->DomainOf(var2) + .MultiplicationBy(coeff2) + .RelaxIfTooComplex()); const Domain infeasible_reachable_values = reachable_rhs_superset.IntersectionWith( ReadDomainFromProto(arg).Complement()); @@ -1739,6 +1742,79 @@ void ExpandSomeLinearOfSizeTwo(ConstraintProto* ct, PresolveContext* context) { ct->Clear(); } +// Note that we used to do that at loading time, but we prefer to do that as +// part of the presolve so that all variables are available for sharing between +// subworkers and also are accessible by the linear relaxation. +// +// TODO(user): Note that currently both encoding introduce extra solutions +// if the constraint has some enforcement literal(). We can either fix this by +// supporting enumeration on a subset of variable. Or add extra constraint to +// fix all new Boolean to false if the initial constraint is not enforced. +void ExpandComplexLinearConstraint(int c, ConstraintProto* ct, + PresolveContext* context) { + // TODO(user): We treat the linear of size 1 differently because we need them + // as is to recognize value encoding. Try to still creates needed Boolean now + // so that we can share more between the different workers. Or revisit how + // linear1 are propagated. + if (ct->linear().domain().size() <= 2) return; + if (ct->linear().vars().size() == 1) return; + + const SatParameters& params = context->params(); + if (params.encode_complex_linear_constraint_with_integer()) { + // Integer encoding. + // + // Here we add a slack with domain equal to rhs and transform + // expr \in rhs to expr - slack = 0 + const Domain rhs = ReadDomainFromProto(ct->linear()); + const int slack = context->NewIntVar(rhs); + ct->mutable_linear()->add_vars(slack); + ct->mutable_linear()->add_coeffs(-1); + ct->mutable_linear()->clear_domain(); + ct->mutable_linear()->add_domain(0); + ct->mutable_linear()->add_domain(0); + } else { + // Boolean encoding. + int single_bool; + BoolArgumentProto* clause = nullptr; + if (ct->enforcement_literal().empty() && ct->linear().domain_size() == 4) { + // We cover the special case of no enforcement and two choices by creating + // a single Boolean. + single_bool = context->NewBoolVar(); + } else { + clause = context->working_model->add_constraints()->mutable_bool_or(); + for (const int ref : ct->enforcement_literal()) { + clause->add_literals(NegatedRef(ref)); + } + } + ct->mutable_enforcement_literal()->Clear(); + for (int i = 0; i < ct->linear().domain_size(); i += 2) { + const int64_t lb = ct->linear().domain(i); + const int64_t ub = ct->linear().domain(i + 1); + + int subdomain_literal; + if (clause != nullptr) { + subdomain_literal = context->NewBoolVar(); + clause->add_literals(subdomain_literal); + } else { + subdomain_literal = i == 0 ? single_bool : NegatedRef(single_bool); + } + + // Create a new constraint which is a copy of the original, but with a + // simple sub-domain and enforcement literal. + ConstraintProto* new_ct = context->working_model->add_constraints(); + *new_ct = *ct; + new_ct->add_enforcement_literal(subdomain_literal); + FillDomainInProto(Domain(lb, ub), new_ct->mutable_linear()); + } + ct->Clear(); + } + + context->UpdateRuleStats("linear: expanded complex rhs"); + context->InitializeNewDomains(); + context->UpdateNewConstraintsVariableUsage(); + context->UpdateConstraintVariableUsage(c); +} + } // namespace void ExpandCpModel(PresolveContext* context) { @@ -1756,29 +1832,52 @@ void ExpandCpModel(PresolveContext* context) { context->ClearPrecedenceCache(); // First pass: we look at constraints that may fully encode variables. - for (int i = 0; i < context->working_model->constraints_size(); ++i) { - ConstraintProto* const ct = context->working_model->mutable_constraints(i); + for (int c = 0; c < context->working_model->constraints_size(); ++c) { + ConstraintProto* const ct = context->working_model->mutable_constraints(c); bool skip = false; switch (ct->constraint_case()) { - case ConstraintProto::ConstraintCase::kReservoir: - ExpandReservoir(ct, context); + case ConstraintProto::kLinear: + // If we only do expansion, we do that as part of the main loop. + // This way we don't need to call FinalExpansionForLinearConstraint(). + if (ct->linear().domain().size() > 2 && + !context->params().cp_model_presolve()) { + ExpandComplexLinearConstraint(c, ct, context); + } break; - case ConstraintProto::ConstraintCase::kIntMod: + case ConstraintProto::kReservoir: + if (context->params().expand_reservoir_constraints()) { + for (const LinearExpressionProto& demand_expr : + ct->reservoir().level_changes()) { + if (!context->IsFixed(demand_expr)) { + skip = true; + break; + } + } + if (skip) { + context->UpdateRuleStats( + "reservoir: expansion is not supported with variable level " + "changes"); + } else { + ExpandReservoir(ct, context); + } + } + break; + case ConstraintProto::kIntMod: ExpandIntMod(ct, context); break; - case ConstraintProto::ConstraintCase::kIntProd: + case ConstraintProto::kIntProd: ExpandIntProd(ct, context); break; - case ConstraintProto::ConstraintCase::kElement: + case ConstraintProto::kElement: ExpandElement(ct, context); break; - case ConstraintProto::ConstraintCase::kInverse: + case ConstraintProto::kInverse: ExpandInverse(ct, context); break; - case ConstraintProto::ConstraintCase::kAutomaton: + case ConstraintProto::kAutomaton: ExpandAutomaton(ct, context); break; - case ConstraintProto::ConstraintCase::kTable: + case ConstraintProto::kTable: if (ct->table().negated()) { ExpandNegativeTable(ct, context); } else { @@ -1794,7 +1893,7 @@ void ExpandCpModel(PresolveContext* context) { // Update variable-constraint graph. context->UpdateNewConstraintsVariableUsage(); if (ct->constraint_case() == ConstraintProto::CONSTRAINT_NOT_SET) { - context->UpdateConstraintVariableUsage(i); + context->UpdateConstraintVariableUsage(c); } // Early exit if the model is unsat. @@ -1811,11 +1910,11 @@ void ExpandCpModel(PresolveContext* context) { ConstraintProto* const ct = context->working_model->mutable_constraints(i); bool skip = false; switch (ct->constraint_case()) { - case ConstraintProto::ConstraintCase::kAllDiff: + case ConstraintProto::kAllDiff: ExpandAllDiff(context->params().expand_alldiff_constraints(), ct, context); break; - case ConstraintProto::ConstraintCase::kLinear: + case ConstraintProto::kLinear: ExpandSomeLinearOfSizeTwo(ct, context); break; default: @@ -1856,5 +1955,22 @@ void ExpandCpModel(PresolveContext* context) { context->NotifyThatModelIsExpanded(); } +void FinalExpansionForLinearConstraint(PresolveContext* context) { + if (context->params().disable_constraint_expansion()) return; + if (context->ModelIsUnsat()) return; + for (int c = 0; c < context->working_model->constraints_size(); ++c) { + ConstraintProto* const ct = context->working_model->mutable_constraints(c); + switch (ct->constraint_case()) { + case ConstraintProto::kLinear: + if (ct->linear().domain().size() > 2) { + ExpandComplexLinearConstraint(c, ct, context); + } + break; + default: + break; + } + } +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_expand.h b/ortools/sat/cp_model_expand.h index 93990bda8e..f4c5d6e8ea 100644 --- a/ortools/sat/cp_model_expand.h +++ b/ortools/sat/cp_model_expand.h @@ -28,6 +28,13 @@ namespace sat { // simplification of the model. Furthermore, this expansion is mandatory. void ExpandCpModel(PresolveContext* context); +// Linear constraint with a complex rhs need to be expanded at the end of the +// presolve. We do that at the end, because the presolve is allowed to simplify +// such constraints by updating the rhs. Also the extra variable we create are +// only linked by a few constraints to the rest of the model and should not be +// presolvable. +void FinalExpansionForLinearConstraint(PresolveContext* context); + // Fills and propagates the set of reachable states/labels. void PropagateAutomaton(const AutomatonConstraintProto& proto, const PresolveContext& context, diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 2bde0e2303..80c6e749d7 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -36,6 +36,7 @@ #include "absl/time/clock.h" #include "absl/time/time.h" #include "absl/types/span.h" +#include "ortools/base/check.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" #include "ortools/graph/connected_components.h" diff --git a/ortools/sat/cp_model_loader.cc b/ortools/sat/cp_model_loader.cc index 4cc26e5f20..dfae0e35ef 100644 --- a/ortools/sat/cp_model_loader.cc +++ b/ortools/sat/cp_model_loader.cc @@ -51,6 +51,7 @@ #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/sat_solver.h" #include "ortools/sat/symmetry.h" +#include "ortools/sat/timetable.h" #include "ortools/util/logging.h" #include "ortools/util/sorted_interval_list.h" #include "ortools/util/strong_integers.h" @@ -1230,6 +1231,27 @@ void LoadCumulativeConstraint(const ConstraintProto& ct, Model* m) { m->Add(Cumulative(intervals, demands, capacity)); } +void LoadReservoirConstraint(const ConstraintProto& ct, Model* m) { + auto* mapping = m->GetOrCreate(); + auto* encoder = m->GetOrCreate(); + const std::vector times = + mapping->Affines(ct.reservoir().time_exprs()); + const std::vector level_changes = + mapping->Affines(ct.reservoir().level_changes()); + std::vector presences; + const int size = ct.reservoir().time_exprs().size(); + for (int i = 0; i < size; ++i) { + if (!ct.reservoir().active_literals().empty()) { + presences.push_back(mapping->Literal(ct.reservoir().active_literals(i))); + } else { + presences.push_back(encoder->GetTrueLiteral()); + } + } + AddReservoirConstraint(times, level_changes, presences, + ct.reservoir().min_level(), ct.reservoir().max_level(), + m); +} + void LoadCircuitConstraint(const ConstraintProto& ct, Model* m) { const auto& circuit = ct.circuit(); if (circuit.tails().empty()) return; @@ -1304,6 +1326,9 @@ bool LoadConstraint(const ConstraintProto& ct, Model* m) { case ConstraintProto::ConstraintProto::kCumulative: LoadCumulativeConstraint(ct, m); return true; + case ConstraintProto::ConstraintProto::kReservoir: + LoadReservoirConstraint(ct, m); + return true; case ConstraintProto::ConstraintProto::kCircuit: LoadCircuitConstraint(ct, m); return true; diff --git a/ortools/sat/cp_model_loader.h b/ortools/sat/cp_model_loader.h index 506d73562b..be97562b09 100644 --- a/ortools/sat/cp_model_loader.h +++ b/ortools/sat/cp_model_loader.h @@ -105,6 +105,7 @@ void LoadNoOverlapConstraint(const ConstraintProto& ct, Model* m); void LoadNoOverlap2dConstraint(const ConstraintProto& ct, Model* m); void LoadCumulativeConstraint(const ConstraintProto& ct, Model* m); void LoadCircuitConstraint(const ConstraintProto& ct, Model* m); +void LoadReservoirConstraint(const ConstraintProto& ct, Model* m); void LoadRoutesConstraint(const ConstraintProto& ct, Model* m); void LoadCircuitCoveringConstraint(const ConstraintProto& ct, Model* m); diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index aab670d76b..78844ba470 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -5472,6 +5472,9 @@ bool CpModelPresolver::PresolveReservoir(ConstraintProto* ct) { for (LinearExpressionProto& exp : *(proto.mutable_time_exprs())) { changed |= CanonicalizeLinearExpression(*ct, &exp); } + for (LinearExpressionProto& exp : *(proto.mutable_level_changes())) { + changed |= CanonicalizeLinearExpression(*ct, &exp); + } if (proto.active_literals().empty()) { const int true_literal = context_->GetTrueLiteral(); @@ -5482,7 +5485,8 @@ bool CpModelPresolver::PresolveReservoir(ConstraintProto* ct) { } const auto& demand_is_null = [&](int i) { - return proto.level_changes(i) == 0 || + return (context_->IsFixed(proto.level_changes(i)) && + context_->FixedValue(proto.level_changes(i)) == 0) || context_->LiteralIsFalse(proto.active_literals(i)); }; @@ -5497,13 +5501,15 @@ bool CpModelPresolver::PresolveReservoir(ConstraintProto* ct) { int new_size = 0; for (int i = 0; i < proto.level_changes_size(); ++i) { if (demand_is_null(i)) continue; - proto.set_level_changes(new_size, proto.level_changes(i)); + *proto.mutable_level_changes(new_size) = proto.level_changes(i); *proto.mutable_time_exprs(new_size) = proto.time_exprs(i); proto.set_active_literals(new_size, proto.active_literals(i)); new_size++; } - proto.mutable_level_changes()->Truncate(new_size); + proto.mutable_level_changes()->erase( + proto.mutable_level_changes()->begin() + new_size, + proto.mutable_level_changes()->end()); proto.mutable_time_exprs()->erase( proto.mutable_time_exprs()->begin() + new_size, proto.mutable_time_exprs()->end()); @@ -5513,15 +5519,21 @@ bool CpModelPresolver::PresolveReservoir(ConstraintProto* ct) { "reservoir: remove zero level_changes or inactive events."); } + // The rest of the presolve only applies if all demands are fixed. + for (const LinearExpressionProto& level_change : proto.level_changes()) { + if (!context_->IsFixed(level_change)) return changed; + } + const int num_events = proto.level_changes_size(); - int64_t gcd = - proto.level_changes().empty() ? 0 : std::abs(proto.level_changes(0)); + int64_t gcd = proto.level_changes().empty() + ? 0 + : std::abs(context_->FixedValue(proto.level_changes(0))); int num_positives = 0; int num_negatives = 0; int64_t max_sum_of_positive_level_changes = 0; int64_t min_sum_of_negative_level_changes = 0; for (int i = 0; i < num_events; ++i) { - const int64_t demand = proto.level_changes(i); + const int64_t demand = context_->FixedValue(proto.level_changes(i)); gcd = MathUtil::GCD64(gcd, std::abs(demand)); if (demand > 0) { num_positives++; @@ -5564,7 +5576,7 @@ bool CpModelPresolver::PresolveReservoir(ConstraintProto* ct) { context_->working_model->add_constraints()->mutable_linear(); int64_t fixed_contrib = 0; for (int i = 0; i < proto.level_changes_size(); ++i) { - const int64_t demand = proto.level_changes(i); + const int64_t demand = context_->FixedValue(proto.level_changes(i)); DCHECK_NE(demand, 0); const int active = proto.active_literals(i); @@ -5585,7 +5597,10 @@ bool CpModelPresolver::PresolveReservoir(ConstraintProto* ct) { if (gcd > 1) { for (int i = 0; i < proto.level_changes_size(); ++i) { - proto.set_level_changes(i, proto.level_changes(i) / gcd); + proto.mutable_level_changes(i)->set_offset( + context_->FixedValue(proto.level_changes(i)) / gcd); + proto.mutable_level_changes(i)->clear_vars(); + proto.mutable_level_changes(i)->clear_coeffs(); } // Adjust min and max levels. @@ -9392,6 +9407,12 @@ CpSolverStatus CpModelPresolver::Presolve() { ExpandCpModel(context_); if (context_->ModelIsUnsat()) return InfeasibleStatus(); + // We still write back the canonical objective has we don't deal well + // with uninitialized domain or duplicate variables. + if (context_->working_model->has_objective()) { + context_->WriteObjectiveToProto(); + } + // We need to append all the variable equivalence that are still used! EncodeAllAffineRelations(); if (logger_->LoggingIsEnabled()) context_->LogInfo(); @@ -9549,6 +9570,9 @@ CpSolverStatus CpModelPresolver::Presolve() { context_->WriteObjectiveToProto(); } + // Take care of linear constraint with a complex rhs. + FinalExpansionForLinearConstraint(context_); + // Adds all needed affine relation to context_->working_model. EncodeAllAffineRelations(); if (context_->ModelIsUnsat()) return InfeasibleStatus(); diff --git a/ortools/sat/cp_model_utils.cc b/ortools/sat/cp_model_utils.cc index b2da364dfe..94bceeed4c 100644 --- a/ortools/sat/cp_model_utils.cc +++ b/ortools/sat/cp_model_utils.cc @@ -120,6 +120,10 @@ IndexReferences GetReferencesUsedByConstraint(const ConstraintProto& ct) { for (const LinearExpressionProto& time : ct.reservoir().time_exprs()) { AddIndices(time.vars(), &output.variables); } + for (const LinearExpressionProto& level : + ct.reservoir().level_changes()) { + AddIndices(level.vars(), &output.variables); + } AddIndices(ct.reservoir().active_literals(), &output.literals); break; case ConstraintProto::ConstraintCase::kTable: @@ -289,6 +293,9 @@ void ApplyToAllVariableIndices(const std::function& f, for (int i = 0; i < ct->reservoir().time_exprs_size(); ++i) { APPLY_TO_REPEATED_FIELD(reservoir, time_exprs(i)->mutable_vars); } + for (int i = 0; i < ct->reservoir().level_changes_size(); ++i) { + APPLY_TO_REPEATED_FIELD(reservoir, level_changes(i)->mutable_vars); + } break; case ConstraintProto::ConstraintCase::kTable: APPLY_TO_REPEATED_FIELD(table, vars); diff --git a/ortools/sat/csharp/Constraints.cs b/ortools/sat/csharp/Constraints.cs index 1d8f32f285..33bc636cc0 100644 --- a/ortools/sat/csharp/Constraints.cs +++ b/ortools/sat/csharp/Constraints.cs @@ -324,7 +324,7 @@ public class ReservoirConstraint : Constraint { ReservoirConstraintProto res = Proto.Reservoir; res.TimeExprs.Add(cp_model_.GetLinearExpressionProto(cp_model_.GetLinearExpr(time))); - res.LevelChanges.Add(Convert.ToInt64(level_change)); + res.LevelChanges.Add(cp_model_.GetLinearExpressionProto(cp_model_.GetLinearExpr(level_change))); res.ActiveLiterals.Add(cp_model_.TrueLiteral().GetIndex()); return this; } @@ -343,7 +343,7 @@ public class ReservoirConstraint : Constraint { ReservoirConstraintProto res = Proto.Reservoir; res.TimeExprs.Add(cp_model_.GetLinearExpressionProto(cp_model_.GetLinearExpr(time))); - res.LevelChanges.Add(Convert.ToInt64(level_change)); + res.LevelChanges.Add(cp_model_.GetLinearExpressionProto(cp_model_.GetLinearExpr(level_change))); res.ActiveLiterals.Add(literal.GetIndex()); return this; } diff --git a/ortools/sat/cumulative.cc b/ortools/sat/cumulative.cc index 936a2562f9..c01fb4b76b 100644 --- a/ortools/sat/cumulative.cc +++ b/ortools/sat/cumulative.cc @@ -359,16 +359,16 @@ std::function CumulativeUsingReservoir( integer_trail->UpperBound(capacity).value()); std::vector times; - std::vector deltas; + std::vector deltas; std::vector presences; const int num_tasks = vars.size(); for (int t = 0; t < num_tasks; ++t) { CHECK(integer_trail->IsFixed(demands[t])); times.push_back(intervals->StartVar(vars[t])); - deltas.push_back(integer_trail->LowerBound(demands[t])); + deltas.push_back(demands[t]); times.push_back(intervals->EndVar(vars[t])); - deltas.push_back(-integer_trail->LowerBound(demands[t])); + deltas.push_back(demands[t].Negated()); if (intervals->IsOptional(vars[t])) { presences.push_back(intervals->PresenceLiteral(vars[t])); presences.push_back(intervals->PresenceLiteral(vars[t])); diff --git a/ortools/sat/disjunctive.cc b/ortools/sat/disjunctive.cc index 1f09a5c88e..ea1a263c1e 100644 --- a/ortools/sat/disjunctive.cc +++ b/ortools/sat/disjunctive.cc @@ -205,11 +205,6 @@ void TaskSet::NotifyEntryIsNowLastIfPresent(const Entry& e) { DCHECK(std::is_sorted(sorted_tasks_.begin(), sorted_tasks_.end())); } -void TaskSet::RemoveEntryWithIndex(int index) { - sorted_tasks_.erase(sorted_tasks_.begin() + index); - optimized_restart_ = 0; -} - IntegerValue TaskSet::ComputeEndMin() const { DCHECK(std::is_sorted(sorted_tasks_.begin(), sorted_tasks_.end())); const int size = sorted_tasks_.size(); @@ -456,6 +451,9 @@ bool CombinedDisjunctive::Propagate() { const IntegerValue shifted_smin = helper_->ShiftedStartMin(t); const IntegerValue size_min = helper_->SizeMin(t); for (const int d_index : task_to_disjunctives_[t]) { + // TODO(user): Refactor the code to use the same algo as in + // DisjunctiveDetectablePrecedences, it is superior and do not need + // this function. task_sets_[d_index].NotifyEntryIsNowLastIfPresent( {t, shifted_smin, size_min}); end_mins_[d_index] = task_sets_[d_index].ComputeEndMin(); diff --git a/ortools/sat/disjunctive.h b/ortools/sat/disjunctive.h index 580a1f944a..05ef9313dc 100644 --- a/ortools/sat/disjunctive.h +++ b/ortools/sat/disjunctive.h @@ -72,7 +72,6 @@ class TaskSet { optimized_restart_ = 0; } void AddEntry(const Entry& e); - void RemoveEntryWithIndex(int index); // Same as AddEntry({t, helper->ShiftedStartMin(t), helper->SizeMin(t)}). // This is a minor optimization to not call SizeMin(t) twice. diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index 4c69419fcb..6b63376e0d 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -98,15 +98,16 @@ std::pair IntegerSumLE::ConditionalLb( bool target_var_present_negatively = false; IntegerValue target_coeff; - // Compute the implied_lb excluding "- target_coeff * target". - IntegerValue implied_lb(-upper_bound_); + // Warning: It is important to do the computation like the propagation is + // doing it to be sure we don't have overflow, since this is what we check + // when creating constraints. + IntegerValue implied_lb(0); for (int i = 0; i < vars_.size(); ++i) { const IntegerVariable var = vars_[i]; const IntegerValue coeff = coeffs_[i]; if (var == NegationOf(target_var)) { target_coeff = coeff; target_var_present_negatively = true; - continue; } const IntegerValue lb = integer_trail_->LowerBound(var); @@ -117,10 +118,25 @@ std::pair IntegerSumLE::ConditionalLb( literal_var_present_positively = (var == integer_literal.var); } } + if (!literal_var_present || !target_var_present_negatively) { return {kMinIntegerValue, kMinIntegerValue}; } + // The upper bound on NegationOf(target_var) are lb(-target) + slack / coeff. + // So the lower bound on target_var is ub - slack / coeff. + const IntegerValue slack = upper_bound_ - implied_lb; + const IntegerValue target_lb = integer_trail_->LowerBound(target_var); + const IntegerValue target_ub = integer_trail_->UpperBound(target_var); + if (slack <= 0) { + // TODO(user): If there is a conflict (negative slack) we can be more + // precise. + return {target_ub, target_ub}; + } + + const IntegerValue target_diff = target_ub - target_lb; + const IntegerValue delta = std::min(slack / target_coeff, target_diff); + // A literal means var >= bound. if (literal_var_present_positively) { // We have var_coeff * var in the expression, the literal is var >= bound. @@ -129,8 +145,11 @@ std::pair IntegerSumLE::ConditionalLb( const IntegerValue diff = std::max( IntegerValue(0), integer_literal.bound - integer_trail_->LowerBound(integer_literal.var)); - return {CeilRatio(implied_lb, target_coeff), - CeilRatio(implied_lb + var_coeff * diff, target_coeff)}; + const IntegerValue tighter_slack = + std::max(IntegerValue(0), slack - var_coeff * diff); + const IntegerValue tighter_delta = + std::min(tighter_slack / target_coeff, target_diff); + return {target_ub - delta, target_ub - tighter_delta}; } else { // We have var_coeff * -var in the expression, the literal is var >= bound. // When it is true, it is not relevant as implied_lb used -var >= -ub. @@ -138,45 +157,14 @@ std::pair IntegerSumLE::ConditionalLb( const IntegerValue diff = std::max( IntegerValue(0), integer_trail_->UpperBound(integer_literal.var) - integer_literal.bound + 1); - return {CeilRatio(implied_lb + var_coeff * diff, target_coeff), - CeilRatio(implied_lb, target_coeff)}; + const IntegerValue tighter_slack = + std::max(IntegerValue(0), slack - var_coeff * diff); + const IntegerValue tighter_delta = + std::min(tighter_slack / target_coeff, target_diff); + return {target_ub - tighter_delta, target_ub - delta}; } } -std::vector> -IntegerSumLE::ConditionalLbs(IntegerVariable target_var) const { - bool target_var_present_negatively = false; - IntegerValue target_coeff; - - // Compute the implied_lb excluding "- target_coeff * target". - IntegerValue implied_lb(-upper_bound_); - for (int i = 0; i < vars_.size(); ++i) { - const IntegerVariable var = vars_[i]; - const IntegerValue coeff = coeffs_[i]; - if (var == NegationOf(target_var)) { - target_coeff = coeff; - target_var_present_negatively = true; - continue; - } - implied_lb += coeff * integer_trail_->LowerBound(var); - } - - std::vector> result; - if (!target_var_present_negatively) return result; - - for (int i = 0; i < vars_.size(); ++i) { - const IntegerVariable var = vars_[i]; - const IntegerValue coeff = coeffs_[i]; - if (var == NegationOf(target_var)) continue; - - const IntegerValue lb = integer_trail_->LowerBound(var); - if (lb == integer_trail_->UpperBound(var)) continue; - result.push_back({IntegerLiteral::GreaterOrEqual(var, lb + 1), - CeilRatio(implied_lb + coeff, target_coeff)}); - } - return result; -} - bool IntegerSumLE::Propagate() { // Reified case: If any of the enforcement_literals are false, we ignore the // constraint. @@ -252,6 +240,8 @@ bool IntegerSumLE::Propagate() { for (int i = rev_num_fixed_vars_; i < num_vars; ++i) { if (max_variations_[i] <= slack) continue; + // TODO(user): If the new ub fall into an hole of the variable, we can + // actually relax the reason more by computing a better slack. const IntegerVariable var = vars_[i]; const IntegerValue coeff = coeffs_[i]; const IntegerValue div = slack / coeff; diff --git a/ortools/sat/integer_expr.h b/ortools/sat/integer_expr.h index 271ff32a13..b46fc16628 100644 --- a/ortools/sat/integer_expr.h +++ b/ortools/sat/integer_expr.h @@ -84,12 +84,6 @@ class IntegerSumLE : public PropagatorInterface { std::pair ConditionalLb( IntegerLiteral integer_literal, IntegerVariable target_var) const; - // Experimental. Similar to ConditionalLb(), but returns all interesting - // conditional lower bounds instead of just analyzing one integer literal. - // All the IntegerLiteral will be of the form >= var_lb + 1. - std::vector> ConditionalLbs( - IntegerVariable target_var) const; - private: // Fills integer_reason_ with all the current lower_bounds. The real // explanation may require removing one of them, but as an optimization, we diff --git a/ortools/sat/precedences.h b/ortools/sat/precedences.h index 0520db7518..c28c1e16f5 100644 --- a/ortools/sat/precedences.h +++ b/ortools/sat/precedences.h @@ -457,86 +457,6 @@ inline std::function ConditionalLowerOrEqualWithOffset( }; } -// is_le => (a <= b). -inline std::function ConditionalLowerOrEqual(IntegerVariable a, - IntegerVariable b, - Literal is_le) { - return ConditionalLowerOrEqualWithOffset(a, b, 0, is_le); -} - -// literals => (a <= b). -inline std::function ConditionalLowerOrEqual( - IntegerVariable a, IntegerVariable b, absl::Span literals) { - return [=](Model* model) { - PrecedencesPropagator* p = model->GetOrCreate(); - p->AddPrecedenceWithAllOptions(a, b, IntegerValue(0), - /*offset_var*/ kNoIntegerVariable, literals); - }; -} - -// is_le <=> (a + offset <= b). -inline std::function ReifiedLowerOrEqualWithOffset( - IntegerVariable a, IntegerVariable b, int64_t offset, Literal is_le) { - return [=](Model* model) { - PrecedencesPropagator* p = model->GetOrCreate(); - p->AddConditionalPrecedenceWithOffset(a, b, IntegerValue(offset), is_le); - - // The negation of (a + offset <= b) is (a + offset > b) which can be - // rewritten as (b + 1 - offset <= a). - p->AddConditionalPrecedenceWithOffset(b, a, IntegerValue(1 - offset), - is_le.Negated()); - }; -} - -// is_eq <=> (a == b). -inline std::function ReifiedEquality(IntegerVariable a, - IntegerVariable b, - Literal is_eq) { - return [=](Model* model) { - // We creates two extra Boolean variables in this case. - // - // TODO(user): Avoid creating them if we already have some literal that - // have the same meaning. For instance if a client also wanted to know if - // a <= b, he would have called ReifiedLowerOrEqualWithOffset() directly. - const Literal is_le = Literal(model->Add(NewBooleanVariable()), true); - const Literal is_ge = Literal(model->Add(NewBooleanVariable()), true); - model->Add(ReifiedBoolAnd({is_le, is_ge}, is_eq)); - model->Add(ReifiedLowerOrEqualWithOffset(a, b, 0, is_le)); - model->Add(ReifiedLowerOrEqualWithOffset(b, a, 0, is_ge)); - }; -} - -// is_eq <=> (a + offset == b). -inline std::function ReifiedEqualityWithOffset(IntegerVariable a, - IntegerVariable b, - int64_t offset, - Literal is_eq) { - return [=](Model* model) { - // We creates two extra Boolean variables in this case. - // - // TODO(user): Avoid creating them if we already have some literal that - // have the same meaning. For instance if a client also wanted to know if - // a <= b, he would have called ReifiedLowerOrEqualWithOffset() directly. - const Literal is_le = Literal(model->Add(NewBooleanVariable()), true); - const Literal is_ge = Literal(model->Add(NewBooleanVariable()), true); - model->Add(ReifiedBoolAnd({is_le, is_ge}, is_eq)); - model->Add(ReifiedLowerOrEqualWithOffset(a, b, offset, is_le)); - model->Add(ReifiedLowerOrEqualWithOffset(b, a, -offset, is_ge)); - }; -} - -// a != b. -inline std::function NotEqual(IntegerVariable a, - IntegerVariable b) { - return [=](Model* model) { - // We have two options (is_gt or is_lt) and one must be true. - const Literal is_lt = Literal(model->Add(NewBooleanVariable()), true); - const Literal is_gt = is_lt.Negated(); - model->Add(ConditionalLowerOrEqualWithOffset(a, b, 1, is_lt)); - model->Add(ConditionalLowerOrEqualWithOffset(b, a, 1, is_gt)); - }; -} - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/python/cp_model.py b/ortools/sat/python/cp_model.py index f3a5170868..d209588978 100644 --- a/ortools/sat/python/cp_model.py +++ b/ortools/sat/python/cp_model.py @@ -47,6 +47,7 @@ rather than for solving specific optimization problems. import collections import threading import time +from typing import Optional import warnings from ortools.sat import cp_model_pb2 @@ -1411,7 +1412,8 @@ class CpModel(object): model_ct = self.__model.constraints[ct.Index()] model_ct.reservoir.time_exprs.extend( [self.ParseLinearExpression(x) for x in times]) - model_ct.reservoir.level_changes.extend(level_changes) + model_ct.reservoir.level_changes.extend( + [self.ParseLinearExpression(x) for x in level_changes]) model_ct.reservoir.min_level = min_level model_ct.reservoir.max_level = max_level return ct @@ -1476,7 +1478,8 @@ class CpModel(object): model_ct = self.__model.constraints[ct.Index()] model_ct.reservoir.time_exprs.extend( [self.ParseLinearExpression(x) for x in times]) - model_ct.reservoir.level_changes.extend(level_changes) + model_ct.reservoir.level_changes.extend( + [self.ParseLinearExpression(x) for x in level_changes]) model_ct.reservoir.active_literals.extend( [self.GetOrMakeIndex(x) for x in actives]) model_ct.reservoir.min_level = min_level @@ -2176,7 +2179,7 @@ class CpSolver(object): """ def __init__(self): - self.__solution: cp_model_pb2.CpSolverResponse = None + self.__solution: Optional[cp_model_pb2.CpSolverResponse] = None self.parameters = sat_parameters_pb2.SatParameters() self.log_callback = None self.__solve_wrapper: swig_helper.SolveWrapper = None diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index e4886d1680..8552ff586d 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -24,7 +24,7 @@ option csharp_namespace = "Google.OrTools.Sat"; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 223 +// NEXT TAG: 224 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -487,10 +487,19 @@ message SatParameters { // Permutations (#Variables = #Values) are always expanded. optional bool expand_alldiff_constraints = 170 [default = false]; + // If true, expand the reservoir constraints by creating booleans for all + // possible precedences between event and encoding the constraint. + optional bool expand_reservoir_constraints = 182 [default = true]; + // If true, it disable all constraint expansion. // This should only be used to test the presolve of expanded constraints. optional bool disable_constraint_expansion = 181 [default = false]; + // Linear constraint with a complex right hand side (more than a single + // interval) need to be expanded, there is a couple of way to do that. + optional bool encode_complex_linear_constraint_with_integer = 223 + [default = false]; + // During presolve, we use a maximum clique heuristic to merge together // no-overlap constraints or at most one constraints. This code can be slow, // so we have a limit in place on the number of explored nodes in the diff --git a/ortools/sat/scheduling_cuts.cc b/ortools/sat/scheduling_cuts.cc index 285f3efafb..6f4e9f2b26 100644 --- a/ortools/sat/scheduling_cuts.cc +++ b/ortools/sat/scheduling_cuts.cc @@ -27,6 +27,7 @@ #include "absl/strings/str_cat.h" #include "absl/types/span.h" +#include "ortools/base/check.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" #include "ortools/base/strong_vector.h" diff --git a/ortools/sat/timetable.cc b/ortools/sat/timetable.cc index 79edb1d3b9..71e198bddd 100644 --- a/ortools/sat/timetable.cc +++ b/ortools/sat/timetable.cc @@ -30,25 +30,23 @@ namespace operations_research { namespace sat { void AddReservoirConstraint(std::vector times, - std::vector deltas, + std::vector deltas, std::vector presences, int64_t min_level, int64_t max_level, Model* model) { // We only create a side if it can fail. IntegerValue min_possible(0); IntegerValue max_possible(0); - for (const IntegerValue d : deltas) { - if (d > 0) { - max_possible += d; - } else { - min_possible += d; - } + auto* integer_trail = model->GetOrCreate(); + for (const AffineExpression d : deltas) { + min_possible += std::min(IntegerValue(0), integer_trail->LowerBound(d)); + max_possible += std::max(IntegerValue(0), integer_trail->UpperBound(d)); } if (max_possible > max_level) { model->TakeOwnership(new ReservoirTimeTabling( times, deltas, presences, IntegerValue(max_level), model)); } if (min_possible < min_level) { - for (IntegerValue& ref : deltas) ref = -ref; + for (AffineExpression& ref : deltas) ref = ref.Negated(); model->TakeOwnership(new ReservoirTimeTabling( times, deltas, presences, IntegerValue(-min_level), model)); } @@ -56,7 +54,7 @@ void AddReservoirConstraint(std::vector times, ReservoirTimeTabling::ReservoirTimeTabling( const std::vector& times, - const std::vector& deltas, + const std::vector& deltas, const std::vector& presences, IntegerValue capacity, Model* model) : times_(times), deltas_(deltas), @@ -68,11 +66,12 @@ ReservoirTimeTabling::ReservoirTimeTabling( const int id = watcher->Register(this); const int num_events = times.size(); for (int e = 0; e < num_events; e++) { - if (deltas_[e] > 0) { + watcher->WatchLowerBound(deltas_[e], id); + if (integer_trail_->UpperBound(deltas_[e]) > 0) { watcher->WatchUpperBound(times_[e].var, id); watcher->WatchLiteral(presences_[e], id); } - if (deltas_[e] < 0) { + if (integer_trail_->LowerBound(deltas_[e]) < 0) { watcher->WatchLowerBound(times_[e].var, id); watcher->WatchLiteral(presences_[e].Negated(), id); } @@ -86,11 +85,12 @@ bool ReservoirTimeTabling::Propagate() { for (int e = 0; e < num_events; e++) { if (assignment_.LiteralIsFalse(presences_[e])) continue; - // For positive delta, we can maybe increase the min. - if (deltas_[e] > 0 && !TryToIncreaseMin(e)) return false; + // For positive delta_min, we can maybe increase the min. + const IntegerValue min_d = integer_trail_->LowerBound(deltas_[e]); + if (min_d > 0 && !TryToIncreaseMin(e)) return false; - // For negative delta, we can maybe decrease the max. - if (deltas_[e] < 0 && !TryToDecreaseMax(e)) return false; + // For negative delta_min, we can maybe decrease the max. + if (min_d < 0 && !TryToDecreaseMax(e)) return false; } return true; } @@ -105,15 +105,16 @@ bool ReservoirTimeTabling::BuildProfile() { const int num_events = times_.size(); profile_.emplace_back(kMinIntegerValue, IntegerValue(0)); // Sentinel. for (int e = 0; e < num_events; e++) { - if (deltas_[e] > 0) { + const IntegerValue min_d = integer_trail_->LowerBound(deltas_[e]); + if (min_d > 0) { // Only consider present event for positive delta. if (!assignment_.LiteralIsTrue(presences_[e])) continue; const IntegerValue ub = integer_trail_->UpperBound(times_[e]); - profile_.push_back({ub, deltas_[e]}); - } else if (deltas_[e] < 0) { + profile_.push_back({ub, min_d}); + } else if (min_d < 0) { // Only consider non-absent event for negative delta. if (assignment_.LiteralIsFalse(presences_[e])) continue; - profile_.push_back({integer_trail_->LowerBound(times_[e]), deltas_[e]}); + profile_.push_back({integer_trail_->LowerBound(times_[e]), min_d}); } } profile_.emplace_back(kMaxIntegerValue, IntegerValue(0)); // Sentinel. @@ -135,6 +136,7 @@ bool ReservoirTimeTabling::BuildProfile() { // Conflict? for (const ProfileRectangle& rect : profile_) { if (rect.height <= capacity_) continue; + FillReasonForProfileAtGivenTime(rect.start); return integer_trail_->ReportConflict(literal_reason_, integer_reason_); } @@ -142,6 +144,22 @@ bool ReservoirTimeTabling::BuildProfile() { return true; } +namespace { + +void AddLowerOrEqual(const AffineExpression& expr, IntegerValue bound, + std::vector* reason) { + if (expr.IsConstant()) return; + reason->push_back(expr.LowerOrEqual(bound)); +} + +void AddGreaterOrEqual(const AffineExpression& expr, IntegerValue bound, + std::vector* reason) { + if (expr.IsConstant()) return; + reason->push_back(expr.GreaterOrEqual(bound)); +} + +} // namespace + // TODO(user): Minimize with how high the profile needs to be. We can also // remove from the reason the absence of a negative event provided that the // level zero min of the event is greater than t anyway. @@ -155,16 +173,21 @@ void ReservoirTimeTabling::FillReasonForProfileAtGivenTime( const int num_events = times_.size(); for (int e = 0; e < num_events; e++) { if (e == event_to_ignore) continue; - if (deltas_[e] > 0) { + const IntegerValue min_d = integer_trail_->LowerBound(deltas_[e]); + if (min_d > 0) { if (!assignment_.LiteralIsTrue(presences_[e])) continue; if (integer_trail_->UpperBound(times_[e]) > t) continue; - integer_reason_.push_back(times_[e].LowerOrEqual(t)); + AddGreaterOrEqual(deltas_[e], min_d, &integer_reason_); + AddLowerOrEqual(times_[e], t, &integer_reason_); literal_reason_.push_back(presences_[e].Negated()); - } else if (deltas_[e] < 0) { + } else if (min_d <= 0) { if (assignment_.LiteralIsFalse(presences_[e])) { literal_reason_.push_back(presences_[e]); - } else if (integer_trail_->LowerBound(times_[e]) > t) { - integer_reason_.push_back(times_[e].GreaterOrEqual(t + 1)); + continue; + } + AddGreaterOrEqual(deltas_[e], min_d, &integer_reason_); + if (min_d < 0 && integer_trail_->LowerBound(times_[e]) > t) { + AddGreaterOrEqual(times_[e], t + 1, &integer_reason_); } } } @@ -173,7 +196,8 @@ void ReservoirTimeTabling::FillReasonForProfileAtGivenTime( // Note that a negative event will always be in the profile, even if its // presence is still not settled. bool ReservoirTimeTabling::TryToDecreaseMax(int event) { - CHECK_LT(deltas_[event], 0); + const IntegerValue min_d = integer_trail_->LowerBound(deltas_[event]); + CHECK_LT(min_d, 0); const IntegerValue start = integer_trail_->LowerBound(times_[event]); const IntegerValue end = integer_trail_->UpperBound(times_[event]); @@ -194,7 +218,7 @@ bool ReservoirTimeTabling::TryToDecreaseMax(int event) { bool push = false; IntegerValue new_end = end; for (; profile_[rec_id].start < end; ++rec_id) { - if (profile_[rec_id].height - deltas_[event] > capacity_) { + if (profile_[rec_id].height - min_d > capacity_) { new_end = profile_[rec_id].start; push = true; break; @@ -210,7 +234,7 @@ bool ReservoirTimeTabling::TryToDecreaseMax(int event) { // detected at profile construction, but then, since the bound might have been // updated, better be defensive. if (new_end < start) { - integer_reason_.push_back(times_[event].GreaterOrEqual(new_end + 1)); + AddGreaterOrEqual(times_[event], new_end + 1, &integer_reason_); return integer_trail_->ReportConflict(literal_reason_, integer_reason_); } @@ -229,7 +253,8 @@ bool ReservoirTimeTabling::TryToDecreaseMax(int event) { } bool ReservoirTimeTabling::TryToIncreaseMin(int event) { - CHECK_GT(deltas_[event], 0); + const IntegerValue min_d = integer_trail_->LowerBound(deltas_[event]); + CHECK_GT(min_d, 0); const IntegerValue start = integer_trail_->LowerBound(times_[event]); const IntegerValue end = integer_trail_->UpperBound(times_[event]); @@ -252,7 +277,7 @@ bool ReservoirTimeTabling::TryToIncreaseMin(int event) { bool push = false; IntegerValue new_start = start; - if (profile_[rec_id].height + deltas_[event] > capacity_) { + if (profile_[rec_id].height + min_d > capacity_) { if (!assignment_.LiteralIsTrue(presences_[event])) { // Push to false since it wasn't part of the profile and cannot fit. push = true; @@ -265,7 +290,7 @@ bool ReservoirTimeTabling::TryToIncreaseMin(int event) { } if (!push) { for (; profile_[rec_id].start > start; --rec_id) { - if (profile_[rec_id - 1].height + deltas_[event] > capacity_) { + if (profile_[rec_id - 1].height + min_d > capacity_) { push = true; new_start = profile_[rec_id].start; break; @@ -276,6 +301,7 @@ bool ReservoirTimeTabling::TryToIncreaseMin(int event) { // The reason is simply the capacity at new_start - 1; FillReasonForProfileAtGivenTime(new_start - 1, event); + AddGreaterOrEqual(deltas_[event], min_d, &integer_reason_); return integer_trail_->ConditionalEnqueue( presences_[event], times_[event].GreaterOrEqual(new_start), &literal_reason_, &integer_reason_); diff --git a/ortools/sat/timetable.h b/ortools/sat/timetable.h index 38c365465c..54b2175085 100644 --- a/ortools/sat/timetable.h +++ b/ortools/sat/timetable.h @@ -35,7 +35,7 @@ namespace sat { // This instantiate one or more ReservoirTimeTabling class to perform the // propagation. void AddReservoirConstraint(std::vector times, - std::vector deltas, + std::vector deltas, std::vector presences, int64_t min_level, int64_t max_level, Model* model); @@ -49,7 +49,7 @@ void AddReservoirConstraint(std::vector times, class ReservoirTimeTabling : public PropagatorInterface { public: ReservoirTimeTabling(const std::vector& times, - const std::vector& deltas, + const std::vector& deltas, const std::vector& presences, IntegerValue capacity, Model* model); @@ -88,7 +88,7 @@ class ReservoirTimeTabling : public PropagatorInterface { // Input. std::vector times_; - std::vector deltas_; + std::vector deltas_; std::vector presences_; IntegerValue capacity_; diff --git a/ortools/sat/util.cc b/ortools/sat/util.cc index 6fa305d05d..2a0421eb9c 100644 --- a/ortools/sat/util.cc +++ b/ortools/sat/util.cc @@ -425,6 +425,7 @@ void MaxBoundedSubsetSum::Add(int64_t value) { if (s <= bound_) { sums_.push_back(s); current_max_ = std::max(current_max_, s); + if (current_max_ == bound_) return; // Abort } } return; @@ -454,6 +455,73 @@ void MaxBoundedSubsetSum::Add(int64_t value) { current_max_ = bound_; } +void MaxBoundedSubsetSum::AddChoices(absl::Span choices) { + if (DEBUG_MODE) { + for (const int64_t c : choices) { + DCHECK_GE(c, 0); + } + } + + // The max is already reachable or we aborted. + if (current_max_ == bound_) return; + + // Filter out zero and values greater than bound_. + filtered_values_.clear(); + for (const int64_t c : choices) { + if (c == 0 || c > bound_) continue; + filtered_values_.push_back(c); + } + if (filtered_values_.empty()) return; + if (filtered_values_.size() == 1) { + Add(filtered_values_[0]); + return; + } + + // Mode 1: vector of all possible sums (with duplicates). + if (!sums_.empty() && sums_.size() <= kMaxComplexityPerAdd) { + const int old_size = sums_.size(); + for (int i = 0; i < old_size; ++i) { + for (const int64_t value : filtered_values_) { + const int64_t s = sums_[i] + value; + if (s <= bound_) { + sums_.push_back(s); + current_max_ = std::max(current_max_, s); + if (current_max_ == bound_) return; // Abort + } + } + } + return; + } + + // Mode 2: bitset of all possible sums. + if (bound_ <= kMaxComplexityPerAdd) { + if (!sums_.empty()) { + expanded_sums_.assign(bound_ + 1, false); + for (const int64_t s : sums_) { + expanded_sums_[s] = true; + } + sums_.clear(); + } + + // The reverse order is important to not add the current value twice. + for (int64_t i = bound_ - 1; i >= 0; --i) { + if (expanded_sums_[i]) { + for (const int64_t value : filtered_values_) { + if (i + value <= bound_) { + expanded_sums_[i + value] = true; + current_max_ = std::max(current_max_, i + value); + if (current_max_ == bound_) return; // Abort + } + } + } + } + return; + } + + // Abort. + current_max_ = bound_; +} + BasicKnapsackSolver::Result BasicKnapsackSolver::Solve( const std::vector& domains, const std::vector& coeffs, const std::vector& costs, const Domain& rhs) { diff --git a/ortools/sat/util.h b/ortools/sat/util.h index 1a611fc8df..a36f03480e 100644 --- a/ortools/sat/util.h +++ b/ortools/sat/util.h @@ -206,6 +206,9 @@ class MaxBoundedSubsetSum { // Add a value to the base set for which subset sums will be taken. void Add(int64_t value); + // Add a choice of values to the base set for which subset sums will be taken. + void AddChoices(absl::Span choices); + // Returns an upper bound (inclusive) on the maximum sum <= bound_. // This might return bound_ if we aborted the computation. int64_t CurrentMax() const { return current_max_; } @@ -219,6 +222,7 @@ class MaxBoundedSubsetSum { int64_t current_max_; std::vector sums_; std::vector expanded_sums_; + std::vector filtered_values_; }; // Use Dynamic programming to solve a single knapsack. This is used by the diff --git a/ortools/util/parse_proto.cc b/ortools/util/parse_proto.cc new file mode 100644 index 0000000000..8e4e928774 --- /dev/null +++ b/ortools/util/parse_proto.cc @@ -0,0 +1,91 @@ +// Copyright 2010-2022 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/util/parse_proto.h" + +#include +#include + +#include "absl/strings/str_cat.h" +#include "absl/strings/str_split.h" +#include "google/protobuf/io/tokenizer.h" +#include "google/protobuf/text_format.h" + +namespace operations_research { +namespace { + +class TextFormatErrorCollector : public google::protobuf::io::ErrorCollector { + public: + struct CollectedError { + bool warning; + int line; + int column; + std::string message; + }; + + TextFormatErrorCollector() = default; + ~TextFormatErrorCollector() override = default; + + void AddError(const int line, const int column, + const std::string& message) override { + collected_errors_.push_back({/*warning=*/false, line, column, message}); + } + + void AddWarning(const int line, const int column, + const std::string& message) override { + collected_errors_.push_back({/*warning=*/true, line, column, message}); + } + + // Returns a string listing each collected error. When an error is associated + // with a line number and column number that can be found in `value`, that + // error message is shown in context using a caret (^), like: + // { good_field: 1 error_field: "bad" } + // ^ + std::string RenderErrorMessage(const absl::string_view value) { + std::string message; + std::vector value_lines = absl::StrSplit(value, '\n'); + for (const auto& error : collected_errors_) { + if (error.warning) { + absl::StrAppend(&message, "warning: "); + } + absl::StrAppend(&message, error.message, "\n"); + if (error.line >= 0 && error.line < value_lines.size()) { + const std::string& error_line = value_lines[error.line]; + if (error.column >= 0 && error.column < error_line.size()) { + // If possible, use ^ to point at error.column in error_line. + absl::StrAppend(&message, error_line, "\n"); + absl::StrAppend(&message, std::string(error.column, ' '), "^\n"); + } + } + } + return message; + } + + private: + std::vector collected_errors_; +}; + +} // namespace + +bool ParseTextProtoForFlag(const absl::string_view text, + google::protobuf::Message* const message_out, + std::string* const error_out) { + TextFormatErrorCollector errors; + google::protobuf::TextFormat::Parser parser; + parser.RecordErrorsTo(&errors); + const bool success = parser.ParseFromString(std::string(text), message_out); + *error_out = errors.RenderErrorMessage(text); + return success; +} + +} // namespace operations_research diff --git a/ortools/util/parse_proto.h b/ortools/util/parse_proto.h new file mode 100644 index 0000000000..899c736b5a --- /dev/null +++ b/ortools/util/parse_proto.h @@ -0,0 +1,36 @@ +// Copyright 2010-2022 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_UTIL_PARSE_PROTO_H_ +#define OR_TOOLS_UTIL_PARSE_PROTO_H_ + +#include + +#include "absl/strings/string_view.h" +#include "google/protobuf/message.h" + +namespace operations_research { + +// Tries to parse `text` as a text format proto. On a success, stores the result +// in `message_out` and returns true, otherwise, returns `false` with an +// explanation in `error_out`. +// +// NOTE: this API is optimized for implementing AbslParseFlag(). The error +// message will be multiline and is designed to be easily read when printed. +bool ParseTextProtoForFlag(absl::string_view text, + google::protobuf::Message* message_out, + std::string* error_out); + +} // namespace operations_research + +#endif // OR_TOOLS_UTIL_PARSE_PROTO_H_ diff --git a/ortools/util/saturated_arithmetic.h b/ortools/util/saturated_arithmetic.h index 1bef05eea9..c2ec1c1498 100644 --- a/ortools/util/saturated_arithmetic.h +++ b/ortools/util/saturated_arithmetic.h @@ -168,6 +168,11 @@ inline int64_t CapSub(int64_t x, int64_t y) { // Note(user): -kint64min != kint64max, but kint64max == ~kint64min. inline int64_t CapOpp(int64_t v) { return v == kint64min ? ~v : -v; } +inline int64_t CapAbs(int64_t v) { + return v == kint64min ? std::numeric_limits::max() + : (v < 0 ? -v : v); +} + namespace cap_prod_util { // Returns an unsigned int equal to the absolute value of n, in a way that // will not produce overflows.