diff --git a/src/sat/cp_constraints.cc b/src/sat/cp_constraints.cc index 7cc28e0cfa..911a647bc8 100644 --- a/src/sat/cp_constraints.cc +++ b/src/sat/cp_constraints.cc @@ -795,11 +795,17 @@ OneOfVarMinPropagator::OneOfVarMinPropagator( bool OneOfVarMinPropagator::Propagate() { // Compute the min of the lower-bound for the still possible variables. + // TODO(user): This could be optimized by keping more info from the last + // Propagate() calls. IntegerValue target_min = kMaxIntegerValue; + const IntegerValue current_min = integer_trail_->LowerBound(target_var_); for (int i = 0; i < vars_.size(); ++i) { if (trail_->Assignment().LiteralIsTrue(selectors_[i])) return true; if (trail_->Assignment().LiteralIsFalse(selectors_[i])) continue; target_min = std::min(target_min, integer_trail_->LowerBound(vars_[i])); + + // Abort if we can't get a better bound. + if (target_min <= current_min) return true; } if (target_min == kMaxIntegerValue) { // All false, conflit. @@ -807,8 +813,6 @@ bool OneOfVarMinPropagator::Propagate() { return false; } - // Increase the min of the target var if we got a better bound. - if (target_min <= integer_trail_->LowerBound(target_var_)) return true; literal_reason_.clear(); integer_reason_.clear(); for (int i = 0; i < vars_.size(); ++i) { diff --git a/src/sat/cumulative.cc b/src/sat/cumulative.cc index 7f937343bc..f8c83d8e38 100644 --- a/src/sat/cumulative.cc +++ b/src/sat/cumulative.cc @@ -11,9 +11,10 @@ // See the License for the specific language governing permissions and // limitations under the License. +#include "sat/cumulative.h" + #include -#include "sat/cumulative.h" #include "sat/disjunctive.h" #include "sat/overload_checker.h" #include "sat/sat_solver.h" @@ -29,64 +30,85 @@ std::function Cumulative( const IntegerVariable& capacity) { return [=](Model* model) { if (vars.empty()) return; + IntervalsRepository* intervals = model->GetOrCreate(); + IntegerEncoder* encoder = model->GetOrCreate(); - if (vars.size() == 1) { - if (intervals->IsOptional(vars[0])) { - model->Add(ConditionalLowerOrEqualWithOffset( - demands[0], capacity, 0, intervals->IsPresentLiteral(vars[0]))); - } else { - model->Add(LowerOrEqual(demands[0], capacity)); - } - return; - } + // Redundant constraints to ensure that the resource capacity is high enough + // for each task. Also ensure that no task consumes more resource than what + // is available. This is useful because the subsequent propagators do not + // filter the capacity variable very well. + for (int i = 0; i < demands.size(); ++i) { + if (intervals->MaxSize(vars[i]) == 0) continue; - // Detect a subset of intervals that needs to be in disjunction. - // - // 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 disjunction with [1, 3) for - // instance. We need to uniformize the handling of interval with size zero. - std::vector in_disjunction; - for (int i = 0; i < vars.size(); ++i) { - if (intervals->MinSize(vars[i]) > 0 && - 2 * model->Get(LowerBound(demands[i])) > - model->Get(UpperBound(capacity))) { - in_disjunction.push_back(vars[i]); - } - } - - // Add a disjunctive constraint on the intervals in in_disjunction. Do not - // create the cumulative at all when all intervals must be in disjunction. - // - // TODO(user): Do proper experiments to see how beneficial this is, the - // disjunctive will propagate more but is also using slower algorithms. That - // said, this is more a question of optimizing the disjunctive propagation - // code. - // - // TODO(user): Another "known" idea is to detect pair of tasks that must be - // in disjunction and to create a Boolean to indicate which one is before - // the other. It shouldn't change the propagation, but may result in a - // faster one with smaller explanations, and the solver can also take - // decision on such Boolean. - // - // TODO(user): A better place for stuff like this could be in the presolver - // so that it is easier to disable and play with alternatives. - if (in_disjunction.size() > 1) { - model->Add(Disjunctive(in_disjunction)); - } - if (in_disjunction.size() == vars.size()) { - // We still need to propagate the minimum capacity though. Note that - // this works since we don't have interval of size zero here. - for (int i = 0; i < vars.size(); ++i) { + if (intervals->MinSize(vars[i]) > 0) { if (intervals->IsOptional(vars[i])) { - model->Add(ConditionalLowerOrEqualWithOffset( - demands[i], capacity, 0, intervals->IsPresentLiteral(vars[i]))); + model->Add(ConditionalLowerOrEqual( + demands[i], capacity, intervals->IsPresentLiteral(vars[i]))); } else { model->Add(LowerOrEqual(demands[i], capacity)); } + continue; } - return; + + // At this point, we know that the duration variable is not fixed. + const Literal size_condition = + encoder->CreateAssociatedLiteral(IntegerLiteral::GreaterOrEqual( + intervals->SizeVar(vars[i]), IntegerValue(1))); + + if (intervals->IsOptional(vars[i])) { + const Literal condition = + Literal(model->Add(NewBooleanVariable()), true); + model->Add(ReifiedBoolAnd( + {size_condition, intervals->IsPresentLiteral(vars[i])}, condition)); + model->Add(ConditionalLowerOrEqual(demands[i], capacity, condition)); + } else { + model->Add( + ConditionalLowerOrEqual(demands[i], capacity, size_condition)); + } + } + + if (vars.size() == 1) return; + + const SatParameters& parameters = model->Get()->parameters(); + + // 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()) { + // 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 + // disjunction with [1, 3) for instance. We need to uniformize the + // handling of interval with size zero. + // + // TODO(user): improve the condition (see CL147454185). + std::vector in_disjunction; + for (int i = 0; i < vars.size(); ++i) { + if (intervals->MinSize(vars[i]) > 0 && + 2 * model->Get(LowerBound(demands[i])) > + model->Get(UpperBound(capacity))) { + in_disjunction.push_back(vars[i]); + } + } + + // Add a disjunctive constraint on the intervals in in_disjunction. Do not + // create the cumulative at all when all intervals must be in disjunction. + // + // TODO(user): Do proper experiments to see how beneficial this is, the + // disjunctive will propagate more but is also using slower algorithms. + // That said, this is more a question of optimizing the disjunctive + // propagation code. + // + // TODO(user): Another "known" idea is to detect pair of tasks that must + // be in disjunction and to create a Boolean to indicate which one is + // before the other. It shouldn't change the propagation, but may result + // in a faster one with smaller explanations, and the solver can also take + // decision on such Boolean. + // + // TODO(user): A better place for stuff like this could be in the + // presolver so that it is easier to disable and play with alternatives. + if (in_disjunction.size() > 1) model->Add(Disjunctive(in_disjunction)); + if (in_disjunction.size() == vars.size()) return; } Trail* trail = model->GetOrCreate(); @@ -100,7 +122,6 @@ std::function Cumulative( time_tabling->RegisterWith(model->GetOrCreate()); model->TakeOwnership(time_tabling); - const SatParameters& parameters = model->Get()->parameters(); // Propagator responsible for applying the Overload Checking filtering rule. // It increases the minimum of the capacity variable. @@ -132,6 +153,8 @@ std::function CumulativeTimeDecomposition( return [=](Model* model) { CHECK(model->Get(IsFixed(capacity_var))); + if (vars.empty()) return; + const int num_tasks = vars.size(); IntegerTrail* integer_trail = model->GetOrCreate(); SatSolver* sat_solver = model->GetOrCreate(); @@ -143,14 +166,9 @@ std::function CumulativeTimeDecomposition( std::vector demands; for (int t = 0; t < num_tasks; ++t) { - const IntervalVariable i = vars[t]; - - CHECK(!intervals->IsOptional(i)); - CHECK(intervals->SizeVar(i) == kNoIntegerVariable); + start_vars.push_back(intervals->StartVar(vars[t])); + end_vars.push_back(intervals->EndVar(vars[t])); CHECK(model->Get(IsFixed(demand_vars[t]))); - - start_vars.push_back(intervals->StartVar(i)); - end_vars.push_back(intervals->EndVar(i)); demands.push_back(integer_trail->LowerBound(demand_vars[t])); } @@ -164,37 +182,37 @@ std::function CumulativeTimeDecomposition( const IntegerValue capacity = integer_trail->UpperBound(capacity_var); - for (IntegerValue time = min_start; time <= max_end; ++time) { + for (IntegerValue time = min_start; time < max_end; ++time) { std::vector literals_with_coeff; for (int t = 0; t < num_tasks; ++t) { const IntegerValue start_min = integer_trail->LowerBound(start_vars[t]); const IntegerValue end_max = integer_trail->UpperBound(end_vars[t]); if (end_max <= time || time < start_min || demands[t] == 0) continue; - // True if the task overlap time. - const Literal overlap_true = - Literal(model->Add(NewBooleanVariable()), true); - // True if the task does not start after time. - const Literal start_true = encoder->CreateAssociatedLiteral( - IntegerLiteral::LowerOrEqual(start_vars[t], IntegerValue(time))); - // True if the task ends after time. - const Literal end_true = + // Task t consumes the resource at time if consume_condition is true. + std::vector consume_condition; + const Literal consume = Literal(model->Add(NewBooleanVariable()), true); + + // Task t consumes the resource at time if it is present. + if (intervals->IsOptional(vars[t])) { + consume_condition.push_back(intervals->IsPresentLiteral(vars[t])); + } + + // Task t overlaps time. + consume_condition.push_back(encoder->CreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(start_vars[t], IntegerValue(time)))); + consume_condition.push_back( encoder->CreateAssociatedLiteral(IntegerLiteral::GreaterOrEqual( - end_vars[t], IntegerValue(time + 1))); + end_vars[t], IntegerValue(time + 1)))); - bool sat = true; - // if overlap_true is true, then start_true is true. - sat &= sat_solver->AddBinaryClause(overlap_true.Negated(), start_true); - // if overlap is true, then end_true is true. - sat &= sat_solver->AddBinaryClause(overlap_true.Negated(), end_true); - // if start_true and end_true are true, then overlap_true is true. - sat &= sat_solver->AddTernaryClause(start_true.Negated(), - end_true.Negated(), overlap_true); + model->Add(ReifiedBoolAnd(consume_condition, consume)); - if (!sat) return; + // TODO(user): this is needed because we currently can't create a + // boolean variable if the model is unsat. + if (sat_solver->IsModelUnsat()) return; literals_with_coeff.push_back( - LiteralWithCoeff(overlap_true, Coefficient(demands[t].value()))); + LiteralWithCoeff(consume, Coefficient(demands[t].value()))); } // The profile cannot exceed the capacity at time. sat_solver->AddLinearConstraint(false, Coefficient(0), true, diff --git a/src/sat/cumulative.h b/src/sat/cumulative.h index 4a72db7cfe..feda168b0c 100644 --- a/src/sat/cumulative.h +++ b/src/sat/cumulative.h @@ -26,7 +26,7 @@ namespace sat { // and the capacity variables. // // Each interval represents a task to be scheduled in time such that the task -// consumes the resource during the time range [lb, ub[ where lb and ub +// consumes the resource during the time range [lb, ub) where lb and ub // respectively represent the lower and upper bounds of the corresponding // interval variable. The amount of resource consumed by the task is the value // of its associated demand variable. @@ -46,9 +46,8 @@ std::function Cumulative( // demands and the capacity variables. See the comment of Cumulative() above for // a definition of the constraint. // -// This constraint assumes that an interval cannot be optional and must have a -// fixed size (which can be zero). The demands and the capacity must be assigned -// to a non-negative number. +// This constraint assumes that task demands and the resource capacity must be +// assigned to a non-negative number. std::function CumulativeTimeDecomposition( const std::vector& vars, const std::vector& demand_vars, diff --git a/src/sat/disjunctive.cc b/src/sat/disjunctive.cc index 40500cd350..5f1761e76b 100644 --- a/src/sat/disjunctive.cc +++ b/src/sat/disjunctive.cc @@ -15,7 +15,6 @@ //#include "base/iterator_adaptors.h" #include "sat/sat_solver.h" -#include "util/sort.h" namespace operations_research { namespace sat { @@ -24,10 +23,10 @@ std::function Disjunctive( const std::vector& vars) { return [=](Model* model) { bool is_all_different = true; - IntervalsRepository* intervals = model->GetOrCreate(); + IntervalsRepository* repository = model->GetOrCreate(); for (const IntervalVariable var : vars) { - if (intervals->IsOptional(var) || intervals->MinSize(var) != 1 || - intervals->MaxSize(var) != 1) { + if (repository->IsOptional(var) || repository->MinSize(var) != 1 || + repository->MaxSize(var) != 1) { is_all_different = false; break; } @@ -41,20 +40,64 @@ std::function Disjunctive( return; } - DisjunctiveConstraint* disjunctive = new DisjunctiveConstraint( + SchedulingConstraintHelper* helper = new SchedulingConstraintHelper( vars, model->GetOrCreate(), model->GetOrCreate(), - intervals, model->GetOrCreate()); - disjunctive->RegisterWith(model->GetOrCreate()); - disjunctive->SetParameters(model->GetOrCreate()->parameters()); - model->TakeOwnership(disjunctive); + repository); + model->TakeOwnership(helper); + GenericLiteralWatcher* watcher = + model->GetOrCreate(); + + // We decided to create the propagators in this particular order, but it + // shouldn't matter much because of the different priorities used. + { + // Only one direction is needed by this one. + DisjunctiveOverloadChecker* overload_checker = + new DisjunctiveOverloadChecker(true, helper); + const int id = overload_checker->RegisterWith(watcher); + watcher->SetPropagatorPriority(id, 1); + model->TakeOwnership(overload_checker); + } + for (const bool time_direction : {true, false}) { + DisjunctiveDetectablePrecedences* detectable_precedences = + new DisjunctiveDetectablePrecedences(time_direction, helper); + const int id = detectable_precedences->RegisterWith(watcher); + watcher->SetPropagatorPriority(id, 2); + model->TakeOwnership(detectable_precedences); + } + for (const bool time_direction : {true, false}) { + DisjunctiveNotLast* not_last = + new DisjunctiveNotLast(time_direction, helper); + const int id = not_last->RegisterWith(watcher); + watcher->SetPropagatorPriority(id, 3); + model->TakeOwnership(not_last); + } + for (const bool time_direction : {true, false}) { + DisjunctiveEdgeFinding* edge_finding = + new DisjunctiveEdgeFinding(time_direction, helper); + const int id = edge_finding->RegisterWith(watcher); + watcher->SetPropagatorPriority(id, 4); + model->TakeOwnership(edge_finding); + } + if (model->GetOrCreate() + ->parameters() + .use_precedences_in_disjunctive_constraint()) { + for (const bool time_direction : {true, false}) { + DisjunctivePrecedences* precedences = new DisjunctivePrecedences( + time_direction, helper, model->GetOrCreate(), + model->GetOrCreate()); + const int id = precedences->RegisterWith(watcher); + watcher->SetPropagatorPriority(id, 5); + model->TakeOwnership(precedences); + } + } }; } -std::function DisjunctiveWithBooleanPrecedences( +std::function DisjunctiveWithBooleanPrecedencesOnly( const std::vector& vars) { return [=](Model* model) { SatSolver* sat_solver = model->GetOrCreate(); - IntervalsRepository* intervals = model->GetOrCreate(); + IntervalsRepository* repository = model->GetOrCreate(); PrecedencesPropagator* precedences = model->GetOrCreate(); for (int i = 0; i < vars.size(); ++i) { @@ -62,20 +105,22 @@ std::function DisjunctiveWithBooleanPrecedences( const BooleanVariable boolean_var = sat_solver->NewBooleanVariable(); const Literal i_before_j = Literal(boolean_var, true); const Literal j_before_i = i_before_j.Negated(); - precedences->AddConditionalPrecedence(intervals->EndVar(vars[i]), - intervals->StartVar(vars[j]), + precedences->AddConditionalPrecedence(repository->EndVar(vars[i]), + repository->StartVar(vars[j]), i_before_j); - precedences->AddConditionalPrecedence(intervals->EndVar(vars[j]), - intervals->StartVar(vars[i]), + precedences->AddConditionalPrecedence(repository->EndVar(vars[j]), + repository->StartVar(vars[i]), j_before_i); } } - DisjunctiveConstraint* disjunctive = new DisjunctiveConstraint( - vars, model->GetOrCreate(), model->GetOrCreate(), - intervals, precedences); - disjunctive->RegisterWith(model->GetOrCreate()); - disjunctive->SetParameters(sat_solver->parameters()); - model->TakeOwnership(disjunctive); + }; +} + +std::function DisjunctiveWithBooleanPrecedences( + const std::vector& vars) { + return [=](Model* model) { + model->Add(DisjunctiveWithBooleanPrecedencesOnly(vars)); + model->Add(Disjunctive(vars)); }; } @@ -90,7 +135,7 @@ void TaskSet::AddEntry(const Entry& e) { DCHECK(std::is_sorted(sorted_tasks_.begin(), sorted_tasks_.end())); // If the task is added after optimized_restart_, we know that we don't need - // to scan the task before optimized_restart_ in the next ComputeMinEnd(). + // to scan the task before optimized_restart_ in the next ComputeEndMin(). if (j <= optimized_restart_) optimized_restart_ = 0; } @@ -120,9 +165,9 @@ void TaskSet::RemoveEntryWithIndex(int index) { optimized_restart_ = 0; } -IntegerValue TaskSet::ComputeMinEnd(int task_to_ignore, +IntegerValue TaskSet::ComputeEndMin(int task_to_ignore, int* critical_index) const { - // The order in which we process tasks with the same min-start doesn't matter. + // The order in which we process tasks with the same start-min doesn't matter. DCHECK(std::is_sorted(sorted_tasks_.begin(), sorted_tasks_.end())); bool ignored = false; const int size = sorted_tasks_.size(); @@ -144,410 +189,178 @@ IntegerValue TaskSet::ComputeMinEnd(int task_to_ignore, return min_end; } -DisjunctiveConstraint::DisjunctiveConstraint( - const std::vector& non_overlapping_intervals, - Trail* trail, IntegerTrail* integer_trail, IntervalsRepository* intervals, - PrecedencesPropagator* precedences) - : non_overlapping_intervals_(non_overlapping_intervals), - trail_(trail), - integer_trail_(integer_trail), - intervals_(intervals), - precedences_(precedences), - stats_("DisjunctiveConstraint") {} - -void DisjunctiveConstraint::SwitchToMirrorProblem() { - std::swap(start_vars_, minus_end_vars_); - std::swap(end_vars_, minus_start_vars_); - std::swap(task_by_increasing_min_start_, task_by_decreasing_max_end_); - std::swap(task_by_increasing_min_end_, task_by_decreasing_max_start_); - std::swap(before_, after_); -} - -bool DisjunctiveConstraint::Propagate() { - SCOPED_TIME_STAT(&stats_); - - start_vars_.clear(); - end_vars_.clear(); - minus_end_vars_.clear(); - minus_start_vars_.clear(); - duration_vars_.clear(); - fixed_durations_.clear(); - reason_for_presence_.clear(); - task_is_currently_present_.clear(); - for (const IntervalVariable i : non_overlapping_intervals_) { - if (intervals_->IsOptional(i)) { - const Literal l = intervals_->IsPresentLiteral(i); - if (trail_->Assignment().LiteralIsFalse(l)) continue; - task_is_currently_present_.push_back( - trail_->Assignment().LiteralIsTrue(l)); - reason_for_presence_.push_back(l.Index()); - } else { - task_is_currently_present_.push_back(true); - reason_for_presence_.push_back(kNoLiteralIndex); - } - if (intervals_->SizeVar(i) == kNoIntegerVariable) { - duration_vars_.push_back(kNoIntegerVariable); - fixed_durations_.push_back(intervals_->MinSize(i)); - } else { - duration_vars_.push_back(intervals_->SizeVar(i)); - fixed_durations_.push_back(IntegerValue(0)); - } - start_vars_.push_back(intervals_->StartVar(i)); - end_vars_.push_back(intervals_->EndVar(i)); - minus_start_vars_.push_back(NegationOf(intervals_->StartVar(i))); - minus_end_vars_.push_back(NegationOf(intervals_->EndVar(i))); - } - - // Because this class doesn't add any new precedences, we can compute the set - // of task before (and after) some IntegerVariable only once before we enter - // the propagation loop below. - if (parameters_.use_precedences_in_disjunctive_constraint()) { - precedences_->ComputePrecedences(end_vars_, task_is_currently_present_, - &before_); - precedences_->ComputePrecedences(minus_start_vars_, - task_is_currently_present_, &after_); - } - - // This is variable as it depends on the optional tasks that are present. - const int num_tasks = start_vars_.size(); - DCHECK_EQ(num_tasks, minus_end_vars_.size()); - DCHECK_EQ(num_tasks, duration_vars_.size()); - DCHECK_EQ(num_tasks, fixed_durations_.size()); - DCHECK_EQ(num_tasks, reason_for_presence_.size()); - DCHECK_EQ(num_tasks, task_is_currently_present_.size()); - - // Initialize our "order" vectors with all the task indices. - if (task_by_increasing_min_start_.size() != num_tasks) { - task_by_increasing_min_start_.resize(num_tasks); - task_by_increasing_min_end_.resize(num_tasks); - task_by_decreasing_max_start_.resize(num_tasks); - task_by_decreasing_max_end_.resize(num_tasks); - for (int t = 0; t < num_tasks; ++t) { - task_by_increasing_min_start_[t] = t; - task_by_increasing_min_end_[t] = t; - task_by_decreasing_max_start_[t] = t; - task_by_decreasing_max_end_[t] = t; - } - } - - // Loop until we reach the fixed-point. It should be unique (see Petr Villim - // PhD). - // - // TODO(user): Some of these passes are idempotent, so there is no need to - // call them if their input didn't change! Improve that. - while (true) { - const int64 old_timestamp = integer_trail_->num_enqueues(); - - // This one is symmetric so there is no need to do it backward in time too. - // It only detects conflict, it never propagates anything. - if (!OverloadCheckingPass()) return false; - - if (!DetectablePrecedencePass()) return false; - SwitchToMirrorProblem(); - if (!DetectablePrecedencePass()) return false; - SwitchToMirrorProblem(); - if (old_timestamp != integer_trail_->num_enqueues()) continue; - - if (!NotLastPass()) return false; - SwitchToMirrorProblem(); - if (!NotLastPass()) return false; - SwitchToMirrorProblem(); - if (old_timestamp != integer_trail_->num_enqueues()) continue; - - if (!EdgeFindingPass()) return false; - SwitchToMirrorProblem(); - if (!EdgeFindingPass()) return false; - SwitchToMirrorProblem(); - if (old_timestamp != integer_trail_->num_enqueues()) continue; - - if (parameters_.use_precedences_in_disjunctive_constraint()) { - if (!PrecedencePass()) return false; - SwitchToMirrorProblem(); - if (!PrecedencePass()) return false; - SwitchToMirrorProblem(); - if (old_timestamp != integer_trail_->num_enqueues()) continue; - } - - break; - } - return true; -} - -void DisjunctiveConstraint::RegisterWith(GenericLiteralWatcher* watcher) { - const int id = watcher->Register(this); - for (const IntervalVariable i : non_overlapping_intervals_) { - watcher->WatchLowerBound(intervals_->StartVar(i), id); - watcher->WatchUpperBound(intervals_->EndVar(i), id); - if (intervals_->SizeVar(i) != kNoIntegerVariable) { - watcher->WatchLowerBound(intervals_->SizeVar(i), id); - } - if (intervals_->IsOptional(i)) { - watcher->WatchLiteral(intervals_->IsPresentLiteral(i), id); - } - } -} - -void DisjunctiveConstraint::AddPresenceAndDurationReason(int t) { - if (reason_for_presence_[t] != kNoLiteralIndex) { - literal_reason_.push_back(Literal(reason_for_presence_[t]).Negated()); - } - if (duration_vars_[t] != kNoIntegerVariable) { - integer_reason_.push_back( - integer_trail_->LowerBoundAsLiteral(duration_vars_[t])); - } -} - -void DisjunctiveConstraint::AddMinDurationReason(int t) { - if (duration_vars_[t] != kNoIntegerVariable) { - integer_reason_.push_back( - integer_trail_->LowerBoundAsLiteral(duration_vars_[t])); - } -} - -void DisjunctiveConstraint::AddMinStartReason(int t, IntegerValue lower_bound) { - integer_reason_.push_back( - IntegerLiteral::GreaterOrEqual(start_vars_[t], lower_bound)); -} - -void DisjunctiveConstraint::AddMinEndReason(int t, IntegerValue lower_bound) { - integer_reason_.push_back( - IntegerLiteral::GreaterOrEqual(end_vars_[t], lower_bound)); -} - -void DisjunctiveConstraint::AddMaxEndReason(int t, IntegerValue upper_bound) { - integer_reason_.push_back( - IntegerLiteral::LowerOrEqual(end_vars_[t], upper_bound)); -} - -void DisjunctiveConstraint::AddMaxStartReason(int t, IntegerValue upper_bound) { - integer_reason_.push_back( - IntegerLiteral::LowerOrEqual(start_vars_[t], upper_bound)); -} - -bool DisjunctiveConstraint::IncreaseMinStart(int t, - IntegerValue new_min_start) { - if (!integer_trail_->Enqueue( - IntegerLiteral::GreaterOrEqual(start_vars_[t], new_min_start), - literal_reason_, integer_reason_)) { - return false; - } - - // We propagate right away the new min-end lower-bound we have. - const IntegerValue min_end_lb = new_min_start + MinDuration(t); - if (MinEnd(t) < min_end_lb) { - integer_reason_.clear(); - literal_reason_.clear(); - AddMinStartReason(t, new_min_start); - AddMinDurationReason(t); - if (!integer_trail_->Enqueue( - IntegerLiteral::GreaterOrEqual(end_vars_[t], min_end_lb), - literal_reason_, integer_reason_)) { - return false; - } - } - - return true; -} - -bool DisjunctiveConstraint::DecreaseMaxEnd(int t, IntegerValue new_max_end) { - if (!integer_trail_->Enqueue( - IntegerLiteral::LowerOrEqual(end_vars_[t], new_max_end), - literal_reason_, integer_reason_)) { - return false; - } - - // We propagate right away the new max-start upper-bound we have. - const IntegerValue max_start_ub = new_max_end - MinDuration(t); - if (MaxStart(t) > max_start_ub) { - integer_reason_.clear(); - literal_reason_.clear(); - AddMaxEndReason(t, new_max_end); - AddMinDurationReason(t); - if (!integer_trail_->Enqueue( - IntegerLiteral::LowerOrEqual(start_vars_[t], max_start_ub), - literal_reason_, integer_reason_)) { - return false; - } - } - - return true; -} - -void DisjunctiveConstraint::UpdateTaskByIncreasingMinStart() { - IncrementalSort(task_by_increasing_min_start_.begin(), - task_by_increasing_min_start_.end(), [this](int t1, int t2) { - return MinStart(t1) < MinStart(t2); - }); -} - -void DisjunctiveConstraint::UpdateTaskByIncreasingMinEnd() { - IncrementalSort(task_by_increasing_min_end_.begin(), - task_by_increasing_min_end_.end(), - [this](int t1, int t2) { return MinEnd(t1) < MinEnd(t2); }); -} - -void DisjunctiveConstraint::UpdateTaskByDecreasingMaxStart() { - IncrementalSort(task_by_decreasing_max_start_.begin(), - task_by_decreasing_max_start_.end(), [this](int t1, int t2) { - return MaxStart(t1) > MaxStart(t2); - }); -} - -void DisjunctiveConstraint::UpdateTaskByDecreasingMaxEnd() { - IncrementalSort(task_by_decreasing_max_end_.begin(), - task_by_decreasing_max_end_.end(), - [this](int t1, int t2) { return MaxEnd(t1) > MaxEnd(t2); }); -} - -bool DisjunctiveConstraint::OverloadCheckingPass() { - SCOPED_TIME_STAT(&stats_); - UpdateTaskByDecreasingMaxEnd(); +bool DisjunctiveOverloadChecker::Propagate() { + helper_->SetTimeDirection(time_direction_); task_set_.Clear(); - for (auto it = task_by_decreasing_max_end_.rbegin(); - it != task_by_decreasing_max_end_.rend(); ++it) { - const int t = *it; - if (task_is_currently_present_[t]) { - task_set_.AddEntry({t, MinStart(t), MinDuration(t)}); + for (auto it = helper_->TaskByDecreasingEndMax().rbegin(); + it != helper_->TaskByDecreasingEndMax().rend(); ++it) { + const auto& task_time = *it; + const int t = task_time.task_index; + if (helper_->IsPresent(t)) { + task_set_.AddEntry({t, helper_->StartMin(t), helper_->DurationMin(t)}); } int critical_index = 0; const IntegerValue min_end_of_critical_tasks = - task_set_.ComputeMinEnd(/*task_to_ignore=*/-1, &critical_index); + task_set_.ComputeEndMin(/*task_to_ignore=*/-1, &critical_index); const std::vector& sorted_tasks = task_set_.SortedTasks(); // Check if we can detect that this optional task can't be executed. // TODO(user): This test doesn't detect all the cases. A better way to do // it would be like in the EdgeFindingPass() with the concept of gray tasks. - if (!task_is_currently_present_[t] && !sorted_tasks.empty()) { - // Skip if it was already shown not to be there. - if (trail_->Assignment().LiteralIsFalse( - Literal(reason_for_presence_[t]))) { - continue; - } - + if (!helper_->IsPresent(t) && !helper_->IsAbsent(t) && + !sorted_tasks.empty()) { const IntegerValue critical_start = sorted_tasks[critical_index].min_start; - if (MinEnd(t) <= critical_start) continue; - const IntegerValue new_min_end = - CapAdd(min_end_of_critical_tasks, - MinStart(t) >= critical_start - ? MinDuration(t) - : MinEnd(t) - critical_start /* cannot overflow */); - if (new_min_end > MaxEnd(t)) { + if (helper_->EndMin(t) <= critical_start) continue; + const IntegerValue new_min_end = CapAdd( + min_end_of_critical_tasks, + helper_->StartMin(t) >= critical_start + ? helper_->DurationMin(t) + : helper_->EndMin(t) - critical_start /* cannot overflow */); + if (new_min_end > helper_->EndMax(t)) { // TODO(user): This could be done lazily, like most of the loop to - // compute the reasons in this class. - integer_reason_.clear(); - literal_reason_.clear(); + // compute the reasons in this file. + helper_->ClearReason(); const IntegerValue window_start = - std::min(MinStart(t), sorted_tasks[critical_index].min_start); + std::min(helper_->StartMin(t), critical_start); const IntegerValue window_end = new_min_end - 1; for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; - AddPresenceAndDurationReason(ct); - AddMinStartReason(ct, window_start); - AddMaxEndReason(ct, window_end); + helper_->AddPresenceReason(ct); + helper_->AddDurationMinReason(ct); + helper_->AddStartMinReason(ct, window_start); + helper_->AddEndMaxReason(ct, window_end); } - AddMinDurationReason(t); - AddMinStartReason(t, window_start); - AddMaxEndReason(t, window_end); + helper_->AddDurationMinReason(t); + helper_->AddStartMinReason(t, window_start); + helper_->AddEndMaxReason(t, window_end); - DCHECK_NE(reason_for_presence_[t], kNoLiteralIndex); - integer_trail_->EnqueueLiteral( - Literal(reason_for_presence_[t]).Negated(), literal_reason_, - integer_reason_); + helper_->PushTaskAbsence(t); + continue; } } // Overload checking. - if (min_end_of_critical_tasks > MaxEnd(t)) { - DCHECK(task_is_currently_present_[t]); - literal_reason_.clear(); - integer_reason_.clear(); + if (min_end_of_critical_tasks > helper_->EndMax(t)) { + DCHECK(helper_->IsPresent(t)); + helper_->ClearReason(); const IntegerValue window_start = sorted_tasks[critical_index].min_start; const IntegerValue window_end = min_end_of_critical_tasks - 1; for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; - AddPresenceAndDurationReason(ct); - AddMinStartReason(ct, window_start); - AddMaxEndReason(ct, window_end); + helper_->AddPresenceReason(ct); + helper_->AddDurationMinReason(ct); + helper_->AddStartMinReason(ct, window_start); + helper_->AddEndMaxReason(ct, window_end); } - - return integer_trail_->ReportConflict(literal_reason_, integer_reason_); + return helper_->ReportConflict(); } } return true; } -bool DisjunctiveConstraint::DetectablePrecedencePass() { - SCOPED_TIME_STAT(&stats_); +int DisjunctiveOverloadChecker::RegisterWith(GenericLiteralWatcher* watcher) { + // This propagator reach the fix point in one pass. + const int id = watcher->Register(this); + helper_->WatchAllTasks(id, watcher); + return id; +} - UpdateTaskByIncreasingMinEnd(); - UpdateTaskByDecreasingMaxStart(); +bool DisjunctiveDetectablePrecedences::Propagate() { + helper_->SetTimeDirection(time_direction_); + const auto& task_by_increasing_min_end = helper_->TaskByIncreasingEndMin(); + const auto& task_by_decreasing_max_start = + helper_->TaskByDecreasingStartMax(); - const int num_tasks = start_vars_.size(); + const int num_tasks = helper_->NumTasks(); int queue_index = num_tasks - 1; task_set_.Clear(); - for (const int t : task_by_increasing_min_end_) { - const IntegerValue min_end = MinEnd(t); + for (const auto task_time : task_by_increasing_min_end) { + const int t = task_time.task_index; + const IntegerValue min_end = task_time.time; + + if (helper_->IsAbsent(t)) continue; + while (queue_index >= 0) { - const int to_insert = task_by_decreasing_max_start_[queue_index]; - if (min_end <= MaxStart(to_insert)) break; - if (task_is_currently_present_[to_insert]) { - task_set_.AddEntry( - {to_insert, MinStart(to_insert), MinDuration(to_insert)}); + const int to_insert = + task_by_decreasing_max_start[queue_index].task_index; + if (min_end <= helper_->StartMax(to_insert)) break; + if (helper_->IsPresent(to_insert)) { + task_set_.AddEntry({to_insert, helper_->StartMin(to_insert), + helper_->DurationMin(to_insert)}); } --queue_index; } // task_set_ contains all the tasks that must be executed before t. // They are in "dectable precedence" because their max_start is smaller than - // the min-end of t like so: + // the end-min of t like so: // [(the task t) - // (a task in task_set_)] - // From there, we deduce that the min-start of t is greater or equal to the - // min-end of the critical tasks. + // (a task in task_set_)] + // From there, we deduce that the start-min of t is greater or equal to the + // end-min of the critical tasks. // // Note that this works as well when task_is_currently_present_[t] is false. int critical_index = 0; const IntegerValue min_end_of_critical_tasks = - task_set_.ComputeMinEnd(/*task_to_ignore=*/t, &critical_index); - if (min_end_of_critical_tasks > MinStart(t)) { + task_set_.ComputeEndMin(/*task_to_ignore=*/t, &critical_index); + if (min_end_of_critical_tasks > helper_->StartMin(t)) { const std::vector& sorted_tasks = task_set_.SortedTasks(); - literal_reason_.clear(); - integer_reason_.clear(); + helper_->ClearReason(); // We need: - // - MaxStart(ct) < MinEnd(t) for the detectable precedence. - // - MinStart(ct) > window_start for the min_end_of_critical_tasks reason. + // - StartMax(ct) < EndMin(t) for the detectable precedence. + // - StartMin(ct) > window_start for the min_end_of_critical_tasks reason. const IntegerValue window_start = sorted_tasks[critical_index].min_start; for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; if (ct == t) continue; - DCHECK(task_is_currently_present_[ct]); - DCHECK_GE(MinStart(ct), window_start); - DCHECK_LT(MaxStart(ct), min_end); - AddPresenceAndDurationReason(ct); - AddMinStartReason(ct, window_start); - AddMaxStartReason(ct, min_end - 1); + DCHECK(helper_->IsPresent(ct)); + DCHECK_GE(helper_->StartMin(ct), window_start); + DCHECK_LT(helper_->StartMax(ct), min_end); + helper_->AddPresenceReason(ct); + helper_->AddDurationMinReason(ct); + helper_->AddStartMinReason(ct, window_start); + helper_->AddStartMaxReason(ct, min_end - 1); } - // Add the reason for t (we only need the min-end). - AddMinEndReason(t, MinEnd(t)); + // Add the reason for t (we only need the end-min). + helper_->AddEndMinReason(t, min_end); - // This augment the min-start of t and subsequently it can augment the + // This augment the start-min of t and subsequently it can augment the // next min_end_of_critical_tasks, but our deduction is still valid. - if (!IncreaseMinStart(t, min_end_of_critical_tasks)) return false; + if (!helper_->IncreaseStartMin(t, min_end_of_critical_tasks)) { + return false; + } // We need to reorder t inside task_set_. - task_set_.NotifyEntryIsNowLastIfPresent({t, MinStart(t), MinDuration(t)}); + task_set_.NotifyEntryIsNowLastIfPresent( + {t, helper_->StartMin(t), helper_->DurationMin(t)}); } } return true; } -bool DisjunctiveConstraint::PrecedencePass() { - SCOPED_TIME_STAT(&stats_); - const int num_tasks = start_vars_.size(); +int DisjunctiveDetectablePrecedences::RegisterWith( + GenericLiteralWatcher* watcher) { + const int id = watcher->Register(this); + helper_->WatchAllTasks(id, watcher); + watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); + return id; +} + +bool DisjunctivePrecedences::Propagate() { + helper_->SetTimeDirection(time_direction_); + + std::vector task_is_currently_present(helper_->EndVars().size()); + for (int i = 0; i < helper_->NumTasks(); ++i) { + task_is_currently_present[i] = helper_->IsPresent(i); + } + precedences_->ComputePrecedences(helper_->EndVars(), + task_is_currently_present, &before_); + + const int num_tasks = helper_->NumTasks(); // We don't care about the initial content of this vector. reason_for_beeing_before_.resize(num_tasks, kNoLiteralIndex); @@ -560,35 +373,37 @@ bool DisjunctiveConstraint::PrecedencePass() { for (; i < size && before_[i].var == var; ++i) { const int task = before_[i].index; reason_for_beeing_before_[task] = before_[i].reason; - task_set_.AddUnsortedEntry({task, MinStart(task), MinDuration(task)}); + task_set_.AddUnsortedEntry( + {task, helper_->StartMin(task), helper_->DurationMin(task)}); } if (task_set_.SortedTasks().size() < 2) continue; + if (integer_trail_->IsCurrentlyIgnored(var)) continue; task_set_.Sort(); const IntegerValue min_end = - task_set_.ComputeMinEnd(/*task_to_ignore=*/-1, &critical_index); + task_set_.ComputeEndMin(/*task_to_ignore=*/-1, &critical_index); if (min_end > integer_trail_->LowerBound(var)) { const std::vector& sorted_tasks = task_set_.SortedTasks(); - literal_reason_.clear(); - integer_reason_.clear(); + helper_->ClearReason(); const IntegerValue window_start = sorted_tasks[critical_index].min_start; for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; - DCHECK(task_is_currently_present_[ct]); - DCHECK_GE(MinStart(ct), window_start); - AddPresenceAndDurationReason(ct); - AddMinStartReason(ct, window_start); + DCHECK(helper_->IsPresent(ct)); + DCHECK_GE(helper_->StartMin(ct), window_start); + helper_->AddPresenceReason(ct); + helper_->AddDurationMinReason(ct); + helper_->AddStartMinReason(ct, window_start); if (reason_for_beeing_before_[ct] != kNoLiteralIndex) { - literal_reason_.push_back( + helper_->MutableLiteralReason()->push_back( Literal(reason_for_beeing_before_[ct]).Negated()); } } - // TODO(user): If var is actually a min-start of an interval, we - // could push the min-end and check the interval consistency right away. - if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(var, min_end), - literal_reason_, integer_reason_)) { + // TODO(user): If var is actually a start-min of an interval, we + // could push the end-min and check the interval consistency right away. + if (!helper_->PushIntegerLiteral( + IntegerLiteral::GreaterOrEqual(var, min_end))) { return false; } } @@ -596,28 +411,40 @@ bool DisjunctiveConstraint::PrecedencePass() { return true; } -bool DisjunctiveConstraint::NotLastPass() { - SCOPED_TIME_STAT(&stats_); - UpdateTaskByDecreasingMaxStart(); - UpdateTaskByDecreasingMaxEnd(); +int DisjunctivePrecedences::RegisterWith(GenericLiteralWatcher* watcher) { + // This propagator reach the fixed point in one go. + const int id = watcher->Register(this); + helper_->WatchAllTasks(id, watcher); + return id; +} - const int num_tasks = start_vars_.size(); +bool DisjunctiveNotLast::Propagate() { + helper_->SetTimeDirection(time_direction_); + const auto& task_by_decreasing_max_start = + helper_->TaskByDecreasingStartMax(); + + const int num_tasks = helper_->NumTasks(); int queue_index = num_tasks - 1; task_set_.Clear(); - for (auto it = task_by_decreasing_max_end_.rbegin(); - it != task_by_decreasing_max_end_.rend(); ++it) { - const int t = *it; - // task_set_ contains all the tasks that must start before the max-end of t. - // These are the only candidates that have a chance to decrease the max-end + for (auto it = helper_->TaskByDecreasingEndMax().rbegin(); + it != helper_->TaskByDecreasingEndMax().rend(); ++it) { + const auto& task_time = *it; + const int t = task_time.task_index; + const IntegerValue max_end = task_time.time; + + if (helper_->IsAbsent(t)) continue; + + // task_set_ contains all the tasks that must start before the end-max of t. + // These are the only candidates that have a chance to decrease the end-max // of t. - const IntegerValue max_end = MaxEnd(t); while (queue_index >= 0) { - const int to_insert = task_by_decreasing_max_start_[queue_index]; - if (max_end <= MaxStart(to_insert)) break; - if (task_is_currently_present_[to_insert]) { - task_set_.AddEntry( - {to_insert, MinStart(to_insert), MinDuration(to_insert)}); + const int to_insert = + task_by_decreasing_max_start[queue_index].task_index; + if (max_end <= helper_->StartMax(to_insert)) break; + if (helper_->IsPresent(to_insert)) { + task_set_.AddEntry({to_insert, helper_->StartMin(to_insert), + helper_->DurationMin(to_insert)}); } --queue_index; } @@ -626,25 +453,25 @@ bool DisjunctiveConstraint::NotLastPass() { // (i.e. it cannot be last): // // [(critical tasks) - // | <- t max-start + // | <- t start-max // - // So we can deduce that the max-end of t is smaller than or equal to the - // largest max-start of the critical tasks. + // So we can deduce that the end-max of t is smaller than or equal to the + // largest start-max of the critical tasks. // // Note that this works as well when task_is_currently_present_[t] is false. int critical_index = 0; const IntegerValue min_end_of_critical_tasks = - task_set_.ComputeMinEnd(/*task_to_ignore=*/t, &critical_index); - if (min_end_of_critical_tasks <= MaxStart(t)) continue; + task_set_.ComputeEndMin(/*task_to_ignore=*/t, &critical_index); + if (min_end_of_critical_tasks <= helper_->StartMax(t)) continue; - // Find the largest max-start of the critical tasks (excluding t). The - // max-end for t need to be smaller than or equal to this. + // Find the largest start-max of the critical tasks (excluding t). The + // end-max for t need to be smaller than or equal to this. IntegerValue largest_ct_max_start = kMinIntegerValue; const std::vector& sorted_tasks = task_set_.SortedTasks(); for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; if (t == ct) continue; - const IntegerValue max_start = MaxStart(ct); + const IntegerValue max_start = helper_->StartMax(ct); if (max_start > largest_ct_max_start) { largest_ct_max_start = max_start; } @@ -655,50 +482,68 @@ bool DisjunctiveConstraint::NotLastPass() { DCHECK(largest_ct_max_start == kMinIntegerValue || max_end > largest_ct_max_start); if (max_end > largest_ct_max_start) { - literal_reason_.clear(); - integer_reason_.clear(); + helper_->ClearReason(); const IntegerValue window_start = sorted_tasks[critical_index].min_start; for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; if (ct == t) continue; - AddPresenceAndDurationReason(ct); - AddMinStartReason(ct, window_start); - AddMaxStartReason(ct, largest_ct_max_start); + helper_->AddPresenceReason(ct); + helper_->AddDurationMinReason(ct); + helper_->AddStartMinReason(ct, window_start); + helper_->AddStartMaxReason(ct, largest_ct_max_start); } - // Add the reason for t, we only need the max-start. - AddMaxStartReason(t, min_end_of_critical_tasks - 1); + // Add the reason for t, we only need the start-max. + helper_->AddStartMaxReason(t, min_end_of_critical_tasks - 1); - // Enqueue the new max-end for t. + // Enqueue the new end-max for t. // Note that changing it will not influence the rest of the loop. - if (!DecreaseMaxEnd(t, largest_ct_max_start)) return false; + if (!helper_->DecreaseEndMax(t, largest_ct_max_start)) return false; } } return true; } -bool DisjunctiveConstraint::EdgeFindingPass() { - SCOPED_TIME_STAT(&stats_); - const int num_tasks = start_vars_.size(); +int DisjunctiveNotLast::RegisterWith(GenericLiteralWatcher* watcher) { + const int id = watcher->Register(this); + helper_->WatchAllTasks(id, watcher); + watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); + return id; +} + +bool DisjunctiveEdgeFinding::Propagate() { + helper_->SetTimeDirection(time_direction_); + const int num_tasks = helper_->NumTasks(); // Some task in sorted_tasks_ will be considered "gray". - // When computing the min-end of the sorted task, we will compute it for: + // When computing the end-min of the sorted task, we will compute it for: // - All the non-gray task // - All the non-gray task + at most one gray task. is_gray_.resize(num_tasks, false); // The algorithm is slightly different than the others. We start with a full // sorted_tasks, and remove tasks one by one starting by the one with highest - // max-end. + // end-max. task_set_.Clear(); - UpdateTaskByIncreasingMinStart(); + const auto& task_by_increasing_min_start = + helper_->TaskByIncreasingStartMin(); int num_not_gray = 0; - for (const int t : task_by_increasing_min_start_) { - task_set_.AddOrderedLastEntry({t, MinStart(t), MinDuration(t)}); + for (const auto task_time : task_by_increasing_min_start) { + const int t = task_time.task_index; + const IntegerValue min_start = task_time.time; + + // Even though we completely ignore absent tasks, we still mark them as + // gray to simplify the rest of the code in this function. + if (helper_->IsAbsent(t)) { + is_gray_[t] = true; + continue; + } + + task_set_.AddOrderedLastEntry({t, min_start, helper_->DurationMin(t)}); // We already mark all the non-present task as gray. - if (task_is_currently_present_[t]) { + if (helper_->IsPresent(t)) { ++num_not_gray; is_gray_[t] = false; } else { @@ -706,13 +551,13 @@ bool DisjunctiveConstraint::EdgeFindingPass() { } } - UpdateTaskByDecreasingMaxEnd(); + const auto& task_by_decreasing_max_end = helper_->TaskByDecreasingEndMax(); int decreasing_max_end_index = 0; // At each iteration we either transform a non-gray task into a gray one or // remove a gray task, so this loop is linear in complexity. while (num_not_gray > 0) { - // This is really similar to task_set_.ComputeMinEnd() except that we + // This is really similar to task_set_.ComputeEndMin() except that we // deal properly with the "gray" tasks. int critical_index = -1; int gray_critical_index = -1; @@ -725,19 +570,20 @@ bool DisjunctiveConstraint::EdgeFindingPass() { int gray_task_index = -1; // Set decreasing_max_end_index to the next non-gray task. - while (is_gray_[task_by_decreasing_max_end_[decreasing_max_end_index]]) { + while (is_gray_[task_by_decreasing_max_end[decreasing_max_end_index] + .task_index]) { ++decreasing_max_end_index; CHECK_LT(decreasing_max_end_index, num_tasks); } - const IntegerValue non_gray_max_end = - MaxEnd(task_by_decreasing_max_end_[decreasing_max_end_index]); + const IntegerValue non_gray_max_end = helper_->EndMax( + task_by_decreasing_max_end[decreasing_max_end_index].task_index); const std::vector& sorted_tasks = task_set_.SortedTasks(); for (int i = 0; i < sorted_tasks.size(); ++i) { const TaskSet::Entry& e = sorted_tasks[i]; if (is_gray_[e.task]) { if (e.min_start >= min_end_of_critical_tasks) { - // Is this gray task increasing the min-end by itself? + // Is this gray task increasing the end-min by itself? const IntegerValue candidate = CapAdd(e.min_start, e.min_duration); if (candidate >= min_end_of_critical_tasks_with_gray) { gray_critical_index = gray_task_index = i; @@ -754,7 +600,7 @@ bool DisjunctiveConstraint::EdgeFindingPass() { } } } else { - DCHECK_LE(MaxEnd(e.task), non_gray_max_end); + DCHECK_LE(helper_->EndMax(e.task), non_gray_max_end); // Augment the gray block. We could augment it more if e.min_start >= // min_end_of_critical_tasks_with_gray, but we don't care much about @@ -776,8 +622,7 @@ bool DisjunctiveConstraint::EdgeFindingPass() { // Overload checking. if (min_end_of_critical_tasks > non_gray_max_end) { - literal_reason_.clear(); - integer_reason_.clear(); + helper_->ClearReason(); // We need the reasons for the critical tasks to fall in: const IntegerValue window_start = sorted_tasks[critical_index].min_start; @@ -785,11 +630,12 @@ bool DisjunctiveConstraint::EdgeFindingPass() { for (int i = critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; if (is_gray_[ct]) continue; - AddPresenceAndDurationReason(ct); - AddMinStartReason(ct, window_start); - AddMaxEndReason(ct, window_end); + helper_->AddPresenceReason(ct); + helper_->AddDurationMinReason(ct); + helper_->AddStartMinReason(ct, window_start); + helper_->AddEndMaxReason(ct, window_end); } - return integer_trail_->ReportConflict(literal_reason_, integer_reason_); + return helper_->ReportConflict(); } // Optimization, we can remove all the gray tasks that starts at or after @@ -805,7 +651,7 @@ bool DisjunctiveConstraint::EdgeFindingPass() { // If we have a situation like: // [(critical_task_with_gray_task) // ] - // ^ max-end without the gray task. + // ^ end-max without the gray task. // // Then the gray task must be after all the critical tasks (all the non-gray // tasks in sorted_tasks actually), otherwise there will be no way to @@ -823,9 +669,8 @@ bool DisjunctiveConstraint::EdgeFindingPass() { CHECK_LE(gray_critical_index, critical_index); // Since the gray task is after all the other, we have a new lower bound. - if (min_end_of_critical_tasks > MinStart(gray_task)) { - literal_reason_.clear(); - integer_reason_.clear(); + if (min_end_of_critical_tasks > helper_->StartMin(gray_task)) { + helper_->ClearReason(); const IntegerValue low_start = sorted_tasks[gray_critical_index].min_start; const IntegerValue high_start = sorted_tasks[critical_index].min_start; @@ -833,22 +678,24 @@ bool DisjunctiveConstraint::EdgeFindingPass() { for (int i = gray_critical_index; i < sorted_tasks.size(); ++i) { const int ct = sorted_tasks[i].task; if (is_gray_[ct]) continue; - AddPresenceAndDurationReason(ct); - AddMinStartReason(ct, i >= critical_index ? high_start : low_start); - AddMaxEndReason(ct, window_end); + helper_->AddPresenceReason(ct); + helper_->AddDurationMinReason(ct); + helper_->AddStartMinReason( + ct, i >= critical_index ? high_start : low_start); + helper_->AddEndMaxReason(ct, window_end); } - // Add the reason for the gray_task (we don't need the max-end or + // Add the reason for the gray_task (we don't need the end-max or // presence reason). - AddMinDurationReason(gray_task); - AddMinStartReason(gray_task, low_start); + helper_->AddDurationMinReason(gray_task); + helper_->AddStartMinReason(gray_task, low_start); - // Enqueue the new min-start for gray_task. + // Enqueue the new start-min for gray_task. // // TODO(user): propagate the precedence Boolean here too? I think it // will be more powerful. Even if eventually all these precedence will // become detectable (see Petr Villim PhD). - if (!IncreaseMinStart(gray_task, min_end_of_critical_tasks)) { + if (!helper_->IncreaseStartMin(gray_task, min_end_of_critical_tasks)) { return false; } } @@ -864,25 +711,34 @@ bool DisjunctiveConstraint::EdgeFindingPass() { const IntegerValue threshold = std::max( min_end_of_critical_tasks, min_end_of_critical_tasks_with_gray); do { - DCHECK( - !is_gray_[task_by_decreasing_max_end_[decreasing_max_end_index]]); - is_gray_[task_by_decreasing_max_end_[decreasing_max_end_index]] = true; + DCHECK(!is_gray_[task_by_decreasing_max_end[decreasing_max_end_index] + .task_index]); + is_gray_[task_by_decreasing_max_end[decreasing_max_end_index] + .task_index] = true; --num_not_gray; if (num_not_gray == 0) break; // Set decreasing_max_end_index to the next non-gray task. - while ( - is_gray_[task_by_decreasing_max_end_[decreasing_max_end_index]]) { + while (is_gray_[task_by_decreasing_max_end[decreasing_max_end_index] + .task_index]) { ++decreasing_max_end_index; CHECK_LT(decreasing_max_end_index, num_tasks); } - } while (MaxEnd(task_by_decreasing_max_end_[decreasing_max_end_index]) >= - threshold); + } while ( + helper_->EndMax(task_by_decreasing_max_end[decreasing_max_end_index] + .task_index) >= threshold); } } return true; } +int DisjunctiveEdgeFinding::RegisterWith(GenericLiteralWatcher* watcher) { + const int id = watcher->Register(this); + helper_->WatchAllTasks(id, watcher); + watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id); + return id; +} + } // namespace sat } // namespace operations_research diff --git a/src/sat/disjunctive.h b/src/sat/disjunctive.h index b584ed10e8..10a44ee4d6 100644 --- a/src/sat/disjunctive.h +++ b/src/sat/disjunctive.h @@ -25,17 +25,26 @@ namespace operations_research { namespace sat { // Enforces a disjunctive (or no overlap) constraint on the given interval -// variables. +// variables. The intervals are interpreted as [start, end) and the constraint +// enforces that no time point belongs to two intervals. +// +// TODO(user): This is not completely true for empty intervals (start == end). +// Make sure such intervals are ignored by the constraint. std::function Disjunctive( const std::vector& vars); -// Same as Disjunctive() but also creates a Boolean variable for all the -// possible precedences of the form (task i is before task j). +// Creates Boolean variables for all the possible precedences of the form (task +// i is before task j) and forces that, for each couple of task (i,j), either i +// is before j or j is before i. Do not create any other propagators. +std::function DisjunctiveWithBooleanPrecedencesOnly( + const std::vector& vars); + +// Same as Disjunctive() + DisjunctiveWithBooleanPrecedencesOnly(). std::function DisjunctiveWithBooleanPrecedences( const std::vector& vars); -// Helper class to compute the min-end of a set of tasks given their min-start -// and min-duration. In Petr Vilim's PhD "Global Constraints in Scheduling", +// Helper class to compute the end-min of a set of tasks given their start-min +// and duration-min. In Petr Vilim's PhD "Global Constraints in Scheduling", // this corresponds to his Theta-tree except that we use a O(n) implementation // for most of the function here, not a O(log(n)) one. class TaskSet { @@ -70,26 +79,26 @@ class TaskSet { void AddUnsortedEntry(const Entry& e) { sorted_tasks_.push_back(e); } void Sort() { std::sort(sorted_tasks_.begin(), sorted_tasks_.end()); } - // Returns the min-end for the task in the set. The time profile of the tasks + // Returns the end-min for the task in the set. The time profile of the tasks // packed to the left will always be a set of contiguous tasks separated by // empty space: // // [Bunch of tasks] ... [Bunch of tasks] ... [critical tasks]. // // We call "critical tasks" the last group. These tasks will be solely - // responsible for for the min-end of the whole set. The returned + // responsible for for the end-min of the whole set. The returned // critical_index will be the index of the first critical task in // SortedTasks(). // // A reason for the min end is: - // - The min-duration of all the critical tasks. - // - The fact that all critical tasks have a min-start greater or equal to the + // - The duration-min of all the critical tasks. + // - The fact that all critical tasks have a start-min greater or equal to the // first of them, that is SortedTasks()[critical_index].min_start. // // It is possible to behave like if one task was not in the set by setting // task_to_ignore to the id of this task. This returns 0 if the set is empty // in which case critical_index will be left unchanged. - IntegerValue ComputeMinEnd(int task_to_ignore, int* critical_index) const; + IntegerValue ComputeEndMin(int task_to_ignore, int* critical_index) const; const std::vector& SortedTasks() const { return sorted_tasks_; } private: @@ -98,160 +107,97 @@ class TaskSet { DISALLOW_COPY_AND_ASSIGN(TaskSet); }; -// This class contains the propagation algorithms with explanation for a -// disjunctive constraint. -class DisjunctiveConstraint : public PropagatorInterface { +// ============================================================================ +// Below are many of the known propagation techniques for the disjunctive, each +// implemented in only one time direction and in its own propagator class. The +// Disjunctive() model function above will instanciate the used ones (according +// to the solver parameters) in both time directions. +// +// See Petr Vilim PhD "Global Constraints in Scheduling" for a description of +// some of the algorithm. +// ============================================================================ + +class DisjunctiveOverloadChecker : public PropagatorInterface { public: - // Creates a disjunctive constraint (or no overlap constraint) between the - // given IntervalVariable. - DisjunctiveConstraint( - const std::vector& non_overlapping_intervals, - Trail* trail, IntegerTrail* integer_trail, - IntervalsRepository* task_repository, PrecedencesPropagator* precedences); - - ~DisjunctiveConstraint() final { - IF_STATS_ENABLED(LOG(INFO) << stats_.StatString()); - } - - // The algorithm is quadratic in the number of tasks. + DisjunctiveOverloadChecker(bool time_direction, + SchedulingConstraintHelper* helper) + : time_direction_(time_direction), helper_(helper) {} bool Propagate() final; - - // Registers this constraint with the GenericLiteralWatcher. - void RegisterWith(GenericLiteralWatcher* watcher); - - void SetParameters(const SatParameters& parameters) { - parameters_ = parameters; - } + int RegisterWith(GenericLiteralWatcher* watcher); private: - // Reverses the time and all the relevant quantities. This is needed because - // most passes will propagate different stuff in the reverse time. For - // example, the not-last propagation will become the not-first one in the - // reverse time direction. - void SwitchToMirrorProblem(); + const bool time_direction_; + SchedulingConstraintHelper* helper_; + TaskSet task_set_; +}; - // Helpers for the current bounds on a task time window. - // [(min-duration) ... (min-duration)] - // ^ ^ ^ ^ - // min-start min-end max-start max-end - // - // Note that for tasks with variable durations, we don't necessarily have - // min-duration between the the min-XXX and max-XXX value. - IntegerValue MinDuration(int t) const { - return duration_vars_[t] == kNoIntegerVariable - ? fixed_durations_[t] - : integer_trail_->LowerBound(duration_vars_[t]); - } - IntegerValue MinStart(int t) const { - return integer_trail_->LowerBound(start_vars_[t]); - } - IntegerValue MaxStart(int t) const { - return integer_trail_->UpperBound(start_vars_[t]); - } - IntegerValue MinEnd(int t) const { - return integer_trail_->LowerBound(end_vars_[t]); - } - IntegerValue MaxEnd(int t) const { - return integer_trail_->UpperBound(end_vars_[t]); - } +class DisjunctiveDetectablePrecedences : public PropagatorInterface { + public: + DisjunctiveDetectablePrecedences(bool time_direction, + SchedulingConstraintHelper* helper) + : time_direction_(time_direction), helper_(helper) {} + bool Propagate() final; + int RegisterWith(GenericLiteralWatcher* watcher); - // Helper functions to compute the reason of a propagation. - // Append to literal_reason_ and integer_reason_ the corresponding reason. - void AddPresenceAndDurationReason(int t); - void AddMinDurationReason(int t); - void AddMinStartReason(int t, IntegerValue lower_bound); - void AddMinEndReason(int t, IntegerValue lower_bound); - void AddMaxEndReason(int t, IntegerValue upper_bound); - void AddMaxStartReason(int t, IntegerValue upper_bound); + private: + const bool time_direction_; + SchedulingConstraintHelper* helper_; + TaskSet task_set_; +}; - // Enqueues new bounds of an interval. The reasons (literal_reason_ and - // integer_reason_) must already be filled. Note that we automatically push - // min-end and max-start accordingly, so we maintain the invariants: - // - min-end >= min-start + min-duration - // - max-start <= max-end + min-duration - bool IncreaseMinStart(int t, IntegerValue new_min_start); - bool DecreaseMaxEnd(int t, IntegerValue new_max_end); +class DisjunctiveNotLast : public PropagatorInterface { + public: + DisjunctiveNotLast(bool time_direction, SchedulingConstraintHelper* helper) + : time_direction_(time_direction), helper_(helper) {} + bool Propagate() final; + int RegisterWith(GenericLiteralWatcher* watcher); - // All these passes use the algorithms described in Petr Vilim PhD "Global - // Constraints in Scheduling". Except that we don't use the O(log(n)) balanced - // tree to compute the min end of a set of task but O(n) algorithm. - // - // This allows to integrate right away any changes we make to the bounds and - // to have a simpler code. Since we also need to produces the reason of each - // propagation, it is unclear if tree-based structure will be a lot better - // here (especially because on the kind of problem we can solve n is - // relatively small, like 30). - bool OverloadCheckingPass(); - bool DetectablePrecedencePass(); - bool NotLastPass(); - bool EdgeFindingPass(); + private: + const bool time_direction_; + SchedulingConstraintHelper* helper_; + TaskSet task_set_; +}; - // Exploits the precedences relations of the form "this set of disjoint - // IntervalVariables must be performed before a given IntegerVariable". The - // relations are computed with PrecedencesPropagator::ComputePrecedences(). - bool PrecedencePass(); +class DisjunctiveEdgeFinding : public PropagatorInterface { + public: + DisjunctiveEdgeFinding(bool time_direction, + SchedulingConstraintHelper* helper) + : time_direction_(time_direction), helper_(helper) {} + bool Propagate() final; + int RegisterWith(GenericLiteralWatcher* watcher); - // Sorts the corresponding vectors of tasks. See below. - void UpdateTaskByIncreasingMinStart(); - void UpdateTaskByIncreasingMinEnd(); - void UpdateTaskByDecreasingMaxStart(); - void UpdateTaskByDecreasingMaxEnd(); + private: + const bool time_direction_; + SchedulingConstraintHelper* helper_; + TaskSet task_set_; + std::vector is_gray_; +}; - // The IntegerVariable of each tasks that must be considered for this - // constraint (note that the index is not an IntevalVariable but a new one - // local to this constraint). For optional tasks, we don't list the one that - // are already known to be absent. - std::vector start_vars_; - std::vector end_vars_; - std::vector duration_vars_; - std::vector fixed_durations_; - std::vector reason_for_presence_; +// Exploits the precedences relations of the form "this set of disjoint +// IntervalVariables must be performed before a given IntegerVariable". The +// relations are computed with PrecedencesPropagator::ComputePrecedences(). +class DisjunctivePrecedences : public PropagatorInterface { + public: + DisjunctivePrecedences(bool time_direction, + SchedulingConstraintHelper* helper, + IntegerTrail* integer_trail, + PrecedencesPropagator* precedences) + : time_direction_(time_direction), + helper_(helper), + integer_trail_(integer_trail), + precedences_(precedences) {} + bool Propagate() final; + int RegisterWith(GenericLiteralWatcher* watcher); - // The negation of the start/end variable so that SwitchToMirrorProblem() - // can do its job in O(1) instead of calling NegationOf() on each entry. - std::vector minus_start_vars_; - std::vector minus_end_vars_; - - // This is used by PrecedencePass(). - std::vector reason_for_beeing_before_; - std::vector before_; - std::vector after_; - - // True for the optional tasks that are known to be present and false for the - // one we don't know yet. The one that are known to be absent are not listed. - // The index is the same as for the *_vars vector. - std::vector task_is_currently_present_; - - // Various task orders. The indices are index in the *_vars vectors. We keep - // them separate so that in most cases they are almost sorted and std::sort() - // should be faster. The increasing/decreasing order is choosen so that - // SwitchToMirrorProblem() works in O(1) and permute these vectors - // accordingly. - // - // TODO(user): Using sets or maintaining these vectors sorted at all time - // may be faster. - std::vector task_by_increasing_min_start_; - std::vector task_by_increasing_min_end_; - std::vector task_by_decreasing_max_start_; - std::vector task_by_decreasing_max_end_; - - // Reason vectors. Declared here to avoid costly initialization. - std::vector literal_reason_; - std::vector integer_reason_; - - const std::vector non_overlapping_intervals_; - Trail* trail_; + private: + const bool time_direction_; + SchedulingConstraintHelper* helper_; IntegerTrail* integer_trail_; - IntervalsRepository* intervals_; PrecedencesPropagator* precedences_; TaskSet task_set_; - std::vector is_gray_; - - SatParameters parameters_; - mutable StatsGroup stats_; - - DISALLOW_COPY_AND_ASSIGN(DisjunctiveConstraint); + std::vector reason_for_beeing_before_; + std::vector before_; }; } // namespace sat diff --git a/src/sat/integer.cc b/src/sat/integer.cc index 9284785a09..2b06897a27 100644 --- a/src/sat/integer.cc +++ b/src/sat/integer.cc @@ -259,7 +259,8 @@ bool IntegerTrail::Propagate(Trail* trail) { // scan the new ones. But then we need a mecanism to detect the new ones. CHECK_EQ(trail->CurrentDecisionLevel(), 0); for (const auto& entry : encoder_->GetFullyEncodedVariables()) { - IntegerVariable var = entry.first; + const IntegerVariable var = entry.first; + if (IsCurrentlyIgnored(var)) continue; // This variable was already added and will be processed below. // Note that this is important, otherwise we may call many times @@ -292,12 +293,16 @@ bool IntegerTrail::Propagate(Trail* trail) { // Bound encoder. for (const IntegerLiteral i_lit : encoder_->GetIntegerLiterals(literal)) { + if (IsCurrentlyIgnored(i_lit.Var())) continue; + // The reason is simply the associated literal. if (!Enqueue(i_lit, {literal.Negated()}, {})) return false; } // Value encoder. - for (IntegerVariable var : watched_min_.Values(literal.Index())) { + for (const IntegerVariable var : watched_min_.Values(literal.Index())) { + DCHECK(!IsOptional(var)) << "Not supported yet"; + // A watched min value just became false. Note that because current_min_ // is also updated by Enqueue(), it may be larger than the watched min. const int min = current_min_.FindOrDie(var); @@ -381,12 +386,15 @@ IntegerVariable IntegerTrail::AddIntegerVariable(IntegerValue lower_bound, CHECK_EQ(vars_.size(), integer_trail_.size()); const IntegerVariable i(vars_.size()); - is_empty_literals_.push_back(kNoLiteralIndex); + is_ignored_literals_.push_back(kNoLiteralIndex); vars_.push_back({lower_bound, static_cast(integer_trail_.size())}); integer_trail_.push_back({lower_bound, i.value()}); + // TODO(user): the is_ignored_literals_ Booleans are currently always the same + // for a variable and its negation. So it may be better not to store it twice + // so that we don't have to be careful when setting them. CHECK_EQ(NegationOf(i).value(), vars_.size()); - is_empty_literals_.push_back(kNoLiteralIndex); + is_ignored_literals_.push_back(kNoLiteralIndex); vars_.push_back({-upper_bound, static_cast(integer_trail_.size())}); integer_trail_.push_back({-upper_bound, NegationOf(i).value()}); @@ -567,8 +575,11 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit, const std::vector& literal_reason, const std::vector& integer_reason) { DCHECK(AllLiteralsAreFalse(literal_reason)); + CHECK(!IsCurrentlyIgnored(i_lit.Var())); // Nothing to do if the bound is not better than the current one. + // TODO(user): Change this to a CHECK? propagator shouldn't try to push such + // bound and waste time explaining it. if (i_lit.bound <= vars_[i_lit.var].current_bound) return true; ++num_enqueues_; @@ -614,6 +625,9 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit, // make the IntegerLiteral bound stronger. const int min_index = FindWithDefault(current_min_, var, -1); if (min_index >= 0) { + DCHECK(!IsOptional(var)) + << "Fully encoded optional variable are not yet supported."; + // Recover the current min, and propagate to false all the values that // are in [min, i_lit.value). All these literals have the same reason, so // we use the "same reason as" mecanism. @@ -650,19 +664,21 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit, } // Check if the integer variable has an empty domain. - bool propagate_is_empty = false; if (i_lit.bound > UpperBound(var)) { - if (!IsOptional(var) || - trail_->Assignment().LiteralIsFalse(Literal(is_empty_literals_[var]))) { + // We relax the upper bound as much as possible to still have a conflict. + const auto ub_reason = IntegerLiteral::LowerOrEqual(var, i_lit.bound - 1); + + if (!IsOptional(var) || trail_->Assignment().LiteralIsFalse( + Literal(is_ignored_literals_[var]))) { std::vector* conflict = trail_->MutableConflict(); *conflict = literal_reason; if (IsOptional(var)) { - conflict->push_back(Literal(is_empty_literals_[var])); + conflict->push_back(Literal(is_ignored_literals_[var])); } // This is the same as: // MergeReasonInto(integer_reason, conflict); - // MergeReasonInto({UpperBoundAsLiteral(var)}, conflict); + // MergeReasonInto({ub_reason)}, conflict); // but with just one call to MergeReasonIntoInternal() for speed. Note // that it may also produce a smaller reason overall. DCHECK(tmp_queue_.empty()); @@ -672,19 +688,24 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit, if (trail_index >= size) tmp_queue_.push_back(trail_index); } { - // Add the reason for UpperBoundAsLiteral(var). - const int trail_index = - vars_[NegationOf(var).value()].current_trail_index; + const int trail_index = FindLowestTrailIndexThatExplainBound(ub_reason); if (trail_index >= size) tmp_queue_.push_back(trail_index); } MergeReasonIntoInternal(conflict); return false; } else { - const Literal is_empty = Literal(is_empty_literals_[var]); - if (!trail_->Assignment().LiteralIsTrue(is_empty)) { - // We delay the propagation for some subtle reasons, see below. - propagate_is_empty = true; + // Note(user): We never make the bound of an optional literal cross. We + // used to have a bug where we propagated these bound and their associated + // literal, and we where reaching a conflict while propagating the + // associated literal instead of setting is_ignored below to false. + const Literal is_ignored = Literal(is_ignored_literals_[var]); + if (integer_decision_levels_.empty()) { + trail_->EnqueueWithUnitReason(is_ignored); + } else { + EnqueueLiteral(is_ignored, literal_reason, integer_reason); + bounds_reason_buffer_.push_back(ub_reason); } + return true; } } @@ -705,10 +726,6 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit, // Special case for level zero. if (integer_decision_levels_.empty()) { - if (propagate_is_empty) { - const Literal is_empty = Literal(is_empty_literals_[var]); - trail_->EnqueueWithUnitReason(is_empty); - } vars_[i_lit.var].current_bound = i_lit.bound; integer_trail_[i_lit.var].bound = i_lit.bound; return true; @@ -718,10 +735,10 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit, {/*bound=*/i_lit.bound, /*var=*/i_lit.var, /*prev_trail_index=*/vars_[i_lit.var].current_trail_index, - /*literals_reason_start_index=*/static_cast( - literals_reason_buffer_.size()), - /*dependencies_start_index=*/static_cast( - bounds_reason_buffer_.size())}); + /*literals_reason_start_index=*/ + static_cast(literals_reason_buffer_.size()), + /*dependencies_start_index=*/ + static_cast(bounds_reason_buffer_.size())}); vars_[i_lit.var].current_bound = i_lit.bound; vars_[i_lit.var].current_trail_index = integer_trail_.size() - 1; @@ -736,17 +753,6 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit, bounds_reason_buffer_.insert(bounds_reason_buffer_.end(), integer_reason.begin(), integer_reason.end()); } - - // Subtle: We do that after the IntegerLiteral as been enqueued for 2 reasons: - // 1/ Any associated literal will use the TrailEntry above as a reason. - // 2/ The precendence propagator currently rely on this. - // TODO(user): The point 2 is brittle, try to clean it up. - if (propagate_is_empty) { - const Literal is_empty = Literal(is_empty_literals_[var]); - EnqueueLiteral(is_empty, literal_reason, integer_reason); - bounds_reason_buffer_.push_back(UpperBoundAsLiteral(var)); - } - return true; } @@ -902,8 +908,6 @@ ClauseRef IntegerTrail::Reason(const Trail& trail, int trail_index) const { tmp_queue_.push_back(next_trail_index); } MergeReasonIntoInternal(reason); - - if (DEBUG_MODE && trail.CurrentDecisionLevel() > 0) CHECK(!reason->empty()); return ClauseRef(*reason); } @@ -917,6 +921,18 @@ void IntegerTrail::EnqueueLiteral( return; } + // This may not indicate an incorectness, but just some propagators that + // didn't reach a fixed-point at level zero. + if (DEBUG_MODE && !integer_decision_levels_.empty()) { + std::vector l = literal_reason; + MergeReasonInto(integer_reason, &l); + LOG_IF(WARNING, l.empty()) + << "Propagating a literal with no reason at a positive level!\n" + << "level:" << integer_decision_levels_.size() << " lit:" << literal + << " " << ReasonDebugString(literal_reason, integer_reason) << "\n" + << DebugString(); + } + const int trail_index = trail_->Index(); if (trail_index >= boolean_trail_index_to_integer_one_.size()) { boolean_trail_index_to_integer_one_.resize(trail_index + 1); @@ -924,10 +940,10 @@ void IntegerTrail::EnqueueLiteral( boolean_trail_index_to_integer_one_[trail_index] = integer_trail_.size(); integer_trail_.push_back({/*bound=*/IntegerValue(0), /*var=*/-1, /*prev_trail_index=*/-1, - /*literals_reason_start_index=*/static_cast( - literals_reason_buffer_.size()), - /*dependencies_start_index=*/static_cast( - bounds_reason_buffer_.size())}); + /*literals_reason_start_index=*/ + static_cast(literals_reason_buffer_.size()), + /*dependencies_start_index=*/ + static_cast(bounds_reason_buffer_.size())}); literals_reason_buffer_.insert(literals_reason_buffer_.end(), literal_reason.begin(), literal_reason.end()); bounds_reason_buffer_.insert(bounds_reason_buffer_.end(), @@ -1141,55 +1157,162 @@ void GenericLiteralWatcher::RegisterReversibleInt(int id, int* rev) { id_to_reversible_ints_[id].push_back(rev); } -SatSolver::Status SolveIntegerProblemWithLazyEncoding( - const std::vector& assumptions, - const std::vector& variables_to_finalize, Model* model) { - TimeLimit* time_limit = model->Mutable(); +// TODO(user): the complexity of this heuristic and the one below is ok when +// use_fixed_search() is false because it is not executed often, however, we do +// a linear scan for each search decision when use_fixed_search() is true, which +// seems bad. Improve. +std::function FirstUnassignedVarAtItsMinHeuristic( + const std::vector& vars, Model* model) { IntegerEncoder* const integer_encoder = model->GetOrCreate(); IntegerTrail* const integer_trail = model->GetOrCreate(); - SatSolver* const solver = model->GetOrCreate(); - SatSolver::Status status = SatSolver::Status::MODEL_UNSAT; - for (;;) { - 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. - status = solver->SolveWithTimeLimit(time_limit); - } else { - status = - solver->ResetAndSolveWithGivenAssumptions(assumptions, time_limit); + return [integer_encoder, integer_trail, /*copy*/ vars]() { + for (const IntegerVariable var : vars) { + // Note that there is no point trying to fix a currently ignored variable. + if (integer_trail->IsCurrentlyIgnored(var)) continue; + const IntegerValue lb = integer_trail->LowerBound(var); + if (lb < integer_trail->UpperBound(var)) { + return integer_encoder + ->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(var, lb)) + .Index(); + } } - if (status != SatSolver::MODEL_SAT) break; + return kNoLiteralIndex; + }; +} - // Look for an integer variable whose domain is not singleton. - // - // Heuristic: we take the first non-fixed variable in the given order and - // try to fix it to its lower bound. - // - // TODO(user): consider better heuristics. Maybe we can use some kind of - // IntegerVariable activity or something from the CP world. We also tried - // to take the one with minimum lb amongst the given variables, which work - // well on some problem but not on other. +std::function UnassignedVarWithLowestMinAtItsMinHeuristic( + const std::vector& vars, Model* model) { + IntegerEncoder* const integer_encoder = model->GetOrCreate(); + IntegerTrail* const integer_trail = model->GetOrCreate(); + return [integer_encoder, integer_trail, /*copy */ vars]() { IntegerVariable candidate = kNoIntegerVariable; - for (const IntegerVariable variable : variables_to_finalize) { - // Note that we use < and not != for optional variable whose domain may - // be empty. - const IntegerValue lb = integer_trail->LowerBound(variable); - if (lb < integer_trail->UpperBound(variable)) { - candidate = variable; - break; + IntegerValue candidate_lb; + for (const IntegerVariable var : vars) { + if (integer_trail->IsCurrentlyIgnored(var)) continue; + const IntegerValue lb = integer_trail->LowerBound(var); + if (lb < integer_trail->UpperBound(var) && + (candidate == kNoIntegerVariable || lb < candidate_lb)) { + candidate = var; + candidate_lb = lb; + } + } + if (candidate == kNoIntegerVariable) return kNoLiteralIndex; + return integer_encoder + ->GetOrCreateAssociatedLiteral( + IntegerLiteral::LowerOrEqual(candidate, candidate_lb)) + .Index(); + }; +} + +SatSolver::Status SolveIntegerProblemWithLazyEncoding( + const std::vector& assumptions, + const std::function& next_decision, Model* model) { + // TODO(user): Improve the time-limit situation, especially when + // use_fixed_search() is true where the deterministic time is not updated. + TimeLimit* time_limit = model->Mutable(); + if (time_limit == nullptr) { + model->SetSingleton(TimeLimit::Infinite()); + time_limit = model->Mutable(); + } + SatSolver* const solver = model->GetOrCreate(); + 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()) { + // If we are not in fixed search, let the SAT solver do a full search to + // 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 { + // TODO(user): We actually don't want to reset the solver from one + // loop to the next as it is less efficient. + const SatSolver::Status status = + solver->ResetAndSolveWithGivenAssumptions(assumptions, time_limit); + if (status != SatSolver::MODEL_SAT) return status; } } - if (candidate == kNoIntegerVariable) break; + // Look for the "best" non-instanciated integer variable. + // If all variables are instanciated, we have a solution. + const LiteralIndex decision = + next_decision == nullptr ? kNoLiteralIndex : next_decision(); + if (decision == kNoLiteralIndex) return SatSolver::MODEL_SAT; - // This add a literal with good polarity. It is very important that the - // decision heuristic assign it to false! Otherwise, our heuristic is not - // good. - integer_encoder->CreateAssociatedLiteral(IntegerLiteral::GreaterOrEqual( - candidate, integer_trail->LowerBound(candidate) + 1)); + // Always try to Enqueue this literal as the next decision so we bypass + // any default polarity heuristic the solver may use on this literal. + solver->EnqueueDecisionAndBackjumpOnConflict(Literal(decision)); + if (solver->IsModelUnsat()) return SatSolver::MODEL_UNSAT; } - return status; + return SatSolver::Status::LIMIT_REACHED; +} + +// Shortcut when there is no assumption, and we consider all variables in +// order for the search decision. +SatSolver::Status SolveIntegerProblemWithLazyEncoding(Model* model) { + const IntegerVariable num_vars = + model->GetOrCreate()->NumIntegerVariables(); + std::vector all_variables; + for (IntegerVariable var(0); var < num_vars; ++var) { + all_variables.push_back(var); + } + return SolveIntegerProblemWithLazyEncoding( + {}, FirstUnassignedVarAtItsMinHeuristic(all_variables, model), model); +} + +// This is really close to ExcludeCurrentSolutionAndBacktrack(). +std::function +ExcludeCurrentSolutionWithoutIgnoredVariableAndBacktrack() { + return [=](Model* model) { + SatSolver* sat_solver = model->GetOrCreate(); + IntegerTrail* integer_trail = model->GetOrCreate(); + IntegerEncoder* encoder = model->GetOrCreate(); + + const int current_level = sat_solver->CurrentDecisionLevel(); + std::vector clause_to_exclude_solution; + clause_to_exclude_solution.reserve(current_level); + for (int i = 0; i < current_level; ++i) { + bool include_decision = true; + const Literal decision = sat_solver->Decisions()[i].literal; + + // Tests if this decision is associated to a bound of an ignored variable + // in the current assignment. + const InlinedIntegerLiteralVector& associated_literals = + encoder->GetIntegerLiterals(decision); + for (const IntegerLiteral bound : associated_literals) { + if (integer_trail->IsCurrentlyIgnored(bound.Var())) { + // In this case we replace the decision (which is a bound on an + // ignored variable) with the fact that the integer variable was + // ignored. This works because the only impact a bound of an ignored + // variable can have on the rest of the model is through the + // is_ignored literal. + clause_to_exclude_solution.push_back( + integer_trail->IsIgnoredLiteral(bound.Var()).Negated()); + include_decision = false; + } + } + + if (include_decision) { + clause_to_exclude_solution.push_back(decision.Negated()); + } + } + + // Note that it is okay to add duplicates literals in ClauseConstraint(), + // the clause will be preprocessed correctly. + sat_solver->Backtrack(0); + model->Add(ClauseConstraint(clause_to_exclude_solution)); + }; } } // namespace sat diff --git a/src/sat/integer.h b/src/sat/integer.h index 4852983103..1688aa1426 100644 --- a/src/sat/integer.h +++ b/src/sat/integer.h @@ -18,6 +18,7 @@ #include #include "base/port.h" +#include "base/logging.h" #include "base/join.h" #include "base/int_type.h" #include "base/map_util.h" @@ -119,7 +120,12 @@ struct IntegerLiteral { return var != o.var || bound == o.bound; } - std::string DebugString() const { return StrCat("I", var, ">=", bound.value()); } + IntegerVariable Var() const { return IntegerVariable(var); } + IntegerValue Bound() const { return bound; } + + std::string DebugString() const { + return StrCat("I", var, ">=", bound.value()); + } private: friend class IntegerEncoder; @@ -403,16 +409,24 @@ class IntegerTrail : public SatPropagator { } // The domain [lb, ub] of an "optional" variable is allowed to be empty (i.e. - // ub < lb) iff the given is_empty literal is true. + // ub < lb) if the given is_ignored literal is true. Moreover, if is_ignored + // is true, then the bound of such variable should NOT impact any non-ignored + // variable in any way (but the reverse is not true). bool IsOptional(IntegerVariable i) const { - return is_empty_literals_[i] != kNoLiteralIndex; + return is_ignored_literals_[i] != kNoLiteralIndex; } - LiteralIndex GetIsEmptyLiteral(IntegerVariable i) const { - return is_empty_literals_[i]; + bool IsCurrentlyIgnored(IntegerVariable i) const { + const LiteralIndex is_ignored_literal = is_ignored_literals_[i]; + return is_ignored_literal != kNoLiteralIndex && + trail_->Assignment().LiteralIsTrue(Literal(is_ignored_literal)); } - void MarkIntegerVariableAsOptional(IntegerVariable i, Literal is_non_empty) { - is_empty_literals_[i] = is_non_empty.NegatedIndex(); - is_empty_literals_[NegationOf(i)] = is_non_empty.NegatedIndex(); + Literal IsIgnoredLiteral(IntegerVariable i) const { + DCHECK(IsOptional(i)); + return Literal(is_ignored_literals_[i]); + } + void MarkIntegerVariableAsOptional(IntegerVariable i, Literal is_considered) { + is_ignored_literals_[i] = is_considered.NegatedIndex(); + is_ignored_literals_[NegationOf(i)] = is_considered.NegatedIndex(); } // Returns the current lower/upper bound of the given integer variable. @@ -576,8 +590,8 @@ class IntegerTrail : public SatPropagator { std::vector literals_reason_buffer_; mutable std::vector bounds_reason_buffer_; - // The "is_empty" literal of the optional variables or kNoLiteralIndex. - ITIVector is_empty_literals_; + // The "is_ignored" literal of the optional variables or kNoLiteralIndex. + ITIVector is_ignored_literals_; // Data used to support the propagation of fully encoded variable. We keep // for each variable the index in encoder_.GetDomainEncoding() of the first @@ -653,7 +667,6 @@ class PropagatorInterface { // indices of the literals modified after the registration will be present. virtual bool IncrementalPropagate(const std::vector& watch_indices) { LOG(FATAL) << "Not implemented."; - return false; } }; @@ -882,6 +895,7 @@ inline std::function ConstantIntegerVariable( inline std::function NewIntegerVariable(int64 lb, int64 ub) { return [=](Model* model) { + CHECK_LE(lb, ub); return model->GetOrCreate()->AddIntegerVariable( IntegerValue(lb), IntegerValue(ub)); }; @@ -894,15 +908,44 @@ inline std::function NewIntegerVariable( }; } +// Constraints might not accept Literals as input, forcing the user to pass +// 0-1 integer views of a literal. +// This class contains all such literal views of a given model, so that asking +// for a view of a literal will always return the same IntegerVariable. +class LiteralViews { + public: + explicit LiteralViews(Model* model) : model_(model) {} + + static LiteralViews* CreateInModel(Model* model) { + return new LiteralViews(model); + } + + IntegerVariable GetIntegerView(const Literal lit) { + const LiteralIndex index = lit.Index(); + + if (!ContainsKey(literal_index_to_integer_, index)) { + const IntegerVariable int_var = model_->Add(NewIntegerVariable(0, 1)); + model_->GetOrCreate() + ->FullyEncodeVariableUsingGivenLiterals( + int_var, {lit.Negated(), lit}, + {IntegerValue(0), IntegerValue(1)}); + literal_index_to_integer_[index] = int_var; + } + + return literal_index_to_integer_[index]; + } + + private: + hash_map literal_index_to_integer_; + Model* model_; +}; + // Creates a 0-1 integer variable "view" of the given literal. It will have a // value of 1 when the literal is true, and 0 when the literal is false. inline std::function NewIntegerVariableFromLiteral( Literal lit) { return [=](Model* model) { - const IntegerVariable int_var = model->Add(NewIntegerVariable(0, 1)); - model->GetOrCreate()->FullyEncodeVariableUsingGivenLiterals( - int_var, {lit.Negated(), lit}, {IntegerValue(0), IntegerValue(1)}); - return int_var; + return model->GetOrCreate()->GetIntegerView(lit); }; } @@ -929,7 +972,7 @@ inline std::function IsFixed(IntegerVariable v) { inline std::function Value(IntegerVariable v) { return [=](const Model& model) { const IntegerTrail* trail = model.Get(); - CHECK_EQ(trail->LowerBound(v), trail->UpperBound(v)); + CHECK_EQ(trail->LowerBound(v), trail->UpperBound(v)) << v; return trail->LowerBound(v).value(); }; } @@ -983,6 +1026,30 @@ inline std::function Equality(IntegerLiteral i, Literal l) { }; } +// TODO(user): This is one of the rare case where it is better to use Equality() +// rather than two Implications(). Maybe we should modify our internal +// implementation to use half-reified encoding? that is do not propagate the +// direction integer-bound => literal, but just literal => integer-bound? This +// is the same as using different underlying variable for an integer literal and +// its negation. +inline std::function Implication(Literal l, IntegerLiteral i) { + return [=](Model* model) { + IntegerTrail* integer_trail = model->GetOrCreate(); + if (i.Bound() <= integer_trail->LowerBound(i.Var())) { + // Always true! nothing to do. + } else if (i.Bound() > integer_trail->UpperBound(i.Var())) { + // Always false. + model->Add(ClauseConstraint({l.Negated()})); + } else { + // TODO(user): Double check what happen when we associate a trivially + // true or false literal. This applies to Equality() too. + IntegerEncoder* encoder = model->GetOrCreate(); + const Literal current = encoder->GetOrCreateAssociatedLiteral(i); + model->Add(Implication(l, current)); + } + }; +} + // in_interval <=> v in [lb, ub]. inline std::function ReifiedInInterval(IntegerVariable v, int64 lb, int64 ub, @@ -1031,13 +1098,46 @@ FullyEncodeVariable(IntegerVariable var) { } // A wrapper around SatSolver::Solve that handles integer variable with lazy -// encoding. Repeatedly calls SatSolver::Solve on the model, and instantiates -// literals for non-fixed variables in vars_with_lazy_encoding, until either all -// variables are fixed to a single value, or the model is proved UNSAT. Returns -// the status of the last call to SatSolver::Solve(). +// encoding. Repeatedly calls SatSolver::Solve() on the model until the given +// next_decision() function return kNoLiteralIndex or the model is proved to +// be UNSAT. +// +// Returns the status of the last call to SatSolver::Solve(). +// +// Note that the next_decision() function must always return an unassigned +// literal or kNoLiteralIndex to end the search. SatSolver::Status SolveIntegerProblemWithLazyEncoding( const std::vector& assumptions, - const std::vector& vars_with_lazy_encoding, Model* model); + const std::function& next_decision, Model* model); + +// Shortcut for SolveIntegerProblemWithLazyEncoding() when there is no +// assumption and we consider all variables in their index order for the next +// search decision. +SatSolver::Status SolveIntegerProblemWithLazyEncoding(Model* model); + +// Decision heuristic for SolveIntegerProblemWithLazyEncoding(). Returns a +// function that will return the literal corresponding to the fact that the +// first currently non-fixed variable value is <= its min. The function will +// return kNoLiteralIndex if all the given variables are fixed. +// +// Note that this function will create the associated literal if needed. +std::function FirstUnassignedVarAtItsMinHeuristic( + const std::vector& vars, Model* model); + +// Decision heuristic for SolveIntegerProblemWithLazyEncoding(). Like +// FirstUnassignedVarAtItsMinHeuristic() but the function will return the +// literal corresponding to the fact that the currently non-assigned variable +// with the lowest min has a value <= this min. +std::function UnassignedVarWithLowestMinAtItsMinHeuristic( + const std::vector& vars, Model* model); + +// Same as ExcludeCurrentSolutionAndBacktrack() but this version works for an +// integer problem with optional variables. The issue is that an optional +// variable that is ignored can basically take any value, and we don't really +// want to enumerate them. This function should exclude all solutions where +// only the ignored variable values change. +std::function +ExcludeCurrentSolutionWithoutIgnoredVariableAndBacktrack(); } // namespace sat } // namespace operations_research diff --git a/src/sat/integer_expr.h b/src/sat/integer_expr.h index aff29b819b..afd2d686b9 100644 --- a/src/sat/integer_expr.h +++ b/src/sat/integer_expr.h @@ -173,7 +173,7 @@ inline std::function WeightedSumLowerOrEqual( const std::vector& vars, const VectorInt& coefficients, int64 upper_bound) { // Special cases. - CHECK_GE(vars.size(), 1) << "Should be encoded differently."; + CHECK_GE(vars.size(), 1); if (vars.size() == 1) { CHECK_NE(coefficients[0], 0); if (coefficients[0] > 0) { @@ -237,8 +237,16 @@ inline std::function ConditionalWeightedSumLowerOrEqual( // Special cases. CHECK_GE(vars.size(), 1); if (vars.size() == 1) { - LOG(INFO) << "ConditionalWeightedSumLowerOrEqual() with 1 term. Should " - "have been presolved!"; + CHECK_NE(coefficients[0], 0); + if (coefficients[0] > 0) { + return Implication( + is_le, IntegerLiteral::LowerOrEqual( + vars[0], IntegerValue(upper_bound / coefficients[0]))); + } else { + return Implication( + is_le, IntegerLiteral::GreaterOrEqual( + vars[0], IntegerValue(upper_bound / coefficients[0]))); + } } if (vars.size() == 2 && (coefficients[0] == 1 || coefficients[0] == -1) && (coefficients[1] == 1 || coefficients[1] == -1)) { diff --git a/src/sat/intervals.cc b/src/sat/intervals.cc index bac63299c2..31d52b72b7 100644 --- a/src/sat/intervals.cc +++ b/src/sat/intervals.cc @@ -13,6 +13,8 @@ #include "sat/intervals.h" +#include "util/sort.h" + namespace operations_research { namespace sat { @@ -57,5 +59,194 @@ IntervalVariable IntervalsRepository::CreateInterval(IntegerVariable start, return i; } +SchedulingConstraintHelper::SchedulingConstraintHelper( + const std::vector& tasks, Trail* trail, + IntegerTrail* integer_trail, IntervalsRepository* repository) + : trail_(trail), + integer_trail_(integer_trail), + current_time_direction_(true) { + start_vars_.clear(); + end_vars_.clear(); + minus_end_vars_.clear(); + minus_start_vars_.clear(); + duration_vars_.clear(); + fixed_durations_.clear(); + reason_for_presence_.clear(); + for (const IntervalVariable i : tasks) { + if (repository->IsOptional(i)) { + reason_for_presence_.push_back(repository->IsPresentLiteral(i).Index()); + } else { + reason_for_presence_.push_back(kNoLiteralIndex); + } + if (repository->SizeVar(i) == kNoIntegerVariable) { + duration_vars_.push_back(kNoIntegerVariable); + fixed_durations_.push_back(repository->MinSize(i)); + } else { + duration_vars_.push_back(repository->SizeVar(i)); + fixed_durations_.push_back(IntegerValue(0)); + } + start_vars_.push_back(repository->StartVar(i)); + end_vars_.push_back(repository->EndVar(i)); + minus_start_vars_.push_back(NegationOf(repository->StartVar(i))); + minus_end_vars_.push_back(NegationOf(repository->EndVar(i))); + } + + const int num_tasks = start_vars_.size(); + task_by_increasing_min_start_.resize(num_tasks); + task_by_increasing_min_end_.resize(num_tasks); + task_by_decreasing_max_start_.resize(num_tasks); + task_by_decreasing_max_end_.resize(num_tasks); + for (int t = 0; t < num_tasks; ++t) { + task_by_increasing_min_start_[t].task_index = t; + task_by_increasing_min_end_[t].task_index = t; + task_by_decreasing_max_start_[t].task_index = t; + task_by_decreasing_max_end_[t].task_index = t; + } +} + +void SchedulingConstraintHelper::SetTimeDirection(bool is_forward) { + if (current_time_direction_ == is_forward) return; + current_time_direction_ = is_forward; + + std::swap(start_vars_, minus_end_vars_); + std::swap(end_vars_, minus_start_vars_); + std::swap(task_by_increasing_min_start_, task_by_decreasing_max_end_); + std::swap(task_by_increasing_min_end_, task_by_decreasing_max_start_); +} + +const std::vector& +SchedulingConstraintHelper::TaskByIncreasingStartMin() { + const int num_tasks = NumTasks(); + for (int i = 0; i < num_tasks; ++i) { + TaskTime& ref = task_by_increasing_min_start_[i]; + ref.time = StartMin(ref.task_index); + } + IncrementalSort(task_by_increasing_min_start_.begin(), + task_by_increasing_min_start_.end()); + return task_by_increasing_min_start_; +} + +const std::vector& +SchedulingConstraintHelper::TaskByIncreasingEndMin() { + const int num_tasks = NumTasks(); + for (int i = 0; i < num_tasks; ++i) { + TaskTime& ref = task_by_increasing_min_end_[i]; + ref.time = EndMin(ref.task_index); + } + IncrementalSort(task_by_increasing_min_end_.begin(), + task_by_increasing_min_end_.end()); + return task_by_increasing_min_end_; +} + +const std::vector& +SchedulingConstraintHelper::TaskByDecreasingStartMax() { + const int num_tasks = NumTasks(); + for (int i = 0; i < num_tasks; ++i) { + TaskTime& ref = task_by_decreasing_max_start_[i]; + ref.time = StartMax(ref.task_index); + } + IncrementalSort(task_by_decreasing_max_start_.begin(), + task_by_decreasing_max_start_.end(), + std::greater()); + return task_by_decreasing_max_start_; +} + +const std::vector& +SchedulingConstraintHelper::TaskByDecreasingEndMax() { + const int num_tasks = NumTasks(); + for (int i = 0; i < num_tasks; ++i) { + TaskTime& ref = task_by_decreasing_max_end_[i]; + ref.time = EndMax(ref.task_index); + } + IncrementalSort(task_by_decreasing_max_end_.begin(), + task_by_decreasing_max_end_.end(), std::greater()); + return task_by_decreasing_max_end_; +} + +bool SchedulingConstraintHelper::PushIntegerLiteral(IntegerLiteral bound) { + return integer_trail_->Enqueue(bound, literal_reason_, integer_reason_); +} + +bool SchedulingConstraintHelper::IncreaseStartMin(int t, + IntegerValue new_min_start) { + if (!integer_trail_->Enqueue( + IntegerLiteral::GreaterOrEqual(start_vars_[t], new_min_start), + literal_reason_, integer_reason_)) { + return false; + } + + // Skip if the task is now absent. + if (IsAbsent(t)) return true; + + // We propagate right away the new end-min lower-bound we have. + const IntegerValue min_end_lb = new_min_start + DurationMin(t); + if (EndMin(t) < min_end_lb) { + ClearReason(); + AddStartMinReason(t, new_min_start); + AddDurationMinReason(t); + if (!integer_trail_->Enqueue( + IntegerLiteral::GreaterOrEqual(end_vars_[t], min_end_lb), + literal_reason_, integer_reason_)) { + return false; + } + } + + return true; +} + +bool SchedulingConstraintHelper::DecreaseEndMax(int t, + IntegerValue new_max_end) { + if (!integer_trail_->Enqueue( + IntegerLiteral::LowerOrEqual(end_vars_[t], new_max_end), + literal_reason_, integer_reason_)) { + return false; + } + + // Skip if the task is now absent. + if (IsAbsent(t)) return true; + + // We propagate right away the new start-max upper-bound we have. + const IntegerValue max_start_ub = new_max_end - DurationMin(t); + if (StartMax(t) > max_start_ub) { + ClearReason(); + AddEndMaxReason(t, new_max_end); + AddDurationMinReason(t); + if (!integer_trail_->Enqueue( + IntegerLiteral::LowerOrEqual(start_vars_[t], max_start_ub), + literal_reason_, integer_reason_)) { + return false; + } + } + + return true; +} + +void SchedulingConstraintHelper::PushTaskAbsence(int t) { + DCHECK_NE(reason_for_presence_[t], kNoLiteralIndex); + DCHECK(!IsPresent(t)); + DCHECK(!IsAbsent(t)); + integer_trail_->EnqueueLiteral(Literal(reason_for_presence_[t]).Negated(), + literal_reason_, integer_reason_); +} + +bool SchedulingConstraintHelper::ReportConflict() { + return integer_trail_->ReportConflict(literal_reason_, integer_reason_); +} + +void SchedulingConstraintHelper::WatchAllTasks( + int id, GenericLiteralWatcher* watcher) const { + const int num_tasks = start_vars_.size(); + for (int t = 0; t < num_tasks; ++t) { + watcher->WatchIntegerVariable(start_vars_[t], id); + watcher->WatchIntegerVariable(end_vars_[t], id); + if (duration_vars_[t] != kNoIntegerVariable) { + watcher->WatchLowerBound(duration_vars_[t], id); + } + if (!IsPresent(t) && !IsAbsent(t)) { + watcher->WatchLiteral(Literal(reason_for_presence_[t]), id); + } + } +} + } // namespace sat } // namespace operations_research diff --git a/src/sat/intervals.h b/src/sat/intervals.h index 5b7ad6e983..30fd1316be 100644 --- a/src/sat/intervals.h +++ b/src/sat/intervals.h @@ -109,6 +109,204 @@ class IntervalsRepository { DISALLOW_COPY_AND_ASSIGN(IntervalsRepository); }; +// Helper class shared by the propagators that manage a given list of tasks. +// +// One of the main advantage of this class is that it allows to share the +// vectors of tasks sorted by various criteria between propagator for a faster +// code. +class SchedulingConstraintHelper { + public: + // All the functions below refer to a task by its index t in the tasks + // vector given at construction. + SchedulingConstraintHelper(const std::vector& tasks, + Trail* trail, IntegerTrail* integer_trail, + IntervalsRepository* intervals_repository); + + // Returns the number of task. + int NumTasks() const { return start_vars_.size(); } + + // Sets the time direction to either forward/backward. This will impact all + // the functions below. + void SetTimeDirection(bool is_forward); + + // Helpers for the current bounds on the current task time window. + // [(duration-min) ... (duration-min)] + // ^ ^ ^ ^ + // start-min end-min start-max end-max + // + // Note that for tasks with variable durations, we don't necessarily have + // duration-min between the the XXX-min and XXX-max value. + IntegerValue DurationMin(int t) const; + IntegerValue StartMin(int t) const; + IntegerValue StartMax(int t) const; + IntegerValue EndMin(int t) const; + IntegerValue EndMax(int t) const; + + // Returns true if the corresponding fact is known for sure. A normal task is + // always present. For optional task for which the presence is still unknown, + // both of these function will return false. + bool IsPresent(int t) const; + bool IsAbsent(int t) const; + + // Sorts and returns the tasks in corresponding order at the time of the call. + // Note that we do not mean strictly-increasing/strictly-decreasing, there + // will be duplicate time values in these vectors. + // + // TODO(user): we could merge the first loop of IncrementalSort() with the + // loop that fill TaskTime.time at each call. + struct TaskTime { + int task_index; + IntegerValue time; + bool operator<(TaskTime other) const { return time < other.time; } + bool operator>(TaskTime other) const { return time > other.time; } + }; + const std::vector& TaskByIncreasingStartMin(); + const std::vector& TaskByIncreasingEndMin(); + const std::vector& TaskByDecreasingStartMax(); + const std::vector& TaskByDecreasingEndMax(); + + // Functions to clear and then set the current reason. + void ClearReason(); + void AddPresenceReason(int t); + void AddDurationMinReason(int t); + void AddStartMinReason(int t, IntegerValue lower_bound); + void AddStartMaxReason(int t, IntegerValue upper_bound); + void AddEndMinReason(int t, IntegerValue lower_bound); + void AddEndMaxReason(int t, IntegerValue upper_bound); + + // It is also possible to directly manipulates the underlying reason vectors + // that will be used when pushing something. + std::vector* MutableLiteralReason() { return &literal_reason_; } + std::vector* MutableIntegerReason() { + return &integer_reason_; + } + + // Push something using the current reason. Note that IncreaseStartMin() will + // also increase the end-min, and DecreaseEndMax() will also decrease the + // start-max. + bool IncreaseStartMin(int t, IntegerValue new_min_start); + bool DecreaseEndMax(int t, IntegerValue new_max_end); + void PushTaskAbsence(int t); + bool PushIntegerLiteral(IntegerLiteral bound); + bool ReportConflict(); + + // Returns the underlying integer variables. + const std::vector& StartVars() const { return start_vars_; } + const std::vector& EndVars() const { return end_vars_; } + + // Registers the given propagator id to be called if any of the tasks + // in this class change. + void WatchAllTasks(int id, GenericLiteralWatcher* watcher) const; + + private: + Trail* trail_; + IntegerTrail* integer_trail_; + + // All the underlying variables of the tasks. + // The vectors are indexed by the task index t. + std::vector start_vars_; + std::vector end_vars_; + std::vector duration_vars_; + std::vector fixed_durations_; + std::vector reason_for_presence_; + + // The current direction of time, true for forward, false for backward. + bool current_time_direction_; + + // The negation of the start/end variable so that SetTimeDirection() + // can do its job in O(1) instead of calling NegationOf() on each entry. + std::vector minus_start_vars_; + std::vector minus_end_vars_; + + // Sorted vectors returned by the TasksBy*() functions. + std::vector task_by_increasing_min_start_; + std::vector task_by_increasing_min_end_; + std::vector task_by_decreasing_max_start_; + std::vector task_by_decreasing_max_end_; + + // Reason vectors. + std::vector literal_reason_; + std::vector integer_reason_; +}; + +// ============================================================================= +// SchedulingConstraintHelper inlined functions. +// ============================================================================= + +inline IntegerValue SchedulingConstraintHelper::DurationMin(int t) const { + return duration_vars_[t] == kNoIntegerVariable + ? fixed_durations_[t] + : integer_trail_->LowerBound(duration_vars_[t]); +} + +inline IntegerValue SchedulingConstraintHelper::StartMin(int t) const { + return integer_trail_->LowerBound(start_vars_[t]); +} + +inline IntegerValue SchedulingConstraintHelper::StartMax(int t) const { + return integer_trail_->UpperBound(start_vars_[t]); +} + +inline IntegerValue SchedulingConstraintHelper::EndMin(int t) const { + return integer_trail_->LowerBound(end_vars_[t]); +} + +inline IntegerValue SchedulingConstraintHelper::EndMax(int t) const { + return integer_trail_->UpperBound(end_vars_[t]); +} + +inline bool SchedulingConstraintHelper::IsPresent(int t) const { + if (reason_for_presence_[t] == kNoLiteralIndex) return true; + return trail_->Assignment().LiteralIsTrue(Literal(reason_for_presence_[t])); +} + +inline bool SchedulingConstraintHelper::IsAbsent(int t) const { + if (reason_for_presence_[t] == kNoLiteralIndex) return false; + return trail_->Assignment().LiteralIsFalse(Literal(reason_for_presence_[t])); +} + +inline void SchedulingConstraintHelper::ClearReason() { + integer_reason_.clear(); + literal_reason_.clear(); +} + +inline void SchedulingConstraintHelper::AddPresenceReason(int t) { + if (reason_for_presence_[t] != kNoLiteralIndex) { + literal_reason_.push_back(Literal(reason_for_presence_[t]).Negated()); + } +} + +inline void SchedulingConstraintHelper::AddDurationMinReason(int t) { + if (duration_vars_[t] != kNoIntegerVariable) { + integer_reason_.push_back( + integer_trail_->LowerBoundAsLiteral(duration_vars_[t])); + } +} + +inline void SchedulingConstraintHelper::AddStartMinReason( + int t, IntegerValue lower_bound) { + integer_reason_.push_back( + IntegerLiteral::GreaterOrEqual(start_vars_[t], lower_bound)); +} + +inline void SchedulingConstraintHelper::AddStartMaxReason( + int t, IntegerValue upper_bound) { + integer_reason_.push_back( + IntegerLiteral::LowerOrEqual(start_vars_[t], upper_bound)); +} + +inline void SchedulingConstraintHelper::AddEndMinReason( + int t, IntegerValue lower_bound) { + integer_reason_.push_back( + IntegerLiteral::GreaterOrEqual(end_vars_[t], lower_bound)); +} + +inline void SchedulingConstraintHelper::AddEndMaxReason( + int t, IntegerValue upper_bound) { + integer_reason_.push_back( + IntegerLiteral::LowerOrEqual(end_vars_[t], upper_bound)); +} + // ============================================================================= // Model based functions. // ============================================================================= diff --git a/src/sat/linear_programming_constraint.cc b/src/sat/linear_programming_constraint.cc new file mode 100644 index 0000000000..48b72501f4 --- /dev/null +++ b/src/sat/linear_programming_constraint.cc @@ -0,0 +1,382 @@ +// Copyright 2010-2014 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "sat/linear_programming_constraint.h" + +#include "base/commandlineflags.h" +#include "util/time_limit.h" + +// TODO(user): remove the option once we know which algo work best. +DEFINE_bool(lp_constraint_use_dual_ray, true, + "If true, use the dual simplex and exploit the dual ray when the " + "problem is DUAL_UNBOUNDED as a reason rather than " + "solving a custom feasibility LP first."); + +namespace operations_research { +namespace sat { + +LinearProgrammingConstraint::LinearProgrammingConstraint( + IntegerTrail* integer_trail) + : integer_trail_(integer_trail) { + if (!FLAGS_lp_constraint_use_dual_ray) { + // The violation_sum_ variable will be the sum of constraints' violation. + violation_sum_constraint_ = lp_data_.CreateNewConstraint(); + lp_data_.SetConstraintBounds(violation_sum_constraint_, 0.0, 0.0); + violation_sum_ = lp_data_.CreateNewVariable(); + lp_data_.SetVariableBounds(violation_sum_, 0.0, + std::numeric_limits::infinity()); + lp_data_.SetCoefficient(violation_sum_constraint_, violation_sum_, -1.0); + } + + // Tweak the default parameters to make the solve incremental. + glop::GlopParameters parameters; + parameters.set_use_dual_simplex(true); + simplex_.SetParameters(parameters); +} + +LinearProgrammingConstraint::ConstraintIndex +LinearProgrammingConstraint::CreateNewConstraint(double lb, double ub) { + DCHECK(!lp_constraint_is_registered_); + const ConstraintIndex ct = lp_data_.CreateNewConstraint(); + lp_data_.SetConstraintBounds(ct, lb, ub); + return ct; +} + +glop::ColIndex LinearProgrammingConstraint::GetOrCreateMirrorVariable( + IntegerVariable ivar) { + const auto it = integer_variable_to_index_.find(ivar); + if (it == integer_variable_to_index_.end()) { + integer_variable_to_index_[ivar] = integer_variables_.size(); + integer_variables_.push_back(ivar); + mirror_lp_variables_.push_back(lp_data_.CreateNewVariable()); + lp_solution_.push_back(std::numeric_limits::infinity()); + } + return mirror_lp_variables_[integer_variable_to_index_[ivar]]; +} + +void LinearProgrammingConstraint::SetCoefficient(ConstraintIndex ct, + IntegerVariable ivar, + double coefficient) { + CHECK(!lp_constraint_is_registered_); + glop::ColIndex cvar = GetOrCreateMirrorVariable(ivar); + lp_data_.SetCoefficient(ct, cvar, coefficient); +} + +void LinearProgrammingConstraint::SetObjective(IntegerVariable ivar, + bool is_minimization) { + CHECK(!lp_constraint_is_registered_); + CHECK(!objective_is_defined_) << "Objective was set more than once."; + objective_is_defined_ = true; + objective_cp_ = ivar; + objective_lp_ = GetOrCreateMirrorVariable(ivar); + objective_is_minimization_ = is_minimization; +} + +void LinearProgrammingConstraint::RegisterWith(GenericLiteralWatcher* watcher) { + DCHECK(!lp_constraint_is_registered_); + lp_constraint_is_registered_ = true; + + lp_data_.Scale(&scaler_); + + // Note that we set the objective AFTER the scaling. + if (objective_is_defined_) { + lp_data_.SetObjectiveCoefficient(objective_lp_, 1.0); + lp_data_.SetMaximizationProblem(!objective_is_minimization_); + } + + if (!FLAGS_lp_constraint_use_dual_ray) { + // Add all the individual violation variables. Note that it is important + // to do that AFTER the scaling so that each constraint is considered on the + // same footing regarding a violation. + // + // Note that scaler_.col_scale() will returns a value of 1.0 for these new + // variables. + // + // TODO(user): See if it is possible to reuse the feasibility code of the + // simplex that do not need to create these extra variables. + // + // TODO(user): It might be better (smaller reasons) to to check the maximum + // of the individual constraint violation rather than the sum. + const double infinity = std::numeric_limits::infinity(); + for (glop::RowIndex row(0); row < lp_data_.num_constraints(); ++row) { + if (row == violation_sum_constraint_) continue; + const glop::Fractional lb = lp_data_.constraint_lower_bounds()[row]; + const glop::Fractional ub = lp_data_.constraint_upper_bounds()[row]; + if (lb != -infinity) { + const glop::ColIndex violation_lb = lp_data_.CreateNewVariable(); + lp_data_.SetVariableBounds(violation_lb, 0.0, infinity); + lp_data_.SetCoefficient(violation_sum_constraint_, violation_lb, 1.0); + lp_data_.SetCoefficient(row, violation_lb, 1.0); + } + if (ub != infinity) { + const glop::ColIndex violation_ub = lp_data_.CreateNewVariable(); + lp_data_.SetVariableBounds(violation_ub, 0.0, infinity); + lp_data_.SetCoefficient(violation_sum_constraint_, violation_ub, 1.0); + lp_data_.SetCoefficient(row, violation_ub, -1.0); + } + } + } + + lp_data_.AddSlackVariablesWhereNecessary(false); + + const int watcher_id = watcher->Register(this); + const int num_vars = integer_variables_.size(); + for (int i = 0; i < num_vars; i++) { + watcher->WatchIntegerVariable(integer_variables_[i], watcher_id, i); + } + watcher->SetPropagatorPriority(watcher_id, 2); +} + +// Check whether the change breaks the current LP solution. +// Call Propagate() only if it does. +bool LinearProgrammingConstraint::IncrementalPropagate( + const std::vector& watch_indices) { + for (const int index : watch_indices) { + const double lb = static_cast( + integer_trail_->LowerBound(integer_variables_[index]).value()); + const double ub = static_cast( + integer_trail_->UpperBound(integer_variables_[index]).value()); + const double value = lp_solution_[index]; + if (value < lb - kEpsilon || value > ub + kEpsilon) return Propagate(); + } + return true; +} + +glop::Fractional LinearProgrammingConstraint::GetVariableValueAtCpScale( + glop::ColIndex var) { + return simplex_.GetVariableValue(var) / scaler_.col_scale(var); +} + +bool LinearProgrammingConstraint::Propagate() { + // Copy CP state into LP state. + const int num_vars = integer_variables_.size(); + for (int i = 0; i < num_vars; i++) { + const IntegerVariable cp_var = integer_variables_[i]; + const double lb = + static_cast(integer_trail_->LowerBound(cp_var).value()); + const double ub = + static_cast(integer_trail_->UpperBound(cp_var).value()); + const double factor = scaler_.col_scale(mirror_lp_variables_[i]); + lp_data_.SetVariableBounds(mirror_lp_variables_[i], lb * factor, + ub * factor); + } + + if (!FLAGS_lp_constraint_use_dual_ray) { + if (objective_is_defined_) { + lp_data_.SetObjectiveCoefficient(objective_lp_, 0.0); + } + lp_data_.SetObjectiveCoefficient(violation_sum_, 1.0); + lp_data_.SetVariableBounds(violation_sum_, 0.0, + std::numeric_limits::infinity()); + lp_data_.SetMaximizationProblem(false); + + // Feasibility deductions. + const auto status = simplex_.Solve(lp_data_, TimeLimit::Infinite().get()); + CHECK(status.ok()) << "LinearProgrammingConstraint encountered an error: " + << status.error_message(); + CHECK_EQ(simplex_.GetProblemStatus(), glop::ProblemStatus::OPTIMAL) + << "simplex Solve() should return optimal, but it returned " + << simplex_.GetProblemStatus(); + + if (simplex_.GetVariableValue(violation_sum_) > kEpsilon) { // infeasible. + FillIntegerReason(1.0); + return integer_trail_->ReportConflict(integer_reason_); + } + + // Reduced cost strengthening for feasibility. + ReducedCostStrengtheningDeductions(1.0, 0.0); + if (!deductions_.empty()) { + FillIntegerReason(1.0); + for (const IntegerLiteral deduction : deductions_) { + if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) { + return false; + } + } + } + + // Revert to the real problem objective and save current solution. + lp_data_.SetVariableBounds(violation_sum_, 0.0, 0.0); + lp_data_.SetObjectiveCoefficient(violation_sum_, 0.0); + if (objective_is_defined_) { + lp_data_.SetObjectiveCoefficient(objective_lp_, 1.0); + lp_data_.SetMaximizationProblem(!objective_is_minimization_); + } + for (int i = 0; i < num_vars; i++) { + lp_solution_[i] = GetVariableValueAtCpScale(mirror_lp_variables_[i]); + } + + // We currently ignore the objective and return right away when we don't + // use the dual ray as an infeasibility reason. + return true; + } + + const auto status = simplex_.Solve(lp_data_, TimeLimit::Infinite().get()); + CHECK(status.ok()) << "LinearProgrammingConstraint encountered an error: " + << status.error_message(); + + // A dual-unbounded problem is infeasible. We use the dual ray reason. + if (simplex_.GetProblemStatus() == glop::ProblemStatus::DUAL_UNBOUNDED) { + FillDualRayReason(); + return integer_trail_->ReportConflict(integer_reason_); + } + + // Optimality deductions if problem has an objective. + if (objective_is_defined_ && + simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) { + const double objective_cp_lb = + static_cast(integer_trail_->LowerBound(objective_cp_).value()); + const double objective_cp_ub = + static_cast(integer_trail_->UpperBound(objective_cp_).value()); + + // Try to filter optimal objective value. + const double objective_value = GetVariableValueAtCpScale(objective_lp_); + if (objective_is_minimization_) { + const double new_lb = std::ceil(objective_value - kEpsilon); + if (objective_cp_lb < new_lb) { + const IntegerValue new_int_lb(static_cast(new_lb)); + FillIntegerReason(1.0); + const IntegerLiteral deduction = + IntegerLiteral::GreaterOrEqual(objective_cp_, new_int_lb); + if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) { + return false; + } + } + } else { + const double new_ub = std::floor(objective_value + kEpsilon); + if (objective_cp_ub > new_ub) { + const IntegerValue new_int_ub(static_cast(new_ub)); + FillIntegerReason(-1.0); + const IntegerLiteral deduction = + IntegerLiteral::LowerOrEqual(objective_cp_, new_int_ub); + if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) { + return false; + } + } + } + + // Reduced cost strengthening. + const double objective_slack = objective_is_minimization_ + ? objective_cp_ub - objective_value + : objective_value - objective_cp_lb; + const double objective_direction = objective_is_minimization_ ? 1.0 : -1.0; + ReducedCostStrengtheningDeductions( + objective_direction, + objective_slack * scaler_.col_scale(objective_lp_)); + + if (!deductions_.empty()) { + FillIntegerReason(objective_direction); + + // Add the opposite bound of the variable used for strengthening. + const IntegerLiteral opposite_bound = + objective_is_minimization_ + ? integer_trail_->UpperBoundAsLiteral(objective_cp_) + : integer_trail_->LowerBoundAsLiteral(objective_cp_); + integer_reason_.push_back(opposite_bound); + + for (const IntegerLiteral deduction : deductions_) { + if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) { + return false; + } + } + } + } + + // Copy current LP solution. + for (int i = 0; i < num_vars; i++) { + lp_solution_[i] = GetVariableValueAtCpScale(mirror_lp_variables_[i]); + } + + return true; +} + +void LinearProgrammingConstraint::FillIntegerReason(double direction) { + integer_reason_.clear(); + + const int num_vars = integer_variables_.size(); + for (int i = 0; i < num_vars; i++) { + // TODO(user): try to extend the bounds that are put in the + // explanation of feasibility: can we compute bounds of variables for which + // the objective would still not be low/high enough for the problem to be + // feasible? If the violation minimum is 10 and a variable has rc 1, + // then decreasing it by 9 would still leave the problem infeasible. + // Using this could allow to generalize some explanations. + const double rc = + simplex_.GetReducedCost(mirror_lp_variables_[i]) * direction; + if (rc > kEpsilon) { + integer_reason_.push_back( + integer_trail_->LowerBoundAsLiteral(integer_variables_[i])); + } else if (rc < -kEpsilon) { + integer_reason_.push_back( + integer_trail_->UpperBoundAsLiteral(integer_variables_[i])); + } + } +} + +void LinearProgrammingConstraint::FillDualRayReason() { + integer_reason_.clear(); + + const int num_vars = integer_variables_.size(); + for (int i = 0; i < num_vars; i++) { + // TODO(user): Like for FillIntegerReason(), the bounds could be + // extended here. Actually, the "dual ray cost updates" is the reduced cost + // of an optimal solution if we were optimizing one direction of one basic + // variable. The simplex_ interface would need to be slightly extended to + // retrieve the basis column in question and the variable values though. + const double rc = + simplex_.GetDualRayRowCombination()[mirror_lp_variables_[i]]; + if (rc > kEpsilon) { + integer_reason_.push_back( + integer_trail_->LowerBoundAsLiteral(integer_variables_[i])); + } else if (rc < -kEpsilon) { + integer_reason_.push_back( + integer_trail_->UpperBoundAsLiteral(integer_variables_[i])); + } + } +} + +void LinearProgrammingConstraint::ReducedCostStrengtheningDeductions( + double direction, double lp_objective_delta) { + deductions_.clear(); + + const int num_vars = integer_variables_.size(); + for (int i = 0; i < num_vars; i++) { + const IntegerVariable cp_var = integer_variables_[i]; + const glop::ColIndex lp_var = mirror_lp_variables_[i]; + const double rc = simplex_.GetReducedCost(lp_var) * direction; + const double value = simplex_.GetVariableValue(lp_var); + const double lp_other_bound = value + lp_objective_delta / rc; + const double cp_other_bound = lp_other_bound / scaler_.col_scale(lp_var); + + if (rc > kEpsilon) { + const double ub = + static_cast(integer_trail_->UpperBound(cp_var).value()); + const double new_ub = std::floor(cp_other_bound + kEpsilon); + if (new_ub < ub) { + const IntegerValue new_ub_int(static_cast(new_ub)); + deductions_.push_back(IntegerLiteral::LowerOrEqual(cp_var, new_ub_int)); + } + } else if (rc < -kEpsilon) { + const double lb = + static_cast(integer_trail_->LowerBound(cp_var).value()); + const double new_lb = std::ceil(cp_other_bound - kEpsilon); + if (new_lb > lb) { + const IntegerValue new_lb_int(static_cast(new_lb)); + deductions_.push_back( + IntegerLiteral::GreaterOrEqual(cp_var, new_lb_int)); + } + } + } +} + +} // namespace sat +} // namespace operations_research diff --git a/src/sat/linear_programming_constraint.h b/src/sat/linear_programming_constraint.h new file mode 100644 index 0000000000..8b372c5828 --- /dev/null +++ b/src/sat/linear_programming_constraint.h @@ -0,0 +1,178 @@ +// Copyright 2010-2014 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_SAT_LINEAR_PROGRAMMING_CONSTRAINT_H_ +#define OR_TOOLS_SAT_LINEAR_PROGRAMMING_CONSTRAINT_H_ + +#include +#include + +#include "glop/revised_simplex.h" +#include "lp_data/lp_data.h" +#include "sat/integer.h" +#include "sat/model.h" +#include "sat/sat_base.h" + +namespace operations_research { +namespace sat { + +// A SAT constraint that enforces a set of linear inequality constraints on +// integer variables using an LP solver. +// +// The propagator uses glop's revised simplex for feasibility and propagation. +// It uses the Reduced Cost Strengthening technique, a classic in mixed integer +// programming, for instance see the thesis of Tobias Achterberg, +// "Constraint Integer Programming", sections 7.7 and 8.8, algorithm 7.11. +// http://nbn-resolving.de/urn:nbn:de:0297-zib-11129 +// Feasibility propagation is done with this technique by setting total +// constraint violation as an objective, i.e. by dualizing all constraints. +// Per-constraint bounds propagation is NOT done by this constraint, +// it should be done by redundant constraints, as reduced cost propagation +// may miss some filtering. +// +// Workflow: create a LinearProgrammingConstraint instance, make linear +// inequality constraints, call RegisterWith() to finalize the set of linear +// constraints. A linear constraint a x + b y + c z <= k, with x y z +// IntegerVariables, can be created by calling: +// auto ct = lp->CreateNewConstraint(-std::numeric_limits::infinity(), +// k); +// lp->SetCoefficient(ct, x, a); +// lp->SetCoefficient(ct, y, b); +// lp->SetCoefficient(ct, z, c); +// +// Note that this constraint works with double floating-point numbers, so one +// could be worried that it may filter too much in case of precision issues. +// However, the underlying LP solver reports infeasibility only if the problem +// is still infeasible by relaxing the bounds by some small relative value. +// Thus the constraint will tend to filter less than it could, not the opposite. +// +// TODO(user): Work with scaled version of the model, maybe by using +// LPSolver instead of RevisedSimplex. +class LinearProgrammingConstraint : public PropagatorInterface { + public: + typedef glop::RowIndex ConstraintIndex; + + explicit LinearProgrammingConstraint(IntegerTrail* integer_trail); + + // Creates a LinearProgrammingConstraint for templated GetOrCreate idiom. + static LinearProgrammingConstraint* CreateInModel(Model* model) { + IntegerTrail* trail = model->GetOrCreate(); + LinearProgrammingConstraint* constraint = + new LinearProgrammingConstraint(trail); + model->TakeOwnership(constraint); + return constraint; + } + + // User API, see header description. + ConstraintIndex CreateNewConstraint(double lb, double ub); + + void SetCoefficient(ConstraintIndex ct, IntegerVariable ivar, + double coefficient); + + // TODO(user): Allow Literals to appear in linear constraints. + // TODO(user): Calling SetCoefficient() twice on the same + // (constraint, variable) pair will overwrite coefficients where accumulating + // them might be desired, this is a common mistake, change API. + + // Objective may or may not be defined. It can be defined only once, + // must be exactly one IntegerVariable, and can be either + // minimized (is_minimization = true) or maximized (is_minimization = false). + // TODO(user): change API for always minimization, so that + // maximization(var) = minimization(Negation(var)). + void SetObjective(IntegerVariable ivar, bool is_minimization); + + // PropagatorInterface API. + bool Propagate() override; + + bool IncrementalPropagate(const std::vector& watch_indices) override; + + void RegisterWith(GenericLiteralWatcher* watcher); + + private: + // Generates a set of IntegerLiterals explaining why the best solution can not + // be improved using reduced costs. This is used to generate explanations for + // both infeasibility and bounds deductions. The direction variable should be + // 1.0 if the last Solve() was a minimization, -1.0 if it was a maximization. + void FillIntegerReason(double direction); + + // Same as FillIntegerReason() but for the case of a DUAL_UNBOUNDED problem. + // This exploit the dual ray as a reason for the primal infeasiblity. + void FillDualRayReason(); + + // Fills the deductions vector with reduced cost deductions that can be made + // from the current state of the LP solver. This should be called after + // Solve(): if the optimization was a minimization, the direction variable + // should be 1.0 and lp_objective_delta the objective's upper bound minus the + // optimal; if the optimization was a maximization, direction should be -1.0 + // and lp_objective_delta the optimal minus the objective's lower bound. + void ReducedCostStrengtheningDeductions(double direction, + double lp_objective_delta); + + // Gets or creates an LP variable that mirrors a CP variable. + // TODO(user): only accept positive variables to prevent having different + // LP variables for the same CP variable. + glop::ColIndex GetOrCreateMirrorVariable(IntegerVariable ivar); + + // Returns the variable value on the same scale as the CP variable value. + glop::Fractional GetVariableValueAtCpScale(glop::ColIndex var); + + // TODO(user): use solver's precision epsilon. + constexpr static double kEpsilon = 1e-6; + + // Underlying LP solver API. + glop::LinearProgram lp_data_; + glop::RevisedSimplex simplex_; + + // For the scaling. + glop::SparseMatrixScaler scaler_; + + // violation_sum_ is used to simulate phase I of the simplex and be able to + // do reduced cost strengthening on problem feasibility by using the sum of + // constraint violations as an optimization objective. + glop::ColIndex violation_sum_; + ConstraintIndex violation_sum_constraint_; + + // Structures used for mirroring IntegerVariables inside the underlying LP + // solver: integer_variables_[i] is mirrored by mirror_lp_variables_[i], + // both are used in IncrementalPropagate() and Propagate() calls; + // integer_variable_to_index_ is used to find which mirroring variable's + // coefficient must be modified on SetCoefficient(). + hash_map integer_variable_to_index_; + std::vector integer_variables_; + std::vector mirror_lp_variables_; + + // We need to remember what to optimize if an objective is given, because + // then we will switch the objective between feasibility and optimization. + bool objective_is_defined_ = false; + IntegerVariable objective_cp_; + glop::ColIndex objective_lp_; + bool objective_is_minimization_; + + // Structures for propagators. + IntegerTrail* integer_trail_; + std::vector integer_reason_; + std::vector deductions_; + + // Last solution found by a call to the underlying LP solver. + // On IncrementalPropagate(), if the bound updates do not invalidate this + // solution, Propagate() will not find domain reductions, no need to call it. + std::vector lp_solution_; + + // Linear constraints cannot be created or modified after this is registered. + bool lp_constraint_is_registered_ = false; +}; + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_LINEAR_PROGRAMMING_CONSTRAINT_H_ diff --git a/src/sat/lp_utils.cc b/src/sat/lp_utils.cc index d3129fdc17..9ef006b18e 100644 --- a/src/sat/lp_utils.cc +++ b/src/sat/lp_utils.cc @@ -210,7 +210,7 @@ void ConvertBooleanProblemToLinearProgram(const LinearBooleanProblem& problem, lp->Clear(); for (int i = 0; i < problem.num_variables(); ++i) { const ColIndex col = lp->CreateNewVariable(); - lp->SetVariableIntegrality(col, true); + lp->SetVariableType(col, glop::LinearProgram::VariableType::INTEGER); lp->SetVariableBounds(col, 0.0, 1.0); } diff --git a/src/sat/lp_utils.h b/src/sat/lp_utils.h index fab702b575..f1fc086d1c 100644 --- a/src/sat/lp_utils.h +++ b/src/sat/lp_utils.h @@ -16,9 +16,9 @@ #ifndef OR_TOOLS_SAT_LP_UTILS_H_ #define OR_TOOLS_SAT_LP_UTILS_H_ -#include "sat/boolean_problem.pb.h" #include "linear_solver/linear_solver.pb.h" #include "lp_data/lp_data.h" +#include "sat/boolean_problem.pb.h" #include "sat/sat_solver.h" namespace operations_research { diff --git a/src/sat/optimization.cc b/src/sat/optimization.cc index 08071b7a3e..ff43558ca3 100644 --- a/src/sat/optimization.cc +++ b/src/sat/optimization.cc @@ -1045,7 +1045,7 @@ void LogSolveInfo(SatSolver::Status result, const SatSolver& sat_solver, SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( bool log_info, IntegerVariable objective_var, - const std::vector& var_for_lazy_encoding, + const std::function& next_decision, const std::function& feasible_solution_observer, Model* model) { // Timing. @@ -1066,7 +1066,7 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( IntegerValue objective(kint64max); while (true) { result = SolveIntegerProblemWithLazyEncoding(/*assumptions=*/{}, - var_for_lazy_encoding, model); + next_decision, model); if (result != SatSolver::MODEL_SAT) break; // The objective is the current lower bound of the objective_var. @@ -1074,7 +1074,9 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( // We have a solution! model_is_feasible = true; - feasible_solution_observer(*model); + if (feasible_solution_observer != nullptr) { + feasible_solution_observer(*model); + } // Restrict the objective. sat_solver->Backtrack(0); @@ -1107,7 +1109,7 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( bool log_info, const std::vector& literals, const std::vector& int64_coeffs, - const std::vector& var_for_lazy_encoding, + const std::function& next_decision, const std::function& feasible_solution_observer, Model* model) { // Timing. @@ -1184,8 +1186,8 @@ SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( } // Solve under the assumptions. - result = SolveIntegerProblemWithLazyEncoding(assumptions, - var_for_lazy_encoding, model); + result = + SolveIntegerProblemWithLazyEncoding(assumptions, next_decision, model); if (result == SatSolver::MODEL_SAT) { // Extract the new solution and save it if it is the best found so far. Coefficient obj(0); @@ -1195,7 +1197,9 @@ SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( } } if (obj + offset < upper_bound) { - feasible_solution_observer(*model); + if (feasible_solution_observer != nullptr) { + feasible_solution_observer(*model); + } upper_bound = obj + offset; if (log_info) LOG(INFO) << "c ub:" << upper_bound; } diff --git a/src/sat/optimization.h b/src/sat/optimization.h index 6af1f279b3..15587d67a8 100644 --- a/src/sat/optimization.h +++ b/src/sat/optimization.h @@ -121,26 +121,13 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScan( const std::function& feasible_solution_observer, Model* model); -// Same as MinimizeIntegerVariableWithLinearScan() but as long as the domain of -// the variables in var_for_lazy_encoding is not a singleton when the problem is -// solved to SAT (i.e all the Boolean variables are assigned), we add a new -// literal that can constrain the non-singleton variable with the lowest lower -// bound to its lower bound. We exploit the default polarity of the solver -// to make sure the first decision when resolving from the current state will be -// to force this variable to its lower bound. -// -// Note(user): For now, we pass the set of var_for_lazy_encoding which must -// contain a set of "decision" variables from which everything else can be -// deduced. On simple model like the jobshop, passing all the variables there -// have the same effect (but we don't want to pass their negation, or we need -// to change the heuristic to not fix variable to their upper bound). -// -// TODO(user): If it is too much work to provide this set, we could add a more -// complex heuristic to decide the next IntegerLiteral that will be associated -// to a decision in the SatSolver search. +// Same as MinimizeIntegerVariableWithLinearScan() but keep solving the problem +// as long as next_decision() do not return kNoLiteralIndex and hence lazily +// encode new variables. See the doc of SolveIntegerProblemWithLazyEncoding() +// for more details. SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( bool log_info, IntegerVariable objective_var, - const std::vector& var_for_lazy_encoding, + const std::function& next_decision, const std::function& feasible_solution_observer, Model* model); @@ -153,7 +140,7 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding( SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding( bool log_info, const std::vector& literals, const std::vector& coeffs, - const std::vector& var_for_lazy_encoding, + const std::function& next_decision, const std::function& feasible_solution_observer, Model* model); diff --git a/src/sat/precedences.cc b/src/sat/precedences.cc index 2aa6ae9c92..3335ee618d 100644 --- a/src/sat/precedences.cc +++ b/src/sat/precedences.cc @@ -316,6 +316,8 @@ void PrecedencesPropagator::AddArc(IntegerVariable tail, IntegerVariable head, bool PrecedencesPropagator::ArcShouldPropagate(const ArcInfo& arc, const Trail& trail) const { + if (integer_trail_->IsCurrentlyIgnored(arc.head_var)) return false; + const LiteralIndex index = OptionalLiteralOf(arc.tail_var); if (index == kNoLiteralIndex) return true; if (trail.Assignment().LiteralIsFalse(Literal(index))) return false; diff --git a/src/sat/precedences.h b/src/sat/precedences.h index 18a7fd360d..6d2b57227f 100644 --- a/src/sat/precedences.h +++ b/src/sat/precedences.h @@ -420,6 +420,13 @@ inline std::function ConditionalLowerOrEqualWithOffset( }; } +// is_le => (a <= b). +inline std::function ConditionalLowerOrEqual(IntegerVariable a, + IntegerVariable b, + Literal is_le) { + return ConditionalLowerOrEqualWithOffset(a, b, 0, is_le); +} + // is_le <=> (a + offset <= b). inline std::function ReifiedLowerOrEqualWithOffset( IntegerVariable a, IntegerVariable b, int64 offset, Literal is_le) { @@ -452,6 +459,25 @@ inline std::function ReifiedEquality(IntegerVariable a, }; } +// is_eq <=> (a + offset == b). +inline std::function ReifiedEqualityWithOffset(IntegerVariable a, + IntegerVariable b, + int64 offset, + Literal is_eq) { + return [=](Model* model) { + // We creates two extra Boolean variables in this case. + // + // TODO(user): Avoid creating them if we already have some literal that + // have the same meaning. For instance if a client also wanted to know if + // a <= b, he would have called ReifiedLowerOrEqualWithOffset() directly. + const Literal is_le = Literal(model->Add(NewBooleanVariable()), true); + const Literal is_ge = Literal(model->Add(NewBooleanVariable()), true); + model->Add(ReifiedBoolAnd({is_le, is_ge}, is_eq)); + model->Add(ReifiedLowerOrEqualWithOffset(a, b, offset, is_le)); + model->Add(ReifiedLowerOrEqualWithOffset(b, a, -offset, is_ge)); + }; +} + // a != b. inline std::function NotEqual(IntegerVariable a, IntegerVariable b) { diff --git a/src/sat/sat_parameters.proto b/src/sat/sat_parameters.proto index 770afe42b7..50483ff3f3 100644 --- a/src/sat/sat_parameters.proto +++ b/src/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: 80 +// NEXT TAG: 83 message SatParameters { // ========================================================================== // Branching and polarity @@ -452,4 +452,31 @@ message SatParameters { // 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]; + + // When this is true, the cumulative constraint is reinforced with propagators + // from the disjunctive constraint to improve the inference on a set of tasks + // that are disjunctive at the root of the problem. This additional level + // supplements the default level of reasoning. + // + // Propagators of the cumulative constraint will not be used at all if all the + // tasks are disjunctive at root node. + // + // 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 = false]; + + // When this is true, then we extract all the linear constraints of the + // CpModel and linearly encode some other constraints in an extra "linear + // programming" constraint that will run a LP solver to detect infeasibility + // and propagates the bounds on the objective. + optional bool use_global_lp_constraint = 81 [default = true]; + + // If true then all decisions taken by the solver are made using a fixed order + // as specified in the API or in the cp_model.proto search_order field. + // + // TODO(user): This is not working in all situation yet. The time limit is not + // really enforced properly, we don't support assumptions, and all the + // decisions variables must be integers. + optional bool use_fixed_search = 82 [default = false]; } diff --git a/src/sat/sat_solver.cc b/src/sat/sat_solver.cc index 69c4321801..f660198fef 100644 --- a/src/sat/sat_solver.cc +++ b/src/sat/sat_solver.cc @@ -2044,7 +2044,6 @@ void SatSolver::ComputeFirstUIPConflict( clause_to_expand = ClauseRef(); } else { clause_to_expand = trail_->Reason(literal.Variable()); - DCHECK(!clause_to_expand.IsEmpty()); } sat_clause = ReasonClauseOrNull(literal.Variable()); diff --git a/src/sat/sat_solver.h b/src/sat/sat_solver.h index eaef993e87..33c58577c9 100644 --- a/src/sat/sat_solver.h +++ b/src/sat/sat_solver.h @@ -501,6 +501,9 @@ class SatSolver { void EnqueueNewDecision(Literal literal); // Returns true if everything has been propagated. + // + // TODO(user): This test is fast but not exhaustive, especially regarding the + // integer propagators. Fix. bool PropagationIsDone() const; // Update the propagators_ list with the relevant propagators. diff --git a/src/sat/table.cc b/src/sat/table.cc index 7372dae592..fd7b25300e 100644 --- a/src/sat/table.cc +++ b/src/sat/table.cc @@ -185,6 +185,52 @@ std::function TableConstraint( }; } +std::function LiteralTableConstraint( + const std::vector>& literal_tuples, + const std::vector& line_literals) { + return [=](Model* model) { + CHECK_EQ(literal_tuples.size(), line_literals.size()); + const int num_tuples = line_literals.size(); + if (num_tuples == 0) return; + const int tuple_size = literal_tuples[0].size(); + if (tuple_size == 0) return; + for (int i = 1; i < num_tuples; ++i) { + CHECK_EQ(tuple_size, literal_tuples[i].size()); + } + + hash_map> line_literals_per_literal; + for (int i = 0; i < num_tuples; ++i) { + const LiteralIndex selected_index = line_literals[i].Index(); + for (const Literal l : literal_tuples[i]) { + line_literals_per_literal[l.Index()].push_back(selected_index); + } + } + + // line_literals[i] == true => literal_tuples[i][j] == true. + // literal_tuples[i][j] == false => line_literals[i] == false. + for (int i = 0; i < num_tuples; ++i) { + const Literal line_is_selected = line_literals[i]; + for (const Literal lit : literal_tuples[i]) { + model->Add(Implication(line_is_selected, lit)); + } + } + + // Exactly one selected literal is true. + model->Add(ExactlyOneConstraint(line_literals)); + + // If all selected literals of the lines containing a literal are false, + // then the literal is false. + for (const auto& p : line_literals_per_literal) { + std::vector clause; + for (const auto& index : p.second) { + clause.push_back(Literal(index)); + } + clause.push_back(Literal(p.first).Negated()); + model->Add(ClauseConstraint(clause)); + } + }; +} + std::function TransitionConstraint( const std::vector& vars, const std::vector>& automata, int64 initial_state, diff --git a/src/sat/table.h b/src/sat/table.h index 1bffc3f672..77df507c43 100644 --- a/src/sat/table.h +++ b/src/sat/table.h @@ -27,6 +27,14 @@ std::function TableConstraint( const std::vector& vars, const std::vector>& tuples); +// Enforces that exactly one literal in line_literals is true, and that +// all literals in the corresponding line of the literal_tuples matrix are true. +// This constraint assumes that exactly one literal per column of the +// literal_tuples matrix is true. +std::function LiteralTableConstraint( + const std::vector>& literal_tuples, + const std::vector& line_literals); + // Given an automata defined by a set of 3-tuples: // (state, transition_with_value_as_label, next_state) // this accepts the sequences of vars.size() variables that are recognized by diff --git a/src/sat/timetable.cc b/src/sat/timetable.cc index 53e2b6dda9..9662a34859 100644 --- a/src/sat/timetable.cc +++ b/src/sat/timetable.cc @@ -11,12 +11,14 @@ // See the License for the specific language governing permissions and // limitations under the License. +#include "sat/timetable.h" + #include #include "sat/overload_checker.h" #include "sat/precedences.h" #include "sat/sat_solver.h" -#include "sat/timetable.h" +#include "util/sort.h" namespace operations_research { namespace sat { @@ -33,18 +35,13 @@ TimeTablingPerTask::TimeTablingPerTask( trail_(trail), integer_trail_(integer_trail), intervals_repository_(intervals_repository) { + // Cached domains. start_min_.resize(num_tasks_); start_max_.resize(num_tasks_); end_min_.resize(num_tasks_); end_max_.resize(num_tasks_); duration_min_.resize(num_tasks_); demand_min_.resize(num_tasks_); - scp_.reserve(num_tasks_); - ecp_.reserve(num_tasks_); - // 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. - profile_.reserve(2 * num_tasks_ + 4); // Collect the variables. start_vars_.resize(num_tasks_); end_vars_.resize(num_tasks_); @@ -54,7 +51,16 @@ TimeTablingPerTask::TimeTablingPerTask( start_vars_[t] = intervals_repository->StartVar(i); end_vars_[t] = intervals_repository->EndVar(i); duration_vars_[t] = intervals_repository->SizeVar(i); + by_start_max_.push_back(TaskTime(t, IntegerValue(0))); + by_end_min_.push_back(TaskTime(t, IntegerValue(0))); } + // 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. + profile_.reserve(2 * num_tasks_ + 4); + in_profile_.resize(num_tasks_); + scp_.reserve(num_tasks_); + ecp_.reserve(num_tasks_); } void TimeTablingPerTask::RegisterWith(GenericLiteralWatcher* watcher) { @@ -99,11 +105,7 @@ bool TimeTablingPerTask::Propagate() { while (profile_changed_) { profile_changed_ = false; - // Rebuild compulsory part events. - // ------------------------------- - scp_.clear(); - ecp_.clear(); - // TODO(user): exclude extreme tasks during the search. + // Cache the variable bounds. for (int t = 0; t < num_tasks_; ++t) { start_min_[t] = StartMin(t); start_max_[t] = StartMax(t); @@ -111,97 +113,15 @@ bool TimeTablingPerTask::Propagate() { end_max_[t] = EndMax(t); demand_min_[t] = DemandMin(t); duration_min_[t] = DurationMin(t); - if (start_max_[t] < end_min_[t] && IsPresent(t)) { - scp_.push_back(Event(start_max_[t], t)); - ecp_.push_back(Event(end_min_[t], t)); - } } - // No filtering is possible. - if (scp_.empty()) return true; - - // TODO(user): use insertion sort with a fall back to std::sort. - // Sort compulsory part events. - std::sort(scp_.begin(), scp_.end()); - std::sort(ecp_.begin(), ecp_.end()); - - // Build the profile. - // ------------------ - profile_.clear(); - - // Add a sentinel to simplify the algorithm. - profile_.push_back( - ProfileRectangle(kMinIntegerValue, kMinIntegerValue, IntegerValue(0))); - - // Start and height of the currently built profile rectange. - IntegerValue current_start = kMinIntegerValue; - IntegerValue current_height = IntegerValue(0); - - // Start and height of the highest profile rectangle. - IntegerValue max_height_start = kMinIntegerValue; - IntegerValue max_height = IntegerValue(0); - - // Next scp and ecp events to be processed. - int next_scp = 0; - int next_ecp = 0; - - while (next_ecp < ecp_.size()) { - const IntegerValue old_height = current_height; - - // Next time point. - IntegerValue t = ecp_[next_ecp].time; - if (next_scp < ecp_.size()) { - t = std::min(t, scp_[next_scp].time); - } - - // Process the starting compulsory parts. - while (next_scp < scp_.size() && scp_[next_scp].time == t) { - current_height += demand_min_[scp_[next_scp].task_id]; - next_scp++; - } - - // Process the ending compulsory parts. - while (next_ecp < ecp_.size() && ecp_[next_ecp].time == t) { - current_height -= demand_min_[ecp_[next_ecp].task_id]; - next_ecp++; - } - - // Insert a new profile rectangle if any. - if (current_height != old_height) { - profile_.push_back(ProfileRectangle(current_start, t, old_height)); - if (current_height > max_height) { - max_height = current_height; - max_height_start = t; - } - current_start = t; - } - } - // Build the last profile rectangle. - DCHECK_EQ(current_height, 0); - profile_.push_back( - ProfileRectangle(current_start, kMaxIntegerValue, IntegerValue(0))); - - // Add a sentinel to simplify the algorithm. - profile_.push_back( - ProfileRectangle(kMaxIntegerValue, kMaxIntegerValue, IntegerValue(0))); - - // Filter the capacity variable. - // ----------------------------- - if (max_height > CapacityMin()) { - reason_.clear(); - literal_reason_.clear(); - ExplainProfileHeight(max_height_start); - if (!integer_trail_->Enqueue( - IntegerLiteral::GreaterOrEqual(capacity_var_, max_height), - literal_reason_, reason_)) { - return false; - } - } + // This can fail if the timetable exceeds the resource capacity. + if (!BuildTimeTable()) return false; // Update the start and end variables. // ----------------------------------- // Tasks with a lower or equal demand will not be pushed. - const IntegerValue min_demand = CapacityMax() - max_height; + const IntegerValue min_demand = CapacityMax() - profile_max_height_; for (int t = 0; t < num_tasks_; ++t) { // The task cannot be pushed. @@ -222,6 +142,110 @@ bool TimeTablingPerTask::Propagate() { return true; } +bool TimeTablingPerTask::BuildTimeTable() { + // Update the value of the events to sort. + for (int i = 0; i < num_tasks_; ++i) { + by_start_max_[i].time = start_max_[by_start_max_[i].task_id]; + by_end_min_[i].time = end_min_[by_end_min_[i].task_id]; + } + // Likely to be already or almost sorted. + IncrementalSort(by_start_max_.begin(), by_start_max_.end()); + IncrementalSort(by_end_min_.begin(), by_end_min_.end()); + + scp_.clear(); + ecp_.clear(); + + // Build start of compulsory part events. + for (int i = 0; i < num_tasks_; ++i) { + const int t = by_start_max_[i].task_id; + in_profile_[t] = IsPresent(t) && start_max_[t] < end_min_[t]; + if (in_profile_[t]) scp_.push_back(by_start_max_[i]); + } + + // Build end of compulsory part events. + for (int i = 0; i < num_tasks_; ++i) { + const int t = by_end_min_[i].task_id; + if (in_profile_[t]) ecp_.push_back(by_end_min_[i]); + } + + DCHECK_EQ(scp_.size(), ecp_.size()); + + // Build the profile. + // ------------------ + profile_.clear(); + + // Start and height of the highest profile rectangle. + profile_max_height_ = kMinIntegerValue; + IntegerValue max_height_start = kMinIntegerValue; + + // Add a sentinel to simplify the algorithm. + profile_.push_back( + ProfileRectangle(kMinIntegerValue, kMinIntegerValue, IntegerValue(0))); + + // Start and height of the currently built profile rectange. + IntegerValue current_start = kMinIntegerValue; + IntegerValue current_height = IntegerValue(0); + + // Next scp and ecp events to be processed. + int next_scp = 0; + int next_ecp = 0; + + while (next_ecp < ecp_.size()) { + const IntegerValue old_height = current_height; + + // Next time point. + IntegerValue t = ecp_[next_ecp].time; + if (next_scp < ecp_.size()) { + t = std::min(t, scp_[next_scp].time); + } + + // Process the starting compulsory parts. + while (next_scp < scp_.size() && scp_[next_scp].time == t) { + current_height += demand_min_[scp_[next_scp].task_id]; + next_scp++; + } + + // Process the ending compulsory parts. + while (next_ecp < ecp_.size() && ecp_[next_ecp].time == t) { + current_height -= demand_min_[ecp_[next_ecp].task_id]; + next_ecp++; + } + + // Insert a new profile rectangle if any. + if (current_height != old_height) { + profile_.push_back(ProfileRectangle(current_start, t, old_height)); + if (current_height > profile_max_height_) { + profile_max_height_ = current_height; + max_height_start = t; + } + current_start = t; + } + } + // Build the last profile rectangle. + DCHECK_EQ(current_height, 0); + profile_.push_back( + ProfileRectangle(current_start, kMaxIntegerValue, IntegerValue(0))); + + // Add a sentinel to simplify the algorithm. + profile_.push_back( + ProfileRectangle(kMaxIntegerValue, kMaxIntegerValue, IntegerValue(0))); + + // Filter the capacity variable. + // ----------------------------- + if (profile_max_height_ > CapacityMin()) { + reason_.clear(); + literal_reason_.clear(); + ExplainProfileHeight(max_height_start); + if (!integer_trail_->Enqueue( + IntegerLiteral::GreaterOrEqual(capacity_var_, profile_max_height_), + literal_reason_, reason_)) { + return false; + } + } + + return true; +} + bool TimeTablingPerTask::SweepTaskRight(int task_id) { // Find the profile rectangle that overlaps the start min of the task. // The sentinel prevents out of bound exceptions. @@ -246,8 +270,8 @@ bool TimeTablingPerTask::SweepTaskRight(int task_id) { } // If the task cannot be scheduled after the conflicting profile rectangle, // we explain all the intermediate pushes to schedule the task to its - // start max. Scheduling the task to its start max may result in a capacity - // overload that will be detected once the profile is rebuilt. + // start max. If the task is not in the profile, then remove the task if it + // is optional or explain the overload. if (s_max < profile_[rec_id].end) { while (end_min_[task_id] < s_max) { if (!UpdateStartingTime(task_id, end_min_[task_id])) return false; @@ -255,6 +279,8 @@ bool TimeTablingPerTask::SweepTaskRight(int task_id) { if (start_min_[task_id] < s_max) { if (!UpdateStartingTime(task_id, s_max)) return false; } + // This function rely on the updated start_min above to build its reason. + if (!in_profile_[task_id]) return OverloadOrRemove(task_id, s_max); profile_changed_ = true; return true; } @@ -297,8 +323,8 @@ bool TimeTablingPerTask::SweepTaskLeft(int task_id) { } // If the task cannot be scheduled before the conflicting profile rectangle, // we explain all the intermediate pushes to schedule the task to its - // end min. Scheduling the task to its end min may result in a capacity - // overload that will be detected once the profile is rebuilt. + // end min. If the task is not in the profile, then remove the task if it is + // optional or explain the overload. if (profile_[rec_id].start < e_min) { while (e_min < start_max_[task_id]) { if (!UpdateEndingTime(task_id, start_max_[task_id])) return false; @@ -306,6 +332,8 @@ bool TimeTablingPerTask::SweepTaskLeft(int task_id) { if (e_min < end_max_[task_id]) { if (!UpdateEndingTime(task_id, e_min)) return false; } + // This function rely on the updated end_max above to build its reason. + if (!in_profile_[task_id]) return OverloadOrRemove(task_id, e_min - 1); profile_changed_ = true; return true; } @@ -428,16 +456,46 @@ void TimeTablingPerTask::AddPresenceReasonIfNeeded(int task_id) { } } -void TimeTablingPerTask::ExplainProfileHeight(IntegerValue time) { +IntegerValue TimeTablingPerTask::ExplainProfileHeight(IntegerValue time) { + IntegerValue height = IntegerValue(0); for (int t = 0; t < num_tasks_; ++t) { // Tasks need to overlap the time point, i.e., start_max <= time < end_min. if (start_max_[t] <= time && time < end_min_[t] && IsPresent(t)) { + height += demand_min_[t]; reason_.push_back(integer_trail_->LowerBoundAsLiteral(demand_vars_[t])); reason_.push_back(IntegerLiteral::LowerOrEqual(start_vars_[t], time)); reason_.push_back(IntegerLiteral::GreaterOrEqual(end_vars_[t], time + 1)); AddPresenceReasonIfNeeded(t); } } + return height; +} + +bool TimeTablingPerTask::OverloadOrRemove(int task_id, IntegerValue time) { + reason_.clear(); + literal_reason_.clear(); + reason_.push_back(integer_trail_->UpperBoundAsLiteral(capacity_var_)); + + // Explain the resource overload if the task cannot be removed. + const IntegerValue height = ExplainProfileHeight(time); + if (IsPresent(task_id)) { + return integer_trail_->Enqueue( + IntegerLiteral::GreaterOrEqual(capacity_var_, height), literal_reason_, + reason_); + } + + // Remove the task to prevent the overload. + reason_.push_back(integer_trail_->LowerBoundAsLiteral(demand_vars_[task_id])); + reason_.push_back(IntegerLiteral::LowerOrEqual(start_vars_[task_id], time)); + reason_.push_back( + IntegerLiteral::GreaterOrEqual(end_vars_[task_id], time + 1)); + + integer_trail_->EnqueueLiteral( + intervals_repository_->IsPresentLiteral(interval_vars_[task_id]) + .Negated(), + literal_reason_, reason_); + + return true; } } // namespace sat diff --git a/src/sat/timetable.h b/src/sat/timetable.h index c199a4955c..b1e1ba6d3a 100644 --- a/src/sat/timetable.h +++ b/src/sat/timetable.h @@ -37,11 +37,11 @@ class TimeTablingPerTask : public PropagatorInterface { void RegisterWith(GenericLiteralWatcher* watcher); private: - struct Event { - /* const */ IntegerValue time; + struct TaskTime { /* const */ int task_id; - Event(IntegerValue time, int task_id) : time(time), task_id(task_id) {} - bool operator<(Event other) const { return time < other.time; } + IntegerValue time; + TaskTime(int task_id, IntegerValue time) : task_id(task_id), time(time) {} + bool operator<(TaskTime other) const { return time < other.time; } }; struct ProfileRectangle { @@ -52,11 +52,15 @@ class TimeTablingPerTask : public PropagatorInterface { : start(start), end(end), height(height) {} }; - // Increase the start min of task task_id. This function may call + // Builds the time table and increases the lower bound of the capacity + // variable accordingly. + bool BuildTimeTable(); + + // Increases the start min of task task_id. This function may call // UpdateStartingTime(). bool SweepTaskRight(int task_id); - // Decrease the end max of task task_id. This function may call + // Decreases the end max of task task_id. This function may call // UpdateEndingTime(). bool SweepTaskLeft(int task_id); @@ -68,9 +72,13 @@ class TimeTablingPerTask : public PropagatorInterface { // reason_ with the corresponding reason. bool UpdateEndingTime(int task_id, IntegerValue new_end); + // Explains the resource overload at time or removes task_id if it is + // optional. + bool OverloadOrRemove(int task_id, IntegerValue time); + // Fills reason_ with the reason that explains the height of the profile at - // the given time point. - void ExplainProfileHeight(IntegerValue time); + // the given time point. Also return the height of the profile at time. + IntegerValue ExplainProfileHeight(IntegerValue time); IntegerValue StartMin(int task_id) const { return integer_trail_->LowerBound(start_vars_[task_id]); @@ -139,13 +147,20 @@ class TimeTablingPerTask : public PropagatorInterface { std::vector duration_min_; std::vector demand_min_; - // Events that represent the start of a compulsory part. - std::vector scp_; - // Events that represent the end of a compulsory part. - std::vector ecp_; + // Tasks sorted by start max (resp. end min). + std::vector by_start_max_; + std::vector by_end_min_; + + // Start (resp. end) of the compulsory parts used to build the profile. + std::vector scp_; + std::vector ecp_; // Optimistic profile of the resource consumption over time. std::vector profile_; + IntegerValue profile_max_height_; + + // True if the corresponding task is part of the profile. + std::vector in_profile_; // True if the last call of the propagator has filtered the domain of a task // and changed the shape of the profile. diff --git a/src/sat/timetable_edgefinding.cc b/src/sat/timetable_edgefinding.cc index 6c74f6db0c..cfd8881233 100644 --- a/src/sat/timetable_edgefinding.cc +++ b/src/sat/timetable_edgefinding.cc @@ -132,21 +132,29 @@ void TimeTableEdgeFinding::SwitchToMirrorProblem() { } void TimeTableEdgeFinding::BuildTimeTable() { - std::vector scp; - std::vector ecp; + scp_.clear(); + ecp_.clear(); - // Build the sorted mandatory parts. + // Build start of compulsory part events. for (int i = 0; i < num_tasks_; ++i) { - const IntegerValue s_max = by_start_max_[i].time; - if (s_max < end_min_[by_start_max_[i].task_id]) { - scp.push_back(by_start_max_[i]); - } - const IntegerValue e_min = by_end_min_[i].time; - if (start_max_[by_end_min_[i].task_id] < e_min) { - ecp.push_back(by_end_min_[i]); + const int t = by_start_max_[i].task_id; + if (!IsPresent(t)) continue; + if (start_max_[t] < end_min_[t]) { + scp_.push_back(by_start_max_[i]); } } + // Build end of compulsory part events. + for (int i = 0; i < num_tasks_; ++i) { + const int t = by_end_min_[i].task_id; + if (!IsPresent(t)) continue; + if (start_max_[t] < end_min_[t]) { + ecp_.push_back(by_end_min_[i]); + } + } + + DCHECK_EQ(scp_.size(), ecp_.size()); + IntegerValue height = IntegerValue(0); IntegerValue energy = IntegerValue(0); IntegerValue previous_time = by_start_min_[0].time; @@ -163,11 +171,11 @@ void TimeTableEdgeFinding::BuildTimeTable() { if (index_smin < num_tasks_) { time = std::min(time, by_start_min_[index_smin].time); } - if (index_scp < scp.size()) { - time = std::min(time, scp[index_scp].time); + if (index_scp < scp_.size()) { + time = std::min(time, scp_[index_scp].time); } - if (index_ecp < ecp.size()) { - time = std::min(time, ecp[index_ecp].time); + if (index_ecp < ecp_.size()) { + time = std::min(time, ecp_[index_ecp].time); } // Total amount of energy contained in the timetable until time. @@ -188,14 +196,14 @@ void TimeTableEdgeFinding::BuildTimeTable() { } // Process the starting compulsory parts. - while (index_scp < scp.size() && scp[index_scp].time == time) { - height += demand_min_[scp[index_scp].task_id]; + while (index_scp < scp_.size() && scp_[index_scp].time == time) { + height += demand_min_[scp_[index_scp].task_id]; index_scp++; } // Process the ending compulsory parts. - while (index_ecp < ecp.size() && ecp[index_ecp].time == time) { - height -= demand_min_[ecp[index_ecp].task_id]; + while (index_ecp < ecp_.size() && ecp_[index_ecp].time == time) { + height -= demand_min_[ecp_[index_ecp].task_id]; index_ecp++; } diff --git a/src/sat/timetable_edgefinding.h b/src/sat/timetable_edgefinding.h index 2eb10d6caf..01161548ef 100644 --- a/src/sat/timetable_edgefinding.h +++ b/src/sat/timetable_edgefinding.h @@ -176,6 +176,10 @@ class TimeTableEdgeFinding : public PropagatorInterface { std::vector by_end_min_; std::vector by_end_max_; + // Start (resp. end) of the compulsory parts used to build the profile. + std::vector scp_; + std::vector ecp_; + // Energy of the free parts. std::vector energy_free_;