Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
2d_orthogonal_packing.cc
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
15
16#include <algorithm>
17#include <cstdint>
18#include <limits>
19#include <optional>
20#include <string>
21#include <tuple>
22#include <utility>
23#include <vector>
24
25#include "absl/log/check.h"
26#include "absl/log/log.h"
27#include "absl/log/vlog_is_on.h"
28#include "absl/numeric/bits.h"
29#include "absl/random/distributions.h"
30#include "absl/types/span.h"
34#include "ortools/sat/util.h"
35#include "ortools/util/bitset.h"
36
37namespace operations_research {
38namespace sat {
39
42 if (!VLOG_IS_ON(1)) return;
43 std::vector<std::pair<std::string, int64_t>> stats;
44 stats.push_back(
45 {"OrthogonalPackingInfeasibilityDetector/called", num_calls_});
46 stats.push_back(
47 {"OrthogonalPackingInfeasibilityDetector/conflicts", num_conflicts_});
48 stats.push_back({"OrthogonalPackingInfeasibilityDetector/dff0_conflicts",
49 num_conflicts_dff0_});
50 stats.push_back({"OrthogonalPackingInfeasibilityDetector/dff2_conflicts",
51 num_conflicts_dff2_});
52 stats.push_back({"OrthogonalPackingInfeasibilityDetector/trivial_conflicts",
53 num_trivial_conflicts_});
54 stats.push_back({"OrthogonalPackingInfeasibilityDetector/conflicts_two_items",
55 num_conflicts_two_items_});
56 stats.push_back({"OrthogonalPackingInfeasibilityDetector/no_energy_conflict",
57 num_scheduling_possible_});
58 stats.push_back({"OrthogonalPackingInfeasibilityDetector/brute_force_calls",
59 num_brute_force_calls_});
60 stats.push_back(
61 {"OrthogonalPackingInfeasibilityDetector/brute_force_conflicts",
62 num_brute_force_conflicts_});
63 stats.push_back(
64 {"OrthogonalPackingInfeasibilityDetector/brute_force_relaxations",
65 num_brute_force_relaxation_});
66
67 shared_stats_->AddStats(stats);
68}
69
70namespace {
71std::optional<std::pair<int, int>> FindPairwiseConflict(
72 absl::Span<const IntegerValue> sizes_x,
73 absl::Span<const IntegerValue> sizes_y,
74 std::pair<IntegerValue, IntegerValue> bounding_box_size,
75 absl::Span<const int> index_by_decreasing_x_size,
76 absl::Span<const int> index_by_decreasing_y_size) {
77 // Look for pairwise incompatible pairs by using the logic such conflict can
78 // only happen between a "tall" item a "wide" item.
79 int x_idx = 0;
80 int y_idx = 0;
81 while (x_idx < index_by_decreasing_x_size.size() &&
82 y_idx < index_by_decreasing_y_size.size()) {
83 if (index_by_decreasing_x_size[x_idx] ==
84 index_by_decreasing_y_size[y_idx]) {
85 if (sizes_x[index_by_decreasing_x_size[x_idx]] >
86 sizes_y[index_by_decreasing_x_size[x_idx]]) {
87 y_idx++;
88 } else {
89 x_idx++;
90 }
91 continue;
92 }
93 const bool overlap_on_x = (sizes_x[index_by_decreasing_x_size[x_idx]] +
94 sizes_x[index_by_decreasing_y_size[y_idx]] >
95 bounding_box_size.first);
96 const bool overlap_on_y = (sizes_y[index_by_decreasing_x_size[x_idx]] +
97 sizes_y[index_by_decreasing_y_size[y_idx]] >
98 bounding_box_size.second);
99 if (overlap_on_x && overlap_on_y) {
100 return std::make_pair(index_by_decreasing_x_size[x_idx],
101 index_by_decreasing_y_size[y_idx]);
102 } else if (overlap_on_x) {
103 x_idx++;
104 } else if (overlap_on_y) {
105 y_idx++;
106 } else {
107 y_idx++;
108 }
109 }
110 return std::nullopt;
111}
112
113IntegerValue RoundingLowestInverse(IntegerValue y, IntegerValue c_k,
114 IntegerValue max_x, IntegerValue k) {
115 DCHECK_GE(y, 0);
116 DCHECK_LE(y, 2 * c_k);
117 IntegerValue ret = std::numeric_limits<IntegerValue>::max();
118
119 // Are we in the case 2 * x == max_x_?
120 if (y <= c_k && (max_x.value() & 1) == 0) {
121 const IntegerValue inverse_mid = max_x / 2;
122 ret = std::min(ret, inverse_mid);
123 if (y == c_k && y.value() & 1) {
124 // This is the only valid case for odd x.
125 return ret;
126 }
127 }
128
129 // The "perfect odd" case is handled above, round up y to an even value.
130 y += y.value() & 1;
131
132 // Check the case 2 * x > max_x_.
133 const IntegerValue inverse_high = max_x - k * (c_k - y / 2);
134 if (2 * inverse_high > max_x) {
135 // We have an inverse in this domain, let's find its minimum value (when
136 // the division rounds down the most) but don't let it go outside the
137 // domain.
138 const IntegerValue lowest_inverse_high =
139 std::max(max_x / 2 + 1, inverse_high - k + 1);
140 ret = std::min(ret, lowest_inverse_high);
141 }
142
143 // Check the case 2 * x < max_x_.
144 const IntegerValue inverse_low = k * y / 2;
145 if (2 * inverse_low < max_x) {
146 ret = std::min(ret, inverse_low);
147 }
148 return ret;
149}
150} // namespace
151
152IntegerValue RoundingDualFeasibleFunction::LowestInverse(IntegerValue y) const {
153 return RoundingLowestInverse(y, c_k_, max_x_, k_);
154}
155
157 IntegerValue y) const {
158 return RoundingLowestInverse(y, c_k_, max_x_, IntegerValue(1) << log2_k_);
159}
160
161// Check for conflict using the `f_0^k` dual feasible function (see
162// documentation for DualFeasibleFunctionF0). This function tries all possible
163// values of the `k` parameter and returns the best conflict found (according to
164// OrthogonalPackingResult::IsBetterThan) if any.
165//
166// The current implementation is a bit more general than a simple check using
167// f_0 described above. This implementation can take a function g(x) that is
168// non-decreasing and satisfy g(0)=0 and it will check for conflict using
169// g(f_0^k(x)) for all values of k, but without recomputing g(x) `k` times. This
170// is handy if g() is a DFF that is slow to compute. g(x) is described by the
171// vector g_x[i] = g(sizes_x[i]) and the variable g_max = g(x_bb_size).
172//
173// The algorithm is the same if we swap the x and y dimension.
174OrthogonalPackingResult OrthogonalPackingInfeasibilityDetector::GetDffConflict(
175 absl::Span<const IntegerValue> sizes_x,
176 absl::Span<const IntegerValue> sizes_y,
177 absl::Span<const int> index_by_decreasing_x_size,
178 absl::Span<const IntegerValue> g_x, IntegerValue g_max,
179 IntegerValue x_bb_size, IntegerValue total_energy, IntegerValue bb_area,
180 IntegerValue* best_k) {
181 // If we found a conflict for a k parameter, which is rare, recompute the
182 // total used energy consumed by the items to find the minimal set of
183 // conflicting items.
184 int num_items = sizes_x.size();
185 auto build_result = [&sizes_x, &sizes_y, num_items, &x_bb_size, &bb_area,
186 &g_max, &g_x](const IntegerValue k) {
187 std::vector<std::pair<int, IntegerValue>> index_to_energy;
188 index_to_energy.reserve(num_items);
189 for (int i = 0; i < num_items; i++) {
190 IntegerValue point_value;
191 if (sizes_x[i] > x_bb_size - k) {
192 point_value = g_max;
193 } else if (sizes_x[i] < k) {
194 continue;
195 } else {
196 point_value = g_x[i];
197 }
198 index_to_energy.push_back({i, point_value * sizes_y[i]});
199 }
200 std::sort(index_to_energy.begin(), index_to_energy.end(),
201 [](const std::pair<int, IntegerValue>& a,
202 const std::pair<int, IntegerValue>& b) {
203 return a.second > b.second;
204 });
205 IntegerValue recomputed_energy = 0;
206 for (int i = 0; i < index_to_energy.size(); i++) {
207 recomputed_energy += index_to_energy[i].second;
208 if (recomputed_energy > bb_area) {
209 OrthogonalPackingResult result(
211 result.conflict_type_ = OrthogonalPackingResult::ConflictType::DFF_F0;
212 result.items_participating_on_conflict_.resize(i + 1);
213 for (int j = 0; j <= i; j++) {
214 const int index = index_to_energy[j].first;
215 result.items_participating_on_conflict_[j] = {
216 .index = index,
217 .size_x = sizes_x[index],
218 .size_y = sizes_y[index]};
219 }
220 result.slack_ = 0;
221 return result;
222 }
223 }
224 LOG(FATAL) << "build_result called with no conflict";
225 };
226
227 // One thing we use in this implementation is that not all values of k are
228 // interesting: what can cause an energy conflict is increasing the size of
229 // the large items, removing the small ones makes it less constrained and we
230 // do it only to preserve correctness. Thus, it is enough to check the values
231 // of k that are just small enough to enlarge a large item. That means that
232 // large items and small ones are not symmetric with respect to what values of
233 // k are important.
234 IntegerValue current_energy = total_energy;
235 OrthogonalPackingResult best_result;
236 if (current_energy > bb_area) {
237 best_result = build_result(0);
238 *best_k = 0;
239 }
240 // We keep an index on the largest item yet-to-be enlarged and the smallest
241 // one yet-to-be removed.
242 int removing_item_index = index_by_decreasing_x_size.size() - 1;
243 int enlarging_item_index = 0;
244 while (enlarging_item_index < index_by_decreasing_x_size.size()) {
245 int index = index_by_decreasing_x_size[enlarging_item_index];
246 IntegerValue size = sizes_x[index];
247 // Note that since `size_x` is decreasing, we test increasingly large
248 // values of k. Also note that a item with size `k` cannot fit alongside
249 // an item with size `size_x`, but smaller ones can.
250 const IntegerValue k = x_bb_size - size + 1;
251 if (2 * k > x_bb_size) {
252 break;
253 }
254 // First, add the area contribution of enlarging the all the items of size
255 // exactly size_x. All larger items were already enlarged in the previous
256 // iterations.
257 do {
258 index = index_by_decreasing_x_size[enlarging_item_index];
259 size = sizes_x[index];
260 current_energy += (g_max - g_x[index]) * sizes_y[index];
261 enlarging_item_index++;
262 } while (enlarging_item_index < index_by_decreasing_x_size.size() &&
263 sizes_x[index_by_decreasing_x_size[enlarging_item_index]] == size);
264
265 // Now remove the area contribution of removing all the items smaller than
266 // k that were not removed before.
267 while (removing_item_index >= 0 &&
268 sizes_x[index_by_decreasing_x_size[removing_item_index]] < k) {
269 const int remove_idx = index_by_decreasing_x_size[removing_item_index];
270 current_energy -= g_x[remove_idx] * sizes_y[remove_idx];
271 removing_item_index--;
272 }
273
274 if (current_energy > bb_area) {
275 OrthogonalPackingResult current_result = build_result(k);
276 if (current_result.IsBetterThan(best_result)) {
277 best_result = current_result;
278 *best_k = k;
279 }
280 }
281 }
282 return best_result;
283}
284
285namespace {
286
287// Tries a simple heuristic to find a solution for the Resource-Constrained
288// Project Scheduling Problem (RCPSP). The RCPSP can be mapped to a
289// 2d bin packing where one dimension (say, x) is chosen to represent the time,
290// and every item is cut into items with size_x = 1 that must remain consecutive
291// in the x-axis but do not need to be aligned on the y axis. This is often
292// called the cumulative relaxation of the 2d bin packing problem.
293//
294// Bin-packing solution RCPSP solution
295// --------------- ---------------
296// | ********** | | ***** |
297// | ********** | | ***** |
298// | ##### | | **#####*** |
299// | ##### | | **#####*** |
300// --------------- ---------------
301//
302// One interesting property is if we find an energy conflict using a
303// superadditive function it means the problem is infeasible both interpreted as
304// a 2d bin packing and as a RCPSP problem. In practice, that means that if we
305// find a RCPSP solution for a 2d bin packing problem, there is no point on
306// using Maximal DFFs to search for energy conflicts.
307//
308// Returns true if it found a feasible solution to the RCPSP problem.
309bool FindHeuristicSchedulingSolution(
310 absl::Span<const IntegerValue> sizes,
311 absl::Span<const IntegerValue> demands,
312 absl::Span<const int> heuristic_order, IntegerValue global_end_max,
313 IntegerValue capacity_max,
314 std::vector<std::pair<IntegerValue, IntegerValue>>& profile,
315 std::vector<std::pair<IntegerValue, IntegerValue>>& new_profile) {
316 // The profile (and new profile) is a set of (time, capa_left) pairs, ordered
317 // by increasing time and capa_left.
318 profile.clear();
319 profile.emplace_back(kMinIntegerValue, capacity_max);
320 profile.emplace_back(kMaxIntegerValue, capacity_max);
321 IntegerValue start_of_previous_task = kMinIntegerValue;
322 for (int i = 0; i < heuristic_order.size(); i++) {
323 const IntegerValue event_size = sizes[heuristic_order[i]];
324 const IntegerValue event_demand = demands[heuristic_order[i]];
325 const IntegerValue event_start_min = 0;
326 const IntegerValue event_start_max = global_end_max - event_size;
327 const IntegerValue start_min =
328 std::max(event_start_min, start_of_previous_task);
329
330 // Iterate on the profile to find the step that contains start_min.
331 // Then push until we find a step with enough capacity.
332 int current = 0;
333 while (profile[current + 1].first <= start_min ||
334 profile[current].second < event_demand) {
335 ++current;
336 }
337
338 const IntegerValue actual_start =
339 std::max(start_min, profile[current].first);
340 start_of_previous_task = actual_start;
341
342 // Compatible with the event.start_max ?
343 if (actual_start > event_start_max) return false;
344
345 const IntegerValue actual_end = actual_start + event_size;
346
347 // No need to update the profile on the last loop.
348 if (i == heuristic_order.size() - 1) break;
349
350 // Update the profile.
351 new_profile.clear();
352 new_profile.push_back(
353 {actual_start, profile[current].second - event_demand});
354 ++current;
355
356 while (profile[current].first < actual_end) {
357 new_profile.push_back(
358 {profile[current].first, profile[current].second - event_demand});
359 ++current;
360 }
361
362 if (profile[current].first > actual_end) {
363 new_profile.push_back(
364 {actual_end, new_profile.back().second + event_demand});
365 }
366 while (current < profile.size()) {
367 new_profile.push_back(profile[current]);
368 ++current;
369 }
370 profile.swap(new_profile);
371 }
372 return true;
373}
374
375} // namespace
376
377// We want to find the minimum set of values of `k` that would always find a
378// conflict if there is a `k` for which it exists. In the literature it is
379// often implied (but not stated) that it is sufficient to test the values of
380// `k` that correspond to the size of an item. This is not true. To find the
381// minimum set of values of `k` we look for all values of `k` that are
382// "extreme": ie., the rounding on the division truncates the most (or the
383// least) amount, depending on the sign it appears in the formula.
384//
385// To find these extreme values, we look for all local minima of the energy
386// slack after applying the DFF (we multiply by `k` for convenience):
387// k * f_k(H) * W - sum_i k * f_k(h_i) * w_i
388// If this value ever becomes negative for a value of `k`, it must happen in a
389// local minimum. Then we use the fact that
390// k * floor(x / k) = x - x % k
391// and that x%k has a local minimum when k=x/i and a local maximum when k=1+x/i
392// for every integer i. The final finer point in the calculation is
393// realizing that if
394// sum_{i, h_i > H/2} w_i > W
395// then you have more "large" objects than it fits in the box, and you will
396// have a conflict using the DFF f_0 for l=H/2. So we can safely ignore this
397// case for the more expensive DFF f_2 calculation.
398void OrthogonalPackingInfeasibilityDetector::GetAllCandidatesForKForDff2(
399 absl::Span<const IntegerValue> sizes, IntegerValue bb_size,
400 IntegerValue sqrt_bb_size, Bitset64<IntegerValue>& candidates) {
401 // x_bb_size is less than 65536, so this fits in only 4kib.
402 candidates.ClearAndResize(bb_size / 2 + 2);
403
404 // `sqrt_bb_size` is lower than 256.
405 for (IntegerValue i = 2; i <= sqrt_bb_size; i++) {
406 candidates.Set(i);
407 }
408 for (int i = 1; i <= sqrt_bb_size; i++) {
409 const ::util::math::ConstantDivisor<uint16_t> div(i);
410 if (i > 1) {
411 candidates.Set(bb_size.value() / div);
412 }
413 for (int k = 0; k < sizes.size(); k++) {
414 IntegerValue size = sizes[k];
415 if (2 * size > bb_size && size < bb_size) {
416 candidates.Set((bb_size.value() - size.value() + 1) / div);
417 } else if (2 * size < bb_size) {
418 candidates.Set(size.value() / div);
419 }
420 }
421 }
422
423 // Remove some bogus candidates added by the logic above.
424 candidates.Clear(0);
425 candidates.Clear(1);
426
427 // Apply the nice result described on [1]: if we are testing the DFF
428 // f_2^k(f_0^l(x)) for all values of `l`, the only values of `k` greater than
429 // C/4 we need to test are {C/4+1, C/3+1}.
430 //
431 // In the same reference there is a proof that this way of composing f_0 and
432 // f_2 cover all possible ways of composing the two functions, including
433 // composing several times each.
434 //
435 // [1] F. Clautiaux, PhD thesis, hal/tel-00749411.
436 candidates.Resize(bb_size / 4 + 1); // Erase all >= C/4
437 candidates.Resize(bb_size / 3 + 2); // Make room for the two special values
438 candidates.Set(bb_size / 4 + 1);
439 if (bb_size > 3) {
440 candidates.Set(bb_size / 3 + 1);
441 }
442}
443
444// Check for conflict all combinations of the two Dual Feasible Functions f_0
445// (see documentation for GetDffConflict()) and f_2 (see documentation for
446// RoundingDualFeasibleFunction). More precisely, check whether there exist `l`
447// and `k` so that
448//
449// sum_i f_2^k(f_0^l(sizes_x[i])) * sizes_y[i] > f_2^k(f_0^l(x_bb_size)) *
450// y_bb_size
451//
452// The function returns the smallest subset of items enough to make the
453// inequality above true or an empty vector if impossible.
455OrthogonalPackingInfeasibilityDetector::CheckFeasibilityWithDualFunction2(
456 absl::Span<const IntegerValue> sizes_x,
457 absl::Span<const IntegerValue> sizes_y,
458 absl::Span<const int> index_by_decreasing_x_size, IntegerValue x_bb_size,
459 IntegerValue y_bb_size, int max_number_of_parameters_to_check) {
460 if (x_bb_size == 1) {
461 return OrthogonalPackingResult();
462 }
463 std::vector<IntegerValue> sizes_x_rescaled;
464 if (x_bb_size >= std::numeric_limits<uint16_t>::max()) {
465 // To do fast division we want our sizes to fit in a uint16_t. The simplest
466 // way of doing that is to just first apply this DFF with the right
467 // power-of-two value of the parameter.
468 const int log2_k =
469 absl::bit_width(static_cast<uint64_t>(x_bb_size.value() + 1)) - 16 + 1;
470 const RoundingDualFeasibleFunctionPowerOfTwo dff(x_bb_size, log2_k);
471 sizes_x_rescaled.resize(sizes_x.size());
472 for (int i = 0; i < sizes_x.size(); i++) {
473 sizes_x_rescaled[i] = dff(sizes_x[i]);
474 }
475 x_bb_size = dff(x_bb_size);
476 CHECK_LT(x_bb_size, std::numeric_limits<uint16_t>::max());
477 sizes_x = sizes_x_rescaled;
478 }
479
480 Bitset64<IntegerValue> candidates;
481 const IntegerValue sqrt_bb_size = FloorSquareRoot(x_bb_size.value());
482 int num_items = sizes_x.size();
483 const IntegerValue max_possible_number_of_parameters =
484 std::min(x_bb_size / 4 + 1, sqrt_bb_size * num_items);
485 if (5ull * max_number_of_parameters_to_check <
486 max_possible_number_of_parameters) {
487 // There are many more possible values than what we want to sample. It is
488 // not worth to pay the price of computing all optimal values to drop most
489 // of them, so let's just pick it randomly.
490 candidates.Resize(x_bb_size / 4 + 1);
491 int num_candidates = 0;
492 while (num_candidates < max_number_of_parameters_to_check) {
493 const IntegerValue pick =
494 absl::Uniform(random_, 1, x_bb_size.value() / 4);
495 if (!candidates.IsSet(pick)) {
496 candidates.Set(pick);
497 num_candidates++;
498 }
499 }
500 } else {
501 GetAllCandidatesForKForDff2(sizes_x, x_bb_size, sqrt_bb_size, candidates);
502
503 if (max_number_of_parameters_to_check < max_possible_number_of_parameters) {
504 // We might have produced too many candidates. Let's count them and if it
505 // is the case, sample them.
506 int count = 0;
507 for (auto it = candidates.begin(); it != candidates.end(); ++it) {
508 count++;
509 }
510 if (count > max_number_of_parameters_to_check) {
511 std::vector<IntegerValue> sampled_candidates(
512 max_number_of_parameters_to_check);
513 std::sample(candidates.begin(), candidates.end(),
514 sampled_candidates.begin(),
515 max_number_of_parameters_to_check, random_);
516 candidates.ClearAll();
517 for (const IntegerValue k : sampled_candidates) {
518 candidates.Set(k);
519 }
520 }
521 }
522 }
523 OrthogonalPackingResult best_result;
524
525 // Finally run our small loop to look for the conflict!
526 std::vector<IntegerValue> modified_sizes(num_items);
527 for (const IntegerValue k : candidates) {
528 const RoundingDualFeasibleFunction dff(x_bb_size, k);
529 IntegerValue energy = 0;
530 for (int i = 0; i < num_items; i++) {
531 modified_sizes[i] = dff(sizes_x[i]);
532 energy += modified_sizes[i] * sizes_y[i];
533 }
534 const IntegerValue modified_x_bb_size = dff(x_bb_size);
535 IntegerValue dff0_k;
536 auto dff0_res =
537 GetDffConflict(sizes_x, sizes_y, index_by_decreasing_x_size,
538 modified_sizes, modified_x_bb_size, x_bb_size, energy,
539 modified_x_bb_size * y_bb_size, &dff0_k);
540 if (dff0_res.result_ != OrthogonalPackingResult::Status::INFEASIBLE) {
541 continue;
542 }
543 DFFComposedF2F0 composed_dff(x_bb_size, dff0_k, k);
544 dff0_res.conflict_type_ = OrthogonalPackingResult::ConflictType::DFF_F2;
545 for (auto& item : dff0_res.items_participating_on_conflict_) {
546 item.size_x =
547 composed_dff.LowestInverse(composed_dff(sizes_x[item.index]));
548
549 // The new size should contribute by the same amount to the energy and
550 // correspond to smaller items.
551 DCHECK_EQ(composed_dff(item.size_x), composed_dff(sizes_x[item.index]));
552 DCHECK_LE(item.size_x, sizes_x[item.index]);
553
554 item.size_y = sizes_y[item.index];
555 }
556 if (dff0_res.IsBetterThan(best_result)) {
557 best_result = dff0_res;
558 }
559 }
560
561 return best_result;
562}
563
564bool OrthogonalPackingInfeasibilityDetector::RelaxConflictWithBruteForce(
566 std::pair<IntegerValue, IntegerValue> bounding_box_size,
567 int brute_force_threshold) {
568 const int num_items_originally =
569 result.items_participating_on_conflict_.size();
570 if (num_items_originally > 2 * brute_force_threshold) {
571 // Don't even try on problems too big.
572 return false;
573 }
574 std::vector<IntegerValue> sizes_x;
575 std::vector<IntegerValue> sizes_y;
576 std::vector<int> indexes;
577 std::vector<bool> to_be_removed(num_items_originally, false);
578
579 sizes_x.reserve(num_items_originally - 1);
580 sizes_y.reserve(num_items_originally - 1);
581 for (int i = 0; i < num_items_originally; i++) {
582 sizes_x.clear();
583 sizes_y.clear();
584 // Look for a conflict using all non-removed items but the i-th one.
585 for (int j = 0; j < num_items_originally; j++) {
586 if (i == j || to_be_removed[j]) {
587 continue;
588 }
589 sizes_x.push_back(result.items_participating_on_conflict_[j].size_x);
590 sizes_y.push_back(result.items_participating_on_conflict_[j].size_y);
591 }
593 sizes_x, sizes_y, bounding_box_size, brute_force_threshold);
595 // We still have a conflict if we remove the i-th item!
596 to_be_removed[i] = true;
597 }
598 }
599 if (!std::any_of(to_be_removed.begin(), to_be_removed.end(),
600 [](bool b) { return b; })) {
601 return false;
602 }
603 OrthogonalPackingResult original = result;
604 result.slack_ = 0;
605 result.conflict_type_ = OrthogonalPackingResult::ConflictType::BRUTE_FORCE;
607 result.items_participating_on_conflict_.clear();
608 for (int i = 0; i < num_items_originally; i++) {
609 if (to_be_removed[i]) {
610 continue;
611 }
612 result.items_participating_on_conflict_.push_back(
613 original.items_participating_on_conflict_[i]);
614 }
615 return true;
616}
617
619OrthogonalPackingInfeasibilityDetector::TestFeasibilityImpl(
620 absl::Span<const IntegerValue> sizes_x,
621 absl::Span<const IntegerValue> sizes_y,
622 std::pair<IntegerValue, IntegerValue> bounding_box_size,
623 const OrthogonalPackingOptions& options) {
624 using ConflictType = OrthogonalPackingResult::ConflictType;
625
626 const int num_items = sizes_x.size();
627 DCHECK_EQ(num_items, sizes_y.size());
628 const IntegerValue bb_area =
629 bounding_box_size.first * bounding_box_size.second;
630 IntegerValue total_energy = 0;
631
632 auto make_item = [sizes_x, sizes_y](int i) {
633 return OrthogonalPackingResult::Item{
634 .index = i, .size_x = sizes_x[i], .size_y = sizes_y[i]};
635 };
636
637 index_by_decreasing_x_size_.resize(num_items);
638 index_by_decreasing_y_size_.resize(num_items);
639 for (int i = 0; i < num_items; i++) {
640 total_energy += sizes_x[i] * sizes_y[i];
641 index_by_decreasing_x_size_[i] = i;
642 index_by_decreasing_y_size_[i] = i;
643 if (sizes_x[i] > bounding_box_size.first ||
644 sizes_y[i] > bounding_box_size.second) {
645 OrthogonalPackingResult result(
647 result.conflict_type_ = ConflictType::TRIVIAL;
648 result.items_participating_on_conflict_ = {make_item(i)};
649 return result;
650 }
651 }
652
653 if (num_items <= 1) {
654 return OrthogonalPackingResult(OrthogonalPackingResult::Status::FEASIBLE);
655 }
656
657 std::sort(index_by_decreasing_x_size_.begin(),
658 index_by_decreasing_x_size_.end(),
659 [&sizes_x, &sizes_y](int a, int b) {
660 // Break ties with y-size
661 return std::tie(sizes_x[a], sizes_y[a]) >
662 std::tie(sizes_x[b], sizes_y[b]);
663 });
664 std::sort(index_by_decreasing_y_size_.begin(),
665 index_by_decreasing_y_size_.end(),
666 [&sizes_y, &sizes_x](int a, int b) {
667 return std::tie(sizes_y[a], sizes_x[a]) >
668 std::tie(sizes_y[b], sizes_x[b]);
669 });
670
671 // First look for pairwise incompatible pairs.
672 if (options.use_pairwise) {
673 if (auto pair = FindPairwiseConflict(sizes_x, sizes_y, bounding_box_size,
674 index_by_decreasing_x_size_,
675 index_by_decreasing_y_size_);
676 pair.has_value()) {
677 OrthogonalPackingResult result(
679 result.conflict_type_ = ConflictType::PAIRWISE;
680 result.items_participating_on_conflict_ = {
681 make_item(pair.value().first), make_item(pair.value().second)};
682 return result;
683 }
684 if (num_items == 2) {
685 return OrthogonalPackingResult(OrthogonalPackingResult::Status::FEASIBLE);
686 }
687 }
688
689 OrthogonalPackingResult result(OrthogonalPackingResult::Status::UNKNOWN);
690 if (total_energy > bb_area) {
691 result.conflict_type_ = ConflictType::TRIVIAL;
693 std::vector<std::pair<int, IntegerValue>> index_to_energy;
694 index_to_energy.reserve(num_items);
695 for (int i = 0; i < num_items; i++) {
696 index_to_energy.push_back({i, sizes_x[i] * sizes_y[i]});
697 }
698 std::sort(index_to_energy.begin(), index_to_energy.end(),
699 [](const std::pair<int, IntegerValue>& a,
700 const std::pair<int, IntegerValue>& b) {
701 return a.second > b.second;
702 });
703 IntegerValue recomputed_energy = 0;
704 for (int i = 0; i < index_to_energy.size(); i++) {
705 recomputed_energy += index_to_energy[i].second;
706 if (recomputed_energy > bb_area) {
707 result.items_participating_on_conflict_.resize(i + 1);
708 for (int j = 0; j <= i; j++) {
709 result.items_participating_on_conflict_[j] =
710 make_item(index_to_energy[j].first);
711 }
712 result.slack_ = recomputed_energy - bb_area - 1;
713 break;
714 }
715 }
716 }
717
718 const int minimum_conflict_size = options.use_pairwise ? 3 : 2;
719 if (result.items_participating_on_conflict_.size() == minimum_conflict_size) {
720 return result;
721 }
722
723 if (options.use_dff_f0) {
724 // If there is no pairwise incompatible pairs, this DFF cannot find a
725 // conflict by enlarging a item on both x and y directions: this would
726 // create an item as long as the whole box and another item as high as the
727 // whole box, which is obviously incompatible, and this incompatibility
728 // would be present already before enlarging the items since it is a DFF. So
729 // it is enough to test making items wide or high, but no need to try both.
730 IntegerValue best_k;
731 auto conflict =
732 GetDffConflict(sizes_x, sizes_y, index_by_decreasing_x_size_, sizes_x,
733 bounding_box_size.first, bounding_box_size.first,
734 total_energy, bb_area, &best_k);
735 if (conflict.IsBetterThan(result)) {
736 result = conflict;
737 }
738
739 conflict =
740 GetDffConflict(sizes_y, sizes_x, index_by_decreasing_y_size_, sizes_y,
741 bounding_box_size.second, bounding_box_size.second,
742 total_energy, bb_area, &best_k);
743 for (auto& item : conflict.items_participating_on_conflict_) {
744 std::swap(item.size_x, item.size_y);
745 }
746 if (conflict.IsBetterThan(result)) {
747 result = conflict;
748 }
749 }
750
751 if (result.items_participating_on_conflict_.size() == minimum_conflict_size) {
752 return result;
753 }
754
755 bool found_scheduling_solution = false;
756 if (options.use_dff_f2) {
757 // Checking for conflicts using f_2 is expensive, so first try a quick
758 // algorithm to check if there is no conflict to be found. See the comments
759 // on top of FindHeuristicSchedulingSolution().
760 if (FindHeuristicSchedulingSolution(
761 sizes_x, sizes_y, index_by_decreasing_x_size_,
762 bounding_box_size.first, bounding_box_size.second,
763 scheduling_profile_, new_scheduling_profile_) ||
764 FindHeuristicSchedulingSolution(
765 sizes_y, sizes_x, index_by_decreasing_y_size_,
766 bounding_box_size.second, bounding_box_size.first,
767 scheduling_profile_, new_scheduling_profile_)) {
768 num_scheduling_possible_++;
769 CHECK(result.result_ != OrthogonalPackingResult::Status::INFEASIBLE);
770 found_scheduling_solution = true;
771 }
772 }
773
774 if (!found_scheduling_solution && options.use_dff_f2) {
775 // We only check for conflicts applying this DFF on heights and widths, but
776 // not on both, which would be too expensive if done naively.
777 auto conflict = CheckFeasibilityWithDualFunction2(
778 sizes_x, sizes_y, index_by_decreasing_x_size_, bounding_box_size.first,
779 bounding_box_size.second,
780 options.dff2_max_number_of_parameters_to_check);
781 if (conflict.IsBetterThan(result)) {
782 result = conflict;
783 }
784
785 if (result.items_participating_on_conflict_.size() ==
786 minimum_conflict_size) {
787 return result;
788 }
789 conflict = CheckFeasibilityWithDualFunction2(
790 sizes_y, sizes_x, index_by_decreasing_y_size_, bounding_box_size.second,
791 bounding_box_size.first,
792 options.dff2_max_number_of_parameters_to_check);
793 for (auto& item : conflict.items_participating_on_conflict_) {
794 std::swap(item.size_x, item.size_y);
795 }
796 if (conflict.IsBetterThan(result)) {
797 result = conflict;
798 }
799 }
800
801 if (result.result_ == OrthogonalPackingResult::Status::UNKNOWN) {
803 sizes_x, sizes_y, bounding_box_size, options.brute_force_threshold);
804 num_brute_force_calls_ +=
807 result.conflict_type_ = ConflictType::BRUTE_FORCE;
809 result.items_participating_on_conflict_.resize(num_items);
810 for (int i = 0; i < num_items; i++) {
811 result.items_participating_on_conflict_[i] = make_item(i);
812 }
815 }
816 }
817
818 if (result.result_ == OrthogonalPackingResult::Status::INFEASIBLE) {
819 num_brute_force_relaxation_ += RelaxConflictWithBruteForce(
820 result, bounding_box_size, options.brute_force_threshold);
821 }
822
823 return result;
824}
825
827 absl::Span<const IntegerValue> sizes_x,
828 absl::Span<const IntegerValue> sizes_y,
829 std::pair<IntegerValue, IntegerValue> bounding_box_size,
830 const OrthogonalPackingOptions& options) {
831 using ConflictType = OrthogonalPackingResult::ConflictType;
832
833 num_calls_++;
835 TestFeasibilityImpl(sizes_x, sizes_y, bounding_box_size, options);
836
837 if (result.result_ == OrthogonalPackingResult::Status::INFEASIBLE) {
838 num_conflicts_++;
839 switch (result.conflict_type_) {
840 case ConflictType::DFF_F0:
841 num_conflicts_dff0_++;
842 break;
843 case ConflictType::DFF_F2:
844 num_conflicts_dff2_++;
845 break;
846 case ConflictType::PAIRWISE:
847 num_conflicts_two_items_++;
848 break;
849 case ConflictType::TRIVIAL:
850 // The total area of the items was larger than the area of the box.
851 num_trivial_conflicts_++;
852 break;
853 case ConflictType::BRUTE_FORCE:
854 num_brute_force_conflicts_++;
855 break;
856 case ConflictType::NO_CONFLICT:
857 LOG(FATAL) << "Should never happen";
858 break;
859 }
860 }
861 return result;
862}
863
865 int i, Coord coord, IntegerValue lower_bound) {
866 Item& item = items_participating_on_conflict_[i];
867 IntegerValue& size = coord == Coord::kCoordX ? item.size_x : item.size_y;
868 const IntegerValue orthogonal_size =
869 coord == Coord::kCoordX ? item.size_y : item.size_x;
870
871 if (size <= lower_bound || orthogonal_size > slack_) {
872 return false;
873 }
874 const IntegerValue new_size =
875 std::max(lower_bound, size - slack_ / orthogonal_size);
876 slack_ -= (size - new_size) * orthogonal_size;
877 DCHECK_NE(size, new_size);
878 DCHECK_GE(slack_, 0);
879 size = new_size;
880 return true;
881}
882
883} // namespace sat
884} // namespace operations_research
OrthogonalPackingResult TestFeasibility(absl::Span< const IntegerValue > sizes_x, absl::Span< const IntegerValue > sizes_y, std::pair< IntegerValue, IntegerValue > bounding_box_size, const OrthogonalPackingOptions &options=OrthogonalPackingOptions())
bool TryUseSlackToReduceItemSize(int i, Coord coord, IntegerValue lower_bound=0)
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
int64_t FloorSquareRoot(int64_t a)
The argument must be non-negative.
Definition util.cc:300
BruteForceResult BruteForceOrthogonalPacking(absl::Span< const IntegerValue > sizes_x, absl::Span< const IntegerValue > sizes_y, std::pair< IntegerValue, IntegerValue > bounding_box_size, int max_complexity)
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue.value())
In SWIG mode, we don't want anything besides these top-level includes.
Select next search node to expand Select next item_i to add this new search node to the search Generate a new search node where item_i is not in the knapsack Check validity of this new partial solution(using propagators) - If valid