improvements in sat: split sub-propagators in scheduling constraints; add linear_programming_constraint in sat using glop

This commit is contained in:
Laurent Perron
2017-03-28 16:11:06 +02:00
parent 0b93b70d41
commit 3a6a5550d5
27 changed files with 2137 additions and 947 deletions

View File

@@ -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) {

View File

@@ -11,9 +11,10 @@
// See the License for the specific language governing permissions and
// limitations under the License.
#include "sat/cumulative.h"
#include <algorithm>
#include "sat/cumulative.h"
#include "sat/disjunctive.h"
#include "sat/overload_checker.h"
#include "sat/sat_solver.h"
@@ -29,64 +30,85 @@ std::function<void(Model*)> Cumulative(
const IntegerVariable& capacity) {
return [=](Model* model) {
if (vars.empty()) return;
IntervalsRepository* intervals = model->GetOrCreate<IntervalsRepository>();
IntegerEncoder* encoder = model->GetOrCreate<IntegerEncoder>();
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<IntervalVariable> 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<SatSolver>()->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<IntervalVariable> 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<Trail>();
@@ -100,7 +122,6 @@ std::function<void(Model*)> Cumulative(
time_tabling->RegisterWith(model->GetOrCreate<GenericLiteralWatcher>());
model->TakeOwnership(time_tabling);
const SatParameters& parameters = model->Get<SatSolver>()->parameters();
// Propagator responsible for applying the Overload Checking filtering rule.
// It increases the minimum of the capacity variable.
@@ -132,6 +153,8 @@ std::function<void(Model*)> 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<IntegerTrail>();
SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
@@ -143,14 +166,9 @@ std::function<void(Model*)> CumulativeTimeDecomposition(
std::vector<IntegerValue> 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<void(Model*)> 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<LiteralWithCoeff> 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<Literal> 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,

View File

@@ -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<void(Model*)> 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<void(Model*)> CumulativeTimeDecomposition(
const std::vector<IntervalVariable>& vars,
const std::vector<IntegerVariable>& demand_vars,

File diff suppressed because it is too large Load Diff

View File

@@ -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<void(Model*)> Disjunctive(
const std::vector<IntervalVariable>& 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<void(Model*)> DisjunctiveWithBooleanPrecedencesOnly(
const std::vector<IntervalVariable>& vars);
// Same as Disjunctive() + DisjunctiveWithBooleanPrecedencesOnly().
std::function<void(Model*)> DisjunctiveWithBooleanPrecedences(
const std::vector<IntervalVariable>& 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<Entry>& 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<IntervalVariable>& 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<bool> 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<IntegerVariable> start_vars_;
std::vector<IntegerVariable> end_vars_;
std::vector<IntegerVariable> duration_vars_;
std::vector<IntegerValue> fixed_durations_;
std::vector<LiteralIndex> 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<IntegerVariable> minus_start_vars_;
std::vector<IntegerVariable> minus_end_vars_;
// This is used by PrecedencePass().
std::vector<LiteralIndex> reason_for_beeing_before_;
std::vector<PrecedencesPropagator::IntegerPrecedences> before_;
std::vector<PrecedencesPropagator::IntegerPrecedences> 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<bool> 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<int> task_by_increasing_min_start_;
std::vector<int> task_by_increasing_min_end_;
std::vector<int> task_by_decreasing_max_start_;
std::vector<int> task_by_decreasing_max_end_;
// Reason vectors. Declared here to avoid costly initialization.
std::vector<Literal> literal_reason_;
std::vector<IntegerLiteral> integer_reason_;
const std::vector<IntervalVariable> non_overlapping_intervals_;
Trail* trail_;
private:
const bool time_direction_;
SchedulingConstraintHelper* helper_;
IntegerTrail* integer_trail_;
IntervalsRepository* intervals_;
PrecedencesPropagator* precedences_;
TaskSet task_set_;
std::vector<bool> is_gray_;
SatParameters parameters_;
mutable StatsGroup stats_;
DISALLOW_COPY_AND_ASSIGN(DisjunctiveConstraint);
std::vector<LiteralIndex> reason_for_beeing_before_;
std::vector<PrecedencesPropagator::IntegerPrecedences> before_;
};
} // namespace sat

View File

@@ -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<int>(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<int>(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>& literal_reason,
const std::vector<IntegerLiteral>& 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<Literal>* 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<int32>(
literals_reason_buffer_.size()),
/*dependencies_start_index=*/static_cast<int32>(
bounds_reason_buffer_.size())});
/*literals_reason_start_index=*/
static_cast<int32>(literals_reason_buffer_.size()),
/*dependencies_start_index=*/
static_cast<int32>(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<Literal> 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<int32>(
literals_reason_buffer_.size()),
/*dependencies_start_index=*/static_cast<int32>(
bounds_reason_buffer_.size())});
/*literals_reason_start_index=*/
static_cast<int32>(literals_reason_buffer_.size()),
/*dependencies_start_index=*/
static_cast<int32>(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<Literal>& assumptions,
const std::vector<IntegerVariable>& variables_to_finalize, Model* model) {
TimeLimit* time_limit = model->Mutable<TimeLimit>();
// 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<LiteralIndex()> FirstUnassignedVarAtItsMinHeuristic(
const std::vector<IntegerVariable>& vars, Model* model) {
IntegerEncoder* const integer_encoder = model->GetOrCreate<IntegerEncoder>();
IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
SatSolver* const solver = model->GetOrCreate<SatSolver>();
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<LiteralIndex()> UnassignedVarWithLowestMinAtItsMinHeuristic(
const std::vector<IntegerVariable>& vars, Model* model) {
IntegerEncoder* const integer_encoder = model->GetOrCreate<IntegerEncoder>();
IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
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<Literal>& assumptions,
const std::function<LiteralIndex()>& 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<TimeLimit>();
if (time_limit == nullptr) {
model->SetSingleton(TimeLimit::Infinite());
time_limit = model->Mutable<TimeLimit>();
}
SatSolver* const solver = model->GetOrCreate<SatSolver>();
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<IntegerTrail>()->NumIntegerVariables();
std::vector<IntegerVariable> 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<void(Model*)>
ExcludeCurrentSolutionWithoutIgnoredVariableAndBacktrack() {
return [=](Model* model) {
SatSolver* sat_solver = model->GetOrCreate<SatSolver>();
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
IntegerEncoder* encoder = model->GetOrCreate<IntegerEncoder>();
const int current_level = sat_solver->CurrentDecisionLevel();
std::vector<Literal> 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

View File

@@ -18,6 +18,7 @@
#include <set>
#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<Literal> literals_reason_buffer_;
mutable std::vector<IntegerLiteral> bounds_reason_buffer_;
// The "is_empty" literal of the optional variables or kNoLiteralIndex.
ITIVector<IntegerVariable, LiteralIndex> is_empty_literals_;
// The "is_ignored" literal of the optional variables or kNoLiteralIndex.
ITIVector<IntegerVariable, LiteralIndex> 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<int>& watch_indices) {
LOG(FATAL) << "Not implemented.";
return false;
}
};
@@ -882,6 +895,7 @@ inline std::function<IntegerVariable(Model*)> ConstantIntegerVariable(
inline std::function<IntegerVariable(Model*)> NewIntegerVariable(int64 lb,
int64 ub) {
return [=](Model* model) {
CHECK_LE(lb, ub);
return model->GetOrCreate<IntegerTrail>()->AddIntegerVariable(
IntegerValue(lb), IntegerValue(ub));
};
@@ -894,15 +908,44 @@ inline std::function<IntegerVariable(Model*)> 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<IntegerEncoder>()
->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<LiteralIndex, IntegerVariable> 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<IntegerVariable(Model*)> NewIntegerVariableFromLiteral(
Literal lit) {
return [=](Model* model) {
const IntegerVariable int_var = model->Add(NewIntegerVariable(0, 1));
model->GetOrCreate<IntegerEncoder>()->FullyEncodeVariableUsingGivenLiterals(
int_var, {lit.Negated(), lit}, {IntegerValue(0), IntegerValue(1)});
return int_var;
return model->GetOrCreate<LiteralViews>()->GetIntegerView(lit);
};
}
@@ -929,7 +972,7 @@ inline std::function<bool(const Model&)> IsFixed(IntegerVariable v) {
inline std::function<int64(const Model&)> Value(IntegerVariable v) {
return [=](const Model& model) {
const IntegerTrail* trail = model.Get<IntegerTrail>();
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<void(Model*)> 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<void(Model*)> Implication(Literal l, IntegerLiteral i) {
return [=](Model* model) {
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
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<IntegerEncoder>();
const Literal current = encoder->GetOrCreateAssociatedLiteral(i);
model->Add(Implication(l, current));
}
};
}
// in_interval <=> v in [lb, ub].
inline std::function<void(Model*)> 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<Literal>& assumptions,
const std::vector<IntegerVariable>& vars_with_lazy_encoding, Model* model);
const std::function<LiteralIndex()>& 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<LiteralIndex()> FirstUnassignedVarAtItsMinHeuristic(
const std::vector<IntegerVariable>& 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<LiteralIndex()> UnassignedVarWithLowestMinAtItsMinHeuristic(
const std::vector<IntegerVariable>& 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<void(Model*)>
ExcludeCurrentSolutionWithoutIgnoredVariableAndBacktrack();
} // namespace sat
} // namespace operations_research

View File

@@ -173,7 +173,7 @@ inline std::function<void(Model*)> WeightedSumLowerOrEqual(
const std::vector<IntegerVariable>& 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<void(Model*)> 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)) {

View File

@@ -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<IntervalVariable>& 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::TaskTime>&
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::TaskTime>&
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::TaskTime>&
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<TaskTime>());
return task_by_decreasing_max_start_;
}
const std::vector<SchedulingConstraintHelper::TaskTime>&
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<TaskTime>());
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

View File

@@ -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<IntervalVariable>& 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<TaskTime>& TaskByIncreasingStartMin();
const std::vector<TaskTime>& TaskByIncreasingEndMin();
const std::vector<TaskTime>& TaskByDecreasingStartMax();
const std::vector<TaskTime>& 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<Literal>* MutableLiteralReason() { return &literal_reason_; }
std::vector<IntegerLiteral>* 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<IntegerVariable>& StartVars() const { return start_vars_; }
const std::vector<IntegerVariable>& 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<IntegerVariable> start_vars_;
std::vector<IntegerVariable> end_vars_;
std::vector<IntegerVariable> duration_vars_;
std::vector<IntegerValue> fixed_durations_;
std::vector<LiteralIndex> 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<IntegerVariable> minus_start_vars_;
std::vector<IntegerVariable> minus_end_vars_;
// Sorted vectors returned by the TasksBy*() functions.
std::vector<TaskTime> task_by_increasing_min_start_;
std::vector<TaskTime> task_by_increasing_min_end_;
std::vector<TaskTime> task_by_decreasing_max_start_;
std::vector<TaskTime> task_by_decreasing_max_end_;
// Reason vectors.
std::vector<Literal> literal_reason_;
std::vector<IntegerLiteral> 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.
// =============================================================================

View File

@@ -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<double>::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<double>::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<double>::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<int>& watch_indices) {
for (const int index : watch_indices) {
const double lb = static_cast<double>(
integer_trail_->LowerBound(integer_variables_[index]).value());
const double ub = static_cast<double>(
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<double>(integer_trail_->LowerBound(cp_var).value());
const double ub =
static_cast<double>(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<double>::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<double>(integer_trail_->LowerBound(objective_cp_).value());
const double objective_cp_ub =
static_cast<double>(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<IntegerValue>(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<IntegerValue>(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<double>(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<IntegerValue>(new_ub));
deductions_.push_back(IntegerLiteral::LowerOrEqual(cp_var, new_ub_int));
}
} else if (rc < -kEpsilon) {
const double lb =
static_cast<double>(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<IntegerValue>(new_lb));
deductions_.push_back(
IntegerLiteral::GreaterOrEqual(cp_var, new_lb_int));
}
}
}
}
} // namespace sat
} // namespace operations_research

View File

@@ -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 <limits>
#include <vector>
#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<double>::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<IntegerTrail>();
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<int>& 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<IntegerVariable, int> integer_variable_to_index_;
std::vector<IntegerVariable> integer_variables_;
std::vector<glop::ColIndex> 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<IntegerLiteral> integer_reason_;
std::vector<IntegerLiteral> 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<double> 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_

View File

@@ -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);
}

View File

@@ -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 {

View File

@@ -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<IntegerVariable>& var_for_lazy_encoding,
const std::function<LiteralIndex()>& next_decision,
const std::function<void(const Model&)>& 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<Literal>& literals,
const std::vector<int64>& int64_coeffs,
const std::vector<IntegerVariable>& var_for_lazy_encoding,
const std::function<LiteralIndex()>& next_decision,
const std::function<void(const Model&)>& 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;
}

View File

@@ -121,26 +121,13 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScan(
const std::function<void(const Model&)>& 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<IntegerVariable>& var_for_lazy_encoding,
const std::function<LiteralIndex()>& next_decision,
const std::function<void(const Model&)>& feasible_solution_observer,
Model* model);
@@ -153,7 +140,7 @@ SatSolver::Status MinimizeIntegerVariableWithLinearScanAndLazyEncoding(
SatSolver::Status MinimizeWeightedLiteralSumWithCoreAndLazyEncoding(
bool log_info, const std::vector<Literal>& literals,
const std::vector<int64>& coeffs,
const std::vector<IntegerVariable>& var_for_lazy_encoding,
const std::function<LiteralIndex()>& next_decision,
const std::function<void(const Model&)>& feasible_solution_observer,
Model* model);

View File

@@ -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;

View File

@@ -420,6 +420,13 @@ inline std::function<void(Model*)> ConditionalLowerOrEqualWithOffset(
};
}
// is_le => (a <= b).
inline std::function<void(Model*)> ConditionalLowerOrEqual(IntegerVariable a,
IntegerVariable b,
Literal is_le) {
return ConditionalLowerOrEqualWithOffset(a, b, 0, is_le);
}
// is_le <=> (a + offset <= b).
inline std::function<void(Model*)> ReifiedLowerOrEqualWithOffset(
IntegerVariable a, IntegerVariable b, int64 offset, Literal is_le) {
@@ -452,6 +459,25 @@ inline std::function<void(Model*)> ReifiedEquality(IntegerVariable a,
};
}
// is_eq <=> (a + offset == b).
inline std::function<void(Model*)> 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<void(Model*)> NotEqual(IntegerVariable a,
IntegerVariable b) {

View File

@@ -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];
}

View File

@@ -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());

View File

@@ -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.

View File

@@ -185,6 +185,52 @@ std::function<void(Model*)> TableConstraint(
};
}
std::function<void(Model*)> LiteralTableConstraint(
const std::vector<std::vector<Literal>>& literal_tuples,
const std::vector<Literal>& 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<LiteralIndex, std::vector<LiteralIndex>> 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<Literal> 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<void(Model*)> TransitionConstraint(
const std::vector<IntegerVariable>& vars,
const std::vector<std::vector<int64>>& automata, int64 initial_state,

View File

@@ -27,6 +27,14 @@ std::function<void(Model*)> TableConstraint(
const std::vector<IntegerVariable>& vars,
const std::vector<std::vector<int64>>& 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<void(Model*)> LiteralTableConstraint(
const std::vector<std::vector<Literal>>& literal_tuples,
const std::vector<Literal>& 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

View File

@@ -11,12 +11,14 @@
// See the License for the specific language governing permissions and
// limitations under the License.
#include "sat/timetable.h"
#include <algorithm>
#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

View File

@@ -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<IntegerValue> duration_min_;
std::vector<IntegerValue> demand_min_;
// Events that represent the start of a compulsory part.
std::vector<Event> scp_;
// Events that represent the end of a compulsory part.
std::vector<Event> ecp_;
// Tasks sorted by start max (resp. end min).
std::vector<TaskTime> by_start_max_;
std::vector<TaskTime> by_end_min_;
// Start (resp. end) of the compulsory parts used to build the profile.
std::vector<TaskTime> scp_;
std::vector<TaskTime> ecp_;
// Optimistic profile of the resource consumption over time.
std::vector<ProfileRectangle> profile_;
IntegerValue profile_max_height_;
// True if the corresponding task is part of the profile.
std::vector<bool> in_profile_;
// True if the last call of the propagator has filtered the domain of a task
// and changed the shape of the profile.

View File

@@ -132,21 +132,29 @@ void TimeTableEdgeFinding::SwitchToMirrorProblem() {
}
void TimeTableEdgeFinding::BuildTimeTable() {
std::vector<TaskTime> scp;
std::vector<TaskTime> 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++;
}

View File

@@ -176,6 +176,10 @@ class TimeTableEdgeFinding : public PropagatorInterface {
std::vector<TaskTime> by_end_min_;
std::vector<TaskTime> by_end_max_;
// Start (resp. end) of the compulsory parts used to build the profile.
std::vector<TaskTime> scp_;
std::vector<TaskTime> ecp_;
// Energy of the free parts.
std::vector<IntegerValue> energy_free_;