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 
25 #include "ortools/base/stl_util.h"
27 #include "ortools/sat/integer.h"
28 #include "ortools/sat/intervals.h"
30 #include "ortools/sat/sat_base.h"
32 
33 namespace operations_research {
34 namespace sat {
35 
36 namespace {
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.
40 const double kMinCutViolation = 1e-4;
41 
42 // Returns the lp value of a Literal.
43 double 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.
59 LinearConstraint 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 
87 bool 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 
98 bool 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 
117 bool 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.
137 int 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 
160 bool 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 
318 double 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;
551  knapsack_solver.set_solution_upper_bound_threshold(
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.
616 IntegerValue 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 
624 std::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 
936  const auto f = GetSuperAdditiveRoundingFunction(
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 
1424 CutGenerator 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 
1491 void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint(
1493  LinearConstraint* cut) {
1494  ProcessUpperBoundedConstraintWithSlackCreation(
1495  /*substitute_only_inner_variables=*/false, IntegerVariable(0), lp_values,
1496  cut, nullptr);
1497 }
1498 
1500 ImpliedBoundsProcessor::GetCachedImpliedBoundInfo(IntegerVariable var) {
1501  auto it = cache_.find(var);
1502  if (it != cache_.end()) return it->second;
1503  return BestImpliedBoundInfo();
1504 }
1505 
1507 ImpliedBoundsProcessor::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.
1575 void 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 
1584 void 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 
1725 bool 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 
1777 namespace {
1778 
1779 void 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 
1852 namespace {
1853 // Returns max((w2i - w1i)*Li, (w2i - w1i)*Ui).
1854 IntegerValue 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
1869 IntegerValue 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).
1891 double 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) {
1938  absl::StrongVector<IntegerVariable, int> variable_partition(
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 
2017 std::function<void(const absl::StrongVector<IntegerVariable, double>&,
2018  LinearConstraintManager*)>
2019 GenerateCumulativeCut(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 
2204  SchedulingConstraintHelper* helper =
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 
2223  SchedulingConstraintHelper* helper =
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 
2335  SchedulingConstraintHelper* helper =
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 
2352  SchedulingConstraintHelper* helper =
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
void STLSortAndRemoveDuplicates(T *v, const LessFunc &less_func)
Definition: stl_util.h:58
bool ContainsKey(const Collection &collection, const Key &key)
Definition: map_util.h:170
static double ToDouble(double f)
Definition: lp_types.h:68
void ConvertToKnapsackForm(const LinearConstraint &constraint, std::vector< LinearConstraint > *knapsack_constraints, IntegerTrail *integer_trail)
Definition: cuts.cc:388
IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor)
Definition: integer.h:90
bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue *result)
Definition: integer.h:110
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
LinearConstraint GetPreprocessedLinearConstraint(const LinearConstraint &constraint, const absl::StrongVector< IntegerVariable, double > &lp_values, const IntegerTrail &integer_trail)
Definition: cuts.cc:250
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
CutGenerator CreateNoOverlapPrecedenceCutGenerator(const std::vector< IntervalVariable > &intervals, Model *model)
Definition: cuts.cc:2348
IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor)
Definition: integer.h:81
const LiteralIndex kNoLiteralIndex(-1)
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)
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
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
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
const IntegerVariable kNoIntegerVariable(-1)
void MakeAllCoefficientsPositive(LinearConstraint *constraint)
std::function< IntegerVariable(Model *)> NewIntegerVariableFromLiteral(Literal lit)
Definition: integer.h:1444
IntegerVariable PositiveVariable(IntegerVariable i)
Definition: integer.h:134
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
void MakeAllVariablesPositive(LinearConstraint *constraint)
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
Definition: integer.cc:27
IntegerValue GetCoefficientOfPositiveVar(const IntegerVariable var, const LinearExpression &expr)
CutGenerator CreateKnapsackCoverCutGenerator(const std::vector< LinearConstraint > &base_constraints, const std::vector< IntegerVariable > &vars, Model *model)
Definition: cuts.cc:437
void AddIntegerVariableFromIntervals(SchedulingConstraintHelper *helper, Model *model, std::vector< IntegerVariable > *vars)
Definition: cuts.cc:1984
bool ConstraintIsTriviallyTrue(const LinearConstraint &constraint, const IntegerTrail &integer_trail)
Definition: cuts.cc:274
std::function< void(Model *)> GreaterOrEqual(IntegerVariable v, int64 lb)
Definition: integer.h:1495
void CleanTermsAndFillConstraint(std::vector< std::pair< IntegerVariable, IntegerValue >> *terms, LinearConstraint *constraint)
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
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)
double ToDouble(IntegerValue value)
Definition: integer.h:69
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapSub(int64 x, int64 y)
int64 SumOfKMinValueInDomain(const Domain &domain, int k)
int64 FloorRatio(int64 value, int64 positive_coeff)
int64 CapAdd(int64 x, int64 y)
int64 CapProd(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)