OR-Tools  8.2
cuts.cc
Go to the documentation of this file.
1// Copyright 2010-2018 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
14#include "ortools/sat/cuts.h"
15
16#include <algorithm>
17#include <cmath>
18#include <functional>
19#include <memory>
20#include <utility>
21#include <vector>
22
27#include "ortools/sat/integer.h"
32
33namespace operations_research {
34namespace sat {
35
36namespace {
37
38// Minimum amount of violation of the cut constraint by the solution. This
39// is needed to avoid numerical issues and adding cuts with minor effect.
40const double kMinCutViolation = 1e-4;
41
42// Returns the lp value of a Literal.
43double GetLiteralLpValue(
44 const Literal lit,
46 const IntegerEncoder* encoder) {
47 const IntegerVariable direct_view = encoder->GetLiteralView(lit);
48 if (direct_view != kNoIntegerVariable) {
49 return lp_values[direct_view];
50 }
51 const IntegerVariable opposite_view = encoder->GetLiteralView(lit.Negated());
52 DCHECK_NE(opposite_view, kNoIntegerVariable);
53 return 1.0 - lp_values[opposite_view];
54}
55
56// Returns a constraint that disallow all given variables to be at their current
57// upper bound. The arguments must form a non-trival constraint of the form
58// sum terms (coeff * var) <= upper_bound.
59LinearConstraint GenerateKnapsackCutForCover(
60 const std::vector<IntegerVariable>& vars,
61 const std::vector<IntegerValue>& coeffs, const IntegerValue upper_bound,
62 const IntegerTrail& integer_trail) {
63 CHECK_EQ(vars.size(), coeffs.size());
64 CHECK_GT(vars.size(), 0);
65 LinearConstraint cut;
66 IntegerValue cut_upper_bound = IntegerValue(0);
67 IntegerValue max_coeff = coeffs[0];
68 // slack = \sum_{i}(coeffs[i] * upper_bound[i]) - upper_bound.
69 IntegerValue slack = -upper_bound;
70 for (int i = 0; i < vars.size(); ++i) {
71 const IntegerValue var_upper_bound =
72 integer_trail.LevelZeroUpperBound(vars[i]);
73 cut_upper_bound += var_upper_bound;
74 cut.vars.push_back(vars[i]);
75 cut.coeffs.push_back(IntegerValue(1));
76 max_coeff = std::max(max_coeff, coeffs[i]);
77 slack += coeffs[i] * var_upper_bound;
78 }
79 CHECK_GT(slack, 0.0) << "Invalid cover for knapsack cut.";
80 cut_upper_bound -= CeilRatio(slack, max_coeff);
81 cut.lb = kMinIntegerValue;
82 cut.ub = cut_upper_bound;
83 VLOG(2) << "Generated Knapsack Constraint:" << cut.DebugString();
84 return cut;
85}
86
87bool SolutionSatisfiesConstraint(
88 const LinearConstraint& constraint,
90 const double activity = ComputeActivity(constraint, lp_values);
91 const double tolerance = 1e-6;
92 return (activity <= constraint.ub.value() + tolerance &&
93 activity >= constraint.lb.value() - tolerance)
94 ? true
95 : false;
96}
97
98bool SmallRangeAndAllCoefficientsMagnitudeAreTheSame(
99 const LinearConstraint& constraint, IntegerTrail* integer_trail) {
100 if (constraint.vars.empty()) return true;
101
102 const int64 magnitude = std::abs(constraint.coeffs[0].value());
103 for (int i = 1; i < constraint.coeffs.size(); ++i) {
104 const IntegerVariable var = constraint.vars[i];
105 if (integer_trail->LevelZeroUpperBound(var) -
106 integer_trail->LevelZeroLowerBound(var) >
107 1) {
108 return false;
109 }
110 if (std::abs(constraint.coeffs[i].value()) != magnitude) {
111 return false;
112 }
113 }
114 return true;
115}
116
117bool AllVarsTakeIntegerValue(
118 const std::vector<IntegerVariable> vars,
120 for (IntegerVariable var : vars) {
121 if (std::abs(lp_values[var] - std::round(lp_values[var])) > 1e-6) {
122 return false;
123 }
124 }
125 return true;
126}
127
128// Returns smallest cover size for the given constraint taking into account
129// level zero bounds. Smallest Cover size is computed as follows.
130// 1. Compute the upper bound if all variables are shifted to have zero lower
131// bound.
132// 2. Sort all terms (coefficient * shifted upper bound) in non decreasing
133// order.
134// 3. Add terms in cover until term sum is smaller or equal to upper bound.
135// 4. Add the last item which violates the upper bound. This forms the smallest
136// cover. Return the size of this cover.
137int GetSmallestCoverSize(const LinearConstraint& constraint,
138 const IntegerTrail& integer_trail) {
139 IntegerValue ub = constraint.ub;
140 std::vector<IntegerValue> sorted_terms;
141 for (int i = 0; i < constraint.vars.size(); ++i) {
142 const IntegerValue coeff = constraint.coeffs[i];
143 const IntegerVariable var = constraint.vars[i];
144 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var);
145 const IntegerValue var_lb = integer_trail.LevelZeroLowerBound(var);
146 ub -= var_lb * coeff;
147 sorted_terms.push_back(coeff * (var_ub - var_lb));
148 }
149 std::sort(sorted_terms.begin(), sorted_terms.end(),
150 std::greater<IntegerValue>());
151 int smallest_cover_size = 0;
152 IntegerValue sorted_term_sum = IntegerValue(0);
153 while (sorted_term_sum <= ub &&
154 smallest_cover_size < constraint.vars.size()) {
155 sorted_term_sum += sorted_terms[smallest_cover_size++];
156 }
157 return smallest_cover_size;
158}
159
160bool ConstraintIsEligibleForLifting(const LinearConstraint& constraint,
161 const IntegerTrail& integer_trail) {
162 for (const IntegerVariable var : constraint.vars) {
163 if (integer_trail.LevelZeroLowerBound(var) != IntegerValue(0) ||
164 integer_trail.LevelZeroUpperBound(var) != IntegerValue(1)) {
165 return false;
166 }
167 }
168 return true;
169}
170} // namespace
171
173 const LinearConstraint& constraint,
175 const std::vector<IntegerValue>& cut_vars_original_coefficients,
176 const IntegerTrail& integer_trail, TimeLimit* time_limit,
177 LinearConstraint* cut) {
178 std::set<IntegerVariable> vars_in_cut;
179 for (IntegerVariable var : cut->vars) {
180 vars_in_cut.insert(var);
181 }
182
183 std::vector<std::pair<IntegerValue, IntegerVariable>> non_zero_vars;
184 std::vector<std::pair<IntegerValue, IntegerVariable>> zero_vars;
185 for (int i = 0; i < constraint.vars.size(); ++i) {
186 const IntegerVariable var = constraint.vars[i];
187 if (integer_trail.LevelZeroLowerBound(var) != IntegerValue(0) ||
188 integer_trail.LevelZeroUpperBound(var) != IntegerValue(1)) {
189 continue;
190 }
191 if (vars_in_cut.find(var) != vars_in_cut.end()) continue;
192 const IntegerValue coeff = constraint.coeffs[i];
193 if (lp_values[var] <= 1e-6) {
194 zero_vars.push_back({coeff, var});
195 } else {
196 non_zero_vars.push_back({coeff, var});
197 }
198 }
199
200 // Decide lifting sequence (nonzeros, zeros in nonincreasing order
201 // of coefficient ).
202 std::sort(non_zero_vars.rbegin(), non_zero_vars.rend());
203 std::sort(zero_vars.rbegin(), zero_vars.rend());
204
205 std::vector<std::pair<IntegerValue, IntegerVariable>> lifting_sequence(
206 std::move(non_zero_vars));
207
208 lifting_sequence.insert(lifting_sequence.end(), zero_vars.begin(),
209 zero_vars.end());
210
211 // Form Knapsack.
212 std::vector<double> lifting_profits;
213 std::vector<double> lifting_weights;
214 for (int i = 0; i < cut->vars.size(); ++i) {
215 lifting_profits.push_back(cut->coeffs[i].value());
216 lifting_weights.push_back(cut_vars_original_coefficients[i].value());
217 }
218
219 // Lift the cut.
220 bool is_lifted = false;
221 bool is_solution_optimal = false;
222 KnapsackSolverForCuts knapsack_solver("Knapsack cut lifter");
223 for (auto entry : lifting_sequence) {
224 is_solution_optimal = false;
225 const IntegerValue var_original_coeff = entry.first;
226 const IntegerVariable var = entry.second;
227 const IntegerValue lifting_capacity = constraint.ub - entry.first;
228 if (lifting_capacity <= IntegerValue(0)) continue;
229 knapsack_solver.Init(lifting_profits, lifting_weights,
230 lifting_capacity.value());
231 knapsack_solver.set_node_limit(100);
232 // NOTE: Since all profits and weights are integer, solution of
233 // knapsack is also integer.
234 // TODO(user): Use an integer solver or heuristic.
235 knapsack_solver.Solve(time_limit, &is_solution_optimal);
236 const double knapsack_upper_bound =
237 std::round(knapsack_solver.GetUpperBound());
238 const IntegerValue cut_coeff = cut->ub - knapsack_upper_bound;
239 if (cut_coeff > IntegerValue(0)) {
240 is_lifted = true;
241 cut->vars.push_back(var);
242 cut->coeffs.push_back(cut_coeff);
243 lifting_profits.push_back(cut_coeff.value());
244 lifting_weights.push_back(var_original_coeff.value());
245 }
246 }
247 return is_lifted;
248}
249
251 const LinearConstraint& constraint,
253 const IntegerTrail& integer_trail) {
254 IntegerValue ub = constraint.ub;
255 LinearConstraint constraint_with_left_vars;
256 for (int i = 0; i < constraint.vars.size(); ++i) {
257 const IntegerVariable var = constraint.vars[i];
258 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var);
259 const IntegerValue coeff = constraint.coeffs[i];
260 if (var_ub.value() - lp_values[var] <= 1.0 - kMinCutViolation) {
261 constraint_with_left_vars.vars.push_back(var);
262 constraint_with_left_vars.coeffs.push_back(coeff);
263 } else {
264 // Variable not in cut
265 const IntegerValue var_lb = integer_trail.LevelZeroLowerBound(var);
266 ub -= coeff * var_lb;
267 }
268 }
269 constraint_with_left_vars.ub = ub;
270 constraint_with_left_vars.lb = constraint.lb;
271 return constraint_with_left_vars;
272}
273
275 const IntegerTrail& integer_trail) {
276 IntegerValue term_sum = IntegerValue(0);
277 for (int i = 0; i < constraint.vars.size(); ++i) {
278 const IntegerVariable var = constraint.vars[i];
279 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var);
280 const IntegerValue coeff = constraint.coeffs[i];
281 term_sum += coeff * var_ub;
282 }
283 if (term_sum <= constraint.ub) {
284 VLOG(2) << "Filtered by cover filter";
285 return true;
286 }
287 return false;
288}
289
291 const LinearConstraint& preprocessed_constraint,
293 const IntegerTrail& integer_trail) {
294 std::vector<double> variable_upper_bound_distances;
295 for (const IntegerVariable var : preprocessed_constraint.vars) {
296 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var);
297 variable_upper_bound_distances.push_back(var_ub.value() - lp_values[var]);
298 }
299 // Compute the min cover size.
300 const int smallest_cover_size =
301 GetSmallestCoverSize(preprocessed_constraint, integer_trail);
302
303 std::nth_element(
304 variable_upper_bound_distances.begin(),
305 variable_upper_bound_distances.begin() + smallest_cover_size - 1,
306 variable_upper_bound_distances.end());
307 double cut_lower_bound = 0.0;
308 for (int i = 0; i < smallest_cover_size; ++i) {
309 cut_lower_bound += variable_upper_bound_distances[i];
310 }
311 if (cut_lower_bound >= 1.0 - kMinCutViolation) {
312 VLOG(2) << "Filtered by kappa heuristic";
313 return true;
314 }
315 return false;
316}
317
318double GetKnapsackUpperBound(std::vector<KnapsackItem> items,
319 const double capacity) {
320 // Sort items by value by weight ratio.
321 std::sort(items.begin(), items.end(), std::greater<KnapsackItem>());
322 double left_capacity = capacity;
323 double profit = 0.0;
324 for (const KnapsackItem item : items) {
325 if (item.weight <= left_capacity) {
326 profit += item.profit;
327 left_capacity -= item.weight;
328 } else {
329 profit += (left_capacity / item.weight) * item.profit;
330 break;
331 }
332 }
333 return profit;
334}
335
337 const LinearConstraint& constraint,
339 const IntegerTrail& integer_trail) {
340 std::vector<KnapsackItem> items;
341 double capacity = -constraint.ub.value() - 1.0;
342 double sum_variable_profit = 0;
343 for (int i = 0; i < constraint.vars.size(); ++i) {
344 const IntegerVariable var = constraint.vars[i];
345 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var);
346 const IntegerValue var_lb = integer_trail.LevelZeroLowerBound(var);
347 const IntegerValue coeff = constraint.coeffs[i];
348 KnapsackItem item;
349 item.profit = var_ub.value() - lp_values[var];
350 item.weight = (coeff * (var_ub - var_lb)).value();
351 items.push_back(item);
352 capacity += (coeff * var_ub).value();
353 sum_variable_profit += item.profit;
354 }
355
356 // Return early if the required upper bound is negative since all the profits
357 // are non negative.
358 if (sum_variable_profit - 1.0 + kMinCutViolation < 0.0) return false;
359
360 // Get the knapsack upper bound.
361 const double knapsack_upper_bound =
362 GetKnapsackUpperBound(std::move(items), capacity);
363 if (knapsack_upper_bound < sum_variable_profit - 1.0 + kMinCutViolation) {
364 VLOG(2) << "Filtered by knapsack upper bound";
365 return true;
366 }
367 return false;
368}
369
371 const LinearConstraint& preprocessed_constraint,
373 const IntegerTrail& integer_trail) {
374 if (ConstraintIsTriviallyTrue(preprocessed_constraint, integer_trail)) {
375 return false;
376 }
377 if (CanBeFilteredUsingCutLowerBound(preprocessed_constraint, lp_values,
378 integer_trail)) {
379 return false;
380 }
381 if (CanBeFilteredUsingKnapsackUpperBound(preprocessed_constraint, lp_values,
382 integer_trail)) {
383 return false;
384 }
385 return true;
386}
387
389 std::vector<LinearConstraint>* knapsack_constraints,
390 IntegerTrail* integer_trail) {
391 // If all coefficient are the same, the generated knapsack cuts cannot be
392 // stronger than the constraint itself. However, when we substitute variables
393 // using the implication graph, this is not longer true. So we only skip
394 // constraints with same coeff and no substitutions.
395 if (SmallRangeAndAllCoefficientsMagnitudeAreTheSame(constraint,
396 integer_trail)) {
397 return;
398 }
399 if (constraint.ub < kMaxIntegerValue) {
400 LinearConstraint canonical_knapsack_form;
401
402 // Negate the variables with negative coefficients.
403 for (int i = 0; i < constraint.vars.size(); ++i) {
404 const IntegerVariable var = constraint.vars[i];
405 const IntegerValue coeff = constraint.coeffs[i];
406 if (coeff > IntegerValue(0)) {
407 canonical_knapsack_form.AddTerm(var, coeff);
408 } else {
409 canonical_knapsack_form.AddTerm(NegationOf(var), -coeff);
410 }
411 }
412 canonical_knapsack_form.ub = constraint.ub;
413 canonical_knapsack_form.lb = kMinIntegerValue;
414 knapsack_constraints->push_back(canonical_knapsack_form);
415 }
416
417 if (constraint.lb > kMinIntegerValue) {
418 LinearConstraint canonical_knapsack_form;
419
420 // Negate the variables with positive coefficients.
421 for (int i = 0; i < constraint.vars.size(); ++i) {
422 const IntegerVariable var = constraint.vars[i];
423 const IntegerValue coeff = constraint.coeffs[i];
424 if (coeff > IntegerValue(0)) {
425 canonical_knapsack_form.AddTerm(NegationOf(var), coeff);
426 } else {
427 canonical_knapsack_form.AddTerm(var, -coeff);
428 }
429 }
430 canonical_knapsack_form.ub = -constraint.lb;
431 canonical_knapsack_form.lb = kMinIntegerValue;
432 knapsack_constraints->push_back(canonical_knapsack_form);
433 }
434}
435
436// TODO(user): Move the cut generator into a class and reuse variables.
438 const std::vector<LinearConstraint>& base_constraints,
439 const std::vector<IntegerVariable>& vars, Model* model) {
440 CutGenerator result;
441 result.vars = vars;
442
443 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
444 std::vector<LinearConstraint> knapsack_constraints;
445 for (const LinearConstraint& constraint : base_constraints) {
446 // There is often a lot of small linear base constraints and it doesn't seem
447 // super useful to generate cuts for constraints of size 2. Any valid cut
448 // of size 1 should be already infered by the propagation.
449 //
450 // TODO(user): The case of size 2 is a bit less clear. investigate more if
451 // it is useful.
452 if (constraint.vars.size() <= 2) continue;
453
454 ConvertToKnapsackForm(constraint, &knapsack_constraints, integer_trail);
455 }
456 VLOG(1) << "#knapsack constraints: " << knapsack_constraints.size();
457
458 // Note(user): for Knapsack cuts, it seems always advantageous to replace a
459 // variable X by a TIGHT lower bound of the form "coeff * binary + lb". This
460 // will not change "covers" but can only result in more violation by the
461 // current LP solution.
462 ImpliedBoundsProcessor implied_bounds_processor(
463 vars, integer_trail, model->GetOrCreate<ImpliedBounds>());
464
465 // TODO(user): do not add generator if there are no knapsack constraints.
466 result.generate_cuts = [implied_bounds_processor, knapsack_constraints, vars,
467 model, integer_trail](
469 lp_values,
470 LinearConstraintManager* manager) mutable {
471 // TODO(user): When we use implied-bound substitution, we might still infer
472 // an interesting cut even if all variables are integer. See if we still
473 // want to skip all such constraints.
474 if (AllVarsTakeIntegerValue(vars, lp_values)) return;
475
476 KnapsackSolverForCuts knapsack_solver(
477 "Knapsack on demand cover cut generator");
478 int64 skipped_constraints = 0;
479 LinearConstraint mutable_constraint;
480
481 // Iterate through all knapsack constraints.
482 implied_bounds_processor.ClearCache();
483 for (const LinearConstraint& constraint : knapsack_constraints) {
484 if (model->GetOrCreate<TimeLimit>()->LimitReached()) break;
485 VLOG(2) << "Processing constraint: " << constraint.DebugString();
486
487 mutable_constraint = constraint;
488 implied_bounds_processor.ProcessUpperBoundedConstraint(
489 lp_values, &mutable_constraint);
490 MakeAllCoefficientsPositive(&mutable_constraint);
491
492 const LinearConstraint preprocessed_constraint =
493 GetPreprocessedLinearConstraint(mutable_constraint, lp_values,
494 *integer_trail);
495 if (preprocessed_constraint.vars.empty()) continue;
496
497 if (!CanFormValidKnapsackCover(preprocessed_constraint, lp_values,
498 *integer_trail)) {
499 skipped_constraints++;
500 continue;
501 }
502
503 // Profits are (upper_bounds[i] - lp_values[i]) for knapsack variables.
504 std::vector<double> profits;
505 profits.reserve(preprocessed_constraint.vars.size());
506
507 // Weights are (coeffs[i] * (upper_bound[i] - lower_bound[i])).
508 std::vector<double> weights;
509 weights.reserve(preprocessed_constraint.vars.size());
510
511 double capacity = -preprocessed_constraint.ub.value() - 1.0;
512
513 // Compute and store the sum of variable profits. This is the constant
514 // part of the objective of the problem we are trying to solve. Hence
515 // this part is not supplied to the knapsack_solver and is subtracted
516 // when we receive the knapsack solution.
517 double sum_variable_profit = 0;
518
519 // Compute the profits, the weights and the capacity for the knapsack
520 // instance.
521 for (int i = 0; i < preprocessed_constraint.vars.size(); ++i) {
522 const IntegerVariable var = preprocessed_constraint.vars[i];
523 const double coefficient = preprocessed_constraint.coeffs[i].value();
524 const double var_ub = ToDouble(integer_trail->LevelZeroUpperBound(var));
525 const double var_lb = ToDouble(integer_trail->LevelZeroLowerBound(var));
526 const double variable_profit = var_ub - lp_values[var];
527 profits.push_back(variable_profit);
528
529 sum_variable_profit += variable_profit;
530
531 const double weight = coefficient * (var_ub - var_lb);
532 weights.push_back(weight);
533 capacity += weight + coefficient * var_lb;
534 }
535 if (capacity < 0.0) continue;
536
537 std::vector<IntegerVariable> cut_vars;
538 std::vector<IntegerValue> cut_vars_original_coefficients;
539
540 VLOG(2) << "Knapsack size: " << profits.size();
541 knapsack_solver.Init(profits, weights, capacity);
542
543 // Set the time limit for the knapsack solver.
544 const double time_limit_for_knapsack_solver =
545 model->GetOrCreate<TimeLimit>()->GetTimeLeft();
546
547 // Solve the instance and subtract the constant part to compute the
548 // sum_of_distance_to_ub_for_vars_in_cover.
549 // TODO(user): Consider solving the instance approximately.
550 bool is_solution_optimal = false;
552 sum_variable_profit - 1.0 + kMinCutViolation);
553 // TODO(user): Consider providing lower bound threshold as
554 // sum_variable_profit - 1.0 + kMinCutViolation.
555 // TODO(user): Set node limit for knapsack solver.
556 auto time_limit_for_solver =
557 absl::make_unique<TimeLimit>(time_limit_for_knapsack_solver);
558 const double sum_of_distance_to_ub_for_vars_in_cover =
559 sum_variable_profit -
560 knapsack_solver.Solve(time_limit_for_solver.get(),
561 &is_solution_optimal);
562 if (is_solution_optimal) {
563 VLOG(2) << "Knapsack Optimal solution found yay !";
564 }
565 if (time_limit_for_solver->LimitReached()) {
566 VLOG(1) << "Knapsack Solver run out of time limit.";
567 }
568 if (sum_of_distance_to_ub_for_vars_in_cover < 1.0 - kMinCutViolation) {
569 // Constraint is eligible for the cover.
570
571 IntegerValue constraint_ub_for_cut = preprocessed_constraint.ub;
572 std::set<IntegerVariable> vars_in_cut;
573 for (int i = 0; i < preprocessed_constraint.vars.size(); ++i) {
574 const IntegerVariable var = preprocessed_constraint.vars[i];
575 const IntegerValue coefficient = preprocessed_constraint.coeffs[i];
576 if (!knapsack_solver.best_solution(i)) {
577 cut_vars.push_back(var);
578 cut_vars_original_coefficients.push_back(coefficient);
579 vars_in_cut.insert(var);
580 } else {
581 const IntegerValue var_lb = integer_trail->LevelZeroLowerBound(var);
582 constraint_ub_for_cut -= coefficient * var_lb;
583 }
584 }
585 LinearConstraint cut = GenerateKnapsackCutForCover(
586 cut_vars, cut_vars_original_coefficients, constraint_ub_for_cut,
587 *integer_trail);
588
589 // Check if the constraint has only binary variables.
590 bool is_lifted = false;
591 if (ConstraintIsEligibleForLifting(cut, *integer_trail)) {
592 if (LiftKnapsackCut(mutable_constraint, lp_values,
593 cut_vars_original_coefficients, *integer_trail,
594 model->GetOrCreate<TimeLimit>(), &cut)) {
595 is_lifted = true;
596 }
597 }
598
599 CHECK(!SolutionSatisfiesConstraint(cut, lp_values));
600 manager->AddCut(cut, is_lifted ? "LiftedKnapsack" : "Knapsack",
601 lp_values);
602 }
603 }
604 if (skipped_constraints > 0) {
605 VLOG(2) << "Skipped constraints: " << skipped_constraints;
606 }
607 };
608
609 return result;
610}
611
612// Compute the larger t <= max_t such that t * rhs_remainder >= divisor / 2.
613//
614// This is just a separate function as it is slightly faster to compute the
615// result only once.
616IntegerValue GetFactorT(IntegerValue rhs_remainder, IntegerValue divisor,
617 IntegerValue max_t) {
618 DCHECK_GE(max_t, 1);
619 return rhs_remainder == 0
620 ? max_t
621 : std::min(max_t, CeilRatio(divisor / 2, rhs_remainder));
622}
623
624std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(
625 IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue t,
626 IntegerValue max_scaling) {
627 DCHECK_GE(max_scaling, 1);
628
629 // Adjust after the multiplication by t.
630 rhs_remainder *= t;
631 DCHECK_LT(rhs_remainder, divisor);
632
633 // Make sure we don't have an integer overflow below. Note that we assume that
634 // divisor and the maximum coeff magnitude are not too different (maybe a
635 // factor 1000 at most) so that the final result will never overflow.
636 max_scaling = std::min(max_scaling, kint64max / divisor);
637
638 const IntegerValue size = divisor - rhs_remainder;
639 if (max_scaling == 1 || size == 1) {
640 // TODO(user): Use everywhere a two step computation to avoid overflow?
641 // First divide by divisor, then multiply by t. For now, we limit t so that
642 // we never have an overflow instead.
643 return [t, divisor](IntegerValue coeff) {
644 return FloorRatio(t * coeff, divisor);
645 };
646 } else if (size <= max_scaling) {
647 return [size, rhs_remainder, t, divisor](IntegerValue coeff) {
648 const IntegerValue ratio = FloorRatio(t * coeff, divisor);
649 const IntegerValue remainder = t * coeff - ratio * divisor;
650 const IntegerValue diff = remainder - rhs_remainder;
651 return size * ratio + std::max(IntegerValue(0), diff);
652 };
653 } else if (max_scaling.value() * rhs_remainder.value() < divisor) {
654 // Because of our max_t limitation, the rhs_remainder might stay small.
655 //
656 // If it is "too small" we cannot use the code below because it will not be
657 // valid. So we just divide divisor into max_scaling bucket. The
658 // rhs_remainder will be in the bucket 0.
659 //
660 // Note(user): This seems the same as just increasing t, modulo integer
661 // overflows. Maybe we should just always do the computation like this so
662 // that we can use larger t even if coeff is close to kint64max.
663 return [t, divisor, max_scaling](IntegerValue coeff) {
664 const IntegerValue ratio = FloorRatio(t * coeff, divisor);
665 const IntegerValue remainder = t * coeff - ratio * divisor;
666 const IntegerValue bucket = FloorRatio(remainder * max_scaling, divisor);
667 return max_scaling * ratio + bucket;
668 };
669 } else {
670 // We divide (size = divisor - rhs_remainder) into (max_scaling - 1) buckets
671 // and increase the function by 1 / max_scaling for each of them.
672 //
673 // Note that for different values of max_scaling, we get a family of
674 // functions that do not dominate each others. So potentially, a max scaling
675 // as low as 2 could lead to the better cut (this is exactly the Letchford &
676 // Lodi function).
677 //
678 // Another intersting fact, is that if we want to compute the maximum alpha
679 // for a constraint with 2 terms like:
680 // divisor * Y + (ratio * divisor + remainder) * X
681 // <= rhs_ratio * divisor + rhs_remainder
682 // so that we have the cut:
683 // Y + (ratio + alpha) * X <= rhs_ratio
684 // This is the same as computing the maximum alpha such that for all integer
685 // X > 0 we have CeilRatio(alpha * divisor * X, divisor)
686 // <= CeilRatio(remainder * X - rhs_remainder, divisor).
687 // We can prove that this alpha is of the form (n - 1) / n, and it will
688 // be reached by such function for a max_scaling of n.
689 //
690 // TODO(user): This function is not always maximal when
691 // size % (max_scaling - 1) == 0. Improve?
692 return [size, rhs_remainder, t, divisor, max_scaling](IntegerValue coeff) {
693 const IntegerValue ratio = FloorRatio(t * coeff, divisor);
694 const IntegerValue remainder = t * coeff - ratio * divisor;
695 const IntegerValue diff = remainder - rhs_remainder;
696 const IntegerValue bucket =
697 diff > 0 ? CeilRatio(diff * (max_scaling - 1), size)
698 : IntegerValue(0);
699 return max_scaling * ratio + bucket;
700 };
701 }
702}
703
704// TODO(user): This has been optimized a bit, but we can probably do even better
705// as it still takes around 25% percent of the run time when all the cuts are on
706// for the opm*mps.gz problems and others.
708 RoundingOptions options, const std::vector<double>& lp_values,
709 const std::vector<IntegerValue>& lower_bounds,
710 const std::vector<IntegerValue>& upper_bounds,
711 ImpliedBoundsProcessor* ib_processor, LinearConstraint* cut) {
712 const int size = lp_values.size();
713 if (size == 0) return;
714 CHECK_EQ(lower_bounds.size(), size);
715 CHECK_EQ(upper_bounds.size(), size);
716 CHECK_EQ(cut->vars.size(), size);
717 CHECK_EQ(cut->coeffs.size(), size);
719
720 // To optimize the computation of the best divisor below, we only need to
721 // look at the indices with a shifted lp value that is not close to zero.
722 //
723 // TODO(user): sort by decreasing lp_values so that our early abort test in
724 // the critical loop below has more chance of returning early? I tried but it
725 // didn't seems to change much though.
726 relevant_indices_.clear();
727 relevant_lp_values_.clear();
728 relevant_coeffs_.clear();
729 relevant_bound_diffs_.clear();
730 divisors_.clear();
731 adjusted_coeffs_.clear();
732
733 // Compute the maximum magnitude for non-fixed variables.
734 IntegerValue max_magnitude(0);
735 for (int i = 0; i < size; ++i) {
736 if (lower_bounds[i] == upper_bounds[i]) continue;
737 const IntegerValue magnitude = IntTypeAbs(cut->coeffs[i]);
738 max_magnitude = std::max(max_magnitude, magnitude);
739 }
740
741 // Shift each variable using its lower/upper bound so that no variable can
742 // change sign. We eventually do a change of variable to its negation so
743 // that all variable are non-negative.
744 bool overflow = false;
745 change_sign_at_postprocessing_.assign(size, false);
746 for (int i = 0; i < size; ++i) {
747 if (cut->coeffs[i] == 0) continue;
748 const IntegerValue magnitude = IntTypeAbs(cut->coeffs[i]);
749
750 // We might change them below.
751 IntegerValue lb = lower_bounds[i];
752 double lp_value = lp_values[i];
753
754 const IntegerValue ub = upper_bounds[i];
755 const IntegerValue bound_diff =
756 IntegerValue(CapSub(ub.value(), lb.value()));
757
758 // Note that since we use ToDouble() this code works fine with lb/ub at
759 // min/max integer value.
760 //
761 // TODO(user): Experiments with different heuristics. Other solver also
762 // seems to try a bunch of possibilities in a "postprocess" phase once
763 // the divisor is chosen. Try that.
764 {
765 // when the magnitude of the entry become smaller and smaller we bias
766 // towards a positive coefficient. This is because after rounding this
767 // will likely become zero instead of -divisor and we need the lp value
768 // to be really close to its bound to compensate.
769 const double lb_dist = std::abs(lp_value - ToDouble(lb));
770 const double ub_dist = std::abs(lp_value - ToDouble(ub));
771 const double bias =
772 std::max(1.0, 0.1 * ToDouble(max_magnitude) / ToDouble(magnitude));
773 if ((bias * lb_dist > ub_dist && cut->coeffs[i] < 0) ||
774 (lb_dist > bias * ub_dist && cut->coeffs[i] > 0)) {
775 change_sign_at_postprocessing_[i] = true;
776 cut->coeffs[i] = -cut->coeffs[i];
777 lp_value = -lp_value;
778 lb = -ub;
779 }
780 }
781
782 // Always shift to lb.
783 // coeff * X = coeff * (X - shift) + coeff * shift.
784 lp_value -= ToDouble(lb);
785 if (!AddProductTo(-cut->coeffs[i], lb, &cut->ub)) {
786 overflow = true;
787 break;
788 }
789
790 // Deal with fixed variable, no need to shift back in this case, we can
791 // just remove the term.
792 if (bound_diff == 0) {
793 cut->coeffs[i] = IntegerValue(0);
794 lp_value = 0.0;
795 }
796
797 if (std::abs(lp_value) > 1e-2) {
798 relevant_coeffs_.push_back(cut->coeffs[i]);
799 relevant_indices_.push_back(i);
800 relevant_lp_values_.push_back(lp_value);
801 relevant_bound_diffs_.push_back(bound_diff);
802 divisors_.push_back(magnitude);
803 }
804 }
805
806 // TODO(user): Maybe this shouldn't be called on such constraint.
807 if (relevant_coeffs_.empty()) {
808 VLOG(2) << "Issue, nothing to cut.";
809 *cut = LinearConstraint(IntegerValue(0), IntegerValue(0));
810 return;
811 }
812 CHECK_NE(max_magnitude, 0);
813
814 // Our heuristic will try to generate a few different cuts, and we will keep
815 // the most violated one scaled by the l2 norm of the relevant position.
816 //
817 // TODO(user): Experiment for the best value of this initial violation
818 // threshold. Note also that we use the l2 norm on the restricted position
819 // here. Maybe we should change that? On that note, the L2 norm usage seems a
820 // bit weird to me since it grows with the number of term in the cut. And
821 // often, we already have a good cut, and we make it stronger by adding extra
822 // terms that do not change its activity.
823 //
824 // The discussion above only concern the best_scaled_violation initial value.
825 // The remainder_threshold allows to not consider cuts for which the final
826 // efficacity is clearly lower than 1e-3 (it is a bound, so we could generate
827 // cuts with a lower efficacity than this).
828 double best_scaled_violation = 0.01;
829 const IntegerValue remainder_threshold(max_magnitude / 1000);
830
831 // The cut->ub might have grown quite a bit with the bound substitution, so
832 // we need to include it too since we will apply the rounding function on it.
833 max_magnitude = std::max(max_magnitude, IntTypeAbs(cut->ub));
834
835 // Make sure that when we multiply the rhs or the coefficient by a factor t,
836 // we do not have an integer overflow. Actually, we need a bit more room
837 // because we might round down a value to the next multiple of
838 // max_magnitude.
839 const IntegerValue threshold = kMaxIntegerValue / 2;
840 if (overflow || max_magnitude >= threshold) {
841 VLOG(2) << "Issue, overflow.";
842 *cut = LinearConstraint(IntegerValue(0), IntegerValue(0));
843 return;
844 }
845 const IntegerValue max_t = threshold / max_magnitude;
846
847 // There is no point trying twice the same divisor or a divisor that is too
848 // small. Note that we use a higher threshold than the remainder_threshold
849 // because we can boost the remainder thanks to our adjusting heuristic below
850 // and also because this allows to have cuts with a small range of
851 // coefficients.
852 //
853 // TODO(user): Note that the std::sort() is visible in some cpu profile.
854 {
855 int new_size = 0;
856 const IntegerValue divisor_threshold = max_magnitude / 10;
857 for (int i = 0; i < divisors_.size(); ++i) {
858 if (divisors_[i] <= divisor_threshold) continue;
859 divisors_[new_size++] = divisors_[i];
860 }
861 divisors_.resize(new_size);
862 }
863 gtl::STLSortAndRemoveDuplicates(&divisors_, std::greater<IntegerValue>());
864
865 // TODO(user): Avoid quadratic algorithm? Note that we are quadratic in
866 // relevant_indices not the full cut->coeffs.size(), but this is still too
867 // much on some problems.
868 IntegerValue best_divisor(0);
869 for (const IntegerValue divisor : divisors_) {
870 // Skip if we don't have the potential to generate a good enough cut.
871 const IntegerValue initial_rhs_remainder =
872 cut->ub - FloorRatio(cut->ub, divisor) * divisor;
873 if (initial_rhs_remainder <= remainder_threshold) continue;
874
875 IntegerValue temp_ub = cut->ub;
876 adjusted_coeffs_.clear();
877
878 // We will adjust coefficient that are just under an exact multiple of
879 // divisor to an exact multiple. This is meant to get rid of small errors
880 // that appears due to rounding error in our exact computation of the
881 // initial constraint given to this class.
882 //
883 // Each adjustement will cause the initial_rhs_remainder to increase, and we
884 // do not want to increase it above divisor. Our threshold below guarantees
885 // this. Note that the higher the rhs_remainder becomes, the more the
886 // function f() has a chance to reduce the violation, so it is not always a
887 // good idea to use all the slack we have between initial_rhs_remainder and
888 // divisor.
889 //
890 // TODO(user): If possible, it might be better to complement these
891 // variables. Even if the adjusted lp_values end up larger, if we loose less
892 // when taking f(), then we will have a better violation.
893 const IntegerValue adjust_threshold =
894 (divisor - initial_rhs_remainder - 1) / IntegerValue(size);
895 if (adjust_threshold > 0) {
896 // Even before we finish the adjust, we can have a lower bound on the
897 // activily loss using this divisor, and so we can abort early. This is
898 // similar to what is done below in the function.
899 bool early_abort = false;
900 double loss_lb = 0.0;
901 const double threshold = ToDouble(initial_rhs_remainder);
902
903 for (int i = 0; i < relevant_coeffs_.size(); ++i) {
904 // Compute the difference of coeff with the next multiple of divisor.
905 const IntegerValue coeff = relevant_coeffs_[i];
906 const IntegerValue remainder =
907 CeilRatio(coeff, divisor) * divisor - coeff;
908
909 if (divisor - remainder <= initial_rhs_remainder) {
910 // We do not know exactly f() yet, but it will always round to the
911 // floor of the division by divisor in this case.
912 loss_lb += ToDouble(divisor - remainder) * relevant_lp_values_[i];
913 if (loss_lb >= threshold) {
914 early_abort = true;
915 break;
916 }
917 }
918
919 // Adjust coeff of the form k * divisor - epsilon.
920 const IntegerValue diff = relevant_bound_diffs_[i];
921 if (remainder > 0 && remainder <= adjust_threshold &&
922 CapProd(diff.value(), remainder.value()) <= adjust_threshold) {
923 temp_ub += remainder * diff;
924 adjusted_coeffs_.push_back({i, coeff + remainder});
925 }
926 }
927
928 if (early_abort) continue;
929 }
930
931 // Create the super-additive function f().
932 const IntegerValue rhs_remainder =
933 temp_ub - FloorRatio(temp_ub, divisor) * divisor;
934 if (rhs_remainder == 0) continue;
935
937 rhs_remainder, divisor, GetFactorT(rhs_remainder, divisor, max_t),
938 options.max_scaling);
939
940 // As we round coefficients, we will compute the loss compared to the
941 // current scaled constraint activity. As soon as this loss crosses the
942 // slack, then we known that there is no violation and we can abort early.
943 //
944 // TODO(user): modulo the scaling, we could compute the exact threshold
945 // using our current best cut. Note that we also have to account the change
946 // in slack due to the adjust code above.
947 const double scaling = ToDouble(f(divisor)) / ToDouble(divisor);
948 const double threshold = scaling * ToDouble(rhs_remainder);
949 double loss = 0.0;
950
951 // Apply f() to the cut and compute the cut violation. Note that it is
952 // okay to just look at the relevant indices since the other have a lp
953 // value which is almost zero. Doing it like this is faster, and even if
954 // the max_magnitude might be off it should still be relevant enough.
955 double violation = -ToDouble(f(temp_ub));
956 double l2_norm = 0.0;
957 bool early_abort = false;
958 int adjusted_coeffs_index = 0;
959 for (int i = 0; i < relevant_coeffs_.size(); ++i) {
960 IntegerValue coeff = relevant_coeffs_[i];
961
962 // Adjust coeff according to our previous computation if needed.
963 if (adjusted_coeffs_index < adjusted_coeffs_.size() &&
964 adjusted_coeffs_[adjusted_coeffs_index].first == i) {
965 coeff = adjusted_coeffs_[adjusted_coeffs_index].second;
966 adjusted_coeffs_index++;
967 }
968
969 if (coeff == 0) continue;
970 const IntegerValue new_coeff = f(coeff);
971 const double new_coeff_double = ToDouble(new_coeff);
972 const double lp_value = relevant_lp_values_[i];
973
974 l2_norm += new_coeff_double * new_coeff_double;
975 violation += new_coeff_double * lp_value;
976 loss += (scaling * ToDouble(coeff) - new_coeff_double) * lp_value;
977 if (loss >= threshold) {
978 early_abort = true;
979 break;
980 }
981 }
982 if (early_abort) continue;
983
984 // Here we scale by the L2 norm over the "relevant" positions. This seems
985 // to work slighly better in practice.
986 violation /= sqrt(l2_norm);
987 if (violation > best_scaled_violation) {
988 best_scaled_violation = violation;
989 best_divisor = divisor;
990 }
991 }
992
993 if (best_divisor == 0) {
994 *cut = LinearConstraint(IntegerValue(0), IntegerValue(0));
995 return;
996 }
997
998 // Adjust coefficients.
999 //
1000 // TODO(user): It might make sense to also adjust the one with a small LP
1001 // value, but then the cut will be slighlty different than the one we computed
1002 // above. Try with and without maybe?
1003 const IntegerValue initial_rhs_remainder =
1004 cut->ub - FloorRatio(cut->ub, best_divisor) * best_divisor;
1005 const IntegerValue adjust_threshold =
1006 (best_divisor - initial_rhs_remainder - 1) / IntegerValue(size);
1007 if (adjust_threshold > 0) {
1008 for (int i = 0; i < relevant_indices_.size(); ++i) {
1009 const int index = relevant_indices_[i];
1010 const IntegerValue diff = relevant_bound_diffs_[i];
1011 if (diff > adjust_threshold) continue;
1012
1013 // Adjust coeff of the form k * best_divisor - epsilon.
1014 const IntegerValue coeff = cut->coeffs[index];
1015 const IntegerValue remainder =
1016 CeilRatio(coeff, best_divisor) * best_divisor - coeff;
1017 if (CapProd(diff.value(), remainder.value()) <= adjust_threshold) {
1018 cut->ub += remainder * diff;
1019 cut->coeffs[index] += remainder;
1020 }
1021 }
1022 }
1023
1024 // Create the super-additive function f().
1025 //
1026 // TODO(user): Try out different rounding function and keep the best. We can
1027 // change max_t and max_scaling. It might not be easy to choose which cut is
1028 // the best, but we can at least know for sure if one dominate the other
1029 // completely. That is, if for all coeff f(coeff)/f(divisor) is greater than
1030 // or equal to the same value for another function f.
1031 const IntegerValue rhs_remainder =
1032 cut->ub - FloorRatio(cut->ub, best_divisor) * best_divisor;
1033 IntegerValue factor_t = GetFactorT(rhs_remainder, best_divisor, max_t);
1034 auto f = GetSuperAdditiveRoundingFunction(rhs_remainder, best_divisor,
1035 factor_t, options.max_scaling);
1036
1037 // Look amongst all our possible function f() for one that dominate greedily
1038 // our current best one. Note that we prefer lower scaling factor since that
1039 // result in a cut with lower coefficients.
1040 remainders_.clear();
1041 for (int i = 0; i < size; ++i) {
1042 const IntegerValue coeff = cut->coeffs[i];
1043 const IntegerValue r =
1044 coeff - FloorRatio(coeff, best_divisor) * best_divisor;
1045 if (r > rhs_remainder) remainders_.push_back(r);
1046 }
1047 gtl::STLSortAndRemoveDuplicates(&remainders_);
1048 if (remainders_.size() <= 100) {
1049 best_rs_.clear();
1050 for (const IntegerValue r : remainders_) {
1051 best_rs_.push_back(f(r));
1052 }
1053 IntegerValue best_d = f(best_divisor);
1054
1055 // Note that the complexity seems high 100 * 2 * options.max_scaling, but
1056 // this only run on cuts that are already efficient and the inner loop tend
1057 // to abort quickly. I didn't see this code in the cpu profile so far.
1058 for (const IntegerValue t :
1059 {IntegerValue(1), GetFactorT(rhs_remainder, best_divisor, max_t)}) {
1060 for (IntegerValue s(2); s <= options.max_scaling; ++s) {
1061 const auto g =
1062 GetSuperAdditiveRoundingFunction(rhs_remainder, best_divisor, t, s);
1063 int num_strictly_better = 0;
1064 rs_.clear();
1065 const IntegerValue d = g(best_divisor);
1066 for (int i = 0; i < best_rs_.size(); ++i) {
1067 const IntegerValue temp = g(remainders_[i]);
1068 if (temp * best_d < best_rs_[i] * d) break;
1069 if (temp * best_d > best_rs_[i] * d) num_strictly_better++;
1070 rs_.push_back(temp);
1071 }
1072 if (rs_.size() == best_rs_.size() && num_strictly_better > 0) {
1073 f = g;
1074 factor_t = t;
1075 best_rs_ = rs_;
1076 best_d = d;
1077 }
1078 }
1079 }
1080 }
1081
1082 // Starts to apply f() to the cut. We only apply it to the rhs here, the
1083 // coefficient will be done after the potential lifting of some Booleans.
1084 cut->ub = f(cut->ub);
1085 tmp_terms_.clear();
1086
1087 // Lift some implied bounds Booleans. Note that we will add them after
1088 // "size" so they will be ignored in the second loop below.
1089 num_lifted_booleans_ = 0;
1090 if (ib_processor != nullptr) {
1091 for (int i = 0; i < size; ++i) {
1092 const IntegerValue coeff = cut->coeffs[i];
1093 if (coeff == 0) continue;
1094
1095 IntegerVariable var = cut->vars[i];
1096 if (change_sign_at_postprocessing_[i]) {
1097 var = NegationOf(var);
1098 }
1099
1101 ib_processor->GetCachedImpliedBoundInfo(var);
1102
1103 // Avoid overflow.
1104 if (CapProd(CapProd(std::abs(coeff.value()), factor_t.value()),
1105 info.bound_diff.value()) == kint64max) {
1106 continue;
1107 }
1108
1109 // Because X = bound_diff * B + S
1110 // We can replace coeff * X by the expression before applying f:
1111 // = f(coeff * bound_diff) * B + f(coeff) * [X - bound_diff * B]
1112 // = f(coeff) * X + (f(coeff * bound_diff) - f(coeff) * bound_diff] B
1113 // So we can "lift" B into the cut.
1114 const IntegerValue coeff_b =
1115 f(coeff * info.bound_diff) - f(coeff) * info.bound_diff;
1116 CHECK_GE(coeff_b, 0);
1117 if (coeff_b == 0) continue;
1118
1119 ++num_lifted_booleans_;
1120 if (info.is_positive) {
1121 tmp_terms_.push_back({info.bool_var, coeff_b});
1122 } else {
1123 tmp_terms_.push_back({info.bool_var, -coeff_b});
1124 cut->ub = CapAdd(-coeff_b.value(), cut->ub.value());
1125 }
1126 }
1127 }
1128
1129 // Apply f() to the cut.
1130 //
1131 // Remove the bound shifts so the constraint is expressed in the original
1132 // variables.
1133 for (int i = 0; i < size; ++i) {
1134 IntegerValue coeff = cut->coeffs[i];
1135 if (coeff == 0) continue;
1136 coeff = f(coeff);
1137 if (coeff == 0) continue;
1138 if (change_sign_at_postprocessing_[i]) {
1139 cut->ub = IntegerValue(
1140 CapAdd((coeff * -upper_bounds[i]).value(), cut->ub.value()));
1141 tmp_terms_.push_back({cut->vars[i], -coeff});
1142 } else {
1143 cut->ub = IntegerValue(
1144 CapAdd((coeff * lower_bounds[i]).value(), cut->ub.value()));
1145 tmp_terms_.push_back({cut->vars[i], coeff});
1146 }
1147 }
1148
1149 // Basic post-processing.
1150 CleanTermsAndFillConstraint(&tmp_terms_, cut);
1151 RemoveZeroTerms(cut);
1152 DivideByGCD(cut);
1153}
1154
1156 const LinearConstraint base_ct, const std::vector<double>& lp_values,
1157 const std::vector<IntegerValue>& lower_bounds,
1158 const std::vector<IntegerValue>& upper_bounds) {
1159 const int base_size = lp_values.size();
1160
1161 // Fill terms with a rewrite of the base constraint where all coeffs &
1162 // variables are positive by using either (X - LB) or (UB - X) as new
1163 // variables.
1164 terms_.clear();
1165 IntegerValue rhs = base_ct.ub;
1166 IntegerValue sum_of_diff(0);
1167 IntegerValue max_base_magnitude(0);
1168 for (int i = 0; i < base_size; ++i) {
1169 const IntegerValue coeff = base_ct.coeffs[i];
1170 const IntegerValue positive_coeff = IntTypeAbs(coeff);
1171 max_base_magnitude = std::max(max_base_magnitude, positive_coeff);
1172 const IntegerValue bound_diff = upper_bounds[i] - lower_bounds[i];
1173 if (!AddProductTo(positive_coeff, bound_diff, &sum_of_diff)) {
1174 return false;
1175 }
1176 const IntegerValue diff = positive_coeff * bound_diff;
1177 if (coeff > 0) {
1178 if (!AddProductTo(-coeff, lower_bounds[i], &rhs)) return false;
1179 terms_.push_back(
1180 {i, ToDouble(upper_bounds[i]) - lp_values[i], positive_coeff, diff});
1181 } else {
1182 if (!AddProductTo(-coeff, upper_bounds[i], &rhs)) return false;
1183 terms_.push_back(
1184 {i, lp_values[i] - ToDouble(lower_bounds[i]), positive_coeff, diff});
1185 }
1186 }
1187
1188 // Try a simple cover heuristic.
1189 // Look for violated CUT of the form: sum (UB - X) or (X - LB) >= 1.
1190 double activity = 0.0;
1191 int new_size = 0;
1192 std::sort(terms_.begin(), terms_.end(), [](const Term& a, const Term& b) {
1193 if (a.dist_to_max_value == b.dist_to_max_value) {
1194 // Prefer low coefficients if the distance is the same.
1195 return a.positive_coeff < b.positive_coeff;
1196 }
1197 return a.dist_to_max_value < b.dist_to_max_value;
1198 });
1199 for (int i = 0; i < terms_.size(); ++i) {
1200 const Term& term = terms_[i];
1201 activity += term.dist_to_max_value;
1202
1203 // As an heuristic we select all the term so that the sum of distance
1204 // to the upper bound is <= 1.0. If the corresponding rhs is negative, then
1205 // we will have a cut of violation at least 0.0. Note that this violation
1206 // can be improved by the lifting.
1207 //
1208 // TODO(user): experiment with different threshold (even greater than one).
1209 // Or come up with an algo that incorporate the lifting into the heuristic.
1210 if (activity > 1.0) {
1211 new_size = i; // before this entry.
1212 break;
1213 }
1214
1215 rhs -= term.diff;
1216 }
1217
1218 // If the rhs is now negative, we have a cut.
1219 //
1220 // Note(user): past this point, now that a given "base" cover has been chosen,
1221 // we basically compute the cut (of the form sum X <= bound) with the maximum
1222 // possible violation. Note also that we lift as much as possible, so we don't
1223 // necessarilly optimize for the cut efficacity though. But we do get a
1224 // stronger cut.
1225 if (rhs >= 0) return false;
1226 if (new_size == 0) return false;
1227
1228 // Transform to a minimal cover. We want to greedily remove the largest coeff
1229 // first, so we have more chance for the "lifting" below which can increase
1230 // the cut violation. If the coeff are the same, we prefer to remove high
1231 // distance from upper bound first.
1232 //
1233 // We compute the cut at the same time.
1234 terms_.resize(new_size);
1235 std::sort(terms_.begin(), terms_.end(), [](const Term& a, const Term& b) {
1236 if (a.positive_coeff == b.positive_coeff) {
1237 return a.dist_to_max_value > b.dist_to_max_value;
1238 }
1239 return a.positive_coeff > b.positive_coeff;
1240 });
1241 in_cut_.assign(base_ct.vars.size(), false);
1242 cut_.ClearTerms();
1243 cut_.lb = kMinIntegerValue;
1244 cut_.ub = IntegerValue(-1);
1245 IntegerValue max_coeff(0);
1246 for (const Term term : terms_) {
1247 if (term.diff + rhs < 0) {
1248 rhs += term.diff;
1249 continue;
1250 }
1251 in_cut_[term.index] = true;
1252 max_coeff = std::max(max_coeff, term.positive_coeff);
1253 cut_.vars.push_back(base_ct.vars[term.index]);
1254 if (base_ct.coeffs[term.index] > 0) {
1255 cut_.coeffs.push_back(IntegerValue(1));
1256 cut_.ub += upper_bounds[term.index];
1257 } else {
1258 cut_.coeffs.push_back(IntegerValue(-1));
1259 cut_.ub -= lower_bounds[term.index];
1260 }
1261 }
1262
1263 // In case the max_coeff variable is not binary, it might be possible to
1264 // tighten the cut a bit more.
1265 //
1266 // Note(user): I never observed this on the miplib so far.
1267 if (max_coeff == 0) return true;
1268 if (max_coeff < -rhs) {
1269 const IntegerValue m = FloorRatio(-rhs - 1, max_coeff);
1270 rhs += max_coeff * m;
1271 cut_.ub -= m;
1272 }
1273 CHECK_LT(rhs, 0);
1274
1275 // Lift all at once the variables not used in the cover.
1276 //
1277 // We have a cut of the form sum_i X_i <= b that we will lift into
1278 // sum_i scaling X_i + sum f(base_coeff_j) X_j <= b * scaling.
1279 //
1280 // Using the super additivity of f() and how we construct it,
1281 // we know that: sum_j base_coeff_j X_j <= N * max_coeff + (max_coeff - slack)
1282 // implies that: sum_j f(base_coeff_j) X_j <= N * scaling.
1283 //
1284 // 1/ cut > b -(N+1) => original sum + (N+1) * max_coeff >= rhs + slack
1285 // 2/ rewrite 1/ as : scaling * cut >= scaling * b - scaling * N => ...
1286 // 3/ lift > N * scaling => lift_sum > N * max_coeff + (max_coeff - slack)
1287 // And adding 2/ + 3/ we prove what we want:
1288 // cut * scaling + lift > b * scaling => original_sum + lift_sum > rhs.
1289 const IntegerValue slack = -rhs;
1290 const IntegerValue remainder = max_coeff - slack;
1291 max_base_magnitude = std::max(max_base_magnitude, IntTypeAbs(cut_.ub));
1292 const IntegerValue max_scaling(std::min(
1293 IntegerValue(60), FloorRatio(kMaxIntegerValue, max_base_magnitude)));
1294 const auto f = GetSuperAdditiveRoundingFunction(remainder, max_coeff,
1295 IntegerValue(1), max_scaling);
1296
1297 const IntegerValue scaling = f(max_coeff);
1298 if (scaling > 1) {
1299 for (int i = 0; i < cut_.coeffs.size(); ++i) cut_.coeffs[i] *= scaling;
1300 cut_.ub *= scaling;
1301 }
1302
1303 num_lifting_ = 0;
1304 for (int i = 0; i < base_size; ++i) {
1305 if (in_cut_[i]) continue;
1306 const IntegerValue positive_coeff = IntTypeAbs(base_ct.coeffs[i]);
1307 const IntegerValue new_coeff = f(positive_coeff);
1308 if (new_coeff == 0) continue;
1309
1310 ++num_lifting_;
1311 if (base_ct.coeffs[i] > 0) {
1312 // Add new_coeff * (X - LB)
1313 cut_.coeffs.push_back(new_coeff);
1314 cut_.vars.push_back(base_ct.vars[i]);
1315 cut_.ub += lower_bounds[i] * new_coeff;
1316 } else {
1317 // Add new_coeff * (UB - X)
1318 cut_.coeffs.push_back(-new_coeff);
1319 cut_.vars.push_back(base_ct.vars[i]);
1320 cut_.ub -= upper_bounds[i] * new_coeff;
1321 }
1322 }
1323
1324 if (scaling > 1) DivideByGCD(&cut_);
1325 return true;
1326}
1327
1329 IntegerVariable x,
1330 IntegerVariable y,
1331 Model* model) {
1332 CutGenerator result;
1333 result.vars = {z, x, y};
1334
1335 IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>();
1336 result.generate_cuts =
1337 [z, x, y, integer_trail](
1339 LinearConstraintManager* manager) {
1340 const int64 x_lb = integer_trail->LevelZeroLowerBound(x).value();
1341 const int64 x_ub = integer_trail->LevelZeroUpperBound(x).value();
1342 const int64 y_lb = integer_trail->LevelZeroLowerBound(y).value();
1343 const int64 y_ub = integer_trail->LevelZeroUpperBound(y).value();
1344
1345 // TODO(user): Compute a better bound (int_max / 4 ?).
1346 const int64 kMaxSafeInteger = (int64{1} << 53) - 1;
1347
1348 if (CapProd(x_ub, y_ub) >= kMaxSafeInteger) {
1349 VLOG(3) << "Potential overflow in PositiveMultiplicationCutGenerator";
1350 return;
1351 }
1352
1353 const double x_lp_value = lp_values[x];
1354 const double y_lp_value = lp_values[y];
1355 const double z_lp_value = lp_values[z];
1356
1357 // TODO(user): As the bounds change monotonically, these cuts
1358 // dominate any previous one. try to keep a reference to the cut and
1359 // replace it. Alternatively, add an API for a level-zero bound change
1360 // callback.
1361
1362 // Cut -z + x_coeff * x + y_coeff* y <= rhs
1363 auto try_add_above_cut = [manager, z_lp_value, x_lp_value, y_lp_value,
1364 x, y, z, &lp_values](
1365 int64 x_coeff, int64 y_coeff, int64 rhs) {
1366 if (-z_lp_value + x_lp_value * x_coeff + y_lp_value * y_coeff >=
1367 rhs + kMinCutViolation) {
1368 LinearConstraint cut;
1369 cut.vars.push_back(z);
1370 cut.coeffs.push_back(IntegerValue(-1));
1371 if (x_coeff != 0) {
1372 cut.vars.push_back(x);
1373 cut.coeffs.push_back(IntegerValue(x_coeff));
1374 }
1375 if (y_coeff != 0) {
1376 cut.vars.push_back(y);
1377 cut.coeffs.push_back(IntegerValue(y_coeff));
1378 }
1379 cut.lb = kMinIntegerValue;
1380 cut.ub = IntegerValue(rhs);
1381 manager->AddCut(cut, "PositiveProduct", lp_values);
1382 }
1383 };
1384
1385 // Cut -z + x_coeff * x + y_coeff* y >= rhs
1386 auto try_add_below_cut = [manager, z_lp_value, x_lp_value, y_lp_value,
1387 x, y, z, &lp_values](
1388 int64 x_coeff, int64 y_coeff, int64 rhs) {
1389 if (-z_lp_value + x_lp_value * x_coeff + y_lp_value * y_coeff <=
1390 rhs - kMinCutViolation) {
1391 LinearConstraint cut;
1392 cut.vars.push_back(z);
1393 cut.coeffs.push_back(IntegerValue(-1));
1394 if (x_coeff != 0) {
1395 cut.vars.push_back(x);
1396 cut.coeffs.push_back(IntegerValue(x_coeff));
1397 }
1398 if (y_coeff != 0) {
1399 cut.vars.push_back(y);
1400 cut.coeffs.push_back(IntegerValue(y_coeff));
1401 }
1402 cut.lb = IntegerValue(rhs);
1403 cut.ub = kMaxIntegerValue;
1404 manager->AddCut(cut, "PositiveProduct", lp_values);
1405 }
1406 };
1407
1408 // McCormick relaxation of bilinear constraints. These 4 cuts are the
1409 // exact facets of the x * y polyhedron for a bounded x and y.
1410 //
1411 // Each cut correspond to plane that contains two of the line
1412 // (x=x_lb), (x=x_ub), (y=y_lb), (y=y_ub). The easiest to
1413 // understand them is to draw the x*y curves and see the 4
1414 // planes that correspond to the convex hull of the graph.
1415 try_add_above_cut(y_lb, x_lb, x_lb * y_lb);
1416 try_add_above_cut(y_ub, x_ub, x_ub * y_ub);
1417 try_add_below_cut(y_ub, x_lb, x_lb * y_ub);
1418 try_add_below_cut(y_lb, x_ub, x_ub * y_lb);
1419 };
1420
1421 return result;
1422}
1423
1424CutGenerator CreateSquareCutGenerator(IntegerVariable y, IntegerVariable x,
1425 Model* model) {
1426 CutGenerator result;
1427 result.vars = {y, x};
1428
1429 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
1430 result.generate_cuts =
1431 [y, x, integer_trail](
1433 LinearConstraintManager* manager) {
1434 const int64 x_ub = integer_trail->LevelZeroUpperBound(x).value();
1435 const int64 x_lb = integer_trail->LevelZeroLowerBound(x).value();
1436
1437 if (x_lb == x_ub) return;
1438
1439 // Check for potential overflows.
1440 if (x_ub > (int64{1} << 31)) return;
1441 DCHECK_GE(x_lb, 0);
1442
1443 const double y_lp_value = lp_values[y];
1444 const double x_lp_value = lp_values[x];
1445
1446 // First cut: target should be below the line:
1447 // (x_lb, x_lb ^ 2) to (x_ub, x_ub ^ 2).
1448 // The slope of that line is (ub^2 - lb^2) / (ub - lb) = ub + lb.
1449 const int64 y_lb = x_lb * x_lb;
1450 const int64 above_slope = x_ub + x_lb;
1451 const double max_lp_y = y_lb + above_slope * (x_lp_value - x_lb);
1452 if (y_lp_value >= max_lp_y + kMinCutViolation) {
1453 // cut: y <= (x_lb + x_ub) * x - x_lb * x_ub
1454 LinearConstraint above_cut;
1455 above_cut.vars.push_back(y);
1456 above_cut.coeffs.push_back(IntegerValue(1));
1457 above_cut.vars.push_back(x);
1458 above_cut.coeffs.push_back(IntegerValue(-above_slope));
1459 above_cut.lb = kMinIntegerValue;
1460 above_cut.ub = IntegerValue(-x_lb * x_ub);
1461 manager->AddCut(above_cut, "SquareUpper", lp_values);
1462 }
1463
1464 // Second cut: target should be above all the lines
1465 // (value, value ^ 2) to (value + 1, (value + 1) ^ 2)
1466 // The slope of that line is 2 * value + 1
1467 //
1468 // Note that we only add one of these cuts. The one for x_lp_value in
1469 // [value, value + 1].
1470 const int64 x_floor = static_cast<int64>(std::floor(x_lp_value));
1471 const int64 below_slope = 2 * x_floor + 1;
1472 const double min_lp_y =
1473 below_slope * x_lp_value - x_floor - x_floor * x_floor;
1474 if (min_lp_y >= y_lp_value + kMinCutViolation) {
1475 // cut: y >= below_slope * (x - x_floor) + x_floor ^ 2
1476 // : y >= below_slope * x - x_floor ^ 2 - x_floor
1477 LinearConstraint below_cut;
1478 below_cut.vars.push_back(y);
1479 below_cut.coeffs.push_back(IntegerValue(1));
1480 below_cut.vars.push_back(x);
1481 below_cut.coeffs.push_back(-IntegerValue(below_slope));
1482 below_cut.lb = IntegerValue(-x_floor - x_floor * x_floor);
1483 below_cut.ub = kMaxIntegerValue;
1484 manager->AddCut(below_cut, "SquareLower", lp_values);
1485 }
1486 };
1487
1488 return result;
1489}
1490
1491void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint(
1493 LinearConstraint* cut) {
1494 ProcessUpperBoundedConstraintWithSlackCreation(
1495 /*substitute_only_inner_variables=*/false, IntegerVariable(0), lp_values,
1496 cut, nullptr);
1497}
1498
1500ImpliedBoundsProcessor::GetCachedImpliedBoundInfo(IntegerVariable var) {
1501 auto it = cache_.find(var);
1502 if (it != cache_.end()) return it->second;
1503 return BestImpliedBoundInfo();
1504}
1505
1507ImpliedBoundsProcessor::ComputeBestImpliedBound(
1508 IntegerVariable var,
1510 auto it = cache_.find(var);
1511 if (it != cache_.end()) return it->second;
1512 BestImpliedBoundInfo result;
1513 const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var);
1514 for (const ImpliedBoundEntry& entry :
1515 implied_bounds_->GetImpliedBounds(var)) {
1516 // Only process entries with a Boolean variable currently part of the LP
1517 // we are considering for this cut.
1518 //
1519 // TODO(user): the more we use cuts, the less it make sense to have a
1520 // lot of small independent LPs.
1521 if (!lp_vars_.contains(PositiveVariable(entry.literal_view))) {
1522 continue;
1523 }
1524
1525 // The equation is X = lb + diff * Bool + Slack where Bool is in [0, 1]
1526 // and slack in [0, ub - lb].
1527 const IntegerValue diff = entry.lower_bound - lb;
1528 CHECK_GE(diff, 0);
1529 const double bool_lp_value = entry.is_positive
1530 ? lp_values[entry.literal_view]
1531 : 1.0 - lp_values[entry.literal_view];
1532 const double slack_lp_value =
1533 lp_values[var] - ToDouble(lb) - bool_lp_value * ToDouble(diff);
1534
1535 // If the implied bound equation is not respected, we just add it
1536 // to implied_bound_cuts, and skip the entry for now.
1537 if (slack_lp_value < -1e-4) {
1538 LinearConstraint ib_cut;
1539 ib_cut.lb = kMinIntegerValue;
1540 std::vector<std::pair<IntegerVariable, IntegerValue>> terms;
1541 if (entry.is_positive) {
1542 // X >= Indicator * (bound - lb) + lb
1543 terms.push_back({entry.literal_view, diff});
1544 terms.push_back({var, IntegerValue(-1)});
1545 ib_cut.ub = -lb;
1546 } else {
1547 // X >= -Indicator * (bound - lb) + bound
1548 terms.push_back({entry.literal_view, -diff});
1549 terms.push_back({var, IntegerValue(-1)});
1550 ib_cut.ub = -entry.lower_bound;
1551 }
1552 CleanTermsAndFillConstraint(&terms, &ib_cut);
1553 ib_cut_pool_.AddCut(std::move(ib_cut), "IB", lp_values);
1554 continue;
1555 }
1556
1557 // We look for tight implied bounds, and amongst the tightest one, we
1558 // prefer larger coefficient in front of the Boolean.
1559 if (slack_lp_value + 1e-4 < result.slack_lp_value ||
1560 (slack_lp_value < result.slack_lp_value + 1e-4 &&
1561 diff > result.bound_diff)) {
1562 result.bool_lp_value = bool_lp_value;
1563 result.slack_lp_value = slack_lp_value;
1564
1565 result.bound_diff = diff;
1566 result.is_positive = entry.is_positive;
1567 result.bool_var = entry.literal_view;
1568 }
1569 }
1570 cache_[var] = result;
1571 return result;
1572}
1573
1574// TODO(user): restrict to a subset of the variables to not spend too much time.
1575void ImpliedBoundsProcessor::SeparateSomeImpliedBoundCuts(
1577 for (const IntegerVariable var :
1578 implied_bounds_->VariablesWithImpliedBounds()) {
1579 if (!lp_vars_.contains(PositiveVariable(var))) continue;
1580 ComputeBestImpliedBound(var, lp_values);
1581 }
1582}
1583
1584void ImpliedBoundsProcessor::ProcessUpperBoundedConstraintWithSlackCreation(
1585 bool substitute_only_inner_variables, IntegerVariable first_slack,
1587 LinearConstraint* cut, std::vector<SlackInfo>* slack_infos) {
1588 tmp_terms_.clear();
1589 IntegerValue new_ub = cut->ub;
1590 bool changed = false;
1591
1592 // TODO(user): we could relax a bit this test.
1593 int64 overflow_detection = 0;
1594
1595 const int size = cut->vars.size();
1596 for (int i = 0; i < size; ++i) {
1597 IntegerVariable var = cut->vars[i];
1598 IntegerValue coeff = cut->coeffs[i];
1599
1600 // Starts by positive coefficient.
1601 // TODO(user): Not clear this is best.
1602 if (coeff < 0) {
1603 coeff = -coeff;
1604 var = NegationOf(var);
1605 }
1606
1607 // Find the best implied bound to use.
1608 // TODO(user): We could also use implied upper bound, that is try with
1609 // NegationOf(var).
1610 const BestImpliedBoundInfo info = ComputeBestImpliedBound(var, lp_values);
1611 {
1612 // This make sure the implied bound for NegationOf(var) is "cached" so
1613 // that GetCachedImpliedBoundInfo() will work. It will also add any
1614 // relevant implied bound cut.
1615 //
1616 // TODO(user): this is a bit hacky. Find a cleaner way.
1617 ComputeBestImpliedBound(NegationOf(var), lp_values);
1618 }
1619
1620 const int old_size = tmp_terms_.size();
1621
1622 // Shall we keep the original term ?
1623 bool keep_term = false;
1624 if (info.bool_var == kNoIntegerVariable) keep_term = true;
1625 if (CapProd(std::abs(coeff.value()), info.bound_diff.value()) ==
1626 kint64max) {
1627 keep_term = true;
1628 }
1629
1630 // TODO(user): On some problem, not replacing the variable at their bound
1631 // by an implied bounds seems beneficial. This is especially the case on
1632 // g200x740.mps.gz
1633 //
1634 // Note that in ComputeCut() the variable with an LP value at the bound do
1635 // not contribute to the cut efficacity (no loss) but do contribute to the
1636 // various heuristic based on the coefficient magnitude.
1637 if (substitute_only_inner_variables) {
1638 const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var);
1639 const IntegerValue ub = integer_trail_->LevelZeroUpperBound(var);
1640 if (lp_values[var] - ToDouble(lb) < 1e-2) keep_term = true;
1641 if (ToDouble(ub) - lp_values[var] < 1e-2) keep_term = true;
1642 }
1643
1644 // This is when we do not add slack.
1645 if (slack_infos == nullptr) {
1646 // We do not want to loose anything, so we only replace if the slack lp is
1647 // zero.
1648 if (info.slack_lp_value > 1e-6) keep_term = true;
1649 }
1650
1651 if (keep_term) {
1652 tmp_terms_.push_back({var, coeff});
1653 } else {
1654 // Substitute.
1655 const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var);
1656 const IntegerValue ub = integer_trail_->LevelZeroUpperBound(var);
1657
1658 SlackInfo slack_info;
1659 slack_info.lp_value = info.slack_lp_value;
1660 slack_info.lb = 0;
1661 slack_info.ub = ub - lb;
1662
1663 if (info.is_positive) {
1664 // X = Indicator * diff + lb + Slack
1665 tmp_terms_.push_back({info.bool_var, coeff * info.bound_diff});
1666 if (!AddProductTo(-coeff, lb, &new_ub)) {
1667 VLOG(2) << "Overflow";
1668 return;
1669 }
1670 if (slack_infos != nullptr) {
1671 tmp_terms_.push_back({first_slack, coeff});
1672 first_slack += 2;
1673
1674 // slack = X - Indicator * info.bound_diff - lb;
1675 slack_info.terms.push_back({var, IntegerValue(1)});
1676 slack_info.terms.push_back({info.bool_var, -info.bound_diff});
1677 slack_info.offset = -lb;
1678 slack_infos->push_back(slack_info);
1679 }
1680 } else {
1681 // X = (1 - Indicator) * (diff) + lb + Slack
1682 // X = -Indicator * (diff) + lb + diff + Slack
1683 tmp_terms_.push_back({info.bool_var, -coeff * info.bound_diff});
1684 if (!AddProductTo(-coeff, lb + info.bound_diff, &new_ub)) {
1685 VLOG(2) << "Overflow";
1686 return;
1687 }
1688 if (slack_infos != nullptr) {
1689 tmp_terms_.push_back({first_slack, coeff});
1690 first_slack += 2;
1691
1692 // slack = X + Indicator * info.bound_diff - lb - diff;
1693 slack_info.terms.push_back({var, IntegerValue(1)});
1694 slack_info.terms.push_back({info.bool_var, +info.bound_diff});
1695 slack_info.offset = -lb - info.bound_diff;
1696 slack_infos->push_back(slack_info);
1697 }
1698 }
1699 changed = true;
1700 }
1701
1702 // Add all the new terms coefficient to the overflow detection to avoid
1703 // issue when merging terms refering to the same variable.
1704 for (int i = old_size; i < tmp_terms_.size(); ++i) {
1705 overflow_detection =
1706 CapAdd(overflow_detection, std::abs(tmp_terms_[i].second.value()));
1707 }
1708 }
1709
1710 if (overflow_detection >= kMaxIntegerValue) {
1711 VLOG(2) << "Overflow";
1712 return;
1713 }
1714 if (!changed) return;
1715
1716 // Update the cut.
1717 //
1718 // Note that because of our overflow_detection variable, there should be
1719 // no integer overflow when we merge identical terms.
1720 cut->lb = kMinIntegerValue; // Not relevant.
1721 cut->ub = new_ub;
1722 CleanTermsAndFillConstraint(&tmp_terms_, cut);
1723}
1724
1725bool ImpliedBoundsProcessor::DebugSlack(IntegerVariable first_slack,
1726 const LinearConstraint& initial_cut,
1727 const LinearConstraint& cut,
1728 const std::vector<SlackInfo>& info) {
1729 tmp_terms_.clear();
1730 IntegerValue new_ub = cut.ub;
1731 for (int i = 0; i < cut.vars.size(); ++i) {
1732 // Simple copy for non-slack variables.
1733 if (cut.vars[i] < first_slack) {
1734 tmp_terms_.push_back({cut.vars[i], cut.coeffs[i]});
1735 continue;
1736 }
1737
1738 // Replace slack by its definition.
1739 const IntegerValue multiplier = cut.coeffs[i];
1740 const int index = (cut.vars[i].value() - first_slack.value()) / 2;
1741 for (const std::pair<IntegerVariable, IntegerValue>& term :
1742 info[index].terms) {
1743 tmp_terms_.push_back({term.first, term.second * multiplier});
1744 }
1745 new_ub -= multiplier * info[index].offset;
1746 }
1747
1748 LinearConstraint tmp_cut;
1749 tmp_cut.lb = kMinIntegerValue; // Not relevant.
1750 tmp_cut.ub = new_ub;
1751 CleanTermsAndFillConstraint(&tmp_terms_, &tmp_cut);
1752 MakeAllVariablesPositive(&tmp_cut);
1753
1754 // We need to canonicalize the initial_cut too for comparison. Note that we
1755 // only use this for debug, so we don't care too much about the memory and
1756 // extra time.
1757 // TODO(user): Expose CanonicalizeConstraint() from the manager.
1758 LinearConstraint tmp_copy;
1759 tmp_terms_.clear();
1760 for (int i = 0; i < initial_cut.vars.size(); ++i) {
1761 tmp_terms_.push_back({initial_cut.vars[i], initial_cut.coeffs[i]});
1762 }
1763 tmp_copy.lb = kMinIntegerValue; // Not relevant.
1764 tmp_copy.ub = new_ub;
1765 CleanTermsAndFillConstraint(&tmp_terms_, &tmp_copy);
1766 MakeAllVariablesPositive(&tmp_copy);
1767
1768 if (tmp_cut == tmp_copy) return true;
1769
1770 LOG(INFO) << first_slack;
1771 LOG(INFO) << tmp_copy.DebugString();
1772 LOG(INFO) << cut.DebugString();
1773 LOG(INFO) << tmp_cut.DebugString();
1774 return false;
1775}
1776
1777namespace {
1778
1779void TryToGenerateAllDiffCut(
1780 const std::vector<std::pair<double, IntegerVariable>>& sorted_vars_lp,
1781 const IntegerTrail& integer_trail,
1783 LinearConstraintManager* manager) {
1784 Domain current_union;
1785 std::vector<IntegerVariable> current_set_vars;
1786 double sum = 0.0;
1787 for (auto value_var : sorted_vars_lp) {
1788 sum += value_var.first;
1789 const IntegerVariable var = value_var.second;
1790 // TODO(user): The union of the domain of the variable being considered
1791 // does not give the tightest bounds, try to get better bounds.
1792 current_union =
1793 current_union.UnionWith(integer_trail.InitialVariableDomain(var));
1794 current_set_vars.push_back(var);
1795 const int64 required_min_sum =
1796 SumOfKMinValueInDomain(current_union, current_set_vars.size());
1797 const int64 required_max_sum =
1798 SumOfKMaxValueInDomain(current_union, current_set_vars.size());
1799 if (sum < required_min_sum || sum > required_max_sum) {
1800 LinearConstraint cut;
1801 for (IntegerVariable var : current_set_vars) {
1802 cut.AddTerm(var, IntegerValue(1));
1803 }
1804 cut.lb = IntegerValue(required_min_sum);
1805 cut.ub = IntegerValue(required_max_sum);
1806 manager->AddCut(cut, "all_diff", lp_values);
1807 // NOTE: We can extend the current set but it is more helpful to generate
1808 // the cut on a different set of variables so we reset the counters.
1809 sum = 0.0;
1810 current_set_vars.clear();
1811 current_union = Domain();
1812 }
1813 }
1814}
1815
1816} // namespace
1817
1819 const std::vector<IntegerVariable>& vars, Model* model) {
1820 CutGenerator result;
1821 result.vars = vars;
1822 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
1823 Trail* trail = model->GetOrCreate<Trail>();
1824 result.generate_cuts =
1825 [vars, integer_trail, trail](
1827 LinearConstraintManager* manager) {
1828 // These cuts work at all levels but the generator adds too many cuts on
1829 // some instances and degrade the performance so we only use it at level
1830 // 0.
1831 if (trail->CurrentDecisionLevel() > 0) return;
1832 std::vector<std::pair<double, IntegerVariable>> sorted_vars;
1833 for (const IntegerVariable var : vars) {
1834 if (integer_trail->LevelZeroLowerBound(var) ==
1835 integer_trail->LevelZeroUpperBound(var)) {
1836 continue;
1837 }
1838 sorted_vars.push_back(std::make_pair(lp_values[var], var));
1839 }
1840 std::sort(sorted_vars.begin(), sorted_vars.end());
1841 TryToGenerateAllDiffCut(sorted_vars, *integer_trail, lp_values,
1842 manager);
1843 // Other direction.
1844 std::reverse(sorted_vars.begin(), sorted_vars.end());
1845 TryToGenerateAllDiffCut(sorted_vars, *integer_trail, lp_values,
1846 manager);
1847 };
1848 VLOG(1) << "Created all_diff cut generator of size: " << vars.size();
1849 return result;
1850}
1851
1852namespace {
1853// Returns max((w2i - w1i)*Li, (w2i - w1i)*Ui).
1854IntegerValue MaxCornerDifference(const IntegerVariable var,
1855 const IntegerValue w1_i,
1856 const IntegerValue w2_i,
1857 const IntegerTrail& integer_trail) {
1858 const IntegerValue lb = integer_trail.LevelZeroLowerBound(var);
1859 const IntegerValue ub = integer_trail.LevelZeroUpperBound(var);
1860 return std::max((w2_i - w1_i) * lb, (w2_i - w1_i) * ub);
1861}
1862
1863// This is the coefficient of zk in the cut, where k = max_index.
1864// MPlusCoefficient_ki = max((wki - wI(i)i) * Li,
1865// (wki - wI(i)i) * Ui)
1866// = max corner difference for variable i,
1867// target expr I(i), max expr k.
1868// The coefficient of zk is Sum(i=1..n)(MPlusCoefficient_ki) + bk
1869IntegerValue MPlusCoefficient(
1870 const std::vector<IntegerVariable>& x_vars,
1871 const std::vector<LinearExpression>& exprs,
1872 const absl::StrongVector<IntegerVariable, int>& variable_partition,
1873 const int max_index, const IntegerTrail& integer_trail) {
1874 IntegerValue coeff = exprs[max_index].offset;
1875 // TODO(user): This algo is quadratic since GetCoefficientOfPositiveVar()
1876 // is linear. This can be optimized (better complexity) if needed.
1877 for (const IntegerVariable var : x_vars) {
1878 const int target_index = variable_partition[var];
1879 if (max_index != target_index) {
1880 coeff += MaxCornerDifference(
1881 var, GetCoefficientOfPositiveVar(var, exprs[target_index]),
1882 GetCoefficientOfPositiveVar(var, exprs[max_index]), integer_trail);
1883 }
1884 }
1885 return coeff;
1886}
1887
1888// Compute the value of
1889// rhs = wI(i)i * xi + Sum(k=1..d)(MPlusCoefficient_ki * zk)
1890// for variable xi for given target index I(i).
1891double ComputeContribution(
1892 const IntegerVariable xi_var, const std::vector<IntegerVariable>& z_vars,
1893 const std::vector<LinearExpression>& exprs,
1895 const IntegerTrail& integer_trail, const int target_index) {
1896 CHECK_GE(target_index, 0);
1897 CHECK_LT(target_index, exprs.size());
1898 const LinearExpression& target_expr = exprs[target_index];
1899 const double xi_value = lp_values[xi_var];
1900 const IntegerValue wt_i = GetCoefficientOfPositiveVar(xi_var, target_expr);
1901 double contrib = wt_i.value() * xi_value;
1902 for (int expr_index = 0; expr_index < exprs.size(); ++expr_index) {
1903 if (expr_index == target_index) continue;
1904 const LinearExpression& max_expr = exprs[expr_index];
1905 const double z_max_value = lp_values[z_vars[expr_index]];
1906 const IntegerValue corner_value = MaxCornerDifference(
1907 xi_var, wt_i, GetCoefficientOfPositiveVar(xi_var, max_expr),
1908 integer_trail);
1909 contrib += corner_value.value() * z_max_value;
1910 }
1911 return contrib;
1912}
1913} // namespace
1914
1916 const IntegerVariable target, const std::vector<LinearExpression>& exprs,
1917 const std::vector<IntegerVariable>& z_vars, Model* model) {
1918 CutGenerator result;
1919 std::vector<IntegerVariable> x_vars;
1920 result.vars = {target};
1921 const int num_exprs = exprs.size();
1922 for (int i = 0; i < num_exprs; ++i) {
1923 result.vars.push_back(z_vars[i]);
1924 x_vars.insert(x_vars.end(), exprs[i].vars.begin(), exprs[i].vars.end());
1925 }
1927 // All expressions should only contain positive variables.
1928 DCHECK(std::all_of(x_vars.begin(), x_vars.end(), [](IntegerVariable var) {
1929 return VariableIsPositive(var);
1930 }));
1931 result.vars.insert(result.vars.end(), x_vars.begin(), x_vars.end());
1932
1933 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
1934 result.generate_cuts =
1935 [x_vars, z_vars, target, num_exprs, exprs, integer_trail, model](
1937 LinearConstraintManager* manager) {
1939 lp_values.size(), -1);
1940 absl::StrongVector<IntegerVariable, double> variable_partition_contrib(
1941 lp_values.size(), std::numeric_limits<double>::infinity());
1942 for (int expr_index = 0; expr_index < num_exprs; ++expr_index) {
1943 for (const IntegerVariable var : x_vars) {
1944 const double contribution = ComputeContribution(
1945 var, z_vars, exprs, lp_values, *integer_trail, expr_index);
1946 const double prev_contribution = variable_partition_contrib[var];
1947 if (contribution < prev_contribution) {
1948 variable_partition[var] = expr_index;
1949 variable_partition_contrib[var] = contribution;
1950 }
1951 }
1952 }
1953
1954 LinearConstraintBuilder cut(model, /*lb=*/IntegerValue(0),
1955 /*ub=*/kMaxIntegerValue);
1956 double violation = lp_values[target];
1957 cut.AddTerm(target, IntegerValue(-1));
1958
1959 for (const IntegerVariable xi_var : x_vars) {
1960 const int input_index = variable_partition[xi_var];
1961 const LinearExpression& expr = exprs[input_index];
1962 const IntegerValue coeff = GetCoefficientOfPositiveVar(xi_var, expr);
1963 if (coeff != IntegerValue(0)) {
1964 cut.AddTerm(xi_var, coeff);
1965 }
1966 violation -= coeff.value() * lp_values[xi_var];
1967 }
1968 for (int expr_index = 0; expr_index < num_exprs; ++expr_index) {
1969 const IntegerVariable z_var = z_vars[expr_index];
1970 const IntegerValue z_coeff = MPlusCoefficient(
1971 x_vars, exprs, variable_partition, expr_index, *integer_trail);
1972 if (z_coeff != IntegerValue(0)) {
1973 cut.AddTerm(z_var, z_coeff);
1974 }
1975 violation -= z_coeff.value() * lp_values[z_var];
1976 }
1977 if (violation > 1e-2) {
1978 manager->AddCut(cut.Build(), "LinMax", lp_values);
1979 }
1980 };
1981 return result;
1982}
1983
1985 Model* model,
1986 std::vector<IntegerVariable>* vars) {
1987 IntegerEncoder* encoder = model->GetOrCreate<IntegerEncoder>();
1988 for (int t = 0; t < helper->NumTasks(); ++t) {
1989 if (helper->Starts()[t].var != kNoIntegerVariable) {
1990 vars->push_back(helper->Starts()[t].var);
1991 }
1992 if (helper->Sizes()[t].var != kNoIntegerVariable) {
1993 vars->push_back(helper->Sizes()[t].var);
1994 }
1995 if (helper->Ends()[t].var != kNoIntegerVariable) {
1996 vars->push_back(helper->Ends()[t].var);
1997 }
1998 if (helper->IsOptional(t) && !helper->IsAbsent(t) &&
1999 !helper->IsPresent(t)) {
2000 const Literal l = helper->PresenceLiteral(t);
2001 if (encoder->GetLiteralView(l) == kNoIntegerVariable &&
2002 encoder->GetLiteralView(l.Negated()) == kNoIntegerVariable) {
2004 }
2005 const IntegerVariable direct_view = encoder->GetLiteralView(l);
2006 if (direct_view != kNoIntegerVariable) {
2007 vars->push_back(direct_view);
2008 } else {
2009 vars->push_back(encoder->GetLiteralView(l.Negated()));
2010 DCHECK_NE(vars->back(), kNoIntegerVariable);
2011 }
2012 }
2013 }
2015}
2016
2017std::function<void(const absl::StrongVector<IntegerVariable, double>&,
2018 LinearConstraintManager*)>
2019GenerateCumulativeCut(const std::string& cut_name,
2021 const std::vector<IntegerVariable>& demands,
2023 Trail* trail = model->GetOrCreate<Trail>();
2024 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
2025 IntegerEncoder* encoder = model->GetOrCreate<IntegerEncoder>();
2026
2027 return [capacity, demands, trail, integer_trail, helper, model, cut_name,
2028 encoder](const absl::StrongVector<IntegerVariable, double>& lp_values,
2029 LinearConstraintManager* manager) {
2030 if (trail->CurrentDecisionLevel() > 0) return;
2031
2032 const auto demand_is_fixed = [integer_trail, &demands](int i) {
2033 return demands.empty() || integer_trail->IsFixed(demands[i]);
2034 };
2035 const auto demand_min = [integer_trail, &demands](int i) {
2036 return demands.empty() ? IntegerValue(1)
2037 : integer_trail->LowerBound(demands[i]);
2038 };
2039 const auto demand_max = [integer_trail, &demands](int i) {
2040 return demands.empty() ? IntegerValue(1)
2041 : integer_trail->UpperBound(demands[i]);
2042 };
2043
2044 std::vector<int> active_intervals;
2045 for (int i = 0; i < helper->NumTasks(); ++i) {
2046 if (!helper->IsAbsent(i) && demand_max(i) > 0 && helper->SizeMin(i) > 0) {
2047 active_intervals.push_back(i);
2048 }
2049 }
2050
2051 if (active_intervals.size() < 2) return;
2052
2053 std::sort(active_intervals.begin(), active_intervals.end(),
2054 [helper](int a, int b) {
2055 return helper->StartMin(a) < helper->StartMin(b) ||
2056 (helper->StartMin(a) == helper->StartMin(b) &&
2057 helper->EndMax(a) < helper->EndMax(b));
2058 });
2059
2060 const IntegerValue capacity_max = integer_trail->UpperBound(capacity);
2061 IntegerValue processed_start = kMinIntegerValue;
2062 for (int i1 = 0; i1 + 1 < active_intervals.size(); ++i1) {
2063 const int start_index = active_intervals[i1];
2064 DCHECK(!helper->IsAbsent(start_index));
2065
2066 // We want maximal cuts. For any start_min value, we only need to create
2067 // cuts starting from the first interval having this start_min value.
2068 if (helper->StartMin(start_index) == processed_start) {
2069 continue;
2070 } else {
2071 processed_start = helper->StartMin(start_index);
2072 }
2073
2074 // For each start time, we will keep the most violated cut generated while
2075 // scanning the residual tasks.
2076 int end_index_of_max_violation = -1;
2077 double max_relative_violation = 1.01;
2078 IntegerValue span_of_max_violation(0);
2079
2080 // Accumulate intervals and check for potential cuts.
2081 double energy_lp = 0.0;
2082 IntegerValue min_of_starts = kMaxIntegerValue;
2083 IntegerValue max_of_ends = kMinIntegerValue;
2084
2085 // We sort all tasks (start_min(task) >= start_min(start_index) by
2086 // increasing end max.
2087 std::vector<int> residual_tasks(active_intervals.begin() + i1,
2088 active_intervals.end());
2089 std::sort(
2090 residual_tasks.begin(), residual_tasks.end(),
2091 [&](int a, int b) { return helper->EndMax(a) < helper->EndMax(b); });
2092
2093 // Let's process residual tasks and evaluate the cut violation of the cut
2094 // at each step. We follow the same structure as the cut creation code
2095 // below.
2096 for (int i2 = 0; i2 < residual_tasks.size(); ++i2) {
2097 const int t = residual_tasks[i2];
2098 if (helper->IsPresent(t)) {
2099 if (demand_is_fixed(t)) {
2100 if (helper->SizeIsFixed(t)) {
2101 energy_lp += ToDouble(helper->SizeMin(t) * demand_min(t));
2102 } else {
2103 energy_lp += ToDouble(demand_min(t)) *
2104 helper->Sizes()[t].LpValue(lp_values);
2105 }
2106 } else if (helper->SizeIsFixed(t)) {
2107 DCHECK(!demands.empty());
2108 energy_lp += lp_values[demands[t]] * ToDouble(helper->SizeMin(t));
2109 } else { // demand and size are not fixed.
2110 DCHECK(!demands.empty());
2111 energy_lp +=
2112 ToDouble(demand_min(t)) * helper->Sizes()[t].LpValue(lp_values);
2113 energy_lp += lp_values[demands[t]] * ToDouble(helper->SizeMin(t));
2114 energy_lp -= ToDouble(demand_min(t) * helper->SizeMin(t));
2115 }
2116 } else {
2117 energy_lp += GetLiteralLpValue(helper->PresenceLiteral(t), lp_values,
2118 encoder) *
2119 ToDouble(helper->SizeMin(t) * demand_min(t));
2120 }
2121
2122 min_of_starts = std::min(min_of_starts, helper->StartMin(t));
2123 max_of_ends = std::max(max_of_ends, helper->EndMax(t));
2124
2125 // Compute the violation of the potential cut.
2126 const double relative_violation =
2127 energy_lp / ToDouble((max_of_ends - min_of_starts) * capacity_max);
2128 if (relative_violation > max_relative_violation) {
2129 end_index_of_max_violation = i2;
2130 max_relative_violation = relative_violation;
2131 span_of_max_violation = max_of_ends - min_of_starts;
2132 }
2133 }
2134
2135 if (end_index_of_max_violation == -1) continue;
2136
2137 // A maximal violated cut has been found.
2138 bool cut_generated = true;
2139 bool has_opt_cuts = false;
2140 bool has_quadratic_cuts = false;
2141
2142 LinearConstraintBuilder cut(model, kMinIntegerValue, IntegerValue(0));
2143
2144 // Build the cut.
2145 cut.AddTerm(capacity, -span_of_max_violation);
2146 for (int i2 = 0; i2 <= end_index_of_max_violation; ++i2) {
2147 const int t = residual_tasks[i2];
2148 if (helper->IsPresent(t)) {
2149 if (demand_is_fixed(t)) {
2150 if (helper->SizeIsFixed(t)) {
2151 cut.AddConstant(helper->SizeMin(t) * demand_min(t));
2152 } else {
2153 cut.AddTerm(helper->Sizes()[t], demand_min(t));
2154 }
2155 } else if (helper->SizeIsFixed(t)) {
2156 DCHECK(!demands.empty());
2157 cut.AddTerm(demands[t], helper->SizeMin(t));
2158 } else { // demand and size are not fixed.
2159 DCHECK(!demands.empty());
2160 // We use McCormick equation.
2161 // demand * size = (demand_min + delta_d) * (min_size +
2162 // delta_s) =
2163 // demand_min * min_size + delta_d * min_size +
2164 // delta_s * demand_min + delta_s * delta_d
2165 // which is >= (by ignoring the quatratic term)
2166 // demand_min * size + min_size * demand - demand_min *
2167 // min_size
2168 cut.AddTerm(helper->Sizes()[t], demand_min(t));
2169 cut.AddTerm(demands[t], helper->SizeMin(t));
2170 // Substract the energy counted twice.
2171 cut.AddConstant(-helper->SizeMin(t) * demand_min(t));
2172 has_quadratic_cuts = true;
2173 }
2174 } else {
2175 has_opt_cuts = true;
2176 if (!helper->SizeIsFixed(t) || !demand_is_fixed(t)) {
2177 has_quadratic_cuts = true;
2178 }
2179 if (!cut.AddLiteralTerm(helper->PresenceLiteral(t),
2180 helper->SizeMin(t) * demand_min(t))) {
2181 cut_generated = false;
2182 break;
2183 }
2184 }
2185 }
2186
2187 if (cut_generated) {
2188 std::string full_name = cut_name;
2189 if (has_opt_cuts) full_name.append("_opt");
2190 if (has_quadratic_cuts) full_name.append("_quad");
2191
2192 manager->AddCut(cut.Build(), cut_name, lp_values);
2193 }
2194 }
2195 };
2196}
2197
2199 const std::vector<IntervalVariable>& intervals,
2200 const IntegerVariable capacity, const std::vector<IntegerVariable>& demands,
2201 Model* model) {
2202 CutGenerator result;
2203
2205 new SchedulingConstraintHelper(intervals, model);
2206 model->TakeOwnership(helper);
2207
2208 result.vars = demands;
2209 result.vars.push_back(capacity);
2210 AddIntegerVariableFromIntervals(helper, model, &result.vars);
2211
2213 "CumulativeEnergy", helper, demands, AffineExpression(capacity), model);
2214 return result;
2215}
2216
2218 const std::vector<IntervalVariable>& intervals,
2219 const IntegerVariable capacity, const std::vector<IntegerVariable>& demands,
2220 Model* model) {
2221 CutGenerator result;
2222
2224 new SchedulingConstraintHelper(intervals, model);
2225 model->TakeOwnership(helper);
2226
2227 result.vars = demands;
2228 result.vars.push_back(capacity);
2229 AddIntegerVariableFromIntervals(helper, model, &result.vars);
2230
2231 struct Event {
2232 int interval_index;
2233 IntegerValue time;
2234 bool positive;
2235 IntegerVariable demand;
2236 };
2237
2238 Trail* trail = model->GetOrCreate<Trail>();
2239 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
2240
2241 result.generate_cuts =
2242 [helper, capacity, demands, trail, integer_trail, model](
2244 LinearConstraintManager* manager) {
2245 if (trail->CurrentDecisionLevel() > 0) return;
2246
2247 std::vector<Event> events;
2248 // Iterate through the intervals. If start_max < end_min, the demand
2249 // is mandatory.
2250 for (int i = 0; i < helper->NumTasks(); ++i) {
2251 if (helper->IsAbsent(i)) continue;
2252
2253 const IntegerValue start_max = helper->StartMax(i);
2254 const IntegerValue end_min = helper->EndMin(i);
2255
2256 if (start_max >= end_min) continue;
2257
2258 Event e1;
2259 e1.interval_index = i;
2260 e1.time = start_max;
2261 e1.demand = demands[i];
2262 e1.positive = true;
2263
2264 Event e2 = e1;
2265 e2.time = end_min;
2266 e2.positive = false;
2267 events.push_back(e1);
2268 events.push_back(e2);
2269 }
2270
2271 // Sort events by time.
2272 // It is also important that all positive event with the same time as
2273 // negative events appear after for the correctness of the algo below.
2274 std::sort(events.begin(), events.end(),
2275 [](const Event i, const Event j) {
2276 if (i.time == j.time) {
2277 if (i.positive == j.positive) {
2278 return i.interval_index < j.interval_index;
2279 }
2280 return !i.positive;
2281 }
2282 return i.time < j.time;
2283 });
2284
2285 std::vector<Event> cut_events;
2286 bool added_positive_event = false;
2287 for (const Event& e : events) {
2288 if (e.positive) {
2289 added_positive_event = true;
2290 cut_events.push_back(e);
2291 continue;
2292 }
2293 if (added_positive_event && cut_events.size() > 1) {
2294 // Create cut.
2295 bool cut_generated = true;
2297 IntegerValue(0));
2298 cut.AddTerm(capacity, IntegerValue(-1));
2299 for (const Event& cut_event : cut_events) {
2300 if (helper->IsPresent(cut_event.interval_index)) {
2301 cut.AddTerm(cut_event.demand, IntegerValue(1));
2302 } else {
2303 cut_generated &= cut.AddLiteralTerm(
2304 helper->PresenceLiteral(cut_event.interval_index),
2305 integer_trail->LowerBound(cut_event.demand));
2306 if (!cut_generated) break;
2307 }
2308 }
2309 if (cut_generated) {
2310 // Violation of the cut is checked by AddCut so we don't check
2311 // it here.
2312 manager->AddCut(cut.Build(), "Cumulative", lp_values);
2313 }
2314 }
2315 // Remove the event.
2316 int new_size = 0;
2317 for (int i = 0; i < cut_events.size(); ++i) {
2318 if (cut_events[i].interval_index == e.interval_index) {
2319 continue;
2320 }
2321 cut_events[new_size] = cut_events[i];
2322 new_size++;
2323 }
2324 cut_events.resize(new_size);
2325 added_positive_event = false;
2326 }
2327 };
2328 return result;
2329}
2330
2332 const std::vector<IntervalVariable>& intervals, Model* model) {
2333 CutGenerator result;
2334
2336 new SchedulingConstraintHelper(intervals, model);
2337 model->TakeOwnership(helper);
2338
2339 AddIntegerVariableFromIntervals(helper, model, &result.vars);
2340
2342 "NoOverlapEnergy", helper,
2343 /*demands=*/{},
2344 /*capacity=*/AffineExpression(IntegerValue(1)), model);
2345 return result;
2346}
2347
2349 const std::vector<IntervalVariable>& intervals, Model* model) {
2350 CutGenerator result;
2351
2353 new SchedulingConstraintHelper(intervals, model);
2354 model->TakeOwnership(helper);
2355
2356 AddIntegerVariableFromIntervals(helper, model, &result.vars);
2357
2358 Trail* trail = model->GetOrCreate<Trail>();
2359
2360 result.generate_cuts =
2361 [trail, helper, model](
2363 LinearConstraintManager* manager) {
2364 if (trail->CurrentDecisionLevel() > 0) return;
2365
2366 // TODO(user): We can do much better in term of complexity:
2367 // Sort all tasks by min start time, loop other them 1 by 1,
2368 // start scanning their successors and stop when the start time of the
2369 // successor is >= duration min of the task.
2370
2371 // TODO(user): each time we go back to level zero, we will generate
2372 // the same cuts over and over again. It is okay because AddCut() will
2373 // not add duplicate cuts, but it might not be the most efficient way.
2374 for (int index1 = 0; index1 < helper->NumTasks(); ++index1) {
2375 if (!helper->IsPresent(index1)) continue;
2376 for (int index2 = index1 + 1; index2 < helper->NumTasks(); ++index2) {
2377 if (!helper->IsPresent(index2)) continue;
2378
2379 // Encode only the interesting pairs.
2380 if (helper->EndMax(index1) <= helper->StartMin(index2) ||
2381 helper->EndMax(index2) <= helper->StartMin(index1)) {
2382 continue;
2383 }
2384
2385 const bool interval_1_can_precede_2 =
2386 helper->EndMin(index1) <= helper->StartMax(index2);
2387 const bool interval_2_can_precede_1 =
2388 helper->EndMin(index2) <= helper->StartMax(index1);
2389
2390 if (interval_1_can_precede_2 && !interval_2_can_precede_1) {
2391 // interval1.end <= interval2.start
2393 IntegerValue(0));
2394 cut.AddTerm(helper->Ends()[index1], IntegerValue(1));
2395 cut.AddTerm(helper->Starts()[index2], IntegerValue(-1));
2396 manager->AddCut(cut.Build(), "NoOverlapPrecedence", lp_values);
2397 } else if (interval_2_can_precede_1 && !interval_1_can_precede_2) {
2398 // interval2.end <= interval1.start
2400 IntegerValue(0));
2401 cut.AddTerm(helper->Ends()[index2], IntegerValue(1));
2402 cut.AddTerm(helper->Starts()[index1], IntegerValue(-1));
2403 manager->AddCut(cut.Build(), "NoOverlapPrecedence", lp_values);
2404 }
2405 }
2406 }
2407 };
2408 return result;
2409}
2410
2412 const std::vector<IntegerVariable>& base_variables, Model* model) {
2413 // Filter base_variables to only keep the one with a literal view, and
2414 // do the conversion.
2415 std::vector<IntegerVariable> variables;
2416 std::vector<Literal> literals;
2417 absl::flat_hash_map<LiteralIndex, IntegerVariable> positive_map;
2418 absl::flat_hash_map<LiteralIndex, IntegerVariable> negative_map;
2419 auto* integer_trail = model->GetOrCreate<IntegerTrail>();
2420 auto* encoder = model->GetOrCreate<IntegerEncoder>();
2421 for (const IntegerVariable var : base_variables) {
2422 if (integer_trail->LowerBound(var) != IntegerValue(0)) continue;
2423 if (integer_trail->UpperBound(var) != IntegerValue(1)) continue;
2424 const LiteralIndex literal_index = encoder->GetAssociatedLiteral(
2425 IntegerLiteral::GreaterOrEqual(var, IntegerValue(1)));
2426 if (literal_index != kNoLiteralIndex) {
2427 variables.push_back(var);
2428 literals.push_back(Literal(literal_index));
2429 positive_map[literal_index] = var;
2430 negative_map[Literal(literal_index).NegatedIndex()] = var;
2431 }
2432 }
2433 CutGenerator result;
2434 result.vars = variables;
2435 auto* implication_graph = model->GetOrCreate<BinaryImplicationGraph>();
2436 result.generate_cuts =
2437 [variables, literals, implication_graph, positive_map, negative_map,
2439 LinearConstraintManager* manager) {
2440 std::vector<double> packed_values;
2441 for (int i = 0; i < literals.size(); ++i) {
2442 packed_values.push_back(lp_values[variables[i]]);
2443 }
2444 const std::vector<std::vector<Literal>> at_most_ones =
2445 implication_graph->GenerateAtMostOnesWithLargeWeight(literals,
2446 packed_values);
2447
2448 for (const std::vector<Literal>& at_most_one : at_most_ones) {
2449 // We need to express such "at most one" in term of the initial
2450 // variables, so we do not use the
2451 // LinearConstraintBuilder::AddLiteralTerm() here.
2452 LinearConstraintBuilder builder(model, IntegerValue(kint64min),
2453 IntegerValue(1));
2454 for (const Literal l : at_most_one) {
2455 if (ContainsKey(positive_map, l.Index())) {
2456 builder.AddTerm(positive_map.at(l.Index()), IntegerValue(1));
2457 } else {
2458 // Add 1 - X to the linear constraint.
2459 builder.AddTerm(negative_map.at(l.Index()), IntegerValue(-1));
2460 builder.AddConstant(IntegerValue(1));
2461 }
2462 }
2463
2464 manager->AddCut(builder.Build(), "clique", lp_values);
2465 }
2466 };
2467 return result;
2468}
2469
2470} // namespace sat
2471} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define CHECK(condition)
Definition: base/logging.h:495
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define CHECK_LT(val1, val2)
Definition: base/logging.h:700
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define CHECK_GE(val1, val2)
Definition: base/logging.h:701
#define CHECK_GT(val1, val2)
Definition: base/logging.h:702
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define CHECK_NE(val1, val2)
Definition: base/logging.h:698
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define VLOG(verboselevel)
Definition: base/logging.h:978
We call domain any subset of Int64 = [kint64min, kint64max].
Domain UnionWith(const Domain &domain) const
Returns the union of D and domain.
void Init(const std::vector< double > &profits, const std::vector< double > &weights, const double capacity)
double Solve(TimeLimit *time_limit, bool *is_solution_optimal)
void set_solution_upper_bound_threshold(const double solution_upper_bound_threshold)
A simple class to enforce both an elapsed time limit and a deterministic time limit in the same threa...
Definition: time_limit.h:105
bool LimitReached()
Returns true when the external limit is true, or the deterministic time is over the deterministic lim...
Definition: time_limit.h:532
bool TrySimpleKnapsack(const LinearConstraint base_ct, const std::vector< double > &lp_values, const std::vector< IntegerValue > &lower_bounds, const std::vector< IntegerValue > &upper_bounds)
Definition: cuts.cc:1155
void ProcessUpperBoundedConstraint(const absl::StrongVector< IntegerVariable, double > &lp_values, LinearConstraint *cut)
Definition: cuts.cc:1491
BestImpliedBoundInfo GetCachedImpliedBoundInfo(IntegerVariable var)
Definition: cuts.cc:1500
const IntegerVariable GetLiteralView(Literal lit) const
Definition: integer.h:420
void ComputeCut(RoundingOptions options, const std::vector< double > &lp_values, const std::vector< IntegerValue > &lower_bounds, const std::vector< IntegerValue > &upper_bounds, ImpliedBoundsProcessor *ib_processor, LinearConstraint *cut)
Definition: cuts.cc:707
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
Definition: integer.h:1355
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
Definition: integer.h:1350
const Domain & InitialVariableDomain(IntegerVariable var) const
Definition: integer.cc:644
ABSL_MUST_USE_RESULT bool AddLiteralTerm(Literal lit, IntegerValue coeff)
void AddTerm(IntegerVariable var, IntegerValue coeff)
bool AddCut(LinearConstraint ct, std::string type_name, const absl::StrongVector< IntegerVariable, double > &lp_solution, std::string extra_info="")
LiteralIndex NegatedIndex() const
Definition: sat_base.h:85
Class that owns everything related to a particular optimization model.
Definition: sat/model.h:38
const std::vector< AffineExpression > & Starts() const
Definition: intervals.h:319
const std::vector< AffineExpression > & Sizes() const
Definition: intervals.h:321
const std::vector< AffineExpression > & Ends() const
Definition: intervals.h:320
SharedTimeLimit * time_limit
int64 value
IntVar * var
Definition: expr_array.cc:1858
GRBmodel * model
static const int64 kint64max
int64_t int64
static const int64 kint64min
const int INFO
Definition: log_severity.h:31
bool ContainsKey(const Collection &collection, const Key &key)
Definition: map_util.h:170
void STLSortAndRemoveDuplicates(T *v, const LessFunc &less_func)
Definition: stl_util.h:58
static double ToDouble(double f)
Definition: lp_types.h:68
std::function< void(const absl::StrongVector< IntegerVariable, double > &, LinearConstraintManager *)> GenerateCumulativeCut(const std::string &cut_name, SchedulingConstraintHelper *helper, const std::vector< IntegerVariable > &demands, AffineExpression capacity, Model *model)
Definition: cuts.cc:2019
void ConvertToKnapsackForm(const LinearConstraint &constraint, std::vector< LinearConstraint > *knapsack_constraints, IntegerTrail *integer_trail)
Definition: cuts.cc:388
void CleanTermsAndFillConstraint(std::vector< std::pair< IntegerVariable, IntegerValue > > *terms, LinearConstraint *constraint)
bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue *result)
Definition: integer.h:110
LinearConstraint GetPreprocessedLinearConstraint(const LinearConstraint &constraint, const absl::StrongVector< IntegerVariable, double > &lp_values, const IntegerTrail &integer_trail)
Definition: cuts.cc:250
CutGenerator CreateNoOverlapPrecedenceCutGenerator(const std::vector< IntervalVariable > &intervals, Model *model)
Definition: cuts.cc:2348
bool CanFormValidKnapsackCover(const LinearConstraint &preprocessed_constraint, const absl::StrongVector< IntegerVariable, double > &lp_values, const IntegerTrail &integer_trail)
Definition: cuts.cc:370
CutGenerator CreateCumulativeCutGenerator(const std::vector< IntervalVariable > &intervals, const IntegerVariable capacity, const std::vector< IntegerVariable > &demands, Model *model)
Definition: cuts.cc:2198
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue)
const LiteralIndex kNoLiteralIndex(-1)
CutGenerator CreateNoOverlapCutGenerator(const std::vector< IntervalVariable > &intervals, Model *model)
Definition: cuts.cc:2331
void RemoveZeroTerms(LinearConstraint *constraint)
IntegerValue GetFactorT(IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue max_t)
Definition: cuts.cc:616
double GetKnapsackUpperBound(std::vector< KnapsackItem > items, const double capacity)
Definition: cuts.cc:318
CutGenerator CreateOverlappingCumulativeCutGenerator(const std::vector< IntervalVariable > &intervals, const IntegerVariable capacity, const std::vector< IntegerVariable > &demands, Model *model)
Definition: cuts.cc:2217
std::function< IntegerVariable(Model *)> NewIntegerVariableFromLiteral(Literal lit)
Definition: integer.h:1444
CutGenerator CreateSquareCutGenerator(IntegerVariable y, IntegerVariable x, Model *model)
Definition: cuts.cc:1424
bool CanBeFilteredUsingCutLowerBound(const LinearConstraint &preprocessed_constraint, const absl::StrongVector< IntegerVariable, double > &lp_values, const IntegerTrail &integer_trail)
Definition: cuts.cc:290
void MakeAllCoefficientsPositive(LinearConstraint *constraint)
const IntegerVariable kNoIntegerVariable(-1)
IntegerVariable PositiveVariable(IntegerVariable i)
Definition: integer.h:134
std::function< void(Model *)> GreaterOrEqual(IntegerVariable v, int64 lb)
Definition: integer.h:1495
CutGenerator CreateLinMaxCutGenerator(const IntegerVariable target, const std::vector< LinearExpression > &exprs, const std::vector< IntegerVariable > &z_vars, Model *model)
Definition: cuts.cc:1915
CutGenerator CreateAllDifferentCutGenerator(const std::vector< IntegerVariable > &vars, Model *model)
Definition: cuts.cc:1818
bool CanBeFilteredUsingKnapsackUpperBound(const LinearConstraint &constraint, const absl::StrongVector< IntegerVariable, double > &lp_values, const IntegerTrail &integer_trail)
Definition: cuts.cc:336
std::function< IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction(IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue t, IntegerValue max_scaling)
Definition: cuts.cc:624
IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor)
Definition: integer.h:90
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
void MakeAllVariablesPositive(LinearConstraint *constraint)
double ToDouble(IntegerValue value)
Definition: integer.h:69
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
Definition: integer.cc:27
IntegerValue GetCoefficientOfPositiveVar(const IntegerVariable var, const LinearExpression &expr)
IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor)
Definition: integer.h:81
CutGenerator CreateKnapsackCoverCutGenerator(const std::vector< LinearConstraint > &base_constraints, const std::vector< IntegerVariable > &vars, Model *model)
Definition: cuts.cc:437
bool ConstraintIsTriviallyTrue(const LinearConstraint &constraint, const IntegerTrail &integer_trail)
Definition: cuts.cc:274
bool LiftKnapsackCut(const LinearConstraint &constraint, const absl::StrongVector< IntegerVariable, double > &lp_values, const std::vector< IntegerValue > &cut_vars_original_coefficients, const IntegerTrail &integer_trail, TimeLimit *time_limit, LinearConstraint *cut)
Definition: cuts.cc:172
CutGenerator CreatePositiveMultiplicationCutGenerator(IntegerVariable z, IntegerVariable x, IntegerVariable y, Model *model)
Definition: cuts.cc:1328
void AddIntegerVariableFromIntervals(SchedulingConstraintHelper *helper, Model *model, std::vector< IntegerVariable > *vars)
Definition: cuts.cc:1984
CutGenerator CreateCliqueCutGenerator(const std::vector< IntegerVariable > &base_variables, Model *model)
Definition: cuts.cc:2411
void DivideByGCD(LinearConstraint *constraint)
double ComputeActivity(const LinearConstraint &constraint, const absl::StrongVector< IntegerVariable, double > &values)
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapAdd(int64 x, int64 y)
int64 FloorRatio(int64 value, int64 positive_coeff)
int64 SumOfKMinValueInDomain(const Domain &domain, int k)
int64 CapProd(int64 x, int64 y)
int64 CapSub(int64 x, int64 y)
int64 SumOfKMaxValueInDomain(const Domain &domain, int k)
int index
Definition: pack.cc:508
int64 weight
Definition: pack.cc:509
int64 time
Definition: resource.cc:1683
int64 demand
Definition: resource.cc:123
Fractional ratio
int64 capacity
int64 coefficient
std::vector< double > lower_bounds
std::vector< double > upper_bounds
Rev< int64 > end_min
Rev< int64 > start_max
std::vector< IntegerVariable > vars
Definition: cuts.h:42
std::function< void(const absl::StrongVector< IntegerVariable, double > &lp_values, LinearConstraintManager *manager)> generate_cuts
Definition: cuts.h:46
std::vector< std::pair< IntegerVariable, IntegerValue > > terms
Definition: cuts.h:80
std::vector< IntegerVariable > vars
void AddTerm(IntegerVariable var, IntegerValue coeff)