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

[CP-SAT] cleanup diffn propagation; remove dead code

Browse Source
This commit is contained in:
Laurent Perron
2024-01-31 14:44:02 +01:00
parent d25ac6a27d
commit 8d4b32967c
13 changed files with 609 additions and 374 deletions
Download Patch File Download Diff File Expand all files Collapse all files

283
ortools/sat/2d_orthogonal_packing.cc
View File

@@ -23,9 +23,9 @@
#include "absl/log/check.h"
#include "absl/numeric/bits.h"
#include "absl/random/distributions.h"
#include "absl/types/span.h"
#include "ortools/base/logging.h"
#include "ortools/base/vlog_is_on.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/util.h"
#include "ortools/util/bitset.h"
@@ -220,47 +220,19 @@ std::vector<int> GetDffConflict(
return result;
}
// Check for conflict all combinations of the two Dual Feasible Functions f_0
// (see documentation for GetDffConflict()) and f_2 (see documentation for
// RoundingDualFeasibleFunction). More precisely, check whether there exist `l`
// and `k` so that
//
// sum_i f_2^k(f_0^l(sizes_x[i])) * sizes_y[i] > f_2^k(f_0^l(x_bb_size)) *
// y_bb_size
//
// The function returns the smallest subset of items enough to make the
// inequality above true or an empty vector if impossible.
std::vector<int> CheckFeasibilityWithDualFunction2(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
absl::Span<const int> index_by_decreasing_x_size, IntegerValue x_bb_size,
IntegerValue y_bb_size) {
if (x_bb_size == 1) {
return {};
}
std::vector<IntegerValue> sizes_x_rescaled;
if (x_bb_size >= std::numeric_limits<uint16_t>::max()) {
// To do fast division we want our sizes to fit in a uint16_t. The simplest
// way of doing that is to just first apply this DFF with the right
// power-of-two value of the parameter.
const int log2_k =
absl::bit_width(static_cast<uint64_t>(x_bb_size.value() + 1)) - 16 + 1;
const RoundingDualFeasibleFunctionPowerOfTwo dff(x_bb_size, log2_k);
sizes_x_rescaled.resize(sizes_x.size());
for (int i = 0; i < sizes_x.size(); i++) {
sizes_x_rescaled[i] = dff(sizes_x[i]);
}
x_bb_size = dff(x_bb_size);
CHECK_LT(x_bb_size, std::numeric_limits<uint16_t>::max());
sizes_x = sizes_x_rescaled;
}
} // namespace
// We now want to find the minimum set of values of `k` that would always
// find a conflict if there is a `k` for which it exists. In the literature
// it is often implied (but not stated) that it is sufficient to test the
// values of `k` that correspond to the size of an item. This is not true. To
// find the minimum set of values of `k` we look for all values of `k` that
// are "extreme": ie., the rounding on the division truncates the most (or the
void OrthogonalPackingInfeasibilityDetector::GetAllCandidatesForKForDff2(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y, IntegerValue x_bb_size,
IntegerValue sqrt_bb_size, IntegerValue y_bb_size,
Bitset64<IntegerValue>& candidates) {
// We want to find the minimum set of values of `k` that would always find a
// conflict if there is a `k` for which it exists. In the literature it is
// often implied (but not stated) that it is sufficient to test the values of
// `k` that correspond to the size of an item. This is not true. To find the
// minimum set of values of `k` we look for all values of `k` that are
// "extreme": ie., the rounding on the division truncates the most (or the
// least) amount, depending on the sign it appears in the formula.
IntegerValue sum_widths = 0;
for (int i = 0; i < sizes_x.size(); i++) {
@@ -273,10 +245,9 @@ std::vector<int> CheckFeasibilityWithDualFunction2(
}
const IntegerValue round_up = sum_widths > 2 * y_bb_size ? 0 : 1;
// x_bb_size is less than 65536, so this fits in only 4kib.
Bitset64<IntegerValue> candidates(x_bb_size / 2 + 2);
candidates.ClearAndResize(x_bb_size / 2 + 2);
// `sqrt_bb_size` is lower than 256.
const IntegerValue sqrt_bb_size = FloorSquareRoot(x_bb_size.value());
for (IntegerValue i = 2; i <= sqrt_bb_size; i++) {
candidates.Set(i);
}
@@ -296,9 +267,6 @@ std::vector<int> CheckFeasibilityWithDualFunction2(
}
}
std::vector<int> result;
int num_items = sizes_x.size();
// Remove some bogus candidates added by the logic above.
candidates.Clear(0);
candidates.Clear(1);
@@ -318,6 +286,89 @@ std::vector<int> CheckFeasibilityWithDualFunction2(
if (x_bb_size > 3) {
candidates.Set(x_bb_size / 3 + 1);
}
}
// Check for conflict all combinations of the two Dual Feasible Functions f_0
// (see documentation for GetDffConflict()) and f_2 (see documentation for
// RoundingDualFeasibleFunction). More precisely, check whether there exist `l`
// and `k` so that
//
// sum_i f_2^k(f_0^l(sizes_x[i])) * sizes_y[i] > f_2^k(f_0^l(x_bb_size)) *
// y_bb_size
//
// The function returns the smallest subset of items enough to make the
// inequality above true or an empty vector if impossible.
std::vector<int>
OrthogonalPackingInfeasibilityDetector::CheckFeasibilityWithDualFunction2(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
absl::Span<const int> index_by_decreasing_x_size, IntegerValue x_bb_size,
IntegerValue y_bb_size, int max_number_of_parameters_to_check) {
if (x_bb_size == 1) {
return {};
}
std::vector<IntegerValue> sizes_x_rescaled;
if (x_bb_size >= std::numeric_limits<uint16_t>::max()) {
// To do fast division we want our sizes to fit in a uint16_t. The simplest
// way of doing that is to just first apply this DFF with the right
// power-of-two value of the parameter.
const int log2_k =
absl::bit_width(static_cast<uint64_t>(x_bb_size.value() + 1)) - 16 + 1;
const RoundingDualFeasibleFunctionPowerOfTwo dff(x_bb_size, log2_k);
sizes_x_rescaled.resize(sizes_x.size());
for (int i = 0; i < sizes_x.size(); i++) {
sizes_x_rescaled[i] = dff(sizes_x[i]);
}
x_bb_size = dff(x_bb_size);
CHECK_LT(x_bb_size, std::numeric_limits<uint16_t>::max());
sizes_x = sizes_x_rescaled;
}
Bitset64<IntegerValue> candidates;
const IntegerValue sqrt_bb_size = FloorSquareRoot(x_bb_size.value());
int num_items = sizes_x.size();
const IntegerValue max_possible_number_of_parameters =
std::min(x_bb_size / 4 + 1, sqrt_bb_size * num_items);
if (5ull * max_number_of_parameters_to_check <
max_possible_number_of_parameters) {
// There are many more possible values than what we want to sample. It is
// not worth to pay the price of computing all optimal values to drop most
// of them, so let's just pick it randomly.
candidates.Resize(x_bb_size / 4 + 1);
int num_candidates = 0;
while (num_candidates < max_number_of_parameters_to_check) {
const IntegerValue pick =
absl::Uniform(random_, 1, x_bb_size.value() / 4);
if (!candidates.IsSet(pick)) {
candidates.Set(pick);
num_candidates++;
}
}
} else {
GetAllCandidatesForKForDff2(sizes_x, sizes_y, x_bb_size, sqrt_bb_size,
y_bb_size, candidates);
if (max_number_of_parameters_to_check < max_possible_number_of_parameters) {
// We might have produced too many candidates. Let's count them and if it
// is the case, sample them.
int count = 0;
for (auto it = candidates.begin(); it != candidates.end(); ++it) {
count++;
}
if (count > max_number_of_parameters_to_check) {
std::vector<IntegerValue> sampled_candidates(
max_number_of_parameters_to_check);
std::sample(candidates.begin(), candidates.end(),
sampled_candidates.begin(),
max_number_of_parameters_to_check, random_);
candidates.ClearAll();
for (const IntegerValue k : sampled_candidates) {
candidates.Set(k);
}
}
}
}
std::vector<int> result;
// Finally run our small loop to look for the conflict!
std::vector<IntegerValue> modified_sizes(num_items);
@@ -341,21 +392,13 @@ std::vector<int> CheckFeasibilityWithDualFunction2(
return result;
}
} // namespace
OrthogonalPackingInfeasibilityDetector::Result
OrthogonalPackingInfeasibilityDetector::TestFeasibility(
OrthogonalPackingResult
OrthogonalPackingInfeasibilityDetector::TestFeasibilityImpl(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
std::pair<IntegerValue, IntegerValue> bounding_box_size) {
enum class ConflictType {
NO_CONFLICT,
TRIVIAL,
DFF_F0,
DFF_F2,
};
ConflictType conflict_type = ConflictType::NO_CONFLICT;
num_calls_++;
std::pair<IntegerValue, IntegerValue> bounding_box_size,
const OrthogonalPackingOptions& options) {
using ConflictType = OrthogonalPackingResult::ConflictType;
const int num_items = sizes_x.size();
DCHECK_EQ(num_items, sizes_y.size());
@@ -363,6 +406,11 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
bounding_box_size.first * bounding_box_size.second;
IntegerValue total_energy = 0;
auto make_item = [sizes_x, sizes_y](int i) {
return OrthogonalPackingResult::Item{
.index = i, .size_x = sizes_x[i], .size_y = sizes_y[i]};
};
index_by_decreasing_x_size_.resize(num_items);
index_by_decreasing_y_size_.resize(num_items);
for (int i = 0; i < num_items; i++) {
@@ -371,14 +419,16 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
index_by_decreasing_y_size_[i] = i;
if (sizes_x[i] > bounding_box_size.first ||
sizes_y[i] > bounding_box_size.second) {
num_conflicts_++;
return Result{.result = Status::INFEASIBLE,
.minimum_unfeasible_subproblem = {i}};
OrthogonalPackingResult result(
OrthogonalPackingResult::Status::INFEASIBLE);
result.conflict_type_ = ConflictType::TRIVIAL;
result.items_participating_on_conflict_ = {make_item(i)};
return result;
}
}
if (num_items <= 1) {
return Result{.result = Status::FEASIBLE};
return OrthogonalPackingResult(OrthogonalPackingResult::Status::FEASIBLE);
}
std::sort(index_by_decreasing_x_size_.begin(),
@@ -389,25 +439,27 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
[&sizes_y](int a, int b) { return sizes_y[a] > sizes_y[b]; });
// First look for pairwise incompatible pairs.
if (options_.use_pairwise) {
if (options.use_pairwise) {
if (auto pair = FindPairwiseConflict(sizes_x, sizes_y, bounding_box_size,
index_by_decreasing_x_size_,
index_by_decreasing_y_size_);
pair.has_value()) {
num_conflicts_++;
num_conflicts_two_items_++;
return Result{.result = Status::INFEASIBLE,
.minimum_unfeasible_subproblem = {pair.value().first,
pair.value().second}};
OrthogonalPackingResult result(
OrthogonalPackingResult::Status::INFEASIBLE);
result.conflict_type_ = ConflictType::PAIRWISE;
result.items_participating_on_conflict_ = {
make_item(pair.value().first), make_item(pair.value().second)};
return result;
}
if (num_items == 2) {
return Result{.result = Status::FEASIBLE};
return OrthogonalPackingResult(OrthogonalPackingResult::Status::FEASIBLE);
}
}
std::vector<int> result;
OrthogonalPackingResult result(OrthogonalPackingResult::Status::UNKNOWN);
if (total_energy > bb_area) {
conflict_type = ConflictType::TRIVIAL;
result.conflict_type_ = ConflictType::TRIVIAL;
result.result_ = OrthogonalPackingResult::Status::INFEASIBLE;
std::vector<std::pair<int, IntegerValue>> index_to_energy;
index_to_energy.reserve(num_items);
for (int i = 0; i < num_items; i++) {
@@ -422,26 +474,39 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
for (int i = 0; i < index_to_energy.size(); i++) {
recomputed_energy += index_to_energy[i].second;
if (recomputed_energy > bb_area) {
result.resize(i + 1);
result.items_participating_on_conflict_.resize(i + 1);
for (int j = 0; j <= i; j++) {
result[j] = index_to_energy[j].first;
result.items_participating_on_conflict_[j] =
make_item(index_to_energy[j].first);
}
result.slack_ = recomputed_energy - bb_area - 1;
break;
}
}
}
auto set_conflict_if_better = [&result, &conflict_type](
const int minimum_conflict_size = options.use_pairwise ? 3 : 2;
if (result.items_participating_on_conflict_.size() == minimum_conflict_size) {
return result;
}
auto set_conflict_if_better = [&result, &make_item](
const std::vector<int>& conflict,
ConflictType type) {
if (!conflict.empty() &&
(result.empty() || conflict.size() < result.size())) {
conflict_type = type;
result = conflict;
(result.result_ != OrthogonalPackingResult::Status::INFEASIBLE ||
conflict.size() < result.items_participating_on_conflict_.size())) {
result.result_ = OrthogonalPackingResult::Status::INFEASIBLE;
result.conflict_type_ = type;
result.items_participating_on_conflict_.clear();
for (int i : conflict) {
result.items_participating_on_conflict_.push_back(make_item(i));
}
result.slack_ = 0; // Only supported for trivial for now.
}
};
if (options_.use_dff_f0) {
if (options.use_dff_f0) {
// If there is no pairwise incompatible pairs, this DFF cannot find a
// conflict by enlarging a item on both x and y directions: this would
// create an item as long as the whole box and another item as high as the
@@ -460,29 +525,56 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
set_conflict_if_better(conflict, ConflictType::DFF_F0);
}
if (options_.use_dff_f2) {
if (result.items_participating_on_conflict_.size() == minimum_conflict_size) {
return result;
}
if (options.use_dff_f2) {
// We only check for conflicts applying this DFF on heights and widths, but
// not on both, which would be too expensive if done naively.
auto conflict = CheckFeasibilityWithDualFunction2(
sizes_x, sizes_y, index_by_decreasing_x_size_, bounding_box_size.first,
bounding_box_size.second);
bounding_box_size.second,
options.dff2_max_number_of_parameters_to_check);
set_conflict_if_better(conflict, ConflictType::DFF_F2);
if (result.items_participating_on_conflict_.size() ==
minimum_conflict_size) {
return result;
}
conflict = CheckFeasibilityWithDualFunction2(
sizes_y, sizes_x, index_by_decreasing_y_size_, bounding_box_size.second,
bounding_box_size.first);
bounding_box_size.first,
options.dff2_max_number_of_parameters_to_check);
set_conflict_if_better(conflict, ConflictType::DFF_F2);
}
if (!result.empty()) {
return result;
}
OrthogonalPackingResult OrthogonalPackingInfeasibilityDetector::TestFeasibility(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
std::pair<IntegerValue, IntegerValue> bounding_box_size,
const OrthogonalPackingOptions& options) {
using ConflictType = OrthogonalPackingResult::ConflictType;
num_calls_++;
OrthogonalPackingResult result =
TestFeasibilityImpl(sizes_x, sizes_y, bounding_box_size, options);
if (result.result_ == OrthogonalPackingResult::Status::INFEASIBLE) {
num_conflicts_++;
switch (conflict_type) {
switch (result.conflict_type_) {
case ConflictType::DFF_F0:
num_conflicts_dff0_++;
break;
case ConflictType::DFF_F2:
num_conflicts_dff2_++;
break;
case ConflictType::PAIRWISE:
num_conflicts_two_items_++;
break;
case ConflictType::TRIVIAL:
// The total area of the items was larger than the area of the box.
num_trivial_conflicts_++;
@@ -491,10 +583,27 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
LOG(FATAL) << "Should never happen";
break;
}
return Result{.result = Status::INFEASIBLE,
.minimum_unfeasible_subproblem = result};
}
return Result{.result = Status::UNKNOWN};
return result;
}
bool OrthogonalPackingResult::TryUseSlackToReduceItemSize(
int i, Coord coord, IntegerValue lower_bound) {
Item& item = items_participating_on_conflict_[i];
IntegerValue& size = coord == Coord::kCoordX ? item.size_x : item.size_y;
const IntegerValue orthogonal_size =
coord == Coord::kCoordX ? item.size_y : item.size_x;
if (size <= lower_bound || orthogonal_size > slack_) {
return false;
}
const IntegerValue new_size =
std::max(lower_bound, size - slack_ / orthogonal_size);
slack_ -= (size - new_size) * orthogonal_size;
DCHECK_NE(size, new_size);
DCHECK_GE(slack_, 0);
size = new_size;
return true;
}
} // namespace sat

102
ortools/sat/2d_orthogonal_packing.h
View File

@@ -20,6 +20,7 @@
#include <vector>
#include "absl/log/check.h"
#include "absl/random/bit_gen_ref.h"
#include "absl/types/span.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/synchronization.h"
@@ -31,6 +32,60 @@ struct OrthogonalPackingOptions {
bool use_pairwise = true;
bool use_dff_f0 = true;
bool use_dff_f2 = true;
int dff2_max_number_of_parameters_to_check = std::numeric_limits<int>::max();
};
class OrthogonalPackingResult {
public:
explicit OrthogonalPackingResult() : result_(Status::UNKNOWN) {}
enum class Status {
INFEASIBLE,
FEASIBLE,
UNKNOWN,
};
Status GetResult() const { return result_; }
struct Item {
// Index of the item on the original sizes_x/sizes_y input.
int index;
// New size for item of index `i` which is smaller or equal to the initial
// size. The subproblem remain infeasible if every item is shrinked to its
// new size.
IntegerValue size_x;
IntegerValue size_y;
};
const std::vector<Item>& GetItemsParticipatingOnConflict() const {
return items_participating_on_conflict_;
}
bool HasSlack() const { return slack_ > IntegerValue(0); }
enum class Coord { kCoordX, kCoordY };
// Use an eventual slack to reduce the size of item corresponding to the
// `i`-th element on GetItemsParticipatingOnConflict(). It will not use any
// slack to reduce it beyond lower_bound. This is a no-op if HasSlack() is
// false.
bool TryUseSlackToReduceItemSize(int i, Coord coord,
IntegerValue lower_bound = 0);
private:
friend class OrthogonalPackingInfeasibilityDetector;
explicit OrthogonalPackingResult(Status result) : result_(result) {}
enum class ConflictType {
NO_CONFLICT,
TRIVIAL,
PAIRWISE,
DFF_F0,
DFF_F2,
};
Status result_;
ConflictType conflict_type_ = ConflictType::NO_CONFLICT;
IntegerValue slack_ = 0;
std::vector<Item> items_participating_on_conflict_;
};
// Class for solving the orthogonal packing problem when it can be done
@@ -38,29 +93,43 @@ struct OrthogonalPackingOptions {
class OrthogonalPackingInfeasibilityDetector {
public:
explicit OrthogonalPackingInfeasibilityDetector(
SharedStatistics* shared_stats,
const OrthogonalPackingOptions& options = OrthogonalPackingOptions())
: options_(options), shared_stats_(shared_stats) {}
absl::BitGenRef random, SharedStatistics* shared_stats)
: random_(random), shared_stats_(shared_stats) {}
~OrthogonalPackingInfeasibilityDetector();
enum class Status {
INFEASIBLE,
FEASIBLE,
UNKNOWN,
};
struct Result {
Status result;
std::vector<int> minimum_unfeasible_subproblem;
};
Result TestFeasibility(
OrthogonalPackingResult TestFeasibility(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
std::pair<IntegerValue, IntegerValue> bounding_box_size);
std::pair<IntegerValue, IntegerValue> bounding_box_size,
const OrthogonalPackingOptions& options = OrthogonalPackingOptions());
private:
const OrthogonalPackingOptions options_;
OrthogonalPackingResult TestFeasibilityImpl(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
std::pair<IntegerValue, IntegerValue> bounding_box_size,
const OrthogonalPackingOptions& options);
// Returns a minimum set of values of the parameter `k` of f_2^k that is
// sufficient to find a conflict. This function runs in
// O(num_items * sqrt(bb_size)) operations.
// All sizes must be positive values less than UINT16_MAX.
// The returned bitset will contain less elements than
// min(sqrt_bb_size * num_items, x_bb_size/4+1).
void GetAllCandidatesForKForDff2(absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
IntegerValue x_bb_size,
IntegerValue sqrt_bb_size,
IntegerValue y_bb_size,
Bitset64<IntegerValue>& candidates);
std::vector<int> CheckFeasibilityWithDualFunction2(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
absl::Span<const int> index_by_decreasing_x_size, IntegerValue x_bb_size,
IntegerValue y_bb_size, int max_number_of_parameters_to_check);
std::vector<int> index_by_decreasing_x_size_;
std::vector<int> index_by_decreasing_y_size_;
@@ -71,6 +140,7 @@ class OrthogonalPackingInfeasibilityDetector {
int64_t num_conflicts_dff2_ = 0;
int64_t num_conflicts_dff0_ = 0;
absl::BitGenRef random_;
SharedStatistics* shared_stats_;
};

5
ortools/sat/BUILD.bazel
View File

@@ -266,6 +266,7 @@ cc_library(
":cp_model_utils",
":integer",
":integer_search",
":linear_propagation",
":model",
":sat_base",
":sat_parameters_cc_proto",
@@ -1823,6 +1824,8 @@ cc_library(
"@com_google_absl//absl/log",
"@com_google_absl//absl/log:check",
"@com_google_absl//absl/numeric:bits",
"@com_google_absl//absl/random:bit_gen_ref",
"@com_google_absl//absl/random:distributions",
"@com_google_absl//absl/types:span",
],
)
@@ -1930,7 +1933,6 @@ cc_library(
":subsolver",
":synchronization",
":util",
"//ortools/base",
"//ortools/base:stl_util",
"//ortools/graph:connected_components",
"//ortools/util:adaptative_parameter_value",
@@ -1943,6 +1945,7 @@ cc_library(
"@com_google_absl//absl/base:core_headers",
"@com_google_absl//absl/container:flat_hash_map",
"@com_google_absl//absl/container:flat_hash_set",
"@com_google_absl//absl/log",
"@com_google_absl//absl/log:check",
"@com_google_absl//absl/meta:type_traits",
"@com_google_absl//absl/random:bit_gen_ref",

1
ortools/sat/cp_model_lns.cc
View File

@@ -38,7 +38,6 @@
#include "absl/types/span.h"
#include "ortools/base/logging.h"
#include "ortools/base/stl_util.h"
#include "ortools/base/vlog_is_on.h"
#include "ortools/graph/connected_components.h"
#include "ortools/sat/cp_model.pb.h"
#include "ortools/sat/cp_model_mapping.h"

1
ortools/sat/cp_model_presolve.cc
View File

@@ -5438,6 +5438,7 @@ bool CpModelPresolver::PresolveNoOverlap2D(int /*c*/, ConstraintProto* ct) {
return RemoveConstraint(ct);
}
// TODO(user): handle this case. See issue #4068.
if (!has_zero_sized_interval && (x_constant || y_constant)) {
context_->UpdateRuleStats(
"no_overlap_2d: a dimension is constant, splitting into many "

14
ortools/sat/cp_model_search.cc
View File

@@ -34,6 +34,7 @@
#include "ortools/sat/cp_model_utils.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/integer_search.h"
#include "ortools/sat/linear_propagation.h"
#include "ortools/sat/model.h"
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_parameters.pb.h"
@@ -327,12 +328,25 @@ std::function<BooleanOrIntegerLiteral()> ConstructHeuristicSearchStrategy(
if (ModelHasSchedulingConstraints(cp_model_proto)) {
std::vector<std::function<BooleanOrIntegerLiteral()>> heuristics;
const auto& params = *model->GetOrCreate<SatParameters>();
bool possible_new_constraints = false;
if (params.use_dynamic_precedence_in_disjunctive()) {
possible_new_constraints = true;
heuristics.push_back(DisjunctivePrecedenceSearchHeuristic(model));
}
if (params.use_dynamic_precedence_in_cumulative()) {
possible_new_constraints = true;
heuristics.push_back(CumulativePrecedenceSearchHeuristic(model));
}
// Tricky: we need to create this at level zero in case there are no linear
// constraint in the model at the beginning.
//
// TODO(user): Alternatively, support creation of SatPropagator at positive
// level.
if (possible_new_constraints && params.new_linear_propagation()) {
model->GetOrCreate<LinearPropagator>();
}
heuristics.push_back(SchedulingSearchHeuristic(model));
return SequentialSearch(std::move(heuristics));
}

289
ortools/sat/diffn.cc
View File

@@ -246,17 +246,14 @@ NonOverlappingRectanglesEnergyPropagator::
{"NonOverlappingRectanglesEnergyPropagator/called", num_calls_});
stats.push_back(
{"NonOverlappingRectanglesEnergyPropagator/conflicts", num_conflicts_});
stats.push_back(
{"NonOverlappingRectanglesEnergyPropagator/multiple_conflicts",
num_multiple_conflicts_});
stats.push_back(
{"NonOverlappingRectanglesEnergyPropagator/conflicts_two_boxes",
num_conflicts_two_boxes_});
stats.push_back({"NonOverlappingRectanglesEnergyPropagator/refined",
num_refined_conflicts_});
stats.push_back(
{"NonOverlappingRectanglesEnergyPropagator/orthogonal_packing_conflict",
num_orthogonal_packing_conflict_});
{"NonOverlappingRectanglesEnergyPropagator/conflicts_with_slack",
num_conflicts_with_slack_});
shared_stats_->AddStats(stats);
}
@@ -267,7 +264,6 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
if (!x_.SynchronizeAndSetTimeDirection(true)) return false;
if (!y_.SynchronizeAndSetTimeDirection(true)) return false;
num_calls_++;
Rectangle bounding_box = {.x_min = std::numeric_limits<IntegerValue>::max(),
.x_max = std::numeric_limits<IntegerValue>::min(),
.y_min = std::numeric_limits<IntegerValue>::max(),
@@ -310,6 +306,8 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
return true;
}
num_calls_++;
std::optional<Conflict> best_conflict =
FindConflict(std::move(active_box_ranges));
if (!best_conflict.has_value()) {
@@ -321,18 +319,30 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
// search for a best explanation. This is specially the case since we only
// want to re-run it over the items that participate in the conflict, so it is
// a much smaller problem.
IntegerValue best_explanation_size = best_conflict->items.size();
IntegerValue best_explanation_size =
best_conflict->opp_result.GetItemsParticipatingOnConflict().size();
bool refined = false;
while (true) {
const std::optional<Conflict> conflict = FindConflict(best_conflict->items);
std::vector<RectangleInRange> items_participating_in_conflict;
items_participating_in_conflict.reserve(
best_conflict->items_for_opp.size());
for (const auto& item :
best_conflict->opp_result.GetItemsParticipatingOnConflict()) {
items_participating_in_conflict.push_back(
best_conflict->items_for_opp[item.index]);
}
std::optional<Conflict> conflict =
FindConflict(items_participating_in_conflict);
if (!conflict.has_value()) break;
// We prefer an explanation with the least number of boxes.
if (conflict->items.size() >= best_explanation_size) {
if (conflict->opp_result.GetItemsParticipatingOnConflict().size() >=
best_explanation_size) {
// The new explanation isn't better than the old one. Stop trying.
break;
}
best_explanation_size = conflict->items.size();
best_conflict = conflict;
best_explanation_size =
conflict->opp_result.GetItemsParticipatingOnConflict().size();
best_conflict = std::move(conflict);
refined = true;
}
@@ -342,9 +352,6 @@ bool NonOverlappingRectanglesEnergyPropagator::Propagate() {
if (best_explanation_size == 2) {
num_conflicts_two_boxes_++;
}
if (best_conflict->opp_problem_size > 0) {
num_orthogonal_packing_conflict_++;
}
BuildAndReportEnergyTooLarge(generalized_explanation);
return false;
}
@@ -359,28 +366,13 @@ NonOverlappingRectanglesEnergyPropagator::FindConflict(
rectangles_with_too_much_energy.candidates.empty()) {
return std::nullopt;
}
num_multiple_conflicts_ +=
rectangles_with_too_much_energy.conflicts.size() > 1;
std::vector<RectangleInRange> best_explanation;
Rectangle best_rectangle;
for (const Rectangle& r : rectangles_with_too_much_energy.conflicts) {
const std::vector<RectangleInRange> range_for_explanation =
GetEnergyConflictForRectangle(r, active_box_ranges);
CheckPropagationIsValid(range_for_explanation, r);
if (best_explanation.empty() ||
best_explanation.size() > range_for_explanation.size()) {
best_explanation = range_for_explanation;
best_rectangle = r;
}
}
Conflict best_conflict;
// Now try with a DFF on the candidates.
int opp_problem_size = 0;
// Also try the DFF on known conflicts, hopefully we will get a smaller
// explanation.
// Sample 10 rectangles for looking for a conflict with DFF. This avoids
// making this heuristic too costly.
// Sample 10 rectangles (at least five among the ones for which we already
// detected an energy overflow), extract an orthogonal packing subproblem for
// each and look for conflict. Sampling avoids making this heuristic too
// costly.
constexpr int kSampleSize = 10;
absl::InlinedVector<Rectangle, kSampleSize> sampled_rectangles;
std::sample(rectangles_with_too_much_energy.conflicts.begin(),
@@ -430,158 +422,87 @@ NonOverlappingRectanglesEnergyPropagator::FindConflict(
// our rectangle is the bounding box, and we need to fit inside it a set of
// items corresponding to the minimum intersection of the original items
// with this bounding box.
const bool use_expensive_heuristics =
(sizes_x.size() < 10) ||
// Propagating a large conflict is expensive, so it's worth running an
// expensive heuristic to make it smaller.
best_conflict.opp_result.GetResult() ==
OrthogonalPackingResult::Status::INFEASIBLE ||
!rectangles_with_too_much_energy.conflicts.empty();
const auto opp_result = orthogonal_packing_checker_.TestFeasibility(
sizes_x, sizes_y, {r.SizeX(), r.SizeY()});
if (opp_result.result ==
OrthogonalPackingInfeasibilityDetector::Status::INFEASIBLE &&
(best_explanation.empty() ||
best_explanation.size() >
opp_result.minimum_unfeasible_subproblem.size())) {
best_explanation.clear();
for (int idx : opp_result.minimum_unfeasible_subproblem) {
best_explanation.push_back(filtered_items[idx]);
}
best_rectangle = r;
opp_problem_size = static_cast<int>(active_box_ranges.size());
sizes_x, sizes_y, {r.SizeX(), r.SizeY()},
OrthogonalPackingOptions{.use_pairwise = true,
.use_dff_f0 = true,
.use_dff_f2 = true,
.dff2_max_number_of_parameters_to_check =
use_expensive_heuristics ? 25 : 3});
if (opp_result.GetResult() == OrthogonalPackingResult::Status::INFEASIBLE &&
(best_conflict.opp_result.GetResult() !=
OrthogonalPackingResult::Status::INFEASIBLE ||
best_conflict.opp_result.GetItemsParticipatingOnConflict().size() >
opp_result.GetItemsParticipatingOnConflict().size())) {
best_conflict.items_for_opp = filtered_items;
best_conflict.opp_result = opp_result;
best_conflict.rectangle_with_too_much_energy = r;
}
// Use the fact that our rectangles are ordered in shrinking order to remove
// all items that will never contribute again.
filtered_items.swap(active_box_ranges);
}
if (!best_explanation.empty()) {
return Conflict{.items = std::move(best_explanation),
.rectangle_with_too_much_energy = best_rectangle,
.opp_problem_size = opp_problem_size};
if (best_conflict.opp_result.GetResult() ==
OrthogonalPackingResult::Status::INFEASIBLE) {
return best_conflict;
} else {
return std::nullopt;
}
return std::nullopt;
}
std::vector<RectangleInRange>
NonOverlappingRectanglesEnergyPropagator::GeneralizeExplanation(
const Conflict& conflict) {
const Rectangle& rectangle = conflict.rectangle_with_too_much_energy;
IntegerValue available_energy = rectangle.Area();
IntegerValue used_energy = 0;
for (const RectangleInRange& box : conflict.items) {
IntegerValue energy = box.GetMinimumIntersectionArea(rectangle);
used_energy += energy;
}
std::vector<RectangleInRange> ranges_for_explanation(conflict.items);
OrthogonalPackingResult relaxed_result = conflict.opp_result;
IntegerValue slack = used_energy - available_energy - 1;
if (conflict.opp_problem_size > 0) {
VLOG_EVERY_N_SEC(2, 3)
<< "Rectangle with energy overflow after solving OPP: " << rectangle
<< " using " << ranges_for_explanation.size() << "/"
<< conflict.opp_problem_size << " items.";
} else {
VLOG_EVERY_N_SEC(2, 3) << "Rectangle with energy overflow: " << rectangle
<< " (" << used_energy << "/" << available_energy
<< "), with " << ranges_for_explanation.size()
<< " boxes inside.";
}
for (auto& range : ranges_for_explanation) {
Rectangle overlap = range.GetMinimumIntersection(rectangle);
DCHECK_GT(overlap.SizeX(), 0);
DCHECK_GT(overlap.SizeY(), 0);
range = RectangleInRange::BiggestWithMinIntersection(
rectangle, range, overlap.SizeX(), overlap.SizeY());
// Now use the energy slack to make the intervals even bigger.
if (overlap.SizeX() < slack) {
IntegerValue y_slack = slack / overlap.SizeX();
// We know we can remove an area as big as y_slack * x_size. That means
// we can reduce the start of an interval of y_slack or increase the end
// by the same amount. But there is no point wasting the slack on a
// interval larger than the zero-level interval, so we check for that.
const IntegerValue y_min_slack = std::max(
IntegerValue(0),
std::min(y_slack, range.bounding_area.y_min -
y_.LevelZeroStartMin(range.box_index)));
if (y_min_slack > 0) {
range.bounding_area.y_min -= y_min_slack;
slack -= y_min_slack * overlap.SizeX();
y_slack -= y_min_slack;
}
const IntegerValue y_max_slack =
std::min(y_slack, y_.LevelZeroStartMax(range.box_index) -
range.bounding_area.y_max);
if (y_max_slack > 0) {
range.bounding_area.y_max += y_max_slack;
slack -= y_max_slack * overlap.SizeX();
}
// Recompute the overlap, since we changed the item dimensions.
overlap = range.GetMinimumIntersection(rectangle);
}
if (overlap.SizeY() < slack) {
IntegerValue x_slack = slack / overlap.SizeY();
const IntegerValue x_min_slack = std::max(
IntegerValue(0),
std::min(x_slack, range.bounding_area.x_min -
x_.LevelZeroStartMax(range.box_index)));
if (x_min_slack > 0) {
range.bounding_area.x_min -= x_min_slack;
slack -= x_min_slack * overlap.SizeY();
x_slack -= x_min_slack;
}
const IntegerValue x_max_slack =
std::min(x_slack, x_.LevelZeroStartMax(range.box_index) -
range.bounding_area.x_max);
if (x_max_slack > 0) {
range.bounding_area.x_max += x_max_slack;
slack -= x_max_slack * overlap.SizeY();
}
}
}
if (conflict.opp_problem_size == 0) {
CHECK_GT(used_energy, available_energy);
CHECK_GT(available_energy, 0);
}
return ranges_for_explanation;
}
std::vector<RectangleInRange>
NonOverlappingRectanglesEnergyPropagator::GetEnergyConflictForRectangle(
const Rectangle& rectangle,
const std::vector<RectangleInRange>& active_box_ranges) {
struct OverlapPerBox {
IntegerValue energy;
IntegerValue overlap_x_size;
IntegerValue overlap_y_size;
int index;
};
// First select the minimum amount of elements that would be enough to exceed
// the energy.
std::vector<OverlapPerBox> energy_per_box(active_box_ranges.size());
for (int i = 0; i < active_box_ranges.size(); ++i) {
const Rectangle overlap =
active_box_ranges[i].GetMinimumIntersection(rectangle);
OverlapPerBox& box = energy_per_box[i];
box.overlap_x_size = overlap.SizeX();
box.overlap_y_size = overlap.SizeY();
box.index = i;
box.energy = box.overlap_x_size * box.overlap_y_size;
}
std::sort(energy_per_box.begin(), energy_per_box.end(),
[](const OverlapPerBox& a, const OverlapPerBox& b) {
return a.energy > b.energy;
});
const IntegerValue available_energy = rectangle.Area();
IntegerValue used_energy = 0;
std::vector<RectangleInRange> ranges_for_explanation;
ranges_for_explanation.reserve(energy_per_box.size());
for (int i = 0; i < energy_per_box.size(); ++i) {
const OverlapPerBox& box = energy_per_box[i];
used_energy += box.energy;
CHECK_GT(box.energy, 0);
ranges_for_explanation.push_back(active_box_ranges[box.index]);
if (used_energy > available_energy) {
// Use the potential slack to have a stronger reason.
int used_slack_count = 0;
const auto& items = relaxed_result.GetItemsParticipatingOnConflict();
for (int i = 0; i < items.size(); ++i) {
if (!relaxed_result.HasSlack()) {
break;
}
const RectangleInRange& range = conflict.items_for_opp[items[i].index];
const RectangleInRange item_in_zero_level_range = {
.bounding_area = {.x_min = x_.LevelZeroStartMin(range.box_index),
.x_max = x_.LevelZeroStartMax(range.box_index) +
range.x_size,
.y_min = y_.LevelZeroStartMin(range.box_index),
.y_max = y_.LevelZeroStartMax(range.box_index) +
range.y_size},
.x_size = range.x_size,
.y_size = range.y_size};
// There is no point trying to intersect less the item with the rectangle
// than it would on zero-level. So do not waste the slack with it.
const Rectangle max_overlap =
item_in_zero_level_range.GetMinimumIntersection(rectangle);
used_slack_count += relaxed_result.TryUseSlackToReduceItemSize(
i, OrthogonalPackingResult::Coord::kCoordX, max_overlap.SizeX());
used_slack_count += relaxed_result.TryUseSlackToReduceItemSize(
i, OrthogonalPackingResult::Coord::kCoordY, max_overlap.SizeY());
}
num_conflicts_with_slack_ += (used_slack_count > 0);
VLOG_EVERY_N_SEC(2, 3)
<< "Found a conflict on the OPP sub-problem of rectangle: " << rectangle
<< " using "
<< conflict.opp_result.GetItemsParticipatingOnConflict().size() << "/"
<< conflict.items_for_opp.size() << " items.";
std::vector<RectangleInRange> ranges_for_explanation;
ranges_for_explanation.reserve(conflict.items_for_opp.size());
for (const auto& item : relaxed_result.GetItemsParticipatingOnConflict()) {
const RectangleInRange& range = conflict.items_for_opp[item.index];
ranges_for_explanation.push_back(
RectangleInRange::BiggestWithMinIntersection(rectangle, range,
item.size_x, item.size_y));
}
return ranges_for_explanation;
}
@@ -596,34 +517,6 @@ int NonOverlappingRectanglesEnergyPropagator::RegisterWith(
return id;
}
void NonOverlappingRectanglesEnergyPropagator::CheckPropagationIsValid(
const std::vector<RectangleInRange>& ranges,
const Rectangle& rectangle_with_too_much_energy) {
const IntegerValue available_energy = rectangle_with_too_much_energy.Area();
IntegerValue used_energy = 0;
for (const auto& range : ranges) {
const int b = range.box_index;
// Check that we are explaining intervals at least as large as the current.
CHECK_GE(x_.StartMin(b), range.bounding_area.x_min);
CHECK_GE(range.bounding_area.x_max - range.x_size, x_.StartMax(b));
CHECK_GE(y_.StartMin(b), range.bounding_area.y_min);
CHECK_GE(range.bounding_area.y_max - range.y_size, y_.StartMax(b));
// Each one of the boxes-in-range that we found on the cut does intersect
// the rectangle we found.
const auto intersection =
range.GetMinimumIntersection(rectangle_with_too_much_energy);
CHECK_GT(intersection.Area(), 0);
// It cannot intersect more than the size of the object.
CHECK_GE(x_.SizeMin(b), intersection.SizeX());
CHECK_GE(y_.SizeMin(b), intersection.SizeY());
used_energy += intersection.Area();
}
// We have more area inside the rectangle than the area of the rectangle.
CHECK_GT(used_energy, available_energy);
}
bool NonOverlappingRectanglesEnergyPropagator::BuildAndReportEnergyTooLarge(
const std::vector<RectangleInRange>& ranges) {
if (ranges.size() == 2) {

17
ortools/sat/diffn.h
View File

@@ -45,7 +45,7 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface {
y_(*y),
random_(model->GetOrCreate<ModelRandomGenerator>()),
shared_stats_(model->GetOrCreate<SharedStatistics>()),
orthogonal_packing_checker_(shared_stats_) {}
orthogonal_packing_checker_(*random_, shared_stats_) {}
~NonOverlappingRectanglesEnergyPropagator() override;
@@ -54,11 +54,10 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface {
private:
struct Conflict {
std::vector<RectangleInRange> items;
// The Orthogonal Packing subproblem we used.
std::vector<RectangleInRange> items_for_opp;
Rectangle rectangle_with_too_much_energy;
// Size of the Orthogonal Packing Problem we used. Can be smaller than
// `items.size()` if we found a OPP conflict that didn't need all items.
int opp_problem_size;
OrthogonalPackingResult opp_result;
};
std::optional<Conflict> FindConflict(
std::vector<RectangleInRange> active_box_ranges);
@@ -67,11 +66,6 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface {
bool BuildAndReportEnergyTooLarge(
const std::vector<RectangleInRange>& ranges);
void CheckPropagationIsValid(const std::vector<RectangleInRange>& ranges,
const Rectangle& rectangle_with_too_much_energy);
std::vector<RectangleInRange> GetEnergyConflictForRectangle(
const Rectangle& rectangle,
const std::vector<RectangleInRange>& active_box_ranges);
SchedulingConstraintHelper& x_;
SchedulingConstraintHelper& y_;
@@ -81,10 +75,9 @@ class NonOverlappingRectanglesEnergyPropagator : public PropagatorInterface {
int64_t num_calls_ = 0;
int64_t num_conflicts_ = 0;
int64_t num_multiple_conflicts_ = 0;
int64_t num_conflicts_two_boxes_ = 0;
int64_t num_refined_conflicts_ = 0;
int64_t num_orthogonal_packing_conflict_ = 0;
int64_t num_conflicts_with_slack_ = 0;
NonOverlappingRectanglesEnergyPropagator(
const NonOverlappingRectanglesEnergyPropagator&) = delete;

30
ortools/sat/integer_expr.h
View File

@@ -455,10 +455,19 @@ inline std::function<void(Model*)> WeightedSumLowerOrEqual(
}
if (params.new_linear_propagation()) {
model->GetOrCreate<LinearPropagator>()->AddConstraint(
const bool ok = model->GetOrCreate<LinearPropagator>()->AddConstraint(
{}, vars,
std::vector<IntegerValue>(coefficients.begin(), coefficients.end()),
IntegerValue(upper_bound));
if (!ok) {
auto* sat_solver = model->GetOrCreate<SatSolver>();
if (sat_solver->CurrentDecisionLevel() == 0) {
sat_solver->NotifyThatModelIsUnsat();
} else {
LOG(FATAL) << "We currently do not support adding conflicting "
"constraint at positive level.";
}
}
} else {
IntegerSumLE* constraint = new IntegerSumLE(
{}, vars,
@@ -544,8 +553,8 @@ inline std::function<void(Model*)> ConditionalWeightedSumLowerOrEqual(
for (int i = 0; i < vars.size(); ++i) {
expression_min +=
coefficients[i] * (coefficients[i] >= 0
? integer_trail->LowerBound(vars[i])
: integer_trail->UpperBound(vars[i]));
? integer_trail->LevelZeroLowerBound(vars[i])
: integer_trail->LevelZeroUpperBound(vars[i]));
}
if (expression_min == upper_bound) {
// Tricky: as we create integer literal, we might propagate stuff and
@@ -554,12 +563,12 @@ inline std::function<void(Model*)> ConditionalWeightedSumLowerOrEqual(
IntegerValue non_cached_min;
for (int i = 0; i < vars.size(); ++i) {
if (coefficients[i] > 0) {
const IntegerValue lb = integer_trail->LowerBound(vars[i]);
const IntegerValue lb = integer_trail->LevelZeroLowerBound(vars[i]);
non_cached_min += coefficients[i] * lb;
model->Add(Implication(enforcement_literals,
IntegerLiteral::LowerOrEqual(vars[i], lb)));
} else if (coefficients[i] < 0) {
const IntegerValue ub = integer_trail->UpperBound(vars[i]);
const IntegerValue ub = integer_trail->LevelZeroUpperBound(vars[i]);
non_cached_min += coefficients[i] * ub;
model->Add(Implication(enforcement_literals,
IntegerLiteral::GreaterOrEqual(vars[i], ub)));
@@ -574,10 +583,19 @@ inline std::function<void(Model*)> ConditionalWeightedSumLowerOrEqual(
}
} else {
if (params.new_linear_propagation()) {
model->GetOrCreate<LinearPropagator>()->AddConstraint(
const bool ok = model->GetOrCreate<LinearPropagator>()->AddConstraint(
enforcement_literals, vars,
std::vector<IntegerValue>(coefficients.begin(), coefficients.end()),
IntegerValue(upper_bound));
if (!ok) {
auto* sat_solver = model->GetOrCreate<SatSolver>();
if (sat_solver->CurrentDecisionLevel() == 0) {
sat_solver->NotifyThatModelIsUnsat();
} else {
LOG(FATAL) << "We currently do not support adding conflicting "
"constraint at positive level.";
}
}
} else {
IntegerSumLE* constraint = new IntegerSumLE(
enforcement_literals, vars,

10
ortools/sat/linear_constraint_manager.cc
View File

@@ -447,8 +447,6 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
{CeilRatio(threshold, IntegerValue(2)), threshold - min_magnitude,
std::min(threshold_lb, threshold_ub)});
if (max_magnitude > second_threshold) {
term_changed = true;
++num_coeff_strenghtening_;
const int num_terms = ct->num_terms;
for (int i = 0; i < num_terms; ++i) {
// In all cases, we reason on a transformed constraint where the term
@@ -460,18 +458,26 @@ bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var);
const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var);
if (coeff > threshold) {
term_changed = true;
++num_coeff_strenghtening_;
ct->coeffs[i] = threshold;
ct->ub -= (coeff - threshold) * ub;
ct->lb -= (coeff - threshold) * lb;
} else if (coeff > second_threshold && coeff < threshold) {
term_changed = true;
++num_coeff_strenghtening_;
ct->coeffs[i] = second_threshold;
ct->ub -= (coeff - second_threshold) * ub;
ct->lb -= (coeff - second_threshold) * lb;
} else if (coeff < -threshold) {
term_changed = true;
++num_coeff_strenghtening_;
ct->coeffs[i] = -threshold;
ct->ub -= (coeff + threshold) * lb;
ct->lb -= (coeff + threshold) * ub;
} else if (coeff < -second_threshold && coeff > -threshold) {
term_changed = true;
++num_coeff_strenghtening_;
ct->coeffs[i] = -second_threshold;
ct->ub -= (coeff + second_threshold) * lb;
ct->lb -= (coeff + second_threshold) * ub;

157
ortools/sat/linear_propagation.cc
View File

@@ -46,12 +46,24 @@ namespace sat {
void CustomFifoQueue::IncreaseSize(int n) {
CHECK_GE(n, pos_.size());
CHECK_LE(left_, right_);
pos_.resize(n, -1);
tmp_positions_.resize(n, 0);
// We need 1 more space since we can add at most n element.
if (n + 1 < queue_.size()) return;
const int old_size = queue_.size();
queue_.resize(n + 1, 0);
queue_.resize(queue_.capacity(), 0);
// We need to reconstruct the queue in this case.
if (left_ > right_) {
const int diff = queue_.size() - old_size;
for (int i = old_size - 1; i >= left_; --i) {
pos_[queue_[i]] = i + diff;
queue_[i + diff] = queue_[i];
}
left_ += diff;
}
}
int CustomFifoQueue::Pop() {
@@ -230,10 +242,11 @@ void EnforcementPropagator::Untrail(const Trail& /*trail*/, int trail_index) {
EnforcementId EnforcementPropagator::Register(
absl::Span<const Literal> enforcement,
std::function<void(EnforcementStatus)> callback) {
CHECK_EQ(trail_.CurrentDecisionLevel(), 0);
int num_true = 0;
int num_false = 0;
bool is_always_false = false;
temp_literals_.clear();
const int level = trail_.CurrentDecisionLevel();
for (const Literal l : enforcement) {
// Make sure we always have enough room for the literal and its negation.
const int size = std::max(l.Index().value(), l.NegatedIndex().value()) + 1;
@@ -241,23 +254,25 @@ EnforcementId EnforcementPropagator::Register(
watcher_.resize(size);
}
if (assignment_.LiteralIsTrue(l)) {
if (level == 0 || trail_.Info(l.Variable()).level == 0) continue;
++num_true;
continue;
}
if (assignment_.LiteralIsFalse(l)) {
} else if (assignment_.LiteralIsFalse(l)) {
if (level == 0 || trail_.Info(l.Variable()).level == 0) {
is_always_false = true;
break;
}
++num_false;
continue;
}
temp_literals_.push_back(l);
}
gtl::STLSortAndRemoveDuplicates(&temp_literals_);
// Return special indices if never/always enforced.
if (num_false > 0) {
if (is_always_false) {
if (callback != nullptr) callback(EnforcementStatus::IS_FALSE);
return EnforcementId(-1);
}
if (num_true == enforcement.size()) {
if (temp_literals_.empty()) {
if (callback != nullptr) callback(EnforcementStatus::IS_ENFORCED);
return EnforcementId(-1);
}
@@ -267,15 +282,56 @@ EnforcementId EnforcementPropagator::Register(
CHECK(!temp_literals_.empty());
buffer_.insert(buffer_.end(), temp_literals_.begin(), temp_literals_.end());
starts_.push_back(buffer_.size()); // Sentinel.
statuses_.push_back(EnforcementStatus::CANNOT_PROPAGATE);
starts_.push_back(buffer_.size()); // Sentinel/next-start.
// The default status at level zero.
statuses_.push_back(temp_literals_.size() == 1
? EnforcementStatus::CAN_PROPAGATE
: EnforcementStatus::CANNOT_PROPAGATE);
if (temp_literals_.size() == 1) {
watcher_[temp_literals_[0].Index()].push_back(id);
ChangeStatus(id, EnforcementStatus::CAN_PROPAGATE);
} else {
watcher_[temp_literals_[0].Index()].push_back(id);
watcher_[temp_literals_[1].Index()].push_back(id);
// Make sure we watch correct literals.
const auto span = GetSpan(id);
int num_not_true = 0;
for (int i = 0; i < span.size(); ++i) {
if (assignment_.LiteralIsTrue(span[i])) continue;
std::swap(span[num_not_true], span[i]);
++num_not_true;
if (num_not_true == 2) break;
}
// We need to watch one of the literals at highest level.
if (num_not_true == 1) {
int max_level = trail_.Info(span[1].Variable()).level;
for (int i = 2; i < span.size(); ++i) {
const int level = trail_.Info(span[i].Variable()).level;
if (level > max_level) {
max_level = level;
std::swap(span[1], span[i]);
}
}
}
watcher_[span[0].Index()].push_back(id);
watcher_[span[1].Index()].push_back(id);
}
// Change status, call callback and set up untrail if the status is different
// from EnforcementStatus::CANNOT_PROPAGATE.
if (num_false > 0) {
ChangeStatus(id, EnforcementStatus::IS_FALSE);
} else if (num_true == temp_literals_.size()) {
ChangeStatus(id, EnforcementStatus::IS_ENFORCED);
} else if (num_true + 1 == temp_literals_.size()) {
ChangeStatus(id, EnforcementStatus::CAN_PROPAGATE);
// Because this is the default status, we still need to call the callback.
if (temp_literals_.size() == 1) {
if (callbacks_[id] != nullptr) {
callbacks_[id](EnforcementStatus::CAN_PROPAGATE);
}
}
}
return id;
}
@@ -390,6 +446,22 @@ void EnforcementPropagator::ChangeStatus(EnforcementId id,
if (callbacks_[id] != nullptr) callbacks_[id](new_status);
}
EnforcementStatus EnforcementPropagator::DebugStatus(EnforcementId id) {
if (id < 0) return EnforcementStatus::IS_ENFORCED;
int num_true = 0;
for (const Literal l : GetSpan(id)) {
if (assignment_.LiteralIsFalse(l)) {
return EnforcementStatus::IS_FALSE;
}
if (assignment_.LiteralIsTrue(l)) ++num_true;
}
const int size = GetSpan(id).size();
if (num_true == size) return EnforcementStatus::IS_ENFORCED;
if (num_true + 1 == size) return EnforcementStatus::CAN_PROPAGATE;
return EnforcementStatus::CANNOT_PROPAGATE;
}
LinearPropagator::LinearPropagator(Model* model)
: trail_(model->GetOrCreate<Trail>()),
integer_trail_(model->GetOrCreate<IntegerTrail>()),
@@ -438,13 +510,14 @@ void LinearPropagator::SetLevel(int level) {
while (!propagation_queue_.empty()) {
in_queue_[propagation_queue_.Pop()] = false;
}
for (int i = rev_at_false_size_; i < in_queue_and_at_false_.size(); ++i) {
in_queue_[in_queue_and_at_false_[i]] = false;
for (int i = rev_unenforced_size_; i < unenforced_constraints_.size();
++i) {
in_queue_[unenforced_constraints_[i]] = false;
}
in_queue_and_at_false_.resize(rev_at_false_size_);
unenforced_constraints_.resize(rev_unenforced_size_);
} else if (level > previous_level_) {
rev_at_false_size_ = in_queue_and_at_false_.size();
rev_int_repository_->SaveState(&rev_at_false_size_);
rev_unenforced_size_ = unenforced_constraints_.size();
rev_int_repository_->SaveState(&rev_unenforced_size_);
}
previous_level_ = level;
@@ -468,12 +541,12 @@ bool LinearPropagator::Propagate() {
AddWatchedToQueue(var);
}
// TODO(user): Abort this propagator as soon as a Boolean is propagated ? so
// that we always finish the Boolean propagation first. This can happen when
// we push a bound that has associated Booleans. The idea is to resume from
// our current state when we are called again. Note however that we have to
// clear the propagated_by_ info has other propagator might have pushed the
// same variable further.
// We abort this propagator as soon as a Boolean is propagated, so that we
// always finish the Boolean propagation first. This can happen when we push a
// bound that has associated Booleans or push enforcement to false. The idea
// is to resume from our current state when we are called again. Note however
// that we have to clear the propagated_by_ info (as done above) has other
// propagator might have pushed the same variable further.
//
// Empty FIFO queue.
const int saved_index = trail_->Index();
@@ -497,13 +570,15 @@ bool LinearPropagator::Propagate() {
}
// Adds a new constraint to the propagator.
void LinearPropagator::AddConstraint(
bool LinearPropagator::AddConstraint(
absl::Span<const Literal> enforcement_literals,
absl::Span<const IntegerVariable> vars,
absl::Span<const IntegerValue> coeffs, IntegerValue upper_bound) {
if (vars.empty()) return;
for (const Literal l : enforcement_literals) {
if (trail_->Assignment().LiteralIsFalse(l)) return;
if (vars.empty()) return true;
if (trail_->CurrentDecisionLevel() == 0) {
for (const Literal l : enforcement_literals) {
if (trail_->Assignment().LiteralIsFalse(l)) return true;
}
}
// Make sure max_variations_ is of correct size.
@@ -551,12 +626,13 @@ void LinearPropagator::AddConstraint(
propagation_queue_.IncreaseSize(in_queue_.size());
if (!enforcement_literals.empty()) {
infos_.back().enf_status = EnforcementStatus::CANNOT_PROPAGATE;
infos_.back().enf_status =
static_cast<int>(EnforcementStatus::CANNOT_PROPAGATE);
infos_.back().enf_id = enforcement_propagator_->Register(
enforcement_literals, [this, id](EnforcementStatus status) {
infos_[id].enf_status = status;
infos_[id].enf_status = static_cast<int>(status);
// TODO(user): With some care, when we cannot propagate or the
// constraint is not enforced, we could live in_queue_[] at true but
// constraint is not enforced, we could leave in_queue_[] at true but
// not put the constraint in the queue.
if (status == EnforcementStatus::CAN_PROPAGATE ||
status == EnforcementStatus::IS_ENFORCED) {
@@ -567,7 +643,7 @@ void LinearPropagator::AddConstraint(
} else {
AddToQueueIfNeeded(id);
infos_.back().enf_id = -1;
infos_.back().enf_status = EnforcementStatus::IS_ENFORCED;
infos_.back().enf_status = static_cast<int>(EnforcementStatus::IS_ENFORCED);
}
for (const IntegerVariable var : GetVariables(infos_[id])) {
@@ -593,6 +669,9 @@ void LinearPropagator::AddConstraint(
watcher_->WatchLowerBound(var, watcher_id_);
}
}
// Propagate this new constraint.
return PropagateOneConstraint(id);
}
absl::Span<IntegerValue> LinearPropagator::GetCoeffs(
@@ -635,14 +714,16 @@ bool LinearPropagator::PropagateOneConstraint(int id) {
// Skip constraint not enforced or that cannot propagate if false.
ConstraintInfo& info = infos_[id];
if (info.enf_status == EnforcementStatus::IS_FALSE ||
info.enf_status == EnforcementStatus::CANNOT_PROPAGATE) {
const EnforcementStatus enf_status = EnforcementStatus(info.enf_status);
DCHECK_EQ(enf_status, enforcement_propagator_->DebugStatus(info.enf_id));
if (enf_status == EnforcementStatus::IS_FALSE ||
enf_status == EnforcementStatus::CANNOT_PROPAGATE) {
DCHECK(!in_queue_[id]);
if (info.enf_status == EnforcementStatus::IS_FALSE) {
if (enf_status == EnforcementStatus::IS_FALSE) {
// We mark this constraint as in the queue but will never inspect it
// again until we backtrack over this time.
in_queue_[id] = true;
in_queue_and_at_false_.push_back(id);
unenforced_constraints_.push_back(id);
}
++num_ignored_;
return true;
@@ -706,7 +787,9 @@ bool LinearPropagator::PropagateOneConstraint(int id) {
}
// We can only propagate more if all the enforcement literals are true.
if (info.enf_status != EnforcementStatus::IS_ENFORCED) return true;
if (info.enf_status != static_cast<int>(EnforcementStatus::IS_ENFORCED)) {
return true;
}
// The lower bound of all the variables except one can be used to update the
// upper bound of the last one.

24
ortools/sat/linear_propagation.h
View File

@@ -48,8 +48,6 @@ class CustomFifoQueue {
public:
CustomFifoQueue() = default;
// Note that this requires the queue to be empty or to never have been poped
// before.
void IncreaseSize(int n);
int Pop();
@@ -83,7 +81,7 @@ class CustomFifoQueue {
// An enforced constraint can be in one of these 4 states.
// Note that we rely on the integer encoding to take 2 bits for optimization.
enum EnforcementStatus {
enum class EnforcementStatus {
// One enforcement literal is false.
IS_FALSE = 0,
// More than two literals are unassigned.
@@ -128,6 +126,10 @@ class EnforcementPropagator : SatPropagator {
EnforcementStatus Status(EnforcementId id) const { return statuses_[id]; }
// Recompute the status from the current assignment.
// This should only used in DCHECK().
EnforcementStatus DebugStatus(EnforcementId id);
private:
absl::Span<Literal> GetSpan(EnforcementId id);
absl::Span<const Literal> GetSpan(EnforcementId id) const;
@@ -181,7 +183,10 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface {
void SetLevel(int level) final;
// Adds a new constraint to the propagator.
void AddConstraint(absl::Span<const Literal> enforcement_literals,
// We support adding constraint at a positive level:
// - This will push new propagation right away.
// - This will returns false if the constraint is currently a conflict.
bool AddConstraint(absl::Span<const Literal> enforcement_literals,
absl::Span<const IntegerVariable> vars,
absl::Span<const IntegerValue> coeffs,
IntegerValue upper_bound);
@@ -244,8 +249,9 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface {
// IntegerVariable that changed since the last clear.
SparseBitset<IntegerVariable> modified_vars_;
// Per constraint info used during propagation.
std::vector<ConstraintInfo> infos_;
// Per constraint info used during propagation. Note that we keep pointer for
// the rev_size/rhs there, so we do need a deque.
std::deque<ConstraintInfo> infos_;
// Buffer of the constraints data.
//
@@ -267,8 +273,10 @@ class LinearPropagator : public PropagatorInterface, ReversibleInterface {
std::vector<bool> in_queue_;
CustomFifoQueue propagation_queue_;
int rev_at_false_size_ = 0;
std::vector<int> in_queue_and_at_false_;
// Unenforced constraints are marked as "in_queue_" but not actually added
// to the propagation_queue_.
int rev_unenforced_size_ = 0;
std::vector<int> unenforced_constraints_;
// Watchers.
absl::StrongVector<IntegerVariable, bool> is_watched_;

50
ortools/util/bitset.h
View File

@@ -19,10 +19,14 @@
#include <string.h>
#include <algorithm>
#include <cstddef>
#include <cstdint>
#include <iterator>
#include <string>
#include <tuple>
#include <vector>
#include "absl/log/check.h"
#include "ortools/base/logging.h"
#include "ortools/base/types.h"
@@ -412,6 +416,8 @@ inline uint64_t TwoBitsFromPos64(uint64_t pos) {
template <typename IndexType = int64_t>
class Bitset64 {
public:
using value_type = IndexType;
// When speed matter, caching the base pointer like this to access this class
// in a read only mode help.
class ConstView {
@@ -592,9 +598,21 @@ class Bitset64 {
//
// IMPORTANT: Because the iterator "caches" the current uint64_t bucket, this
// will probably not do what you want if Bitset64 is modified while iterating.
class EndIterator {};
class Iterator {
public:
// Make this iterator a std::forward_iterator, so it works with std::sample,
// std::max_element, etc.
Iterator() : data_(nullptr), size_(0) {}
Iterator(Iterator&& other) = default;
Iterator(const Iterator& other) = default;
Iterator& operator=(const Iterator& other) = default;
using difference_type = std::ptrdiff_t;
using iterator_category = std::forward_iterator_tag;
using value_type = IndexType;
using size_type = std::size_t;
using reference = value_type&;
using pointer = value_type*;
explicit Iterator(const Bitset64& bitset)
: data_(bitset.data_.data()), size_(bitset.data_.size()) {
if (!bitset.data_.empty()) {
@@ -603,17 +621,34 @@ class Bitset64 {
}
}
bool operator!=(const EndIterator&) const { return size_ != 0; }
static Iterator EndIterator(const Bitset64& bitset) {
return Iterator(bitset.data_.data());
}
bool operator==(const Iterator& other) const { return !(*this == other); }
bool operator!=(const Iterator& other) const {
if (other.size_ == 0) {
return size_ != 0;
}
return std::tie(index_, current_) !=
std::tie(other.index_, other.current_);
}
IndexType operator*() const { return IndexType(index_); }
void operator++() {
Iterator operator++(int) {
Iterator other = *this;
(*this)++;
return other;
}
Iterator& operator++() {
int bucket = BitOffset64(index_);
while (current_ == 0) {
bucket++;
if (bucket == size_) {
size_ = 0;
return;
return *this;
}
current_ = data_[bucket];
}
@@ -621,10 +656,13 @@ class Bitset64 {
// Computes the index and clear the least significant bit of current_.
index_ = BitShift64(bucket) | LeastSignificantBitPosition64(current_);
current_ &= current_ - 1;
return *this;
}
private:
const uint64_t* const data_;
explicit Iterator(const uint64_t* data) : data_(data), size_(0) {}
const uint64_t* data_;
int size_;
int index_ = 0;
uint64_t current_ = 0;
@@ -633,7 +671,7 @@ class Bitset64 {
// Allows range-based "for" loop on the non-zero positions:
// for (const IndexType index : bitset) {}
Iterator begin() const { return Iterator(*this); }
EndIterator end() const { return EndIterator(); }
Iterator end() const { return Iterator::EndIterator(*this); }
// Cryptic function! This is just an optimized version of a given piece of
// code and has probably little general use.