Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
diffn_util.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 <stddef.h>
17
18#include <algorithm>
19#include <array>
20#include <cmath>
21#include <cstddef>
22#include <cstdint>
23#include <limits>
24#include <optional>
25#include <ostream>
26#include <sstream>
27#include <string>
28#include <string_view>
29#include <tuple>
30#include <type_traits>
31#include <utility>
32#include <vector>
33
34#include "absl/algorithm/container.h"
35#include "absl/container/btree_set.h"
36#include "absl/container/flat_hash_map.h"
37#include "absl/container/flat_hash_set.h"
38#include "absl/container/inlined_vector.h"
39#include "absl/log/check.h"
40#include "absl/log/log.h"
41#include "absl/log/vlog_is_on.h"
42#include "absl/random/bit_gen_ref.h"
43#include "absl/types/optional.h"
44#include "absl/types/span.h"
51#include "ortools/sat/util.h"
55
56namespace operations_research {
57namespace sat {
58
59bool Rectangle::IsDisjoint(const Rectangle& other) const {
60 return x_min >= other.x_max || other.x_min >= x_max || y_min >= other.y_max ||
61 other.y_min >= y_max;
62}
63
64absl::InlinedVector<Rectangle, 4> Rectangle::RegionDifference(
65 const Rectangle& other) const {
66 const Rectangle intersect = Intersect(other);
67 if (intersect.SizeX() == 0) {
68 return {*this};
69 }
70
71 //-------------------
72 //| | 4 | |
73 //| |---------| |
74 //| 1 | other | 2 |
75 //| |---------| |
76 //| | 3 | |
77 //-------------------
78 absl::InlinedVector<Rectangle, 4> result;
79 if (x_min < intersect.x_min) {
80 // Piece 1
81 result.push_back({.x_min = x_min,
82 .x_max = intersect.x_min,
83 .y_min = y_min,
84 .y_max = y_max});
85 }
86 if (x_max > intersect.x_max) {
87 // Piece 2
88 result.push_back({.x_min = intersect.x_max,
89 .x_max = x_max,
90 .y_min = y_min,
91 .y_max = y_max});
92 }
93 if (y_min < intersect.y_min) {
94 // Piece 3
95 result.push_back({.x_min = intersect.x_min,
96 .x_max = intersect.x_max,
97 .y_min = y_min,
98 .y_max = intersect.y_min});
99 }
100 if (y_max > intersect.y_max) {
101 // Piece 4
102 result.push_back({.x_min = intersect.x_min,
103 .x_max = intersect.x_max,
104 .y_min = intersect.y_max,
105 .y_max = y_max});
106 }
107
108 return result;
109}
110
112 absl::Span<const Rectangle> rectangles) {
113 if (rectangles.empty()) return {};
114
115 std::vector<std::pair<int, int>> intersections =
117 const int num_intersections = intersections.size();
118 intersections.reserve(num_intersections * 2 + 1);
119 for (int i = 0; i < num_intersections; ++i) {
120 intersections.push_back({intersections[i].second, intersections[i].first});
121 }
122
124 view.ResetFromPairs(intersections, /*minimum_num_nodes=*/rectangles.size());
125 CompactVectorVector<int> components;
126 FindStronglyConnectedComponents(static_cast<int>(rectangles.size()), view,
127 &components);
129 for (int i = 0; i < components.size(); ++i) {
130 absl::Span<const int> component = components[i];
131 result.Add({});
132 for (const int r : component) {
133 result.AppendToLastVector(r);
134 }
135 }
136 return result;
137}
138
139bool ReportEnergyConflict(Rectangle bounding_box, absl::Span<const int> boxes,
142 x->ClearReason();
143 y->ClearReason();
144 IntegerValue total_energy(0);
145 for (const int b : boxes) {
146 const IntegerValue x_min = x->ShiftedStartMin(b);
147 const IntegerValue x_max = x->ShiftedEndMax(b);
148 if (x_min < bounding_box.x_min || x_max > bounding_box.x_max) continue;
149 const IntegerValue y_min = y->ShiftedStartMin(b);
150 const IntegerValue y_max = y->ShiftedEndMax(b);
151 if (y_min < bounding_box.y_min || y_max > bounding_box.y_max) continue;
152
153 x->AddEnergyMinInIntervalReason(b, bounding_box.x_min, bounding_box.x_max);
154 y->AddEnergyMinInIntervalReason(b, bounding_box.y_min, bounding_box.y_max);
155
156 x->AddPresenceReason(b);
158
159 total_energy += x->SizeMin(b) * y->SizeMin(b);
160
161 // We abort early if a subset of boxes is enough.
162 // TODO(user): Also relax the box if possible.
163 if (total_energy > bounding_box.Area()) break;
164 }
165
166 CHECK_GT(total_energy, bounding_box.Area());
167 x->ImportOtherReasons(*y);
168 return x->ReportConflict();
169}
170
171bool BoxesAreInEnergyConflict(absl::Span<const Rectangle> rectangles,
172 absl::Span<const IntegerValue> energies,
173 absl::Span<const int> boxes,
174 Rectangle* conflict) {
175 // First consider all relevant intervals along the x axis.
176 std::vector<IntegerValue> x_starts;
177 std::vector<TaskTime> boxes_by_increasing_x_max;
178 for (const int b : boxes) {
179 x_starts.push_back(rectangles[b].x_min);
180 boxes_by_increasing_x_max.push_back({b, rectangles[b].x_max});
181 }
183 std::sort(boxes_by_increasing_x_max.begin(), boxes_by_increasing_x_max.end());
184
185 std::vector<IntegerValue> y_starts;
186 std::vector<IntegerValue> energy_sum;
187 std::vector<TaskTime> boxes_by_increasing_y_max;
188
189 std::vector<std::vector<int>> stripes(x_starts.size());
190 for (int i = 0; i < boxes_by_increasing_x_max.size(); ++i) {
191 const int b = boxes_by_increasing_x_max[i].task_index;
192 const IntegerValue x_min = rectangles[b].x_min;
193 const IntegerValue x_max = rectangles[b].x_max;
194 for (int j = 0; j < x_starts.size(); ++j) {
195 if (x_starts[j] > x_min) break;
196 stripes[j].push_back(b);
197
198 // Redo the same on the y coordinate for the current x interval
199 // which is [starts[j], x_max].
200 y_starts.clear();
201 boxes_by_increasing_y_max.clear();
202 for (const int b : stripes[j]) {
203 y_starts.push_back(rectangles[b].y_min);
204 boxes_by_increasing_y_max.push_back({b, rectangles[b].y_max});
205 }
207 std::sort(boxes_by_increasing_y_max.begin(),
208 boxes_by_increasing_y_max.end());
209
210 const IntegerValue x_size = x_max - x_starts[j];
211 energy_sum.assign(y_starts.size(), IntegerValue(0));
212 for (int i = 0; i < boxes_by_increasing_y_max.size(); ++i) {
213 const int b = boxes_by_increasing_y_max[i].task_index;
214 const IntegerValue y_min = rectangles[b].y_min;
215 const IntegerValue y_max = rectangles[b].y_max;
216 for (int j = 0; j < y_starts.size(); ++j) {
217 if (y_starts[j] > y_min) break;
218 energy_sum[j] += energies[b];
219 if (energy_sum[j] > x_size * (y_max - y_starts[j])) {
220 if (conflict != nullptr) {
221 *conflict = rectangles[b];
222 for (int k = 0; k < i; ++k) {
223 const int task_index = boxes_by_increasing_y_max[k].task_index;
224 if (rectangles[task_index].y_min >= y_starts[j]) {
225 conflict->GrowToInclude(rectangles[task_index]);
226 }
227 }
228 }
229 return true;
230 }
231 }
232 }
233 }
234 }
235 return false;
236}
237
238bool AnalyzeIntervals(bool transpose, absl::Span<const int> local_boxes,
239 absl::Span<const Rectangle> rectangles,
240 absl::Span<const IntegerValue> rectangle_energies,
241 IntegerValue* x_threshold, IntegerValue* y_threshold,
242 Rectangle* conflict) {
243 // First, we compute the possible x_min values (removing duplicates).
244 // We also sort the relevant tasks by their x_max.
245 //
246 // TODO(user): If the number of unique x_max is smaller than the number of
247 // unique x_min, it is better to do it the other way around.
248 std::vector<IntegerValue> starts;
249 std::vector<TaskTime> task_by_increasing_x_max;
250 for (const int t : local_boxes) {
251 const IntegerValue x_min =
252 transpose ? rectangles[t].y_min : rectangles[t].x_min;
253 const IntegerValue x_max =
254 transpose ? rectangles[t].y_max : rectangles[t].x_max;
255 starts.push_back(x_min);
256 task_by_increasing_x_max.push_back({t, x_max});
257 }
259
260 // Note that for the same end_max, the order change our heuristic to
261 // evaluate the max_conflict_height.
262 std::sort(task_by_increasing_x_max.begin(), task_by_increasing_x_max.end());
263
264 // The maximum y dimension of a bounding area for which there is a potential
265 // conflict.
266 IntegerValue max_conflict_height(0);
267
268 // This is currently only used for logging.
269 absl::flat_hash_set<std::pair<IntegerValue, IntegerValue>> stripes;
270
271 // All quantities at index j correspond to the interval [starts[j], x_max].
272 std::vector<IntegerValue> energies(starts.size(), IntegerValue(0));
273 std::vector<IntegerValue> y_mins(starts.size(), kMaxIntegerValue);
274 std::vector<IntegerValue> y_maxs(starts.size(), -kMaxIntegerValue);
275 std::vector<IntegerValue> energy_at_max_y(starts.size(), IntegerValue(0));
276 std::vector<IntegerValue> energy_at_min_y(starts.size(), IntegerValue(0));
277
278 // Sentinel.
279 starts.push_back(kMaxIntegerValue);
280
281 // Iterate over all boxes by increasing x_max values.
282 int first_j = 0;
283 const IntegerValue threshold = transpose ? *y_threshold : *x_threshold;
284 for (int i = 0; i < task_by_increasing_x_max.size(); ++i) {
285 const int t = task_by_increasing_x_max[i].task_index;
286
287 const IntegerValue energy = rectangle_energies[t];
288 IntegerValue x_min = rectangles[t].x_min;
289 IntegerValue x_max = rectangles[t].x_max;
290 IntegerValue y_min = rectangles[t].y_min;
291 IntegerValue y_max = rectangles[t].y_max;
292 if (transpose) {
293 std::swap(x_min, y_min);
294 std::swap(x_max, y_max);
295 }
296
297 // Add this box contribution to all the [starts[j], x_max] intervals.
298 while (first_j + 1 < starts.size() && x_max - starts[first_j] > threshold) {
299 ++first_j;
300 }
301 for (int j = first_j; starts[j] <= x_min; ++j) {
302 const IntegerValue old_energy_at_max = energy_at_max_y[j];
303 const IntegerValue old_energy_at_min = energy_at_min_y[j];
304
305 energies[j] += energy;
306
307 const bool is_disjoint = y_min >= y_maxs[j] || y_max <= y_mins[j];
308
309 if (y_min <= y_mins[j]) {
310 if (y_min < y_mins[j]) {
311 y_mins[j] = y_min;
312 energy_at_min_y[j] = energy;
313 } else {
314 energy_at_min_y[j] += energy;
315 }
316 }
317
318 if (y_max >= y_maxs[j]) {
319 if (y_max > y_maxs[j]) {
320 y_maxs[j] = y_max;
321 energy_at_max_y[j] = energy;
322 } else {
323 energy_at_max_y[j] += energy;
324 }
325 }
326
327 // If the new box is disjoint in y from the ones added so far, there
328 // cannot be a new conflict involving this box, so we skip until we add
329 // new boxes.
330 if (is_disjoint) continue;
331
332 const IntegerValue width = x_max - starts[j];
333 IntegerValue conflict_height = CeilRatio(energies[j], width) - 1;
334 if (y_max - y_min > conflict_height) continue;
335 if (conflict_height >= y_maxs[j] - y_mins[j]) {
336 // We have a conflict.
337 if (conflict != nullptr) {
338 *conflict = rectangles[t];
339 for (int k = 0; k < i; ++k) {
340 const int task_index = task_by_increasing_x_max[k].task_index;
341 const IntegerValue task_x_min = transpose
342 ? rectangles[task_index].y_min
343 : rectangles[task_index].x_min;
344 if (task_x_min < starts[j]) continue;
345 conflict->GrowToInclude(rectangles[task_index]);
346 }
347 }
348 return false;
349 }
350
351 // Because we currently do not have a conflict involving the new box, the
352 // only way to have one is to remove enough energy to reduce the y domain.
353 IntegerValue can_remove = std::min(old_energy_at_min, old_energy_at_max);
354 if (old_energy_at_min < old_energy_at_max) {
355 if (y_maxs[j] - y_min >=
356 CeilRatio(energies[j] - old_energy_at_min, width)) {
357 // In this case, we need to remove at least old_energy_at_max to have
358 // a conflict.
359 can_remove = old_energy_at_max;
360 }
361 } else if (old_energy_at_max < old_energy_at_min) {
362 if (y_max - y_mins[j] >=
363 CeilRatio(energies[j] - old_energy_at_max, width)) {
364 can_remove = old_energy_at_min;
365 }
366 }
367 conflict_height = CeilRatio(energies[j] - can_remove, width) - 1;
368
369 // If the new box height is above the conflict_height, do not count
370 // it now. We only need to consider conflict involving the new box.
371 if (y_max - y_min > conflict_height) continue;
372
373 if (VLOG_IS_ON(2)) stripes.insert({starts[j], x_max});
374 max_conflict_height = std::max(max_conflict_height, conflict_height);
375 }
376 }
377
378 VLOG(2) << " num_starts: " << starts.size() - 1 << "/" << local_boxes.size()
379 << " conflict_height: " << max_conflict_height
380 << " num_stripes:" << stripes.size() << " (<= " << threshold << ")";
381
382 if (transpose) {
383 *x_threshold = std::min(*x_threshold, max_conflict_height);
384 } else {
385 *y_threshold = std::min(*y_threshold, max_conflict_height);
386 }
387 return true;
388}
389
390absl::Span<int> FilterBoxesAndRandomize(
391 absl::Span<const Rectangle> cached_rectangles, absl::Span<int> boxes,
392 IntegerValue threshold_x, IntegerValue threshold_y,
393 absl::BitGenRef random) {
394 size_t new_size = 0;
395 for (const int b : boxes) {
396 const Rectangle& dim = cached_rectangles[b];
397 if (dim.x_max - dim.x_min > threshold_x) continue;
398 if (dim.y_max - dim.y_min > threshold_y) continue;
399 boxes[new_size++] = b;
400 }
401 if (new_size == 0) return {};
402 std::shuffle(&boxes[0], &boxes[0] + new_size, random);
403 return {&boxes[0], new_size};
404}
405
406absl::Span<int> FilterBoxesThatAreTooLarge(
407 absl::Span<const Rectangle> cached_rectangles,
408 absl::Span<const IntegerValue> energies, absl::Span<int> boxes) {
409 // Sort the boxes by increasing area.
410 std::sort(boxes.begin(), boxes.end(), [cached_rectangles](int a, int b) {
411 return cached_rectangles[a].Area() < cached_rectangles[b].Area();
412 });
413
414 IntegerValue total_energy(0);
415 for (const int box : boxes) total_energy += energies[box];
416
417 // Remove all the large boxes until we have one with area smaller than the
418 // energy of the boxes below.
419 int new_size = boxes.size();
420 while (new_size > 0 &&
421 cached_rectangles[boxes[new_size - 1]].Area() >= total_energy) {
422 --new_size;
423 total_energy -= energies[boxes[new_size]];
424 }
425 return boxes.subspan(0, new_size);
426}
427
428void ConstructOverlappingSets(absl::Span<IndexedInterval> intervals,
429 CompactVectorVector<int>* result,
430 absl::Span<const int> order) {
431 result->clear();
432 DCHECK(std::is_sorted(intervals.begin(), intervals.end(),
434 IntegerValue min_end_in_set = kMaxIntegerValue;
435
436 // We do a line sweep. The "current" subset crossing the "line" at
437 // (time, time + 1) will be in (*intervals)[start_index, end_index) at the end
438 // of the loop block.
439 int start_index = 0;
440 const int size = intervals.size();
441 for (int end_index = 0; end_index < size;) {
442 const IntegerValue time = intervals[end_index].start;
443
444 // First, if there is some deletion, we will push the "old" set to the
445 // result before updating it. Otherwise, we will have a superset later, so
446 // we just continue for now.
447 if (min_end_in_set <= time) {
448 // Push the current set to result first if its size is > 1.
449 if (start_index + 1 < end_index) {
450 result->Add({});
451 for (int i = start_index; i < end_index; ++i) {
452 result->AppendToLastVector(intervals[i].index);
453 }
454 }
455
456 // Update the set. Note that we keep the order.
457 min_end_in_set = kMaxIntegerValue;
458 int new_start = end_index;
459 for (int i = end_index; --i >= start_index;) {
460 if (intervals[i].end > time) {
461 min_end_in_set = std::min(min_end_in_set, intervals[i].end);
462 intervals[--new_start] = intervals[i];
463 }
464 }
465 start_index = new_start;
466 }
467
468 // Add all the new intervals starting exactly at "time".
469 // Note that we always add at least one here.
470 const int old_end = end_index;
471 while (end_index < size && intervals[end_index].start == time) {
472 min_end_in_set = std::min(min_end_in_set, intervals[end_index].end);
473 ++end_index;
474 }
475
476 // If order is not empty, make sure we maintain the order.
477 // TODO(user): we could only do that when we push a new set.
478 if (!order.empty() && end_index > old_end) {
479 std::sort(intervals.data() + old_end, intervals.data() + end_index,
480 [order](const IndexedInterval& a, const IndexedInterval& b) {
481 return order[a.index] < order[b.index];
482 });
483 std::inplace_merge(
484 intervals.data() + start_index, intervals.data() + old_end,
485 intervals.data() + end_index,
486 [order](const IndexedInterval& a, const IndexedInterval& b) {
487 return order[a.index] < order[b.index];
488 });
489 }
490 }
491
492 // Push final set.
493 if (start_index + 1 < size) {
494 result->Add({});
495 for (int i = start_index; i < size; ++i) {
496 result->AppendToLastVector(intervals[i].index);
497 }
498 }
499}
500
502 std::vector<IndexedInterval>* intervals,
503 std::vector<std::vector<int>>* components) {
504 components->clear();
505 if (intervals->empty()) return;
506 if (intervals->size() == 1) {
507 components->push_back({intervals->front().index});
508 return;
509 }
510
511 // For correctness, ComparatorByStart is enough, but in unit tests we want to
512 // verify this function against another implementation, and fully defined
513 // sorting with tie-breaking makes that much easier.
514 // If that becomes a performance bottleneck:
515 // - One may want to sort the list outside of this function, and simply
516 // have this function DCHECK that it's sorted by start.
517 // - One may use stable_sort() with ComparatorByStart().
518 std::sort(intervals->begin(), intervals->end(),
520
521 IntegerValue end_max_so_far = (*intervals)[0].end;
522 components->push_back({(*intervals)[0].index});
523 for (int i = 1; i < intervals->size(); ++i) {
524 const IndexedInterval& interval = (*intervals)[i];
525 if (interval.start >= end_max_so_far) {
526 components->push_back({interval.index});
527 } else {
528 components->back().push_back(interval.index);
529 }
530 end_max_so_far = std::max(end_max_so_far, interval.end);
531 }
532}
533
534std::vector<int> GetIntervalArticulationPoints(
535 std::vector<IndexedInterval>* intervals) {
536 std::vector<int> articulation_points;
537 if (intervals->size() < 3) return articulation_points; // Empty.
538 if (DEBUG_MODE) {
539 for (const IndexedInterval& interval : *intervals) {
540 DCHECK_LT(interval.start, interval.end);
541 }
542 }
543
544 std::sort(intervals->begin(), intervals->end(),
545 IndexedInterval::ComparatorByStart());
546
547 IntegerValue end_max_so_far = (*intervals)[0].end;
548 int index_of_max = 0;
549 IntegerValue prev_end_max = kMinIntegerValue; // Initialized as a sentinel.
550 for (int i = 1; i < intervals->size(); ++i) {
551 const IndexedInterval& interval = (*intervals)[i];
552 if (interval.start >= end_max_so_far) {
553 // New connected component.
554 end_max_so_far = interval.end;
555 index_of_max = i;
556 prev_end_max = kMinIntegerValue;
557 continue;
558 }
559 // Still the same connected component. Was the previous "max" an
560 // articulation point ?
561 if (prev_end_max != kMinIntegerValue && interval.start >= prev_end_max) {
562 // We might be re-inserting the same articulation point: guard against it.
563 if (articulation_points.empty() ||
564 articulation_points.back() != index_of_max) {
565 articulation_points.push_back(index_of_max);
566 }
567 }
568 // Update the max end.
569 if (interval.end > end_max_so_far) {
570 prev_end_max = end_max_so_far;
571 end_max_so_far = interval.end;
572 index_of_max = i;
573 } else if (interval.end > prev_end_max) {
574 prev_end_max = interval.end;
575 }
576 }
577 // Convert articulation point indices to IndexedInterval.index.
578 for (int& index : articulation_points) index = (*intervals)[index].index;
579 return articulation_points;
580}
581
582namespace {
583bool IsZeroOrPowerOfTwo(int value) { return (value & (value - 1)) == 0; }
584
585void AppendPairwiseRestriction(const ItemWithVariableSize& item1,
586 const ItemWithVariableSize& item2,
587 std::vector<PairwiseRestriction>* result) {
588 const int state =
589 // box1 can be left of box2.
590 (item1.x.end_min <= item2.x.start_max) +
591 // box1 can be right of box2.
592 2 * (item2.x.end_min <= item1.x.start_max) +
593 // box1 can be below box2.
594 4 * (item1.y.end_min <= item2.y.start_max) +
595 // box1 can be up of box2.
596 8 * (item2.y.end_min <= item1.y.start_max);
597
598 if (!IsZeroOrPowerOfTwo(state)) {
599 return;
600 }
601
603 switch (state) {
604 case 0:
605 // Conflict. The two boxes must overlap in both dimensions.
607 break;
608 case 1:
609 // box2 can only be after box1 on x.
610 if (item1.x.end_min > item2.x.start_min ||
611 item2.x.start_max < item1.x.end_max) {
612 type =
614 break;
615 }
616 return;
617 case 2:
618 // box1 can only be after box2 on x.
619 if (item2.x.end_min > item1.x.start_min ||
620 item1.x.start_max < item2.x.end_max) {
621 type =
623 break;
624 }
625 return;
626 case 4:
627 // box2 can only be after box1 on y.
628 if (item1.y.end_min > item2.y.start_min ||
629 item2.y.start_max < item1.y.end_max) {
631 break;
632 }
633 return;
634 case 8:
635 // box1 can only be after box2 on y.
636 if (item2.y.end_min > item1.y.start_min ||
637 item1.y.start_max < item2.y.end_max) {
639 break;
640 }
641 return;
642 default:
643 return;
644 }
645 result->push_back(
646 {.first_index = item1.index, .second_index = item2.index, .type = type});
647}
648} // namespace
649
650void AppendPairwiseRestrictions(absl::Span<const ItemWithVariableSize> items,
651 std::vector<PairwiseRestriction>* result) {
652 for (int i1 = 0; i1 + 1 < items.size(); ++i1) {
653 for (int i2 = i1 + 1; i2 < items.size(); ++i2) {
654 AppendPairwiseRestriction(items[i1], items[i2], result);
655 }
656 }
657}
658
660 absl::Span<const ItemWithVariableSize> items,
661 absl::Span<const ItemWithVariableSize> other_items,
662 std::vector<PairwiseRestriction>* result) {
663 for (int i1 = 0; i1 < items.size(); ++i1) {
664 for (int i2 = 0; i2 < other_items.size(); ++i2) {
665 AppendPairwiseRestriction(items[i1], other_items[i2], result);
666 }
667 }
668}
669
671 events_.clear();
672 num_rectangles_added_ = 0;
673}
674
675void CapacityProfile::AddRectangle(IntegerValue x_min, IntegerValue x_max,
676 IntegerValue y_min, IntegerValue y_max) {
677 DCHECK_LE(x_min, x_max);
678 if (x_min == x_max) return;
679
680 events_.push_back(
681 StartRectangleEvent(num_rectangles_added_, x_min, y_min, y_max));
682 events_.push_back(EndRectangleEvent(num_rectangles_added_, x_max));
683 ++num_rectangles_added_;
684}
685
686void CapacityProfile::AddMandatoryConsumption(IntegerValue x_min,
687 IntegerValue x_max,
688 IntegerValue y_height) {
689 DCHECK_LE(x_min, x_max);
690 if (x_min == x_max) return;
691
692 events_.push_back(ChangeMandatoryProfileEvent(x_min, y_height));
693 events_.push_back(ChangeMandatoryProfileEvent(x_max, -y_height));
694}
695
697 std::vector<CapacityProfile::Rectangle>* result) {
698 std::sort(events_.begin(), events_.end());
699 IntegerPriorityQueue<QueueElement> min_pq(num_rectangles_added_);
700 IntegerPriorityQueue<QueueElement> max_pq(num_rectangles_added_);
701 IntegerValue mandatory_capacity(0);
702
703 result->clear();
704
705 result->push_back({kMinIntegerValue, IntegerValue(0)});
706
707 for (int i = 0; i < events_.size();) {
708 const IntegerValue current_time = events_[i].time;
709 for (; i < events_.size(); ++i) {
710 const Event& event = events_[i];
711 if (event.time != current_time) break;
712
713 switch (events_[i].type) {
714 case START_RECTANGLE: {
715 min_pq.Add({event.index, -event.y_min});
716 max_pq.Add({event.index, event.y_max});
717 break;
718 }
719 case END_RECTANGLE: {
720 min_pq.Remove(event.index);
721 max_pq.Remove(event.index);
722 break;
723 }
724 case CHANGE_MANDATORY_PROFILE: {
725 mandatory_capacity += event.y_min;
726 break;
727 }
728 }
729 }
730
731 DCHECK(!max_pq.IsEmpty() || mandatory_capacity == 0);
732 const IntegerValue new_height =
733 max_pq.IsEmpty()
734 ? IntegerValue(0)
735 : max_pq.Top().value + min_pq.Top().value - mandatory_capacity;
736 if (new_height != result->back().height) {
737 result->push_back({current_time, new_height});
738 }
739 }
740}
741
742IntegerValue CapacityProfile::GetBoundingArea() {
743 std::sort(events_.begin(), events_.end());
744 IntegerPriorityQueue<QueueElement> min_pq(num_rectangles_added_);
745 IntegerPriorityQueue<QueueElement> max_pq(num_rectangles_added_);
746
747 IntegerValue area(0);
748 IntegerValue previous_time = kMinIntegerValue;
749 IntegerValue previous_height(0);
750
751 for (int i = 0; i < events_.size();) {
752 const IntegerValue current_time = events_[i].time;
753 for (; i < events_.size(); ++i) {
754 const Event& event = events_[i];
755 if (event.time != current_time) break;
756
757 switch (event.type) {
758 case START_RECTANGLE: {
759 min_pq.Add({event.index, -event.y_min});
760 max_pq.Add({event.index, event.y_max});
761 break;
762 }
763 case END_RECTANGLE: {
764 min_pq.Remove(event.index);
765 max_pq.Remove(event.index);
766 break;
767 }
768 case CHANGE_MANDATORY_PROFILE: {
769 break;
770 }
771 }
772 }
773 const IntegerValue new_height =
774 max_pq.IsEmpty() ? IntegerValue(0)
775 : max_pq.Top().value + min_pq.Top().value;
776 if (previous_height != 0) {
777 area += previous_height * (current_time - previous_time);
778 }
779 previous_time = current_time;
780 previous_height = new_height;
781 }
782 return area;
783}
784
785IntegerValue Smallest1DIntersection(IntegerValue range_min,
786 IntegerValue range_max, IntegerValue size,
787 IntegerValue interval_min,
788 IntegerValue interval_max) {
789 // If the item is on the left of the range, we get the intersection between
790 // [range_min, range_min + size] and [interval_min, interval_max].
791 const IntegerValue overlap_on_left =
792 std::min(range_min + size, interval_max) -
793 std::max(range_min, interval_min);
794
795 // If the item is on the right of the range, we get the intersection between
796 // [range_max - size, range_max] and [interval_min, interval_max].
797 const IntegerValue overlap_on_right =
798 std::min(range_max, interval_max) -
799 std::max(range_max - size, interval_min);
800
801 return std::max(IntegerValue(0), std::min(overlap_on_left, overlap_on_right));
802}
803
804ProbingRectangle::ProbingRectangle(
805 const std::vector<RectangleInRange>& intervals)
806 : intervals_(intervals) {
807 minimum_energy_ = 0;
808 if (intervals_.empty()) {
809 return;
810 }
811 interval_points_sorted_by_x_.reserve(intervals_.size() * 4 + 2);
812 interval_points_sorted_by_y_.reserve(intervals_.size() * 4 + 2);
813
814 Rectangle bounding_box = {.x_min = std::numeric_limits<IntegerValue>::max(),
815 .x_max = std::numeric_limits<IntegerValue>::min(),
816 .y_min = std::numeric_limits<IntegerValue>::max(),
817 .y_max = std::numeric_limits<IntegerValue>::min()};
818
819 for (int i = 0; i < intervals_.size(); ++i) {
820 const RectangleInRange& interval = intervals_[i];
821 minimum_energy_ += interval.x_size * interval.y_size;
822
823 bounding_box.x_min =
824 std::min(bounding_box.x_min, interval.bounding_area.x_min);
825 bounding_box.x_max =
826 std::max(bounding_box.x_max, interval.bounding_area.x_max);
827 bounding_box.y_min =
828 std::min(bounding_box.y_min, interval.bounding_area.y_min);
829 bounding_box.y_max =
830 std::max(bounding_box.y_max, interval.bounding_area.y_max);
831
832 interval_points_sorted_by_x_.push_back({interval.bounding_area.x_min, i});
833 interval_points_sorted_by_x_.push_back(
834 {interval.bounding_area.x_min + interval.x_size, i});
835 interval_points_sorted_by_x_.push_back(
836 {interval.bounding_area.x_max - interval.x_size, i});
837 interval_points_sorted_by_x_.push_back({interval.bounding_area.x_max, i});
838
839 interval_points_sorted_by_y_.push_back({interval.bounding_area.y_min, i});
840 interval_points_sorted_by_y_.push_back(
841 {interval.bounding_area.y_min + interval.y_size, i});
842 interval_points_sorted_by_y_.push_back(
843 {interval.bounding_area.y_max - interval.y_size, i});
844 interval_points_sorted_by_y_.push_back({interval.bounding_area.y_max, i});
845 }
846
847 full_energy_ = minimum_energy_;
848 // Add four bogus points in the extremities so we can delegate setting up all
849 // internal state to Shrink().
850 interval_points_sorted_by_x_.push_back({bounding_box.x_min - 1, -1});
851 interval_points_sorted_by_x_.push_back({bounding_box.x_max + 1, -1});
852 interval_points_sorted_by_y_.push_back({bounding_box.y_min - 1, -1});
853 interval_points_sorted_by_y_.push_back({bounding_box.y_max + 1, -1});
854
855 auto comparator = [](const IntervalPoint& a, const IntervalPoint& b) {
856 return std::tie(a.value, a.index) < std::tie(b.value, b.index);
857 };
858 gtl::STLSortAndRemoveDuplicates(&interval_points_sorted_by_x_, comparator);
859 gtl::STLSortAndRemoveDuplicates(&interval_points_sorted_by_y_, comparator);
860
861 grouped_intervals_sorted_by_x_.reserve(interval_points_sorted_by_x_.size());
862 grouped_intervals_sorted_by_y_.reserve(interval_points_sorted_by_y_.size());
863 int i = 0;
864 while (i < interval_points_sorted_by_x_.size()) {
865 int idx_begin = i;
866 while (i < interval_points_sorted_by_x_.size() &&
867 interval_points_sorted_by_x_[i].value ==
868 interval_points_sorted_by_x_[idx_begin].value) {
869 i++;
870 }
871 grouped_intervals_sorted_by_x_.push_back(
872 {interval_points_sorted_by_x_[idx_begin].value,
873 absl::Span<IntervalPoint>(interval_points_sorted_by_x_)
874 .subspan(idx_begin, i - idx_begin)});
875 }
876
877 i = 0;
878 while (i < interval_points_sorted_by_y_.size()) {
879 int idx_begin = i;
880 while (i < interval_points_sorted_by_y_.size() &&
881 interval_points_sorted_by_y_[i].value ==
882 interval_points_sorted_by_y_[idx_begin].value) {
883 i++;
884 }
885 grouped_intervals_sorted_by_y_.push_back(
886 {interval_points_sorted_by_y_[idx_begin].value,
887 absl::Span<IntervalPoint>(interval_points_sorted_by_y_)
888 .subspan(idx_begin, i - idx_begin)});
889 }
890
891 Reset();
892}
893
895 indexes_[Edge::LEFT] = 0;
896 indexes_[Edge::RIGHT] = grouped_intervals_sorted_by_x_.size() - 1;
897 indexes_[Edge::BOTTOM] = 0;
898 indexes_[Edge::TOP] = grouped_intervals_sorted_by_y_.size() - 1;
899
900 next_indexes_[Edge::LEFT] = 1;
901 next_indexes_[Edge::RIGHT] = grouped_intervals_sorted_by_x_.size() - 2;
902 next_indexes_[Edge::BOTTOM] = 1;
903 next_indexes_[Edge::TOP] = grouped_intervals_sorted_by_y_.size() - 2;
904
905 minimum_energy_ = full_energy_;
906 ranges_touching_both_boundaries_[0].clear();
907 ranges_touching_both_boundaries_[1].clear();
908
909 for (int i = 0; i < 4; ++i) {
910 corner_count_[i] = 0;
911 intersect_length_[i] = 0;
912 cached_delta_energy_[i] = 0;
913 }
914
915 // Remove the four bogus points we added.
916 Shrink(Edge::LEFT);
917 Shrink(Edge::BOTTOM);
918 Shrink(Edge::RIGHT);
919 Shrink(Edge::TOP);
920}
921
922Rectangle ProbingRectangle::GetCurrentRectangle() const {
923 return {
924 .x_min = grouped_intervals_sorted_by_x_[indexes_[Edge::LEFT]].coordinate,
925 .x_max = grouped_intervals_sorted_by_x_[indexes_[Edge::RIGHT]].coordinate,
926 .y_min =
927 grouped_intervals_sorted_by_y_[indexes_[Edge::BOTTOM]].coordinate,
928 .y_max = grouped_intervals_sorted_by_y_[indexes_[Edge::TOP]].coordinate};
929}
930
931namespace {
932// Intersects `rectangle` with the largest rectangle that must intersect with
933// the range in some way. To visualize this largest rectangle, imagine the four
934// possible extreme positions for the item in range (the four corners). This
935// rectangle is the one defined by the interior points of each position. This
936// don't use IsDisjoint() because it also works when the rectangle would be
937// malformed (it's bounding box less than twice the size).
938bool CanConsumeEnergy(const Rectangle& rectangle,
939 const RectangleInRange& item) {
940 return (rectangle.x_max > item.bounding_area.x_max - item.x_size) &&
941 (rectangle.y_max > item.bounding_area.y_max - item.y_size) &&
942 (rectangle.x_min < item.bounding_area.x_min + item.x_size) &&
943 (rectangle.y_min < item.bounding_area.y_min + item.y_size);
944}
945
946std::array<bool, 4> GetPossibleEdgeIntersection(const Rectangle& rectangle,
947 const RectangleInRange& range) {
948 std::array<bool, 4> result;
949 using Edge = ProbingRectangle::Edge;
950 result[Edge::LEFT] = rectangle.x_min >= range.bounding_area.x_min;
951 result[Edge::BOTTOM] = rectangle.y_min >= range.bounding_area.y_min;
952 result[Edge::RIGHT] = rectangle.x_max <= range.bounding_area.x_max;
953 result[Edge::TOP] = rectangle.y_max <= range.bounding_area.y_max;
954 return result;
955}
956
957} // namespace
958
959// NOMUTANTS -- This is a sanity check, it is hard to corrupt the data in an
960// unit test to check it will fail.
961void ProbingRectangle::ValidateInvariants() const {
962 const Rectangle current_rectangle = GetCurrentRectangle();
964 IntegerValue intersect_length[4] = {0, 0, 0, 0};
965 IntegerValue corner_count[4] = {0, 0, 0, 0};
966 IntegerValue energy = 0;
967 CHECK_LE(next_indexes_[Edge::LEFT], indexes_[Edge::RIGHT]);
968 CHECK_LE(next_indexes_[Edge::BOTTOM], indexes_[Edge::TOP]);
969 CHECK_GE(next_indexes_[Edge::TOP], indexes_[Edge::BOTTOM]);
970 CHECK_GE(next_indexes_[Edge::RIGHT], indexes_[Edge::LEFT]);
971
972 for (int interval_idx = 0; interval_idx < intervals_.size(); interval_idx++) {
973 const RectangleInRange& range = intervals_[interval_idx];
974
975 const Rectangle min_intersect =
976 range.GetMinimumIntersection(current_rectangle);
977 CHECK_LE(min_intersect.SizeX(), range.x_size);
978 CHECK_LE(min_intersect.SizeY(), range.y_size);
979 energy += min_intersect.Area();
980
981 std::array<bool, 4> touching_boundary = {false, false, false, false};
982 CHECK_EQ(CanConsumeEnergy(current_rectangle, range) &&
983 current_rectangle.Area() != 0,
984 range.GetMinimumIntersectionArea(current_rectangle) != 0);
985 if (CanConsumeEnergy(current_rectangle, range)) {
986 touching_boundary = GetPossibleEdgeIntersection(current_rectangle, range);
987 }
988
989 CHECK_EQ(
990 touching_boundary[Edge::LEFT] && touching_boundary[Edge::RIGHT],
991 ranges_touching_both_boundaries_[Direction::LEFT_AND_RIGHT].contains(
992 interval_idx));
993 CHECK_EQ(
994 touching_boundary[Edge::TOP] && touching_boundary[Edge::BOTTOM],
995 ranges_touching_both_boundaries_[Direction::TOP_AND_BOTTOM].contains(
996 interval_idx));
997
998 if (touching_boundary[Edge::LEFT] && !touching_boundary[Edge::RIGHT]) {
999 intersect_length[Edge::LEFT] += Smallest1DIntersection(
1000 range.bounding_area.y_min, range.bounding_area.y_max, range.y_size,
1001 current_rectangle.y_min, current_rectangle.y_max);
1002 }
1003
1004 if (touching_boundary[Edge::RIGHT] && !touching_boundary[Edge::LEFT]) {
1005 intersect_length[Edge::RIGHT] += Smallest1DIntersection(
1006 range.bounding_area.y_min, range.bounding_area.y_max, range.y_size,
1007 current_rectangle.y_min, current_rectangle.y_max);
1008 }
1009
1010 if (touching_boundary[Edge::TOP] && !touching_boundary[Edge::BOTTOM]) {
1011 intersect_length[Edge::TOP] += Smallest1DIntersection(
1012 range.bounding_area.x_min, range.bounding_area.x_max, range.x_size,
1013 current_rectangle.x_min, current_rectangle.x_max);
1014 }
1015
1016 if (touching_boundary[Edge::BOTTOM] && !touching_boundary[Edge::TOP]) {
1017 intersect_length[Edge::BOTTOM] += Smallest1DIntersection(
1018 range.bounding_area.x_min, range.bounding_area.x_max, range.x_size,
1019 current_rectangle.x_min, current_rectangle.x_max);
1020 }
1021
1022 if ((touching_boundary[Edge::LEFT] && touching_boundary[Edge::RIGHT]) ||
1023 (touching_boundary[Edge::TOP] && touching_boundary[Edge::BOTTOM])) {
1024 // We account separately for the problematic items that touches both
1025 // sides.
1026 continue;
1027 }
1028 if (touching_boundary[Edge::BOTTOM] && touching_boundary[Edge::LEFT]) {
1029 corner_count[RectangleInRange::Corner::BOTTOM_LEFT]++;
1030 }
1031 if (touching_boundary[Edge::BOTTOM] && touching_boundary[Edge::RIGHT]) {
1032 corner_count[RectangleInRange::Corner::BOTTOM_RIGHT]++;
1033 }
1034 if (touching_boundary[Edge::TOP] && touching_boundary[Edge::LEFT]) {
1035 corner_count[RectangleInRange::Corner::TOP_LEFT]++;
1036 }
1037 if (touching_boundary[Edge::TOP] && touching_boundary[Edge::RIGHT]) {
1038 corner_count[RectangleInRange::Corner::TOP_RIGHT]++;
1039 }
1040 }
1041
1042 CHECK_EQ(energy, minimum_energy_);
1043 for (int i = 0; i < 4; i++) {
1044 CHECK_EQ(intersect_length[i], intersect_length_[i]);
1045 CHECK_EQ(corner_count[i], corner_count_[i]);
1046 }
1047}
1048
1049namespace {
1050
1051struct EdgeInfo {
1052 using Edge = ProbingRectangle::Edge;
1053 using Direction = ProbingRectangle::Direction;
1054 using Corner = RectangleInRange::Corner;
1055
1056 Edge opposite_edge;
1057
1058 struct OrthogonalInfo {
1059 Edge edge;
1060 Corner adjacent_corner;
1061 };
1062
1063 Direction shrink_direction;
1064 Direction orthogonal_shrink_direction;
1065 // Lower coordinate one first (ie., BOTTOM before TOP, LEFT before RIGHT).
1066 OrthogonalInfo orthogonal_edges[2];
1067};
1068
1069struct EdgeInfoHolder {
1070 using Edge = ProbingRectangle::Edge;
1071 using Direction = ProbingRectangle::Direction;
1072 using Corner = RectangleInRange::Corner;
1073
1074 static constexpr EdgeInfo kLeft = {
1075 .opposite_edge = Edge::RIGHT,
1076 .shrink_direction = Direction::LEFT_AND_RIGHT,
1077 .orthogonal_shrink_direction = Direction::TOP_AND_BOTTOM,
1078 .orthogonal_edges = {
1079 {.edge = Edge::BOTTOM, .adjacent_corner = Corner::BOTTOM_LEFT},
1080 {.edge = Edge::TOP, .adjacent_corner = Corner::TOP_LEFT}}};
1081
1082 static constexpr EdgeInfo kRight = {
1083 .opposite_edge = Edge::LEFT,
1084 .shrink_direction = Direction::LEFT_AND_RIGHT,
1085 .orthogonal_shrink_direction = Direction::TOP_AND_BOTTOM,
1086 .orthogonal_edges = {
1087 {.edge = Edge::BOTTOM, .adjacent_corner = Corner::BOTTOM_RIGHT},
1088 {.edge = Edge::TOP, .adjacent_corner = Corner::TOP_RIGHT}}};
1089 static constexpr EdgeInfo kBottom = {
1090 .opposite_edge = Edge::TOP,
1091 .shrink_direction = Direction::TOP_AND_BOTTOM,
1092 .orthogonal_shrink_direction = Direction::LEFT_AND_RIGHT,
1093 .orthogonal_edges = {
1094 {.edge = Edge::LEFT, .adjacent_corner = Corner::BOTTOM_LEFT},
1095 {.edge = Edge::RIGHT, .adjacent_corner = Corner::BOTTOM_RIGHT}}};
1096 static constexpr EdgeInfo kTop = {
1097 .opposite_edge = Edge::BOTTOM,
1098 .shrink_direction = Direction::TOP_AND_BOTTOM,
1099 .orthogonal_shrink_direction = Direction::LEFT_AND_RIGHT,
1100 .orthogonal_edges = {
1101 {.edge = Edge::LEFT, .adjacent_corner = Corner::TOP_LEFT},
1102 {.edge = Edge::RIGHT, .adjacent_corner = Corner::TOP_RIGHT}}};
1103};
1104
1105constexpr const EdgeInfo& GetEdgeInfo(ProbingRectangle::Edge edge) {
1106 using Edge = ProbingRectangle::Edge;
1107 switch (edge) {
1108 case Edge::LEFT:
1109 return EdgeInfoHolder::kLeft;
1110 case Edge::RIGHT:
1111 return EdgeInfoHolder::kRight;
1112 case Edge::BOTTOM:
1113 return EdgeInfoHolder::kBottom;
1114 case Edge::TOP:
1115 return EdgeInfoHolder::kTop;
1116 }
1117}
1118
1119IntegerValue GetSmallest1DIntersection(ProbingRectangle::Direction direction,
1120 const RectangleInRange& range,
1121 const Rectangle& rectangle) {
1122 switch (direction) {
1123 case ProbingRectangle::Direction::LEFT_AND_RIGHT:
1124 return Smallest1DIntersection(range.bounding_area.x_min,
1125 range.bounding_area.x_max, range.x_size,
1126 rectangle.x_min, rectangle.x_max);
1127 case ProbingRectangle::Direction::TOP_AND_BOTTOM:
1128 return Smallest1DIntersection(range.bounding_area.y_min,
1129 range.bounding_area.y_max, range.y_size,
1130 rectangle.y_min, rectangle.y_max);
1131 }
1132}
1133
1134} // namespace
1135
1136template <ProbingRectangle::Edge edge>
1137void ProbingRectangle::ShrinkImpl() {
1138 constexpr EdgeInfo e = GetEdgeInfo(edge);
1139
1140 bool update_next_index[4] = {false, false, false, false};
1141 update_next_index[edge] = true;
1142
1143 IntegerValue step_1d_size;
1144 minimum_energy_ -= GetShrinkDeltaEnergy(edge);
1145 const std::vector<PointsForCoordinate>& sorted_intervals =
1146 e.shrink_direction == Direction::LEFT_AND_RIGHT
1147 ? grouped_intervals_sorted_by_x_
1148 : grouped_intervals_sorted_by_y_;
1149
1150 const Rectangle prev_rectangle = GetCurrentRectangle();
1151 indexes_[edge] = next_indexes_[edge];
1152 const Rectangle current_rectangle = GetCurrentRectangle();
1153
1154 switch (edge) {
1155 case Edge::LEFT:
1156 step_1d_size = current_rectangle.x_min - prev_rectangle.x_min;
1157 next_indexes_[edge] =
1158 std::min(indexes_[edge] + 1, indexes_[e.opposite_edge]);
1159 next_indexes_[e.opposite_edge] =
1160 std::max(indexes_[edge], next_indexes_[e.opposite_edge]);
1161 break;
1162 case Edge::BOTTOM:
1163 step_1d_size = current_rectangle.y_min - prev_rectangle.y_min;
1164 next_indexes_[edge] =
1165 std::min(indexes_[edge] + 1, indexes_[e.opposite_edge]);
1166 next_indexes_[e.opposite_edge] =
1167 std::max(indexes_[edge], next_indexes_[e.opposite_edge]);
1168 break;
1169 case Edge::RIGHT:
1170 step_1d_size = prev_rectangle.x_max - current_rectangle.x_max;
1171 next_indexes_[edge] =
1172 std::max(indexes_[edge] - 1, indexes_[e.opposite_edge]);
1173 next_indexes_[e.opposite_edge] =
1174 std::min(indexes_[edge], next_indexes_[e.opposite_edge]);
1175 break;
1176 case Edge::TOP:
1177 step_1d_size = prev_rectangle.y_max - current_rectangle.y_max;
1178 next_indexes_[edge] =
1179 std::max(indexes_[edge] - 1, indexes_[e.opposite_edge]);
1180 next_indexes_[e.opposite_edge] =
1181 std::min(indexes_[edge], next_indexes_[e.opposite_edge]);
1182 break;
1183 }
1184
1185 absl::Span<ProbingRectangle::IntervalPoint> items_touching_coordinate =
1186 sorted_intervals[indexes_[edge]].items_touching_coordinate;
1187
1188 IntegerValue delta_corner_count[4] = {0, 0, 0, 0};
1189 for (const auto& item : items_touching_coordinate) {
1190 const RectangleInRange& range = intervals_[item.index];
1191 if (!CanConsumeEnergy(prev_rectangle, range)) {
1192 // This item is out of our area of interest, skip.
1193 continue;
1194 }
1195
1196 const std::array<bool, 4> touching_boundary_before =
1197 GetPossibleEdgeIntersection(prev_rectangle, range);
1198 const std::array<bool, 4> touching_boundary_after =
1199 CanConsumeEnergy(current_rectangle, range)
1200 ? GetPossibleEdgeIntersection(current_rectangle, range)
1201 : std::array<bool, 4>({false, false, false, false});
1202
1203 bool remove_corner[4] = {false, false, false, false};
1204 auto erase_item = [this, &prev_rectangle, &range, &touching_boundary_before,
1205 &remove_corner](Edge edge_to_erase) {
1206 const EdgeInfo& erase_info = GetEdgeInfo(edge_to_erase);
1207 intersect_length_[edge_to_erase] -= GetSmallest1DIntersection(
1208 erase_info.orthogonal_shrink_direction, range, prev_rectangle);
1209
1210 if (touching_boundary_before[erase_info.orthogonal_edges[0].edge] &&
1211 touching_boundary_before[erase_info.orthogonal_edges[1].edge]) {
1212 // Ignore touching both corners
1213 return;
1214 }
1215 for (const auto [orthogonal_edge, corner] : erase_info.orthogonal_edges) {
1216 if (touching_boundary_before[orthogonal_edge]) {
1217 remove_corner[corner] = true;
1218 }
1219 }
1220 };
1221
1222 if (touching_boundary_after[edge] && !touching_boundary_before[edge]) {
1223 if (touching_boundary_before[e.opposite_edge]) {
1224 ranges_touching_both_boundaries_[e.shrink_direction].insert(item.index);
1225 erase_item(e.opposite_edge);
1226 } else {
1227 // Do the opposite of remove_item().
1228 intersect_length_[edge] += GetSmallest1DIntersection(
1229 e.orthogonal_shrink_direction, range, prev_rectangle);
1230 // Update the corner count unless it is touching both.
1231 if (!touching_boundary_before[e.orthogonal_edges[0].edge] ||
1232 !touching_boundary_before[e.orthogonal_edges[1].edge]) {
1233 for (const auto [orthogonal_edge, corner] : e.orthogonal_edges) {
1234 if (touching_boundary_before[orthogonal_edge]) {
1235 delta_corner_count[corner]++;
1236 }
1237 }
1238 }
1239 }
1240 }
1241
1242 for (int i = 0; i < 4; i++) {
1243 const Edge edge_to_update = (Edge)i;
1244 const EdgeInfo& info = GetEdgeInfo(edge_to_update);
1245 const bool remove_edge = touching_boundary_before[edge_to_update] &&
1246 !touching_boundary_after[edge_to_update];
1247 if (!remove_edge) {
1248 continue;
1249 }
1250
1251 update_next_index[edge_to_update] = true;
1252
1253 if (touching_boundary_before[info.opposite_edge]) {
1254 ranges_touching_both_boundaries_[info.shrink_direction].erase(
1255 item.index);
1256 } else {
1257 erase_item(edge_to_update);
1258 }
1259 }
1260
1261 for (int i = 0; i < 4; i++) {
1262 corner_count_[i] -= remove_corner[i];
1263 }
1264 }
1265
1266 // Update the intersection length for items touching both sides.
1267 for (const int idx : ranges_touching_both_boundaries_[e.shrink_direction]) {
1268 const RectangleInRange& range = intervals_[idx];
1269 const std::array<bool, 2> touching_corner =
1270 (e.shrink_direction == Direction::LEFT_AND_RIGHT)
1271 ? std::array<bool, 2>(
1272 {current_rectangle.y_min >= range.bounding_area.y_min,
1273 current_rectangle.y_max <= range.bounding_area.y_max})
1274 : std::array<bool, 2>(
1275 {current_rectangle.x_min >= range.bounding_area.x_min,
1276 current_rectangle.x_max <= range.bounding_area.x_max});
1277 if (touching_corner[0] == touching_corner[1]) {
1278 // Either it is not touching neither corners (so no length to update) or
1279 // it is touching both corners, which will be handled by the "both
1280 // sides" code and should not contribute to intersect_length_.
1281 continue;
1282 }
1283
1284 const IntegerValue incr =
1285 GetSmallest1DIntersection(e.shrink_direction, range, prev_rectangle) -
1286 GetSmallest1DIntersection(e.shrink_direction, range, current_rectangle);
1287 for (int i = 0; i < 2; i++) {
1288 if (touching_corner[i]) {
1289 intersect_length_[e.orthogonal_edges[i].edge] -= incr;
1290 }
1291 }
1292 }
1293
1294 for (const auto [orthogonal_edge, corner] : e.orthogonal_edges) {
1295 intersect_length_[orthogonal_edge] -= corner_count_[corner] * step_1d_size;
1296 }
1297
1298 for (int i = 0; i < 4; i++) {
1299 corner_count_[i] += delta_corner_count[i];
1300 }
1301
1302 auto points_consume_energy =
1303 [this,
1304 &current_rectangle](absl::Span<ProbingRectangle::IntervalPoint> points) {
1305 for (const auto& item : points) {
1306 const RectangleInRange& range = intervals_[item.index];
1307 if (CanConsumeEnergy(current_rectangle, range)) {
1308 return true;
1309 }
1310 }
1311 return false;
1312 };
1313 if (update_next_index[Edge::LEFT]) {
1314 for (; next_indexes_[Edge::LEFT] < indexes_[Edge::RIGHT];
1315 ++next_indexes_[Edge::LEFT]) {
1316 if (points_consume_energy(
1317 grouped_intervals_sorted_by_x_[next_indexes_[Edge::LEFT]]
1318 .items_touching_coordinate)) {
1319 break;
1320 }
1321 }
1322 }
1323 if (update_next_index[Edge::BOTTOM]) {
1324 for (; next_indexes_[Edge::BOTTOM] < indexes_[Edge::TOP];
1325 ++next_indexes_[Edge::BOTTOM]) {
1326 if (points_consume_energy(
1327 grouped_intervals_sorted_by_y_[next_indexes_[Edge::BOTTOM]]
1328 .items_touching_coordinate)) {
1329 break;
1330 }
1331 }
1332 }
1333 if (update_next_index[Edge::RIGHT]) {
1334 for (; next_indexes_[Edge::RIGHT] > indexes_[Edge::LEFT];
1335 --next_indexes_[Edge::RIGHT]) {
1336 if (points_consume_energy(
1337 grouped_intervals_sorted_by_x_[next_indexes_[Edge::RIGHT]]
1338 .items_touching_coordinate)) {
1339 break;
1340 }
1341 }
1342 }
1343 if (update_next_index[Edge::TOP]) {
1344 for (; next_indexes_[Edge::TOP] > indexes_[Edge::BOTTOM];
1345 --next_indexes_[Edge::TOP]) {
1346 if (points_consume_energy(
1347 grouped_intervals_sorted_by_y_[next_indexes_[Edge::TOP]]
1348 .items_touching_coordinate)) {
1349 break;
1350 }
1351 }
1352 }
1353
1354 probe_area_ = current_rectangle.Area();
1355 CacheShrinkDeltaEnergy(0);
1356 CacheShrinkDeltaEnergy(1);
1357}
1358
1359void ProbingRectangle::Shrink(Edge edge) {
1360 switch (edge) {
1362 ShrinkImpl<Edge::LEFT>();
1363 return;
1364 case Edge::BOTTOM:
1365 ShrinkImpl<Edge::BOTTOM>();
1366 return;
1367 case Edge::RIGHT:
1368 ShrinkImpl<Edge::RIGHT>();
1369 return;
1370 case Edge::TOP:
1371 ShrinkImpl<Edge::TOP>();
1372 return;
1373 }
1374}
1375
1376IntegerValue ProbingRectangle::GetShrinkDeltaArea(Edge edge) const {
1377 const Rectangle current_rectangle = GetCurrentRectangle();
1378 const std::vector<PointsForCoordinate>& sorted_intervals =
1379 (edge == Edge::LEFT || edge == Edge::RIGHT)
1380 ? grouped_intervals_sorted_by_x_
1381 : grouped_intervals_sorted_by_y_;
1382 const IntegerValue coordinate =
1383 sorted_intervals[next_indexes_[edge]].coordinate;
1384 switch (edge) {
1385 case Edge::LEFT:
1386 return (coordinate - current_rectangle.x_min) * current_rectangle.SizeY();
1387 case Edge::BOTTOM:
1388 return (coordinate - current_rectangle.y_min) * current_rectangle.SizeX();
1389 case Edge::RIGHT:
1390 return (current_rectangle.x_max - coordinate) * current_rectangle.SizeY();
1391 case Edge::TOP:
1392 return (current_rectangle.y_max - coordinate) * current_rectangle.SizeX();
1393 }
1394}
1395
1396void ProbingRectangle::CacheShrinkDeltaEnergy(int dimension) {
1397 const Rectangle current_rectangle = GetCurrentRectangle();
1398 Rectangle next_rectangle_up = current_rectangle;
1399 Rectangle next_rectangle_down = current_rectangle;
1400 IntegerValue step_1d_size_up, step_1d_size_down;
1401 IntegerValue units_crossed_up, units_crossed_down;
1402 IntegerValue* delta_energy_up_ptr;
1403 IntegerValue* delta_energy_down_ptr;
1404
1405 if (dimension == 0) {
1406 // CanShrink(Edge::RIGHT) and CanShrink(Edge::LEFT) are equivalent
1407 if (!CanShrink(Edge::LEFT)) {
1408 cached_delta_energy_[Edge::LEFT] = 0;
1409 cached_delta_energy_[Edge::RIGHT] = 0;
1410 return;
1411 }
1412
1413 next_rectangle_up.x_min =
1414 grouped_intervals_sorted_by_x_[next_indexes_[Edge::LEFT]].coordinate;
1415 next_rectangle_down.x_max =
1416 grouped_intervals_sorted_by_x_[next_indexes_[Edge::RIGHT]].coordinate;
1417
1418 step_1d_size_up = next_rectangle_up.x_min - current_rectangle.x_min;
1419 step_1d_size_down = current_rectangle.x_max - next_rectangle_down.x_max;
1420 units_crossed_up = intersect_length_[Edge::LEFT];
1421 units_crossed_down = intersect_length_[Edge::RIGHT];
1422 delta_energy_up_ptr = &cached_delta_energy_[Edge::LEFT];
1423 delta_energy_down_ptr = &cached_delta_energy_[Edge::RIGHT];
1424 } else {
1425 if (!CanShrink(Edge::TOP)) {
1426 cached_delta_energy_[Edge::BOTTOM] = 0;
1427 cached_delta_energy_[Edge::TOP] = 0;
1428 return;
1429 }
1430
1431 next_rectangle_up.y_min =
1432 grouped_intervals_sorted_by_y_[next_indexes_[Edge::BOTTOM]].coordinate;
1433 next_rectangle_down.y_max =
1434 grouped_intervals_sorted_by_y_[next_indexes_[Edge::TOP]].coordinate;
1435
1436 step_1d_size_up = next_rectangle_up.y_min - current_rectangle.y_min;
1437 step_1d_size_down = current_rectangle.y_max - next_rectangle_down.y_max;
1438 units_crossed_up = intersect_length_[Edge::BOTTOM];
1439 units_crossed_down = intersect_length_[Edge::TOP];
1440 delta_energy_up_ptr = &cached_delta_energy_[Edge::BOTTOM];
1441 delta_energy_down_ptr = &cached_delta_energy_[Edge::TOP];
1442 }
1443 IntegerValue delta_energy_up = 0;
1444 IntegerValue delta_energy_down = 0;
1445
1446 // Note that the non-deterministic iteration order is fine here.
1447 for (const int idx : ranges_touching_both_boundaries_[dimension]) {
1448 const RectangleInRange& range = intervals_[idx];
1449 const IntegerValue curr_x = Smallest1DIntersection(
1450 range.bounding_area.x_min, range.bounding_area.x_max, range.x_size,
1451 current_rectangle.x_min, current_rectangle.x_max);
1452 const IntegerValue curr_y = Smallest1DIntersection(
1453 range.bounding_area.y_min, range.bounding_area.y_max, range.y_size,
1454 current_rectangle.y_min, current_rectangle.y_max);
1455 const IntegerValue curr = curr_x * curr_y;
1456 delta_energy_up += curr;
1457 delta_energy_down += curr;
1458
1459 if (dimension == 0) {
1460 const IntegerValue up_x = Smallest1DIntersection(
1461 range.bounding_area.x_min, range.bounding_area.x_max, range.x_size,
1462 next_rectangle_up.x_min, next_rectangle_up.x_max);
1463 const IntegerValue down_x = Smallest1DIntersection(
1464 range.bounding_area.x_min, range.bounding_area.x_max, range.x_size,
1465 next_rectangle_down.x_min, next_rectangle_down.x_max);
1466
1467 delta_energy_up -= curr_y * up_x;
1468 delta_energy_down -= curr_y * down_x;
1469 } else {
1470 const IntegerValue up_y = Smallest1DIntersection(
1471 range.bounding_area.y_min, range.bounding_area.y_max, range.y_size,
1472 next_rectangle_up.y_min, next_rectangle_up.y_max);
1473 const IntegerValue down_y = Smallest1DIntersection(
1474 range.bounding_area.y_min, range.bounding_area.y_max, range.y_size,
1475 next_rectangle_down.y_min, next_rectangle_down.y_max);
1476
1477 delta_energy_up -= curr_x * up_y;
1478 delta_energy_down -= curr_x * down_y;
1479 }
1480 }
1481 delta_energy_up += units_crossed_up * step_1d_size_up;
1482 delta_energy_down += units_crossed_down * step_1d_size_down;
1483 *delta_energy_up_ptr = delta_energy_up;
1484 *delta_energy_down_ptr = delta_energy_down;
1485}
1486
1487bool ProbingRectangle::CanShrink(Edge edge) const {
1488 switch (edge) {
1490 case Edge::RIGHT:
1491 return (next_indexes_[Edge::RIGHT] > indexes_[Edge::LEFT]);
1492 case Edge::BOTTOM:
1493 case Edge::TOP:
1494 return (indexes_[Edge::TOP] > next_indexes_[Edge::BOTTOM]);
1495 }
1496}
1497
1498namespace {
1499std::vector<double> GetExpTable() {
1500 std::vector<double> table(101);
1501 for (int i = 0; i <= 100; ++i) {
1502 table[i] = std::exp(-(i - 50) / 5.0);
1503 }
1504 return table;
1505}
1506
1507} // namespace
1508
1509FindRectanglesResult FindRectanglesWithEnergyConflictMC(
1510 const std::vector<RectangleInRange>& intervals, absl::BitGenRef random,
1511 double temperature, double candidate_energy_usage_factor) {
1512 FindRectanglesResult result;
1513 ProbingRectangle ranges(intervals);
1514
1515 static const std::vector<double>* cached_probabilities =
1516 new std::vector<double>(GetExpTable());
1517
1518 const double inv_temp = 1.0 / temperature;
1519 absl::InlinedVector<ProbingRectangle::Edge, 4> candidates;
1520 absl::InlinedVector<IntegerValue, 4> energy_deltas;
1521 absl::InlinedVector<double, 4> weights;
1522 while (!ranges.IsMinimal()) {
1523 const IntegerValue rect_area = ranges.GetCurrentRectangleArea();
1524 const IntegerValue min_energy = ranges.GetMinimumEnergy();
1525 if (min_energy > rect_area) {
1526 result.conflicts.push_back(ranges.GetCurrentRectangle());
1527 } else if (min_energy.value() >
1528 candidate_energy_usage_factor * rect_area.value()) {
1529 result.candidates.push_back(ranges.GetCurrentRectangle());
1530 }
1531 if (min_energy == 0) {
1532 break;
1533 }
1534 candidates.clear();
1535 energy_deltas.clear();
1536
1537 for (int border_idx = 0; border_idx < 4; ++border_idx) {
1538 const ProbingRectangle::Edge border =
1539 static_cast<ProbingRectangle::Edge>(border_idx);
1540 if (!ranges.CanShrink(border)) {
1541 continue;
1542 }
1543 candidates.push_back(border);
1544 const IntegerValue delta_area = ranges.GetShrinkDeltaArea(border);
1545 const IntegerValue delta_energy = ranges.GetShrinkDeltaEnergy(border);
1546 energy_deltas.push_back(delta_energy - delta_area);
1547 }
1548 const IntegerValue min_energy_delta =
1549 *std::min_element(energy_deltas.begin(), energy_deltas.end());
1550 weights.clear();
1551 for (const IntegerValue delta_slack : energy_deltas) {
1552 const int64_t table_lookup =
1553 std::max((int64_t)0,
1554 std::min((int64_t)((delta_slack - min_energy_delta).value() *
1555 5 * inv_temp +
1556 50),
1557 (int64_t)100));
1558 weights.push_back((*cached_probabilities)[table_lookup]);
1559 }
1560 // Pick a change with a probability proportional to exp(- delta_E / Temp)
1561 ranges.Shrink(candidates[WeightedPick(weights, random)]);
1562 }
1563 if (ranges.GetMinimumEnergy() > ranges.GetCurrentRectangleArea()) {
1564 result.conflicts.push_back(ranges.GetCurrentRectangle());
1565 }
1566 return result;
1567}
1568
1569std::string RenderDot(std::optional<Rectangle> bb,
1570 absl::Span<const Rectangle> solution,
1571 std::string_view extra_dot_payload) {
1572 const std::vector<std::string> colors = {"#0000ff80", "#ee00ee80",
1573 "#ff000080", "#eeee0080",
1574 "#00ff0080", "#00eeee80"};
1575 std::stringstream ss;
1576 ss << "digraph {\n";
1577 ss << " graph [ bgcolor=lightgray ]\n";
1578 ss << " node [style=filled shape=box]\n";
1579 if (bb.has_value()) {
1580 ss << " bb [fillcolor=\"grey\" pos=\"" << 2 * bb->x_min + bb->SizeX()
1581 << "," << 2 * bb->y_min + bb->SizeY() << "!\" width=" << 2 * bb->SizeX()
1582 << " height=" << 2 * bb->SizeY() << "]\n";
1583 }
1584 for (int i = 0; i < solution.size(); ++i) {
1585 ss << " " << i << " [fillcolor=\"" << colors[i % colors.size()]
1586 << "\" pos=\"" << 2 * solution[i].x_min + solution[i].SizeX() << ","
1587 << 2 * solution[i].y_min + solution[i].SizeY()
1588 << "!\" width=" << 2 * solution[i].SizeX()
1589 << " height=" << 2 * solution[i].SizeY() << "]\n";
1590 }
1591 ss << extra_dot_payload;
1592 ss << "}\n";
1593 return ss.str();
1594}
1595
1596std::vector<Rectangle> FindEmptySpaces(
1597 const Rectangle& bounding_box, std::vector<Rectangle> ocupied_rectangles) {
1598 // Sorting is not necessary for correctness but makes it faster.
1599 std::sort(ocupied_rectangles.begin(), ocupied_rectangles.end(),
1600 [](const Rectangle& a, const Rectangle& b) {
1601 return std::tuple(a.x_min, -a.x_max, a.y_min) <
1602 std::tuple(b.x_min, -b.x_max, b.y_min);
1603 });
1604 return PavedRegionDifference({bounding_box}, ocupied_rectangles);
1605}
1606
1607std::vector<Rectangle> PavedRegionDifference(
1608 std::vector<Rectangle> original_region,
1609 absl::Span<const Rectangle> area_to_remove) {
1610 std::vector<Rectangle> new_area_to_cover;
1611 new_area_to_cover.reserve(original_region.size());
1612 for (const Rectangle& rectangle : area_to_remove) {
1613 new_area_to_cover.clear();
1614 for (const Rectangle& r : original_region) {
1615 const auto& new_rectangles = r.RegionDifference(rectangle);
1616 new_area_to_cover.insert(new_area_to_cover.end(), new_rectangles.begin(),
1617 new_rectangles.end());
1618 }
1619 original_region.swap(new_area_to_cover);
1620 if (original_region.empty()) break;
1621 }
1622 return original_region;
1623}
1624
1625// Each node in the tree will hold either a single box that is covering the
1626// whole interval represented by the node, or, if no such box exists, a superset
1627// of all the connected components of boxes that are overlapping the interval.
1628// It is a superset and not the exact set of connected components because we
1629// don't delete nodes that became stale, as explained in the class comment
1630// below.
1631struct BinaryTreeNode {
1632 // Contains exactly one element if occupying_box_index != -1.
1633 absl::flat_hash_set<int> connected_components_descendants;
1634 // Hold the x_max of the box that is currently occupying this node (if any) to
1635 // know when it is stale.
1637 // -1 if not occupied.
1639};
1641// A data structure to store which boxes are overlapping the current sweep line.
1642// This uses a binary tree in a slight non-standard way: in a typical use of a
1643// binary tree the actual values are stored in the leaves and the intermediate
1644// nodes are there just to make finding the right leaf efficient. Here we do the
1645// opposite: the values are stored as high up in the tree as possible.
1646// For example, for a tree of size 8 a box that occupies the y interval [0, 7]
1647// will be stored as a single node at the root. In the same tree, a box that
1648// occupies [3, 7] will be stored with the nodes representing the [3, 4), [4, 6)
1649// and [6, 8) intervals. There is no difference on what is stored in the
1650// intermediate nodes or on the leaves. When the sweep line moves, we don't
1651// update the existing nodes on the tree. Thus, some nodes will become stale and
1652// will represent boxes that no longer overlap the sweep line. Those stale nodes
1653// get removed lazily.
1654struct SweepLineIntervalTree {
1655 explicit SweepLineIntervalTree(int max_y, int num_boxes)
1656 : tree(LeafIndex(max_y + 1)), tree_nodes(tree.StorageSize()) {
1657 union_find.SetNumberOfNodes(num_boxes);
1658 }
1659
1660 // Recompute the connected components of a given node, by simply setting it to
1661 // {self} + left.connected_components + right.connected_components.
1662 void RecomputeConnectedComponents(TreeNodeIndex idx) {
1663 BinaryTreeNode& node = tree_nodes[idx];
1664 if (node.occupying_box_index != -1) {
1665 node.connected_components_descendants = {
1666 union_find.FindRoot(node.occupying_box_index)};
1667 return;
1668 }
1669 node.connected_components_descendants.clear();
1670 if (tree.IsLeaf(idx)) return;
1671 for (const TreeNodeIndex child_idx :
1672 {tree.LeftChild(idx), tree.RightChild(idx)}) {
1673 // The order is non-deterministic, but since this is doing the union of
1674 // hash sets the result is deterministic.
1675 for (const int c :
1676 tree_nodes[child_idx].connected_components_descendants) {
1677 node.connected_components_descendants.insert(union_find.FindRoot(c));
1678 }
1679 }
1680 }
1681
1682 // We don't have global deletion method in this class, but this method
1683 // checks if a single interval is fully to the left of the sweep line and
1684 // removes it if so, also updating its connected components.
1685 void RemoveNodeIfXMaxLowerOrEqual(TreeNodeIndex idx, int x_threshold) {
1686 BinaryTreeNode& node = tree_nodes[idx];
1687 if (node.occupying_box_index == -1) {
1688 // Node is already empty.
1689 return;
1690 }
1691 if (node.occupying_box_x_max > x_threshold) {
1692 // Node is still overlapping the sweep line.
1693 return;
1694 }
1695 node.occupying_box_index = -1;
1696 RecomputeConnectedComponents(idx);
1697 }
1698
1699 void UpdateChildrenIntersecting(TreeNodeIndex idx, int sweep_line_x_pos,
1700 int component_index,
1701 std::vector<int>* new_connections) {
1702 if (tree.IsLeaf(idx)) return;
1703 for (const TreeNodeIndex child_idx :
1704 {tree.LeftChild(idx), tree.RightChild(idx)}) {
1705 RemoveNodeIfXMaxLowerOrEqual(child_idx, sweep_line_x_pos);
1706 BinaryTreeNode& child_node = tree_nodes[child_idx];
1707 if (child_node.occupying_box_index != -1) {
1708 if (union_find.AddEdge(child_node.occupying_box_index,
1709 component_index)) {
1710 new_connections->push_back(child_node.occupying_box_index);
1711 }
1712 // No need to recurse here: we already connected everything on this
1713 // branch to the new box.
1714 continue;
1715 }
1716 const bool had_different_component =
1717 absl::c_any_of(child_node.connected_components_descendants,
1718 [this, component_index](const int c) {
1719 return !union_find.Connected(c, component_index);
1720 });
1721 // Since everything is intersecting the current box, all descendants
1722 // must be in one single component.
1723 child_node.connected_components_descendants = {component_index};
1724
1725 // Only go down on the tree if we have below either:
1726 // - a different component to connect.
1727 // - a box to remove that is in a different component.
1728 // In any case, we will visit O(c + d) terminals, where c is the number of
1729 // components we are connecting and d is the number of boxes that we will
1730 // delete. Since a box can only be deleted log N times (one per interval
1731 // it was cut into) and we can only connect O(N) components in total, the
1732 // amortized cost of a call to UpdateChildrenIntersecting is O((log N)^2).
1733 if (had_different_component) {
1734 UpdateChildrenIntersecting(child_idx, sweep_line_x_pos, component_index,
1735 new_connections);
1736 }
1737 }
1738 }
1739
1740 bool UpdateParents(TreeNodeIndex node, int sweep_line_x_pos,
1741 int component_index, std::vector<int>* new_connections) {
1742 if (node == tree.Root()) return false;
1743 for (TreeNodeIndex parent = tree.Parent(node); parent != tree.Root();
1744 parent = tree.Parent(parent)) {
1745 RemoveNodeIfXMaxLowerOrEqual(parent, sweep_line_x_pos);
1746 BinaryTreeNode& parent_value = tree_nodes[parent];
1747 if (parent_value.occupying_box_index != -1) {
1748 if (union_find.AddEdge(parent_value.occupying_box_index,
1749 component_index)) {
1750 new_connections->push_back(parent_value.occupying_box_index);
1751 return true;
1752 }
1753 }
1754 parent_value.connected_components_descendants.insert(component_index);
1755 }
1756 return false;
1757 }
1758
1759 // Add a new box to the sweep line. This will store it in the tree (split in
1760 // log N intervals) check if it connects to one or more existing connected
1761 // components, and for each case it does, add the box that it is overlapping
1762 // to new_connections.
1763 void AddInterval(TreeNodeIndex idx, int sweep_line_x_pos, int box_index,
1764 int x_max, std::vector<int>* new_connections) {
1765 RemoveNodeIfXMaxLowerOrEqual(idx, sweep_line_x_pos);
1766 int cur_box_component = union_find.FindRoot(box_index);
1767 BinaryTreeNode& node = tree_nodes[idx];
1768 if (node.occupying_box_index == -1) {
1769 node.connected_components_descendants = {box_index};
1770 node.occupying_box_index = box_index;
1771 node.occupying_box_x_max = x_max;
1772 const bool had_occupied_parent = UpdateParents(
1773 idx, sweep_line_x_pos, cur_box_component, new_connections);
1774 // We can only be connecting children if it is not already connect via
1775 // something above on the tree.
1776 if (!had_occupied_parent) {
1777 UpdateChildrenIntersecting(idx, sweep_line_x_pos, cur_box_component,
1778 new_connections);
1779 }
1780 } else {
1781 // We have already something fully occupying this interval.
1782 if (union_find.AddEdge(node.occupying_box_index, cur_box_component)) {
1783 new_connections->push_back(node.occupying_box_index);
1784 cur_box_component = union_find.FindRoot(cur_box_component);
1785 }
1786 node.connected_components_descendants = {cur_box_component};
1787 if (node.occupying_box_x_max < x_max) {
1788 // Replace the existing box by the new one.
1789 node.occupying_box_index = box_index;
1790 node.occupying_box_x_max = x_max;
1791 }
1792 }
1793 }
1794
1795 FixedShapeBinaryTree tree;
1796 util_intops::StrongVector<TreeNodeIndex, BinaryTreeNode> tree_nodes;
1800struct Rectangle32 {
1801 int x_min;
1808// Requires that rectangles are sorted by x_min and that sizes on both
1809// dimensions are > 0.
1810std::vector<std::pair<int, int>> FindPartialRectangleIntersectionsImpl(
1811 absl::Span<const Rectangle32> rectangles, int32_t y_max) {
1812 // We are going to use a sweep line algorithm to find the intersections.
1813 // First, we sort the rectangles by their x coordinates, then consider a sweep
1814 // line that goes from the left to the right. See the comment on the
1815 // SweepLineIntervalTree class for more details about what we store for each
1816 // line.
1817 DCHECK(std::is_sorted(rectangles.begin(), rectangles.end(),
1818 [](const Rectangle32& a, const Rectangle32& b) {
1819 return a.x_min < b.x_min;
1820 }));
1821
1822 SweepLineIntervalTree interval_tree(y_max, rectangles.size());
1823
1824 std::vector<TreeNodeIndex> interval_pieces;
1825 std::vector<int> new_connections;
1826 std::vector<std::pair<int, int>> arcs;
1827 for (int rectangle_index = 0; rectangle_index < rectangles.size();
1828 ++rectangle_index) {
1829 DCHECK_LT(rectangles[rectangle_index].x_min,
1830 rectangles[rectangle_index].x_max);
1831 DCHECK_LT(rectangles[rectangle_index].y_min,
1832 rectangles[rectangle_index].y_max);
1833 const int sweep_line_x_pos = rectangles[rectangle_index].x_min;
1834 const Rectangle32& r = rectangles[rectangle_index];
1835 interval_pieces.clear();
1836 interval_tree.tree.PartitionIntervalIntoNodes(
1837 LeafIndex(r.y_min), LeafIndex(r.y_max - 1), &interval_pieces);
1838 new_connections.clear();
1839 for (const TreeNodeIndex& node : interval_pieces) {
1840 interval_tree.AddInterval(node, sweep_line_x_pos, rectangle_index,
1841 r.x_max, &new_connections);
1842 }
1843 for (const int new_connection : new_connections) {
1844 arcs.push_back({rectangles[new_connection].index, r.index});
1845 }
1846 }
1847 return arcs;
1848}
1849
1850namespace {
1851struct PostProcessedResult {
1852 std::vector<Rectangle32> rectangles_sorted_by_x_min;
1853 std::pair<int32_t, int32_t> bounding_box; // Always starting at (0,0).
1854};
1855
1856PostProcessedResult ConvertToRectangle32WithNonZeroSizes(
1857 absl::Span<const Rectangle> rectangles) {
1858 // This function is a preprocessing function for algorithms that find overlap
1859 // between rectangles. It does the following:
1860 // - It converts the arbitrary int64_t coordinates into a small integer by
1861 // sorting the possible values and assigning them consecutive integers.
1862 // - It grows zero size intervals to make them size one. This simplifies
1863 // things considerably, since it is hard to reason about degenerated
1864 // rectangles in the general algorithm.
1865 //
1866 // Note that the last point need to be done with care. Imagine the following
1867 // example:
1868 // +----------+
1869 // | |
1870 // | +--------------+
1871 // | | |
1872 // | | p,q r |
1873 // | +----*-----*-+-+
1874 // | | |
1875 // | | |
1876 // | | |
1877 // | +------------+
1878 // | |
1879 // | |
1880 // +----------+
1881 // Where p,q and r are points (ie, boxes of size 0x0) and p and q have the
1882 // same coordinates. We replace them by the following:
1883 // +----------+
1884 // | |
1885 // | +----------------------+
1886 // | | |
1887 // | | |
1888 // | +----+-+---------------+
1889 // | | |p|
1890 // | | +-+-+
1891 // | | |q|
1892 // | | +-+ +-+
1893 // | | |r|
1894 // | +--------------+-+---+
1895 // | | |
1896 // | | |
1897 // | | |
1898 // | +--------------------+
1899 // | |
1900 // | |
1901 // +----------+
1902 //
1903 // That is a pretty radical deformation of the original shape, but it retains
1904 // the property of whether a pair of rectangles intersect or not.
1905
1906 if (rectangles.empty()) return {};
1907
1908 enum class Event {
1909 kEnd = 0,
1910 kPoint = 1,
1911 kBegin = 2,
1912 };
1913 std::vector<std::tuple<IntegerValue, Event, int>> x_events;
1914 std::vector<std::tuple<IntegerValue, Event, int>> y_events;
1915 x_events.reserve(rectangles.size() * 2);
1916 y_events.reserve(rectangles.size() * 2);
1917 for (int i = 0; i < rectangles.size(); ++i) {
1918 const Rectangle& r = rectangles[i];
1919 DCHECK_GE(r.SizeX(), 0);
1920 DCHECK_GE(r.SizeY(), 0);
1921 if (r.SizeX() == 0) {
1922 x_events.push_back({r.x_min, Event::kPoint, i});
1923 } else {
1924 x_events.push_back({r.x_min, Event::kBegin, i});
1925 x_events.push_back({r.x_max, Event::kEnd, i});
1926 }
1927 if (r.SizeY() == 0) {
1928 y_events.push_back({r.y_min, Event::kPoint, i});
1929 } else {
1930 y_events.push_back({r.y_min, Event::kBegin, i});
1931 y_events.push_back({r.y_max, Event::kEnd, i});
1932 }
1933 }
1934 std::sort(y_events.begin(), y_events.end());
1935
1936 std::vector<Rectangle32> rectangles32;
1937 rectangles32.resize(rectangles.size());
1938 IntegerValue prev_y = 0;
1939 Event prev_event = Event::kEnd;
1940 int cur_index = -1;
1941 for (const auto [y, event, index] : y_events) {
1942 if ((prev_event != event && prev_event != Event::kEnd) || prev_y != y ||
1943 event == Event::kPoint || cur_index == -1) {
1944 ++cur_index;
1945 }
1946
1947 switch (event) {
1948 case Event::kBegin:
1949 rectangles32[index].y_min = cur_index;
1950 rectangles32[index].index = index;
1951 break;
1952 case Event::kEnd:
1953 rectangles32[index].y_max = cur_index;
1954 break;
1955 case Event::kPoint:
1956 rectangles32[index].y_min = cur_index;
1957 rectangles32[index].y_max = cur_index + 1;
1958 rectangles32[index].index = index;
1959 break;
1960 }
1961 prev_event = event;
1962 prev_y = y;
1963 }
1964 const int max_y_index = cur_index + 1;
1965
1966 gtl::STLClearObject(&y_events);
1967
1968 std::sort(x_events.begin(), x_events.end());
1969 IntegerValue prev_x = 0;
1970 prev_event = Event::kEnd;
1971 cur_index = -1;
1972 for (const auto [x, event, index] : x_events) {
1973 if ((prev_event != event && prev_event != Event::kEnd) || prev_x != x ||
1974 event == Event::kPoint || cur_index == -1) {
1975 ++cur_index;
1976 }
1977
1978 switch (event) {
1979 case Event::kBegin:
1980 rectangles32[index].x_min = cur_index;
1981 break;
1982 case Event::kEnd:
1983 rectangles32[index].x_max = cur_index;
1984 break;
1985 case Event::kPoint:
1986 rectangles32[index].x_min = cur_index;
1987 rectangles32[index].x_max = cur_index + 1;
1988 break;
1989 }
1990 prev_event = event;
1991 prev_x = x;
1992 }
1993 const int max_x_index = cur_index + 1;
1994
1995 std::vector<Rectangle32> sorted_rectangles32;
1996 sorted_rectangles32.reserve(rectangles.size());
1997 for (const auto [x, event, index] : x_events) {
1998 if (event == Event::kBegin || event == Event::kPoint) {
1999 sorted_rectangles32.push_back(rectangles32[index]);
2000 }
2001 }
2002
2003 return {sorted_rectangles32, {max_x_index, max_y_index}};
2004}
2005template <typename RectangleT>
2006std::optional<std::pair<int, int>> FindOneIntersectionIfPresentImpl(
2007 absl::Span<const RectangleT> rectangles) {
2008 using CoordinateType = std::decay_t<decltype(RectangleT::x_min)>;
2009 DCHECK(absl::c_is_sorted(rectangles,
2010 [](const RectangleT& a, const RectangleT& b) {
2011 return a.x_min < b.x_min;
2012 }));
2013
2014 // Set of box intersection the sweep line. We only store y_min, other
2015 // coordinates can be accessed via rectangles[index].coordinate.
2016 struct Element {
2017 mutable int index;
2018 CoordinateType y_min;
2019 bool operator<(const Element& other) const { return y_min < other.y_min; }
2020 };
2021
2022 // Note: To use btree_set that has no iterator stability, we have to be
2023 // a bit careful below.
2024 absl::btree_set<Element> interval_set;
2025
2026 for (int i = 0; i < rectangles.size(); ++i) {
2027 const CoordinateType x = rectangles[i].x_min;
2028 const CoordinateType y_min = rectangles[i].y_min;
2029 const CoordinateType y_max = rectangles[i].y_max;
2030
2031 // TODO(user): We can handle that, but it require some changes below.
2032 DCHECK_LE(y_min, y_max);
2033
2034 // Try to add this rectangle to the set, if there is an intersection, lazily
2035 // remove it if its x_max is already passed, otherwise report the
2036 // intersection.
2037 auto [it, inserted] = interval_set.insert({i, y_min});
2038 if (!inserted) {
2039 if (rectangles[it->index].x_max <= x) {
2040 // We just replace if the rectangle at position i is stale.
2041 it->index = i;
2042 } else {
2043 // Intersection.
2044 return {{it->index, i}};
2045 }
2046 } else {
2047 // If there was no element at position y_min, we need to test if the
2048 // interval before is stale or if it overlap with the new one.
2049 if (it != interval_set.begin()) {
2050 auto it_before = it;
2051 --it_before;
2052
2053 // Lazy erase stale entry.
2054 if (rectangles[it_before->index].x_max <= x) {
2055 // For absl::btree_set we don't have iterator stability, so we do need
2056 // to re-assign 'it' to the element just after the one we erased.
2057 it = interval_set.erase(it_before);
2058 } else {
2059 DCHECK_LE(it_before->y_min, y_min);
2060 const CoordinateType y_max_before =
2061 rectangles[it_before->index].y_max;
2062 if (y_max_before > y_min) {
2063 // Intersection.
2064 return {{it_before->index, i}};
2065 }
2066 }
2067 }
2068 }
2069
2070 // We handled the part before, now we need to deal with the interval that
2071 // starts after y_min.
2072 ++it;
2073 while (it != interval_set.end()) {
2074 // Lazy erase stale entry.
2075 if (rectangles[it->index].x_max <= x) {
2076 it = interval_set.erase(it);
2077 continue;
2078 }
2079
2080 DCHECK_LE(y_min, it->y_min);
2081 if (y_max > it->y_min) {
2082 // Intersection.
2083 return {{it->index, i}};
2084 }
2085 break;
2086 }
2087 }
2088
2089 return {};
2090}
2091
2092} // namespace
2093
2094std::vector<std::pair<int, int>> FindPartialRectangleIntersections(
2095 absl::Span<const Rectangle> rectangles) {
2096 auto postprocessed = ConvertToRectangle32WithNonZeroSizes(rectangles);
2098 postprocessed.rectangles_sorted_by_x_min,
2099 postprocessed.bounding_box.second);
2100}
2101
2102std::optional<std::pair<int, int>> FindOneIntersectionIfPresent(
2103 absl::Span<const Rectangle> rectangles) {
2104 return FindOneIntersectionIfPresentImpl(rectangles);
2105}
2106
2107std::optional<std::pair<int, int>> FindOneIntersectionIfPresentWithZeroArea(
2108 absl::Span<const Rectangle> rectangles) {
2109 auto postprocessed = ConvertToRectangle32WithNonZeroSizes(rectangles);
2110 std::optional<std::pair<int, int>> result = FindOneIntersectionIfPresentImpl(
2111 absl::MakeConstSpan(postprocessed.rectangles_sorted_by_x_min));
2112 if (!result.has_value()) return {};
2113 return {{postprocessed.rectangles_sorted_by_x_min[result->first].index,
2114 postprocessed.rectangles_sorted_by_x_min[result->second].index}};
2115}
2116
2117} // namespace sat
2118} // namespace operations_research
A connected components finder that only works on dense ints.
void AddMandatoryConsumption(IntegerValue x_min, IntegerValue x_max, IntegerValue y_height)
void AddRectangle(IntegerValue x_min, IntegerValue x_max, IntegerValue y_min, IntegerValue y_max)
Adds a rectangle to the current shape.
void BuildResidualCapacityProfile(std::vector< Rectangle > *result)
void ResetFromPairs(const Collection &pairs, int minimum_num_nodes=0)
Similar to ResetFromFlatMapping().
Definition util.h:895
IntegerValue GetShrinkDeltaEnergy(Edge edge) const
How much of GetMinimumEnergy() will change if Shrink() is called.
Definition diffn_util.h:599
void Reset()
Reset to the bounding box of the whole set.
IntegerValue GetShrinkDeltaArea(Edge edge) const
How much of GetCurrentRectangleArea() will change if Shrink() is called.
void ClearReason()
Functions to clear and then set the current reason.
void AddEnergyMinInIntervalReason(int t, IntegerValue min, IntegerValue max)
void STLSortAndRemoveDuplicates(T *v, const LessFunc &less_func)
Definition stl_util.h:55
void STLClearObject(T *obj)
Definition stl_util.h:120
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor)
CompactVectorVector< int > GetOverlappingRectangleComponents(absl::Span< const Rectangle > rectangles)
absl::Span< int > FilterBoxesAndRandomize(absl::Span< const Rectangle > cached_rectangles, absl::Span< int > boxes, IntegerValue threshold_x, IntegerValue threshold_y, absl::BitGenRef random)
bool AnalyzeIntervals(bool transpose, absl::Span< const int > local_boxes, absl::Span< const Rectangle > rectangles, absl::Span< const IntegerValue > rectangle_energies, IntegerValue *x_threshold, IntegerValue *y_threshold, Rectangle *conflict)
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue.value())
absl::Span< int > FilterBoxesThatAreTooLarge(absl::Span< const Rectangle > cached_rectangles, absl::Span< const IntegerValue > energies, absl::Span< int > boxes)
std::vector< int > GetIntervalArticulationPoints(std::vector< IndexedInterval > *intervals)
void AppendPairwiseRestrictions(absl::Span< const ItemWithVariableSize > items, std::vector< PairwiseRestriction > *result)
void ConstructOverlappingSets(absl::Span< IndexedInterval > intervals, CompactVectorVector< int > *result, absl::Span< const int > order)
int WeightedPick(absl::Span< const double > input, absl::BitGenRef random)
Definition util.cc:426
bool ReportEnergyConflict(Rectangle bounding_box, absl::Span< const int > boxes, SchedulingConstraintHelper *x, SchedulingConstraintHelper *y)
std::vector< std::pair< int, int > > FindPartialRectangleIntersectionsImpl(absl::Span< const Rectangle32 > rectangles, int32_t y_max)
FindRectanglesResult FindRectanglesWithEnergyConflictMC(const std::vector< RectangleInRange > &intervals, absl::BitGenRef random, double temperature, double candidate_energy_usage_factor)
std::vector< Rectangle > PavedRegionDifference(std::vector< Rectangle > original_region, absl::Span< const Rectangle > area_to_remove)
bool BoxesAreInEnergyConflict(absl::Span< const Rectangle > rectangles, absl::Span< const IntegerValue > energies, absl::Span< const int > boxes, Rectangle *conflict)
std::vector< std::pair< int, int > > FindPartialRectangleIntersections(absl::Span< const Rectangle > rectangles)
IntegerValue Smallest1DIntersection(IntegerValue range_min, IntegerValue range_max, IntegerValue size, IntegerValue interval_min, IntegerValue interval_max)
void GetOverlappingIntervalComponents(std::vector< IndexedInterval > *intervals, std::vector< std::vector< int > > *components)
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
ClosedInterval::Iterator end(ClosedInterval interval)
const bool DEBUG_MODE
Definition radix_sort.h:266
void FindStronglyConnectedComponents(NodeIndex num_nodes, const Graph &graph, SccOutput *components)
Simple wrapper function for most usage.
int occupying_box_index
-1 if not occupied.
absl::flat_hash_set< int > connected_components_descendants
Contains exactly one element if occupying_box_index != -1.
void GrowToInclude(const Rectangle &other)
Definition diffn_util.h:50
absl::InlinedVector< Rectangle, 4 > RegionDifference(const Rectangle &other) const
Definition diffn_util.cc:64
bool IsDisjoint(const Rectangle &other) const
Definition diffn_util.cc:59
Rectangle Intersect(const Rectangle &other) const
The methods below are not meant to be used with zero-area rectangles.
Definition diffn_util.h:105
DenseConnectedComponentsFinder union_find
void RemoveNodeIfXMaxLowerOrEqual(TreeNodeIndex idx, int x_threshold)
void UpdateChildrenIntersecting(TreeNodeIndex idx, int sweep_line_x_pos, int component_index, std::vector< int > *new_connections)
util_intops::StrongVector< TreeNodeIndex, BinaryTreeNode > tree_nodes