Google OR-Tools v9.12
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
diffn_cuts.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 <cmath>
18#include <string>
19#include <tuple>
20#include <utility>
21#include <vector>
22
23#include "absl/base/attributes.h"
24#include "absl/log/check.h"
25#include "absl/strings/str_cat.h"
26#include "absl/strings/string_view.h"
27#include "absl/types/span.h"
30#include "ortools/sat/cuts.h"
33#include "ortools/sat/integer.h"
37#include "ortools/sat/model.h"
41#include "ortools/sat/util.h"
43
44namespace operations_research {
45namespace sat {
46
47namespace {
48
49// Minimum amount of violation of the cut constraint by the solution. This
50// is needed to avoid numerical issues and adding cuts with minor effect.
51const double kMinCutViolation = 1e-4;
52
53} // namespace
54
56 const SchedulingConstraintHelper* x_helper)
57 : x_start_min(x_helper->StartMin(t)),
58 x_start_max(x_helper->StartMax(t)),
59 x_end_min(x_helper->EndMin(t)),
60 x_end_max(x_helper->EndMax(t)),
61 x_size_min(x_helper->SizeMin(t)) {}
62
65 : DiffnBaseEvent(t, x_helper) {}
66
67 // We need this for linearizing the energy in some cases.
69
70 // If set, this event is optional and its presence is controlled by this.
72
73 // A linear expression which is a valid lower bound on the total energy of
74 // this event. We also cache the activity of the expression to not recompute
75 // it all the time.
78
79 // True if linearized_energy is not exact and a McCormick relaxation.
80 bool energy_is_quadratic = false;
81
82 // Used to minimize the increase on the y axis for rectangles.
83 double y_spread = 0.0;
84
85 // The actual value of the presence literal of the interval(s) is checked
86 // when the event is created. A value of kNoLiteralIndex indicates that either
87 // the interval was not optional, or that its presence literal is true at
88 // level zero.
90
91 // Computes the mandatory minimal overlap of the interval with the time window
92 // [start, end].
93 IntegerValue GetMinOverlap(IntegerValue start, IntegerValue end) const {
94 return std::max(std::min({x_end_min - start, end - x_start_max, x_size_min,
95 end - start}),
96 IntegerValue(0));
97 }
98
99 // This method expects all the other fields to have been filled before.
100 // It must be called before the EnergyEvent is used.
101 ABSL_MUST_USE_RESULT bool FillEnergyLp(
102 AffineExpression x_size,
104 Model* model) {
105 LinearConstraintBuilder tmp_energy(model);
106 if (IsPresent()) {
107 if (!decomposed_energy.empty()) {
108 if (!tmp_energy.AddDecomposedProduct(decomposed_energy)) return false;
109 } else {
110 tmp_energy.AddQuadraticLowerBound(x_size, y_size,
111 model->GetOrCreate<IntegerTrail>(),
113 }
114 } else {
116 energy_min)) {
117 return false;
118 }
119 }
120 linearized_energy = tmp_energy.BuildExpression();
122 return true;
123 }
124
125 std::string DebugString() const {
126 return absl::StrCat(
127 "DiffnEnergyEvent(x_start_min = ", x_start_min.value(),
128 ", x_start_max = ", x_start_max.value(),
129 ", x_end_min = ", x_end_min.value(),
130 ", x_end_max = ", x_end_max.value(), ", y_min = ", y_min.value(),
131 ", y_max = ", y_max.value(), ", y_size = ", y_size.DebugString(),
132 ", energy = ",
133 decomposed_energy.empty()
134 ? "{}"
135 : absl::StrCat(decomposed_energy.size(), " terms"),
136 ", presence_literal_index = ", presence_literal_index.value(), ")");
137 }
138};
139
141 absl::Span<const std::vector<LiteralValueValue>> energies,
142 absl::Span<const int> rectangles, absl::string_view cut_name, Model* model,
145 SchedulingDemandHelper* y_demands_helper) {
146 std::vector<DiffnEnergyEvent> events;
147 const auto& lp_values = manager->LpValues();
148 for (const int rect : rectangles) {
149 if (y_helper->SizeMax(rect) == 0 || x_helper->SizeMax(rect) == 0) {
150 continue;
151 }
152
153 DiffnEnergyEvent e(rect, x_helper);
154 e.y_min = y_helper->StartMin(rect);
155 e.y_max = y_helper->EndMax(rect);
156 e.y_size = y_helper->Sizes()[rect];
157 e.decomposed_energy = energies[rect];
159 x_helper->IsPresent(rect)
160 ? (y_helper->IsPresent(rect)
162 : y_helper->PresenceLiteral(rect).Index())
163 : x_helper->PresenceLiteral(rect).Index();
164 e.y_size_min = y_helper->SizeMin(rect);
165 e.energy_min = y_demands_helper->EnergyMin(rect);
166 e.energy_is_quadratic = y_demands_helper->EnergyIsQuadratic(rect);
167
168 // We can always skip events.
169 if (!e.FillEnergyLp(x_helper->Sizes()[rect], lp_values, model)) continue;
170 events.push_back(e);
171 }
172
173 if (events.empty()) return;
174
175 // Compute y_spread.
176 double average_d = 0.0;
177 for (const auto& e : events) {
178 average_d += ToDouble(e.y_min + e.y_max);
179 }
180 const double average = average_d / 2.0 / static_cast<double>(events.size());
181 for (auto& e : events) {
182 e.y_spread = std::abs(ToDouble(e.y_max) - average) +
183 std::abs(average - ToDouble(e.y_min));
184 }
185
186 TopNCuts top_n_cuts(5);
187
188 std::sort(events.begin(), events.end(),
189 [](const DiffnEnergyEvent& a, const DiffnEnergyEvent& b) {
190 return std::tie(a.x_start_min, a.y_spread, a.x_end_max) <
191 std::tie(b.x_start_min, b.y_spread, b.x_end_max);
192 });
193
194 // The sum of all energies can be used to stop iterating early.
195 double sum_of_all_energies = 0.0;
196 for (const auto& e : events) {
197 sum_of_all_energies += e.linearized_energy_lp_value;
198 }
199
200 CapacityProfile capacity_profile;
201 for (int i1 = 0; i1 + 1 < events.size(); ++i1) {
202 // For each start time, we will keep the most violated cut generated while
203 // scanning the residual intervals.
204 int max_violation_end_index = -1;
205 double max_relative_violation = 1.0 + kMinCutViolation;
206 IntegerValue max_violation_window_start(0);
207 IntegerValue max_violation_window_end(0);
208 IntegerValue max_violation_y_min(0);
209 IntegerValue max_violation_y_max(0);
210 IntegerValue max_violation_area(0);
211 bool max_violation_use_precise_area = false;
212
213 // Accumulate intervals, areas, energies and check for potential cuts.
214 double energy_lp = 0.0;
215 IntegerValue window_min = kMaxIntegerValue;
216 IntegerValue window_max = kMinIntegerValue;
217 IntegerValue y_min = kMaxIntegerValue;
218 IntegerValue y_max = kMinIntegerValue;
219 capacity_profile.Clear();
220
221 // We sort all tasks (x_start_min(task) >= x_start_min(start_index) by
222 // increasing end max.
223 std::vector<DiffnEnergyEvent> residual_events(events.begin() + i1,
224 events.end());
225 std::sort(residual_events.begin(), residual_events.end(),
226 [](const DiffnEnergyEvent& a, const DiffnEnergyEvent& b) {
227 return std::tie(a.x_end_max, a.y_spread) <
228 std::tie(b.x_end_max, b.y_spread);
229 });
230 // Let's process residual tasks and evaluate the violation of the cut at
231 // each step. We follow the same structure as the cut creation code below.
232 for (int i2 = 0; i2 < residual_events.size(); ++i2) {
233 const DiffnEnergyEvent& e = residual_events[i2];
234 energy_lp += e.linearized_energy_lp_value;
235 window_min = std::min(window_min, e.x_start_min);
236 window_max = std::max(window_max, e.x_end_max);
237 y_min = std::min(y_min, e.y_min);
238 y_max = std::max(y_max, e.y_max);
239 capacity_profile.AddRectangle(e.x_start_min, e.x_end_max, e.y_min,
240 e.y_max);
241
242 // Dominance rule. If the next interval also fits in
243 // [window_min, window_max]*[y_min, y_max], the cut will be stronger with
244 // the next interval/rectangle.
245 if (i2 + 1 < residual_events.size() &&
246 residual_events[i2 + 1].x_start_min >= window_min &&
247 residual_events[i2 + 1].x_end_max <= window_max &&
248 residual_events[i2 + 1].y_min >= y_min &&
249 residual_events[i2 + 1].y_max <= y_max) {
250 continue;
251 }
252
253 // Checks the current area vs the sum of all energies.
254 // The area is capacity_profile.GetBoundingArea().
255 // We can compare it to the bounding box area:
256 // (window_max - window_min) * (y_max - y_min).
257 bool use_precise_area = false;
258 IntegerValue precise_area(0);
259 double area_lp = 0.0;
260 const IntegerValue bbox_area =
261 (window_max - window_min) * (y_max - y_min);
262 precise_area = capacity_profile.GetBoundingArea();
263 use_precise_area = precise_area < bbox_area;
264 area_lp = ToDouble(std::min(precise_area, bbox_area));
265
266 if (area_lp >= sum_of_all_energies) {
267 break;
268 }
269
270 // Compute the violation of the potential cut.
271 const double relative_violation = energy_lp / area_lp;
272 if (relative_violation > max_relative_violation) {
273 max_violation_end_index = i2;
274 max_relative_violation = relative_violation;
275 max_violation_window_start = window_min;
276 max_violation_window_end = window_max;
277 max_violation_y_min = y_min;
278 max_violation_y_max = y_max;
279 max_violation_area = std::min(precise_area, bbox_area);
280 max_violation_use_precise_area = use_precise_area;
281 }
282 }
283
284 if (max_violation_end_index == -1) continue;
285
286 // A maximal violated cut has been found.
287 // Build it and add it to the pool.
288 bool add_opt_to_name = false;
289 bool add_quadratic_to_name = false;
290 bool add_energy_to_name = false;
291 LinearConstraintBuilder cut(model, kMinIntegerValue, max_violation_area);
292 for (int i2 = 0; i2 <= max_violation_end_index; ++i2) {
293 const DiffnEnergyEvent& event = residual_events[i2];
294 cut.AddLinearExpression(event.linearized_energy);
295 if (!event.IsPresent()) add_opt_to_name = true;
296 if (event.energy_is_quadratic) add_quadratic_to_name = true;
297 if (event.energy_min > event.x_size_min * event.y_size_min) {
298 add_energy_to_name = true;
299 }
300 }
301 std::string full_name(cut_name);
302 if (add_opt_to_name) full_name.append("_optional");
303 if (add_quadratic_to_name) full_name.append("_quadratic");
304 if (add_energy_to_name) full_name.append("_energy");
305 if (max_violation_use_precise_area) full_name.append("_precise");
306 top_n_cuts.AddCut(cut.Build(), full_name, lp_values);
307 }
308 top_n_cuts.TransferToManager(manager);
309}
310
312 NoOverlap2DConstraintHelper* helper, Model* model) {
313 CutGenerator result;
314 result.only_run_at_level_zero = true;
316 &helper->x_helper(), model, &result.vars,
319 &helper->y_helper(), model, &result.vars,
322 ProductDecomposer* product_decomposer =
324 result.generate_cuts = [helper, model, product_decomposer](
325 LinearConstraintManager* manager) {
326 const int num_rectangles = helper->NumBoxes();
327 std::vector<std::vector<LiteralValueValue>> energies(num_rectangles);
328 // TODO(user): We could compute this once and for all in the helper.
329 for (int i = 0; i < num_rectangles; ++i) {
330 energies[i] = product_decomposer->TryToDecompose(
331 helper->x_helper().Sizes()[i], helper->y_helper().Sizes()[i]);
332 }
333 if (!helper->SynchronizeAndSetDirection(true, true, false)) return false;
334 SchedulingDemandHelper* x_demands_helper = &helper->x_demands_helper();
335 SchedulingDemandHelper* y_demands_helper = &helper->y_demands_helper();
336 if (!x_demands_helper->CacheAllEnergyValues()) return true;
337 if (!y_demands_helper->CacheAllEnergyValues()) return true;
338
339 std::vector<int> rectangles;
340 rectangles.reserve(num_rectangles);
341 for (const auto& component :
343 for (const int rect : component) {
344 rectangles.clear();
345 if (helper->IsAbsent(rect)) continue;
346 // We do not consider rectangles controlled by 2 different unassigned
347 // enforcement literals.
348 if (!helper->x_helper().IsPresent(rect) &&
349 !helper->y_helper().IsPresent(rect) &&
350 helper->x_helper().PresenceLiteral(rect) !=
351 helper->y_helper().PresenceLiteral(rect)) {
352 continue;
353 }
354
355 rectangles.push_back(rect);
356 }
357
358 if (rectangles.size() <= 1) continue;
359
360 GenerateNoOverlap2dEnergyCut(energies, rectangles, "NoOverlap2dXEnergy",
361 model, manager, &helper->x_helper(),
362 &helper->y_helper(), y_demands_helper);
363 GenerateNoOverlap2dEnergyCut(energies, rectangles, "NoOverlap2dYEnergy",
364 model, manager, &helper->y_helper(),
365 &helper->x_helper(), x_demands_helper);
366 }
367 return true;
368 };
369 return result;
370}
371
373 : DiffnBaseEvent(t, x_helper) {}
374
375std::string DiffnCtEvent::DebugString() const {
376 return absl::StrCat("DiffnCtEvent(x_end = ", x_end.DebugString(),
377 ", x_start_min = ", x_start_min.value(),
378 ", x_start_max = ", x_start_max.value(),
379 ", x_size_min = ", x_size_min.value(),
380 ", x_lp_end = ", x_lp_end, ", y_min = ", y_min.value(),
381 ", y_max = ", y_max.value(),
382 ", y_size_min = ", y_size_min.value(),
383 ", energy_min = ", energy_min.value(),
384 ", use_energy = ", use_energy, ", lifted = ", lifted);
385}
386
387// We generate the cut from the Smith's rule from:
388// M. Queyranne, Structure of a simple scheduling polyhedron,
389// Mathematical Programming 58 (1993), 263–285
390//
391// The original cut is:
392// sum(end_min_i * duration_min_i) >=
393// (sum(duration_min_i^2) + sum(duration_min_i)^2) / 2
394//
395// Let's build a figure where each horizontal rectangle represent a task. It
396// ends at the end of the task, and its height is the duration of the task.
397// For a given order, we pack each rectangle to the left while not overlapping,
398// that is one rectangle starts when the previous one ends.
399//
400// e1
401// -----
402// :\ | s1
403// : \| e2
404// -------------
405// :\ |
406// : \ | s2
407// : \| e3
408// ----------------
409// : \| s3
410// ----------------
411//
412// We can notice that the total area is independent of the order of tasks.
413// The first term of the rhs is the area above the diagonal.
414// The second term of the rhs is the area below the diagonal.
415//
416// We apply the following changes (see the code for cumulative constraints):
417// - we strengthen this cuts by noticing that if all tasks starts after S,
418// then replacing end_min_i by (end_min_i - S) is still valid.
419// - we lift rectangles that start before the start of the sequence, but must
420// overlap with it.
421// - we apply the same transformation that was applied to the cumulative
422// constraint to use the no_overlap cut in the no_overlap_2d setting.
423// - we use a limited complexity subset-sum to compute reachable capacity
424// - we look at a set of intervals starting after a given start_min, sorted by
425// relative (end_lp - start_min).
426void GenerateNoOvelap2dCompletionTimeCuts(absl::string_view cut_name,
427 std::vector<DiffnCtEvent> events,
428 bool use_lifting, Model* model,
429 LinearConstraintManager* manager) {
430 TopNCuts top_n_cuts(5);
431
432 // Sort by start min to bucketize by start_min.
433 std::sort(events.begin(), events.end(),
434 [](const DiffnCtEvent& e1, const DiffnCtEvent& e2) {
435 return std::tie(e1.x_start_min, e1.y_size_min, e1.x_lp_end) <
436 std::tie(e2.x_start_min, e2.y_size_min, e2.x_lp_end);
437 });
438 for (int start = 0; start + 1 < events.size(); ++start) {
439 // Skip to the next bucket (of start_min).
440 if (start > 0 &&
441 events[start].x_start_min == events[start - 1].x_start_min) {
442 continue;
443 }
444
445 const IntegerValue sequence_start_min = events[start].x_start_min;
446 std::vector<DiffnCtEvent> residual_tasks(events.begin() + start,
447 events.end());
448
449 // We look at event that start before sequence_start_min, but are forced
450 // to cross this time point. In that case, we replace this event by a
451 // truncated event starting at sequence_start_min. To do this, we reduce
452 // the size_min, align the start_min with the sequence_start_min, and
453 // scale the energy down accordingly.
454 if (use_lifting) {
455 for (int before = 0; before < start; ++before) {
456 if (events[before].x_start_min + events[before].x_size_min >
457 sequence_start_min) {
458 // Build the vector of energies as the vector of sizes.
459 DiffnCtEvent event = events[before]; // Copy.
460 event.lifted = true;
461 event.energy_min = ComputeEnergyMinInWindow(
462 event.x_start_min, event.x_start_max, event.x_end_min,
463 event.x_end_max, event.x_size_min, event.y_size_min,
464 event.decomposed_energy, sequence_start_min, event.x_end_max);
465 event.x_size_min =
466 event.x_size_min + event.x_start_min - sequence_start_min;
467 event.x_start_min = sequence_start_min;
468 if (event.energy_min > event.x_size_min * event.y_size_min) {
469 event.use_energy = true;
470 }
471 DCHECK_GE(event.energy_min, event.x_size_min * event.y_size_min);
472 if (event.energy_min <= 0) continue;
473 residual_tasks.push_back(event);
474 }
475 }
476 }
477
478 std::sort(residual_tasks.begin(), residual_tasks.end(),
479 [](const DiffnCtEvent& e1, const DiffnCtEvent& e2) {
480 return e1.x_lp_end < e2.x_lp_end;
481 });
482
483 // Best cut so far for this loop.
484 int best_end = -1;
485 double best_efficacy = 0.01;
486 IntegerValue best_min_rhs = 0;
487 bool best_use_subset_sum = false;
488
489 // Used in the first term of the rhs of the equation.
490 IntegerValue sum_event_areas = 0;
491 // Used in the second term of the rhs of the equation.
492 IntegerValue sum_energy = 0;
493 // For normalization.
494 IntegerValue sum_square_energy = 0;
495
496 double lp_contrib = 0.0;
497 IntegerValue current_start_min(kMaxIntegerValue);
498 IntegerValue y_min_of_subset = kMaxIntegerValue;
499 IntegerValue y_max_of_subset = kMinIntegerValue;
500 IntegerValue sum_of_y_size_min = 0;
501
502 bool use_dp = true;
504
505 for (int i = 0; i < residual_tasks.size(); ++i) {
506 const DiffnCtEvent& event = residual_tasks[i];
507 DCHECK_GE(event.x_start_min, sequence_start_min);
508 // Make sure we do not overflow.
509 if (!AddTo(event.energy_min, &sum_energy)) break;
510 if (!AddProductTo(event.energy_min, event.x_size_min, &sum_event_areas)) {
511 break;
512 }
513 if (!AddSquareTo(event.energy_min, &sum_square_energy)) break;
514 if (!AddTo(event.y_size_min, &sum_of_y_size_min)) break;
515
516 lp_contrib += event.x_lp_end * ToDouble(event.energy_min);
517 current_start_min = std::min(current_start_min, event.x_start_min);
518
519 // For the capacity, we use the worse |y_max - y_min| and if all the tasks
520 // so far have a fixed demand with a gcd > 1, we can round it down.
521 y_min_of_subset = std::min(y_min_of_subset, event.y_min);
522 y_max_of_subset = std::max(y_max_of_subset, event.y_max);
523 if (!event.y_size_is_fixed) use_dp = false;
524 if (use_dp) {
525 if (i == 0) {
526 dp.Reset((y_max_of_subset - y_min_of_subset).value());
527 } else {
528 // TODO(user): Can we increase the bound dynamically ?
529 if (y_max_of_subset - y_min_of_subset > dp.Bound()) {
530 use_dp = false;
531 }
532 }
533 }
534 if (use_dp) {
535 dp.Add(event.y_size_min.value());
536 }
537
538 const IntegerValue reachable_capacity =
539 use_dp ? IntegerValue(dp.CurrentMax())
540 : y_max_of_subset - y_min_of_subset;
541
542 // If we have not reached capacity, there can be no cuts on ends.
543 if (sum_of_y_size_min <= reachable_capacity) continue;
544
545 // Do we have a violated cut ?
546 const IntegerValue square_sum_energy = CapProdI(sum_energy, sum_energy);
547 if (AtMinOrMaxInt64I(square_sum_energy)) break;
548 const IntegerValue rhs_second_term =
549 CeilRatio(square_sum_energy, reachable_capacity);
550
551 IntegerValue min_rhs = CapAddI(sum_event_areas, rhs_second_term);
552 if (AtMinOrMaxInt64I(min_rhs)) break;
553 min_rhs = CeilRatio(min_rhs, 2);
554
555 // shift contribution by current_start_min.
556 if (!AddProductTo(sum_energy, current_start_min, &min_rhs)) break;
557
558 // The efficacy of the cut is the normalized violation of the above
559 // equation. We will normalize by the sqrt of the sum of squared energies.
560 const double efficacy = (ToDouble(min_rhs) - lp_contrib) /
561 std::sqrt(ToDouble(sum_square_energy));
562
563 // For a given start time, we only keep the best cut.
564 // The reason is that if the cut is strongly violated, we can get a
565 // sequence of violated cuts as we add more tasks. These new cuts will
566 // be less violated, but will not bring anything useful to the LP
567 // relaxation. At the same time, this sequence of cuts can push out
568 // other cuts from a disjoint set of tasks.
569 if (efficacy > best_efficacy) {
570 best_efficacy = efficacy;
571 best_end = i;
572 best_min_rhs = min_rhs;
573 best_use_subset_sum =
574 reachable_capacity < y_max_of_subset - y_min_of_subset;
575 }
576 }
577 if (best_end != -1) {
578 LinearConstraintBuilder cut(model, best_min_rhs, kMaxIntegerValue);
579 bool is_lifted = false;
580 bool add_energy_to_name = false;
581 for (int i = 0; i <= best_end; ++i) {
582 const DiffnCtEvent& event = residual_tasks[i];
583 is_lifted |= event.lifted;
584 add_energy_to_name |= event.use_energy;
585 cut.AddTerm(event.x_end, event.energy_min);
586 }
587 std::string full_name(cut_name);
588 if (is_lifted) full_name.append("_lifted");
589 if (add_energy_to_name) full_name.append("_energy");
590 if (best_use_subset_sum) full_name.append("_subsetsum");
591 top_n_cuts.AddCut(cut.Build(), full_name, manager->LpValues());
592 }
593 }
594 top_n_cuts.TransferToManager(manager);
595}
596
598 NoOverlap2DConstraintHelper* helper, Model* model) {
599 CutGenerator result;
600 result.only_run_at_level_zero = true;
602 &helper->x_helper(), model, &result.vars,
605 &helper->y_helper(), model, &result.vars,
608
609 auto* product_decomposer = model->GetOrCreate<ProductDecomposer>();
610 result.generate_cuts = [helper, product_decomposer,
611 model](LinearConstraintManager* manager) {
612 if (!helper->SynchronizeAndSetDirection()) {
613 return false;
614 }
615
616 const int num_rectangles = helper->NumBoxes();
617 std::vector<int> rectangles;
618 rectangles.reserve(num_rectangles);
619 const SchedulingConstraintHelper* x_helper = &helper->x_helper();
620 const SchedulingConstraintHelper* y_helper = &helper->y_helper();
621 for (const auto& component :
622 helper->connected_components().AsVectorOfSpan()) {
623 rectangles.clear();
624 if (component.size() <= 1) continue;
625 for (int rect : component) {
626 if (!helper->IsPresent(rect)) continue;
627 if (x_helper->SizeMin(rect) == 0 || y_helper->SizeMin(rect) == 0) {
628 continue;
629 }
630 rectangles.push_back(rect);
631 }
632
633 auto generate_cuts = [product_decomposer, manager, model, helper,
634 &rectangles](absl::string_view cut_name) {
635 std::vector<DiffnCtEvent> events;
636
637 const SchedulingConstraintHelper* x_helper = &helper->x_helper();
638 const SchedulingConstraintHelper* y_helper = &helper->y_helper();
639 const auto& lp_values = manager->LpValues();
640 for (const int rect : rectangles) {
641 DiffnCtEvent event(rect, x_helper);
642 event.x_end = x_helper->Ends()[rect];
643 event.x_lp_end = event.x_end.LpValue(lp_values);
644 event.y_min = y_helper->StartMin(rect);
645 event.y_max = y_helper->EndMax(rect);
646 event.y_size_min = y_helper->SizeMin(rect);
647
648 // TODO(user): Use improved energy from demands helper.
649 event.energy_min = event.x_size_min * event.y_size_min;
650 event.decomposed_energy = product_decomposer->TryToDecompose(
651 x_helper->Sizes()[rect], y_helper->Sizes()[rect]);
652 events.push_back(event);
653 }
654
655 GenerateNoOvelap2dCompletionTimeCuts(cut_name, std::move(events),
656 /*use_lifting=*/true, model,
657 manager);
658 };
659
660 if (!helper->SynchronizeAndSetDirection(true, true, false)) {
661 return false;
662 }
663 generate_cuts("NoOverlap2dXCompletionTime");
664 if (!helper->SynchronizeAndSetDirection(true, true, true)) {
665 return false;
666 }
667 generate_cuts("NoOverlap2dYCompletionTime");
668 if (!helper->SynchronizeAndSetDirection(false, false, false)) {
669 return false;
670 }
671 generate_cuts("NoOverlap2dXCompletionTime_mirror");
672 if (!helper->SynchronizeAndSetDirection(false, false, true)) {
673 return false;
674 }
675 generate_cuts("NoOverlap2dYCompletionTime_mirror");
676 }
677 return true;
678 };
679 return result;
680}
681
682} // namespace sat
683} // namespace operations_research
void AddRectangle(IntegerValue x_min, IntegerValue x_max, IntegerValue y_min, IntegerValue y_max)
Adds a rectangle to the current shape.
std::vector< absl::Span< const V > > AsVectorOfSpan() const
Definition util.h:816
void AddQuadraticLowerBound(AffineExpression left, AffineExpression right, IntegerTrail *integer_trail, bool *is_quadratic=nullptr)
ABSL_MUST_USE_RESULT bool AddLiteralTerm(Literal lit, IntegerValue coeff=IntegerValue(1))
void AddLinearExpression(const LinearExpression &expr)
ABSL_MUST_USE_RESULT bool AddDecomposedProduct(absl::Span< const LiteralValueValue > product)
const util_intops::StrongVector< IntegerVariable, double > & LpValues()
To simplify CutGenerator api.
LiteralIndex Index() const
Definition sat_base.h:91
const CompactVectorVector< int > & connected_components() const
bool SynchronizeAndSetDirection(bool x_is_forward_after_swap=true, bool y_is_forward_after_swap=true, bool swap_x_and_y=false)
Helper class to express a product as a linear constraint.
std::vector< LiteralValueValue > TryToDecompose(const AffineExpression &left, const AffineExpression &right)
absl::Span< const AffineExpression > Sizes() const
void TransferToManager(LinearConstraintManager *manager)
Empty the local pool and add all its content to the manager.
void AddCut(LinearConstraint ct, absl::string_view name, const util_intops::StrongVector< IntegerVariable, double > &lp_solution)
Adds a cut to the local pool.
void STLSortAndRemoveDuplicates(T *v, const LessFunc &less_func)
Definition stl_util.h:58
bool AddSquareTo(IntegerValue a, IntegerValue *result)
Computes result += a * a, and return false iff there is an overflow.
bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue *result)
Computes result += a * b, and return false iff there is an overflow.
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
bool AddTo(IntegerValue a, IntegerValue *result)
IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor)
CutGenerator CreateNoOverlap2dEnergyCutGenerator(NoOverlap2DConstraintHelper *helper, Model *model)
const LiteralIndex kNoLiteralIndex(-1)
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue.value())
void GenerateNoOverlap2dEnergyCut(absl::Span< const std::vector< LiteralValueValue > > energies, absl::Span< const int > rectangles, absl::string_view cut_name, Model *model, LinearConstraintManager *manager, SchedulingConstraintHelper *x_helper, SchedulingConstraintHelper *y_helper, SchedulingDemandHelper *y_demands_helper)
void AddIntegerVariableFromIntervals(const SchedulingConstraintHelper *helper, Model *model, std::vector< IntegerVariable > *vars, int mask)
IntegerValue CapAddI(IntegerValue a, IntegerValue b)
CutGenerator CreateNoOverlap2dCompletionTimeCutGenerator(NoOverlap2DConstraintHelper *helper, Model *model)
bool AtMinOrMaxInt64I(IntegerValue t)
IntegerValue CapProdI(IntegerValue a, IntegerValue b)
Overflows and saturated arithmetic.
double ToDouble(IntegerValue value)
IntegerValue ComputeEnergyMinInWindow(IntegerValue start_min, IntegerValue start_max, IntegerValue end_min, IntegerValue end_max, IntegerValue size_min, IntegerValue demand_min, absl::Span< const LiteralValueValue > filtered_energy, IntegerValue window_start, IntegerValue window_end)
In SWIG mode, we don't want anything besides these top-level includes.
absl::AnyInvocable< bool(LinearConstraintManager *manager)> generate_cuts
Definition cuts.h:59
std::vector< IntegerVariable > vars
Definition cuts.h:58
std::vector< LiteralValueValue > decomposed_energy
Definition diffn_cuts.h:74
IntegerValue x_start_min
Cache of the intervals bound on the x direction.
Definition diffn_cuts.h:59
DiffnBaseEvent(int t, const SchedulingConstraintHelper *x_helper)
Definition diffn_cuts.cc:55
IntegerValue y_min
Useful for no_overlap_2d or cumulative.
Definition diffn_cuts.h:65
IntegerValue energy_min
The energy min of this event.
Definition diffn_cuts.h:70
DiffnCtEvent(int t, const SchedulingConstraintHelper *x_helper)
AffineExpression x_end
The lp value of the end of the x interval.
Definition diffn_cuts.h:86
AffineExpression y_size
We need this for linearizing the energy in some cases.
Definition diffn_cuts.cc:68
DiffnEnergyEvent(int t, const SchedulingConstraintHelper *x_helper)
Definition diffn_cuts.cc:64
bool energy_is_quadratic
True if linearized_energy is not exact and a McCormick relaxation.
Definition diffn_cuts.cc:80
LiteralIndex presence_literal_index
If set, this event is optional and its presence is controlled by this.
Definition diffn_cuts.cc:71
double y_spread
Used to minimize the increase on the y axis for rectangles.
Definition diffn_cuts.cc:83
IntegerValue GetMinOverlap(IntegerValue start, IntegerValue end) const
Definition diffn_cuts.cc:93
ABSL_MUST_USE_RESULT bool FillEnergyLp(AffineExpression x_size, const util_intops::StrongVector< IntegerVariable, double > &lp_values, Model *model)