always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

[CP-SAT] add basic scheduling lns; improve diffn constraint; rewrite alldiff; polish branching code

Browse Source
This commit is contained in:
Laurent Perron
2019-02-26 14:42:30 +01:00
parent 0e3968b3eb
commit cfbf1db09b
11 changed files with 603 additions and 423 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
ortools/lp_data/sparse_vector.h
View File

@@ -503,8 +503,8 @@ SparseVector<IndexType, IteratorType>::SparseVector(const SparseVector& other) {
}
template <typename IndexType, typename IteratorType>
SparseVector<IndexType, IteratorType>& SparseVector<IndexType, IteratorType>::
operator=(const SparseVector& other) {
SparseVector<IndexType, IteratorType>&
SparseVector<IndexType, IteratorType>::operator=(const SparseVector& other) {
PopulateFromSparseVector(other);
return *this;
}

223
ortools/sat/all_different.cc
View File

@@ -56,7 +56,8 @@ std::function<void(Model*)> AllDifferentBinary(
std::function<void(Model*)> AllDifferentOnBounds(
const std::vector<IntegerVariable>& vars) {
return [=](Model* model) {
AllDifferentBoundsPropagator* constraint = new AllDifferentBoundsPropagator(
if (vars.empty()) return;
auto* constraint = new AllDifferentBoundsPropagator(
vars, model->GetOrCreate<IntegerTrail>());
constraint->RegisterWith(model->GetOrCreate<GenericLiteralWatcher>());
model->TakeOwnership(constraint);
@@ -407,14 +408,22 @@ bool AllDifferentConstraint::Propagate() {
AllDifferentBoundsPropagator::AllDifferentBoundsPropagator(
const std::vector<IntegerVariable>& vars, IntegerTrail* integer_trail)
: vars_(vars), integer_trail_(integer_trail), num_calls_(0) {
: integer_trail_(integer_trail) {
CHECK(!vars.empty());
// We need +2 for sentinels.
const int capacity = vars.size() + 2;
index_to_start_index_.resize(capacity);
index_to_end_index_.resize(capacity);
index_to_var_.resize(capacity, kNoIntegerVariable);
for (int i = 0; i < vars.size(); ++i) {
negated_vars_.push_back(NegationOf(vars_[i]));
vars_.push_back({vars[i]});
negated_vars_.push_back({NegationOf(vars[i])});
}
}
bool AllDifferentBoundsPropagator::Propagate() {
if (vars_.empty()) return true;
if (!PropagateLowerBounds()) return false;
// Note that it is not required to swap back vars_ and negated_vars_.
@@ -425,77 +434,183 @@ bool AllDifferentBoundsPropagator::Propagate() {
return result;
}
// TODO(user): we could gain by pushing all the new bound at the end, so that
// we just have to sort to_insert_ once.
void AllDifferentBoundsPropagator::FillHallReason(IntegerValue hall_lb,
IntegerValue hall_ub) {
for (auto entry : to_insert_) {
value_to_variable_[entry.first] = entry.second;
}
to_insert_.clear();
integer_reason_.clear();
for (int64 v = hall_lb.value(); v <= hall_ub; ++v) {
const IntegerVariable var = gtl::FindOrDie(value_to_variable_, v);
const int limit = GetIndex(hall_ub);
for (int i = GetIndex(hall_lb); i <= limit; ++i) {
const IntegerVariable var = index_to_var_[i];
integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, hall_lb));
integer_reason_.push_back(IntegerLiteral::LowerOrEqual(var, hall_ub));
}
}
bool AllDifferentBoundsPropagator::PropagateLowerBounds() {
++num_calls_;
critical_intervals_.clear();
hall_starts_.clear();
hall_ends_.clear();
to_insert_.clear();
if (num_calls_ % 20 == 0) {
// We don't really need to clear this, but we do from time to time to
// save memory (in case the variable domains are huge). This optimization
// helps a bit.
value_to_variable_.clear();
int AllDifferentBoundsPropagator::FindStartIndexAndCompressPath(int index) {
// First, walk the pointer until we find one pointing to itself.
int start_index = index;
while (true) {
const int next = index_to_start_index_[start_index];
if (start_index == next) break;
start_index = next;
}
// Loop over the variables by increasing ub.
IncrementalSort(
vars_.begin(), vars_.end(), [this](IntegerVariable a, IntegerVariable b) {
return integer_trail_->UpperBound(a) < integer_trail_->UpperBound(b);
});
for (const IntegerVariable var : vars_) {
const IntegerValue lb = integer_trail_->LowerBound(var);
// Second, redo the same thing and make everyone point to the representative.
while (true) {
const int next = index_to_start_index_[index];
if (start_index == next) break;
index_to_start_index_[index] = start_index;
index = next;
}
return start_index;
}
// Check if lb is in an Hall interval, and push it if this is the case.
const int hall_index =
std::lower_bound(hall_ends_.begin(), hall_ends_.end(), lb) -
hall_ends_.begin();
if (hall_index < hall_ends_.size() && hall_starts_[hall_index] <= lb) {
const IntegerValue hs = hall_starts_[hall_index];
const IntegerValue he = hall_ends_[hall_index];
FillHallReason(hs, he);
integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, hs));
if (!integer_trail_->Enqueue(IntegerLiteral::GreaterOrEqual(var, he + 1),
/*literal_reason=*/{}, integer_reason_)) {
bool AllDifferentBoundsPropagator::PropagateLowerBounds() {
// Start by filling the cached bounds and sorting by increasing lb.
for (VarValue& entry : vars_) {
entry.lb = integer_trail_->LowerBound(entry.var);
entry.ub = integer_trail_->UpperBound(entry.var);
}
IncrementalSort(vars_.begin(), vars_.end(),
[](VarValue a, VarValue b) { return a.lb < b.lb; });
// We will split the variable in vars sorted by lb in contiguous subset with
// index of the form [start, start + num_in_window).
int start = 0;
int num_in_window = 1;
// Minimum lower bound in the current window.
IntegerValue min_lb = vars_.front().lb;
const int size = vars_.size();
for (int i = 1; i < size; ++i) {
const IntegerValue lb = vars_[i].lb;
// If the lower bounds of all the other variables is greater, then it can
// never fall into a potential hall interval formed by the variable in the
// current window, so we can split the problem into independent parts.
if (lb <= min_lb + IntegerValue(num_in_window - 1)) {
++num_in_window;
continue;
}
// Process the current window.
if (num_in_window > 1) {
absl::Span<VarValue> window(&vars_[start], num_in_window);
if (!PropagateLowerBoundsInternal(min_lb, window)) {
return false;
}
}
// Updates critical_intervals_. Note that we use the old lb, but that
// doesn't change the value of newly_covered. This block is what takes the
// most time.
int64 newly_covered;
const auto it =
critical_intervals_.GrowRightByOne(lb.value(), &newly_covered);
to_insert_.push_back({newly_covered, var});
const IntegerValue end(it->end);
// Start of the next window.
start = i;
num_in_window = 1;
min_lb = lb;
}
// Take care of the last window.
if (num_in_window > 1) {
absl::Span<VarValue> window(&vars_[start], num_in_window);
return PropagateLowerBoundsInternal(min_lb, window);
}
return true;
}
bool AllDifferentBoundsPropagator::PropagateLowerBoundsInternal(
IntegerValue min_lb, absl::Span<VarValue> vars) {
hall_starts_.clear();
hall_ends_.clear();
// All cached lb in vars will be in [min_lb, min_lb + vars_.size()).
// Make sure we change our base_ so that GetIndex() fit in our buffers.
base_ = min_lb - IntegerValue(1);
// Sparse cleaning of value_to_nodes_.
for (const int i : indices_to_clear_) {
index_to_var_[i] = kNoIntegerVariable;
}
indices_to_clear_.clear();
// Sort vars by increasing ub.
std::sort(vars.begin(), vars.end(),
[](VarValue a, VarValue b) { return a.ub < b.ub; });
for (const VarValue entry : vars) {
const IntegerVariable var = entry.var;
// Note that it is important to use the cache to make sure GetIndex() is
// not out of bound in case integer_trail_->LowerBound() changed when we
// pushed something.
const IntegerValue lb = entry.lb;
const int lb_index = GetIndex(lb);
const bool value_is_covered = PointIsPresent(lb_index);
// Check if lb is in an Hall interval, and push it if this is the case.
if (value_is_covered) {
const int hall_index =
std::lower_bound(hall_ends_.begin(), hall_ends_.end(), lb) -
hall_ends_.begin();
if (hall_index < hall_ends_.size() && hall_starts_[hall_index] <= lb) {
const IntegerValue hs = hall_starts_[hall_index];
const IntegerValue he = hall_ends_[hall_index];
FillHallReason(hs, he);
integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, hs));
if (!integer_trail_->Enqueue(
IntegerLiteral::GreaterOrEqual(var, he + 1),
/*literal_reason=*/{}, integer_reason_)) {
return false;
}
}
}
// Update our internal representation of the non-consecutive intervals.
//
// If lb is not used, we add a node there, otherwise we add it to the
// right of the interval that contains lb. In both cases, if there is an
// interval to the left (resp. right) we merge them.
int new_index = lb_index;
int start_index = lb_index;
int end_index = lb_index;
if (value_is_covered) {
start_index = FindStartIndexAndCompressPath(new_index);
new_index = index_to_end_index_[start_index] + 1;
end_index = new_index;
} else {
if (PointIsPresent(new_index - 1)) {
start_index = FindStartIndexAndCompressPath(new_index - 1);
}
}
if (PointIsPresent(new_index + 1)) {
end_index = index_to_end_index_[new_index + 1];
index_to_start_index_[new_index + 1] = start_index;
}
// Update the end of the representative.
index_to_end_index_[start_index] = end_index;
// This is the only place where we "add" a new node.
{
index_to_start_index_[new_index] = start_index;
index_to_var_[new_index] = var;
indices_to_clear_.push_back(new_index);
}
// We cannot have a conflict, because it should have beend detected before
// by pushing an interval lower bound past its upper bound.
//
// TODO(user): Not 100% clear since pushing can have side-effect, maybe we
// should just report the conflict if it happens!
const IntegerValue end = GetValue(end_index);
DCHECK_LE(end, integer_trail_->UpperBound(var));
// If we have a new Hall interval, add it to the set. Note that it will
// always be last, and if it overlaps some previous Hall intervals, it
// always overlaps them fully.
if (end == integer_trail_->UpperBound(var)) {
const IntegerValue start(it->start);
//
// Note: It is okay to not use entry.ub here if we want to fetch the last
// value, but in practice it shouldn't really change when we push a
// lower_bound and it is faster to use the cached entry.
if (end == entry.ub) {
const IntegerValue start = GetValue(start_index);
while (!hall_starts_.empty() && start <= hall_starts_.back()) {
hall_starts_.pop_back();
hall_ends_.pop_back();
@@ -511,8 +626,8 @@ bool AllDifferentBoundsPropagator::PropagateLowerBounds() {
void AllDifferentBoundsPropagator::RegisterWith(
GenericLiteralWatcher* watcher) {
const int id = watcher->Register(this);
for (const IntegerVariable& var : vars_) {
watcher->WatchIntegerVariable(var, id);
for (const VarValue entry : vars_) {
watcher->WatchIntegerVariable(entry.var, id);
}
watcher->NotifyThatPropagatorMayNotReachFixedPointInOnePass(id);
}

81
ortools/sat/all_different.h
View File

@@ -24,7 +24,6 @@
#include "ortools/sat/integer.h"
#include "ortools/sat/model.h"
#include "ortools/sat/sat_base.h"
#include "ortools/util/sorted_interval_list.h"
namespace operations_research {
namespace sat {
@@ -140,14 +139,11 @@ class AllDifferentConstraint : PropagatorInterface {
// deduce that the domain of the other variables cannot contains such Hall
// interval.
//
// We use a "simple" O(n log n) algorithm.
// We use a "fast" O(n log n) algorithm.
//
// TODO(user): implement the faster algorithm described in:
// TODO(user): It might be difficult to find something faster than what is
// implemented here. Some related reference:
// https://cs.uwaterloo.ca/~vanbeek/Publications/ijcai03_TR.pdf
// Note that the algorithms are similar, the gain comes by replacing our
// SortedDisjointIntervalList with a more customized class for our operations.
// It is even possible to get an O(n) complexity if the values of the bounds are
// in a range of size O(n).
class AllDifferentBoundsPropagator : public PropagatorInterface {
public:
AllDifferentBoundsPropagator(const std::vector<IntegerVariable>& vars,
@@ -157,36 +153,71 @@ class AllDifferentBoundsPropagator : public PropagatorInterface {
void RegisterWith(GenericLiteralWatcher* watcher);
private:
// We locally cache the lb/ub for faster sorting and to guarantee some
// invariant when we push bounds.
struct VarValue {
IntegerVariable var;
IntegerValue lb;
IntegerValue ub;
};
// Fills integer_reason_ with the reason why we have the given hall interval.
void FillHallReason(IntegerValue hall_lb, IntegerValue hall_ub);
// Do half the job of Propagate().
// Do half the job of Propagate(). This will split the variable into
// independent subset, and call PropagateLowerBoundsInternal() on each of
// them.
bool PropagateLowerBounds();
bool PropagateLowerBoundsInternal(IntegerValue min_lb,
absl::Span<VarValue> vars);
// Internally, we will maintain a set of non-consecutive integer intervals of
// the form [start, end]. Each point (i.e. IntegerValue) of such interval will
// be associated to an unique variable and via an union-find algorithm point
// to its start. The end only make sense for representative.
//
// TODO(user): Because we don't use rank, we have a worst case complexity of
// O(n log n). We could try a normal Union-find data structure, but then we
// also have to maintain a start vector.
//
// Note that during the execution of the algorithm we start from empty
// intervals and finish with a set of points of size num_vars.
//
// The list of all points are maintained in the dense vectors index_to_*_
// where we have remapped values to indices (with GetIndex()) to make sure it
// always fall into the correct range.
int FindStartIndexAndCompressPath(int index);
int GetIndex(IntegerValue value) const {
DCHECK_GE(value, base_);
DCHECK_LT(value - base_, index_to_start_index_.size());
return (value - base_).value();
}
IntegerValue GetValue(int index) const { return base_ + IntegerValue(index); }
bool PointIsPresent(int index) const {
return index_to_var_[index] != kNoIntegerVariable;
}
std::vector<IntegerVariable> vars_;
std::vector<IntegerVariable> negated_vars_;
IntegerTrail* integer_trail_;
// The sets of "critical" intervals. This has the same meaning as in the
// disjunctive constraint.
SortedDisjointIntervalList critical_intervals_;
// These vector will be either sorted by lb or by ub.
std::vector<VarValue> vars_;
std::vector<VarValue> negated_vars_;
// The list of Hall intervalls detected so far, sorted.
std::vector<IntegerValue> hall_starts_;
std::vector<IntegerValue> hall_ends_;
// Members needed for explaining the propagation.
//
// The IntegerVariable in an hall interval [lb, ub] are the variables with key
// in [lb, ub] in this map. Note(user): if the set of bounds is small, we
// could use a std::vector here. The O(ub - lb) to create the reason is fine
// since this is the size of the reason.
//
// Optimization: we only insert the entry in the map lazily when the reason
// is needed.
int64 num_calls_;
std::vector<std::pair<int64, IntegerVariable>> to_insert_;
absl::flat_hash_map<int64, IntegerVariable> value_to_variable_;
// Non-consecutive intervals related data-structures.
IntegerValue base_;
std::vector<int> indices_to_clear_;
std::vector<int> index_to_start_index_;
std::vector<int> index_to_end_index_;
std::vector<IntegerVariable> index_to_var_;
// Temporary integer reason.
std::vector<IntegerLiteral> integer_reason_;
DISALLOW_COPY_AND_ASSIGN(AllDifferentBoundsPropagator);

120
ortools/sat/cp_model_lns.cc
View File

@@ -61,6 +61,18 @@ NeighborhoodGeneratorHelper::NeighborhoodGeneratorHelper(
}
}
}
type_to_constraints_.clear();
const int num_constraints = model_proto_.constraints_size();
for (int c = 0; c < num_constraints; ++c) {
const int type = model_proto_.constraints(c).constraint_case();
CHECK_GE(type, 0) << "Negative constraint case ??";
CHECK_LT(type, 10000) << "Large constraint case ??";
if (type >= type_to_constraints_.size()) {
type_to_constraints_.resize(type + 1);
}
type_to_constraints_[type].push_back(c);
}
}
bool NeighborhoodGeneratorHelper::IsActive(int var) const {
@@ -111,20 +123,26 @@ CpModelProto NeighborhoodGeneratorHelper::RelaxGivenVariables(
return FixGivenVariables(initial_solution, fixed_variables);
}
CpModelProto SimpleNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, int64 seed,
double difficulty) const {
namespace {
void GetRandomSubset(int seed, double relative_size, std::vector<int>* base) {
random_engine_t random;
random.seed(seed);
// TODO(user): we could generate this more efficiently than using random
// shuffle.
std::vector<int> fixed_variables = helper_.ActiveVariables();
std::shuffle(base->begin(), base->end(), random);
const int target_size = std::ceil(relative_size * base->size());
base->resize(target_size);
}
std::shuffle(fixed_variables.begin(), fixed_variables.end(), random);
const int target_size =
std::ceil((1.0 - difficulty) * fixed_variables.size());
fixed_variables.resize(target_size);
} // namespace
CpModelProto SimpleNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, int64 seed,
double difficulty) const {
std::vector<int> fixed_variables = helper_.ActiveVariables();
GetRandomSubset(seed, 1.0 - difficulty, &fixed_variables);
return helper_.FixGivenVariables(initial_solution, fixed_variables);
}
@@ -242,5 +260,91 @@ CpModelProto ConstraintGraphNeighborhoodGenerator::Generate(
return helper_.RelaxGivenVariables(initial_solution, relaxed_variables);
}
CpModelProto SchedulingNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, int64 seed,
double difficulty) const {
std::set<int> intervals_to_relax;
{
const auto span = helper_.TypeToConstraints(ConstraintProto::kInterval);
std::vector<int> v(span.begin(), span.end());
GetRandomSubset(seed, difficulty, &v);
intervals_to_relax.insert(v.begin(), v.end());
}
CpModelProto copy = helper_.ModelProto();
// Fix the presence/absence of non-relaxed intervals.
for (const int i : helper_.TypeToConstraints(ConstraintProto::kInterval)) {
if (intervals_to_relax.count(i)) continue;
const ConstraintProto& interval_ct = copy.constraints(i);
if (interval_ct.enforcement_literal().empty()) continue;
CHECK_EQ(interval_ct.enforcement_literal().size(), 1);
const int enforcement_ref = interval_ct.enforcement_literal(0);
const int enforcement_var = PositiveRef(enforcement_ref);
const int value = initial_solution.solution(enforcement_var);
// Fix the value.
copy.mutable_variables(enforcement_var)->clear_domain();
copy.mutable_variables(enforcement_var)->add_domain(value);
copy.mutable_variables(enforcement_var)->add_domain(value);
// If the interval is ignored, skip for the loop below as there is no
// point adding precedence on it.
if (RefIsPositive(enforcement_ref) == (value == 0)) {
intervals_to_relax.insert(i);
}
}
for (const int c : helper_.TypeToConstraints(ConstraintProto::kNoOverlap)) {
// Sort all non-relaxed intervals of this constraint by current start time.
std::vector<std::pair<int64, int>> start_interval_pairs;
for (const int i : copy.constraints(c).no_overlap().intervals()) {
if (intervals_to_relax.count(i)) continue;
const ConstraintProto& interval_ct = copy.constraints(i);
// TODO(user): we ignore size zero for now.
const int size_var = interval_ct.interval().size();
if (initial_solution.solution(size_var) == 0) continue;
const int start_var = interval_ct.interval().start();
const int64 start_value = initial_solution.solution(start_var);
start_interval_pairs.push_back({start_value, i});
}
std::sort(start_interval_pairs.begin(), start_interval_pairs.end());
// Add precedence between the remaining intervals, forcing their order.
for (int i = 0; i + 1 < start_interval_pairs.size(); ++i) {
const int before_var =
copy.constraints(start_interval_pairs[i].second).interval().end();
const int after_var = copy.constraints(start_interval_pairs[i + 1].second)
.interval()
.start();
CHECK_LE(initial_solution.solution(before_var),
initial_solution.solution(after_var));
LinearConstraintProto* linear = copy.add_constraints()->mutable_linear();
linear->add_domain(kint64min);
linear->add_domain(0);
linear->add_vars(before_var);
linear->add_coeffs(1);
linear->add_vars(after_var);
linear->add_coeffs(-1);
}
}
// Set the current solution as a hint.
//
// TODO(user): Move to common function?
copy.clear_solution_hint();
for (int var = 0; var < copy.variables_size(); ++var) {
copy.mutable_solution_hint()->add_vars(var);
copy.mutable_solution_hint()->add_values(initial_solution.solution(var));
}
return copy;
}
} // namespace sat
} // namespace operations_research

26
ortools/sat/cp_model_lns.h
View File

@@ -16,6 +16,7 @@
#include <vector>
#include "absl/types/span.h"
#include "ortools/base/integral_types.h"
#include "ortools/sat/cp_model.pb.h"
@@ -60,6 +61,13 @@ class NeighborhoodGeneratorHelper {
return var_to_constraint_;
}
// Returns all the constraints indices of a given type.
const absl::Span<const int> TypeToConstraints(
ConstraintProto::ConstraintCase type) const {
if (type >= type_to_constraints_.size()) return {};
return absl::MakeSpan(type_to_constraints_[type]);
}
// The initial problem.
const CpModelProto& ModelProto() const { return model_proto_; }
@@ -69,6 +77,9 @@ class NeighborhoodGeneratorHelper {
const CpModelProto& model_proto_;
// Constraints by types.
std::vector<std::vector<int>> type_to_constraints_;
// Variable-Constraint graph.
std::vector<std::vector<int>> constraint_to_var_;
std::vector<std::vector<int>> var_to_constraint_;
@@ -150,6 +161,21 @@ class ConstraintGraphNeighborhoodGenerator : public NeighborhoodGenerator {
double difficulty) const final;
};
// Only make sense for scheduling problem. This select a random set of interval
// of the problem according to the difficulty. Then, for each no_overlap
// constraints, it adds strict relation order between the non-relaxed intervals.
//
// TODO(user): Also deal with cumulative constraint.
class SchedulingNeighborhoodGenerator : public NeighborhoodGenerator {
public:
explicit SchedulingNeighborhoodGenerator(
NeighborhoodGeneratorHelper const* helper, const std::string& name)
: NeighborhoodGenerator(name, helper) {}
CpModelProto Generate(const CpSolverResponse& initial_solution, int64 seed,
double difficulty) const final;
};
} // namespace sat
} // namespace operations_research

32
ortools/sat/cp_model_solver.cc
View File

@@ -497,9 +497,9 @@ std::string CpSolverResponseStats(const CpSolverResponse& response) {
absl::StrAppend(&result, "\nobjective: NA");
absl::StrAppend(&result, "\nbest_bound: NA");
} else {
absl::StrAppendFormat(&result, "\nobjective: %g",
absl::StrAppendFormat(&result, "\nobjective: %.9g",
response.objective_value());
absl::StrAppendFormat(&result, "\nbest_bound: %g",
absl::StrAppendFormat(&result, "\nbest_bound: %.9g",
response.best_objective_bound());
}
@@ -513,10 +513,10 @@ std::string CpSolverResponseStats(const CpSolverResponse& response) {
"\npropagations: ", response.num_binary_propagations());
absl::StrAppend(
&result, "\ninteger_propagations: ", response.num_integer_propagations());
absl::StrAppendFormat(&result, "\nwalltime: %g", response.wall_time());
absl::StrAppendFormat(&result, "\nusertime: %g", response.user_time());
absl::StrAppendFormat(&result, "\ndeterministic_time: %g",
response.deterministic_time());
absl::StrAppend(&result, "\nwalltime: ", response.wall_time());
absl::StrAppend(&result, "\nusertime: ", response.user_time());
absl::StrAppend(&result,
"\ndeterministic_time: ", response.deterministic_time());
absl::StrAppend(&result, "\n");
return result;
}
@@ -1900,8 +1900,16 @@ CpSolverResponse SolveCpModelWithLNS(
const bool focus_on_decision_variables =
parameters->lns_focus_on_decision_variables();
// TODO(user): Find a way to propagate the level zero bounds from the other
// worker inside this base LNS problem.
const NeighborhoodGeneratorHelper helper(&model_proto,
focus_on_decision_variables);
// For now we will just alternate between our possible neighborhoods.
NeighborhoodGeneratorHelper helper(&model_proto, focus_on_decision_variables);
//
// TODO(user): select with higher probability the one that seems to work best?
// We can evaluate this compared to the best solution at the time of the
// neighborhood creation.
std::vector<std::unique_ptr<NeighborhoodGenerator>> generators;
generators.push_back(
absl::make_unique<SimpleNeighborhoodGenerator>(&helper, "rnd_lns"));
@@ -1909,9 +1917,17 @@ CpSolverResponse SolveCpModelWithLNS(
&helper, "var_lns"));
generators.push_back(absl::make_unique<ConstraintGraphNeighborhoodGenerator>(
&helper, "cst_lns"));
if (!helper.TypeToConstraints(ConstraintProto::kNoOverlap).empty()) {
generators.push_back(absl::make_unique<SchedulingNeighborhoodGenerator>(
&helper, "scheduling_lns"));
}
// The "optimal" difficulties do not have to be the same for different
// generators. TODO(user): move this inside the generator API?
// generators.
//
// TODO(user): This should be shared across LNS threads! create a thread
// safe class that accept signals of the form (difficulty and result) and
// properly update the difficulties.
std::vector<AdaptiveParameterValue> difficulties(generators.size(),
AdaptiveParameterValue(0.5));

260
ortools/sat/diffn.cc
View File

@@ -38,18 +38,6 @@ NonOverlappingRectanglesPropagator::NonOverlappingRectanglesPropagator(
strict_(strict),
integer_trail_(integer_trail) {
CHECK_GT(num_boxes_, 0);
neighbors_.resize(num_boxes_ * (num_boxes_ - 1));
neighbors_begins_.resize(num_boxes_);
neighbors_ends_.resize(num_boxes_);
for (int i = 0; i < num_boxes_; ++i) {
const int begin = i * (num_boxes_ - 1);
neighbors_begins_[i] = begin;
neighbors_ends_[i] = begin + (num_boxes_ - 1);
for (int j = 0; j < num_boxes_; ++j) {
if (j == i) continue;
neighbors_[begin + (j > i ? j - 1 : j)] = j;
}
}
}
NonOverlappingRectanglesPropagator::~NonOverlappingRectanglesPropagator() {}
@@ -68,7 +56,7 @@ bool NonOverlappingRectanglesPropagator::Propagate() {
y_.SetTimeDirection(true);
for (int box = 0; box < num_boxes_; ++box) {
if (!strict_ && cached_areas_[box] == 0) continue;
if (cached_areas_[box] == 0) continue;
RETURN_IF_FALSE(FailWhenEnergyIsTooLarge(box));
}
@@ -80,9 +68,6 @@ void NonOverlappingRectanglesPropagator::RegisterWith(
const int id = watcher->Register(this);
x_.WatchAllTasks(id, watcher);
y_.WatchAllTasks(id, watcher);
for (int t = 0; t < num_boxes_; ++t) {
watcher->RegisterReversibleInt(id, &neighbors_ends_[t]);
}
watcher->SetPropagatorPriority(id, 3);
}
@@ -106,57 +91,40 @@ IntegerValue DistanceToBoundingInterval(IntegerValue min_a, IntegerValue max_a,
} // namespace
void NonOverlappingRectanglesPropagator::UpdateNeighbors(int box) {
tmp_removed_.clear();
void NonOverlappingRectanglesPropagator::SortNeighbors(int box) {
cached_distance_to_bounding_box_.resize(num_boxes_);
neighbors_.clear();
const IntegerValue box_x_min = x_.StartMin(box);
const IntegerValue box_x_max = x_.EndMax(box);
const IntegerValue box_y_min = y_.StartMin(box);
const IntegerValue box_y_max = y_.EndMax(box);
int new_index = neighbors_begins_[box];
const int end = neighbors_ends_[box];
for (int i = new_index; i < end; ++i) {
const int other = neighbors_[i];
for (int other = 0; other < num_boxes_; ++other) {
if (other == box) continue;
if (cached_areas_[other] == 0) continue;
const IntegerValue other_x_min = x_.StartMin(other);
const IntegerValue other_x_max = x_.EndMax(other);
if (IntervalAreDisjointForSure(box_x_min, box_x_max, other_x_min,
other_x_max)) {
tmp_removed_.push_back(other);
continue;
}
const IntegerValue other_y_min = y_.StartMin(other);
const IntegerValue other_y_max = y_.EndMax(other);
if (IntervalAreDisjointForSure(box_y_min, box_y_max, other_y_min,
other_y_max)) {
tmp_removed_.push_back(other);
continue;
}
neighbors_[new_index++] = other;
neighbors_.push_back(other);
cached_distance_to_bounding_box_[other] =
std::max(DistanceToBoundingInterval(box_x_min, box_x_max, other_x_min,
other_x_max),
DistanceToBoundingInterval(box_y_min, box_y_max, other_y_min,
other_y_max));
}
neighbors_ends_[box] = new_index;
for (int i = 0; i < tmp_removed_.size();) {
neighbors_[new_index++] = tmp_removed_[i++];
}
IncrementalSort(neighbors_.begin() + neighbors_begins_[box],
neighbors_.begin() + neighbors_ends_[box],
[this](int i, int j) {
return cached_distance_to_bounding_box_[i] <
cached_distance_to_bounding_box_[j];
});
IncrementalSort(neighbors_.begin(), neighbors_.begin(), [this](int i, int j) {
return cached_distance_to_bounding_box_[i] <
cached_distance_to_bounding_box_[j];
});
}
bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) {
// Note that we only consider the smallest dimension of each boxes here.
UpdateNeighbors(box);
SortNeighbors(box);
IntegerValue area_min_x = x_.StartMin(box);
IntegerValue area_max_x = x_.EndMax(box);
IntegerValue area_min_y = y_.StartMin(box);
@@ -165,9 +133,7 @@ bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) {
IntegerValue sum_of_areas = cached_areas_[box];
IntegerValue total_sum_of_areas = sum_of_areas;
const int end = neighbors_ends_[box];
for (int i = neighbors_begins_[box]; i < end; ++i) {
const int other = neighbors_[i];
for (const int other : neighbors_) {
total_sum_of_areas += cached_areas_[other];
}
@@ -182,9 +148,9 @@ bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) {
// TODO(user): Is there a better order, maybe sort by distance
// with the current box.
for (int i = neighbors_begins_[box]; i < end; ++i) {
for (int i = 0; i < neighbors_.size(); ++i) {
const int other = neighbors_[i];
if (cached_areas_[other] == 0) continue;
CHECK_GT(cached_areas_[other], 0);
// Update Bounding box.
area_min_x = std::min(area_min_x, x_.StartMin(other));
@@ -205,10 +171,8 @@ bool NonOverlappingRectanglesPropagator::FailWhenEnergyIsTooLarge(int box) {
x_.ClearReason();
y_.ClearReason();
add_box_energy_in_rectangle_reason(box);
for (int j = neighbors_begins_[box]; j <= i; ++j) {
const int other = neighbors_[j];
if (cached_areas_[other] == 0) continue;
add_box_energy_in_rectangle_reason(other);
for (int j = 0; j <= i; ++j) {
add_box_energy_in_rectangle_reason(neighbors_[j]);
}
x_.ImportOtherReasons(y_);
return x_.ReportConflict();
@@ -290,7 +254,6 @@ bool CheckOverload(bool time_direction, IntegerValue other_time,
other->AddEndMinReason(task, other_time + 1);
}
}
// LOG(INFO) << "Overload: " << num_events - critical_event;
helper->ImportOtherReasons(*other);
return helper->ReportConflict();
}
@@ -298,12 +261,30 @@ bool CheckOverload(bool time_direction, IntegerValue other_time,
return true;
}
void AddOtherReasons(const std::vector<int>& tasks, IntegerValue other_time,
int main_task, SchedulingConstraintHelper* other) {
other->ClearReason();
bool main_task_seen = false;
for (const int task : tasks) {
other->AddStartMaxReason(task, other_time);
other->AddEndMinReason(task, other_time + 1);
if (task == main_task) {
main_task_seen = true;
}
}
if (!main_task_seen) {
other->AddStartMaxReason(main_task, other_time);
other->AddEndMinReason(main_task, other_time + 1);
}
}
bool DetectPrecedences(bool time_direction, IntegerValue other_time,
const absl::flat_hash_set<int>& boxes,
const absl::flat_hash_set<int>& active_boxes,
SchedulingConstraintHelper* helper,
SchedulingConstraintHelper* other) {
helper->SetTimeDirection(time_direction);
other->SetTimeDirection(true);
const auto& task_by_increasing_end_min = helper->TaskByIncreasingEndMin();
const auto& task_by_decreasing_start_max = helper->TaskByDecreasingStartMax();
@@ -314,7 +295,7 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time,
const int t = task_time.task_index;
const IntegerValue end_min = task_time.time;
if (helper->IsAbsent(t) || !gtl::ContainsKey(boxes, t)) continue;
if (helper->IsAbsent(t) || !gtl::ContainsKey(active_boxes, t)) continue;
while (queue_index >= 0) {
const auto to_insert = task_by_decreasing_start_max[queue_index];
@@ -322,7 +303,7 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time,
const IntegerValue start_max = to_insert.time;
if (end_min <= start_max) break;
if (gtl::ContainsKey(boxes, task_index)) {
if (gtl::ContainsKey(active_boxes, task_index)) {
CHECK(helper->IsPresent(task_index));
task_set.AddEntry({task_index, helper->ShiftedStartMin(task_index),
helper->DurationMin(task_index)});
@@ -343,44 +324,22 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time,
if (end_min_of_critical_tasks > helper->StartMin(t)) {
const std::vector<TaskSet::Entry>& sorted_tasks = task_set.SortedTasks();
helper->ClearReason();
other->ClearReason();
// We need:
// - StartMax(ct) < EndMin(t) for the detectable precedence.
// - StartMin(ct) > window_start for the end_min_of_critical_tasks reason.
const IntegerValue window_start = sorted_tasks[critical_index].start_min;
const int num_conflicting_tasks = sorted_tasks.size() - critical_index;
if (num_conflicting_tasks < 4) {
for (int i = critical_index; i < sorted_tasks.size(); ++i) {
const int ct = sorted_tasks[i].task;
if (ct == t) continue;
other->AddReasonForBeingBefore(t, ct);
other->AddReasonForBeingBefore(ct, t);
}
for (int i = critical_index; i + 1 < sorted_tasks.size(); ++i) {
const int ct1 = sorted_tasks[i].task;
if (ct1 == t) continue;
for (int j = i + 1; j < sorted_tasks.size(); ++j) {
const int ct2 = sorted_tasks[j].task;
if (ct2 == t) continue;
other->AddReasonForBeingBefore(ct1, ct2);
other->AddReasonForBeingBefore(ct2, ct1);
}
}
} else {
for (int i = critical_index; i < sorted_tasks.size(); ++i) {
const int ct = sorted_tasks[i].task;
if (ct == t) continue;
other->AddStartMaxReason(ct, other_time);
other->AddEndMinReason(ct, other_time + 1);
}
other->AddStartMaxReason(t, other_time);
other->AddEndMinReason(t, other_time + 1);
std::vector<int> tasks(sorted_tasks.size() - critical_index);
for (int i = critical_index; i < sorted_tasks.size(); ++i) {
tasks[i - critical_index] = sorted_tasks[i].task;
}
AddOtherReasons(tasks, other_time, t, other);
for (int i = critical_index; i < sorted_tasks.size(); ++i) {
const int ct = sorted_tasks[i].task;
if (ct == t) continue;
CHECK(gtl::ContainsKey(boxes, ct));
CHECK(gtl::ContainsKey(active_boxes, ct));
CHECK(helper->IsPresent(ct));
helper->AddEnergyAfterReason(ct, sorted_tasks[i].duration_min,
window_start);
@@ -390,7 +349,7 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time,
// Add the reason for t (we only need the end-min).
helper->AddEndMinReason(t, end_min);
// Import other reasons.
// Import reasons from the 'other' dimension.
helper->ImportOtherReasons(*other);
// This augment the start-min of t and subsequently it can augment the
@@ -409,7 +368,6 @@ bool DetectPrecedences(bool time_direction, IntegerValue other_time,
}
return true;
}
} // namespace
bool NonOverlappingRectanglesPropagator::PropagateOnProjections() {
@@ -472,7 +430,26 @@ bool NonOverlappingRectanglesPropagator::FindMandatoryBoxesOnOneDimension(
}
for (const auto& it : mandatory_boxes) {
RETURN_IF_FALSE(PropagateMandatoryBoxesOnOneDimension(it.first, it.second,
// Compute the 'canonical' line to use when explaining that boxes overlap
// on the 'other' dimension.
IntegerValue lb(kint64min);
IntegerValue ub(kint64max);
for (const int task : it.second) {
lb = std::max(lb, other->StartMax(task));
ub = std::min(ub, other->EndMin(task));
}
const IntegerValue span = ub - lb + 1;
IntegerValue selected = lb;
for (int shift = 30; shift >= 0; --shift) {
const IntegerValue mask(static_cast<int64>(1) >> shift);
if (mask <= span) {
selected = (ub / mask) * mask;
break;
}
}
// And propagate.
RETURN_IF_FALSE(PropagateMandatoryBoxesOnOneDimension(selected, it.second,
helper, other));
}
return true;
@@ -485,109 +462,6 @@ bool NonOverlappingRectanglesPropagator::PropagateMandatoryBoxesOnOneDimension(
RETURN_IF_FALSE(CheckOverload(true, event, active_boxes, helper, other));
RETURN_IF_FALSE(DetectPrecedences(true, event, active_boxes, helper, other));
RETURN_IF_FALSE(DetectPrecedences(false, event, active_boxes, helper, other));
// RETURN_IF_FALSE(NaivePush(true, event, active_boxes, helper, other));
// RETURN_IF_FALSE(NaivePush(false, event, active_boxes, helper, other));
return true;
}
bool NonOverlappingRectanglesPropagator::NaivePush(
bool time_direction, IntegerValue other_time,
const absl::flat_hash_set<int>& boxes, SchedulingConstraintHelper* helper,
SchedulingConstraintHelper* other) {
// const auto left_box_before_right_box =
// [](int left, int right, SchedulingConstraintHelper* helper) {
// // left box pushes right box.
// const IntegerValue left_end_min = helper->EndMin(left);
// if (left_end_min > helper->StartMin(right)) {
// // Store reasons state.
// const int literal_size = helper->MutableLiteralReason()->size();
// const int integer_size = helper->MutableIntegerReason()->size();
//
// helper->AddEndMinReason(left, left_end_min);
// if (!helper->IncreaseStartMin(right, left_end_min)) {
// return false;
// }
//
// // Restore the reasons to the state before the increase of the
// start. helper->MutableLiteralReason()->resize(literal_size);
// helper->MutableIntegerReason()->resize(integer_size);
// }
//
// // right box pushes left box.
// const IntegerValue right_start_max = helper->StartMax(right);
// if (right_start_max < helper->EndMax(left)) {
// helper->AddStartMaxReason(right, right_start_max);
// return helper->DecreaseEndMax(left, right_start_max);
// }
//
// return true;
// };
//
// helper->SetTimeDirection(time_direction);
// other->SetTimeDirection(true);
// const auto& task_by_increasing_end_min =
// helper->TaskByIncreasingEndMin(); const auto&
// task_by_decreasing_start_max =
// helper->TaskByDecreasingStartMax();
//
// int queue_index = helper->NumTasks() - 1;
// for (const auto task_time : task_by_increasing_end_min) {
// const int t = task_time.task_index;
// const IntegerValue end_min = task_time.time;
//
// if (!gtl::ContainsKey(boxes, t)) continue;
//
//
// for (int b1 = 0; b1 + 1 < boxes.size(); ++b1) {
// const int box1 = boxes[b1];
// if (!strict_ && cached_areas_[box1] == 0) continue;
// for (int b2 = b1 + 1; b2 < boxes.size(); ++b2) {
// const int box2 = boxes[b2];
// if (!strict_ && cached_areas_[box2] == 0) continue;
// // For each direction and each order, we test if the boxes can be
// // disjoint.
// const int state = (helper->EndMin(box1) <= helper->StartMax(box2)) +
// 2 * (helper->EndMin(box2) <=
// helper->StartMax(box1));
//
// if (state == 3) continue;
// // Clean up reasons.
// helper->ClearReason();
// other->ClearReason();
//
// // This is an "hack" to be able to easily test for none or for one
// // and only one of the conditions below.
// switch (state) {
// case 0: {
// helper->AddReasonForBeingBefore(box1, box2);
// helper->AddReasonForBeingBefore(box2, box1);
// other->AddReasonForBeingBefore(box1, box2);
// other->AddReasonForBeingBefore(box2, box1);
// helper->ImportOtherReasons(*other);
// return helper->ReportConflict();
// }
// case 1: {
// other->AddReasonForBeingBefore(box1, box2);
// other->AddReasonForBeingBefore(box2, box1);
// helper->AddReasonForBeingBefore(box1, box2);
// helper->ImportOtherReasons(*other);
// RETURN_IF_FALSE(left_box_before_right_box(box1, box2, helper));
// break;
// }
// case 2: {
// other->AddReasonForBeingBefore(box1, box2);
// other->AddReasonForBeingBefore(box2, box1);
// helper->AddReasonForBeingBefore(box2, box1);
// helper->ImportOtherReasons(*other);
// RETURN_IF_FALSE(left_box_before_right_box(box2, box1, helper));
// break;
// }
// default: {
// break;
// }
// }
// }
// }
return true;
}

12
ortools/sat/diffn.h
View File

@@ -52,11 +52,7 @@ class NonOverlappingRectanglesPropagator : public PropagatorInterface {
const std::vector<int>& boxes,
SchedulingConstraintHelper* active,
SchedulingConstraintHelper* other);
bool NaivePush(bool time_direction, IntegerValue other_time,
const absl::flat_hash_set<int>& boxes,
SchedulingConstraintHelper* helper,
SchedulingConstraintHelper* other_helper);
void UpdateNeighbors(int box);
void SortNeighbors(int box);
bool FailWhenEnergyIsTooLarge(int box);
const int num_boxes_;
@@ -65,13 +61,7 @@ class NonOverlappingRectanglesPropagator : public PropagatorInterface {
const bool strict_;
IntegerTrail* integer_trail_;
// The neighbors_ of a box will be in
// [neighbors_[begin[box]], neighbors_[end[box]])
std::vector<int> neighbors_;
std::vector<int> neighbors_begins_;
std::vector<int> neighbors_ends_;
std::vector<int> tmp_removed_;
std::vector<IntegerValue> cached_areas_;
std::vector<IntegerValue> cached_distance_to_bounding_box_;

243
ortools/sat/integer_search.cc
View File

@@ -25,25 +25,94 @@
namespace operations_research {
namespace sat {
LiteralIndex AtMinValue(IntegerVariable var, Model* model) {
IntegerEncoder* const integer_encoder = model->GetOrCreate<IntegerEncoder>();
IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
DCHECK(!integer_trail->IsCurrentlyIgnored(var));
const IntegerValue lb = integer_trail->LowerBound(var);
DCHECK_LE(lb, integer_trail->UpperBound(var));
if (lb == integer_trail->UpperBound(var)) return kNoLiteralIndex;
return integer_encoder
->GetOrCreateAssociatedLiteral(IntegerLiteral::LowerOrEqual(var, lb))
.Index();
}
LiteralIndex GreaterOrEqualToMiddleValue(IntegerVariable var, Model* model) {
auto* encoder = model->GetOrCreate<IntegerEncoder>();
auto* trail = model->GetOrCreate<Trail>();
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
const IntegerValue var_lb = integer_trail->LowerBound(var);
const IntegerValue var_ub = integer_trail->UpperBound(var);
CHECK_LT(var_lb, var_ub);
const IntegerValue chosen_value =
var_lb + std::max(IntegerValue(1), (var_ub - var_lb) / IntegerValue(2));
const Literal ge = encoder->GetOrCreateAssociatedLiteral(
IntegerLiteral::GreaterOrEqual(var, chosen_value));
CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable()));
VLOG(2) << "Chosen " << var << " >= " << chosen_value;
return ge.Index();
}
LiteralIndex SplitDomainUsingLpValue(IntegerVariable var, Model* model) {
auto* encoder = model->GetOrCreate<IntegerEncoder>();
auto* trail = model->GetOrCreate<Trail>();
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
auto* lp_dispatcher = model->GetOrCreate<LinearProgrammingDispatcher>();
const IntegerVariable positive_var =
VariableIsPositive(var) ? var : NegationOf(var);
DCHECK(!integer_trail->IsCurrentlyIgnored(positive_var));
LinearProgrammingConstraint* lp =
gtl::FindWithDefault(*lp_dispatcher, positive_var, nullptr);
if (lp == nullptr) return kNoLiteralIndex;
const IntegerValue value = IntegerValue(
static_cast<int64>(std::round(lp->GetSolutionValue(positive_var))));
// Because our lp solution might be from higher up in the tree, it
// is possible that value is now outside the domain of positive_var.
// In this case, we just revert to the current literal.
const IntegerValue lb = integer_trail->LowerBound(positive_var);
const IntegerValue ub = integer_trail->UpperBound(positive_var);
// We try first (<= value), but if this do not reduce the domain we
// try to enqueue (>= value). Note that even for domain with hole,
// since we know that this variable is not fixed, one of the two
// alternative must reduce the domain.
//
// TODO(user): use GetOrCreateLiteralAssociatedToEquality() instead?
// It may replace two decision by only one. However this function
// cannot currently be called during search, but that should be easy
// enough to fix.
if (value >= lb && value < ub) {
const Literal le = encoder->GetOrCreateAssociatedLiteral(
IntegerLiteral::LowerOrEqual(positive_var, value));
CHECK(!trail->Assignment().VariableIsAssigned(le.Variable()));
return le.Index();
}
if (value > lb && value <= ub) {
const Literal ge = encoder->GetOrCreateAssociatedLiteral(
IntegerLiteral::GreaterOrEqual(positive_var, value));
CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable()));
return ge.Index();
}
return kNoLiteralIndex;
}
// TODO(user): the complexity caused by the linear scan in this heuristic and
// the one below is ok when search_branching is set to SAT_SEARCH because it is
// not executed often, but otherwise it is done for each search decision,
// which seems expensive. 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>();
return [integer_encoder, integer_trail, /*copy*/ vars]() {
return [/*copy*/ vars, model]() {
IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
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();
}
const LiteralIndex decision = AtMinValue(var, model);
if (decision != kNoLiteralIndex) return decision;
}
return kNoLiteralIndex;
};
@@ -51,9 +120,8 @@ std::function<LiteralIndex()> FirstUnassignedVarAtItsMinHeuristic(
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]() {
return [/*copy */ vars, model]() {
IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
IntegerVariable candidate = kNoIntegerVariable;
IntegerValue candidate_lb;
for (const IntegerVariable var : vars) {
@@ -66,10 +134,7 @@ std::function<LiteralIndex()> UnassignedVarWithLowestMinAtItsMinHeuristic(
}
}
if (candidate == kNoIntegerVariable) return kNoLiteralIndex;
return integer_encoder
->GetOrCreateAssociatedLiteral(
IntegerLiteral::LowerOrEqual(candidate, candidate_lb))
.Index();
return AtMinValue(candidate, model);
};
}
@@ -96,10 +161,6 @@ std::function<LiteralIndex()> SatSolverHeuristic(Model* model) {
}
std::function<LiteralIndex()> PseudoCost(Model* model) {
auto* encoder = model->GetOrCreate<IntegerEncoder>();
auto* trail = model->GetOrCreate<Trail>();
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
const ObjectiveSynchronizationHelper* helper =
model->Get<ObjectiveSynchronizationHelper>();
const bool has_objective =
@@ -111,29 +172,14 @@ std::function<LiteralIndex()> PseudoCost(Model* model) {
PseudoCosts* pseudo_costs = model->GetOrCreate<PseudoCosts>();
// NOTE: The algorithm to choose the value for a variable is kept outside of
// PseudoCosts class since this is an independent heuristic. For now this is
// only used for pseudo costs however this can be refactored a bit for other
// branching strategies to use.
return [=]() {
const IntegerVariable chosen_var = pseudo_costs->GetBestDecisionVar();
if (chosen_var == kNoIntegerVariable) return kNoLiteralIndex;
const IntegerValue chosen_var_lb = integer_trail->LowerBound(chosen_var);
const IntegerValue chosen_var_ub = integer_trail->UpperBound(chosen_var);
CHECK_LT(chosen_var_lb, chosen_var_ub);
// TODO(user): Experiment with different heuristics for choosing
// value.
const IntegerValue chosen_value =
chosen_var_lb +
std::max(IntegerValue(1),
(chosen_var_ub - chosen_var_lb) / IntegerValue(2));
const Literal ge = encoder->GetOrCreateAssociatedLiteral(
IntegerLiteral::GreaterOrEqual(chosen_var, chosen_value));
CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable()));
VLOG(2) << "Chosen " << chosen_var << " >= " << chosen_value;
return ge.Index();
return GreaterOrEqualToMiddleValue(chosen_var, model);
};
}
@@ -214,16 +260,50 @@ std::function<LiteralIndex()> FollowHint(
};
}
std::function<LiteralIndex()> ExploitLpSolution(
std::function<LiteralIndex()> heuristic, Model* model) {
LiteralIndex ExploitLpSolution(const LiteralIndex decision, Model* model) {
auto* encoder = model->GetOrCreate<IntegerEncoder>();
auto* trail = model->GetOrCreate<Trail>();
auto* integer_trail = model->GetOrCreate<IntegerTrail>();
auto* lp_dispatcher = model->GetOrCreate<LinearProgrammingDispatcher>();
auto* lp_constraints =
model->GetOrCreate<LinearProgrammingConstraintCollection>();
const SatParameters& parameters = *(model->GetOrCreate<SatParameters>());
if (decision == kNoLiteralIndex) {
return decision;
}
bool lp_solution_is_exploitable = true;
// TODO(user,user): When we have more than one LP, their set of variable
// is always disjoint. So we could still change the polarity if the next
// variable we branch on is part of a LP that has a solution.
for (LinearProgrammingConstraint* lp : *lp_constraints) {
if (!lp->HasSolution() ||
!(parameters.exploit_all_lp_solution() || lp->SolutionIsInteger())) {
lp_solution_is_exploitable = false;
break;
}
}
if (lp_solution_is_exploitable) {
for (const IntegerLiteral l :
encoder->GetAllIntegerLiterals(Literal(decision))) {
const IntegerVariable positive_var =
VariableIsPositive(l.var) ? l.var : NegationOf(l.var);
if (integer_trail->IsCurrentlyIgnored(positive_var)) continue;
const LiteralIndex lp_decision =
SplitDomainUsingLpValue(positive_var, model);
if (lp_decision != kNoLiteralIndex) {
return lp_decision;
}
}
}
return decision;
}
std::function<LiteralIndex()> ExploitLpSolution(
std::function<LiteralIndex()> heuristic, Model* model) {
auto* lp_constraints =
model->GetOrCreate<LinearProgrammingConstraintCollection>();
// Use the normal heuristic if the LP(s) do not seem to cover enough variables
// to be relevant.
// TODO(user): Instead, try and depending on the result call it again or not?
@@ -237,86 +317,9 @@ std::function<LiteralIndex()> ExploitLpSolution(
if (num_integer_variables > 2 * num_lp_variables) return heuristic;
}
bool last_decision_followed_lp = false;
int old_level = 0;
double old_obj = 0.0;
return [=]() mutable {
const LiteralIndex decision = heuristic();
if (decision == kNoLiteralIndex) {
if (last_decision_followed_lp) {
VLOG(1) << "Integer LP solution is feasible, level:" << old_level
<< "->" << trail->CurrentDecisionLevel() << " obj:" << old_obj;
}
last_decision_followed_lp = false;
return kNoLiteralIndex;
}
bool lp_solution_is_exploitable = true;
double obj = 0.0;
// TODO(user,user): When we have more than one LP, their set of variable
// is always disjoint. So we could still change the polarity if the next
// variable we branch on is part of a LP that has a solution.
for (LinearProgrammingConstraint* lp : *lp_constraints) {
if (!lp->HasSolution() ||
!(parameters.exploit_all_lp_solution() || lp->SolutionIsInteger())) {
lp_solution_is_exploitable = false;
break;
}
obj += lp->SolutionObjectiveValue();
}
if (lp_solution_is_exploitable) {
if (!last_decision_followed_lp || obj != old_obj) {
old_level = trail->CurrentDecisionLevel();
old_obj = obj;
VLOG(2) << "Integer LP solution at level:" << old_level
<< " obj:" << old_obj;
}
for (const IntegerLiteral l :
encoder->GetAllIntegerLiterals(Literal(decision))) {
const IntegerVariable positive_var =
VariableIsPositive(l.var) ? l.var : NegationOf(l.var);
if (integer_trail->IsCurrentlyIgnored(positive_var)) continue;
LinearProgrammingConstraint* lp =
gtl::FindWithDefault(*lp_dispatcher, positive_var, nullptr);
if (lp != nullptr) {
const IntegerValue value = IntegerValue(static_cast<int64>(
std::round(lp->GetSolutionValue(positive_var))));
// Because our lp solution might be from higher up in the tree, it
// is possible that value is now outside the domain of positive_var.
// In this case, we just revert to the current literal.
const IntegerValue lb = integer_trail->LowerBound(positive_var);
const IntegerValue ub = integer_trail->UpperBound(positive_var);
// We try first (<= value), but if this do not reduce the domain we
// try to enqueue (>= value). Note that even for domain with hole,
// since we know that this variable is not fixed, one of the two
// alternative must reduce the domain.
//
// TODO(user): use GetOrCreateLiteralAssociatedToEquality() instead?
// It may replace two decision by only one. However this function
// cannot currently be called during search, but that should be easy
// enough to fix.
if (value >= lb && value < ub) {
const Literal le = encoder->GetOrCreateAssociatedLiteral(
IntegerLiteral::LowerOrEqual(positive_var, value));
CHECK(!trail->Assignment().VariableIsAssigned(le.Variable()));
last_decision_followed_lp = true;
return le.Index();
}
if (value > lb && value <= ub) {
const Literal ge = encoder->GetOrCreateAssociatedLiteral(
IntegerLiteral::GreaterOrEqual(positive_var, value));
CHECK(!trail->Assignment().VariableIsAssigned(ge.Variable()));
last_decision_followed_lp = true;
return ge.Index();
}
}
}
}
last_decision_followed_lp = false;
return decision;
return ExploitLpSolution(decision, model);
};
}
@@ -598,7 +601,7 @@ SatSolver::Status SolveIntegerProblemWithLazyEncoding(Model* model) {
void LogNewSolution(const std::string& event_or_solution_count,
double time_in_seconds, double obj_lb, double obj_ub,
const std::string& solution_info) {
LOG(INFO) << absl::StrFormat("#%-5s %6.2fs obj:[%g,%g] %s",
LOG(INFO) << absl::StrFormat("#%-5s %6.2fs obj:[%.9g,%.9g] %s",
event_or_solution_count, time_in_seconds, obj_lb,
obj_ub, solution_info);
}

18
ortools/sat/integer_search.h
View File

@@ -17,11 +17,24 @@
#include <vector>
#include "ortools/sat/integer.h"
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_solver.h"
namespace operations_research {
namespace sat {
// Returns decision corresponding to var at its lower bound. Returns
// kNoLiteralIndex if the variable is fixed.
LiteralIndex AtMinValue(IntegerVariable var, Model* model);
// Returns decision corresponding to var <= round(lp_value). If the variable
// does not appear in the LP, this method returns kNoLiteralIndex.
LiteralIndex SplitDomainUsingLpValue(IntegerVariable var, Model* model);
// Returns decision corresponding to var >= lb + max(1, (ub - lb) / 2). It also
// CHECKs that the variable is not fixed.
LiteralIndex GreaterOrEqualToMiddleValue(IntegerVariable var, 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
@@ -75,6 +88,11 @@ std::function<LiteralIndex()> PseudoCost(Model* model);
std::function<LiteralIndex()> ExploitLpSolution(
std::function<LiteralIndex()> heuristic, Model* model);
// Similar to ExploitLpSolution(). Takes LiteralIndex as base decision and
// changes change the returned decision to AtLpValue() of the underlying integer
// variable if LP solution is exploitable.
LiteralIndex ExploitLpSolution(const LiteralIndex decision, Model* model);
// A restart policy that restarts every k failures.
std::function<bool()> RestartEveryKFailures(int k, SatSolver* solver);

7
ortools/util/sorted_interval_list.cc
View File

@@ -139,9 +139,12 @@ bool Domain::IsEmpty() const { return intervals_.empty(); }
int64 Domain::Size() const {
int64 size = 0;
for (const ClosedInterval interval : intervals_) {
size += operations_research::CapAdd(
1, operations_research::CapSub(interval.end, interval.start));
size = operations_research::CapAdd(
size, operations_research::CapSub(interval.end, interval.start));
}
// Because the intervals are closed on both side above, with miss 1 per
// interval.
size = operations_research::CapAdd(size, intervals_.size());
return size;
}