more work on 2d packing

This commit is contained in:
Laurent Perron
2024-01-25 13:26:36 +01:00
parent 3f511acb83
commit d9036bffeb
6 changed files with 403 additions and 67 deletions

View File

@@ -15,15 +15,20 @@
#include <algorithm>
#include <cstdint>
#include <limits>
#include <optional>
#include <string>
#include <utility>
#include <vector>
#include "absl/log/check.h"
#include "absl/numeric/bits.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"
namespace operations_research {
namespace sat {
@@ -36,6 +41,10 @@ OrthogonalPackingInfeasibilityDetector::
{"OrthogonalPackingInfeasibilityDetector/called", num_calls_});
stats.push_back(
{"OrthogonalPackingInfeasibilityDetector/conflicts", num_conflicts_});
stats.push_back({"OrthogonalPackingInfeasibilityDetector/dff0_conflicts",
num_conflicts_dff0_});
stats.push_back({"OrthogonalPackingInfeasibilityDetector/dff2_conflicts",
num_conflicts_dff2_});
stats.push_back({"OrthogonalPackingInfeasibilityDetector/trivial_conflicts",
num_trivial_conflicts_});
stats.push_back({"OrthogonalPackingInfeasibilityDetector/conflicts_two_items",
@@ -87,11 +96,9 @@ std::optional<std::pair<int, int>> FindPairwiseConflict(
return std::nullopt;
}
// Check for conflict using the f_0 dual feasible function described on Côté,
// Jean-François, Mohamed Haouari, and Manuel Iori. "A primal decomposition
// algorithm for the two-dimensional bin packing problem." arXiv preprint
// arXiv:1909.06835 (2019). This function tries all possible values of the `k`
// parameter.
// Check for conflict using the f_0 dual feasible function. See for example [1]
// for some discussion. This function tries all possible values of the `k`
// parameter and returns the conflict needing the least number of items.
//
// The `f_0^k(x)` function for a total bin width `C` and a parameter `k` taking
// values in [0, C/2] is defined as:
@@ -100,30 +107,46 @@ std::optional<std::pair<int, int>> FindPairwiseConflict(
// f_0^k(x) = | x, if k <= x <= C - k,
// \ 0, if x < k.
//
// If a conflict is found, this replaces the conflict in `result` if it detects
// a conflict using less items or if `result` is empty.
// If no conflict is found, returns an empty vector.
//
// The current implementation is a bit more general than a simple check using
// f_0 described above. This implementation can take a function g(x) that is
// non-decreasing and satisfy g(0)=0 and it will check for conflict using
// g(f_0^k(x)) for all values of k, but without recomputing g(x) `k` times. This
// is handy if g() is a DFF that is slow to compute. g(x) is described by the
// vector g_x[i] = g(sizes_x[i]) and the variable g_max = g(x_bb_size).
//
// The algorithm is the same if we swap the x and y dimension.
void GetOrReplaceDffConflict(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 total_energy,
IntegerValue bb_area, std::vector<int>* result) {
//
// [1] Côté, Jean-François, Mohamed Haouari, and Manuel Iori. "A primal
// decomposition algorithm for the two-dimensional bin packing problem." arXiv
// preprint https://arxiv.org/abs/1909.06835 (2019).
std::vector<int> GetDffConflict(
absl::Span<const IntegerValue> sizes_x,
absl::Span<const IntegerValue> sizes_y,
absl::Span<const int> index_by_decreasing_x_size,
absl::Span<const IntegerValue> g_x, IntegerValue g_max,
IntegerValue x_bb_size, IntegerValue total_energy, IntegerValue bb_area) {
std::vector<int> result;
// If we found a conflict for a k parameter, which is rare, recompute the
// total used energy consumed by the items to find the minimal set of
// conflicting items.
int num_items = sizes_x.size();
auto add_result_if_better = [&result, &sizes_x, &sizes_y, num_items,
&x_bb_size, &bb_area](const IntegerValue k) {
&x_bb_size, &bb_area, &g_max,
&g_x](const IntegerValue k) {
std::vector<std::pair<int, IntegerValue>> index_to_energy;
index_to_energy.reserve(num_items);
for (int i = 0; i < num_items; i++) {
IntegerValue x_size = sizes_x[i];
if (x_size > x_bb_size - k) {
x_size = x_bb_size;
} else if (x_size < k) {
IntegerValue point_value;
if (sizes_x[i] > x_bb_size - k) {
point_value = g_max;
} else if (sizes_x[i] < k) {
continue;
} else {
point_value = g_x[i];
}
index_to_energy.push_back({i, x_size * sizes_y[i]});
index_to_energy.push_back({i, point_value * sizes_y[i]});
}
std::sort(index_to_energy.begin(), index_to_energy.end(),
[](const std::pair<int, IntegerValue>& a,
@@ -134,10 +157,10 @@ void GetOrReplaceDffConflict(absl::Span<const IntegerValue> sizes_x,
for (int i = 0; i < index_to_energy.size(); i++) {
recomputed_energy += index_to_energy[i].second;
if (recomputed_energy > bb_area) {
if (result->empty() || i + 1 < result->size()) {
result->resize(i + 1);
if (result.empty() || i + 1 < result.size()) {
result.resize(i + 1);
for (int j = 0; j <= i; j++) {
(*result)[j] = index_to_energy[j].first;
result[j] = index_to_energy[j].first;
}
}
break;
@@ -153,6 +176,9 @@ void GetOrReplaceDffConflict(absl::Span<const IntegerValue> sizes_x,
// large items and small ones are not symmetric with respect to what values of
// k are important.
IntegerValue current_energy = total_energy;
if (current_energy > bb_area) {
add_result_if_better(0);
}
// We keep an index on the largest item yet-to-be enlarged and the smallest
// one yet-to-be removed.
int removing_item_index = index_by_decreasing_x_size.size() - 1;
@@ -160,20 +186,20 @@ void GetOrReplaceDffConflict(absl::Span<const IntegerValue> sizes_x,
while (enlarging_item_index < index_by_decreasing_x_size.size()) {
int index = index_by_decreasing_x_size[enlarging_item_index];
IntegerValue size = sizes_x[index];
if (size <= x_bb_size / 2) {
break;
}
// Note that since `size_x` is decreasing, we test increasingly large
// values of k. Also note that a item with size `k` cannot fit alongside
// an item with size `size_x`, but smaller ones can.
const IntegerValue k = x_bb_size - size + 1;
if (2 * k > x_bb_size) {
break;
}
// First, add the area contribution of enlarging the all the items of size
// exactly size_x. All larger items were already enlarged in the previous
// iterations.
do {
index = index_by_decreasing_x_size[enlarging_item_index];
size = sizes_x[index];
current_energy += (x_bb_size - size) * sizes_y[index];
current_energy += (g_max - g_x[index]) * sizes_y[index];
enlarging_item_index++;
} while (enlarging_item_index < index_by_decreasing_x_size.size() &&
sizes_x[index_by_decreasing_x_size[enlarging_item_index]] == size);
@@ -183,7 +209,7 @@ void GetOrReplaceDffConflict(absl::Span<const IntegerValue> sizes_x,
while (removing_item_index >= 0 &&
sizes_x[index_by_decreasing_x_size[removing_item_index]] < k) {
const int remove_idx = index_by_decreasing_x_size[removing_item_index];
current_energy -= sizes_x[remove_idx] * sizes_y[remove_idx];
current_energy -= g_x[remove_idx] * sizes_y[remove_idx];
removing_item_index--;
}
@@ -191,6 +217,128 @@ void GetOrReplaceDffConflict(absl::Span<const IntegerValue> sizes_x,
add_result_if_better(k);
}
}
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;
}
// 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
// least) amount, depending on the sign it appears in the formula.
IntegerValue sum_widths = 0;
for (int i = 0; i < sizes_x.size(); i++) {
const IntegerValue x_size = sizes_x[i];
if (2 * x_size > x_bb_size) {
sum_widths += 2 * sizes_y[i];
} else if (2 * x_size == x_bb_size) {
sum_widths += sizes_y[i];
}
}
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);
// `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);
}
for (int i = 1; i <= sqrt_bb_size; i++) {
const QuickSmallDivision div(i);
if (i > 1) {
candidates.Set(div.DivideByDivisor(x_bb_size.value() + round_up.value()));
}
for (int k = 0; k < sizes_x.size(); k++) {
IntegerValue x_size = sizes_x[k];
if (2 * x_size > x_bb_size && x_size < x_bb_size) {
candidates.Set(
div.DivideByDivisor(x_bb_size.value() - x_size.value() + 1));
} else if (2 * x_size < x_bb_size) {
candidates.Set(div.DivideByDivisor(x_size.value()));
}
}
}
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);
// Apply the nice result described on [1]: if we are testing the DFF
// f_2^k(f_0^l(x)) for all values of `l`, the only values of `k` greater than
// C/4 we need to test are {C/4+1, C/3+1}.
//
// In the same reference there is a proof that this way of composing f_0 and
// f_2 cover all possible ways of composing the two functions, including
// composing several times each.
//
// [1] F. Clautiaux, PhD thesis, hal/tel-00749411.
candidates.Resize(x_bb_size / 4 + 1); // Erase all >= C/4
candidates.Resize(x_bb_size / 3 + 2); // Make room for the two special values
candidates.Set(x_bb_size / 4 + 1);
if (x_bb_size > 3) {
candidates.Set(x_bb_size / 3 + 1);
}
// Finally run our small loop to look for the conflict!
std::vector<IntegerValue> modified_sizes(num_items);
for (const IntegerValue k : candidates) {
const RoundingDualFeasibleFunction dff(x_bb_size, k);
IntegerValue energy = 0;
for (int i = 0; i < num_items; i++) {
modified_sizes[i] = dff(sizes_x[i]);
energy += modified_sizes[i] * sizes_y[i];
}
const IntegerValue modified_x_bb_size = dff(x_bb_size);
auto dff0_res = GetDffConflict(
sizes_x, sizes_y, index_by_decreasing_x_size, modified_sizes,
modified_x_bb_size, x_bb_size, energy, modified_x_bb_size * y_bb_size);
if (!dff0_res.empty() &&
(result.empty() || dff0_res.size() < result.size())) {
result = dff0_res;
}
}
return result;
}
} // namespace
@@ -200,6 +348,13 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
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_++;
const int num_items = sizes_x.size();
@@ -234,49 +389,108 @@ OrthogonalPackingInfeasibilityDetector::TestFeasibility(
[&sizes_y](int a, int b) { return sizes_y[a] > sizes_y[b]; });
// First look for pairwise incompatible pairs.
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}};
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}};
}
if (num_items == 2) {
return Result{.result = Status::FEASIBLE};
}
}
if (num_items == 2) {
return Result{.result = Status::FEASIBLE};
}
// 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 whole box, which
// is obviously incompatible, and this incompatibility would be present
// already before enlarging the items since it is a DFF. So it is enough to
// test making items wide or high, but no need to try both.
std::vector<int> result;
GetOrReplaceDffConflict(sizes_x, sizes_y, index_by_decreasing_x_size_,
bounding_box_size.first, total_energy, bb_area,
&result);
if (result.size() == 3) {
// Since we already checked pairwise conflict, we won't find anything
// narrower than a three-item conflict.
num_conflicts_++;
return Result{.result = Status::INFEASIBLE,
.minimum_unfeasible_subproblem = result};
if (total_energy > bb_area) {
conflict_type = ConflictType::TRIVIAL;
std::vector<std::pair<int, IntegerValue>> index_to_energy;
index_to_energy.reserve(num_items);
for (int i = 0; i < num_items; i++) {
index_to_energy.push_back({i, sizes_x[i] * sizes_y[i]});
}
std::sort(index_to_energy.begin(), index_to_energy.end(),
[](const std::pair<int, IntegerValue>& a,
const std::pair<int, IntegerValue>& b) {
return a.second > b.second;
});
IntegerValue recomputed_energy = 0;
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);
for (int j = 0; j <= i; j++) {
result[j] = index_to_energy[j].first;
}
break;
}
}
}
GetOrReplaceDffConflict(sizes_y, sizes_x, index_by_decreasing_y_size_,
bounding_box_size.second, total_energy, bb_area,
&result);
// The items have a trivial conflict and we didn't found any tighter one.
if (total_energy > bb_area && result.size() == num_items) {
num_trivial_conflicts_++;
auto set_conflict_if_better = [&result, &conflict_type](
const std::vector<int>& conflict,
ConflictType type) {
if (!conflict.empty() &&
(result.empty() || conflict.size() < result.size())) {
conflict_type = type;
result = conflict;
}
};
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
// whole box, which is obviously incompatible, and this incompatibility
// would be present already before enlarging the items since it is a DFF. So
// it is enough to test making items wide or high, but no need to try both.
auto conflict =
GetDffConflict(sizes_x, sizes_y, index_by_decreasing_x_size_, sizes_x,
bounding_box_size.first, bounding_box_size.first,
total_energy, bb_area);
set_conflict_if_better(conflict, ConflictType::DFF_F0);
conflict = GetDffConflict(sizes_y, sizes_x, index_by_decreasing_y_size_,
sizes_y, bounding_box_size.second,
bounding_box_size.second, total_energy, bb_area);
set_conflict_if_better(conflict, ConflictType::DFF_F0);
}
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);
set_conflict_if_better(conflict, ConflictType::DFF_F2);
conflict = CheckFeasibilityWithDualFunction2(
sizes_y, sizes_x, index_by_decreasing_y_size_, bounding_box_size.second,
bounding_box_size.first);
set_conflict_if_better(conflict, ConflictType::DFF_F2);
}
if (!result.empty()) {
num_conflicts_++;
switch (conflict_type) {
case ConflictType::DFF_F0:
num_conflicts_dff0_++;
break;
case ConflictType::DFF_F2:
num_conflicts_dff2_++;
break;
case ConflictType::TRIVIAL:
// The total area of the items was larger than the area of the box.
num_trivial_conflicts_++;
break;
case ConflictType::NO_CONFLICT:
LOG(FATAL) << "Should never happen";
break;
}
return Result{.result = Status::INFEASIBLE,
.minimum_unfeasible_subproblem = result};
}

View File

@@ -15,9 +15,11 @@
#define OR_TOOLS_SAT_2D_ORTHOGONAL_PACKING_H_
#include <cstdint>
#include <limits>
#include <utility>
#include <vector>
#include "absl/log/check.h"
#include "absl/types/span.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/synchronization.h"
@@ -25,13 +27,20 @@
namespace operations_research {
namespace sat {
struct OrthogonalPackingOptions {
bool use_pairwise = true;
bool use_dff_f0 = true;
bool use_dff_f2 = true;
};
// Class for solving the orthogonal packing problem when it can be done
// efficiently (ie., not applying any heuristic slower than O(N^2)).
class OrthogonalPackingInfeasibilityDetector {
public:
explicit OrthogonalPackingInfeasibilityDetector(
SharedStatistics* shared_stats)
: shared_stats_(shared_stats) {}
SharedStatistics* shared_stats,
const OrthogonalPackingOptions& options = OrthogonalPackingOptions())
: options_(options), shared_stats_(shared_stats) {}
~OrthogonalPackingInfeasibilityDetector();
@@ -51,6 +60,7 @@ class OrthogonalPackingInfeasibilityDetector {
std::pair<IntegerValue, IntegerValue> bounding_box_size);
private:
const OrthogonalPackingOptions options_;
std::vector<int> index_by_decreasing_x_size_;
std::vector<int> index_by_decreasing_y_size_;
@@ -58,10 +68,93 @@ class OrthogonalPackingInfeasibilityDetector {
int64_t num_conflicts_ = 0;
int64_t num_conflicts_two_items_ = 0;
int64_t num_trivial_conflicts_ = 0;
int64_t num_conflicts_dff2_ = 0;
int64_t num_conflicts_dff0_ = 0;
SharedStatistics* shared_stats_;
};
// Dual Feasible Function based on rounding. Called f_2 on [1].
//
// The `f_2^k(x)` function for an integer x in [0, C] and a parameter `k` taking
// values in [0, C/2] is defined as:
//
// / 2 * [ floor(C / k) - floor( (C - x) / k) ], if x > C / 2,
// f_2^k(x) = | floor(C / k), if k = C / 2,
// \ floor(x / k), if x < C / 2.
//
// This function is a Maximal Dual Feasible Function. Ie., it satisfies:
// - f_2 is nondecreasing,
// - f_2 is superadditive, i.e., f_2(x) + f_2(y) <= f_2(x + y),
// - f_2 is symmetric, i.e., f_2(x) + f_2(C - x) = f_2(C),
// - f_2(0) = 0.
//
// [1] Carlier, Jacques, François Clautiaux, and Aziz Moukrim. "New reduction
// procedures and lower bounds for the two-dimensional bin packing problem with
// fixed orientation." Computers & Operations Research 34.8 (2007): 2223-2250.
class RoundingDualFeasibleFunction {
public:
// `max_x` must fit in a uint16_t and `k` in [0, max_x/2].
RoundingDualFeasibleFunction(IntegerValue max_x, IntegerValue k)
: div_(k.value()),
max_x_(max_x),
c_k_(div_.DivideByDivisor(max_x_.value())) {
DCHECK_GT(k, 0);
DCHECK_LE(2 * k, max_x_);
DCHECK_LE(max_x_, std::numeric_limits<uint16_t>::max());
}
// `x` must be in [0, max_x].
IntegerValue operator()(IntegerValue x) const {
DCHECK_GE(x, 0);
DCHECK_LE(x, max_x_);
if (2 * x > max_x_) {
return 2 * (c_k_ - div_.DivideByDivisor(max_x_.value() - x.value()));
} else if (2 * x == max_x_) {
return c_k_;
} else {
return 2 * div_.DivideByDivisor(x.value());
}
}
private:
const QuickSmallDivision div_;
const IntegerValue max_x_;
const IntegerValue c_k_;
};
// Same as above for k = 2^log2_k.
class RoundingDualFeasibleFunctionPowerOfTwo {
public:
RoundingDualFeasibleFunctionPowerOfTwo(IntegerValue max_x,
IntegerValue log2_k)
: log2_k_(log2_k), max_x_(max_x), c_k_(max_x_ >> log2_k_) {
DCHECK_GE(log2_k_, 0);
DCHECK_LT(log2_k_, 63);
DCHECK_LE(2 * (1 << log2_k), max_x_);
DCHECK_LE(max_x_, std::numeric_limits<int64_t>::max() / 2);
}
IntegerValue operator()(IntegerValue x) const {
DCHECK_GE(x, 0);
DCHECK_LE(x, max_x_);
if (2 * x > max_x_) {
return 2 * (c_k_ - ((max_x_ - x) >> log2_k_));
} else if (2 * x == max_x_) {
return c_k_;
} else {
return 2 * (x >> log2_k_);
}
}
private:
const IntegerValue log2_k_;
const IntegerValue max_x_;
const IntegerValue c_k_;
};
} // namespace sat
} // namespace operations_research

View File

@@ -895,6 +895,7 @@ cc_library(
"//ortools/algorithms:sparse_permutation",
"//ortools/base",
"@com_google_absl//absl/log:check",
"@com_google_absl//absl/types:span",
],
)
@@ -1817,7 +1818,11 @@ cc_library(
deps = [
":integer",
":synchronization",
":util",
"//ortools/util:bitset",
"@com_google_absl//absl/log",
"@com_google_absl//absl/log:check",
"@com_google_absl//absl/numeric:bits",
"@com_google_absl//absl/types:span",
],
)

