From 8af04faff5fa214220226d4d1db2c3db5f78123d Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Thu, 16 Jun 2022 07:45:19 +0200 Subject: [PATCH] [CP-SAT] introduce decomposed energy and use it in scheduling propagation and cuts; remove duplicate code in presolve; simplify the name of a few parameters; add SchedulingDemandHelper class; re-enable model status in search log --- ortools/sat/BUILD.bazel | 1 + ortools/sat/cp_model_lns.cc | 5 + ortools/sat/cp_model_presolve.cc | 361 ++++++++--------- ortools/sat/cp_model_presolve.h | 1 - ortools/sat/cp_model_search.cc | 27 +- ortools/sat/cumulative.cc | 23 +- ortools/sat/cumulative_energy.cc | 142 ++----- ortools/sat/cumulative_energy.h | 37 +- ortools/sat/diffn.cc | 13 +- ortools/sat/disjunctive.cc | 7 +- ortools/sat/implied_bounds.cc | 213 +++++----- ortools/sat/implied_bounds.h | 13 +- ortools/sat/integer.cc | 20 + ortools/sat/integer.h | 24 +- ortools/sat/integer_search.cc | 9 +- ortools/sat/intervals.cc | 305 +++++++++++++- ortools/sat/intervals.h | 100 ++++- ortools/sat/lb_tree_search.cc | 5 +- ortools/sat/linear_constraint.cc | 58 +-- ortools/sat/linear_constraint.h | 11 +- ortools/sat/linear_relaxation.cc | 110 +++-- ortools/sat/linear_relaxation.h | 27 +- ortools/sat/lp_utils.cc | 8 +- ortools/sat/model.h | 3 +- ortools/sat/parameters_validation.cc | 1 + ortools/sat/presolve_context.cc | 2 - ortools/sat/presolve_context.h | 11 - ortools/sat/probing.cc | 3 + ortools/sat/probing.h | 10 + ortools/sat/sat_parameters.proto | 24 +- ortools/sat/scheduling_cuts.cc | 579 ++++++++++++++------------- ortools/sat/scheduling_cuts.h | 24 +- ortools/sat/synchronization.cc | 6 +- ortools/sat/synchronization.h | 1 + ortools/sat/timetable.cc | 48 ++- ortools/sat/timetable.h | 23 +- ortools/sat/timetable_edgefinding.cc | 239 ++++++----- ortools/sat/timetable_edgefinding.h | 36 +- 38 files changed, 1456 insertions(+), 1074 deletions(-) diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index c8c58108a1..80152d41aa 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -753,6 +753,7 @@ cc_library( hdrs = ["intervals.h"], deps = [ ":cp_constraints", + ":implied_bounds", ":integer", ":integer_expr", ":model", diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index ee791cf580..fddec14744 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -337,11 +337,16 @@ void NeighborhoodGeneratorHelper::RecomputeHelperData() { absl::StrCat(" compo:", absl::StrJoin(component_sizes, ","), ",..."); } } + + // TODO(user): This is not ideal, as if two reductions appears in a row and + // nothing else is done for a while, we will never see the "latest" size + // in the log until it is reduced again. shared_response_->LogPeriodicMessage( "Model", absl::StrCat("var:", active_variables_.size(), "/", num_variables, " constraints:", simplied_model_proto_.constraints().size(), "/", model_proto_.constraints().size(), compo_message), + parameters_.model_reduction_log_frequency_in_seconds(), &last_logging_time_); } diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 32d0a7a40a..262a73c69d 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -2372,22 +2372,6 @@ void CpModelPresolver::TryToReduceCoefficientsOfLinearConstraint( namespace { -// Return true if the given domain only restrict the values with an upper bound. -bool IsLeConstraint(const Domain& domain, const Domain& all_values) { - return all_values - .IntersectionWith( - Domain(std::numeric_limits::min(), domain.Max())) - .IsIncludedIn(domain); -} - -// Same as IsLeConstraint() but in the other direction. -bool IsGeConstraint(const Domain& domain, const Domain& all_values) { - return all_values - .IntersectionWith( - Domain(domain.Min(), std::numeric_limits::max())) - .IsIncludedIn(domain); -} - // In the equation terms + coeff * var_domain \included rhs, returns true if can // we always fix rhs to its min value for any value in terms. It is okay to // not be as generic as possible here. @@ -2425,138 +2409,8 @@ bool RhsCanBeFixedToMax(int64_t coeff, const Domain& var_domain, return false; } -// Remove from to_clear any entry not in current. -void TakeIntersectionWith(const absl::flat_hash_set& current, - absl::flat_hash_set* to_clear) { - std::vector new_set; - for (const int c : *to_clear) { - if (current.contains(c)) new_set.push_back(c); - } - to_clear->clear(); - for (const int c : new_set) to_clear->insert(c); -} - } // namespace -bool CpModelPresolver::DetectAndProcessOneSidedLinearConstraint( - int c, ConstraintProto* ct) { - if (ct->constraint_case() != ConstraintProto::kLinear) return false; - if (context_->ModelIsUnsat()) return false; - if (context_->keep_all_feasible_solutions) return false; - - // TODO(user): There is a bit of code and effort duplication with - // PropagateDomainsInLinear(). Try to remove that. - Domain implied_rhs(0); - const int num_vars = ct->linear().vars().size(); - for (int i = 0; i < num_vars; ++i) { - const int ref = ct->linear().vars(i); - const int64_t coeff = ct->linear().coeffs(i); - implied_rhs = - implied_rhs - .AdditionWith(context_->DomainOf(ref).MultiplicationBy(coeff)) - .RelaxIfTooComplex(); - } - - // Abort if trivial. - const Domain old_rhs = ReadDomainFromProto(ct->linear()); - if (implied_rhs.IsIncludedIn(old_rhs)) { - context_->UpdateRuleStats("linear: always true"); - return RemoveConstraint(ct); - } - - // Incorporate the implied rhs information. - const Domain rhs = old_rhs.SimplifyUsingImpliedDomain(implied_rhs); - if (rhs.IsEmpty()) { - context_->UpdateRuleStats("linear: infeasible"); - return MarkConstraintAsFalse(ct); - } - if (rhs != old_rhs) { - context_->UpdateRuleStats("linear: simplified rhs"); - } - FillDomainInProto(rhs, ct->mutable_linear()); - - // Detect if it is always good for a term of this constraint to move towards - // its lower (resp. upper) bound. This is the same as saying that this - // constraint only bound in one direction. - const bool is_le_constraint = IsLeConstraint(rhs, implied_rhs); - const bool is_ge_constraint = IsGeConstraint(rhs, implied_rhs); - if (!is_le_constraint && !is_ge_constraint) return false; - CHECK_NE(is_le_constraint, is_ge_constraint); - - // Tricky: If a variable appears in the enforcement literal list, it is - // constraining in the "true" direction. We rely on the constraint to have - // been canonicalized and these literal removed for correctness here. In debug - // mode we still do some check since it is easy to add presolve rules and call - // this with a non-canonicalized constraint. - absl::flat_hash_set enforcement_set; - if (DEBUG_MODE) { - for (const int ref : ct->enforcement_literal()) { - enforcement_set.insert(ref); - } - } - - bool recanonicalize = false; - for (int i = 0; i < num_vars; ++i) { - const int var = ct->linear().vars(i); - const int64_t var_coeff = ct->linear().coeffs(i); - CHECK(RefIsPositive(var)); - - if ((var_coeff > 0) == is_ge_constraint) { - DCHECK(!enforcement_set.contains(var)); - context_->var_to_lb_only_constraints[var].insert(c); - } else { - DCHECK(!enforcement_set.contains(NegatedRef(var))); - context_->var_to_ub_only_constraints[var].insert(c); - } - - // Simple dual fixing: If for any feasible solution, any solution with var - // higher (resp. lower) is also valid, then we can fix that variable to - // its bound if it also moves the objective in the good direction. - const bool is_in_objective = context_->VarToConstraints(var).contains(-1); - const int size = - context_->VarToConstraints(var).size() - (is_in_objective ? 1 : 0); - const int64_t obj_coeff = - is_in_objective ? context_->ObjectiveMap().at(var) : 0; - - // We cannot fix anything if the domain of the objective is excluding - // some objective values. - if (obj_coeff != 0 && context_->ObjectiveDomainIsConstraining()) { - continue; - } - - if (obj_coeff <= 0 && - context_->var_to_lb_only_constraints[var].size() >= size) { - TakeIntersectionWith(context_->VarToConstraints(var), - &(context_->var_to_lb_only_constraints[var])); - if (context_->var_to_lb_only_constraints[var].size() >= size) { - if (!context_->IntersectDomainWith(var, Domain(context_->MaxOf(var)))) { - return false; - } - context_->UpdateRuleStats("linear: dual fixing"); - recanonicalize = true; - continue; - } - } - - if (obj_coeff >= 0 && - context_->var_to_ub_only_constraints[var].size() >= size) { - TakeIntersectionWith(context_->VarToConstraints(var), - &(context_->var_to_ub_only_constraints[var])); - if (context_->var_to_ub_only_constraints[var].size() >= size) { - if (!context_->IntersectDomainWith(var, Domain(context_->MinOf(var)))) { - return false; - } - context_->UpdateRuleStats("linear: dual fixing"); - recanonicalize = true; - continue; - } - } - } - - if (recanonicalize) return CanonicalizeLinear(ct); - return false; -} - // TODO(user): generalize to disjoint amo for even better propagation! // TODO(user): merge code with the fully included case. // TODO(user): Similarly amo and bool_or intersection or amo and enforcement @@ -2834,7 +2688,7 @@ bool CpModelPresolver::PropagateDomainsInLinear(int ct_index, for (int i = 0; i < num_vars; ++i) { const int var = ct->linear().vars(i); const int64_t coeff = ct->linear().coeffs(i); - CHECK(RefIsPositive(var)); + DCHECK(RefIsPositive(var)); term_domains[i] = context_->DomainOf(var).MultiplicationBy(coeff); left_domains[i + 1] = left_domains[i].AdditionWith(term_domains[i]).RelaxIfTooComplex(); @@ -2905,7 +2759,7 @@ bool CpModelPresolver::PropagateDomainsInLinear(int ct_index, // Given a variable that only appear in one constraint and in the // objective, for any feasible solution, it will be always better to move // this singleton variable as much as possible towards its good objective - // direction. Sometime_exprs, we can detect that we will always be able to + // direction. Sometime, we can detect that we will always be able to // do this until the only constraint of this singleton variable is tight. // // When this happens, we can make the constraint an equality. Note that it @@ -5555,9 +5409,167 @@ void CpModelPresolver::Probe() { auto* sat_solver = model.GetOrCreate(); auto* mapping = model.GetOrCreate(); auto* prober = model.GetOrCreate(); + + // Try to detect trivial clauses thanks to implications. + // This can be slow, so we bound the amount of work done. + // + // Idea: If we have l1, l2 in a bool_or and not(l1) => l2, the constraint is + // always true. + // + // Correctness: Note that we always replace a clause with another one that + // subsumes it. So we are correct even if new clauses are learned and used + // for propagation along the way. + // + // TODO(user): Improve the algo? + int64_t work_done = 0; + const int64_t work_limit = 1e8; + if (true) { + const auto& assignment = sat_solver->Assignment(); + prober->SetPropagationCallback([&](Literal decision) { + if (work_done > work_limit) return; + const int decision_var = + mapping->GetProtoVariableFromBooleanVariable(decision.Variable()); + if (decision_var < 0) return; + for (const int c : context_->VarToConstraints(decision_var)) { + ++work_done; + if (c < 0) continue; + const ConstraintProto& ct = context_->working_model->constraints(c); + if (ct.enforcement_literal().size() > 2) { + // Any l for which decision => l can be removed. + // + // If decision => not(l), constraint can never be satisfied. However + // because we don't know if this constraint was part of the + // propagation we replace it by an implication. + // + // TODO(user): remove duplication with code below. + // TODO(user): If decision appear positively, we could potentially + // remove a bunch of terms (all the ones involving variables implied + // by the decision) from the innner constraint, especially in the + // linear case. + int decision_ref; + int false_ref; + bool decision_is_positive = false; + bool has_false_literal = false; + bool simplification_possible = false; + for (const int ref : ct.enforcement_literal()) { + ++work_done; + const Literal lit = mapping->Literal(ref); + if (PositiveRef(ref) == decision_var) { + decision_ref = ref; + decision_is_positive = assignment.LiteralIsTrue(lit); + if (!decision_is_positive) break; + continue; + } + if (assignment.LiteralIsFalse(lit)) { + false_ref = ref; + has_false_literal = true; + } else if (assignment.LiteralIsTrue(lit)) { + // If decision => l, we can remove l from the list. + simplification_possible = true; + } + } + if (!decision_is_positive) continue; + + if (has_false_literal) { + // Reduce to implication. + auto* mutable_ct = context_->working_model->mutable_constraints(c); + mutable_ct->Clear(); + mutable_ct->add_enforcement_literal(decision_ref); + mutable_ct->mutable_bool_and()->add_literals(NegatedRef(false_ref)); + context_->UpdateRuleStats( + "probing: reduced enforced constraint to implication."); + context_->UpdateConstraintVariableUsage(c); + continue; + } + + if (simplification_possible) { + int new_size = 0; + auto* mutable_enforcements = + context_->working_model->mutable_constraints(c) + ->mutable_enforcement_literal(); + for (const int ref : ct.enforcement_literal()) { + if (PositiveRef(ref) != decision_var && + assignment.LiteralIsTrue(mapping->Literal(ref))) { + continue; + } + mutable_enforcements->Set(new_size++, ref); + } + mutable_enforcements->Truncate(new_size); + context_->UpdateRuleStats("probing: simplified enforcement list."); + context_->UpdateConstraintVariableUsage(c); + } + continue; + } + + if (ct.constraint_case() != ConstraintProto::kBoolOr) continue; + if (ct.bool_or().literals().size() <= 2) continue; + + int decision_ref; + int true_ref; + bool decision_is_negative = false; + bool has_true_literal = false; + bool simplification_possible = false; + for (const int ref : ct.bool_or().literals()) { + ++work_done; + const Literal lit = mapping->Literal(ref); + if (PositiveRef(ref) == decision_var) { + decision_ref = ref; + decision_is_negative = assignment.LiteralIsFalse(lit); + if (!decision_is_negative) break; + continue; + } + if (assignment.LiteralIsTrue(lit)) { + true_ref = ref; + has_true_literal = true; + } else if (assignment.LiteralIsFalse(lit)) { + // If not(l1) => not(l2), we can remove l2 from the clause. + simplification_possible = true; + } + } + if (!decision_is_negative) continue; + + if (has_true_literal) { + // This will later be merged with the current implications and removed + // if it is a duplicate. + auto* mutable_bool_or = + context_->working_model->mutable_constraints(c) + ->mutable_bool_or(); + mutable_bool_or->mutable_literals()->Clear(); + mutable_bool_or->add_literals(decision_ref); + mutable_bool_or->add_literals(true_ref); + context_->UpdateRuleStats("probing: bool_or reduced to implication"); + context_->UpdateConstraintVariableUsage(c); + continue; + } + + if (simplification_possible) { + int new_size = 0; + auto* mutable_bool_or = + context_->working_model->mutable_constraints(c) + ->mutable_bool_or(); + for (const int ref : ct.bool_or().literals()) { + if (PositiveRef(ref) != decision_var && + assignment.LiteralIsFalse(mapping->Literal(ref))) { + continue; + } + mutable_bool_or->set_literals(new_size++, ref); + } + mutable_bool_or->mutable_literals()->Truncate(new_size); + context_->UpdateRuleStats("probing: simplified clauses."); + context_->UpdateConstraintVariableUsage(c); + } + } + }); + } + prober->ProbeBooleanVariables(/*deterministic_time_limit=*/1.0); context_->time_limit()->AdvanceDeterministicTime( model.GetOrCreate()->GetElapsedDeterministicTime()); + if (work_done > 0) { + SOLVER_LOG(logger_, + "[Probing] implications and bool_or (work_done=", work_done, + ").", (work_done > work_limit ? " Aborted." : "")); + } if (sat_solver->IsModelUnsat() || !implication_graph->DetectEquivalences()) { return (void)context_->NotifyThatModelIsUnsat("during probing"); } @@ -6258,20 +6270,6 @@ bool CpModelPresolver::PresolveOneConstraint(int c) { } return PresolveIntMod(ct); case ConstraintProto::kLinear: { - // In the presence of affine relation, it is possible that the sign of a - // variable change during canonicalization, and a variable that could - // freely move in one direction can no longer do so. So we make sure we - // always remove c from all the maps before re-inserting it in - // DetectAndProcessOneSidedLinearConstraint(). - // - // TODO(user): The move in only one direction code is redundant with the - // dual bound strengthening code. So maybe we don't need both. - for (const int ref : ct->linear().vars()) { - const int var = PositiveRef(ref); - context_->var_to_lb_only_constraints[var].erase(c); - context_->var_to_ub_only_constraints[var].erase(c); - } - if (CanonicalizeLinear(ct)) { context_->UpdateConstraintVariableUsage(c); } @@ -6301,7 +6299,7 @@ bool CpModelPresolver::PresolveOneConstraint(int c) { context_->UpdateConstraintVariableUsage(c); } - // If we extracted some enforcement, we redo some preoslve. + // If we extracted some enforcement, we redo some presolve. const int old_num_enforcement_literals = ct->enforcement_literal_size(); ExtractEnforcementLiteralFromLinearConstraint(c, ct); if (ct->enforcement_literal_size() > old_num_enforcement_literals) { @@ -6314,14 +6312,6 @@ bool CpModelPresolver::PresolveOneConstraint(int c) { } TryToReduceCoefficientsOfLinearConstraint(c, ct); - - // Note that it is important for this to be last, so that if a constraint - // is marked as being in one direction, no other presolve is applied until - // it is processed again and unmarked at the beginning of this case. - if (DetectAndProcessOneSidedLinearConstraint(c, ct)) { - context_->UpdateConstraintVariableUsage(c); - } - return false; } case ConstraintProto::kInterval: @@ -6843,16 +6833,6 @@ void CpModelPresolver::DetectDuplicateConstraints() { const Domain d = ReadDomainFromProto( context_->working_model->constraints(dup).linear()); if (rep_domain != d) { - // Tricky: we modify the domain so we need to clear this information. - // - // TODO(user): Revisit this algorithm and integrate it with the dual - // reduction. We can either have incrementaly maintained info, or just - // do it with one pass of DualBoundStrengthening. - for (const int var : context_->ConstraintToVars(rep)) { - context_->var_to_lb_only_constraints[var].erase(rep); - context_->var_to_ub_only_constraints[var].erase(rep); - } - context_->UpdateRuleStats("duplicate: merged rhs of linear constraint"); const Domain rhs = rep_domain.IntersectionWith(d); if (rhs.IsEmpty()) { @@ -6895,11 +6875,6 @@ void CpModelPresolver::DetectDuplicateConstraints() { context_->ReadObjectiveFromProto(); } } - - for (const int var : context_->ConstraintToVars(dup)) { - context_->var_to_lb_only_constraints[var].erase(dup); - context_->var_to_ub_only_constraints[var].erase(dup); - } context_->working_model->mutable_constraints(dup)->Clear(); context_->UpdateConstraintVariableUsage(dup); context_->UpdateRuleStats("duplicate: removed constraint"); @@ -7128,10 +7103,6 @@ void CpModelPresolver::DetectDominatedLinearConstraints() { }); for (const int c : constraint_indices_to_clean) { - for (const int var : context_->ConstraintToVars(c)) { - context_->var_to_lb_only_constraints[var].erase(c); - context_->var_to_ub_only_constraints[var].erase(c); - } context_->UpdateConstraintVariableUsage(c); } @@ -8562,10 +8533,12 @@ bool ModelCopy::CopyLinear(const ConstraintProto& ct) { const Domain new_domain = ReadDomainFromProto(ct.linear()).AdditionWith(Domain(-offset)); - if (non_fixed_variables_.empty() && !new_domain.Contains(0)) { - if (ct.enforcement_literal().empty()) { - return false; - } + if (non_fixed_variables_.empty()) { + // Trivial constraint. + if (new_domain.Contains(0)) return true; + + // Constraint is false. + if (ct.enforcement_literal().empty()) return false; temp_literals_.clear(); for (const int literal : ct.enforcement_literal()) { if (context_->LiteralIsTrue(literal)) { diff --git a/ortools/sat/cp_model_presolve.h b/ortools/sat/cp_model_presolve.h index 23c1eb3dce..055a510994 100644 --- a/ortools/sat/cp_model_presolve.h +++ b/ortools/sat/cp_model_presolve.h @@ -156,7 +156,6 @@ class CpModelPresolver { bool AddVarAffineRepresentativeFromLinearEquality(int target_index, ConstraintProto* ct); bool PresolveLinearEqualityWithModulo(ConstraintProto* ct); - bool DetectAndProcessOneSidedLinearConstraint(int c, ConstraintProto* ct); // It can be interesting to know for a given linear constraint that a subset // of its variables are in at most one relation. diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index 5bc32eff69..080fae5d23 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -483,9 +483,9 @@ std::vector GetDiverseSetOfParameters( new_params.set_optimize_with_lb_tree_search(true); new_params.set_linearization_level(2); if (base_params.use_dual_scheduling_heuristics()) { - new_params.set_use_overload_checker_in_cumulative_constraint(true); - new_params.set_use_timetable_edge_finding_in_cumulative_constraint(true); - new_params.set_use_hard_precedences_in_cumulative_constraint(true); + new_params.set_use_overload_checker_in_cumulative(true); + new_params.set_use_timetable_edge_finding_in_cumulative(true); + new_params.set_use_hard_precedences_in_cumulative(true); } // We do not want to change the objective_var lb from outside as it gives @@ -499,12 +499,15 @@ std::vector GetDiverseSetOfParameters( new_params.set_search_branching(SatParameters::AUTOMATIC_SEARCH); new_params.set_use_probing_search(true); if (base_params.use_dual_scheduling_heuristics()) { - new_params.set_use_overload_checker_in_cumulative_constraint(true); - new_params.set_use_timetable_edge_finding_in_cumulative_constraint(true); - new_params.set_use_hard_precedences_in_cumulative_constraint(true); + new_params.set_use_overload_checker_in_cumulative(true); + new_params.set_use_timetable_edge_finding_in_cumulative(true); + new_params.set_use_hard_precedences_in_cumulative(true); } strategies["probing"] = new_params; + new_params.set_linearization_level(1); + strategies["probing_default_lp"] = new_params; + new_params.set_linearization_level(2); strategies["probing_max_lp"] = new_params; } @@ -533,9 +536,9 @@ std::vector GetDiverseSetOfParameters( new_params.set_linearization_level(2); new_params.set_search_branching(SatParameters::LP_SEARCH); if (base_params.use_dual_scheduling_heuristics()) { - new_params.set_use_overload_checker_in_cumulative_constraint(true); - new_params.set_use_timetable_edge_finding_in_cumulative_constraint(true); - new_params.set_use_hard_precedences_in_cumulative_constraint(true); + new_params.set_use_overload_checker_in_cumulative(true); + new_params.set_use_timetable_edge_finding_in_cumulative(true); + new_params.set_use_hard_precedences_in_cumulative(true); } strategies["reduced_costs"] = new_params; } @@ -608,7 +611,6 @@ std::vector GetDiverseSetOfParameters( names.push_back("quick_restart_no_lp"); names.push_back("lb_tree_search"); names.push_back("probing"); - if (base_params.num_workers() > 16) names.push_back("probing_max_lp"); #if !defined(__PORTABLE_PLATFORM__) && defined(USE_SCIP) if (absl::GetFlag(FLAGS_cp_model_use_max_hs)) names.push_back("max_hs"); #endif // !defined(__PORTABLE_PLATFORM__) && defined(USE_SCIP) @@ -618,6 +620,11 @@ std::vector GetDiverseSetOfParameters( } } + // Add subsolvers. + for (const std::string& name : base_params.add_subsolvers()) { + names.push_back(name); + } + // Remove the names that should be ignored. absl::flat_hash_set to_ignore; for (const std::string& name : base_params.ignore_subsolvers()) { diff --git a/ortools/sat/cumulative.cc b/ortools/sat/cumulative.cc index cac0341684..ce9ba618cc 100644 --- a/ortools/sat/cumulative.cc +++ b/ortools/sat/cumulative.cc @@ -87,7 +87,7 @@ std::function Cumulative( // Detect a subset of intervals that needs to be in disjunction and add a // Disjunctive() constraint over them. - if (parameters.use_disjunctive_constraint_in_cumulative_constraint()) { + if (parameters.use_disjunctive_constraint_in_cumulative()) { // TODO(user): We need to exclude intervals that can be of size zero // because the disjunctive do not "ignore" them like the cumulative // does. That is, the interval [2,2) will be assumed to be in @@ -148,6 +148,9 @@ std::function Cumulative( if (helper == nullptr) { helper = intervals->GetOrCreateHelper(vars); } + SchedulingDemandHelper* demands_helper = + new SchedulingDemandHelper(demands, helper, model); + model->TakeOwnership(demands_helper); // For each variables that is after a subset of task ends (i.e. like a // makespan objective), we detect it and add a special constraint to @@ -167,7 +170,7 @@ std::function Cumulative( // TODO(user): We compute this only once, so we should explore the full // precedence graph, not just task in direct precedence. Make sure not to // create to many such constraints though. - if (parameters.use_hard_precedences_in_cumulative_constraint()) { + if (parameters.use_hard_precedences_in_cumulative()) { // The CumulativeIsAfterSubsetConstraint() always reset the helper to the // forward time direction, so it is important to also precompute the // precedence relation using the same direction! This is needed in case @@ -218,8 +221,8 @@ std::function Cumulative( sum_of_demand_max > integer_trail->LevelZeroLowerBound(capacity)) { CumulativeIsAfterSubsetConstraint* constraint = new CumulativeIsAfterSubsetConstraint(var, min_offset, capacity, - demands, subtasks, - integer_trail, helper); + subtasks, helper, + demands_helper, model); constraint->RegisterWith(watcher); model->TakeOwnership(constraint); } @@ -230,22 +233,22 @@ std::function Cumulative( // increases the minimum of the start variables, decrease the maximum of the // end variables, and increase the minimum of the capacity variable. TimeTablingPerTask* time_tabling = - new TimeTablingPerTask(demands, capacity, integer_trail, helper); + new TimeTablingPerTask(capacity, helper, demands_helper, model); time_tabling->RegisterWith(watcher); model->TakeOwnership(time_tabling); // Propagator responsible for applying the Overload Checking filtering rule. // It increases the minimum of the capacity variable. - if (parameters.use_overload_checker_in_cumulative_constraint()) { - AddCumulativeOverloadChecker(demands, capacity, helper, model); + if (parameters.use_overload_checker_in_cumulative()) { + AddCumulativeOverloadChecker(capacity, helper, demands_helper, model); } // Propagator responsible for applying the Timetable Edge finding filtering // rule. It increases the minimum of the start variables and decreases the // maximum of the end variables, - if (parameters.use_timetable_edge_finding_in_cumulative_constraint()) { - TimeTableEdgeFinding* time_table_edge_finding = new TimeTableEdgeFinding( - demands, capacity, helper, integer_trail, model); + if (parameters.use_timetable_edge_finding_in_cumulative()) { + TimeTableEdgeFinding* time_table_edge_finding = + new TimeTableEdgeFinding(capacity, helper, demands_helper, model); time_table_edge_finding->RegisterWith(watcher); model->TakeOwnership(time_table_edge_finding); } diff --git a/ortools/sat/cumulative_energy.cc b/ortools/sat/cumulative_energy.cc index 19c3a1531a..48c709c6b6 100644 --- a/ortools/sat/cumulative_energy.cc +++ b/ortools/sat/cumulative_energy.cc @@ -29,75 +29,26 @@ namespace operations_research { namespace sat { -void AddCumulativeEnergyConstraint(std::vector energies, - AffineExpression capacity, - SchedulingConstraintHelper* helper, - Model* model) { - auto* watcher = model->GetOrCreate(); - auto* integer_trail = model->GetOrCreate(); - - CumulativeEnergyConstraint* constraint = new CumulativeEnergyConstraint( - std::move(energies), capacity, integer_trail, helper); - constraint->RegisterWith(watcher); - model->TakeOwnership(constraint); -} - -void AddCumulativeOverloadChecker(const std::vector& demands, - AffineExpression capacity, +void AddCumulativeOverloadChecker(AffineExpression capacity, SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, Model* model) { auto* watcher = model->GetOrCreate(); - auto* integer_trail = model->GetOrCreate(); - - std::vector energies; - const int num_tasks = helper->NumTasks(); - CHECK_EQ(demands.size(), num_tasks); - for (int t = 0; t < num_tasks; ++t) { - const AffineExpression size = helper->Sizes()[t]; - const AffineExpression demand = demands[t]; - // This will cover the basic cases (constant * constant, constant * affine), - // as well as the case where both demand and size are controlled by the same - // set of literals. This is the case for the RCPSP problems. - const std::optional linearized_energy = - TryToLinearizeProduct(size, demand, model); - if (!linearized_energy.has_value()) { - // The case where both demand and size are variable should be rare. - // - // TODO(user): Handle when needed by creating an intermediate product - // variable equal to demand * size. Note that because of the affine - // expression, we do need some custom code for this. - VLOG(1) << "Overload checker with non-linearizable energy is currently " - "not implemented. Skipping."; - return; - } - energies.emplace_back(linearized_energy.value()); - } - CumulativeEnergyConstraint* constraint = - new CumulativeEnergyConstraint(energies, capacity, integer_trail, helper); + new CumulativeEnergyConstraint(capacity, helper, demands, model); constraint->RegisterWith(watcher); model->TakeOwnership(constraint); } CumulativeEnergyConstraint::CumulativeEnergyConstraint( - std::vector energies, AffineExpression capacity, - IntegerTrail* integer_trail, SchedulingConstraintHelper* helper) - : energies_(std::move(energies)), - capacity_(capacity), - integer_trail_(integer_trail), + AffineExpression capacity, SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, Model* model) + : capacity_(capacity), + integer_trail_(model->GetOrCreate()), helper_(helper), - theta_tree_() { + demands_(demands) { const int num_tasks = helper_->NumTasks(); - CHECK_EQ(energies_.size(), num_tasks); task_to_start_event_.resize(num_tasks); - if (DEBUG_MODE) { - for (const LinearExpression& energy : energies_) { - DCHECK_GE(energy.offset, 0); - for (const IntegerValue coeff : energy.coeffs) { - DCHECK_GE(coeff, 0); - } - } - } } void CumulativeEnergyConstraint::RegisterWith(GenericLiteralWatcher* watcher) { @@ -110,6 +61,7 @@ bool CumulativeEnergyConstraint::Propagate() { // This only uses one time direction, but the helper might be used elsewhere. // TODO(user): just keep the current direction? if (!helper_->SynchronizeAndSetTimeDirection(true)) return false; + demands_->CacheAllEnergyValues(); const IntegerValue capacity_max = integer_trail_->UpperBound(capacity_); // TODO(user): force capacity_max >= 0, fail/remove optionals when 0. @@ -120,7 +72,7 @@ bool CumulativeEnergyConstraint::Propagate() { int num_events = 0; for (const auto task_time : helper_->TaskByIncreasingStartMin()) { const int task = task_time.task_index; - if (helper_->IsAbsent(task) || energies_[task].Max(*integer_trail_) == 0) { + if (helper_->IsAbsent(task) || demands_->EnergyMax(task) == 0) { task_to_start_event_[task] = -1; continue; } @@ -148,14 +100,13 @@ bool CumulativeEnergyConstraint::Propagate() { start_event_is_present_[current_event] = is_present; if (is_present) { tree_has_mandatory_intervals = true; - theta_tree_.AddOrUpdateEvent( - current_event, start_min * capacity_max, - energies_[current_task].Min(*integer_trail_), - energies_[current_task].Max(*integer_trail_)); + theta_tree_.AddOrUpdateEvent(current_event, start_min * capacity_max, + demands_->EnergyMin(current_task), + demands_->EnergyMax(current_task)); } else { - theta_tree_.AddOrUpdateOptionalEvent( - current_event, start_min * capacity_max, - energies_[current_task].Max(*integer_trail_)); + theta_tree_.AddOrUpdateOptionalEvent(current_event, + start_min * capacity_max, + demands_->EnergyMax(current_task)); } } @@ -184,10 +135,7 @@ bool CumulativeEnergyConstraint::Propagate() { if (start_event_is_present_[event]) { const int task = start_event_task_time_[event].task_index; helper_->AddPresenceReason(task); - for (const IntegerVariable var : energies_[task].vars) { - helper_->MutableIntegerReason()->push_back( - integer_trail_->LowerBoundAsLiteral(var)); - } + demands_->AddEnergyMinReason(task); helper_->AddStartMinReason(task, window_start); helper_->AddEndMaxReason(task, window_end); } @@ -235,10 +183,7 @@ bool CumulativeEnergyConstraint::Propagate() { helper_->AddPresenceReason(task); helper_->AddStartMinReason(task, window_start); helper_->AddEndMaxReason(task, window_end); - for (const IntegerVariable var : energies_[task].vars) { - helper_->MutableIntegerReason()->push_back( - integer_trail_->LowerBoundAsLiteral(var)); - } + demands_->AddEnergyMinReason(task); } } if (capacity_.var != kNoIntegerVariable) { @@ -250,27 +195,9 @@ bool CumulativeEnergyConstraint::Propagate() { start_event_task_time_[event_with_new_energy_max].task_index; helper_->AddStartMinReason(task_with_new_energy_max, window_start); helper_->AddEndMaxReason(task_with_new_energy_max, window_end); - - if (new_energy_max < - energies_[task_with_new_energy_max].Min(*integer_trail_)) { - if (helper_->IsOptional(task_with_new_energy_max)) { - return helper_->PushTaskAbsence(task_with_new_energy_max); - } else { - return helper_->ReportConflict(); - } - } else if (energies_[task_with_new_energy_max].vars.size() == 1) { - const LinearExpression& e = energies_[task_with_new_energy_max]; - const AffineExpression affine_energy(e.vars[0], e.coeffs[0], e.offset); - const IntegerLiteral deduction = - affine_energy.LowerOrEqual(new_energy_max); - if (!helper_->PushIntegerLiteralIfTaskPresent(task_with_new_energy_max, - deduction)) { - return false; - } - } else { - // TODO(user): Propagate if possible. - VLOG(3) << "Cumulative energy missed propagation " - << energies_[task_with_new_energy_max].vars.size(); + if (!demands_->DecreaseEnergyMax(task_with_new_energy_max, + new_energy_max)) { + return false; } if (helper_->IsPresent(task_with_new_energy_max)) { @@ -278,8 +205,7 @@ bool CumulativeEnergyConstraint::Propagate() { task_to_start_event_[task_with_new_energy_max], start_event_task_time_[event_with_new_energy_max].time * capacity_max, - energies_[task_with_new_energy_max].Min(*integer_trail_), - new_energy_max); + demands_->EnergyMin(task_with_new_energy_max), new_energy_max); } else { theta_tree_.RemoveEvent(event_with_new_energy_max); } @@ -290,16 +216,15 @@ bool CumulativeEnergyConstraint::Propagate() { CumulativeIsAfterSubsetConstraint::CumulativeIsAfterSubsetConstraint( IntegerVariable var, IntegerValue offset, AffineExpression capacity, - const std::vector demands, - const std::vector subtasks, IntegerTrail* integer_trail, - SchedulingConstraintHelper* helper) + const std::vector subtasks, SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, Model* model) : var_to_push_(var), offset_(offset), capacity_(capacity), - demands_(demands), subtasks_(subtasks), - integer_trail_(integer_trail), - helper_(helper) { + integer_trail_(model->GetOrCreate()), + helper_(helper), + demands_(demands) { is_in_subtasks_.assign(helper->NumTasks(), false); for (const int t : subtasks) is_in_subtasks_[t] = true; } @@ -345,11 +270,11 @@ bool CumulativeIsAfterSubsetConstraint::Propagate() { const int t = profile[i].task; if (!helper_->IsPresent(t) || !is_in_subtasks_[t]) continue; - const IntegerValue demand_min = integer_trail_->LowerBound(demands_[t]); + const IntegerValue demand_min = demands_->DemandMin(t); const IntegerValue delta = profile[i].is_first ? -demand_min : demand_min; profile_height += delta; if (delta > 0) { - if (demands_[t].IsConstant()) { + if (demands_->Demands()[t].IsConstant()) { dp_.Add(delta.value()); } else { dp_.Add(capacity_max.value()); // Abort DP. @@ -387,7 +312,7 @@ bool CumulativeIsAfterSubsetConstraint::Propagate() { const IntegerValue size_min = helper_->SizeMin(t); if (size_min == 0) continue; - const IntegerValue demand_min = integer_trail_->LowerBound(demands_[t]); + const IntegerValue demand_min = demands_->DemandMin(t); if (demand_min == 0) continue; const IntegerValue end_min = helper_->EndMin(t); @@ -396,10 +321,7 @@ bool CumulativeIsAfterSubsetConstraint::Propagate() { helper_->AddEndMinReason(t, std::min(best_time + size_min, end_min)); helper_->AddSizeMinReason(t); helper_->AddPresenceReason(t); - if (demands_[t].var != kNoIntegerVariable) { - helper_->MutableIntegerReason()->push_back( - integer_trail_->LowerBoundAsLiteral(demands_[t].var)); - } + demands_->AddDemandMinReason(t); } if (capacity_.var != kNoIntegerVariable) { helper_->MutableIntegerReason()->push_back( @@ -425,7 +347,7 @@ void CumulativeIsAfterSubsetConstraint::RegisterWith( watcher->WatchLowerBound(helper_->Starts()[t], id); watcher->WatchLowerBound(helper_->Ends()[t], id); watcher->WatchLowerBound(helper_->Sizes()[t], id); - watcher->WatchLowerBound(demands_[t], id); + watcher->WatchLowerBound(demands_->Demands()[t], id); if (!helper_->IsPresent(t) && !helper_->IsAbsent(t)) { watcher->WatchLiteral(helper_->PresenceLiteral(t), id); } diff --git a/ortools/sat/cumulative_energy.h b/ortools/sat/cumulative_energy.h index 5276de343a..2050e570bb 100644 --- a/ortools/sat/cumulative_energy.h +++ b/ortools/sat/cumulative_energy.h @@ -29,6 +29,9 @@ namespace sat { // Enforces the existence of a preemptive schedule where every task is executed // inside its interval, using energy units of the resource during execution. // +// Important: This only uses the energies min/max and not the actual demand +// of a task. It can thus be used in some non-conventional situation. +// // All energy expression are assumed to take a non-negative value; // if the energy of a task is 0, the task can run anywhere. // The schedule never uses more than capacity units of energy at a given time. @@ -36,35 +39,27 @@ namespace sat { // This is mathematically equivalent to making a model with energy(task) // different tasks with demand and size 1, but is much more efficient, // since it uses O(|tasks|) variables instead of O(sum_{task} |energy(task)|). -void AddCumulativeEnergyConstraint(std::vector energies, - AffineExpression capacity, - SchedulingConstraintHelper* helper, - Model* model); - -// Creates a CumulativeEnergyConstraint where the energy of each interval is -// the product of the demands times its size. -// -// TODO(user): This is not ideal if both size and demands are variable. -void AddCumulativeOverloadChecker(const std::vector& demands, - AffineExpression capacity, +void AddCumulativeOverloadChecker(AffineExpression capacity, SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, Model* model); -// Implementation of AddCumulativeEnergyConstraint(). +// Implementation of AddCumulativeOverloadChecker(). class CumulativeEnergyConstraint : public PropagatorInterface { public: - CumulativeEnergyConstraint(std::vector energies, - AffineExpression capacity, - IntegerTrail* integer_trail, - SchedulingConstraintHelper* helper); + CumulativeEnergyConstraint(AffineExpression capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, Model* model); + bool Propagate() final; void RegisterWith(GenericLiteralWatcher* watcher); private: - const std::vector energies_; const AffineExpression capacity_; IntegerTrail* integer_trail_; SchedulingConstraintHelper* helper_; + SchedulingDemandHelper* demands_; + ThetaLambdaTree theta_tree_; // Task characteristics. @@ -85,10 +80,10 @@ class CumulativeIsAfterSubsetConstraint : public PropagatorInterface { public: CumulativeIsAfterSubsetConstraint(IntegerVariable var, IntegerValue offset, AffineExpression capacity, - const std::vector demands, const std::vector subtasks, - IntegerTrail* integer_trail, - SchedulingConstraintHelper* helper); + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, + Model* model); bool Propagate() final; void RegisterWith(GenericLiteralWatcher* watcher); @@ -97,7 +92,6 @@ class CumulativeIsAfterSubsetConstraint : public PropagatorInterface { const IntegerVariable var_to_push_; const IntegerValue offset_; const AffineExpression capacity_; - const std::vector demands_; const std::vector subtasks_; // Computed at construction time, this is const. @@ -109,6 +103,7 @@ class CumulativeIsAfterSubsetConstraint : public PropagatorInterface { IntegerTrail* integer_trail_; SchedulingConstraintHelper* helper_; + SchedulingDemandHelper* demands_; }; } // namespace sat diff --git a/ortools/sat/diffn.cc b/ortools/sat/diffn.cc index 326aadd307..5a7a50936b 100644 --- a/ortools/sat/diffn.cc +++ b/ortools/sat/diffn.cc @@ -112,7 +112,6 @@ void AddDiffnCumulativeRelationOnX(SchedulingConstraintHelper* x, WeightedSumGreaterOrEqual({capacity.var, min_start_var, max_end_var}, coeffs, capacity.constant.value())); - auto* integer_trail = model->GetOrCreate(); auto* watcher = model->GetOrCreate(); const SatParameters* params = model->GetOrCreate(); @@ -121,12 +120,19 @@ void AddDiffnCumulativeRelationOnX(SchedulingConstraintHelper* x, bool add_energetic_relaxation = params->use_energetic_reasoning_in_no_overlap_2d(); + // Needed if we use one of the relaxation below. + SchedulingDemandHelper* demands; + if (add_timetabling_relaxation || add_energetic_relaxation) { + demands = model->TakeOwnership(new SchedulingDemandHelper(sizes, x, model)); + } + // Propagator responsible for applying Timetabling filtering rule. It // increases the minimum of the start variables, decrease the maximum of the // end variables, and increase the minimum of the capacity variable. if (add_timetabling_relaxation) { + DCHECK(demands != nullptr); TimeTablingPerTask* time_tabling = - new TimeTablingPerTask(sizes, capacity, integer_trail, x); + new TimeTablingPerTask(capacity, x, demands, model); time_tabling->RegisterWith(watcher); model->TakeOwnership(time_tabling); } @@ -134,7 +140,8 @@ void AddDiffnCumulativeRelationOnX(SchedulingConstraintHelper* x, // Propagator responsible for applying the Overload Checking filtering rule. // It increases the minimum of the capacity variable. if (add_energetic_relaxation) { - AddCumulativeOverloadChecker(sizes, capacity, x, model); + DCHECK(demands != nullptr); + AddCumulativeOverloadChecker(capacity, x, demands, model); } } diff --git a/ortools/sat/disjunctive.cc b/ortools/sat/disjunctive.cc index ee5e46079f..78ed91f6b2 100644 --- a/ortools/sat/disjunctive.cc +++ b/ortools/sat/disjunctive.cc @@ -73,8 +73,11 @@ std::function Disjunctive( if (/*DISABLES_CODE*/ (false)) { const AffineExpression one(IntegerValue(1)); std::vector demands(vars.size(), one); - TimeTablingPerTask* timetable = new TimeTablingPerTask( - demands, one, model->GetOrCreate(), helper); + SchedulingDemandHelper* demands_helper = model->TakeOwnership( + new SchedulingDemandHelper(demands, helper, model)); + + TimeTablingPerTask* timetable = + new TimeTablingPerTask(one, helper, demands_helper, model); timetable->RegisterWith(watcher); model->TakeOwnership(timetable); return; diff --git a/ortools/sat/implied_bounds.cc b/ortools/sat/implied_bounds.cc index 2ac1310cb7..4a9918edb8 100644 --- a/ortools/sat/implied_bounds.cc +++ b/ortools/sat/implied_bounds.cc @@ -23,7 +23,6 @@ #include #include "absl/container/flat_hash_map.h" -#include "absl/meta/type_traits.h" #include "absl/strings/str_cat.h" #include "ortools/base/logging.h" #include "ortools/base/strong_vector.h" @@ -282,22 +281,26 @@ std::string EncodingStr(const std::vector& enc) { // bounds repository. Because if we can reconcile an encoding, then any of the // literal in the at most one should imply a value on the boolean view use in // the size2 affine. -bool TryToReconcileEncodings( +std::vector TryToReconcileEncodings( const AffineExpression& size2_affine, const AffineExpression& affine, - const std::vector& affine_var_encoding, Model* model, - LinearConstraintBuilder* builder) { + const std::vector& affine_var_encoding, + bool put_affine_left_in_result, Model* model) { IntegerEncoder* integer_encoder = model->GetOrCreate(); IntegerVariable binary = size2_affine.var; - if (!integer_encoder->VariableIsFullyEncoded(binary)) return false; + std::vector terms; + if (!integer_encoder->VariableIsFullyEncoded(binary)) return terms; const std::vector& size2_enc = integer_encoder->FullDomainEncoding(binary); - CHECK_EQ(2, size2_enc.size()); + + // TODO(user): I am not sure how this can happen since size2_affine is + // supposed to be non-fixed. Maybe we miss some propag. Investigate. + if (size2_enc.size() != 2) return terms; + Literal lit0 = size2_enc[0].literal; - IntegerValue value0 = - size2_enc[0].value * size2_affine.coeff + size2_affine.constant; + IntegerValue value0 = size2_affine.ValueAt(size2_enc[0].value); Literal lit1 = size2_enc[1].literal; - IntegerValue value1 = - size2_enc[1].value * size2_affine.coeff + size2_affine.constant; + IntegerValue value1 = size2_affine.ValueAt(size2_enc[1].value); + for (const auto& [unused, candidate_literal] : affine_var_encoding) { if (candidate_literal == lit1) { std::swap(lit0, lit1); @@ -305,32 +308,93 @@ bool TryToReconcileEncodings( } if (candidate_literal != lit0) continue; - // Compute the minimum energy. - IntegerValue min_energy = kMaxIntegerValue; + // Build the decomposition. for (const auto& [value, literal] : affine_var_encoding) { - const IntegerValue energy = literal == lit0 - ? value0 * affine.ValueAt(value) - : value1 * affine.ValueAt(value); - min_energy = std::min(energy, min_energy); + const IntegerValue size_2_value = literal == lit0 ? value0 : value1; + const IntegerValue affine_value = affine.ValueAt(value); + if (put_affine_left_in_result) { + terms.push_back({literal, affine_value, size_2_value}); + } else { + terms.push_back({literal, size_2_value, affine_value}); + } } + break; + } - // Build the energy expression. - builder->Clear(); - builder->AddConstant(min_energy); - for (const auto& [value, literal] : affine_var_encoding) { - const IntegerValue energy = literal == lit0 - ? value0 * affine.ValueAt(value) - : value1 * affine.ValueAt(value); - if (energy > min_energy) { - if (!builder->AddLiteralTerm(literal, energy - min_energy)) { - return false; + return terms; +} + +std::vector TryToDecomposeProduct( + const AffineExpression& left, const AffineExpression& right, Model* model) { + IntegerTrail* integer_trail = model->GetOrCreate(); + if (integer_trail->IsFixed(left) || integer_trail->IsFixed(right)) return {}; + + // Fill in the encodings for the left variable. + ImpliedBounds* implied_bounds = model->GetOrCreate(); + const absl::flat_hash_map>& + left_encodings = implied_bounds->GetElementEncodings(left.var); + + // Fill in the encodings for the right variable. + const absl::flat_hash_map>& + right_encodings = implied_bounds->GetElementEncodings(right.var); + + std::vector compatible_keys; + for (const auto& [index, encoding] : left_encodings) { + if (right_encodings.contains(index)) { + compatible_keys.push_back(index); + } + } + + if (compatible_keys.empty()) { + if (integer_trail->InitialVariableDomain(left.var).Size() == 2) { + for (const auto& [index, right_encoding] : right_encodings) { + const std::vector result = + TryToReconcileEncodings(left, right, right_encoding, + /*put_affine_left_in_result=*/false, model); + if (!result.empty()) { + return result; } } } - return true; + if (integer_trail->InitialVariableDomain(right.var).Size() == 2) { + for (const auto& [index, left_encoding] : left_encodings) { + const std::vector result = + TryToReconcileEncodings(right, left, left_encoding, + /*put_affine_left_in_result=*/true, model); + if (!result.empty()) { + return result; + } + } + } + return {}; } - return false; + if (compatible_keys.size() > 1) { + VLOG(1) << "More than one exactly_one involved in the encoding of the two " + "variables"; + } + + // Select the compatible encoding with the minimum index. + const int min_index = + *std::min_element(compatible_keys.begin(), compatible_keys.end()); + // By construction, encodings follow the order of literals in the exactly_one + // constraint. + const std::vector& left_encoding = + left_encodings.at(min_index); + const std::vector& right_encoding = + right_encodings.at(min_index); + DCHECK_EQ(left_encoding.size(), right_encoding.size()); + + // Build decomposition of the product. + std::vector terms; + for (int i = 0; i < left_encoding.size(); ++i) { + const Literal literal = left_encoding[i].literal; + DCHECK_EQ(literal, right_encoding[i].literal); + terms.push_back({literal, left.ValueAt(left_encoding[i].value), + right.ValueAt(right_encoding[i].value)}); + } + + return terms; } // TODO(user): Experiment with x * x where constants = 0, x is @@ -373,94 +437,27 @@ bool DetectLinearEncodingOfProducts(const AffineExpression& left, return true; } - // Fill in the encodings for the left variable. - ImpliedBounds* implied_bounds = model->GetOrCreate(); - const absl::flat_hash_map>& - left_encodings = implied_bounds->GetElementEncodings(left.var); + const std::vector product = + TryToDecomposeProduct(left, right, model); + if (product.empty()) return false; - // Fill in the encodings for the right variable. - const absl::flat_hash_map>& - right_encodings = implied_bounds->GetElementEncodings(right.var); - - std::vector compatible_keys; - for (const auto& [index, encoding] : left_encodings) { - if (right_encodings.contains(index)) { - compatible_keys.push_back(index); - } + IntegerValue min_coefficient = kMaxIntegerValue; + for (const LiteralValueValue& term : product) { + min_coefficient = + std::min(min_coefficient, term.left_value * term.right_value); } - if (compatible_keys.empty()) { - if (integer_trail->InitialVariableDomain(left.var).Size() == 2) { - for (const auto& [index, right_encoding] : right_encodings) { - if (TryToReconcileEncodings(left, right, right_encoding, model, - builder)) { - return true; - } - } - } - if (integer_trail->InitialVariableDomain(right.var).Size() == 2) { - for (const auto& [index, left_encoding] : left_encodings) { - if (TryToReconcileEncodings(right, left, left_encoding, model, - builder)) { - return true; - } - } - } - return false; - } - - if (compatible_keys.size() > 1) { - VLOG(1) << "More than one exactly_one involved in the encoding of the two " - "variables"; - } - - // Select the compatible encoding with the minimum index. - const int min_index = - *std::min_element(compatible_keys.begin(), compatible_keys.end()); - // By construction, encodings follow the order of literals in the exactly_one - // constraint. - const std::vector& left_encoding = - left_encodings.at(min_index); - const std::vector& right_encoding = - right_encodings.at(min_index); - DCHECK_EQ(left_encoding.size(), right_encoding.size()); - - // Compute the min energy. - IntegerValue min_energy = kMaxIntegerValue; - for (int i = 0; i < left_encoding.size(); ++i) { - const IntegerValue energy = left.ValueAt(left_encoding[i].value) * - right.ValueAt(right_encoding[i].value); - min_energy = std::min(min_energy, energy); - } - - // Build the linear formulation of the energy. - for (int i = 0; i < left_encoding.size(); ++i) { - const IntegerValue energy = left.ValueAt(left_encoding[i].value) * - right.ValueAt(right_encoding[i].value); - if (energy == min_energy) continue; - DCHECK_GT(energy, min_energy); - const Literal lit = left_encoding[i].literal; - DCHECK_EQ(lit, right_encoding[i].literal); - - if (!builder->AddLiteralTerm(lit, energy - min_energy)) { + for (const LiteralValueValue& term : product) { + const IntegerValue coefficient = + term.left_value * term.right_value - min_coefficient; + if (coefficient == 0) continue; + if (!builder->AddLiteralTerm(term.literal, coefficient)) { return false; } } - builder->AddConstant(min_energy); + builder->AddConstant(min_coefficient); return true; } -std::optional TryToLinearizeProduct( - const AffineExpression& left, const AffineExpression& right, Model* model) { - LinearConstraintBuilder builder(model); - if (DetectLinearEncodingOfProducts(left, right, model, &builder)) { - // The expression must only have positive coefficient because we will call - // Min() on it and that function expect it this way. - return CanonicalizeExpr(builder.BuildExpression()); - } else { - return std::optional(); - } -} - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/implied_bounds.h b/ortools/sat/implied_bounds.h index a25220b46d..49ebd3ca63 100644 --- a/ortools/sat/implied_bounds.h +++ b/ortools/sat/implied_bounds.h @@ -205,6 +205,14 @@ class ImpliedBounds { int64_t num_enqueued_in_var_to_bounds_ = 0; }; +// Tries to decompose a product left * right in a list of constant alternative +// left_value * right_value controlled by literals in an exactly one +// relationship. We construct this by using literals from the full encoding or +// element encodings of the variables of the two affine expressions. +// If it fails, it returns an empty vector. +std::vector TryToDecomposeProduct( + const AffineExpression& left, const AffineExpression& right, Model* model); + // Looks at value encodings and detects if the product of two variables can be // linearized. // @@ -217,11 +225,6 @@ bool DetectLinearEncodingOfProducts(const AffineExpression& left, const AffineExpression& right, Model* model, LinearConstraintBuilder* builder); -// Try to linearize left * right and returns the result. If we cannot linearize -// the result will have no value. -std::optional TryToLinearizeProduct( - const AffineExpression& left, const AffineExpression& right, Model* model); - } // namespace sat } // namespace operations_research diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index 778a7b7353..018e9ca38c 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -508,6 +508,26 @@ LiteralIndex IntegerEncoder::SearchForLiteralAtOrBefore( return after_it->second.Index(); } +ABSL_MUST_USE_RESULT bool IntegerEncoder::LiteralOrNegationHasView( + Literal lit, IntegerVariable* view, bool* view_is_direct) const { + const IntegerVariable direct_var = GetLiteralView(lit); + const IntegerVariable opposite_var = GetLiteralView(lit.Negated()); + // If a literal has both views, we want to always keep the same + // representative: the smallest IntegerVariable. + if (direct_var != kNoIntegerVariable && + (opposite_var == kNoIntegerVariable || direct_var <= opposite_var)) { + if (view != nullptr) *view = direct_var; + if (view_is_direct != nullptr) *view_is_direct = true; + return true; + } + if (opposite_var != kNoIntegerVariable) { + if (view != nullptr) *view = opposite_var; + if (view_is_direct != nullptr) *view_is_direct = false; + return true; + } + return false; +} + IntegerTrail::~IntegerTrail() { if (parameters_.log_search_progress() && num_decisions_to_break_loop_ > 0) { VLOG(1) << "Num decisions to break propagation loop: " diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index c61c29b3f7..2054138dff 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -338,6 +338,23 @@ struct ValueLiteralPair { std::ostream& operator<<(std::ostream& os, const ValueLiteralPair& p); +struct LiteralValueValue { + Literal literal; + IntegerValue left_value; + IntegerValue right_value; + + // Used for testing. + bool operator==(const LiteralValueValue& rhs) const { + return literal.Index() == rhs.literal.Index() && + left_value == rhs.left_value && right_value == rhs.right_value; + } + + std::string DebugString() const { + return absl::StrCat("(lit(", literal.Index().value(), ") * ", + left_value.value(), " * ", right_value.value(), ")"); + } +}; + // Each integer variable x will be associated with a set of literals encoding // (x >= v) for some values of v. This class maintains the relationship between // the integer variables and such literals which can be created by a call to @@ -500,10 +517,9 @@ class IntegerEncoder { // If this is true, then a literal can be linearized with an affine expression // involving an integer variable. - const bool LiteralOrNegationHasView(Literal lit) const { - return GetLiteralView(lit) != kNoIntegerVariable || - GetLiteralView(lit.Negated()) != kNoIntegerVariable; - } + ABSL_MUST_USE_RESULT bool LiteralOrNegationHasView( + Literal lit, IntegerVariable* view = nullptr, + bool* view_is_direct = nullptr) const; // Returns a Boolean literal associated with a bound lower than or equal to // the one of the given IntegerLiteral. If the given IntegerLiteral is true, diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index 7122d91da4..c74fc44293 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -1045,6 +1045,10 @@ SatSolver::Status SolveIntegerProblem(Model* model) { return sat_solver->UnsatStatus(); } + // In multi-thread, we really only want to save the LP relaxation for thread + // with high linearization level to avoid to pollute the repository with + // sub-par lp solutions. + // // TODO(user): Experiment more around dynamically changing the // threshold for storing LP solutions in the pool. Alternatively expose // this as parameter so this can be tuned later. @@ -1053,7 +1057,8 @@ SatSolver::Status SolveIntegerProblem(Model* model) { // change. Avoid adding solution that are too deep in the tree (most // variable fixed). Also use a callback rather than having this here, we // don't want this file to depend on cp_model.proto. - if (model->Get() != nullptr) { + if (model->Get() != nullptr && + sat_parameters.linearization_level() >= 2) { num_decisions_since_last_lp_record_++; if (num_decisions_since_last_lp_record_ >= 100) { // NOTE: We can actually record LP solutions more frequently. However @@ -1332,7 +1337,7 @@ void ContinuousProber::LogStatistics() { num_bounds_tried_, " #new_integer_bounds:", shared_bounds_manager_->NumBoundsExported("probing"), ", #new_binary_clauses:", prober_->num_new_binary_clauses()), - &last_logging_time_); + parameters_.log_frequency_in_seconds(), &last_logging_time_); } } // namespace sat diff --git a/ortools/sat/intervals.cc b/ortools/sat/intervals.cc index efeda1de6e..bd16c20f27 100644 --- a/ortools/sat/intervals.cc +++ b/ortools/sat/intervals.cc @@ -22,6 +22,7 @@ #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/strong_vector.h" +#include "ortools/sat/implied_bounds.h" #include "ortools/sat/integer.h" #include "ortools/sat/integer_expr.h" #include "ortools/sat/linear_constraint.h" @@ -519,22 +520,33 @@ bool SchedulingConstraintHelper::PushIntervalBound(int t, IntegerLiteral lit) { return true; } -bool SchedulingConstraintHelper::IncreaseStartMin(int t, - IntegerValue new_start_min) { +bool SchedulingConstraintHelper::IncreaseStartMin(int t, IntegerValue value) { if (starts_[t].var == kNoIntegerVariable) { - if (new_start_min > starts_[t].constant) return PushTaskAbsence(t); + if (value > starts_[t].constant) return PushTaskAbsence(t); return true; } - return PushIntervalBound(t, starts_[t].GreaterOrEqual(new_start_min)); + return PushIntervalBound(t, starts_[t].GreaterOrEqual(value)); } -bool SchedulingConstraintHelper::DecreaseEndMax(int t, - IntegerValue new_end_max) { +bool SchedulingConstraintHelper::IncreaseEndMin(int t, IntegerValue value) { if (ends_[t].var == kNoIntegerVariable) { - if (new_end_max < ends_[t].constant) return PushTaskAbsence(t); + if (value > ends_[t].constant) return PushTaskAbsence(t); return true; } - return PushIntervalBound(t, ends_[t].LowerOrEqual(new_end_max)); + return PushIntervalBound(t, ends_[t].GreaterOrEqual(value)); +} + +bool SchedulingConstraintHelper::DecreaseEndMax(int t, IntegerValue value) { + if (ends_[t].var == kNoIntegerVariable) { + if (value < ends_[t].constant) return PushTaskAbsence(t); + return true; + } + return PushIntervalBound(t, ends_[t].LowerOrEqual(value)); +} + +bool SchedulingConstraintHelper::PushLiteral(Literal l) { + integer_trail_->EnqueueLiteral(l, literal_reason_, integer_reason_); + return true; } bool SchedulingConstraintHelper::PushTaskAbsence(int t) { @@ -632,5 +644,282 @@ IntegerValue SchedulingConstraintHelper::GetMinOverlap(int t, std::min(EndMin(t) - start, end - StartMax(t))); } +SchedulingDemandHelper::SchedulingDemandHelper( + std::vector demands, SchedulingConstraintHelper* helper, + Model* model) + : integer_trail_(model->GetOrCreate()), + assignment_(model->GetOrCreate()->Assignment()), + demands_(std::move(demands)), + helper_(helper) { + const int num_tasks = helper->NumTasks(); + linearized_energies_.resize(num_tasks); + decomposed_energies_.resize(num_tasks); + cached_energies_min_.resize(num_tasks, kMinIntegerValue); + cached_energies_max_.resize(num_tasks, kMaxIntegerValue); + energy_is_quadratic_.resize(num_tasks, false); + + // For the special case were demands is empty. + if (demands_.size() != num_tasks) return; + for (int t = 0; t < num_tasks; ++t) { + const AffineExpression size = helper->Sizes()[t]; + const AffineExpression demand = demands_[t]; + decomposed_energies_[t] = TryToDecomposeProduct(size, demand, model); + } +} + +IntegerValue SchedulingDemandHelper::SimpleEnergyMin(int t) const { + if (demands_.empty()) return kMinIntegerValue; + return DemandMin(t) * helper_->SizeMin(t); +} + +IntegerValue SchedulingDemandHelper::LinearEnergyMin(int t) const { + if (!linearized_energies_[t].has_value()) return kMinIntegerValue; + return linearized_energies_[t]->Min(*integer_trail_); +} + +IntegerValue SchedulingDemandHelper::DecomposedEnergyMin(int t) const { + if (decomposed_energies_[t].empty()) return kMinIntegerValue; + IntegerValue result = kMaxIntegerValue; + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energies_[t]) { + if (assignment_.LiteralIsTrue(lit)) { + return fixed_size * fixed_demand; + } + if (assignment_.LiteralIsFalse(lit)) continue; + result = std::min(result, fixed_size * fixed_demand); + } + DCHECK_NE(result, kMaxIntegerValue); + return result; +} + +IntegerValue SchedulingDemandHelper::SimpleEnergyMax(int t) const { + if (demands_.empty()) return kMaxIntegerValue; + return DemandMax(t) * helper_->SizeMax(t); +} + +IntegerValue SchedulingDemandHelper::LinearEnergyMax(int t) const { + if (!linearized_energies_[t].has_value()) return kMaxIntegerValue; + return linearized_energies_[t]->Max(*integer_trail_); +} + +IntegerValue SchedulingDemandHelper::DecomposedEnergyMax(int t) const { + if (decomposed_energies_[t].empty()) return kMaxIntegerValue; + IntegerValue result = kMinIntegerValue; + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energies_[t]) { + if (assignment_.LiteralIsTrue(lit)) { + return fixed_size * fixed_demand; + } + if (assignment_.LiteralIsFalse(lit)) continue; + result = std::max(result, fixed_size * fixed_demand); + } + DCHECK_NE(result, kMinIntegerValue); + return result; +} + +void SchedulingDemandHelper::CacheAllEnergyValues() { + const int num_tasks = cached_energies_min_.size(); + for (int t = 0; t < num_tasks; ++t) { + cached_energies_min_[t] = std::max( + {SimpleEnergyMin(t), LinearEnergyMin(t), DecomposedEnergyMin(t)}); + CHECK_NE(cached_energies_min_[t], kMinIntegerValue); + energy_is_quadratic_[t] = + decomposed_energies_[t].empty() && !demands_.empty() && + !integer_trail_->IsFixed(demands_[t]) && !helper_->SizeIsFixed(t); + cached_energies_max_[t] = std::min( + {SimpleEnergyMax(t), LinearEnergyMax(t), DecomposedEnergyMax(t)}); + CHECK_NE(cached_energies_min_[t], kMaxIntegerValue); + } +} + +IntegerValue SchedulingDemandHelper::DemandMin(int t) const { + DCHECK_LT(t, demands_.size()); + return integer_trail_->LowerBound(demands_[t]); +} + +IntegerValue SchedulingDemandHelper::DemandMax(int t) const { + DCHECK_LT(t, demands_.size()); + return integer_trail_->UpperBound(demands_[t]); +} + +bool SchedulingDemandHelper::DemandIsFixed(int t) const { + return integer_trail_->IsFixed(demands_[t]); +} + +bool SchedulingDemandHelper::DecreaseEnergyMax(int t, IntegerValue value) { + if (value < EnergyMin(t)) { + if (helper_->IsOptional(t)) { + return helper_->PushTaskAbsence(t); + } else { + return helper_->ReportConflict(); + } + } else if (!decomposed_energies_[t].empty()) { + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energies_[t]) { + if (fixed_size * fixed_demand > value) { + if (assignment_.LiteralIsTrue(lit)) return helper_->ReportConflict(); + if (assignment_.LiteralIsFalse(lit)) continue; + if (!helper_->PushLiteral(lit.Negated())) return false; + } + } + } else if (linearized_energies_[t].has_value() && + linearized_energies_[t]->vars.size() == 1) { + const LinearExpression& e = linearized_energies_[t].value(); + const AffineExpression affine_energy(e.vars[0], e.coeffs[0], e.offset); + const IntegerLiteral deduction = affine_energy.LowerOrEqual(value); + if (!helper_->PushIntegerLiteralIfTaskPresent(t, deduction)) { + return false; + } + } else { + // TODO(user): Propagate if possible. + VLOG(3) << "Cumulative energy missed propagation"; + } + return true; +} + +void SchedulingDemandHelper::AddDemandMinReason(int t) { + DCHECK_LT(t, demands_.size()); + if (demands_[t].var != kNoIntegerVariable) { + helper_->MutableIntegerReason()->push_back( + integer_trail_->LowerBoundAsLiteral(demands_[t].var)); + } +} + +void SchedulingDemandHelper::AddEnergyMinReason(int t) { + // We prefer these reason in order. + const IntegerValue value = cached_energies_min_[t]; + if (DecomposedEnergyMin(t) >= value) { + auto* reason = helper_->MutableLiteralReason(); + const int old_size = reason->size(); + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energies_[t]) { + if (assignment_.LiteralIsTrue(lit)) { + reason->resize(old_size); + reason->push_back(lit.Negated()); + return; + } else if (fixed_size * fixed_demand < value && + assignment_.LiteralIsFalse(lit)) { + reason->push_back(lit); + } + } + } else if (SimpleEnergyMin(t) >= value) { + AddDemandMinReason(t); + helper_->AddSizeMinReason(t); + } else { + DCHECK_GE(LinearEnergyMin(t), value); + for (const IntegerVariable var : linearized_energies_[t]->vars) { + helper_->MutableIntegerReason()->push_back( + integer_trail_->LowerBoundAsLiteral(var)); + } + } +} + +void SchedulingDemandHelper::OverrideLinearizedEnergies( + const std::vector& energies) { + const int num_tasks = energies.size(); + DCHECK_EQ(num_tasks, helper_->NumTasks()); + linearized_energies_.resize(num_tasks); + for (int t = 0; t < num_tasks; ++t) { + linearized_energies_[t] = energies[t]; + if (DEBUG_MODE) { + for (const IntegerValue coeff : linearized_energies_[t]->coeffs) { + DCHECK_GE(coeff, 0); + } + } + } +} + +void SchedulingDemandHelper::OverrideDecomposedEnergies( + const std::vector>& energies) { + DCHECK_EQ(energies.size(), helper_->NumTasks()); + decomposed_energies_ = energies; +} + +IntegerValue SchedulingDemandHelper::EnergyMinInWindow( + int t, IntegerValue window_start, IntegerValue window_end) { + if (window_end <= window_start) return IntegerValue(0); + + // Returns zero if the interval do not necessarily overlap. + const IntegerValue start_max = helper_->StartMax(t); + const IntegerValue end_min = helper_->EndMin(t); + if (end_min <= window_start) return IntegerValue(0); + if (start_max >= window_end) return IntegerValue(0); + const IntegerValue window_size = window_end - window_start; + const IntegerValue simple_energy_min = + DemandMin(t) * std::min({end_min - window_start, window_end - start_max, + helper_->SizeMin(t), window_size}); + if (decomposed_energies_[t].empty()) return simple_energy_min; + + IntegerValue result = kMaxIntegerValue; + const IntegerValue start_min = helper_->StartMin(t); + const IntegerValue end_max = helper_->EndMax(t); + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energies_[t]) { + if (assignment_.LiteralIsTrue(lit)) { + // Both should be identical, so we don't recompute it. + return simple_energy_min; + } + if (assignment_.LiteralIsFalse(lit)) continue; + const IntegerValue alt_em = std::max(end_min, start_min + fixed_size); + const IntegerValue alt_sm = std::min(start_max, end_max - fixed_size); + const IntegerValue energy_min = + fixed_demand * std::min({alt_em - window_start, window_end - alt_sm, + fixed_size, window_size}); + result = std::min(result, energy_min); + } + if (result == kMaxIntegerValue) return simple_energy_min; + return std::max(simple_energy_min, result); +} + +// Since we usually ask way less often for the reason, we redo the computation +// here. +void SchedulingDemandHelper::AddEnergyMinInWindowReason( + int t, IntegerValue window_start, IntegerValue window_end) { + const IntegerValue actual_energy_min = + EnergyMinInWindow(t, window_start, window_end); + if (actual_energy_min == 0) return; + + // Return simple reason right away if there is no decomposition or the simple + // energy is enough. + const IntegerValue start_max = helper_->StartMax(t); + const IntegerValue end_min = helper_->EndMin(t); + const IntegerValue simple_energy_min = + DemandMin(t) * std::min({end_min - window_start, window_end - start_max, + helper_->SizeMin(t)}); + if (simple_energy_min == actual_energy_min) { + AddDemandMinReason(t); + helper_->AddSizeMinReason(t); + helper_->AddStartMaxReason(t, start_max); + helper_->AddEndMinReason(t, end_min); + return; + } + + // TODO(user): only include the one we need? + const IntegerValue start_min = helper_->StartMin(t); + const IntegerValue end_max = helper_->EndMax(t); + DCHECK(!decomposed_energies_[t].empty()); + helper_->AddStartMinReason(t, start_min); + helper_->AddStartMaxReason(t, start_max); + helper_->AddEndMinReason(t, end_min); + helper_->AddEndMaxReason(t, end_max); + + auto* literal_reason = helper_->MutableLiteralReason(); + const int old_size = literal_reason->size(); + + DCHECK(!decomposed_energies_[t].empty()); + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energies_[t]) { + // Should be the same in most cases. + if (assignment_.LiteralIsTrue(lit)) { + literal_reason->resize(old_size); + literal_reason->push_back(lit.Negated()); + return; + } + if (assignment_.LiteralIsFalse(lit)) { + const IntegerValue alt_em = std::max(end_min, start_min + fixed_size); + const IntegerValue alt_sm = std::min(start_max, end_max - fixed_size); + const IntegerValue energy_min = + fixed_demand * + std::min({alt_em - window_start, window_end - alt_sm, fixed_size}); + if (energy_min >= actual_energy_min) continue; + literal_reason->push_back(lit); + } + } +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/intervals.h b/ortools/sat/intervals.h index cefd47b477..614dddfeb4 100644 --- a/ortools/sat/intervals.h +++ b/ortools/sat/intervals.h @@ -358,8 +358,10 @@ class SchedulingConstraintHelper : public PropagatorInterface, // conditionned on its presence. The functions will do the correct thing // depending on whether or not the start_min/end_max are optional variables // whose presence implies the interval presence. - ABSL_MUST_USE_RESULT bool IncreaseStartMin(int t, IntegerValue new_start_min); - ABSL_MUST_USE_RESULT bool DecreaseEndMax(int t, IntegerValue new_end_max); + ABSL_MUST_USE_RESULT bool IncreaseStartMin(int t, IntegerValue value); + ABSL_MUST_USE_RESULT bool IncreaseEndMin(int t, IntegerValue value); + ABSL_MUST_USE_RESULT bool DecreaseEndMax(int t, IntegerValue value); + ABSL_MUST_USE_RESULT bool PushLiteral(Literal l); ABSL_MUST_USE_RESULT bool PushTaskAbsence(int t); ABSL_MUST_USE_RESULT bool PushTaskPresence(int t); ABSL_MUST_USE_RESULT bool PushIntegerLiteral(IntegerLiteral lit); @@ -494,6 +496,100 @@ class SchedulingConstraintHelper : public PropagatorInterface, std::vector already_added_to_other_reasons_; }; +// Helper class for cumulative constraint to wrap demands and expose concept +// like energy. +// +// In a cumulative constraint, an interval always has a size and a demand, but +// it can also have a set of "selector" literals each associated with a fixed +// size / fixed demands. This allows more precise energy estimation. +// +// TODO(user): Cache energy min and reason for the non O(1) cases. +class SchedulingDemandHelper { + public: + // Hack: this can be called with and empty demand vector as long as + // OverrideEnergies() is called to define the energies. + SchedulingDemandHelper(std::vector demands, + SchedulingConstraintHelper* helper, Model* model); + + // When defined, the interval will consume this much demand during its whole + // duration. Some propagator only relies on the "energy" and thus never uses + // this. + IntegerValue DemandMin(int t) const; + IntegerValue DemandMax(int t) const; + bool DemandIsFixed(int t) const; + void AddDemandMinReason(int t); + const std::vector& Demands() const { return demands_; } + + // The "energy" is usually size * demand, but in some non-conventional usage + // it might have a more complex formula. In all case, the energy is assumed + // to be only consumed during the interval duration. + // + // IMPORTANT: One must call CacheAllEnergyValues() for the values to be + // updated. TODO(user): this is error prone, maybe we should revisit. But if + // there is many alternatives, we don't want to rescan the list more than a + // linear number of time per propagation. + // + // TODO(user): Add more complex EnergyMinBefore(time) once we also support + // expressing the interval as a set of alternatives. + void CacheAllEnergyValues(); + IntegerValue EnergyMin(int t) const { return cached_energies_min_[t]; } + IntegerValue EnergyMax(int t) const { return cached_energies_max_[t]; } + bool EnergyIsQuadratic(int t) const { return energy_is_quadratic_[t]; } + void AddEnergyMinReason(int t); + + // Returns the energy min in [start, end]. + // + // Note(user): These functions are not in O(1) if the decomposition is used, + // so we have to be careful in not calling them too often. + IntegerValue EnergyMinInWindow(int t, IntegerValue window_start, + IntegerValue window_end); + void AddEnergyMinInWindowReason(int t, IntegerValue window_start, + IntegerValue window_end); + + // Important: This might not do anything depending on the representation of + // the energy we have. + ABSL_MUST_USE_RESULT bool DecreaseEnergyMax(int t, IntegerValue value); + + // Different optional representation of the energy of an interval. + // Important: first value is size, second value is demand. + const std::vector>& DecomposedEnergies() + const { + return decomposed_energies_; + } + + // Visible for testing. + void OverrideLinearizedEnergies( + const std::vector& energies); + void OverrideDecomposedEnergies( + const std::vector>& energies); + + private: + IntegerValue SimpleEnergyMin(int t) const; + IntegerValue LinearEnergyMin(int t) const; + IntegerValue SimpleEnergyMax(int t) const; + IntegerValue LinearEnergyMax(int t) const; + IntegerValue DecomposedEnergyMin(int t) const; + IntegerValue DecomposedEnergyMax(int t) const; + + IntegerTrail* integer_trail_; + const VariablesAssignment& assignment_; + std::vector demands_; + SchedulingConstraintHelper* helper_; + + // Cached value of the energies, as it can be a bit costly to compute. + std::vector cached_energies_min_; + std::vector cached_energies_max_; + std::vector energy_is_quadratic_; + + // A representation of the energies as a set of alternative. + // If subvector is empty, we don't have this representation. + std::vector> decomposed_energies_; + + // A representation of the energies as a set of linear expression. + // If the optional is not set, we don't have this representation. + std::vector> linearized_energies_; +}; + // ============================================================================= // SchedulingConstraintHelper inlined functions. // ============================================================================= diff --git a/ortools/sat/lb_tree_search.cc b/ortools/sat/lb_tree_search.cc index 94d9709a77..42ab67c639 100644 --- a/ortools/sat/lb_tree_search.cc +++ b/ortools/sat/lb_tree_search.cc @@ -475,8 +475,9 @@ SatSolver::Status LbTreeSearch::Search( if (lb > current_objective_lb_) break; } - shared_response_->LogPeriodicMessage("TreeS", SmallProgressString(), - &last_logging_time_); + shared_response_->LogPeriodicMessage( + "TreeS", SmallProgressString(), + parameters_.log_frequency_in_seconds(), &last_logging_time_); if (n < nodes_.size()) { current_branch_.push_back(n); diff --git a/ortools/sat/linear_constraint.cc b/ortools/sat/linear_constraint.cc index 5d6ed705ed..ca2097066e 100644 --- a/ortools/sat/linear_constraint.cc +++ b/ortools/sat/linear_constraint.cc @@ -80,9 +80,30 @@ void LinearConstraintBuilder::AddLinearExpression(const LinearExpression& expr, offset_ += expr.offset * coeff; } +ABSL_MUST_USE_RESULT bool LinearConstraintBuilder::AddDecomposedProduct( + const std::vector& product) { + if (product.empty()) return true; + + IntegerValue product_min = kMaxIntegerValue; + // TODO(user): Checks the value of literals. + for (const LiteralValueValue& term : product) { + product_min = std::min(product_min, term.left_value * term.right_value); + } + + for (const LiteralValueValue& term : product) { + IntegerValue coeff = term.left_value * term.right_value - product_min; + if (coeff == 0) continue; + if (!AddLiteralTerm(term.literal, coeff)) { + return false; + } + } + AddConstant(product_min); + return true; +} + void LinearConstraintBuilder::AddQuadraticLowerBound( - AffineExpression left, AffineExpression right, - IntegerTrail* integer_trail) { + AffineExpression left, AffineExpression right, IntegerTrail* integer_trail, + bool* is_quadratic) { if (integer_trail->IsFixed(left)) { AddTerm(right, integer_trail->FixedValue(left)); } else if (integer_trail->IsFixed(right)) { @@ -94,6 +115,7 @@ void LinearConstraintBuilder::AddQuadraticLowerBound( AddTerm(right, left_min); // Substract the energy counted twice. AddConstant(-left_min * right_min); + if (is_quadratic != nullptr) *is_quadratic = true; } } @@ -104,31 +126,19 @@ void LinearConstraintBuilder::AddConstant(IntegerValue value) { ABSL_MUST_USE_RESULT bool LinearConstraintBuilder::AddLiteralTerm( Literal lit, IntegerValue coeff) { DCHECK(encoder_ != nullptr); - bool has_direct_view = encoder_->GetLiteralView(lit) != kNoIntegerVariable; - bool has_opposite_view = - encoder_->GetLiteralView(lit.Negated()) != kNoIntegerVariable; + IntegerVariable var = kNoIntegerVariable; + bool view_is_direct = true; + if (!encoder_->LiteralOrNegationHasView(lit, &var, &view_is_direct)) { + return false; + } - // If a literal has both views, we want to always keep the same - // representative: the smallest IntegerVariable. Note that AddTerm() will - // also make sure to use the associated positive variable. - if (has_direct_view && has_opposite_view) { - if (encoder_->GetLiteralView(lit) <= - encoder_->GetLiteralView(lit.Negated())) { - has_opposite_view = false; - } else { - has_direct_view = false; - } - } - if (has_direct_view) { - AddTerm(encoder_->GetLiteralView(lit), coeff); - return true; - } - if (has_opposite_view) { - AddTerm(encoder_->GetLiteralView(lit.Negated()), -coeff); + if (view_is_direct) { + AddTerm(var, coeff); + } else { + AddTerm(var, -coeff); offset_ += coeff; - return true; } - return false; + return true; } LinearConstraint LinearConstraintBuilder::Build() { diff --git a/ortools/sat/linear_constraint.h b/ortools/sat/linear_constraint.h index 325129df05..381199a462 100644 --- a/ortools/sat/linear_constraint.h +++ b/ortools/sat/linear_constraint.h @@ -180,6 +180,14 @@ class LinearConstraintBuilder { void AddLinearExpression(const LinearExpression& expr); void AddLinearExpression(const LinearExpression& expr, IntegerValue coeff); + // Add the corresponding decomposed products (obtained from + // TryToDecomposeProduct). The code assumes all literals to be in an + // exactly_one relation. + // It returns false if one literal does not have an integer view, as it + // actually calls AddLiteralTerm(). + ABSL_MUST_USE_RESULT bool AddDecomposedProduct( + const std::vector& product); + // Add literal * coeff to the constaint. Returns false and do nothing if the // given literal didn't have an integer view. ABSL_MUST_USE_RESULT bool AddLiteralTerm(Literal lit, IntegerValue coeff); @@ -197,7 +205,8 @@ class LinearConstraintBuilder { // expression instead. This would depend on the LP value of the left and // right. void AddQuadraticLowerBound(AffineExpression left, AffineExpression right, - IntegerTrail* integer_trail); + IntegerTrail* integer_trail, + bool* is_quadratic = nullptr); // Clears all added terms and constants. Keeps the original bounds. void Clear() { diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc index 3ed2f102bc..1240818293 100644 --- a/ortools/sat/linear_relaxation.cc +++ b/ortools/sat/linear_relaxation.cc @@ -669,15 +669,12 @@ void AppendNoOverlapRelaxationAndCutGenerator(const ConstraintProto& ct, SchedulingConstraintHelper* helper = repository->GetOrCreateHelper(intervals); if (!helper->SynchronizeAndSetTimeDirection(true)) return; - std::vector> energies; - energies.reserve(helper->NumTasks()); - for (int i = 0; i < helper->NumTasks(); ++i) { - LinearConstraintBuilder e(model); - e.AddTerm(helper->Sizes()[i], one); - energies.push_back(CanonicalizeExpr(e.BuildExpression())); - } - AddCumulativeRelaxation(helper, demands, /*capacity=*/one, energies, model, + SchedulingDemandHelper* demands_helper = + new SchedulingDemandHelper(demands, helper, model); + model->TakeOwnership(demands_helper); + + AddCumulativeRelaxation(/*capacity=*/one, helper, demands_helper, model, relaxation); if (model->GetOrCreate()->linearization_level() > 1) { AddNoOverlapCutGenerator(helper, makespan, model, relaxation); @@ -708,19 +705,15 @@ void AppendCumulativeRelaxationAndCutGenerator(const ConstraintProto& ct, // We try to linearize the energy of each task (size * demand). SchedulingConstraintHelper* helper = repository->GetOrCreateHelper(intervals); if (!helper->SynchronizeAndSetTimeDirection(true)) return; - std::vector> energies; - energies.reserve(helper->NumTasks()); - for (int i = 0; i < helper->NumTasks(); ++i) { - energies.push_back( - TryToLinearizeProduct(demands[i], helper->Sizes()[i], model)); - } + SchedulingDemandHelper* demands_helper = + new SchedulingDemandHelper(demands, helper, model); + model->TakeOwnership(demands_helper); // We can now add the relaxation and the cut generators. - AddCumulativeRelaxation(helper, demands, capacity, energies, model, - relaxation); + AddCumulativeRelaxation(capacity, helper, demands_helper, model, relaxation); if (model->GetOrCreate()->linearization_level() > 1) { - AddCumulativeCutGenerator(helper, demands, capacity, energies, makespan, - model, relaxation); + AddCumulativeCutGenerator(capacity, helper, demands_helper, makespan, model, + relaxation); } } @@ -728,13 +721,12 @@ void AppendCumulativeRelaxationAndCutGenerator(const ConstraintProto& ct, // and add the constraint that the sum of energies of each task must fit in the // capacity * span area. // TODO(user): Exploit the makespan if found. -void AddCumulativeRelaxation( - SchedulingConstraintHelper* helper, - const std::vector& demands, - const AffineExpression& capacity, - const std::vector>& energies, Model* model, - LinearRelaxation* relaxation) { +void AddCumulativeRelaxation(const AffineExpression& capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands_helper, + Model* model, LinearRelaxation* relaxation) { const int num_intervals = helper->NumTasks(); + demands_helper->CacheAllEnergyValues(); std::vector presence_literals; std::vector starts; @@ -763,8 +755,7 @@ void AddCumulativeRelaxation( model->GetOrCreate()->GetTrueLiteral()); } - if (!helper->SizeIsFixed(index) || - (!demands.empty() && !integer_trail->IsFixed(demands[index]))) { + if (!helper->SizeIsFixed(index) || !demands_helper->DemandIsFixed(index)) { num_variable_energies++; } starts.push_back(helper->Starts()[index]); @@ -811,24 +802,21 @@ void AddCumulativeRelaxation( if (helper->IsAbsent(i)) continue; if (!helper->IsOptional(i)) { - if (energies[i].has_value()) { - // The energy is defined if built from a constant value, a linear - // expression, or a linearized product. - lc.AddLinearExpression(energies[i].value()); + const std::vector& product = + demands_helper->DecomposedEnergies()[i]; + if (!product.empty()) { + // The energy is defined if the vector is not empty. + if (!lc.AddDecomposedProduct(product)) return; } else { - // The demand and the size are variable, and their product could not be - // linearized. - lc.AddQuadraticLowerBound(helper->Sizes()[i], demands[i], - integer_trail); + // The energy is not a decomposed product, but it could still be + // constant or linear. If not, a McCormick relaxation will be + // introduced. AddQuadraticLowerBound() supports all cases. + lc.AddQuadraticLowerBound(helper->Sizes()[i], + demands_helper->Demands()[i], integer_trail); } } else { - const IntegerValue product_min = - helper->SizeMin(i) * integer_trail->LowerBound(demands[i]); - const IntegerValue energy_min = energies[i].has_value() - ? energies[i]->Min(*integer_trail) - : IntegerValue(0); - if (!lc.AddLiteralTerm(helper->PresenceLiteral(i), - std::max(energy_min, product_min))) { + const IntegerValue energy_min = demands_helper->EnergyMin(i); + if (!lc.AddLiteralTerm(helper->PresenceLiteral(i), energy_min)) { return; } } @@ -879,10 +867,10 @@ void AppendNoOverlap2dRelaxation(const ConstraintProto& ct, Model* model, for (int i = 0; i < ct.no_overlap_2d().x_intervals_size(); ++i) { if (intervals_repository->IsPresent(x_intervals[i]) && intervals_repository->IsPresent(y_intervals[i])) { - LinearConstraintBuilder linear_energy(model); - if (DetectLinearEncodingOfProducts(x_sizes[i], y_sizes[i], model, - &linear_energy)) { - lc.AddLinearExpression(linear_energy.BuildExpression()); + const std::vector energy = + TryToDecomposeProduct(x_sizes[i], y_sizes[i], model); + if (!energy.empty()) { + if (!lc.AddDecomposedProduct(energy)) return; } else { lc.AddQuadraticLowerBound(x_sizes[i], y_sizes[i], integer_trail); } @@ -1197,6 +1185,10 @@ void TryToLinearizeConstraint(const CpModelProto& model_proto, break; } case ConstraintProto::ConstraintCase::kNoOverlap2D: { + // TODO(user): Use the same pattern as the other 2 scheduling methods: + // - single function + // - generate helpers once + // // Adds an energetic relaxation (sum of areas fits in bounding box). AppendNoOverlap2dRelaxation(ct, model, relaxation); if (linearization_level > 1) { @@ -1340,20 +1332,18 @@ bool IntervalIsVariable(const IntervalVariable interval, return false; } -void AddCumulativeCutGenerator( - SchedulingConstraintHelper* helper, - const std::vector& demands, - const AffineExpression& capacity, - const std::vector>& energies, - std::optional& makespan, Model* m, - LinearRelaxation* relaxation) { +void AddCumulativeCutGenerator(const AffineExpression& capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands_helper, + std::optional& makespan, + Model* m, LinearRelaxation* relaxation) { + relaxation->cut_generators.push_back(CreateCumulativeTimeTableCutGenerator( + helper, demands_helper, capacity, m)); relaxation->cut_generators.push_back( - CreateCumulativeTimeTableCutGenerator(helper, capacity, demands, m)); - relaxation->cut_generators.push_back( - CreateCumulativeCompletionTimeCutGenerator(helper, capacity, demands, - energies, m)); - relaxation->cut_generators.push_back( - CreateCumulativePrecedenceCutGenerator(helper, capacity, demands, m)); + CreateCumulativeCompletionTimeCutGenerator(helper, demands_helper, + capacity, m)); + relaxation->cut_generators.push_back(CreateCumulativePrecedenceCutGenerator( + helper, demands_helper, capacity, m)); // Checks if at least one rectangle has a variable size, is optional, or if // the demand or the capacity are variable. @@ -1365,14 +1355,14 @@ void AddCumulativeCutGenerator( break; } // Checks variable demand. - if (!integer_trail->IsFixed(demands[i])) { + if (!demands_helper->DemandIsFixed(i)) { has_variable_part = true; break; } } if (has_variable_part || !integer_trail->IsFixed(capacity)) { relaxation->cut_generators.push_back(CreateCumulativeEnergyCutGenerator( - helper, capacity, demands, energies, makespan, m)); + helper, demands_helper, capacity, makespan, m)); } } diff --git a/ortools/sat/linear_relaxation.h b/ortools/sat/linear_relaxation.h index 5400bac090..3eca84e0cb 100644 --- a/ortools/sat/linear_relaxation.h +++ b/ortools/sat/linear_relaxation.h @@ -17,13 +17,10 @@ #include #include "ortools/sat/cp_model.pb.h" -#include "ortools/sat/cp_model_mapping.h" #include "ortools/sat/cuts.h" #include "ortools/sat/integer.h" -#include "ortools/sat/linear_constraint.h" -#include "ortools/sat/linear_programming_constraint.h" +#include "ortools/sat/intervals.h" #include "ortools/sat/model.h" -#include "ortools/sat/sat_base.h" namespace operations_research { namespace sat { @@ -178,20 +175,16 @@ void AddNoOverlap2dCutGenerator(const ConstraintProto& ct, Model* m, // Adds linearization of cumulative constraints.The second part adds an // energetic equation linking the duration of all potential tasks to the actual // max span * capacity of the cumulative constraint. -void AddCumulativeRelaxation( - SchedulingConstraintHelper* helper, - const std::vector& demands, - const AffineExpression& capacity, - const std::vector>& energies, Model* model, - LinearRelaxation* relaxation); +void AddCumulativeRelaxation(const AffineExpression& capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, Model* model, + LinearRelaxation* relaxation); -void AddCumulativeCutGenerator( - SchedulingConstraintHelper* helper, - const std::vector& demands, - const AffineExpression& capacity, - const std::vector>& energies, - std::optional& makespan, Model* m, - LinearRelaxation* relaxation); +void AddCumulativeCutGenerator(const AffineExpression& capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, + std::optional& makespan, + Model* m, LinearRelaxation* relaxation); void AddNoOverlapCutGenerator(SchedulingConstraintHelper* helper, const std::optional& makespan, diff --git a/ortools/sat/lp_utils.cc b/ortools/sat/lp_utils.cc index 199f19de18..d8ae326290 100644 --- a/ortools/sat/lp_utils.cc +++ b/ortools/sat/lp_utils.cc @@ -853,9 +853,11 @@ bool ConvertMPModelProtoToCpModelProto(const SatParameters& params, // Notify if a continuous variable has a small domain as this is likely to // make an all integer solution far from a continuous one. - if (!mp_var.is_integer() && cp_var->domain(0) != cp_var->domain(1) && - cp_var->domain(1) - cp_var->domain(0) < kSmallDomainSize) { - ++num_small_domains; + if (!mp_var.is_integer()) { + const double diff = mp_var.upper_bound() - mp_var.lower_bound(); + if (diff > kWantedPrecision && diff < kSmallDomainSize) { + ++num_small_domains; + } } } diff --git a/ortools/sat/model.h b/ortools/sat/model.h index faf0424154..d64a679784 100644 --- a/ortools/sat/model.h +++ b/ortools/sat/model.h @@ -149,8 +149,9 @@ class Model { * It will be destroyed when the model is. */ template - void TakeOwnership(T* t) { + T* TakeOwnership(T* t) { cleanup_list_.emplace_back(new Delete(t)); + return t; } /** diff --git a/ortools/sat/parameters_validation.cc b/ortools/sat/parameters_validation.cc index 0c8b981c26..2935b2cfb9 100644 --- a/ortools/sat/parameters_validation.cc +++ b/ortools/sat/parameters_validation.cc @@ -59,6 +59,7 @@ std::string ValidateParameters(const SatParameters& params) { TEST_NOT_NAN(absolute_gap_limit); TEST_NOT_NAN(relative_gap_limit); TEST_NOT_NAN(log_frequency_in_seconds); + TEST_NOT_NAN(model_reduction_log_frequency_in_seconds); TEST_NOT_NAN(presolve_probing_deterministic_time_limit); TEST_NOT_NAN(merge_no_overlap_work_limit); TEST_NOT_NAN(merge_at_most_one_work_limit); diff --git a/ortools/sat/presolve_context.cc b/ortools/sat/presolve_context.cc index 730b51c638..6a314be040 100644 --- a/ortools/sat/presolve_context.cc +++ b/ortools/sat/presolve_context.cc @@ -1201,8 +1201,6 @@ void PresolveContext::InitializeNewDomains() { var_with_reduced_small_degree.Resize(domains.size()); var_to_constraints_.resize(domains.size()); var_to_num_linear1_.resize(domains.size()); - var_to_ub_only_constraints.resize(domains.size()); - var_to_lb_only_constraints.resize(domains.size()); } void PresolveContext::CanonicalizeDomainOfSizeTwo(int var) { diff --git a/ortools/sat/presolve_context.h b/ortools/sat/presolve_context.h index d98f534cef..80326b078d 100644 --- a/ortools/sat/presolve_context.h +++ b/ortools/sat/presolve_context.h @@ -516,17 +516,6 @@ class PresolveContext { TimeLimit* time_limit() { return time_limit_; } ModelRandomGenerator* random() { return random_; } - // For each variables, list the constraints that just enforce a lower bound - // (resp. upper bound) on that variable. If all the constraints in which a - // variable appear are in the same direction, then we can usually fix a - // variable to one of its bound (modulo its cost). - // - // TODO(user): Keeping these extra vector of hash_set seems inefficient. Come - // up with a better way to detect if a variable is only constrainted in one - // direction. - std::vector> var_to_ub_only_constraints; - std::vector> var_to_lb_only_constraints; - CpModelProto* working_model = nullptr; CpModelProto* mapping_model = nullptr; diff --git a/ortools/sat/probing.cc b/ortools/sat/probing.cc index eb1296dbb8..0b7e4e6ccf 100644 --- a/ortools/sat/probing.cc +++ b/ortools/sat/probing.cc @@ -76,6 +76,9 @@ bool Prober::ProbeOneVariableInternal(BooleanVariable b) { if (sat_solver_->IsModelUnsat()) return false; if (sat_solver_->CurrentDecisionLevel() == 0) continue; + if (trail_.Index() > saved_index) { + if (callback_ != nullptr) callback_(decision); + } implied_bounds_->ProcessIntegerTrail(decision); integer_trail_->AppendNewBounds(&new_integer_bounds_); diff --git a/ortools/sat/probing.h b/ortools/sat/probing.h index 005dabe257..e12cd04901 100644 --- a/ortools/sat/probing.h +++ b/ortools/sat/probing.h @@ -14,6 +14,7 @@ #ifndef OR_TOOLS_SAT_PROBING_H_ #define OR_TOOLS_SAT_PROBING_H_ +#include #include #include #include @@ -86,6 +87,13 @@ class Prober { int num_new_literals_fixed() const { return num_new_literals_fixed_; } int num_new_binary_clauses() const { return num_new_binary_; } + // Register a callback that will be called on each "propagation". + // One can inspect the VariablesAssignment to see what are the inferred + // literals. + void SetPropagationCallback(std::function f) { + callback_ = f; + } + private: bool ProbeOneVariableInternal(BooleanVariable b); @@ -114,6 +122,8 @@ class Prober { int num_new_integer_bounds_ = 0; int num_new_literals_fixed_ = 0; + std::function callback_ = nullptr; + // Logger. SolverLogger* logger_; }; diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 6f39a67271..9381e05d71 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: 218 +// NEXT TAG: 220 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -382,8 +382,9 @@ message SatParameters { // Indicates how much logging should wait before logging periodic search // information from specialized workers (lb_tree_search, probing). - // A value <= 0.0 will disable periodic logs. - optional double log_frequency_in_seconds = 212 [default = 0.0]; + // A value < 0 will disable periodic logs. + optional double log_frequency_in_seconds = 212 [default = -1]; + optional double model_reduction_log_frequency_in_seconds = 218 [default = 5]; // Whether the solver should display per sub-solver search statistics. // This is only useful is log_search_progress is set to true, and if the @@ -582,6 +583,11 @@ message SatParameters { // is an objective, or randomized SAT search for pure satisfiability problems. repeated string subsolvers = 207; + // A convenient way to add more workers types. + // Note that these will be added at the end of the list, so this is only + // useful if you run with many workers. + repeated string add_subsolvers = 219; + // Rather than fully specifying subsolvers, it is often convenient to just // remove the ones that are not useful on a given problem. repeated string ignore_subsolvers = 209; @@ -706,8 +712,7 @@ message SatParameters { // // This always result in better propagation, but it is usually slow, so // depending on the problem, turning this off may lead to a faster solution. - optional bool use_overload_checker_in_cumulative_constraint = 78 - [default = false]; + optional bool use_overload_checker_in_cumulative = 78 [default = false]; // When this is true, the cumulative constraint is reinforced with timetable // edge finding, i.e., an additional level of reasoning based on the @@ -716,13 +721,11 @@ message SatParameters { // // This always result in better propagation, but it is usually slow, so // depending on the problem, turning this off may lead to a faster solution. - optional bool use_timetable_edge_finding_in_cumulative_constraint = 79 - [default = false]; + optional bool use_timetable_edge_finding_in_cumulative = 79 [default = false]; // If true, detect and create constraint for integer variable that are "after" // a set of intervals in the same cumulative constraint. - optional bool use_hard_precedences_in_cumulative_constraint = 215 - [default = false]; + optional bool use_hard_precedences_in_cumulative = 215 [default = false]; // When this is true, the cumulative constraint is reinforced with propagators // from the disjunctive constraint to improve the inference on a set of tasks @@ -734,8 +737,7 @@ message SatParameters { // // This always result in better propagation, but it is usually slow, so // depending on the problem, turning this off may lead to a faster solution. - optional bool use_disjunctive_constraint_in_cumulative_constraint = 80 - [default = true]; + optional bool use_disjunctive_constraint_in_cumulative = 80 [default = true]; // When this is true, the no_overlap_2d constraint is reinforced with // propagators from the cumulative constraints. It consists of ignoring the diff --git a/ortools/sat/scheduling_cuts.cc b/ortools/sat/scheduling_cuts.cc index 166372cc5a..9744bdca83 100644 --- a/ortools/sat/scheduling_cuts.cc +++ b/ortools/sat/scheduling_cuts.cc @@ -49,20 +49,6 @@ namespace { // is needed to avoid numerical issues and adding cuts with minor effect. const double kMinCutViolation = 1e-4; -// Returns the lp value of a Literal. -double GetLiteralLpValue( - const Literal lit, - const absl::StrongVector& lp_values, - const IntegerEncoder* encoder) { - const IntegerVariable direct_view = encoder->GetLiteralView(lit); - if (direct_view != kNoIntegerVariable) { - return lp_values[direct_view]; - } - const IntegerVariable opposite_view = encoder->GetLiteralView(lit.Negated()); - DCHECK_NE(opposite_view, kNoIntegerVariable); - return 1.0 - lp_values[opposite_view]; -} - void AddIntegerVariableFromIntervals(SchedulingConstraintHelper* helper, Model* model, std::vector* vars) { @@ -80,17 +66,11 @@ void AddIntegerVariableFromIntervals(SchedulingConstraintHelper* helper, if (helper->IsOptional(t) && !helper->IsAbsent(t) && !helper->IsPresent(t)) { const Literal l = helper->PresenceLiteral(t); - if (encoder->GetLiteralView(l) == kNoIntegerVariable && - encoder->GetLiteralView(l.Negated()) == kNoIntegerVariable) { - model->Add(NewIntegerVariableFromLiteral(l)); - } - const IntegerVariable direct_view = encoder->GetLiteralView(l); - if (direct_view != kNoIntegerVariable) { - vars->push_back(direct_view); - } else { - vars->push_back(encoder->GetLiteralView(l.Negated())); - DCHECK_NE(vars->back(), kNoIntegerVariable); + IntegerVariable view = kNoIntegerVariable; + if (!encoder->LiteralOrNegationHasView(l, &view)) { + view = model->Add(NewIntegerVariableFromLiteral(l)); } + vars->push_back(view); } } } @@ -98,25 +78,41 @@ void AddIntegerVariableFromIntervals(SchedulingConstraintHelper* helper, } // namespace struct EnergyEvent { + // Cache of the intervals bound on the x direction. IntegerValue x_start_min; IntegerValue x_start_max; IntegerValue x_end_min; IntegerValue x_end_max; + + // Useful for no_overlap_2d. + IntegerValue y_min = IntegerValue(0); + IntegerValue y_max = IntegerValue(0); + + // Sizes in both dimensions. + // We also cache the minimum value to not recompute it. AffineExpression x_size; - IntegerValue y_min = IntegerValue(0); // Useful for no_overlap_2d. - IntegerValue y_max = IntegerValue(0); // Useful for no_overlap_2d. AffineExpression y_size; - - // Energy will be present only if the x_size * y_size could be linearized. - std::optional energy; - LiteralIndex presence_literal_index = kNoLiteralIndex; - - // Caches for MinOf(x_size), MinOf(y_size), and the LP values of the energy - // and the presence literal. IntegerValue x_size_min; IntegerValue y_size_min; - double literal_lp = 1.0; - double energy_lp = 0.0; + + // If set, this event is optional and its presence is controlled by this. + LiteralIndex presence_literal_index = kNoLiteralIndex; + + // The energy min of this event. + IntegerValue energy_min; + + // A linear expression which is a valid lower bound on the total energy of + // this event. We also cache the activity of the expression to not recompute + // it all the time. + LinearExpression linearized_energy; + double linearized_energy_lp_value = 0.0; + + // True if linearized_energy is not exact and a McCormick relaxation. + bool energy_is_quadratic = false; + + // If non empty, a decomposed view of the energy of this event. + // First value in each pair is x_size, second is y_size. + std::vector decomposed_energy; // Used to minimize the increase on the y axis for rectangles. double y_spread = 0.0; @@ -135,6 +131,64 @@ struct EnergyEvent { IntegerValue(0)); } + // This method expects all the other fields to have been filled before. + // It must be called before the EnergyEvent is used. + ABSL_MUST_USE_RESULT bool FillEnergyLp( + const absl::StrongVector& lp_values, + Model* model) { + LinearConstraintBuilder tmp_energy(model); + if (IsPresent()) { + if (!decomposed_energy.empty()) { + if (!tmp_energy.AddDecomposedProduct(decomposed_energy)) return false; + } else { + tmp_energy.AddQuadraticLowerBound(x_size, y_size, + model->GetOrCreate(), + &energy_is_quadratic); + } + } else { + if (!tmp_energy.AddLiteralTerm(Literal(presence_literal_index), + energy_min)) { + return false; + } + } + linearized_energy = tmp_energy.BuildExpression(); + linearized_energy_lp_value = linearized_energy.LpValue(lp_values); + return true; + } + + IntegerValue EnergyMinInWindow(const VariablesAssignment& assignment, + IntegerValue window_start, + IntegerValue window_end) const { + if (window_end <= window_start) return IntegerValue(0); + + // Returns zero if the interval do not necessarily overlap. + if (x_end_min <= window_start) return IntegerValue(0); + if (x_start_max >= window_end) return IntegerValue(0); + const IntegerValue window_size = window_end - window_start; + const IntegerValue simple_energy_min = + y_size_min * + std::min({x_end_min - window_start, window_end - x_start_max, + x_size_min, window_size}); + if (decomposed_energy.empty()) return simple_energy_min; + + IntegerValue result = kMaxIntegerValue; + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energy) { + if (assignment.LiteralIsTrue(lit)) { + // Both should be identical, so we don't recompute it. + return simple_energy_min; + } + if (assignment.LiteralIsFalse(lit)) continue; + const IntegerValue alt_em = std::max(x_end_min, x_start_min + fixed_size); + const IntegerValue alt_sm = std::min(x_start_max, x_end_max - fixed_size); + const IntegerValue energy_min = + fixed_demand * std::min({alt_em - window_start, window_end - alt_sm, + fixed_size, window_size}); + result = std::min(result, energy_min); + } + if (result == kMaxIntegerValue) return simple_energy_min; + return std::max(simple_energy_min, result); + } + std::string DebugString() const { return absl::StrCat( "EnergyEvent(x_start_min = ", x_start_min.value(), @@ -143,7 +197,10 @@ struct EnergyEvent { ", x_end_max = ", x_end_max.value(), ", x_size = ", x_size.DebugString(), ", y_min = ", y_min.value(), ", y_max = ", y_max.value(), ", y_size = ", y_size.DebugString(), - ", energy = ", energy ? energy.value().DebugString() : "{}", + ", energy = ", + decomposed_energy.empty() + ? "{}" + : absl::StrCat(decomposed_energy.size(), " terms"), ", presence_literal_index = ", presence_literal_index.value(), ")"); } }; @@ -169,24 +226,24 @@ void GenerateCumulativeEnergeticCuts( time_points_set.end()); double sum_of_energies_lp = 0.0; - for (const EnergyEvent& event : events) { // Min contribution is zero. - sum_of_energies_lp += event.energy_lp; + for (const EnergyEvent& event : events) { + sum_of_energies_lp += event.linearized_energy_lp_value; } - // Compute the energetic contribution of a task in a given time window, and - // add it either to energy_lp, or to the cut. - // It returns false if it tried to generate the cut, and failed. - IntegerTrail* integer_trail = model->GetOrCreate(); - // Checks the precondition of the code. + IntegerTrail* integer_trail = model->GetOrCreate(); DCHECK(!makespan.has_value() || integer_trail->IsFixed(capacity)); - const auto add_one_event = [integer_trail, &lp_values]( + // Compute the energetic contribution of a task in a given time window, and + // add it to the cut. It returns false if it tried to generate the cut, and + // failed. + const VariablesAssignment& assignment = + model->GetOrCreate()->Assignment(); + const auto add_one_event = [&assignment]( const EnergyEvent& event, IntegerValue window_start, IntegerValue window_end, - double* energy_lp = nullptr, - LinearConstraintBuilder* cut = nullptr, + LinearConstraintBuilder* cut, bool* add_energy_to_name = nullptr, bool* add_quadratic_to_name = nullptr, bool* add_opt_to_name = nullptr, @@ -195,85 +252,60 @@ void GenerateCumulativeEnergeticCuts( return true; // Event can move outside the time window. } - const bool compute_energy_only = energy_lp != nullptr; - if (compute_energy_only) { - DCHECK(cut == nullptr); - DCHECK(add_energy_to_name == nullptr); - DCHECK(add_quadratic_to_name == nullptr); - DCHECK(add_opt_to_name == nullptr); - DCHECK(add_lifted_to_name == nullptr); - } else { - DCHECK(cut != nullptr); - DCHECK(add_energy_to_name != nullptr); - DCHECK(add_quadratic_to_name != nullptr); - DCHECK(add_opt_to_name != nullptr); - DCHECK(add_lifted_to_name != nullptr); - } + DCHECK(cut != nullptr); if (event.x_start_min >= window_start && event.x_end_max <= window_end) { // Event is always contained by the time window. - if (event.IsPresent()) { - if (event.energy) { - // We favor the energy info instead of the McCormick relaxation. - if (compute_energy_only) { - *energy_lp += event.energy.value().LpValue(lp_values); - } else { - cut->AddLinearExpression(event.energy.value()); - *add_energy_to_name = true; - } - } else { // The energy was not linearized. Thus it is quadratic. - if (compute_energy_only) { - *energy_lp += - ToDouble(event.y_size_min) * event.x_size.LpValue(lp_values) + - ToDouble(event.x_size_min) * event.y_size.LpValue(lp_values) - - ToDouble(event.x_size_min * event.y_size_min); - } else { - cut->AddQuadraticLowerBound(event.x_size, event.y_size, - integer_trail); - *add_quadratic_to_name = true; - } - } - } else { // Event is optional. - const IntegerValue min_energy = std::max( - event.x_size_min * event.y_size_min, - event.energy ? event.energy.value().LevelZeroMin(integer_trail) - : IntegerValue(0)); - if (compute_energy_only) { - *energy_lp += event.literal_lp * ToDouble(min_energy); - } else { - if (!cut->AddLiteralTerm(Literal(event.presence_literal_index), - min_energy)) { - return false; - } + cut->AddLinearExpression(event.linearized_energy); - *add_opt_to_name = true; - if (min_energy > event.x_size_min * event.y_size_min) { - *add_energy_to_name = true; - } - } + if (event.energy_is_quadratic && add_quadratic_to_name != nullptr) { + *add_quadratic_to_name = true; + } + if (add_energy_to_name != nullptr && + event.energy_min > event.x_size_min * event.y_size_min) { + *add_energy_to_name = true; + } + if (!event.IsPresent() && add_opt_to_name != nullptr) { + *add_opt_to_name = true; } } else { // The event has a mandatory overlap with the time window. const IntegerValue min_overlap = event.GetMinOverlap(window_start, window_end); if (min_overlap <= 0) return true; + if (add_lifted_to_name != nullptr) *add_lifted_to_name = true; if (event.IsPresent()) { - if (compute_energy_only) { - *energy_lp += event.y_size.LpValue(lp_values) * ToDouble(min_overlap); - } else { + const std::vector& energy = event.decomposed_energy; + if (energy.empty()) { cut->AddTerm(event.y_size, min_overlap); - *add_lifted_to_name = true; + } else { + const IntegerValue window_size = window_end - window_start; + for (const auto [lit, fixed_size, fixed_demand] : energy) { + if (assignment.LiteralIsFalse(lit)) continue; + const IntegerValue alt_em = + std::max(event.x_end_min, event.x_start_min + fixed_size); + const IntegerValue alt_sm = + std::min(event.x_start_max, event.x_end_max - fixed_size); + const IntegerValue energy_min = + fixed_demand * + std::min({alt_em - window_start, window_end - alt_sm, + fixed_size, window_size}); + DCHECK_GT(energy_min, 0); + if (!cut->AddLiteralTerm(lit, energy_min)) return false; + } + if (add_energy_to_name != nullptr) *add_energy_to_name = true; } } else { - if (compute_energy_only) { - *energy_lp += - event.literal_lp * ToDouble(min_overlap * event.y_size_min); - } else { - if (!cut->AddLiteralTerm(Literal(event.presence_literal_index), - min_overlap * event.y_size_min)) { - return false; - } - *add_opt_to_name = true; + if (add_opt_to_name != nullptr) *add_opt_to_name = true; + const IntegerValue min_energy = + event.EnergyMinInWindow(assignment, window_start, window_end); + if (min_energy > event.x_size_min * event.y_size_min && + add_energy_to_name != nullptr) { + *add_energy_to_name = true; + } + if (!cut->AddLiteralTerm(Literal(event.presence_literal_index), + min_energy)) { + return false; } } } @@ -298,6 +330,7 @@ void GenerateCumulativeEnergeticCuts( : std::numeric_limits::infinity(); const int num_time_points = time_points.size(); + LinearConstraintBuilder tmp_energy(model); for (int i = 0; i + 1 < num_time_points; ++i) { const IntegerValue window_start = time_points[i]; // After max_end_min, all tasks can fit before window_start. @@ -319,10 +352,12 @@ void GenerateCumulativeEnergeticCuts( double energy_lp = 0.0; bool energy_correctly_computed = true; for (const EnergyEvent& event : events) { - if (!add_one_event(event, window_start, window_end, &energy_lp)) { + tmp_energy.Clear(); + if (!add_one_event(event, window_start, window_end, &tmp_energy)) { energy_correctly_computed = false; break; // Abort. } + energy_lp += tmp_energy.BuildExpression().LpValue(lp_values); } if (!energy_correctly_computed) continue; @@ -365,10 +400,9 @@ void GenerateCumulativeEnergeticCuts( cut.AddTerm(capacity, window_start - window_end); } for (const EnergyEvent& event : events) { - if (!add_one_event(event, window_start, window_end, - /*energy_lp=*/nullptr, &cut, &add_energy_to_name, - &add_quadratic_to_name, &add_opt_to_name, - &add_lifted_to_name)) { + if (!add_one_event(event, window_start, window_end, &cut, + &add_energy_to_name, &add_quadratic_to_name, + &add_opt_to_name, &add_lifted_to_name)) { cut_generated = false; break; // Exit the event loop. } @@ -388,75 +422,57 @@ void GenerateCumulativeEnergeticCuts( top_n_cuts.TransferToManager(lp_values, manager); } -void AppendVariablesToCumulativeCut( - const AffineExpression& capacity, - const std::vector& demands, IntegerTrail* integer_trail, - CutGenerator* result) { - for (const AffineExpression& demand_expr : demands) { +void AppendVariablesToCumulativeCut(const AffineExpression& capacity, + SchedulingDemandHelper* demands_helper, + Model* model, + std::vector* vars) { + auto* integer_trail = model->GetOrCreate(); + for (const AffineExpression& demand_expr : demands_helper->Demands()) { if (!integer_trail->IsFixed(demand_expr)) { - result->vars.push_back(demand_expr.var); + vars->push_back(demand_expr.var); } } - if (!integer_trail->IsFixed(capacity)) { - result->vars.push_back(capacity.var); + IntegerEncoder* encoder = model->GetOrCreate(); + for (const auto& product : demands_helper->DecomposedEnergies()) { + for (const auto& lit_val_val : product) { + IntegerVariable view = kNoIntegerVariable; + if (!encoder->LiteralOrNegationHasView(lit_val_val.literal, &view)) { + view = model->Add(NewIntegerVariableFromLiteral(lit_val_val.literal)); + } + vars->push_back(view); + } } -} -double ComputeEnergyLp( - const EnergyEvent& e, - const absl::StrongVector& lp_values, - IntegerTrail* integer_trail) { - if (e.presence_literal_index == kNoLiteralIndex) { - if (e.energy) { - return e.energy.value().LpValue(lp_values); - } else { // demand and size are not fixed. - // X * Y >= X * y_min + x_min * Y - x_min * y_min. - return ToDouble(e.y_size_min) * e.x_size.LpValue(lp_values) + - ToDouble(e.x_size_min) * e.y_size.LpValue(lp_values) - - ToDouble(e.x_size_min * e.y_size_min); - } - } else { - const IntegerValue min_energy = - std::max(e.x_size_min * e.y_size_min, - e.energy ? e.energy.value().LevelZeroMin(integer_trail) - : IntegerValue(0)); - return e.literal_lp * ToDouble(min_energy); + if (!integer_trail->IsFixed(capacity)) { + vars->push_back(capacity.var); } } CutGenerator CreateCumulativeEnergyCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, - const std::vector>& energies, + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, const std::optional& makespan, Model* model) { CutGenerator result; Trail* trail = model->GetOrCreate(); IntegerEncoder* encoder = model->GetOrCreate(); - IntegerTrail* integer_trail = model->GetOrCreate(); - AppendVariablesToCumulativeCut(capacity, demands, integer_trail, &result); + AppendVariablesToCumulativeCut(capacity, demands_helper, model, &result.vars); AddIntegerVariableFromIntervals(helper, model, &result.vars); - for (const auto& energy : energies) { - if (!energy.has_value()) continue; - result.vars.insert(result.vars.end(), energy.value().vars.begin(), - energy.value().vars.end()); - } gtl::STLSortAndRemoveDuplicates(&result.vars); result.generate_cuts = - [makespan, capacity, demands, energies, trail, helper, model, encoder]( + [makespan, capacity, demands_helper, trail, helper, model, encoder]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { if (trail->CurrentDecisionLevel() > 0) return true; if (!helper->SynchronizeAndSetTimeDirection(true)) return false; - - IntegerTrail* integer_trail = model->GetOrCreate(); + demands_helper->CacheAllEnergyValues(); std::vector events; for (int i = 0; i < helper->NumTasks(); ++i) { if (helper->IsAbsent(i)) continue; - if (integer_trail->LevelZeroUpperBound(demands[i]) == 0 || - helper->SizeMin(i) == 0) { + // TODO(user): use level 0 bounds ? + if (demands_helper->DemandMax(i) == 0 || helper->SizeMin(i) == 0) { continue; } @@ -466,16 +482,17 @@ CutGenerator CreateCumulativeEnergyCutGenerator( e.x_end_min = helper->EndMin(i); e.x_end_max = helper->EndMax(i); e.x_size = helper->Sizes()[i]; - e.y_size = demands[i]; - e.energy = energies[i]; + e.y_size = demands_helper->Demands()[i]; + e.decomposed_energy = demands_helper->DecomposedEnergies()[i]; e.x_size_min = helper->SizeMin(i); - e.y_size_min = integer_trail->LevelZeroLowerBound(demands[i]); + e.y_size_min = demands_helper->DemandMin(i); + e.energy_min = demands_helper->EnergyMin(i); + e.energy_is_quadratic = demands_helper->EnergyIsQuadratic(i); if (!helper->IsPresent(i)) { e.presence_literal_index = helper->PresenceLiteral(i).Index(); - e.literal_lp = GetLiteralLpValue(Literal(e.presence_literal_index), - lp_values, encoder); } - e.energy_lp = ComputeEnergyLp(e, lp_values, integer_trail); + // We can always skip events. + if (!e.FillEnergyLp(lp_values, model)) continue; events.push_back(e); } @@ -528,17 +545,14 @@ CutGenerator CreateNoOverlapEnergyCutGenerator( e.x_end_max = helper->EndMax(i); e.x_size = helper->Sizes()[i]; e.y_size = IntegerValue(1); - e.energy = sizes[i]; e.x_size_min = helper->SizeMin(i); e.y_size_min = IntegerValue(1); - if (helper->IsPresent(i)) { - e.energy_lp = e.energy->LpValue(lp_values); - } else { + e.energy_min = e.x_size_min; + if (!helper->IsPresent(i)) { e.presence_literal_index = helper->PresenceLiteral(i).Index(); - e.literal_lp = GetLiteralLpValue(Literal(e.presence_literal_index), - lp_values, encoder); - e.energy_lp = e.literal_lp * ToDouble(e.x_size_min); } + // We can always skip events. + if (!e.FillEnergyLp(lp_values, model)) continue; events.push_back(e); } @@ -551,12 +565,12 @@ CutGenerator CreateNoOverlapEnergyCutGenerator( } void GenerateNoOverlap2dEnergyCut( - const std::vector>& energies, + const std::vector>& energies, absl::Span rectangles, const std::string& cut_name, const absl::StrongVector& lp_values, Model* model, - IntegerTrail* integer_trail, IntegerEncoder* encoder, - LinearConstraintManager* manager, SchedulingConstraintHelper* x_helper, - SchedulingConstraintHelper* y_helper) { + IntegerTrail* integer_trail, LinearConstraintManager* manager, + SchedulingConstraintHelper* x_helper, SchedulingConstraintHelper* y_helper, + SchedulingDemandHelper* y_demands_helper) { std::vector events; for (const int rect : rectangles) { if (y_helper->SizeMax(rect) == 0 || x_helper->SizeMax(rect) == 0) { @@ -572,20 +586,20 @@ void GenerateNoOverlap2dEnergyCut( e.y_min = y_helper->StartMin(rect); e.y_max = y_helper->EndMax(rect); e.y_size = y_helper->Sizes()[rect]; - e.energy = energies[rect]; + e.decomposed_energy = energies[rect]; e.presence_literal_index = x_helper->IsPresent(rect) ? (y_helper->IsPresent(rect) ? kNoLiteralIndex : y_helper->PresenceLiteral(rect).Index()) : x_helper->PresenceLiteral(rect).Index(); - if (e.presence_literal_index != kNoLiteralIndex) { - e.literal_lp = GetLiteralLpValue(Literal(e.presence_literal_index), - lp_values, encoder); - } e.x_size_min = x_helper->SizeMin(rect); e.y_size_min = y_helper->SizeMin(rect); - e.energy_lp = ComputeEnergyLp(e, lp_values, integer_trail); + e.energy_min = y_demands_helper->EnergyMin(rect); + e.energy_is_quadratic = y_demands_helper->EnergyIsQuadratic(rect); + + // We can always skip events. + if (!e.FillEnergyLp(lp_values, model)) continue; events.push_back(e); } @@ -613,7 +627,7 @@ void GenerateNoOverlap2dEnergyCut( // The sum of all energies can be used to stop iterating early. double sum_of_all_energies = 0.0; for (const auto& e : events) { - sum_of_all_energies += e.energy_lp; + sum_of_all_energies += e.linearized_energy_lp_value; } CapacityProfile capacity_profile; @@ -649,7 +663,7 @@ void GenerateNoOverlap2dEnergyCut( // each step. We follow the same structure as the cut creation code below. for (int i2 = 0; i2 < residual_events.size(); ++i2) { const EnergyEvent& e = residual_events[i2]; - energy_lp += e.energy_lp; + energy_lp += e.linearized_energy_lp_value; window_min = std::min(window_min, e.x_start_min); window_max = std::max(window_max, e.x_end_max); y_min = std::min(y_min, e.y_min); @@ -702,50 +716,26 @@ void GenerateNoOverlap2dEnergyCut( if (max_violation_end_index == -1) continue; // A maximal violated cut has been found. - bool cut_generated = true; + // Build it and add it to the pool. bool add_opt_to_name = false; bool add_quadratic_to_name = false; bool add_energy_to_name = false; - - // Build the cut. LinearConstraintBuilder cut(model, kMinIntegerValue, max_violation_area); - for (int i2 = 0; i2 <= max_violation_end_index; ++i2) { - const EnergyEvent& e = residual_events[i2]; - if (e.IsPresent()) { - if (e.energy) { - // We favor the energy info instead of the McCormick relaxation. - cut.AddLinearExpression(e.energy.value()); - add_energy_to_name = true; - } else { - cut.AddQuadraticLowerBound(e.x_size, e.y_size, integer_trail); - add_quadratic_to_name = true; - } - } else { - add_opt_to_name = true; - const IntegerValue min_energy = - std::max(e.x_size_min * e.y_size_min, - e.energy ? e.energy.value().LevelZeroMin(integer_trail) - : IntegerValue(0)); - if (min_energy > e.x_size_min * e.y_size_min) { - add_energy_to_name = true; - } - if (!cut.AddLiteralTerm(Literal(e.presence_literal_index), - min_energy)) { - cut_generated = false; - break; - } + const EnergyEvent& event = residual_events[i2]; + cut.AddLinearExpression(event.linearized_energy); + if (!event.IsPresent()) add_opt_to_name = true; + if (event.energy_is_quadratic) add_quadratic_to_name = true; + if (event.energy_min > event.x_size_min * event.y_size_min) { + add_energy_to_name = true; } } - - if (cut_generated) { - std::string full_name = cut_name; - if (add_opt_to_name) full_name.append("_opt"); - if (add_quadratic_to_name) full_name.append("_quadratic"); - if (add_energy_to_name) full_name.append("_energy"); - if (max_violation_use_precise_area) full_name.append("_precise"); - top_n_cuts.AddCut(cut.Build(), full_name, lp_values); - } + std::string full_name = cut_name; + if (add_opt_to_name) full_name.append("_opt"); + if (add_quadratic_to_name) full_name.append("_quadratic"); + if (add_energy_to_name) full_name.append("_energy"); + if (max_violation_use_precise_area) full_name.append("_precise"); + top_n_cuts.AddCut(cut.Build(), full_name, lp_values); } top_n_cuts.TransferToManager(lp_values, manager); } @@ -758,14 +748,14 @@ CutGenerator CreateNoOverlap2dEnergyCutGenerator( model->GetOrCreate(); const int num_rectangles = x_intervals.size(); - std::vector> energies; + std::vector> energies; std::vector x_sizes; std::vector y_sizes; for (int i = 0; i < num_rectangles; ++i) { x_sizes.push_back(intervals_repository->Size(x_intervals[i])); y_sizes.push_back(intervals_repository->Size(y_intervals[i])); energies.push_back( - TryToLinearizeProduct(x_sizes.back(), y_sizes.back(), model)); + TryToDecomposeProduct(x_sizes.back(), y_sizes.back(), model)); } SchedulingConstraintHelper* x_helper = @@ -773,22 +763,32 @@ CutGenerator CreateNoOverlap2dEnergyCutGenerator( SchedulingConstraintHelper* y_helper = model->GetOrCreate()->GetOrCreateHelper(y_intervals); + SchedulingDemandHelper* x_demands_helper = + new SchedulingDemandHelper(x_helper->Sizes(), y_helper, model); + model->TakeOwnership(x_demands_helper); + + SchedulingDemandHelper* y_demands_helper = + new SchedulingDemandHelper(y_helper->Sizes(), x_helper, model); + model->TakeOwnership(y_demands_helper); + AddIntegerVariableFromIntervals(x_helper, model, &result.vars); AddIntegerVariableFromIntervals(y_helper, model, &result.vars); gtl::STLSortAndRemoveDuplicates(&result.vars); Trail* trail = model->GetOrCreate(); IntegerTrail* integer_trail = model->GetOrCreate(); - IntegerEncoder* encoder = model->GetOrCreate(); result.generate_cuts = - [integer_trail, trail, encoder, x_helper, y_helper, model, energies]( - const absl::StrongVector& lp_values, - LinearConstraintManager* manager) { + [integer_trail, trail, x_helper, y_helper, x_demands_helper, + y_demands_helper, model, + energies](const absl::StrongVector& lp_values, + LinearConstraintManager* manager) { if (trail->CurrentDecisionLevel() > 0) return true; if (!x_helper->SynchronizeAndSetTimeDirection(true)) return false; if (!y_helper->SynchronizeAndSetTimeDirection(true)) return false; + x_demands_helper->CacheAllEnergyValues(); + y_demands_helper->CacheAllEnergyValues(); const int num_rectangles = x_helper->NumTasks(); std::vector active_rectangles; @@ -828,10 +828,10 @@ CutGenerator CreateNoOverlap2dEnergyCutGenerator( GenerateNoOverlap2dEnergyCut( energies, rectangles, "NoOverlap2dXEnergy", lp_values, model, - integer_trail, encoder, manager, x_helper, y_helper); + integer_trail, manager, x_helper, y_helper, y_demands_helper); GenerateNoOverlap2dEnergyCut( energies, rectangles, "NoOverlap2dYEnergy", lp_values, model, - integer_trail, encoder, manager, y_helper, x_helper); + integer_trail, manager, y_helper, x_helper, x_demands_helper); } return true; @@ -840,12 +840,12 @@ CutGenerator CreateNoOverlap2dEnergyCutGenerator( } CutGenerator CreateCumulativeTimeTableCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, Model* model) { + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, Model* model) { CutGenerator result; IntegerTrail* integer_trail = model->GetOrCreate(); - AppendVariablesToCumulativeCut(capacity, demands, integer_trail, &result); + AppendVariablesToCumulativeCut(capacity, demands_helper, model, &result.vars); AddIntegerVariableFromIntervals(helper, model, &result.vars); gtl::STLSortAndRemoveDuplicates(&result.vars); @@ -860,7 +860,7 @@ CutGenerator CreateCumulativeTimeTableCutGenerator( Trail* trail = model->GetOrCreate(); result.generate_cuts = - [helper, capacity, demands, trail, integer_trail, model]( + [helper, capacity, demands_helper, trail, integer_trail, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { if (trail->CurrentDecisionLevel() > 0) return true; @@ -880,7 +880,7 @@ CutGenerator CreateCumulativeTimeTableCutGenerator( TimeTableEvent e1; e1.interval_index = i; e1.time = start_max; - e1.demand = demands[i]; + e1.demand = demands_helper->Demands()[i]; e1.positive = true; TimeTableEvent e2 = e1; @@ -1064,12 +1064,12 @@ void GenerateCutsBetweenPairOfNonOverlappingTasks( } CutGenerator CreateCumulativePrecedenceCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, Model* model) { + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, Model* model) { CutGenerator result; IntegerTrail* integer_trail = model->GetOrCreate(); - AppendVariablesToCumulativeCut(capacity, demands, integer_trail, &result); + AppendVariablesToCumulativeCut(capacity, demands_helper, model, &result.vars); AddIntegerVariableFromIntervals(helper, model, &result.vars); gtl::STLSortAndRemoveDuplicates(&result.vars); @@ -1077,7 +1077,7 @@ CutGenerator CreateCumulativePrecedenceCutGenerator( Trail* trail = model->GetOrCreate(); result.generate_cuts = - [trail, integer_trail, helper, demands, capacity, model]( + [trail, integer_trail, helper, demands_helper, capacity, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { if (trail->CurrentDecisionLevel() > 0) return true; @@ -1093,7 +1093,7 @@ CutGenerator CreateCumulativePrecedenceCutGenerator( event.end_min = helper->EndMin(t); event.end_max = helper->EndMax(t); event.end = helper->Ends()[t]; - event.demand_min = integer_trail->LowerBound(demands[t]); + event.demand_min = demands_helper->DemandMin(t); event.duration_min = helper->SizeMin(t); events.push_back(event); } @@ -1171,10 +1171,16 @@ struct CtEvent { // The end max of the y interval. IntegerValue y_end_max; + // The size min of the y interval. + IntegerValue y_size_min; + // The min energy of the task (this is always larger or equal to x_size_min * // y_size_min). IntegerValue energy_min; + // The decomposed energy of the product. + std::vector decomposed_energy; + // Indicates if the events used the optional energy information from the // model. bool use_energy = false; @@ -1189,6 +1195,32 @@ struct CtEvent { // will not use this. IntegerValue fixed_y_size = -1; + IntegerValue EnergyMinAfter(const VariablesAssignment& assignment, + IntegerValue window_start) const { + // Returns zero if the interval do not necessarily overlap. + if (x_start_min + x_size_min <= window_start) return IntegerValue(0); + const IntegerValue size_reduction = window_start - x_start_min; + const IntegerValue simple_energy_min = + y_size_min * (x_size_min - size_reduction); + DCHECK_GT(simple_energy_min, 0); + if (decomposed_energy.empty()) return simple_energy_min; + + IntegerValue result = kMaxIntegerValue; + for (const auto [lit, fixed_size, fixed_demand] : decomposed_energy) { + if (assignment.LiteralIsTrue(lit)) { + // Both should be identical, so we don't recompute it. + return simple_energy_min; + } + if (assignment.LiteralIsFalse(lit)) continue; + const IntegerValue energy_min = + fixed_demand * (fixed_size - size_reduction); + DCHECK_GT(energy_min, 0); + result = std::min(result, energy_min); + } + if (result == kMaxIntegerValue) return simple_energy_min; + return std::max(simple_energy_min, result); + } + std::string DebugString() const { return absl::StrCat("CtEvent(x_end = ", x_end.DebugString(), ", x_start_min = ", x_start_min.value(), @@ -1219,6 +1251,8 @@ void GenerateCompletionTimeCuts( std::vector events, bool use_lifting, Model* model, LinearConstraintManager* manager) { TopNCuts top_n_cuts(5); + const VariablesAssignment& assignment = + model->GetOrCreate()->Assignment(); // Sort by start min to bucketize by start_min. std::sort(events.begin(), events.end(), @@ -1247,21 +1281,15 @@ void GenerateCompletionTimeCuts( // Build the vector of energies as the vector of sizes. CtEvent event = events[before]; // Copy. event.lifted = true; - const IntegerValue old_size_min = event.x_size_min; event.x_size_min = event.x_size_min + event.x_start_min - sequence_start_min; event.x_start_min = sequence_start_min; - // We can rescale the energy min correctly. - // - // Let's take the example of a rectangle of size 2 * 20 that overlaps - // sequence start min by 1, and that can rotate by 90 degrees. - // The energy min is 40, size min is 2, size_max is 20. - // If the rectangle is horizontal, the lifted energy is: - // (20 - 1) * 2 = 38 - // If the rectangle is vertical, the lifted energy is: - // (2 - 1) * 20 = 20 - // The min of the two is always reached when size = size_min. - event.energy_min = event.energy_min * event.x_size_min / old_size_min; + event.energy_min = + event.EnergyMinAfter(assignment, sequence_start_min); + if (event.energy_min > event.x_size_min * event.y_size_min) { + event.use_energy = true; + } + if (event.energy_min <= 0) continue; residual_tasks.push_back(event); } } @@ -1391,12 +1419,13 @@ CutGenerator CreateNoOverlapCompletionTimeCutGenerator( event.x_lp_end = end_expr.LpValue(lp_values); event.y_start_min = IntegerValue(0); event.y_end_max = IntegerValue(1); + event.y_size_min = IntegerValue(1); event.energy_min = size_min; events.push_back(event); } } GenerateCompletionTimeCuts(cut_name, lp_values, std::move(events), - /*use_lifting=*/false, model, manager); + /*use_lifting=*/true, model, manager); }; if (!helper->SynchronizeAndSetTimeDirection(true)) return false; generate_cuts("NoOverlapCompletionTime"); @@ -1408,40 +1437,36 @@ CutGenerator CreateNoOverlapCompletionTimeCutGenerator( } CutGenerator CreateCumulativeCompletionTimeCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, - const std::vector>& energies, - Model* model) { + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, Model* model) { CutGenerator result; IntegerTrail* integer_trail = model->GetOrCreate(); - AppendVariablesToCumulativeCut(capacity, demands, integer_trail, &result); - + AppendVariablesToCumulativeCut(capacity, demands_helper, model, &result.vars); AddIntegerVariableFromIntervals(helper, model, &result.vars); gtl::STLSortAndRemoveDuplicates(&result.vars); Trail* trail = model->GetOrCreate(); result.generate_cuts = - [trail, integer_trail, helper, demands, energies, capacity, model]( + [trail, integer_trail, helper, demands_helper, capacity, model]( const absl::StrongVector& lp_values, LinearConstraintManager* manager) { if (trail->CurrentDecisionLevel() > 0) return true; if (!helper->SynchronizeAndSetTimeDirection(true)) return false; const IntegerValue capacity_max = integer_trail->UpperBound(capacity); - auto generate_cuts = [&lp_values, model, manager, helper, capacity_max, - integer_trail, &demands, - &energies](const std::string& cut_name) { + auto generate_cuts = [&lp_values, model, manager, helper, + demands_helper, capacity_max, + integer_trail](const std::string& cut_name) { std::vector events; for (int index = 0; index < helper->NumTasks(); ++index) { if (!helper->IsPresent(index)) continue; if (helper->SizeMin(index) > 0 && - integer_trail->LowerBound(demands[index]) > 0) { + demands_helper->DemandMin(index) > 0) { const AffineExpression end_expr = helper->Ends()[index]; const IntegerValue size_min = helper->SizeMin(index); - const IntegerValue demand_min = - integer_trail->LowerBound(demands[index]); + const IntegerValue demand_min = demands_helper->DemandMin(index); CtEvent event; event.x_start_min = helper->StartMin(index); @@ -1450,20 +1475,13 @@ CutGenerator CreateCumulativeCompletionTimeCutGenerator( event.x_lp_end = end_expr.LpValue(lp_values); event.y_start_min = IntegerValue(0); event.y_end_max = IntegerValue(capacity_max); - event.energy_min = size_min * demand_min; - if (integer_trail->IsFixed(demands[index])) { + event.y_size_min = demand_min; + event.energy_min = demands_helper->EnergyMin(index); + event.decomposed_energy = + demands_helper->DecomposedEnergies()[index]; + if (demands_helper->DemandIsFixed(index)) { event.fixed_y_size = demand_min; } - - // TODO(user): Investigate and re-enable a correct version. - if (/*DISABLES_CODE*/ (false) && energies[index].has_value()) { - const IntegerValue linearized_energy = - energies[index]->Min(*integer_trail); - if (linearized_energy > event.energy_min) { - event.energy_min = linearized_energy; - event.use_energy = true; - } - } events.push_back(event); } } @@ -1479,6 +1497,7 @@ CutGenerator CreateCumulativeCompletionTimeCutGenerator( return result; } +// TODO(user): Use demands_helper and decomposed energy. CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( const std::vector& x_intervals, const std::vector& y_intervals, Model* model) { @@ -1553,13 +1572,17 @@ CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator( event.x_lp_end = x_end_expr.LpValue(lp_values); event.y_start_min = y_helper->ShiftedStartMin(rect); event.y_end_max = y_helper->ShiftedEndMax(rect); + event.y_size_min = y_helper->SizeMin(rect); + // TODO(user): Use improved energy from demands helper. event.energy_min = x_helper->SizeMin(rect) * y_helper->SizeMin(rect); + event.decomposed_energy = TryToDecomposeProduct( + x_helper->Sizes()[rect], y_helper->Sizes()[rect], model); events.push_back(event); } GenerateCompletionTimeCuts(cut_name, lp_values, std::move(events), - /*use_lifting=*/true, model, manager); + /*use_lifting=*/false, model, manager); }; if (!x_helper->SynchronizeAndSetTimeDirection(true)) return false; diff --git a/ortools/sat/scheduling_cuts.h b/ortools/sat/scheduling_cuts.h index 77da8522e9..c5c339e8cd 100644 --- a/ortools/sat/scheduling_cuts.h +++ b/ortools/sat/scheduling_cuts.h @@ -19,15 +19,9 @@ #include #include -#include "ortools/base/strong_vector.h" #include "ortools/sat/cuts.h" #include "ortools/sat/integer.h" #include "ortools/sat/intervals.h" -#include "ortools/sat/linear_constraint.h" -#include "ortools/sat/linear_constraint_manager.h" -#include "ortools/sat/model.h" -#include "ortools/util/strong_integers.h" -#include "ortools/util/time_limit.h" namespace operations_research { namespace sat { @@ -46,9 +40,8 @@ namespace sat { // // The maximum energy is capacity * span of intervals at level 0. CutGenerator CreateCumulativeEnergyCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, - const std::vector>& energies, + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, const std::optional& makespan, Model* model); // For a given set of intervals and demands, we first compute the mandatory part @@ -63,22 +56,21 @@ CutGenerator CreateCumulativeEnergyCutGenerator( // sum(demands of always present intervals) // + sum(presence_literal * min_of_demand) <= capacity. CutGenerator CreateCumulativeTimeTableCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, Model* model); + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, Model* model); // Completion time cuts for the cumulative constraint. It is a simple relaxation // where we replace a cumulative task with demand k and duration d by a // no_overlap task with duration d * k / capacity_max. CutGenerator CreateCumulativeCompletionTimeCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, - const std::vector>& energies, Model* model); + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, Model* model); // For a given set of intervals in a cumulative constraint, we detect violated // mandatory precedences and create a cut for these. CutGenerator CreateCumulativePrecedenceCutGenerator( - SchedulingConstraintHelper* helper, const AffineExpression& capacity, - const std::vector& demands, Model* model); + SchedulingConstraintHelper* helper, SchedulingDemandHelper* demands_helper, + const AffineExpression& capacity, Model* model); // Completion time cuts for the no_overlap_2d constraint. It actually generates // the completion time cumulative cuts in both axis. diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index cffdaaa57c..130537b9d7 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -165,11 +165,11 @@ void SharedResponseManager::LogMessage(const std::string& prefix, void SharedResponseManager::LogPeriodicMessage(const std::string& prefix, const std::string& message, + double frequency_seconds, absl::Time* last_logging_time) { - const double freq = parameters_.log_frequency_in_seconds(); - if (freq <= 0.0 || last_logging_time == nullptr) return; + if (frequency_seconds < 0.0 || last_logging_time == nullptr) return; const absl::Time now = absl::Now(); - if (now - *last_logging_time < absl::Seconds(freq)) { + if (now - *last_logging_time < absl::Seconds(frequency_seconds)) { return; } diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 2b3ec1078f..bbcd4554df 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -335,6 +335,7 @@ class SharedResponseManager { void LogMessage(const std::string& prefix, const std::string& message); void LogPeriodicMessage(const std::string& prefix, const std::string& message, + double frequency_seconds, absl::Time* last_logging_time); bool LoggingIsEnabled() const { return logger_->LoggingIsEnabled(); } diff --git a/ortools/sat/timetable.cc b/ortools/sat/timetable.cc index 54ed95b198..e186c21c2f 100644 --- a/ortools/sat/timetable.cc +++ b/ortools/sat/timetable.cc @@ -281,14 +281,15 @@ bool ReservoirTimeTabling::TryToIncreaseMin(int event) { &literal_reason_, &integer_reason_); } -TimeTablingPerTask::TimeTablingPerTask( - const std::vector& demands, AffineExpression capacity, - IntegerTrail* integer_trail, SchedulingConstraintHelper* helper) +TimeTablingPerTask::TimeTablingPerTask(AffineExpression capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, + Model* model) : num_tasks_(helper->NumTasks()), - demands_(demands), capacity_(capacity), - integer_trail_(integer_trail), - helper_(helper) { + helper_(helper), + demands_(demands), + integer_trail_(model->GetOrCreate()) { // Each task may create at most two profile rectangles. Such pattern appear if // the profile is shaped like the Hanoi tower. The additional space is for // both extremities and the sentinels. @@ -312,7 +313,7 @@ void TimeTablingPerTask::RegisterWith(GenericLiteralWatcher* watcher) { helper_->WatchAllTasks(id, watcher); watcher->WatchUpperBound(capacity_.var, id); for (int t = 0; t < num_tasks_; t++) { - watcher->WatchLowerBound(demands_[t].var, id); + watcher->WatchLowerBound(demands_->Demands()[t], id); } watcher->RegisterReversibleInt(id, &num_profile_tasks_); @@ -396,14 +397,14 @@ bool TimeTablingPerTask::BuildProfile() { while (next_start >= 0 && by_decreasing_start_max[next_start].time == time) { const int t = by_decreasing_start_max[next_start].task_index; - if (IsInProfile(t)) current_height += DemandMin(t); + if (IsInProfile(t)) current_height += demands_->DemandMin(t); --next_start; } // Process the ending compulsory parts. while (next_end < num_tasks_ && by_end_min[next_end].time == time) { const int t = by_end_min[next_end].task_index; - if (IsInProfile(t)) current_height -= DemandMin(t); + if (IsInProfile(t)) current_height -= demands_->DemandMin(t); ++next_end; } @@ -448,7 +449,7 @@ bool TimeTablingPerTask::SweepAllTasks() { for (const auto& [t, time] : helper_->TaskByIncreasingStartMin()) { if (helper_->IsAbsent(t)) continue; if (helper_->SizeMin(t) == 0) continue; - if (DemandMin(t) <= demand_threshold) continue; + if (demands_->DemandMin(t) <= demand_threshold) continue; if (!SweepTask(t, time, &profile_index)) return false; } @@ -471,7 +472,8 @@ bool TimeTablingPerTask::SweepTask(int task_id, IntegerValue initial_start_min, // A profile rectangle is in conflict with the task if its height exceeds // conflict_height. - const IntegerValue conflict_height = CapacityMax() - DemandMin(task_id); + const IntegerValue conflict_height = + CapacityMax() - demands_->DemandMin(task_id); // Last time point during which task_id was in conflict with a profile // rectangle before being pushed. @@ -507,15 +509,16 @@ bool TimeTablingPerTask::SweepTask(int task_id, IntegerValue initial_start_min, // Since the task is part of the profile, try to lower its demand max // if possible. - const IntegerValue delta = DemandMax(task_id) - DemandMin(task_id); + const IntegerValue delta = + demands_->DemandMax(task_id) - demands_->DemandMin(task_id); if (delta > 0) { const IntegerValue threshold = CapacityMax() - delta; if (profile_[rec_id].start > start_max) --rec_id; for (; profile_[rec_id].start < initial_end_min; ++rec_id) { DCHECK_GT(profile_[rec_id + 1].start, start_max); if (profile_[rec_id].height <= threshold) continue; - const IntegerValue new_max = - CapacityMax() - profile_[rec_id].height + DemandMin(task_id); + const IntegerValue new_max = CapacityMax() - profile_[rec_id].height + + demands_->DemandMin(task_id); // Note that the task_id is already part of the profile reason, so // there is nothing else needed. @@ -523,7 +526,7 @@ bool TimeTablingPerTask::SweepTask(int task_id, IntegerValue initial_start_min, const IntegerValue time = std::max(start_max, profile_[rec_id].start); AddProfileReason(task_id, time, time + 1, CapacityMax()); if (!helper_->PushIntegerLiteralIfTaskPresent( - task_id, demands_[task_id].LowerOrEqual(new_max))) { + task_id, demands_->Demands()[task_id].LowerOrEqual(new_max))) { return false; } } @@ -561,7 +564,8 @@ bool TimeTablingPerTask::UpdateStartingTime(int task_id, IntegerValue left, IntegerValue right) { DCHECK_LT(left, right); helper_->ClearReason(); - AddProfileReason(task_id, left, right, CapacityMax() - DemandMin(task_id)); + AddProfileReason(task_id, left, right, + CapacityMax() - demands_->DemandMin(task_id)); if (capacity_.var != kNoIntegerVariable) { helper_->MutableIntegerReason()->push_back( integer_trail_->UpperBoundAsLiteral(capacity_.var)); @@ -570,10 +574,7 @@ bool TimeTablingPerTask::UpdateStartingTime(int task_id, IntegerValue left, // State of the task to be pushed. helper_->AddEndMinReason(task_id, left + 1); helper_->AddSizeMinReason(task_id); - if (demands_[task_id].var != kNoIntegerVariable) { - helper_->MutableIntegerReason()->push_back( - integer_trail_->LowerBoundAsLiteral(demands_[task_id].var)); - } + demands_->AddDemandMinReason(task_id); // Explain the increase of the minimum start and end times. return helper_->IncreaseStartMin(task_id, right); @@ -618,10 +619,7 @@ void TimeTablingPerTask::AddProfileReason(int task_id, IntegerValue left, // Note that we exclude the demand min for the task we push. // If we push the demand_max, we don't need it. And otherwise the task_id // shouldn't be part of the profice anyway. - if (demands_[t].var != kNoIntegerVariable && t != task_id) { - helper_->MutableIntegerReason()->push_back( - integer_trail_->LowerBoundAsLiteral(demands_[t].var)); - } + if (t != task_id) demands_->AddDemandMinReason(t); if (mode == 0) { helper_->AddStartMaxReason(t, left); @@ -632,7 +630,7 @@ void TimeTablingPerTask::AddProfileReason(int task_id, IntegerValue left, // TODO(user): Improve what task we "exclude" instead of always taking // the last ones? Note however that profile_tasks_ should be in order in // which task have a mandatory part. - sum_of_demand += integer_trail_->LowerBound(demands_[t]); + sum_of_demand += demands_->DemandMin(t); if (sum_of_demand > capacity_threshold) break; } else if (mode == 1) { helper_->AddStartMaxReason(t, start_max <= left ? left : right - 1); diff --git a/ortools/sat/timetable.h b/ortools/sat/timetable.h index 6674321b45..0f7dcc8546 100644 --- a/ortools/sat/timetable.h +++ b/ortools/sat/timetable.h @@ -104,11 +104,14 @@ class ReservoirTimeTabling : public PropagatorInterface { // A strongly quadratic version of Time Tabling filtering. This propagator // is similar to the CumulativeTimeTable propagator of the constraint solver. +// +// TODO(user): Use SchedulingDemandHelper. In particular, if we know the task +// is from a set of fixed alternatives, we might be able to push it more. class TimeTablingPerTask : public PropagatorInterface { public: - TimeTablingPerTask(const std::vector& demands, - AffineExpression capacity, IntegerTrail* integer_trail, - SchedulingConstraintHelper* helper); + TimeTablingPerTask(AffineExpression capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, Model* model); bool Propagate() final; @@ -171,14 +174,6 @@ class TimeTablingPerTask : public PropagatorInterface { return integer_trail_->UpperBound(capacity_); } - IntegerValue DemandMin(int task_id) const { - return integer_trail_->LowerBound(demands_[task_id]); - } - - IntegerValue DemandMax(int task_id) const { - return integer_trail_->UpperBound(demands_[task_id]); - } - // Returns true if the tasks is present and has a mantatory part. bool IsInProfile(int t) const { return positions_in_profile_tasks_[t] < num_profile_tasks_; @@ -187,14 +182,12 @@ class TimeTablingPerTask : public PropagatorInterface { // Number of tasks. const int num_tasks_; - // The demand variables of the tasks. - std::vector demands_; - // Capacity of the resource. const AffineExpression capacity_; - IntegerTrail* integer_trail_; SchedulingConstraintHelper* helper_; + SchedulingDemandHelper* demands_; + IntegerTrail* integer_trail_; // Optimistic profile of the resource consumption over time. std::vector profile_; diff --git a/ortools/sat/timetable_edgefinding.cc b/ortools/sat/timetable_edgefinding.cc index 69e0ae94d7..10624f653d 100644 --- a/ortools/sat/timetable_edgefinding.cc +++ b/ortools/sat/timetable_edgefinding.cc @@ -27,15 +27,15 @@ namespace operations_research { namespace sat { -TimeTableEdgeFinding::TimeTableEdgeFinding( - const std::vector& demands, AffineExpression capacity, - SchedulingConstraintHelper* helper, IntegerTrail* integer_trail, - Model* model) +TimeTableEdgeFinding::TimeTableEdgeFinding(AffineExpression capacity, + SchedulingConstraintHelper* helper, + SchedulingDemandHelper* demands, + Model* model) : num_tasks_(helper->NumTasks()), - demands_(demands), capacity_(capacity), helper_(helper), - integer_trail_(integer_trail) { + demands_(demands), + integer_trail_(model->GetOrCreate()) { // Edge finding structures. mandatory_energy_before_end_max_.resize(num_tasks_); mandatory_energy_before_start_min_.resize(num_tasks_); @@ -43,21 +43,16 @@ TimeTableEdgeFinding::TimeTableEdgeFinding( // Energy of free parts. size_free_.resize(num_tasks_); energy_free_.resize(num_tasks_); - - // Try to linearize the energy. - energies_.resize(num_tasks_); - for (int t = 0; t < num_tasks_; t++) { - energies_[t] = TryToLinearizeProduct(helper->Sizes()[t], demands[t], model); - } } void TimeTableEdgeFinding::RegisterWith(GenericLiteralWatcher* watcher) { const int id = watcher->Register(this); - watcher->WatchUpperBound(capacity_.var, id); + watcher->WatchUpperBound(capacity_, id); helper_->WatchAllTasks(id, watcher); for (int t = 0; t < num_tasks_; t++) { - watcher->WatchLowerBound(demands_[t].var, id); + watcher->WatchLowerBound(demands_->Demands()[t], id); } + watcher->SetPropagatorPriority(id, 3); watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); } @@ -113,7 +108,6 @@ void TimeTableEdgeFinding::BuildTimeTable() { while (index_emax >= 0) { // Next time point. - // TODO(user): could be simplified with a sentinel. IntegerValue time = by_decreasing_end_max[index_emax].time; if (index_smin < num_tasks_) { time = std::min(time, by_start_min[index_smin].time); @@ -145,48 +139,46 @@ void TimeTableEdgeFinding::BuildTimeTable() { // Process the starting compulsory parts. while (index_scp < scp_.size() && scp_[index_scp].time == time) { - height += DemandMin(scp_[index_scp].task_index); + height += demands_->DemandMin(scp_[index_scp].task_index); index_scp++; } // Process the ending compulsory parts. while (index_ecp < ecp_.size() && ecp_[index_ecp].time == time) { - height -= DemandMin(ecp_[index_ecp].task_index); + height -= demands_->DemandMin(ecp_[index_ecp].task_index); index_ecp++; } } } bool TimeTableEdgeFinding::TimeTableEdgeFindingPass() { + demands_->CacheAllEnergyValues(); + // Initialize the data structures and build the free parts. // -------------------------------------------------------- for (int t = 0; t < num_tasks_; ++t) { // If the task has no mandatory part, then its free part is the task itself. const IntegerValue start_max = helper_->StartMax(t); const IntegerValue end_min = helper_->EndMin(t); - const IntegerValue demand_min = DemandMin(t); + const IntegerValue demand_min = demands_->DemandMin(t); IntegerValue mandatory_energy(0); if (start_max >= end_min) { size_free_[t] = helper_->SizeMin(t); } else { - size_free_[t] = helper_->SizeMin(t) + start_max - end_min; - mandatory_energy = (end_min - start_max) * demand_min; + const IntegerValue mandatory_size = end_min - start_max; + size_free_[t] = helper_->SizeMin(t) - mandatory_size; + mandatory_energy = mandatory_size * demand_min; } - const IntegerValue min_energy_of_product = helper_->SizeMin(t) * demand_min; - const IntegerValue min_linear_energy = - energies_[t].has_value() ? energies_[t].value().Min(*integer_trail_) - : IntegerValue(0); - CHECK_EQ(size_free_[t] * DemandMin(t) + mandatory_energy, - min_energy_of_product); - // if (min_linear_energy > min_energy_of_product) { - // LOG(INFO) << "linear: " << min_linear_energy.value() - // << " vs product: " << min_energy_of_product.value(); - // } - energy_free_[t] = - std::max(min_linear_energy, min_energy_of_product) - mandatory_energy; + + const IntegerValue min_energy = demands_->EnergyMin(t); + energy_free_[t] = min_energy - mandatory_energy; + DCHECK_GE(energy_free_[t], 0); } + // TODO(user): Is it possible to have a 'higher' mandatory profile using + // the min energy instead of the demand_min * size_min? How can we incorporate + // this extra energy in the mandatory profile ? BuildTimeTable(); const auto& by_start_min = helper_->TaskByIncreasingStartMin(); @@ -206,39 +198,43 @@ bool TimeTableEdgeFinding::TimeTableEdgeFindingPass() { if (end_task_time.time == previous_end) continue; previous_end = end_task_time.time; - // Energy of the free parts contained in the interval [begin, end). + // Energy of the free parts contained in the interval + // [window_min, window_max]. IntegerValue energy_free_parts = IntegerValue(0); + reason_tasks_fully_included_in_window_.clear(); + reason_tasks_partially_included_in_window_.clear(); // Task that requires the biggest additional amount of energy to be - // scheduled at its minimum start time in the task interval [begin, end). + // scheduled at its minimum start time in the task interval + // [window_min, window_max]. int max_task = -1; IntegerValue free_energy_of_max_task_in_window(0); IntegerValue extra_energy_required_by_max_task = kMinIntegerValue; // Process task by decreasing start min. + const IntegerValue window_max = end_task_time.time; for (const TaskTime begin_task_time : gtl::reversed_view(by_start_min)) { const int begin_task = begin_task_time.task_index; + // The considered time window. Note that we use the "cached" values so + // that our mandatory energy before computation is correct. + const IntegerValue window_min = begin_task_time.time; + + // Not a valid time window. + if (window_max <= window_min) continue; + // TODO(user): consider optional tasks for additional propagation. if (!helper_->IsPresent(begin_task)) continue; if (energy_free_[begin_task] == 0) continue; - // The considered time window. Note that we use the "cached" values so - // that our mandatory energy before computation is correct. - const IntegerValue begin = begin_task_time.time; // Start min. - const IntegerValue end = end_task_time.time; // End max. - - // Not a valid time window. - if (end <= begin) continue; - // We consider two different cases: either the free part overlaps the - // end of the interval (right) or it does not (inside). + // window_max of the interval (right) or it does not (inside). // - // begin end + // window_min window_max // v v // right: ======|=== // - // begin end + // window_min window_max // v v // inside: ========== | // @@ -248,22 +244,30 @@ bool TimeTableEdgeFinding::TimeTableEdgeFindingPass() { // equal to the largest part of the free part that can fit in the task // interval. const IntegerValue end_max = helper_->EndMax(begin_task); - if (end_max <= end) { - // The whole task energy is contained in the task interval. + if (end_max <= window_max) { + // The whole task energy is contained in the window. + reason_tasks_fully_included_in_window_.push_back(begin_task); energy_free_parts += energy_free_[begin_task]; } else { - const IntegerValue demand_min = DemandMin(begin_task); + const IntegerValue demand_min = demands_->DemandMin(begin_task); const IntegerValue extra_energy = - std::min(size_free_[begin_task], (end - begin)) * demand_min; + std::min(size_free_[begin_task], (window_max - window_min)) * + demand_min; // This is not in the paper, but it is almost free for us to account for // the free energy of this task that must be present in the window. const IntegerValue free_energy_in_window = std::max(IntegerValue(0), - size_free_[begin_task] - (end_max - end)) * + size_free_[begin_task] - (end_max - window_max)) * demand_min; + // TODO(user): There is no point setting max_task if its start min + // is already bigger that what we can push. Maybe we can exploit that? if (extra_energy > extra_energy_required_by_max_task) { + if (max_task != -1 && free_energy_of_max_task_in_window > 0) { + reason_tasks_partially_included_in_window_.push_back(max_task); + } + max_task = begin_task; extra_energy_required_by_max_task = extra_energy; @@ -271,7 +275,8 @@ bool TimeTableEdgeFinding::TimeTableEdgeFindingPass() { // new one for later. energy_free_parts += free_energy_of_max_task_in_window; free_energy_of_max_task_in_window = free_energy_in_window; - } else { + } else if (free_energy_in_window > 0) { + reason_tasks_partially_included_in_window_.push_back(begin_task); energy_free_parts += free_energy_in_window; } } @@ -283,15 +288,34 @@ bool TimeTableEdgeFinding::TimeTableEdgeFindingPass() { if (max_task == -1) continue; // Compute the amount of energy available to schedule max_task. - const IntegerValue interval_energy = CapacityMax() * (end - begin); + const IntegerValue window_energy = + CapacityMax() * (window_max - window_min); const IntegerValue energy_mandatory = mandatory_energy_before_end_max_[end_task] - mandatory_energy_before_start_min_[begin_task]; const IntegerValue available_energy = - interval_energy - energy_free_parts - energy_mandatory; + window_energy - energy_free_parts - energy_mandatory; - // Enough energy to schedule max_task at its minimum start time. - if (extra_energy_required_by_max_task <= available_energy) continue; + // Enough energy to schedule max_task at its minimum start time? + // + // TODO(user, fdid): In case of alternatives, for each fixed + // size/demand pair, we can compute a new_start and use the min of them. + if (extra_energy_required_by_max_task <= available_energy) { + // If the test below is true, we know the max_task cannot fully + // fit in the time window, so at least end_min > window_max. + // + // TODO(user): We currently only do that if we are not about to push the + // start as we assume the start push is just stronger. Maybe we should + // do it in more situation? + if (energy_free_[max_task] > available_energy && + helper_->EndMin(max_task) <= window_max) { + FillEnergyInWindowReason(window_min, window_max, max_task); + demands_->AddEnergyMinReason(max_task); + helper_->AddStartMinReason(max_task, window_min); + if (!helper_->IncreaseEndMin(max_task, window_max + 1)) return false; + } + continue; + } // Compute the length of the mandatory subpart of max_task that should be // considered as available. @@ -299,17 +323,27 @@ bool TimeTableEdgeFinding::TimeTableEdgeFindingPass() { // TODO(user): Because this use updated bounds, it might be more than what // we accounted for in the precomputation. This is correct but could be // improved uppon. - const IntegerValue mandatory_in = std::max( - IntegerValue(0), std::min(end, helper_->EndMin(max_task)) - - std::max(begin, helper_->StartMax(max_task))); + const IntegerValue mandatory_size_in_window = + std::max(IntegerValue(0), + std::min(window_max, helper_->EndMin(max_task)) - + std::max(window_min, helper_->StartMax(max_task))); // Compute the new minimum start time of max_task. + const IntegerValue max_free_size_that_fit = + available_energy / demands_->DemandMin(max_task); const IntegerValue new_start = - end - mandatory_in - (available_energy / DemandMin(max_task)); + window_max - mandatory_size_in_window - max_free_size_that_fit; // Push and explain only if the new start is bigger than the current one. if (helper_->StartMin(max_task) < new_start) { - if (!IncreaseStartMin(begin, end, max_task, new_start)) return false; + FillEnergyInWindowReason(window_min, window_max, max_task); + + // Reason needed for task_index. + // We only need start_min and demand_min to push the start. + helper_->AddStartMinReason(max_task, window_min); + demands_->AddDemandMinReason(max_task); + + if (!helper_->IncreaseStartMin(max_task, new_start)) return false; } } } @@ -317,61 +351,62 @@ bool TimeTableEdgeFinding::TimeTableEdgeFindingPass() { return true; } -void TimeTableEdgeFinding::AddEnergyReason(int task_index) { - if (energies_[task_index].has_value()) { - for (const IntegerVariable var : energies_[task_index].value().vars) { - helper_->MutableIntegerReason()->push_back( - integer_trail_->LowerBoundAsLiteral(var)); - } - } else { - // Variables of the task to be pushed. We do not need the end max for this - // task and we only need for it to begin in the time window. - if (demands_[task_index].var != kNoIntegerVariable) { - helper_->MutableIntegerReason()->push_back( - integer_trail_->LowerBoundAsLiteral(demands_[task_index].var)); - } - helper_->AddSizeMinReason(task_index); - } -} - -bool TimeTableEdgeFinding::IncreaseStartMin(IntegerValue begin, - IntegerValue end, int task_index, - IntegerValue new_start) { +void TimeTableEdgeFinding::FillEnergyInWindowReason(IntegerValue window_min, + IntegerValue window_max, + int task_index) { helper_->ClearReason(); - std::vector* mutable_reason = helper_->MutableIntegerReason(); // Capacity of the resource. if (capacity_.var != kNoIntegerVariable) { - mutable_reason->push_back( + helper_->MutableIntegerReason()->push_back( integer_trail_->UpperBoundAsLiteral(capacity_.var)); } - helper_->AddStartMinReason(task_index, begin); - AddEnergyReason(task_index); - - // Task contributing to the energy in the interval. + // Tasks contributing to the mandatory energy in the interval. for (int t = 0; t < num_tasks_; ++t) { if (t == task_index) continue; if (!helper_->IsPresent(t)) continue; - if (helper_->EndMax(t) <= begin) continue; - if (helper_->StartMin(t) >= end) continue; - - // We need the reason for the energy contribution of this interval into - // [begin, end]. - // - // TODO(user): Since we actually do not account fully for this energy, we - // could relax the reason more. - // - // TODO(user): This reason might not be enough in the presence of variable - // size intervals where StartMax and EndMin give rise to more energy - // that just using size min and these bounds. Fix. - helper_->AddStartMinReason(t, std::min(begin, helper_->StartMin(t))); - helper_->AddEndMaxReason(t, std::max(end, helper_->EndMax(t))); + const IntegerValue smax = helper_->StartMax(t); + const IntegerValue emin = helper_->EndMin(t); + if (smax >= emin) continue; + if (emin <= window_min) continue; + if (smax >= window_max) continue; + helper_->AddStartMaxReason(t, std::max(smax, window_min)); + helper_->AddEndMinReason(t, std::min(emin, window_max)); helper_->AddPresenceReason(t); - AddEnergyReason(t); + demands_->AddDemandMinReason(t); } - return helper_->IncreaseStartMin(task_index, new_start); + // Tasks contributing to the free energy in [window_min, window_max]. + // + // TODO(user): If a task appears in both, we could avoid adding twice the + // same things, but the core solver should merge duplicates anyway. + for (const int t : reason_tasks_fully_included_in_window_) { + DCHECK_NE(t, task_index); + DCHECK(helper_->IsPresent(t)); + DCHECK_GT(helper_->EndMax(t), window_min); + DCHECK_LT(helper_->StartMin(t), window_max); + DCHECK_GE(helper_->StartMin(t), window_min); + + helper_->AddStartMinReason(t, helper_->StartMin(t)); + helper_->AddEndMaxReason(t, std::max(window_max, helper_->EndMax(t))); + helper_->AddPresenceReason(t); + demands_->AddEnergyMinReason(t); + } + for (const int t : reason_tasks_partially_included_in_window_) { + DCHECK_NE(t, task_index); + DCHECK(helper_->IsPresent(t)); + DCHECK_GT(helper_->EndMax(t), window_min); + DCHECK_LT(helper_->StartMin(t), window_max); + DCHECK_GE(helper_->StartMin(t), window_min); + + helper_->AddStartMinReason(t, helper_->StartMin(t)); + helper_->AddEndMaxReason(t, std::max(window_max, helper_->EndMax(t))); + helper_->AddPresenceReason(t); + + helper_->AddSizeMinReason(t); + demands_->AddDemandMinReason(t); + } } } // namespace sat diff --git a/ortools/sat/timetable_edgefinding.h b/ortools/sat/timetable_edgefinding.h index 45392654c3..00a26c47cf 100644 --- a/ortools/sat/timetable_edgefinding.h +++ b/ortools/sat/timetable_edgefinding.h @@ -60,10 +60,9 @@ namespace sat { // time if this would cause an overload in one of the task intervals. class TimeTableEdgeFinding : public PropagatorInterface { public: - TimeTableEdgeFinding(const std::vector& demands, - AffineExpression capacity, + TimeTableEdgeFinding(AffineExpression capacity, SchedulingConstraintHelper* helper, - IntegerTrail* integer_trail, Model* model); + SchedulingDemandHelper* demands, Model* model); bool Propagate() final; @@ -83,40 +82,27 @@ class TimeTableEdgeFinding : public PropagatorInterface { // update the end times by calling the SwitchToMirrorProblem method first. bool TimeTableEdgeFindingPass(); - // Increases the start min of task_index with the proper explanation. - bool IncreaseStartMin(IntegerValue begin, IntegerValue end, int task_index, - IntegerValue new_start); - - // Add the minimum energy reason of a task. - void AddEnergyReason(int task_index); - - IntegerValue DemandMin(int task_index) const { - return integer_trail_->LowerBound(demands_[task_index]); - } + // Fills the reason for the energy in [window_min, window_max]. + // We exclude the given task_index mandatory energy and uses + // tasks_contributing_to_free_energy_. + void FillEnergyInWindowReason(IntegerValue window_min, + IntegerValue window_max, int task_index); IntegerValue CapacityMax() const { return integer_trail_->UpperBound(capacity_); } - // Number of tasks. const int num_tasks_; - - // IntervalVariable and IntegerVariable of each tasks that must be considered - // in this constraint. - std::vector demands_; const AffineExpression capacity_; - SchedulingConstraintHelper* helper_; + SchedulingDemandHelper* demands_; IntegerTrail* integer_trail_; - std::vector> energies_; - // Start (resp. end) of the compulsory parts used to build the profile. std::vector scp_; std::vector ecp_; - // Sizes and energy of the free parts. One is just the other times the - // minimum demand. + // Sizes and energy of the free parts. std::vector size_free_; std::vector energy_free_; @@ -125,6 +111,10 @@ class TimeTableEdgeFinding : public PropagatorInterface { std::vector mandatory_energy_before_start_min_; std::vector mandatory_energy_before_end_max_; + // List of task that should participate in the reason. + std::vector reason_tasks_fully_included_in_window_; + std::vector reason_tasks_partially_included_in_window_; + DISALLOW_COPY_AND_ASSIGN(TimeTableEdgeFinding); };