diff --git a/ortools/sat/disjunctive.cc b/ortools/sat/disjunctive.cc index 3208ab56c0..808ca7f18c 100644 --- a/ortools/sat/disjunctive.cc +++ b/ortools/sat/disjunctive.cc @@ -202,8 +202,7 @@ bool DisjunctiveOverloadChecker::Propagate() { // Set up theta tree. start_event_to_task_.clear(); start_event_time_.clear(); - start_event_duration_.clear(); - int num_events = 0; + int num_events_ = 0; for (const auto task_time : helper_->TaskByIncreasingStartMin()) { const int task = task_time.task_index; // TODO(user): Add max energy deduction for variable durations. @@ -214,23 +213,30 @@ bool DisjunctiveOverloadChecker::Propagate() { } start_event_to_task_.push_back(task); start_event_time_.push_back(task_time.time); - start_event_duration_.push_back(helper_->DurationMin(task)); - task_to_start_event_[task] = num_events; - num_events++; + task_to_start_event_[task] = num_events_; + num_events_++; } - theta_tree_.Reset(start_event_time_, start_event_duration_); + const int num_events = num_events_; + start_event_is_present_.assign(num_events, false); + theta_tree_.Reset(num_events); - // Overload checker loop. - const auto& task_by_increasing_max_end = - ::gtl::reversed_view(helper_->TaskByDecreasingEndMax()); - for (const auto task_time : task_by_increasing_max_end) { + // Introduce events by nondecreasing end_max, check for overloads. + for (const auto task_time : + ::gtl::reversed_view(helper_->TaskByDecreasingEndMax())) { const int current_task = task_time.task_index; - if (task_to_start_event_[current_task] == -1) { - continue; - } else if (helper_->IsPresent(current_task)) { - theta_tree_.AddEvent(task_to_start_event_[current_task]); - } else { - theta_tree_.AddOptionalEvent(task_to_start_event_[current_task]); + if (task_to_start_event_[current_task] == -1) continue; + + { + const int current_event = task_to_start_event_[current_task]; + const bool is_present = helper_->IsPresent(current_task); + start_event_is_present_[current_event] = is_present; + // TODO(user): consider reducing max available duration. + const IntegerValue energy_max = helper_->DurationMin(current_task); + const IntegerValue energy_min = is_present ? energy_max : IntegerValue(0); + + theta_tree_.AddOrUpdateEvent(current_event, + start_event_time_[current_event], energy_min, + energy_max); } const IntegerValue current_end = task_time.time; @@ -238,14 +244,13 @@ bool DisjunctiveOverloadChecker::Propagate() { // Explain failure with tasks in critical interval. helper_->ClearReason(); const int critical_event = - theta_tree_.GetRightmostEventWithEnvelopeGreaterOrEqualTo( - current_end + 1); + theta_tree_.GetMaxEventWithEnvelopeGreaterThan(current_end); const IntegerValue window_start = start_event_time_[critical_event]; const IntegerValue window_end = theta_tree_.GetEnvelopeOf(critical_event) - 1; for (int event = critical_event; event < num_events; event++) { - if (theta_tree_.IsPresent(event)) { + if (start_event_is_present_[event]) { const int task = start_event_to_task_[event]; helper_->AddPresenceReason(task); helper_->AddDurationMinReason(task); @@ -262,13 +267,20 @@ bool DisjunctiveOverloadChecker::Propagate() { // TODO(user): This could be done lazily, like most of the loop to // compute the reasons in this file. helper_->ClearReason(); - const int critical_event = - theta_tree_.GetEventRealizingOptionalEnvelope(); + int critical_event; + int optional_event; + IntegerValue unused; + theta_tree_.GetEventsWithOptionalEnvelopeGreaterThan( + current_end, &critical_event, &optional_event, &unused); + const IntegerValue window_start = start_event_time_[critical_event]; - const IntegerValue window_end = theta_tree_.GetOptionalEnvelope() - 1; + const int optional_task = start_event_to_task_[optional_event]; + const IntegerValue window_end = + theta_tree_.GetEnvelopeOf(critical_event) + + helper_->DurationMin(optional_task) - 1; for (int event = critical_event; event < num_events; event++) { - if (theta_tree_.IsPresent(event)) { + if (start_event_is_present_[event]) { const int task = start_event_to_task_[event]; helper_->AddPresenceReason(task); helper_->AddDurationMinReason(task); @@ -277,9 +289,6 @@ bool DisjunctiveOverloadChecker::Propagate() { } } - const int optional_event = - theta_tree_.GetOptionalEventRealizingOptionalEnvelope(); - const int optional_task = start_event_to_task_[optional_event]; helper_->AddDurationMinReason(optional_task); helper_->AddStartMinReason(optional_task, window_start); helper_->AddEndMaxReason(optional_task, window_end); diff --git a/ortools/sat/disjunctive.h b/ortools/sat/disjunctive.h index 1df71ecc75..636fe5ef49 100644 --- a/ortools/sat/disjunctive.h +++ b/ortools/sat/disjunctive.h @@ -136,7 +136,7 @@ class DisjunctiveOverloadChecker : public PropagatorInterface { std::vector task_to_start_event_; std::vector start_event_to_task_; std::vector start_event_time_; - std::vector start_event_duration_; + std::vector start_event_is_present_; }; class DisjunctiveDetectablePrecedences : public PropagatorInterface { diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index 8a2e791251..8ff54b9609 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -1216,13 +1216,16 @@ SatSolver::Status SolveIntegerProblemWithLazyEncoding( time_limit = model->Mutable(); } SatSolver* const solver = model->GetOrCreate(); + if (solver->IsModelUnsat()) return SatSolver::MODEL_UNSAT; + solver->Backtrack(0); + solver->SetAssumptionLevel(0); + const SatParameters parameters = solver->parameters(); if (parameters.use_fixed_search()) { CHECK(assumptions.empty()) << "Not supported yet"; // Note that it is important to do the level-zero propagation if it wasn't // already done because EnqueueDecisionAndBackjumpOnConflict() assumes that // the solver is in a "propagated" state. - solver->Backtrack(0); if (!solver->Propagate()) return SatSolver::MODEL_UNSAT; } while (!time_limit->LimitReached()) { @@ -1230,9 +1233,6 @@ SatSolver::Status SolveIntegerProblemWithLazyEncoding( // instanciate all the already created Booleans. if (!parameters.use_fixed_search()) { if (assumptions.empty()) { - // TODO(user): This one doesn't do Backtrack(0), and doing it seems to - // trigger a bug in research/devtools/compote/scheduler - // :instruction_scheduler_test, investigate. const SatSolver::Status status = solver->SolveWithTimeLimit(time_limit); if (status != SatSolver::MODEL_SAT) return status; } else { diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index 1a90fbf5e8..29896deec0 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -117,7 +117,7 @@ struct IntegerLiteral { return var == o.var && bound == o.bound; } bool operator!=(IntegerLiteral o) const { - return var != o.var || bound == o.bound; + return var != o.var || bound != o.bound; } IntegerVariable Var() const { return IntegerVariable(var); } @@ -139,10 +139,6 @@ struct IntegerLiteral { // Our external API uses IntegerVariable but internally we only use an int for // simplicity. TODO(user): change this? // - // TODO(user): We can't use const because we want to be able to copy a - // std::vector. So instead make them private and provide some - // getters. - // // Note that bound is always in [kMinIntegerValue, kMaxIntegerValue + 1]. /*const*/ int var; /*const*/ IntegerValue bound; @@ -433,6 +429,11 @@ class IntegerTrail : public SatPropagator { IntegerValue LowerBound(IntegerVariable i) const; IntegerValue UpperBound(IntegerVariable i) const; + // Returns the value of the lower bound before the last Enqueue() that changed + // it. Note that PreviousLowerBound() == LowerBound() iff this is the level + // zero bound. + IntegerValue PreviousLowerBound(IntegerVariable i) const; + // Returns the integer literal that represent the current lower/upper bound of // the given integer variable. IntegerLiteral LowerBoundAsLiteral(IntegerVariable i) const; @@ -823,6 +824,12 @@ inline IntegerValue IntegerTrail::LowerBound(IntegerVariable i) const { return vars_[i.value()].current_bound; } +inline IntegerValue IntegerTrail::PreviousLowerBound(IntegerVariable i) const { + const int index = vars_[i.value()].current_trail_index; + if (index < vars_.size()) return LowerBound(i); + return integer_trail_[integer_trail_[index].prev_trail_index].bound; +} + inline IntegerValue IntegerTrail::UpperBound(IntegerVariable i) const { return -vars_[NegationOf(i).value()].current_bound; } diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index ff76e18ab8..b10986b178 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -94,9 +94,27 @@ bool IntegerSumLE::Propagate() { const IntegerValue new_lb = rev_lb_fixed_vars_ + lb_unfixed_vars; // Conflict? - const IntegerValue slack = upper_bound_ - new_lb; + IntegerValue slack = upper_bound_ - new_lb; if (slack < 0) { - FillIntegerReason(); + // Like FillIntegerReason() but try to relax the reason a bit. + // + // TODO(user): if not all the slack is consumed, we could relax it even + // more. It might also be advantageous to relax first the variable with the + // highest "trail index". + integer_reason_.clear(); + for (int i = 0; i < vars_.size(); ++i) { + const IntegerVariable var = vars_[i]; + const IntegerValue lb = integer_trail_->LowerBound(var); + const IntegerValue prev_lb = integer_trail_->PreviousLowerBound(var); + if (lb == prev_lb) continue; // level zero. + const IntegerValue diff = (lb - prev_lb) * coeffs_[i]; + if (slack + diff < 0) { + integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, prev_lb)); + slack += diff; + } else { + integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, lb)); + } + } // Reified case: If the reified literal is unassigned, we set it to false, // otherwise we have a conflict. @@ -132,9 +150,7 @@ bool IntegerSumLE::Propagate() { FillIntegerReason(); } - // We need to remove the entry index from the reason temporarily. This is - // OK with the reversible aspect of integer_reason_ because fixed - // variables will never be swapped. + // We need to remove the entry index from the reason temporarily. IntegerLiteral saved; const int index = index_in_integer_reason_[i]; if (index >= 0) { diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index 236f0f0516..5c6a7f96ff 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -19,6 +19,7 @@ #include "ortools/base/stringprintf.h" #include "google/protobuf/descriptor.h" #include "ortools/sat/encoding.h" +#include "ortools/sat/integer_expr.h" #include "ortools/sat/util.h" namespace operations_research { @@ -1107,6 +1108,278 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( return result; } +SatSolver::Status MinimizeWithCoreAndLazyEncoding( + bool log_info, IntegerVariable objective_var, + const std::vector& variables, + const std::vector& coefficients, + const std::function& next_decision, + const std::function& feasible_solution_observer, + Model* model) { + SatSolver* sat_solver = model->GetOrCreate(); + IntegerTrail* integer_trail = model->GetOrCreate(); + IntegerEncoder* integer_encoder = model->GetOrCreate(); + + // This will be called each time a feasible solution is found. Returns false + // if a conflict was detected while trying to constrain the objective to a + // smaller value. + int num_solutions = 0; + IntegerValue best_objective = integer_trail->UpperBound(objective_var); + const auto process_solution = [&]() { + const IntegerValue objective(model->Get(Value(objective_var))); + if (objective >= best_objective) return true; + + ++num_solutions; + best_objective = objective; + if (feasible_solution_observer != nullptr) { + feasible_solution_observer(*model); + } + + // TODO(user): Investigate if constraining the objective is better. + if (/* DISALBES CODE */ (false)) { + sat_solver->Backtrack(0); + sat_solver->SetAssumptionLevel(0); + if (!integer_trail->Enqueue( + IntegerLiteral::LowerOrEqual(objective_var, objective - 1), {}, + {})) { + return false; + } + } + return true; + }; + + // We express the objective as a linear sum of terms. These will evolve as the + // algorithm progress. + struct ObjectiveTerm { + IntegerVariable var; + IntegerValue weight; + + // These fields are only used for logging/debugging. + int depth; + IntegerValue old_var_lb; + }; + std::vector terms; + CHECK_EQ(variables.size(), coefficients.size()); + for (int i = 0; i < variables.size(); ++i) { + if (coefficients[i] > 0) { + terms.push_back({variables[i], coefficients[i]}); + } else if (coefficients[i] < 0) { + terms.push_back({NegationOf(variables[i]), -coefficients[i]}); + } + terms.back().depth = 0; + } + + // This is used by the "stratified" approach. We will only consider terms with + // a weight not lower than this threshold. The threshold will decrease as the + // algorithm progress. + IntegerValue stratified_threshold = kMaxIntegerValue; + + // TODO(user): The core is returned in the same order as the assumptions, + // so we don't really need this map, we could just do a linear scan to + // recover which node are part of the core. + std::map assumption_to_term_index; + + // Start the algorithm. + int max_depth = 0; + SatSolver::Status result; + for (int iter = 0;; ++iter) { + sat_solver->Backtrack(0); + sat_solver->SetAssumptionLevel(0); + + // We assumes all terms at their lower-bound. + std::vector assumptions; + assumption_to_term_index.clear(); + IntegerValue next_stratified_threshold(0); + IntegerValue implied_objective_lb; + for (int i = 0; i < terms.size(); ++i) { + const ObjectiveTerm term = terms[i]; + const IntegerValue var_lb = integer_trail->LowerBound(term.var); + terms[i].old_var_lb = var_lb; + implied_objective_lb += term.weight * var_lb.value(); + + // TODO(user): These can be simply removed from the list. + if (term.weight == 0) continue; + + // Skip fixed terms. + // We still keep them around for a proper lower-bound computation. + // TODO(user): we could keep an objective offset instead. + if (var_lb == integer_trail->UpperBound(term.var)) continue; + + // Only consider the terms above the threshold. + if (term.weight < stratified_threshold) { + next_stratified_threshold = + std::max(next_stratified_threshold, term.weight); + } else { + assumptions.push_back(integer_encoder->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(term.var, var_lb))); + InsertOrDie(&assumption_to_term_index, assumptions.back().Index(), i); + } + } + + // Update the objective lower bound with our current bound. + // + // Note(user): This is not needed for correctness, but it might cause + // more propagation and is nice to have for reporting/logging purpose. + if (!integer_trail->Enqueue( + IntegerLiteral::GreaterOrEqual(objective_var, implied_objective_lb), + {}, {})) { + result = SatSolver::MODEL_UNSAT; + break; + } + + // No assumptions with the current stratified_threshold? use the new one. + if (assumptions.empty() && next_stratified_threshold > 0) { + stratified_threshold = next_stratified_threshold; + --iter; // "false" iteration, the lower bound does not increase. + continue; + } + + // Display the progress. + const IntegerValue objective_lb = integer_trail->LowerBound(objective_var); + if (log_info) { + LOG(INFO) << StringPrintf( + " iter:%d lb:%lld (%lld) gap:%lld assumptions:%zu strat:%lld " + "depth:%d", + iter, objective_lb.value(), implied_objective_lb.value(), + (best_objective - objective_lb).value(), assumptions.size(), + stratified_threshold.value(), max_depth); + } + + // Solve under the assumptions. + result = + SolveIntegerProblemWithLazyEncoding(assumptions, next_decision, model); + if (result == SatSolver::MODEL_SAT) { + if (!process_solution()) { + result = SatSolver::MODEL_UNSAT; + break; + } + + // If not all assumptions where taken, continue with a lower stratified + // bound. Otherwise we have an optimal solution. + stratified_threshold = next_stratified_threshold; + if (stratified_threshold == 0) break; + --iter; // "false" iteration, the lower bound does not increase. + continue; + } + if (result != SatSolver::ASSUMPTIONS_UNSAT) break; + + // We have a new core. + std::vector core = sat_solver->GetLastIncompatibleDecisions(); + if (sat_solver->parameters().minimize_core()) { + MinimizeCore(sat_solver, &core); + } + CHECK(!core.empty()); + + // This just increase the lower-bound of the corresponding node, which + // should already be done by the solver. + if (core.size() == 1) continue; + + sat_solver->Backtrack(0); + sat_solver->SetAssumptionLevel(0); + + // Compute the min weight of all the terms in the core. The lower bound will + // be increased by that much because at least one assumption in the core + // must be true. This is also why we can start at 1 for new_var_lb. + IntegerValue min_weight = kMaxIntegerValue; + IntegerValue max_weight(0); + IntegerValue new_var_lb(1); + IntegerValue new_var_ub(0); + int new_depth = 0; + for (const Literal lit : core) { + const int index = FindOrDie(assumption_to_term_index, lit.Index()); + min_weight = std::min(min_weight, terms[index].weight); + max_weight = std::max(max_weight, terms[index].weight); + new_depth = std::max(new_depth, terms[index].depth + 1); + new_var_lb += integer_trail->LowerBound(terms[index].var); + new_var_ub += integer_trail->UpperBound(terms[index].var); + CHECK_EQ(terms[index].old_var_lb, + integer_trail->LowerBound(terms[index].var)); + } + max_depth = std::max(max_depth, new_depth); + if (log_info) { + LOG(INFO) << StringPrintf( + " core:%zu weight:[%lld,%lld] domain:[%lld,%lld] depth:%d", + core.size(), min_weight.value(), max_weight.value(), + new_var_lb.value(), new_var_ub.value(), new_depth); + } + + // We will "transfer" min_weight from all the variables of the core + // to a new variable. + const IntegerVariable new_var = + model->Add(NewIntegerVariable(new_var_lb.value(), new_var_ub.value())); + terms.push_back({new_var, min_weight, new_depth}); + + // Sum variables in the core <= new_var. + // TODO(user): Experiment with FixedWeightedSum() instead. + { + std::vector constraint_vars; + std::vector constraint_coeffs; + for (const Literal lit : core) { + const int index = FindOrDie(assumption_to_term_index, lit.Index()); + terms[index].weight -= min_weight; + constraint_vars.push_back(terms[index].var); + constraint_coeffs.push_back(1); + } + constraint_vars.push_back(new_var); + constraint_coeffs.push_back(-1); + model->Add( + WeightedSumLowerOrEqual(constraint_vars, constraint_coeffs, 0)); + } + + // Re-express the objective with the new terms. + // TODO(user): Do more experiments to decide if this is better. + // TODO(user): Experiment with FixedWeightedSum(). + if (/* DISALBES CODE */ (false)) { + std::vector constraint_vars; + std::vector constraint_coeffs; + for (const ObjectiveTerm node : terms) { + if (node.weight == 0) continue; + constraint_vars.push_back(node.var); + constraint_coeffs.push_back(node.weight.value()); + } + constraint_vars.push_back(objective_var); + constraint_coeffs.push_back(-1); + model->Add(FixedWeightedSum(constraint_vars, constraint_coeffs, 0)); + } + + // Find out the true lower bound of new_var. This is called "cover + // optimization" in the max-SAT literature. + // + // TODO(user): Do more experiments to decide if this is better. This + // approach kind of mix the basic linear-scan one with the core based + // approach. + if (/* DISALBES CODE */ (false)) { + IntegerValue best = new_var_ub; + + // Simple linear scan algorithm to find the optimal of new_var. + while (best > new_var_lb) { + const Literal a = integer_encoder->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(new_var, best - 1)); + result = SolveIntegerProblemWithLazyEncoding({a}, next_decision, model); + if (result != SatSolver::MODEL_SAT) break; + best = integer_trail->LowerBound(new_var); + if (!process_solution()) { + result = SatSolver::MODEL_UNSAT; + break; + } + } + if (result == SatSolver::ASSUMPTIONS_UNSAT) { + sat_solver->Backtrack(0); + sat_solver->SetAssumptionLevel(0); + if (!integer_trail->Enqueue( + IntegerLiteral::GreaterOrEqual(new_var, best), {}, {})) { + result = SatSolver::MODEL_UNSAT; + break; + } + } + } + } + + // Returns MODEL_SAT if we found the optimal. + return num_solutions > 0 && result == SatSolver::MODEL_UNSAT + ? SatSolver::MODEL_SAT + : result; +} + SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( bool log_info, const std::vector& literals, const std::vector& int64_coeffs, diff --git a/ortools/sat/optimization.h b/ortools/sat/optimization.h index eb6c9a740c..32a4fecbfc 100644 --- a/ortools/sat/optimization.h +++ b/ortools/sat/optimization.h @@ -131,12 +131,27 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( const std::function& feasible_solution_observer, Model* model); +// Same as MinimizeIntegerVariableWithLinearScanAndLazyEncoding() but use +// a core-based approach instead. The given objective_var must be equal to the +// sum of the given variables using the given coefficients. +// +// TODO(user): It is not needed to have objective_var and the linear objective +// constraint encoded in the model. Remove this preconditions in order to +// improve the solving time. +SatSolver::Status MinimizeWithCoreAndLazyEncoding( + bool log_info, IntegerVariable objective_var, + const std::vector& variables, + const std::vector& coefficients, + const std::function& next_decision, + const std::function& feasible_solution_observer, + Model* model); + // Similar to MinimizeIntegerVariableWithLinearScanAndLazyEncoding() but use // a core based approach. Note that this require the objective to be given as // a weighted sum of literals // -// TODO(user): It should be easy to adapt this to a weighted sum of -// IntegerVariable. +// TODO(user): The function above is more general, remove this one after +// checking that the performances are similar. SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( bool log_info, const std::vector& literals, const std::vector& coeffs, diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 50483ff3f3..be453f5273 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -18,7 +18,7 @@ package operations_research.sat; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 83 +// NEXT TAG: 84 message SatParameters { // ========================================================================== // Branching and polarity @@ -479,4 +479,10 @@ message SatParameters { // really enforced properly, we don't support assumptions, and all the // decisions variables must be integers. optional bool use_fixed_search = 82 [default = false]; + + // The default optimization method is a simple "linear scan", each time trying + // to find a better solution than the previous one. If this is true, then we + // use a core-based approach (like in max-SAT) when we try to increase the + // lower bound instead. + optional bool optimize_with_core = 83 [default = false]; } diff --git a/ortools/sat/sat_solver.cc b/ortools/sat/sat_solver.cc index e12964fcc8..fb251ec310 100644 --- a/ortools/sat/sat_solver.cc +++ b/ortools/sat/sat_solver.cc @@ -349,7 +349,7 @@ void SatSolver::AddLearnedClauseAndEnqueueUnitPropagation( SCOPED_TIME_STAT(&stats_); // Note that we need to output the learned clause before cleaning the clause - // database. This is because we already backtacked and some of the clauses + // database. This is because we already backtracked and some of the clauses // that where needed to infer the conflict may not be "reasons" anymore and // may be deleted. if (drat_writer_ != nullptr) { @@ -544,7 +544,7 @@ bool SatSolver::PropagateAndStopAfterOneConflictResolution() { } } - // A conflict occured, compute a nice reason for this failure. + // A conflict occurred, compute a nice reason for this failure. same_reason_identifier_.Clear(); const int max_trail_index = ComputeMaxTrailIndex(trail_->FailingClause()); ComputeFirstUIPConflict(max_trail_index, &learned_conflict_, @@ -798,7 +798,7 @@ SatSolver::Status SatSolver::ReapplyDecisionsUpTo( *first_propagation_index = std::min(*first_propagation_index, index); if (index == kUnsatTrailIndex) return MODEL_UNSAT; if (current_decision_level_ <= old_level) { - // A conflict occured which backjumped to an earlier decision level. + // A conflict occurred which backjumped to an earlier decision level. // We potentially backjumped over some valid decisions, so we need to // continue the loop and try to re-enqueue them. // @@ -1045,7 +1045,7 @@ SatSolver::Status SatSolver::SolveInternal(TimeLimit* time_limit) { } if (!PropagateAndStopAfterOneConflictResolution()) { - // A conflict occured, continue the loop. + // A conflict occurred, continue the loop. if (is_model_unsat_) return StatusWithLog(MODEL_UNSAT); } else { // We need to reapply any assumptions that are not currently applied. diff --git a/ortools/sat/theta_tree.cc b/ortools/sat/theta_tree.cc index e9f8a56301..f15ffda867 100644 --- a/ortools/sat/theta_tree.cc +++ b/ortools/sat/theta_tree.cc @@ -19,188 +19,162 @@ namespace sat { ThetaLambdaTree::ThetaLambdaTree() {} // Make a tree using the first num_events events of the vectors. -void ThetaLambdaTree::Reset(std::vector event_initial_envelope, - std::vector event_energy) { - const int num_events = event_initial_envelope.size(); - DCHECK_EQ(num_events, event_energy.size()); - DCHECK(std::is_sorted(event_initial_envelope.begin(), - event_initial_envelope.end())); - event_initial_envelope_ = std::move(event_initial_envelope); - event_energy_ = std::move(event_energy); - +void ThetaLambdaTree::Reset(int num_events) { // Use 2^k leaves in the tree, with 2^k >= std::max(num_events, 2). num_events_ = num_events; for (num_leaves_ = 2; num_leaves_ < num_events; num_leaves_ <<= 1) { } const int num_nodes = 2 * num_leaves_; - tree_energy_.assign(num_nodes, IntegerValue(0)); + tree_energy_min_.assign(num_nodes, IntegerValue(0)); tree_energy_opt_.assign(num_nodes, IntegerValue(0)); tree_envelope_.assign(num_nodes, kMinIntegerValue); tree_envelope_opt_.assign(num_nodes, kMinIntegerValue); - is_present_.assign(num_events, false); - is_optional_.assign(num_events, false); } -void ThetaLambdaTree::AddEvent(int event) { +void ThetaLambdaTree::AddOrUpdateEvent(int event, IntegerValue initial_envelope, + IntegerValue energy_min, + IntegerValue energy_max) { DCHECK_LE(0, event); DCHECK_LT(event, num_events_); - DCHECK(!IsPresent(event) && !IsOptional(event)); - is_present_[event] = true; - is_optional_[event] = false; + DCHECK_LE(0, energy_min); + DCHECK_LE(energy_min, energy_max); const int node = num_leaves_ + event; - tree_envelope_[node] = event_initial_envelope_[event] + event_energy_[event]; - tree_energy_[node] = event_energy_[event]; - tree_envelope_opt_[node] = tree_envelope_[node]; - tree_energy_opt_[node] = tree_energy_[node]; - RefreshNode(node); -} - -void ThetaLambdaTree::AddOptionalEvent(int event) { - DCHECK_LE(0, event); - DCHECK_LT(event, num_events_); - DCHECK(!IsPresent(event) && !IsOptional(event)); - is_present_[event] = false; - is_optional_[event] = true; - const int node = num_leaves_ + event; - tree_envelope_[node] = kMinIntegerValue; - tree_energy_[node] = IntegerValue(0); - tree_envelope_opt_[node] = - event_initial_envelope_[event] + event_energy_[event]; - tree_energy_opt_[node] = event_energy_[event]; + tree_envelope_[node] = initial_envelope + energy_min; + tree_energy_min_[node] = energy_min; + tree_envelope_opt_[node] = initial_envelope + energy_max; + tree_energy_opt_[node] = energy_max; RefreshNode(node); } void ThetaLambdaTree::RemoveEvent(int event) { DCHECK_LE(0, event); DCHECK_LT(event, num_events_); - DCHECK(IsPresent(event) || IsOptional(event)); - is_present_[event] = false; - is_optional_[event] = false; const int node = num_leaves_ + event; tree_envelope_[node] = kMinIntegerValue; - tree_energy_[node] = IntegerValue(0); + tree_energy_min_[node] = IntegerValue(0); tree_envelope_opt_[node] = kMinIntegerValue; tree_energy_opt_[node] = IntegerValue(0); RefreshNode(node); } -bool ThetaLambdaTree::IsPresent(int event) const { - DCHECK_LE(0, event); - DCHECK_LT(event, num_events_); - return is_present_[event]; -} -bool ThetaLambdaTree::IsOptional(int event) const { - DCHECK_LE(0, event); - DCHECK_LT(event, num_events_); - return is_optional_[event]; -} - IntegerValue ThetaLambdaTree::GetEnvelope() const { return tree_envelope_[1]; } IntegerValue ThetaLambdaTree::GetOptionalEnvelope() const { return tree_envelope_opt_[1]; } -int ThetaLambdaTree::GetRightmostEventWithEnvelopeGreaterOrEqualTo( +int ThetaLambdaTree::GetMaxEventWithEnvelopeGreaterThan( IntegerValue target_envelope) const { - DCHECK_LE(target_envelope, tree_envelope_[1]); - const int event = GetLeafRealizingEnvelope(1, target_envelope) - num_leaves_; - DCHECK(IsPresent(event)); - return event; + DCHECK_LT(target_envelope, tree_envelope_[1]); + return GetMaxLeafWithEnvelopeGreaterThan(1, target_envelope) - num_leaves_; +} + +void ThetaLambdaTree::GetEventsWithOptionalEnvelopeGreaterThan( + IntegerValue target_envelope, int* critical_event, int* optional_event, + IntegerValue* available_energy) const { + int critical_leaf; + int optional_leaf; + GetLeavesWithOptionalEnvelopeGreaterThan(target_envelope, &critical_leaf, + &optional_leaf, available_energy); + *critical_event = critical_leaf - num_leaves_; + *optional_event = optional_leaf - num_leaves_; } IntegerValue ThetaLambdaTree::GetEnvelopeOf(int event) const { - DCHECK(IsPresent(event)); DCHECK_LE(0, event); DCHECK_LT(event, num_events_); IntegerValue env = tree_envelope_[event + num_leaves_]; for (int node = event + num_leaves_; node > 1; node >>= 1) { const int right = node | 1; - if (right != node) env += tree_energy_[right]; + if (right != node) env += tree_energy_min_[right]; } return env; } -int ThetaLambdaTree::GetOptionalEventRealizingOptionalEnvelope() const { - const int event = GetLeafRealizingOptionalEnvelope(1) - num_leaves_; - DCHECK(IsOptional(event)); - return event; -} - -int ThetaLambdaTree::GetEventRealizingOptionalEnvelope() const { - const int event = GetLeafRealizingOptionalEnvelope(1) - num_leaves_; - DCHECK(IsPresent(event) || IsOptional(event)); - return event; -} - void ThetaLambdaTree::RefreshNode(int node) { do { const int right = node | 1; const int left = right ^ 1; node >>= 1; - tree_energy_[node] = tree_energy_[left] + tree_energy_[right]; - tree_envelope_[node] = std::max(tree_envelope_[right], - tree_envelope_[left] + tree_energy_[right]); + tree_energy_min_[node] = tree_energy_min_[left] + tree_energy_min_[right]; + tree_envelope_[node] = std::max( + tree_envelope_[right], tree_envelope_[left] + tree_energy_min_[right]); tree_energy_opt_[node] = - std::max(tree_energy_[left] + tree_energy_opt_[right], - tree_energy_[right] + tree_energy_opt_[left]); + std::max(tree_energy_min_[left] + tree_energy_opt_[right], + tree_energy_min_[right] + tree_energy_opt_[left]); tree_envelope_opt_[node] = - std::max(tree_envelope_opt_[left] + tree_energy_[right], + std::max(tree_envelope_opt_[left] + tree_energy_min_[right], tree_envelope_[left] + tree_energy_opt_[right]); tree_envelope_opt_[node] = std::max(tree_envelope_opt_[node], tree_envelope_opt_[right]); } while (node > 1); } -int ThetaLambdaTree::GetLeafRealizingEnvelope(int node, - IntegerValue envelope) const { - DCHECK_LE(envelope, tree_envelope_[node]); +int ThetaLambdaTree::GetMaxLeafWithEnvelopeGreaterThan( + int node, IntegerValue target_envelope) const { + DCHECK_LT(target_envelope, tree_envelope_[node]); while (node < num_leaves_) { const int left = node << 1; const int right = left | 1; - if (envelope <= tree_envelope_[right]) { + if (target_envelope < tree_envelope_[right]) { node = right; } else { - envelope -= tree_energy_[right]; + target_envelope -= tree_energy_min_[right]; node = left; } } return node; } -int ThetaLambdaTree::GetOptionalLeafOfMaximalEnergy(int node) const { +int ThetaLambdaTree::GetMaxLeafWithOptionalEnergyGreaterThan( + int node, IntegerValue node_available_energy, + IntegerValue* available_energy) const { + DCHECK_LT(node_available_energy, tree_energy_opt_[node]); while (node < num_leaves_) { const int left = node << 1; const int right = left | 1; - if (tree_energy_opt_[node] == - tree_energy_[left] + tree_energy_opt_[right]) { + const IntegerValue available_energy_right = + node_available_energy - tree_energy_min_[left]; + if (available_energy_right < tree_energy_opt_[right]) { + node_available_energy = available_energy_right; node = right; - } else { // tree_energy_opt_[left] + tree_energy_[right] + } else { // available_energy_left < tree_energy_opt_[left] + node_available_energy -= tree_energy_min_[right]; node = left; } } + *available_energy = node_available_energy; return node; } -template -int ThetaLambdaTree::GetLeafRealizingOptionalEnvelope(int node) const { +void ThetaLambdaTree::GetLeavesWithOptionalEnvelopeGreaterThan( + IntegerValue target_envelope, int* critical_leaf, int* optional_leaf, + IntegerValue* available_energy) const { + DCHECK_LT(target_envelope, tree_envelope_opt_[1]); + int node = 1; while (node < num_leaves_) { const int left = node << 1; const int right = left | 1; - if (tree_envelope_opt_[node] == tree_envelope_opt_[right]) { + if (target_envelope < tree_envelope_opt_[right]) { node = right; - } else if (tree_envelope_opt_[node] == + } else if (target_envelope < tree_envelope_[left] + tree_energy_opt_[right]) { - return initial_envelope_leaf_wanted - ? GetLeafRealizingEnvelope(left, tree_envelope_[left]) - : GetOptionalLeafOfMaximalEnergy(right); - } else { // tree_envelope_opt_[left] + tree_energy_[right] + *critical_leaf = + GetMaxLeafWithEnvelopeGreaterThan(left, tree_envelope_[left] - 1); + *optional_leaf = GetMaxLeafWithOptionalEnergyGreaterThan( + right, target_envelope - tree_envelope_[left], available_energy); + return; + } else { // < tree_envelope_opt_[left] + tree_energy_min_[right] + target_envelope -= tree_energy_min_[right]; node = left; } } - return node; + *critical_leaf = node; + *optional_leaf = node; + *available_energy = + target_envelope - tree_envelope_[node] + tree_energy_min_[node]; } } // namespace sat diff --git a/ortools/sat/theta_tree.h b/ortools/sat/theta_tree.h index 8bf2219c1b..e04370e4b7 100644 --- a/ortools/sat/theta_tree.h +++ b/ortools/sat/theta_tree.h @@ -30,9 +30,7 @@ namespace sat { // right child, etc. To represent num_events events, we use the smallest // possible amount of leaves, num_leaves = 2^k >= num_events. The unused leaves // are filled with dummy values, as if they were absent events. -// The API gives access to events that realize the maximal output. When it does -// so, it always gives the rightmost ones if there are several solutions, in -// an attempt to minimize the number of events used in explanations. +// The API gives access to rightmost events that realize a given envelope. // // See: // _ (0) Petr Vilim's PhD thesis "Global Constraints in Scheduling". @@ -51,31 +49,33 @@ namespace sat { // algorithm, this generalization intends to provide a data structure that can // fit several algorithms. // This tree is based around the notion of events. It has events at its leaves -// that can be present, optional, or absent, internal nodes compute values on -// the set of events that are present/optional. -// Each event has an initial_envelope and an energy. -// The tree maintains: -// _ energy(node) = sum_{leaf \in leaves(node)} energy(leaf) +// that can be present or absent, and present events come with an +// initial_envelope, a minimal and a maximal energy. +// All nodes maintain values on the set of present events under them: +// _ energy_min(node) = sum_{leaf \in leaves(node)} energy_min(leaf) // _ envelope(node) = // max_{leaf \in leaves(node)} // initial_envelope(leaf) + -// sum_{leaf' \in leaves(node), leaf' at or right of leaf} energy(leaf'). +// sum_{leaf' \in leaves(node), leaf' >= leaf} energy_min(leaf'). // // Thus, the envelope of a leaf representing an event, when present, is -// initial_envelope(event) + energy(event). +// initial_envelope(event) + energy_min(event). // -// envelope_opt and energy_opt are similar, but compute the maximum value a node -// could have if one of the optional event were present: -// _ energy_opt(node) = sum_{leaf \in leaves(node)} energy(leaf) -// + max_{leaf \in leaves(node)} energy_opt(leaf) +// envelope_opt and energy_opt are similar, but represent the maximum value +// a node could have if one leaf took its maximum energy: +// _ energy_opt(node) = sum_{leaf \in leaves(node)} energy_min(leaf) +// + max_{leaf \in leaves(node)} +// energy_max(leaf) - energy_min(leaf) // _ envelope_opt(node) = // max_{leaf \in leaves(node)} // initial_envelope(leaf) + -// sum_{leaf' \in leaves(node), leaf' at or right of leaf} energy(leaf') + -// max_{leaf_opt \in leaves(node)} energy_opt(leaf_opt) +// sum_{leaf' \in leaves(node), leaf' >= leaf} energy_min(leaf') + +// max_{leaf_opt \in leaves(node)} +// . (energy_max(leaf_opt) - energy_min(leaf_opt)) // max_{leaf_opt \in leaves(node)} -// initial_envelope(leaf_opt) + energy_opt(leaf_opt) + -// sum_{leaf \in leaves(node), leaf at or right of leaf_opt} energy(leaf) +// initial_envelope(leaf_opt) + +// energy_max(leaf_opt) - energy_min(leaf_opt) + +// sum_{leaf' \in leaves(node), leaf' >= leaf_opt} energy_min(leaf') // // Most articles using theta-tree variants hack Vilim's original theta tree // for the disjunctive resource constraint by manipulating envelope and @@ -96,90 +96,87 @@ namespace sat { // tasks away to reason only on events. class ThetaLambdaTree { public: - // Build a reusable tree. Initialization is done with Reset(). + // Builds a reusable tree. Initialization is done with Reset(). ThetaLambdaTree(); - // Reset() initializes this class; other operations refer to the values - // passed by the last Reset(). - // Instead of allocating and de-allocating trees at every usage, i.e. at - // every Propagate() of the scheduling algorithms that uses it, - // this class allows to keep the same memory for each call. - // Reset() should be called with (vector) attributes for events, - // then the resulting Add/Remove/Get operations will use these attributes. - // Events should be sorted by nondecreasing initial_envelope. - void Reset(std::vector event_initial_envelope, - std::vector event_energy); + // Initializes this class for events in [0, num_events) and makes all of them + // absent. Instead of allocating and de-allocating trees at every usage, i.e. + // at every Propagate() of the scheduling algorithms that uses it, this class + // allows to keep the same memory for each call. + void Reset(int num_events); - // Make event present, compute the new envelopes in O(log n). - // The event must not be already present or optional. - void AddEvent(int event); + // Makes event present and updates its initial envelope and min/max energies. + // This updates the tree in O(log n). + void AddOrUpdateEvent(int event, IntegerValue initial_envelope, + IntegerValue energy_min, IntegerValue energy_max); - // Make event optional, compute the new envelopes in O(log n). - // The event must not be already present or optional. - void AddOptionalEvent(int event); - - // Make event absent, compute the new envelope in O(log n). - // The event must be present or optional. + // Makes event absent, compute the new envelope in O(log n). void RemoveEvent(int event); - // Returns true iff event is present, O(1). - bool IsPresent(int event) const; - - // Returns true iff event is optional, O(1). - bool IsOptional(int event) const; - - // Returns the current envelope of present tasks in O(1). + // Returns the maximum envelope using all the energy_min in O(1). IntegerValue GetEnvelope() const; - // Return the current envelope of present tasks + 1 optional task in O(1). + // Returns the maximum envelope using the energy min of all task but + // one and the energy max of the last one in O(1). IntegerValue GetOptionalEnvelope() const; - // Returns the rightmost event such that - // initial_envelope(event) + sum_{event' at or after event} energy(event'), - // is larger or equal to the given envelope, in O(log n). - // The given envelope can be at most GetEnvelope(). - int GetRightmostEventWithEnvelopeGreaterOrEqualTo( - IntegerValue target_envelope) const; + // Computes the maximum event s.t. GetEnvelopeOf(event) > envelope_max. + // There must be such an event, i.e. GetEnvelope() > envelope_max. + // This finds the maximum event e such that + // initial_envelope(e) + sum_{e' >= e} energy_min(e') > target_envelope. + // This operation is O(log n). + int GetMaxEventWithEnvelopeGreaterThan(IntegerValue target_envelope) const; - // Returns the envelope obtained by using the given event's initial_envelope - // plus the sum of energies of present events at or on the right of given - // event in O(log n). + // Returns initial_envelope(event) + sum_{event' >= event} energy_min(event'), + // in time O(log n). IntegerValue GetEnvelopeOf(int event) const; - // Returns the rightmost optional event responsible for the current optional - // envelope in O(log n). There must be an optional event, - // i.e. GetOptionalEnvelope() > GetEnvelope(). - int GetOptionalEventRealizingOptionalEnvelope() const; - - // Returns the rightmost event responsible for the initial_envelope part of - // the current optional envelope, in O(log n). The event can be present or - // optional. - int GetEventRealizingOptionalEnvelope() const; + // Computes a pair of events (critical_event, optional_event) such that + // if optional_event was at its maximum energy, the envelope of critical_event + // would be greater than target_envelope. + // This assumes that such a pair exists, i.e. GetOptionalEnvelope() + // should be greater than target_envelope. + // More formally, this finds events such that + // initial_envelope(critical_event) + + // sum_{event' >= critical_event} energy_min(event') + + // max_{optional_event >= critical_event} + // (energy_max(optional_event) - energy_min(optional_event)) + // > target envelope. + // For efficiency reasons, this also fills available_energy with the maximum + // value such that the optional envelope of the pair would be target_envelope, + // i.e. target_envelope - GetEnvelopeOf(event) + energy_min(optional_event). + // This operation is O(log n). + void GetEventsWithOptionalEnvelopeGreaterThan( + IntegerValue target_envelope, int* critical_event, int* optional_event, + IntegerValue* available_energy) const; private: // Propagates the change of leaf energies and envelopes towards the root. void RefreshNode(int leaf); - // Returns the rightmost present leaf under node responsible for given - // envelope, see GetRightmostEventWithEnvelopeGreaterOrEqualTo(). - int GetLeafRealizingEnvelope(int node, IntegerValue envelope) const; + // Finds the maximum leaf under node such that + // initial_envelope(leaf) + sum_{leaf' >= leaf} energy_min(leaf') + // > target_envelope. + int GetMaxLeafWithEnvelopeGreaterThan(int node, + IntegerValue target_envelope) const; - // Returns the rightmost optional leaf of maximal energy under node. - // There must be an optional leaf under node. - int GetOptionalLeafOfMaximalEnergy(int node) const; + // Returns the maximum leaf under node whose optional energy would overload + // node. + // Finds the maximum leaf under node such that + // sum_{leaf' under node} energy_min(leaf') + + // energy_max(leaf) - energy_min(leaf) > node_available_energy. + // available_energy will be the energy available for this leaf, + // i.e. node_available_energy - sum_{leaf' under node} energy_min(leaf') + + // energy_min(leaf). + int GetMaxLeafWithOptionalEnergyGreaterThan( + int node, IntegerValue node_available_energy, + IntegerValue* available_energy) const; - // If initial_envelope_leaf_wanted is true, returns the rightmost leaf - // such that the optional envelope is initial_envelope(leaf) plus some - // energies. It can be a present or an optional leaf. - // If initial_envelope_leaf_wanted is false, returns the rightmost optional - // leaf whose energy is in the optional envelope above. - // Notice that with a given set of present and optional events, - // the leaf returned by calling with initial_envelope_leaf_wanted set to true - // is the same or at the left as the leaf returned by calling with the boolean - // set to false. - // See class description for more details. - template - int GetLeafRealizingOptionalEnvelope(int node) const; + // Finds the leaves and energy relevant for + // GetEventsWithOptionalEnvelopeGreaterThan(). + void GetLeavesWithOptionalEnvelopeGreaterThan( + IntegerValue target_envelope, int* critical_leaf, int* optional_leaf, + IntegerValue* available_energy) const; // Number of events of the last Reset(); int num_events_; @@ -188,19 +185,11 @@ class ThetaLambdaTree { // that 2 <= num_leaves_ and num_events_ <= num_leaves_. int num_leaves_; - // Event characteristics: initial_envelope and energy. - std::vector event_initial_envelope_; - std::vector event_energy_; - // Envelopes and energies of nodes. std::vector tree_envelope_; - std::vector tree_energy_; + std::vector tree_energy_min_; std::vector tree_envelope_opt_; std::vector tree_energy_opt_; - - // Stores whether an event is currently present/optional. - std::vector is_present_; - std::vector is_optional_; }; } // namespace sat