View File

@@ -122,6 +122,28 @@ inline bool AtMinOrMaxInt64I(IntegerValue t) {
return AtMinOrMaxInt64(t.value());
}
// Helper for dividing several small integers by the same value. Note that there
// is no point using this class is the divisor is a compile-time constant, since
// the compiler should be smart enough to do this automatically.
// Building a `QuickSmallDivision` object costs an integer division, but each
// call to `DivideByDivisor` will only do an integer multiplication and a shift.
//
// This class always return the exact value of the division for all possible
// values of `dividend` and `divisor`.
class QuickSmallDivision {
public:
explicit QuickSmallDivision(uint16_t divisor)
: inverse_((1ull << 48) / divisor + 1) {}
uint16_t DivideByDivisor(uint16_t dividend) const {
return static_cast<uint16_t>((inverse_ * static_cast<uint64_t>(dividend)) >>
48);
}
private:
uint64_t inverse_;
};
// Returns dividend - FloorRatio(dividend, divisor) * divisor;
//
// This function is around the same speed than the computation above, but it

View File

@@ -19,6 +19,7 @@
#include <vector>
#include "absl/log/check.h"
#include "absl/types/span.h"
#include "ortools/algorithms/dynamic_partition.h"
#include "ortools/algorithms/sparse_permutation.h"
#include "ortools/base/logging.h"
@@ -27,7 +28,7 @@ namespace operations_research {
namespace sat {
std::vector<std::vector<int>> BasicOrbitopeExtraction(
const std::vector<std::unique_ptr<SparsePermutation>>& generators) {
absl::Span<const std::unique_ptr<SparsePermutation>> generators) {
// Count the number of permutations that are compositions of 2-cycle and
// regroup them according to the number of cycles.
std::vector<std::vector<int>> num_cycles_to_2cyclers;
@@ -150,7 +151,7 @@ std::vector<std::vector<int>> BasicOrbitopeExtraction(
}
std::vector<int> GetOrbits(
int n, const std::vector<std::unique_ptr<SparsePermutation>>& generators) {
int n, absl::Span<const std::unique_ptr<SparsePermutation>> generators) {
MergingPartition union_find;
union_find.Reset(n);
for (const std::unique_ptr<SparsePermutation>& perm : generators) {

View File

@@ -17,6 +17,7 @@
#include <memory>
#include <vector>
#include "absl/types/span.h"
#include "ortools/algorithms/sparse_permutation.h"
namespace operations_research {
@@ -45,7 +46,7 @@ namespace sat {
// graph20-20-1rand.mps.gz. I suspect the generators provided by the detection
// code follow our preconditions.
std::vector<std::vector<int>> BasicOrbitopeExtraction(
const std::vector<std::unique_ptr<SparsePermutation>>& generators);
absl::Span<const std::unique_ptr<SparsePermutation>> generators);
// Returns a vector of size n such that
// - orbits[i] == -1 iff i is never touched by the generators (singleton orbit).
@@ -53,7 +54,7 @@ std::vector<std::vector<int>> BasicOrbitopeExtraction(
//
// TODO(user): We could reuse the internal memory if needed.
std::vector<int> GetOrbits(
int n, const std::vector<std::unique_ptr<SparsePermutation>>& generators);
int n, absl::Span<const std::unique_ptr<SparsePermutation>> generators);
// Returns the orbits under the given orbitope action.
// Same results format as in GetOrbits(). Note that here, the orbit index