OR-Tools  8.2
linear_programming_constraint.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 
15 
16 #include <algorithm>
17 #include <cmath>
18 #include <iterator>
19 #include <limits>
20 #include <string>
21 #include <utility>
22 #include <vector>
23 
24 #include "absl/container/flat_hash_map.h"
25 #include "absl/numeric/int128.h"
28 #include "ortools/base/logging.h"
29 #include "ortools/base/map_util.h"
30 #include "ortools/base/mathutil.h"
31 #include "ortools/base/stl_util.h"
33 #include "ortools/glop/parameters.pb.h"
35 #include "ortools/glop/status.h"
39 #include "ortools/sat/integer.h"
42 
43 namespace operations_research {
44 namespace sat {
45 
46 using glop::ColIndex;
47 using glop::Fractional;
48 using glop::RowIndex;
49 
51  if (is_sparse_) {
52  for (const glop::ColIndex col : non_zeros_) {
53  dense_vector_[col] = IntegerValue(0);
54  }
55  dense_vector_.resize(size, IntegerValue(0));
56  } else {
57  dense_vector_.assign(size, IntegerValue(0));
58  }
59  for (const glop::ColIndex col : non_zeros_) {
60  is_zeros_[col] = true;
61  }
62  is_zeros_.resize(size, true);
63  non_zeros_.clear();
64  is_sparse_ = true;
65 }
66 
67 bool ScatteredIntegerVector::Add(glop::ColIndex col, IntegerValue value) {
68  const int64 add = CapAdd(value.value(), dense_vector_[col].value());
69  if (add == kint64min || add == kint64max) return false;
70  dense_vector_[col] = IntegerValue(add);
71  if (is_sparse_ && is_zeros_[col]) {
72  is_zeros_[col] = false;
73  non_zeros_.push_back(col);
74  }
75  return true;
76 }
77 
79  IntegerValue multiplier,
80  const std::vector<std::pair<glop::ColIndex, IntegerValue>>& terms) {
81  const double threshold = 0.1 * static_cast<double>(dense_vector_.size());
82  if (is_sparse_ && static_cast<double>(terms.size()) < threshold) {
83  for (const std::pair<glop::ColIndex, IntegerValue> term : terms) {
84  if (is_zeros_[term.first]) {
85  is_zeros_[term.first] = false;
86  non_zeros_.push_back(term.first);
87  }
88  if (!AddProductTo(multiplier, term.second, &dense_vector_[term.first])) {
89  return false;
90  }
91  }
92  if (static_cast<double>(non_zeros_.size()) < threshold) {
93  is_sparse_ = false;
94  }
95  } else {
96  is_sparse_ = false;
97  for (const std::pair<glop::ColIndex, IntegerValue> term : terms) {
98  if (!AddProductTo(multiplier, term.second, &dense_vector_[term.first])) {
99  return false;
100  }
101  }
102  }
103  return true;
104 }
105 
107  const std::vector<IntegerVariable>& integer_variables,
108  IntegerValue upper_bound, LinearConstraint* result) {
109  result->vars.clear();
110  result->coeffs.clear();
111  if (is_sparse_) {
112  std::sort(non_zeros_.begin(), non_zeros_.end());
113  for (const glop::ColIndex col : non_zeros_) {
114  const IntegerValue coeff = dense_vector_[col];
115  if (coeff == 0) continue;
116  result->vars.push_back(integer_variables[col.value()]);
117  result->coeffs.push_back(coeff);
118  }
119  } else {
120  const int size = dense_vector_.size();
121  for (glop::ColIndex col(0); col < size; ++col) {
122  const IntegerValue coeff = dense_vector_[col];
123  if (coeff == 0) continue;
124  result->vars.push_back(integer_variables[col.value()]);
125  result->coeffs.push_back(coeff);
126  }
127  }
128  result->lb = kMinIntegerValue;
129  result->ub = upper_bound;
130 }
131 
132 std::vector<std::pair<glop::ColIndex, IntegerValue>>
134  std::vector<std::pair<glop::ColIndex, IntegerValue>> result;
135  if (is_sparse_) {
136  std::sort(non_zeros_.begin(), non_zeros_.end());
137  for (const glop::ColIndex col : non_zeros_) {
138  const IntegerValue coeff = dense_vector_[col];
139  if (coeff != 0) result.push_back({col, coeff});
140  }
141  } else {
142  const int size = dense_vector_.size();
143  for (glop::ColIndex col(0); col < size; ++col) {
144  const IntegerValue coeff = dense_vector_[col];
145  if (coeff != 0) result.push_back({col, coeff});
146  }
147  }
148  return result;
149 }
150 
151 // TODO(user): make SatParameters singleton too, otherwise changing them after
152 // a constraint was added will have no effect on this class.
154  : constraint_manager_(model),
155  sat_parameters_(*(model->GetOrCreate<SatParameters>())),
156  model_(model),
157  time_limit_(model->GetOrCreate<TimeLimit>()),
158  integer_trail_(model->GetOrCreate<IntegerTrail>()),
159  trail_(model->GetOrCreate<Trail>()),
160  integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
161  random_(model->GetOrCreate<ModelRandomGenerator>()),
162  implied_bounds_processor_({}, integer_trail_,
163  model->GetOrCreate<ImpliedBounds>()),
164  dispatcher_(model->GetOrCreate<LinearProgrammingDispatcher>()),
165  expanded_lp_solution_(
166  *model->GetOrCreate<LinearProgrammingConstraintLpSolution>()) {
167  // Tweak the default parameters to make the solve incremental.
168  glop::GlopParameters parameters;
169  parameters.set_use_dual_simplex(true);
170  simplex_.SetParameters(parameters);
171  if (sat_parameters_.use_branching_in_lp() ||
172  sat_parameters_.search_branching() == SatParameters::LP_SEARCH) {
173  compute_reduced_cost_averages_ = true;
174  }
175 
176  // Register our local rev int repository.
177  integer_trail_->RegisterReversibleClass(&rc_rev_int_repository_);
178 }
179 
181  VLOG(1) << "Total number of simplex iterations: "
182  << total_num_simplex_iterations_;
183  for (int i = 0; i < num_solves_by_status_.size(); ++i) {
184  if (num_solves_by_status_[i] == 0) continue;
185  VLOG(1) << "#" << glop::ProblemStatus(i) << " : "
186  << num_solves_by_status_[i];
187  }
188 }
189 
191  const LinearConstraint& ct) {
192  DCHECK(!lp_constraint_is_registered_);
193  constraint_manager_.Add(ct);
194 
195  // We still create the mirror variable right away though.
196  //
197  // TODO(user): clean this up? Note that it is important that the variable
198  // in lp_data_ never changes though, so we can restart from the current
199  // lp solution and be incremental (even if the constraints changed).
200  for (const IntegerVariable var : ct.vars) {
201  GetOrCreateMirrorVariable(PositiveVariable(var));
202  }
203 }
204 
205 glop::ColIndex LinearProgrammingConstraint::GetOrCreateMirrorVariable(
206  IntegerVariable positive_variable) {
207  DCHECK(VariableIsPositive(positive_variable));
208  const auto it = mirror_lp_variable_.find(positive_variable);
209  if (it == mirror_lp_variable_.end()) {
210  const glop::ColIndex col(integer_variables_.size());
211  implied_bounds_processor_.AddLpVariable(positive_variable);
212  mirror_lp_variable_[positive_variable] = col;
213  integer_variables_.push_back(positive_variable);
214  lp_solution_.push_back(std::numeric_limits<double>::infinity());
215  lp_reduced_cost_.push_back(0.0);
216  (*dispatcher_)[positive_variable] = this;
217 
218  const int index = std::max(positive_variable.value(),
219  NegationOf(positive_variable).value());
220  if (index >= expanded_lp_solution_.size()) {
221  expanded_lp_solution_.resize(index + 1, 0.0);
222  }
223  return col;
224  }
225  return it->second;
226 }
227 
229  IntegerValue coeff) {
230  CHECK(!lp_constraint_is_registered_);
231  objective_is_defined_ = true;
232  IntegerVariable pos_var = VariableIsPositive(ivar) ? ivar : NegationOf(ivar);
233  if (ivar != pos_var) coeff = -coeff;
234 
235  constraint_manager_.SetObjectiveCoefficient(pos_var, coeff);
236  const glop::ColIndex col = GetOrCreateMirrorVariable(pos_var);
237  integer_objective_.push_back({col, coeff});
238  objective_infinity_norm_ =
239  std::max(objective_infinity_norm_, IntTypeAbs(coeff));
240 }
241 
242 // TODO(user): As the search progress, some variables might get fixed. Exploit
243 // this to reduce the number of variables in the LP and in the
244 // ConstraintManager? We might also detect during the search that two variable
245 // are equivalent.
246 //
247 // TODO(user): On TSP/VRP with a lot of cuts, this can take 20% of the overall
248 // running time. We should be able to almost remove most of this from the
249 // profile by being more incremental (modulo LP scaling).
250 //
251 // TODO(user): A longer term idea for LP with a lot of variables is to not
252 // add all variables to each LP solve and do some "sifting". That can be useful
253 // for TSP for instance where the number of edges is large, but only a small
254 // fraction will be used in the optimal solution.
255 bool LinearProgrammingConstraint::CreateLpFromConstraintManager() {
256  // Fill integer_lp_.
257  integer_lp_.clear();
258  infinity_norms_.clear();
259  const auto& all_constraints = constraint_manager_.AllConstraints();
260  for (const auto index : constraint_manager_.LpConstraints()) {
261  const LinearConstraint& ct = all_constraints[index].constraint;
262 
263  integer_lp_.push_back(LinearConstraintInternal());
264  LinearConstraintInternal& new_ct = integer_lp_.back();
265  new_ct.lb = ct.lb;
266  new_ct.ub = ct.ub;
267  const int size = ct.vars.size();
268  IntegerValue infinity_norm(0);
269  if (ct.lb > ct.ub) {
270  VLOG(1) << "Trivial infeasible bound in an LP constraint";
271  return false;
272  }
273  if (ct.lb > kMinIntegerValue) {
274  infinity_norm = std::max(infinity_norm, IntTypeAbs(ct.lb));
275  }
276  if (ct.ub < kMaxIntegerValue) {
277  infinity_norm = std::max(infinity_norm, IntTypeAbs(ct.ub));
278  }
279  for (int i = 0; i < size; ++i) {
280  // We only use positive variable inside this class.
281  IntegerVariable var = ct.vars[i];
282  IntegerValue coeff = ct.coeffs[i];
283  if (!VariableIsPositive(var)) {
284  var = NegationOf(var);
285  coeff = -coeff;
286  }
287  infinity_norm = std::max(infinity_norm, IntTypeAbs(coeff));
288  new_ct.terms.push_back({GetOrCreateMirrorVariable(var), coeff});
289  }
290  infinity_norms_.push_back(infinity_norm);
291 
292  // Important to keep lp_data_ "clean".
293  std::sort(new_ct.terms.begin(), new_ct.terms.end());
294  }
295 
296  // Copy the integer_lp_ into lp_data_.
297  lp_data_.Clear();
298  for (int i = 0; i < integer_variables_.size(); ++i) {
299  CHECK_EQ(glop::ColIndex(i), lp_data_.CreateNewVariable());
300  }
301 
302  // We remove fixed variables from the objective. This should help the LP
303  // scaling, but also our integer reason computation.
304  int new_size = 0;
305  objective_infinity_norm_ = 0;
306  for (const auto entry : integer_objective_) {
307  const IntegerVariable var = integer_variables_[entry.first.value()];
308  if (integer_trail_->IsFixedAtLevelZero(var)) {
309  integer_objective_offset_ +=
310  entry.second * integer_trail_->LevelZeroLowerBound(var);
311  continue;
312  }
313  objective_infinity_norm_ =
314  std::max(objective_infinity_norm_, IntTypeAbs(entry.second));
315  integer_objective_[new_size++] = entry;
316  lp_data_.SetObjectiveCoefficient(entry.first, ToDouble(entry.second));
317  }
318  objective_infinity_norm_ =
319  std::max(objective_infinity_norm_, IntTypeAbs(integer_objective_offset_));
320  integer_objective_.resize(new_size);
321  lp_data_.SetObjectiveOffset(ToDouble(integer_objective_offset_));
322 
323  for (const LinearConstraintInternal& ct : integer_lp_) {
324  const ConstraintIndex row = lp_data_.CreateNewConstraint();
325  lp_data_.SetConstraintBounds(row, ToDouble(ct.lb), ToDouble(ct.ub));
326  for (const auto& term : ct.terms) {
327  lp_data_.SetCoefficient(row, term.first, ToDouble(term.second));
328  }
329  }
330  lp_data_.NotifyThatColumnsAreClean();
331 
332  // We scale the LP using the level zero bounds that we later override
333  // with the current ones.
334  //
335  // TODO(user): As part of the scaling, we may also want to shift the initial
336  // variable bounds so that each variable contain the value zero in their
337  // domain. Maybe just once and for all at the beginning.
338  const int num_vars = integer_variables_.size();
339  for (int i = 0; i < num_vars; i++) {
340  const IntegerVariable cp_var = integer_variables_[i];
341  const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var));
342  const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var));
343  lp_data_.SetVariableBounds(glop::ColIndex(i), lb, ub);
344  }
345 
346  // TODO(user): As we have an idea of the LP optimal after the first solves,
347  // maybe we can adapt the scaling accordingly.
348  glop::GlopParameters params;
349  params.set_cost_scaling(glop::GlopParameters::MEAN_COST_SCALING);
350  scaler_.Scale(params, &lp_data_);
351  UpdateBoundsOfLpVariables();
352 
353  // Set the information for the step to polish the LP basis. All our variables
354  // are integer, but for now, we just try to minimize the fractionality of the
355  // binary variables.
356  if (sat_parameters_.polish_lp_solution()) {
357  simplex_.ClearIntegralityScales();
358  for (int i = 0; i < num_vars; ++i) {
359  const IntegerVariable cp_var = integer_variables_[i];
360  const IntegerValue lb = integer_trail_->LevelZeroLowerBound(cp_var);
361  const IntegerValue ub = integer_trail_->LevelZeroUpperBound(cp_var);
362  if (lb != 0 || ub != 1) continue;
363  simplex_.SetIntegralityScale(
364  glop::ColIndex(i),
365  1.0 / scaler_.VariableScalingFactor(glop::ColIndex(i)));
366  }
367  }
368 
369  lp_data_.NotifyThatColumnsAreClean();
370  lp_data_.AddSlackVariablesWhereNecessary(false);
371  VLOG(1) << "LP relaxation: " << lp_data_.GetDimensionString() << ". "
372  << constraint_manager_.AllConstraints().size()
373  << " Managed constraints.";
374  return true;
375 }
376 
377 LPSolveInfo LinearProgrammingConstraint::SolveLpForBranching() {
378  LPSolveInfo info;
379  glop::BasisState basis_state = simplex_.GetState();
380 
381  const glop::Status status = simplex_.Solve(lp_data_, time_limit_);
382  total_num_simplex_iterations_ += simplex_.GetNumberOfIterations();
383  simplex_.LoadStateForNextSolve(basis_state);
384  if (!status.ok()) {
385  VLOG(1) << "The LP solver encountered an error: " << status.error_message();
386  info.status = glop::ProblemStatus::ABNORMAL;
387  return info;
388  }
389  info.status = simplex_.GetProblemStatus();
390  if (info.status == glop::ProblemStatus::OPTIMAL ||
391  info.status == glop::ProblemStatus::DUAL_FEASIBLE) {
392  // Record the objective bound.
393  info.lp_objective = simplex_.GetObjectiveValue();
394  info.new_obj_bound = IntegerValue(
395  static_cast<int64>(std::ceil(info.lp_objective - kCpEpsilon)));
396  }
397  return info;
398 }
399 
400 void LinearProgrammingConstraint::FillReducedCostReasonIn(
401  const glop::DenseRow& reduced_costs,
402  std::vector<IntegerLiteral>* integer_reason) {
403  integer_reason->clear();
404  const int num_vars = integer_variables_.size();
405  for (int i = 0; i < num_vars; i++) {
406  const double rc = reduced_costs[glop::ColIndex(i)];
407  if (rc > kLpEpsilon) {
408  integer_reason->push_back(
409  integer_trail_->LowerBoundAsLiteral(integer_variables_[i]));
410  } else if (rc < -kLpEpsilon) {
411  integer_reason->push_back(
412  integer_trail_->UpperBoundAsLiteral(integer_variables_[i]));
413  }
414  }
415 
416  integer_trail_->RemoveLevelZeroBounds(integer_reason);
417 }
418 
419 bool LinearProgrammingConstraint::BranchOnVar(IntegerVariable positive_var) {
420  // From the current LP solution, branch on the given var if fractional.
421  DCHECK(lp_solution_is_set_);
422  const double current_value = GetSolutionValue(positive_var);
423  DCHECK_GT(std::abs(current_value - std::round(current_value)), kCpEpsilon);
424 
425  // Used as empty reason in this method.
426  integer_reason_.clear();
427 
428  bool deductions_were_made = false;
429 
430  UpdateBoundsOfLpVariables();
431 
432  const IntegerValue current_obj_lb = integer_trail_->LowerBound(objective_cp_);
433  // This will try to branch in both direction around the LP value of the
434  // given variable and push any deduction done this way.
435 
436  const glop::ColIndex lp_var = GetOrCreateMirrorVariable(positive_var);
437  const double current_lb = ToDouble(integer_trail_->LowerBound(positive_var));
438  const double current_ub = ToDouble(integer_trail_->UpperBound(positive_var));
439  const double factor = scaler_.VariableScalingFactor(lp_var);
440  if (current_value < current_lb || current_value > current_ub) {
441  return false;
442  }
443 
444  // Form LP1 var <= floor(current_value)
445  const double new_ub = std::floor(current_value);
446  lp_data_.SetVariableBounds(lp_var, current_lb * factor, new_ub * factor);
447 
448  LPSolveInfo lower_branch_info = SolveLpForBranching();
449  if (lower_branch_info.status != glop::ProblemStatus::OPTIMAL &&
450  lower_branch_info.status != glop::ProblemStatus::DUAL_FEASIBLE &&
451  lower_branch_info.status != glop::ProblemStatus::DUAL_UNBOUNDED) {
452  return false;
453  }
454 
455  if (lower_branch_info.status == glop::ProblemStatus::DUAL_UNBOUNDED) {
456  // Push the other branch.
457  const IntegerLiteral deduction = IntegerLiteral::GreaterOrEqual(
458  positive_var, IntegerValue(std::ceil(current_value)));
459  if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) {
460  return false;
461  }
462  deductions_were_made = true;
463  } else if (lower_branch_info.new_obj_bound <= current_obj_lb) {
464  return false;
465  }
466 
467  // Form LP2 var >= ceil(current_value)
468  const double new_lb = std::ceil(current_value);
469  lp_data_.SetVariableBounds(lp_var, new_lb * factor, current_ub * factor);
470 
471  LPSolveInfo upper_branch_info = SolveLpForBranching();
472  if (upper_branch_info.status != glop::ProblemStatus::OPTIMAL &&
473  upper_branch_info.status != glop::ProblemStatus::DUAL_FEASIBLE &&
474  upper_branch_info.status != glop::ProblemStatus::DUAL_UNBOUNDED) {
475  return deductions_were_made;
476  }
477 
478  if (upper_branch_info.status == glop::ProblemStatus::DUAL_UNBOUNDED) {
479  // Push the other branch if not infeasible.
480  if (lower_branch_info.status != glop::ProblemStatus::DUAL_UNBOUNDED) {
481  const IntegerLiteral deduction = IntegerLiteral::LowerOrEqual(
482  positive_var, IntegerValue(std::floor(current_value)));
483  if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) {
484  return deductions_were_made;
485  }
486  deductions_were_made = true;
487  }
488  } else if (upper_branch_info.new_obj_bound <= current_obj_lb) {
489  return deductions_were_made;
490  }
491 
492  IntegerValue approximate_obj_lb = kMinIntegerValue;
493 
494  if (lower_branch_info.status == glop::ProblemStatus::DUAL_UNBOUNDED &&
495  upper_branch_info.status == glop::ProblemStatus::DUAL_UNBOUNDED) {
496  return integer_trail_->ReportConflict(integer_reason_);
497  } else if (lower_branch_info.status == glop::ProblemStatus::DUAL_UNBOUNDED) {
498  approximate_obj_lb = upper_branch_info.new_obj_bound;
499  } else if (upper_branch_info.status == glop::ProblemStatus::DUAL_UNBOUNDED) {
500  approximate_obj_lb = lower_branch_info.new_obj_bound;
501  } else {
502  approximate_obj_lb = std::min(lower_branch_info.new_obj_bound,
503  upper_branch_info.new_obj_bound);
504  }
505 
506  // NOTE: On some problems, the approximate_obj_lb could be inexact which add
507  // some tolerance to CP-SAT where currently there is none.
508  if (approximate_obj_lb <= current_obj_lb) return deductions_were_made;
509 
510  // Push the bound to the trail.
511  const IntegerLiteral deduction =
512  IntegerLiteral::GreaterOrEqual(objective_cp_, approximate_obj_lb);
513  if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) {
514  return deductions_were_made;
515  }
516 
517  return true;
518 }
519 
521  DCHECK(!lp_constraint_is_registered_);
522  lp_constraint_is_registered_ = true;
523  model->GetOrCreate<LinearProgrammingConstraintCollection>()->push_back(this);
524 
525  // Note fdid, this is not really needed by should lead to better cache
526  // locality.
527  std::sort(integer_objective_.begin(), integer_objective_.end());
528 
529  // Set the LP to its initial content.
530  if (!sat_parameters_.add_lp_constraints_lazily()) {
531  constraint_manager_.AddAllConstraintsToLp();
532  }
533  if (!CreateLpFromConstraintManager()) {
534  model->GetOrCreate<SatSolver>()->NotifyThatModelIsUnsat();
535  return;
536  }
537 
538  GenericLiteralWatcher* watcher = model->GetOrCreate<GenericLiteralWatcher>();
539  const int watcher_id = watcher->Register(this);
540  const int num_vars = integer_variables_.size();
541  for (int i = 0; i < num_vars; i++) {
542  watcher->WatchIntegerVariable(integer_variables_[i], watcher_id, i);
543  }
544  if (objective_is_defined_) {
545  watcher->WatchUpperBound(objective_cp_, watcher_id);
546  }
547  watcher->SetPropagatorPriority(watcher_id, 2);
548  watcher->AlwaysCallAtLevelZero(watcher_id);
549 
550  // Registering it with the trail make sure this class is always in sync when
551  // it is used in the decision heuristics.
552  integer_trail_->RegisterReversibleClass(this);
553  watcher->RegisterReversibleInt(watcher_id, &rev_optimal_constraints_size_);
554 }
555 
557  optimal_constraints_.resize(rev_optimal_constraints_size_);
558  if (lp_solution_is_set_ && level < lp_solution_level_) {
559  lp_solution_is_set_ = false;
560  }
561 
562  // Special case for level zero, we "reload" any previously known optimal
563  // solution from that level.
564  //
565  // TODO(user): Keep all optimal solution in the current branch?
566  // TODO(user): Still try to add cuts/constraints though!
567  if (level == 0 && !level_zero_lp_solution_.empty()) {
568  lp_solution_is_set_ = true;
569  lp_solution_ = level_zero_lp_solution_;
570  lp_solution_level_ = 0;
571  for (int i = 0; i < lp_solution_.size(); i++) {
572  expanded_lp_solution_[integer_variables_[i]] = lp_solution_[i];
573  expanded_lp_solution_[NegationOf(integer_variables_[i])] =
574  -lp_solution_[i];
575  }
576  }
577 }
578 
580  for (const IntegerVariable var : generator.vars) {
581  GetOrCreateMirrorVariable(VariableIsPositive(var) ? var : NegationOf(var));
582  }
583  cut_generators_.push_back(std::move(generator));
584 }
585 
587  const std::vector<int>& watch_indices) {
588  if (!lp_solution_is_set_) return Propagate();
589 
590  // At level zero, if there is still a chance to add cuts or lazy constraints,
591  // we re-run the LP.
592  if (trail_->CurrentDecisionLevel() == 0 && !lp_at_level_zero_is_final_) {
593  return Propagate();
594  }
595 
596  // Check whether the change breaks the current LP solution. If it does, call
597  // Propagate() on the current LP.
598  for (const int index : watch_indices) {
599  const double lb =
600  ToDouble(integer_trail_->LowerBound(integer_variables_[index]));
601  const double ub =
602  ToDouble(integer_trail_->UpperBound(integer_variables_[index]));
603  const double value = lp_solution_[index];
604  if (value < lb - kCpEpsilon || value > ub + kCpEpsilon) return Propagate();
605  }
606 
607  // TODO(user): The saved lp solution is still valid given the current variable
608  // bounds, so the LP optimal didn't change. However we might still want to add
609  // new cuts or new lazy constraints?
610  //
611  // TODO(user): Propagate the last optimal_constraint? Note that we need
612  // to be careful since the reversible int in IntegerSumLE are not registered.
613  // However, because we delete "optimalconstraints" on backtrack, we might not
614  // care.
615  return true;
616 }
617 
618 glop::Fractional LinearProgrammingConstraint::GetVariableValueAtCpScale(
619  glop::ColIndex var) {
620  return scaler_.UnscaleVariableValue(var, simplex_.GetVariableValue(var));
621 }
622 
624  IntegerVariable variable) const {
625  return lp_solution_[gtl::FindOrDie(mirror_lp_variable_, variable).value()];
626 }
627 
629  IntegerVariable variable) const {
630  return lp_reduced_cost_[gtl::FindOrDie(mirror_lp_variable_, variable)
631  .value()];
632 }
633 
634 void LinearProgrammingConstraint::UpdateBoundsOfLpVariables() {
635  const int num_vars = integer_variables_.size();
636  for (int i = 0; i < num_vars; i++) {
637  const IntegerVariable cp_var = integer_variables_[i];
638  const double lb = ToDouble(integer_trail_->LowerBound(cp_var));
639  const double ub = ToDouble(integer_trail_->UpperBound(cp_var));
640  const double factor = scaler_.VariableScalingFactor(glop::ColIndex(i));
641  lp_data_.SetVariableBounds(glop::ColIndex(i), lb * factor, ub * factor);
642  }
643 }
644 
645 bool LinearProgrammingConstraint::SolveLp() {
646  if (trail_->CurrentDecisionLevel() == 0) {
647  lp_at_level_zero_is_final_ = false;
648  }
649 
650  const auto status = simplex_.Solve(lp_data_, time_limit_);
651  total_num_simplex_iterations_ += simplex_.GetNumberOfIterations();
652  if (!status.ok()) {
653  VLOG(1) << "The LP solver encountered an error: " << status.error_message();
654  simplex_.ClearStateForNextSolve();
655  return false;
656  }
657  average_degeneracy_.AddData(CalculateDegeneracy());
658  if (average_degeneracy_.CurrentAverage() >= 1000.0) {
659  VLOG(2) << "High average degeneracy: "
660  << average_degeneracy_.CurrentAverage();
661  }
662 
663  const int status_as_int = static_cast<int>(simplex_.GetProblemStatus());
664  if (status_as_int >= num_solves_by_status_.size()) {
665  num_solves_by_status_.resize(status_as_int + 1);
666  }
667  num_solves_by_status_[status_as_int]++;
668  VLOG(2) << "lvl:" << trail_->CurrentDecisionLevel() << " "
669  << simplex_.GetProblemStatus()
670  << " iter:" << simplex_.GetNumberOfIterations()
671  << " obj:" << simplex_.GetObjectiveValue();
672 
674  lp_solution_is_set_ = true;
675  lp_solution_level_ = trail_->CurrentDecisionLevel();
676  const int num_vars = integer_variables_.size();
677  for (int i = 0; i < num_vars; i++) {
678  const glop::Fractional value =
679  GetVariableValueAtCpScale(glop::ColIndex(i));
680  lp_solution_[i] = value;
681  expanded_lp_solution_[integer_variables_[i]] = value;
682  expanded_lp_solution_[NegationOf(integer_variables_[i])] = -value;
683  }
684 
685  if (lp_solution_level_ == 0) {
686  level_zero_lp_solution_ = lp_solution_;
687  }
688  }
689  return true;
690 }
691 
692 bool LinearProgrammingConstraint::AddCutFromConstraints(
693  const std::string& name,
694  const std::vector<std::pair<RowIndex, IntegerValue>>& integer_multipliers) {
695  // This is initialized to a valid linear contraint (by taking linear
696  // combination of the LP rows) and will be transformed into a cut if
697  // possible.
698  //
699  // TODO(user): For CG cuts, Ideally this linear combination should have only
700  // one fractional variable (basis_col). But because of imprecision, we get a
701  // bunch of fractional entry with small coefficient (relative to the one of
702  // basis_col). We try to handle that in IntegerRoundingCut(), but it might be
703  // better to add small multiple of the involved rows to get rid of them.
704  IntegerValue cut_ub;
705  if (!ComputeNewLinearConstraint(integer_multipliers, &tmp_scattered_vector_,
706  &cut_ub)) {
707  VLOG(1) << "Issue, overflow!";
708  return false;
709  }
710 
711  // Important: because we use integer_multipliers below, we cannot just
712  // divide by GCD or call PreventOverflow() here.
713  //
714  // TODO(user): the conversion col_index -> IntegerVariable is slow and could
715  // in principle be removed. Easy for cuts, but not so much for
716  // implied_bounds_processor_. Note that in theory this could allow us to
717  // use Literal directly without the need to have an IntegerVariable for them.
718  tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, cut_ub,
719  &cut_);
720 
721  // Note that the base constraint we use are currently always tight.
722  // It is not a requirement though.
723  if (DEBUG_MODE) {
724  const double norm = ToDouble(ComputeInfinityNorm(cut_));
725  const double activity = ComputeActivity(cut_, expanded_lp_solution_);
726  if (std::abs(activity - ToDouble(cut_.ub)) / norm > 1e-4) {
727  VLOG(1) << "Cut not tight " << activity << " <= " << ToDouble(cut_.ub);
728  return false;
729  }
730  }
731  CHECK(constraint_manager_.DebugCheckConstraint(cut_));
732 
733  // We will create "artificial" variables after this index that will be
734  // substitued back into LP variables afterwards. Also not that we only use
735  // positive variable indices for these new variables, so that algorithm that
736  // take their negation will not mess up the indexing.
737  const IntegerVariable first_new_var(expanded_lp_solution_.size());
738  CHECK_EQ(first_new_var.value() % 2, 0);
739 
740  LinearConstraint copy_in_debug;
741  if (DEBUG_MODE) {
742  copy_in_debug = cut_;
743  }
744 
745  // Unlike for the knapsack cuts, it might not be always beneficial to
746  // process the implied bounds even though it seems to be better in average.
747  //
748  // TODO(user): Perform more experiments, in particular with which bound we use
749  // and if we complement or not before the MIR rounding. Other solvers seems
750  // to try different complementation strategies in a "potprocessing" and we
751  // don't. Try this too.
752  std::vector<ImpliedBoundsProcessor::SlackInfo> ib_slack_infos;
753  implied_bounds_processor_.ProcessUpperBoundedConstraintWithSlackCreation(
754  /*substitute_only_inner_variables=*/false, first_new_var,
755  expanded_lp_solution_, &cut_, &ib_slack_infos);
756  DCHECK(implied_bounds_processor_.DebugSlack(first_new_var, copy_in_debug,
757  cut_, ib_slack_infos));
758 
759  // Fills data for IntegerRoundingCut().
760  //
761  // Note(user): we use the current bound here, so the reasonement will only
762  // produce locally valid cut if we call this at a non-root node. We could
763  // use the level zero bounds if we wanted to generate a globally valid cut
764  // at another level. For now this is only called at level zero anyway.
765  tmp_lp_values_.clear();
766  tmp_var_lbs_.clear();
767  tmp_var_ubs_.clear();
768  for (const IntegerVariable var : cut_.vars) {
769  if (var >= first_new_var) {
771  const auto& info =
772  ib_slack_infos[(var.value() - first_new_var.value()) / 2];
773  tmp_lp_values_.push_back(info.lp_value);
774  tmp_var_lbs_.push_back(info.lb);
775  tmp_var_ubs_.push_back(info.ub);
776  } else {
777  tmp_lp_values_.push_back(expanded_lp_solution_[var]);
778  tmp_var_lbs_.push_back(integer_trail_->LevelZeroLowerBound(var));
779  tmp_var_ubs_.push_back(integer_trail_->LevelZeroUpperBound(var));
780  }
781  }
782 
783  // Add slack.
784  // definition: integer_lp_[row] + slack_row == bound;
785  const IntegerVariable first_slack(first_new_var +
786  IntegerVariable(2 * ib_slack_infos.size()));
787  tmp_slack_rows_.clear();
788  tmp_slack_bounds_.clear();
789  for (const auto pair : integer_multipliers) {
790  const RowIndex row = pair.first;
791  const IntegerValue coeff = pair.second;
792  const auto status = simplex_.GetConstraintStatus(row);
793  if (status == glop::ConstraintStatus::FIXED_VALUE) continue;
794 
795  tmp_lp_values_.push_back(0.0);
796  cut_.vars.push_back(first_slack +
797  2 * IntegerVariable(tmp_slack_rows_.size()));
798  tmp_slack_rows_.push_back(row);
799  cut_.coeffs.push_back(coeff);
800 
801  const IntegerValue diff(
802  CapSub(integer_lp_[row].ub.value(), integer_lp_[row].lb.value()));
803  if (coeff > 0) {
804  tmp_slack_bounds_.push_back(integer_lp_[row].ub);
805  tmp_var_lbs_.push_back(IntegerValue(0));
806  tmp_var_ubs_.push_back(diff);
807  } else {
808  tmp_slack_bounds_.push_back(integer_lp_[row].lb);
809  tmp_var_lbs_.push_back(-diff);
810  tmp_var_ubs_.push_back(IntegerValue(0));
811  }
812  }
813 
814  bool at_least_one_added = false;
815 
816  // Try cover appraoch to find cut.
817  {
818  if (cover_cut_helper_.TrySimpleKnapsack(cut_, tmp_lp_values_, tmp_var_lbs_,
819  tmp_var_ubs_)) {
820  at_least_one_added |= PostprocessAndAddCut(
821  absl::StrCat(name, "_K"), cover_cut_helper_.Info(), first_new_var,
822  first_slack, ib_slack_infos, cover_cut_helper_.mutable_cut());
823  }
824  }
825 
826  // Try integer rounding heuristic to find cut.
827  {
828  RoundingOptions options;
829  options.max_scaling = sat_parameters_.max_integer_rounding_scaling();
830  integer_rounding_cut_helper_.ComputeCut(options, tmp_lp_values_,
831  tmp_var_lbs_, tmp_var_ubs_,
832  &implied_bounds_processor_, &cut_);
833  at_least_one_added |= PostprocessAndAddCut(
834  name,
835  absl::StrCat("num_lifted_booleans=",
836  integer_rounding_cut_helper_.NumLiftedBooleans()),
837  first_new_var, first_slack, ib_slack_infos, &cut_);
838  }
839  return at_least_one_added;
840 }
841 
842 bool LinearProgrammingConstraint::PostprocessAndAddCut(
843  const std::string& name, const std::string& info,
844  IntegerVariable first_new_var, IntegerVariable first_slack,
845  const std::vector<ImpliedBoundsProcessor::SlackInfo>& ib_slack_infos,
846  LinearConstraint* cut) {
847  // Compute the activity. Warning: the cut no longer have the same size so we
848  // cannot use tmp_lp_values_. Note that the substitution below shouldn't
849  // change the activity by definition.
850  double activity = 0.0;
851  for (int i = 0; i < cut->vars.size(); ++i) {
852  if (cut->vars[i] < first_new_var) {
853  activity +=
854  ToDouble(cut->coeffs[i]) * expanded_lp_solution_[cut->vars[i]];
855  }
856  }
857  const double kMinViolation = 1e-4;
858  const double violation = activity - ToDouble(cut->ub);
859  if (violation < kMinViolation) {
860  VLOG(3) << "Bad cut " << activity << " <= " << ToDouble(cut->ub);
861  return false;
862  }
863 
864  // Substitute any slack left.
865  {
866  int num_slack = 0;
867  tmp_scattered_vector_.ClearAndResize(integer_variables_.size());
868  IntegerValue cut_ub = cut->ub;
869  bool overflow = false;
870  for (int i = 0; i < cut->vars.size(); ++i) {
871  const IntegerVariable var = cut->vars[i];
872 
873  // Simple copy for non-slack variables.
874  if (var < first_new_var) {
875  const glop::ColIndex col =
876  gtl::FindOrDie(mirror_lp_variable_, PositiveVariable(var));
877  if (VariableIsPositive(var)) {
878  tmp_scattered_vector_.Add(col, cut->coeffs[i]);
879  } else {
880  tmp_scattered_vector_.Add(col, -cut->coeffs[i]);
881  }
882  continue;
883  }
884 
885  // Replace slack from bound substitution.
886  if (var < first_slack) {
887  const IntegerValue multiplier = cut->coeffs[i];
888  const int index = (var.value() - first_new_var.value()) / 2;
889  CHECK_LT(index, ib_slack_infos.size());
890 
891  std::vector<std::pair<ColIndex, IntegerValue>> terms;
892  for (const std::pair<IntegerVariable, IntegerValue>& term :
893  ib_slack_infos[index].terms) {
894  terms.push_back(
895  {gtl::FindOrDie(mirror_lp_variable_,
896  PositiveVariable(term.first)),
897  VariableIsPositive(term.first) ? term.second : -term.second});
898  }
899  if (!tmp_scattered_vector_.AddLinearExpressionMultiple(multiplier,
900  terms)) {
901  overflow = true;
902  break;
903  }
904  if (!AddProductTo(multiplier, -ib_slack_infos[index].offset, &cut_ub)) {
905  overflow = true;
906  break;
907  }
908  continue;
909  }
910 
911  // Replace slack from LP constraints.
912  ++num_slack;
913  const int slack_index = (var.value() - first_slack.value()) / 2;
914  const glop::RowIndex row = tmp_slack_rows_[slack_index];
915  const IntegerValue multiplier = -cut->coeffs[i];
916  if (!tmp_scattered_vector_.AddLinearExpressionMultiple(
917  multiplier, integer_lp_[row].terms)) {
918  overflow = true;
919  break;
920  }
921 
922  // Update rhs.
923  if (!AddProductTo(multiplier, tmp_slack_bounds_[slack_index], &cut_ub)) {
924  overflow = true;
925  break;
926  }
927  }
928 
929  if (overflow) {
930  VLOG(1) << "Overflow in slack removal.";
931  return false;
932  }
933 
934  VLOG(3) << " num_slack: " << num_slack;
935  tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, cut_ub,
936  cut);
937  }
938 
939  // Display some stats used for investigation of cut generation.
940  const std::string extra_info =
941  absl::StrCat(info, " num_ib_substitutions=", ib_slack_infos.size());
942 
943  const double new_violation =
944  ComputeActivity(*cut, expanded_lp_solution_) - ToDouble(cut_.ub);
945  if (std::abs(violation - new_violation) >= 1e-4) {
946  VLOG(1) << "Violation discrepancy after slack removal. "
947  << " before = " << violation << " after = " << new_violation;
948  }
949 
950  DivideByGCD(cut);
951  return constraint_manager_.AddCut(*cut, name, expanded_lp_solution_,
952  extra_info);
953 }
954 
955 // TODO(user): This can be still too slow on some problems like
956 // 30_70_45_05_100.mps.gz. Not this actual function, but the set of computation
957 // it triggers. We should add heuristics to abort earlier if a cut is not
958 // promising. Or only test a few positions and not all rows.
959 void LinearProgrammingConstraint::AddCGCuts() {
960  const RowIndex num_rows = lp_data_.num_constraints();
961  std::vector<std::pair<RowIndex, double>> lp_multipliers;
962  std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
963  for (RowIndex row(0); row < num_rows; ++row) {
964  ColIndex basis_col = simplex_.GetBasis(row);
965  const Fractional lp_value = GetVariableValueAtCpScale(basis_col);
966 
967  // Only consider fractional basis element. We ignore element that are close
968  // to an integer to reduce the amount of positions we try.
969  //
970  // TODO(user): We could just look at the diff with std::floor() in the hope
971  // that when we are just under an integer, the exact computation below will
972  // also be just under it.
973  if (std::abs(lp_value - std::round(lp_value)) < 0.01) continue;
974 
975  // If this variable is a slack, we ignore it. This is because the
976  // corresponding row is not tight under the given lp values.
977  if (basis_col >= integer_variables_.size()) continue;
978 
979  if (time_limit_->LimitReached()) break;
980 
981  // TODO(user): Avoid code duplication between the sparse/dense path.
982  double magnitude = 0.0;
983  lp_multipliers.clear();
984  const glop::ScatteredRow& lambda = simplex_.GetUnitRowLeftInverse(row);
985  if (lambda.non_zeros.empty()) {
986  for (RowIndex row(0); row < num_rows; ++row) {
987  const double value = lambda.values[glop::RowToColIndex(row)];
988  if (std::abs(value) < kZeroTolerance) continue;
989 
990  // There should be no BASIC status, but they could be imprecision
991  // in the GetUnitRowLeftInverse() code? not sure, so better be safe.
992  const auto status = simplex_.GetConstraintStatus(row);
993  if (status == glop::ConstraintStatus::BASIC) {
994  VLOG(1) << "BASIC row not expected! " << value;
995  continue;
996  }
997 
998  magnitude = std::max(magnitude, std::abs(value));
999  lp_multipliers.push_back({row, value});
1000  }
1001  } else {
1002  for (const ColIndex col : lambda.non_zeros) {
1003  const RowIndex row = glop::ColToRowIndex(col);
1004  const double value = lambda.values[col];
1005  if (std::abs(value) < kZeroTolerance) continue;
1006 
1007  const auto status = simplex_.GetConstraintStatus(row);
1008  if (status == glop::ConstraintStatus::BASIC) {
1009  VLOG(1) << "BASIC row not expected! " << value;
1010  continue;
1011  }
1012 
1013  magnitude = std::max(magnitude, std::abs(value));
1014  lp_multipliers.push_back({row, value});
1015  }
1016  }
1017  if (lp_multipliers.empty()) continue;
1018 
1019  Fractional scaling;
1020  for (int i = 0; i < 2; ++i) {
1021  if (i == 1) {
1022  // Try other sign.
1023  //
1024  // TODO(user): Maybe add an heuristic to know beforehand which sign to
1025  // use?
1026  for (std::pair<RowIndex, double>& p : lp_multipliers) {
1027  p.second = -p.second;
1028  }
1029  }
1030 
1031  // TODO(user): We use a lower value here otherwise we might run into
1032  // overflow while computing the cut. This should be fixable.
1033  integer_multipliers =
1034  ScaleLpMultiplier(/*take_objective_into_account=*/false,
1035  lp_multipliers, &scaling, /*max_pow=*/52);
1036  AddCutFromConstraints("CG", integer_multipliers);
1037  }
1038  }
1039 }
1040 
1041 namespace {
1042 
1043 // For each element of a, adds a random one in b and append the pair to output.
1044 void RandomPick(const std::vector<RowIndex>& a, const std::vector<RowIndex>& b,
1045  ModelRandomGenerator* random,
1046  std::vector<std::pair<RowIndex, RowIndex>>* output) {
1047  if (a.empty() || b.empty()) return;
1048  for (const RowIndex row : a) {
1049  const RowIndex other = b[absl::Uniform<int>(*random, 0, b.size())];
1050  if (other != row) {
1051  output->push_back({row, other});
1052  }
1053  }
1054 }
1055 
1056 template <class ListOfTerms>
1057 IntegerValue GetCoeff(ColIndex col, const ListOfTerms& terms) {
1058  for (const auto& term : terms) {
1059  if (term.first == col) return term.second;
1060  }
1061  return IntegerValue(0);
1062 }
1063 
1064 } // namespace
1065 
1066 void LinearProgrammingConstraint::AddMirCuts() {
1067  // Heuristic to generate MIR_n cuts by combining a small number of rows. This
1068  // works greedily and follow more or less the MIR cut description in the
1069  // literature. We have a current cut, and we add one more row to it while
1070  // eliminating a variable of the current cut whose LP value is far from its
1071  // bound.
1072  //
1073  // A notable difference is that we randomize the variable we eliminate and
1074  // the row we use to do so. We still have weights to indicate our preferred
1075  // choices. This allows to generate different cuts when called again and
1076  // again.
1077  //
1078  // TODO(user): We could combine n rows to make sure we eliminate n variables
1079  // far away from their bounds by solving exactly in integer small linear
1080  // system.
1082  integer_variables_.size(), IntegerValue(0));
1083  SparseBitset<ColIndex> non_zeros_(ColIndex(integer_variables_.size()));
1084 
1085  // We compute all the rows that are tight, these will be used as the base row
1086  // for the MIR_n procedure below.
1087  const RowIndex num_rows = lp_data_.num_constraints();
1088  std::vector<std::pair<RowIndex, IntegerValue>> base_rows;
1089  absl::StrongVector<RowIndex, double> row_weights(num_rows.value(), 0.0);
1090  for (RowIndex row(0); row < num_rows; ++row) {
1091  const auto status = simplex_.GetConstraintStatus(row);
1092  if (status == glop::ConstraintStatus::BASIC) continue;
1093  if (status == glop::ConstraintStatus::FREE) continue;
1094 
1097  base_rows.push_back({row, IntegerValue(1)});
1098  }
1101  base_rows.push_back({row, IntegerValue(-1)});
1102  }
1103 
1104  // For now, we use the dual values for the row "weights".
1105  //
1106  // Note that we use the dual at LP scale so that it make more sense when we
1107  // compare different rows since the LP has been scaled.
1108  //
1109  // TODO(user): In Kati Wolter PhD "Implementation of Cutting Plane
1110  // Separators for Mixed Integer Programs" which describe SCIP's MIR cuts
1111  // implementation (or at least an early version of it), a more complex score
1112  // is used.
1113  //
1114  // Note(user): Because we only consider tight rows under the current lp
1115  // solution (i.e. non-basic rows), most should have a non-zero dual values.
1116  // But there is some degenerate problem where these rows have a really low
1117  // weight (or even zero), and having only weight of exactly zero in
1118  // std::discrete_distribution will result in a crash.
1119  row_weights[row] = std::max(1e-8, std::abs(simplex_.GetDualValue(row)));
1120  }
1121 
1122  std::vector<double> weights;
1124  std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
1125  for (const std::pair<RowIndex, IntegerValue>& entry : base_rows) {
1126  if (time_limit_->LimitReached()) break;
1127 
1128  // First try to generate a cut directly from this base row (MIR1).
1129  //
1130  // Note(user): We abort on success like it seems to be done in the
1131  // literature. Note that we don't succeed that often in generating an
1132  // efficient cut, so I am not sure aborting will make a big difference
1133  // speedwise. We might generate similar cuts though, but hopefully the cut
1134  // management can deal with that.
1135  integer_multipliers = {entry};
1136  if (AddCutFromConstraints("MIR_1", integer_multipliers)) {
1137  continue;
1138  }
1139 
1140  // Cleanup.
1141  for (const ColIndex col : non_zeros_.PositionsSetAtLeastOnce()) {
1142  dense_cut[col] = IntegerValue(0);
1143  }
1144  non_zeros_.SparseClearAll();
1145 
1146  // Copy cut.
1147  const IntegerValue multiplier = entry.second;
1148  for (const std::pair<ColIndex, IntegerValue> term :
1149  integer_lp_[entry.first].terms) {
1150  const ColIndex col = term.first;
1151  const IntegerValue coeff = term.second;
1152  non_zeros_.Set(col);
1153  dense_cut[col] += coeff * multiplier;
1154  }
1155 
1156  used_rows.assign(num_rows.value(), false);
1157  used_rows[entry.first] = true;
1158 
1159  // We will aggregate at most kMaxAggregation more rows.
1160  //
1161  // TODO(user): optim + tune.
1162  const int kMaxAggregation = 5;
1163  for (int i = 0; i < kMaxAggregation; ++i) {
1164  // First pick a variable to eliminate. We currently pick a random one with
1165  // a weight that depend on how far it is from its closest bound.
1166  IntegerValue max_magnitude(0);
1167  weights.clear();
1168  std::vector<ColIndex> col_candidates;
1169  for (const ColIndex col : non_zeros_.PositionsSetAtLeastOnce()) {
1170  if (dense_cut[col] == 0) continue;
1171 
1172  max_magnitude = std::max(max_magnitude, IntTypeAbs(dense_cut[col]));
1173  const int col_degree =
1174  lp_data_.GetSparseColumn(col).num_entries().value();
1175  if (col_degree <= 1) continue;
1177  continue;
1178  }
1179 
1180  const IntegerVariable var = integer_variables_[col.value()];
1181  const double lp_value = expanded_lp_solution_[var];
1182  const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(var));
1183  const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(var));
1184  const double bound_distance = std::min(ub - lp_value, lp_value - lb);
1185  if (bound_distance > 1e-2) {
1186  weights.push_back(bound_distance);
1187  col_candidates.push_back(col);
1188  }
1189  }
1190  if (col_candidates.empty()) break;
1191 
1192  const ColIndex var_to_eliminate =
1193  col_candidates[std::discrete_distribution<>(weights.begin(),
1194  weights.end())(*random_)];
1195 
1196  // What rows can we add to eliminate var_to_eliminate?
1197  std::vector<RowIndex> possible_rows;
1198  weights.clear();
1199  for (const auto entry : lp_data_.GetSparseColumn(var_to_eliminate)) {
1200  const RowIndex row = entry.row();
1201  const auto status = simplex_.GetConstraintStatus(row);
1202  if (status == glop::ConstraintStatus::BASIC) continue;
1203  if (status == glop::ConstraintStatus::FREE) continue;
1204 
1205  // We disallow all the rows that contain a variable that we already
1206  // eliminated (or are about to). This mean that we choose rows that
1207  // form a "triangular" matrix on the position we choose to eliminate.
1208  if (used_rows[row]) continue;
1209  used_rows[row] = true;
1210 
1211  // TODO(user): Instead of using FIXED_VALUE consider also both direction
1212  // when we almost have an equality? that is if the LP constraints bounds
1213  // are close from each others (<1e-6 ?). Initial experiments shows it
1214  // doesn't change much, so I kept this version for now. Note that it
1215  // might just be better to use the side that constrain the current lp
1216  // optimal solution (that we get from the status).
1217  bool add_row = false;
1218  if (status == glop::ConstraintStatus::FIXED_VALUE ||
1220  if (entry.coefficient() > 0.0) {
1221  if (dense_cut[var_to_eliminate] < 0) add_row = true;
1222  } else {
1223  if (dense_cut[var_to_eliminate] > 0) add_row = true;
1224  }
1225  }
1226  if (status == glop::ConstraintStatus::FIXED_VALUE ||
1228  if (entry.coefficient() > 0.0) {
1229  if (dense_cut[var_to_eliminate] > 0) add_row = true;
1230  } else {
1231  if (dense_cut[var_to_eliminate] < 0) add_row = true;
1232  }
1233  }
1234  if (add_row) {
1235  possible_rows.push_back(row);
1236  weights.push_back(row_weights[row]);
1237  }
1238  }
1239  if (possible_rows.empty()) break;
1240 
1241  const RowIndex row_to_combine =
1242  possible_rows[std::discrete_distribution<>(weights.begin(),
1243  weights.end())(*random_)];
1244  const IntegerValue to_combine_coeff =
1245  GetCoeff(var_to_eliminate, integer_lp_[row_to_combine].terms);
1246  CHECK_NE(to_combine_coeff, 0);
1247 
1248  IntegerValue mult1 = -to_combine_coeff;
1249  IntegerValue mult2 = dense_cut[var_to_eliminate];
1250  CHECK_NE(mult2, 0);
1251  if (mult1 < 0) {
1252  mult1 = -mult1;
1253  mult2 = -mult2;
1254  }
1255 
1256  const IntegerValue gcd = IntegerValue(
1257  MathUtil::GCD64(std::abs(mult1.value()), std::abs(mult2.value())));
1258  CHECK_NE(gcd, 0);
1259  mult1 /= gcd;
1260  mult2 /= gcd;
1261 
1262  // Overflow detection.
1263  //
1264  // TODO(user): do that in the possible_rows selection? only problem is
1265  // that we do not have the integer coefficient there...
1266  for (std::pair<RowIndex, IntegerValue>& entry : integer_multipliers) {
1267  max_magnitude = std::max(max_magnitude, entry.second);
1268  }
1269  if (CapAdd(CapProd(max_magnitude.value(), std::abs(mult1.value())),
1270  CapProd(infinity_norms_[row_to_combine].value(),
1271  std::abs(mult2.value()))) == kint64max) {
1272  break;
1273  }
1274 
1275  for (std::pair<RowIndex, IntegerValue>& entry : integer_multipliers) {
1276  entry.second *= mult1;
1277  }
1278  integer_multipliers.push_back({row_to_combine, mult2});
1279 
1280  // TODO(user): Not supper efficient to recombine the rows.
1281  if (AddCutFromConstraints(absl::StrCat("MIR_", i + 2),
1282  integer_multipliers)) {
1283  break;
1284  }
1285 
1286  // Minor optim: the computation below is only needed if we do one more
1287  // iteration.
1288  if (i + 1 == kMaxAggregation) break;
1289 
1290  for (ColIndex col : non_zeros_.PositionsSetAtLeastOnce()) {
1291  dense_cut[col] *= mult1;
1292  }
1293  for (const std::pair<ColIndex, IntegerValue> term :
1294  integer_lp_[row_to_combine].terms) {
1295  const ColIndex col = term.first;
1296  const IntegerValue coeff = term.second;
1297  non_zeros_.Set(col);
1298  dense_cut[col] += coeff * mult2;
1299  }
1300  }
1301  }
1302 }
1303 
1304 void LinearProgrammingConstraint::AddZeroHalfCuts() {
1305  if (time_limit_->LimitReached()) return;
1306 
1307  tmp_lp_values_.clear();
1308  tmp_var_lbs_.clear();
1309  tmp_var_ubs_.clear();
1310  for (const IntegerVariable var : integer_variables_) {
1311  tmp_lp_values_.push_back(expanded_lp_solution_[var]);
1312  tmp_var_lbs_.push_back(integer_trail_->LevelZeroLowerBound(var));
1313  tmp_var_ubs_.push_back(integer_trail_->LevelZeroUpperBound(var));
1314  }
1315 
1316  // TODO(user): See if it make sense to try to use implied bounds there.
1317  zero_half_cut_helper_.ProcessVariables(tmp_lp_values_, tmp_var_lbs_,
1318  tmp_var_ubs_);
1319  for (glop::RowIndex row(0); row < integer_lp_.size(); ++row) {
1320  // Even though we could use non-tight row, for now we prefer to use tight
1321  // ones.
1322  const auto status = simplex_.GetConstraintStatus(row);
1323  if (status == glop::ConstraintStatus::BASIC) continue;
1324  if (status == glop::ConstraintStatus::FREE) continue;
1325 
1326  zero_half_cut_helper_.AddOneConstraint(
1327  row, integer_lp_[row].terms, integer_lp_[row].lb, integer_lp_[row].ub);
1328  }
1329  for (const std::vector<std::pair<RowIndex, IntegerValue>>& multipliers :
1330  zero_half_cut_helper_.InterestingCandidates(random_)) {
1331  if (time_limit_->LimitReached()) break;
1332 
1333  // TODO(user): Make sure that if the resulting linear coefficients are not
1334  // too high, we do try a "divisor" of two and thus try a true zero-half cut
1335  // instead of just using our best MIR heuristic (which might still be better
1336  // though).
1337  AddCutFromConstraints("ZERO_HALF", multipliers);
1338  }
1339 }
1340 
1341 void LinearProgrammingConstraint::UpdateSimplexIterationLimit(
1342  const int64 min_iter, const int64 max_iter) {
1343  if (sat_parameters_.linearization_level() < 2) return;
1344  const int64 num_degenerate_columns = CalculateDegeneracy();
1345  const int64 num_cols = simplex_.GetProblemNumCols().value();
1346  if (num_cols <= 0) {
1347  return;
1348  }
1349  CHECK_GT(num_cols, 0);
1350  const int64 decrease_factor = (10 * num_degenerate_columns) / num_cols;
1352  // We reached here probably because we predicted wrong. We use this as a
1353  // signal to increase the iterations or punish less for degeneracy compare
1354  // to the other part.
1355  if (is_degenerate_) {
1356  next_simplex_iter_ /= std::max(int64{1}, decrease_factor);
1357  } else {
1358  next_simplex_iter_ *= 2;
1359  }
1360  } else if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) {
1361  if (is_degenerate_) {
1362  next_simplex_iter_ /= std::max(int64{1}, 2 * decrease_factor);
1363  } else {
1364  // This is the most common case. We use the size of the problem to
1365  // determine the limit and ignore the previous limit.
1366  next_simplex_iter_ = num_cols / 40;
1367  }
1368  }
1369  next_simplex_iter_ =
1370  std::max(min_iter, std::min(max_iter, next_simplex_iter_));
1371 }
1372 
1374  UpdateBoundsOfLpVariables();
1375 
1376  // TODO(user): It seems the time we loose by not stopping early might be worth
1377  // it because we end up with a better explanation at optimality.
1378  glop::GlopParameters parameters = simplex_.GetParameters();
1379  if (/* DISABLES CODE */ (false) && objective_is_defined_) {
1380  // We put a limit on the dual objective since there is no point increasing
1381  // it past our current objective upper-bound (we will already fail as soon
1382  // as we pass it). Note that this limit is properly transformed using the
1383  // objective scaling factor and offset stored in lp_data_.
1384  //
1385  // Note that we use a bigger epsilon here to be sure that if we abort
1386  // because of this, we will report a conflict.
1387  parameters.set_objective_upper_limit(
1388  static_cast<double>(integer_trail_->UpperBound(objective_cp_).value() +
1389  100.0 * kCpEpsilon));
1390  }
1391 
1392  // Put an iteration limit on the work we do in the simplex for this call. Note
1393  // that because we are "incremental", even if we don't solve it this time we
1394  // will make progress towards a solve in the lower node of the tree search.
1395  if (trail_->CurrentDecisionLevel() == 0) {
1396  // TODO(user): Dynamically change the iteration limit for root node as
1397  // well.
1398  parameters.set_max_number_of_iterations(2000);
1399  } else {
1400  parameters.set_max_number_of_iterations(next_simplex_iter_);
1401  }
1402  if (sat_parameters_.use_exact_lp_reason()) {
1403  parameters.set_change_status_to_imprecise(false);
1404  parameters.set_primal_feasibility_tolerance(1e-7);
1405  parameters.set_dual_feasibility_tolerance(1e-7);
1406  }
1407 
1408  simplex_.SetParameters(parameters);
1410  if (!SolveLp()) return true;
1411 
1412  // Add new constraints to the LP and resolve?
1413  const int max_cuts_rounds =
1414  trail_->CurrentDecisionLevel() == 0
1415  ? sat_parameters_.max_cut_rounds_at_level_zero()
1416  : 1;
1417  int cuts_round = 0;
1418  while (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL &&
1419  cuts_round < max_cuts_rounds) {
1420  // We wait for the first batch of problem constraints to be added before we
1421  // begin to generate cuts.
1422  cuts_round++;
1423  if (!integer_lp_.empty()) {
1424  implied_bounds_processor_.ClearCache();
1425  implied_bounds_processor_.SeparateSomeImpliedBoundCuts(
1426  expanded_lp_solution_);
1427 
1428  // The "generic" cuts are currently part of this class as they are using
1429  // data from the current LP.
1430  //
1431  // TODO(user): Refactor so that they are just normal cut generators?
1432  if (trail_->CurrentDecisionLevel() == 0) {
1433  if (sat_parameters_.add_mir_cuts()) AddMirCuts();
1434  if (sat_parameters_.add_cg_cuts()) AddCGCuts();
1435  if (sat_parameters_.add_zero_half_cuts()) AddZeroHalfCuts();
1436  }
1437 
1438  // Try to add cuts.
1439  if (!cut_generators_.empty() &&
1440  (trail_->CurrentDecisionLevel() == 0 ||
1441  !sat_parameters_.only_add_cuts_at_level_zero())) {
1442  for (const CutGenerator& generator : cut_generators_) {
1443  generator.generate_cuts(expanded_lp_solution_, &constraint_manager_);
1444  }
1445  }
1446 
1447  implied_bounds_processor_.IbCutPool().TransferToManager(
1448  expanded_lp_solution_, &constraint_manager_);
1449  }
1450 
1451  glop::BasisState state = simplex_.GetState();
1452  if (constraint_manager_.ChangeLp(expanded_lp_solution_, &state)) {
1453  simplex_.LoadStateForNextSolve(state);
1454  if (!CreateLpFromConstraintManager()) {
1455  return integer_trail_->ReportConflict({});
1456  }
1457  const double old_obj = simplex_.GetObjectiveValue();
1458  if (!SolveLp()) return true;
1459  if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) {
1460  VLOG(1) << "Relaxation improvement " << old_obj << " -> "
1461  << simplex_.GetObjectiveValue()
1462  << " diff: " << simplex_.GetObjectiveValue() - old_obj
1463  << " level: " << trail_->CurrentDecisionLevel();
1464  }
1465  } else {
1466  if (trail_->CurrentDecisionLevel() == 0) {
1467  lp_at_level_zero_is_final_ = true;
1468  }
1469  break;
1470  }
1471  }
1472 
1473  // A dual-unbounded problem is infeasible. We use the dual ray reason.
1475  if (sat_parameters_.use_exact_lp_reason()) {
1476  if (!FillExactDualRayReason()) return true;
1477  } else {
1478  FillReducedCostReasonIn(simplex_.GetDualRayRowCombination(),
1479  &integer_reason_);
1480  }
1481  return integer_trail_->ReportConflict(integer_reason_);
1482  }
1483 
1484  // TODO(user): Update limits for DUAL_UNBOUNDED status as well.
1485  UpdateSimplexIterationLimit(/*min_iter=*/10, /*max_iter=*/1000);
1486 
1487  // Optimality deductions if problem has an objective.
1488  if (objective_is_defined_ &&
1491  // Try to filter optimal objective value. Note that GetObjectiveValue()
1492  // already take care of the scaling so that it returns an objective in the
1493  // CP world.
1494  const double relaxed_optimal_objective = simplex_.GetObjectiveValue();
1495  const IntegerValue approximate_new_lb(
1496  static_cast<int64>(std::ceil(relaxed_optimal_objective - kCpEpsilon)));
1497 
1498  // TODO(user): Maybe do a bit less computation when we cannot propagate
1499  // anything.
1500  if (sat_parameters_.use_exact_lp_reason()) {
1501  if (!ExactLpReasonning()) return false;
1502 
1503  // Display when the inexact bound would have propagated more.
1504  const IntegerValue propagated_lb =
1505  integer_trail_->LowerBound(objective_cp_);
1506  if (approximate_new_lb > propagated_lb) {
1507  VLOG(2) << "LP objective [ " << ToDouble(propagated_lb) << ", "
1508  << ToDouble(integer_trail_->UpperBound(objective_cp_))
1509  << " ] approx_lb += "
1510  << ToDouble(approximate_new_lb - propagated_lb) << " gap: "
1511  << integer_trail_->UpperBound(objective_cp_) - propagated_lb;
1512  }
1513  } else {
1514  FillReducedCostReasonIn(simplex_.GetReducedCosts(), &integer_reason_);
1515  const double objective_cp_ub =
1516  ToDouble(integer_trail_->UpperBound(objective_cp_));
1517  ReducedCostStrengtheningDeductions(objective_cp_ub -
1518  relaxed_optimal_objective);
1519  if (!deductions_.empty()) {
1520  deductions_reason_ = integer_reason_;
1521  deductions_reason_.push_back(
1522  integer_trail_->UpperBoundAsLiteral(objective_cp_));
1523  }
1524 
1525  // Push new objective lb.
1526  if (approximate_new_lb > integer_trail_->LowerBound(objective_cp_)) {
1527  const IntegerLiteral deduction =
1528  IntegerLiteral::GreaterOrEqual(objective_cp_, approximate_new_lb);
1529  if (!integer_trail_->Enqueue(deduction, {}, integer_reason_)) {
1530  return false;
1531  }
1532  }
1533 
1534  // Push reduced cost strengthening bounds.
1535  if (!deductions_.empty()) {
1536  const int trail_index_with_same_reason = integer_trail_->Index();
1537  for (const IntegerLiteral deduction : deductions_) {
1538  if (!integer_trail_->Enqueue(deduction, {}, deductions_reason_,
1539  trail_index_with_same_reason)) {
1540  return false;
1541  }
1542  }
1543  }
1544  }
1545  }
1546 
1547  // Copy more info about the current solution.
1548  if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL) {
1549  CHECK(lp_solution_is_set_);
1550 
1551  lp_objective_ = simplex_.GetObjectiveValue();
1552  lp_solution_is_integer_ = true;
1553  const int num_vars = integer_variables_.size();
1554  for (int i = 0; i < num_vars; i++) {
1555  lp_reduced_cost_[i] = scaler_.UnscaleReducedCost(
1556  glop::ColIndex(i), simplex_.GetReducedCost(glop::ColIndex(i)));
1557  if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) >
1558  kCpEpsilon) {
1559  lp_solution_is_integer_ = false;
1560  }
1561  }
1562 
1563  if (compute_reduced_cost_averages_) {
1564  UpdateAverageReducedCosts();
1565  }
1566  }
1567 
1568  if (sat_parameters_.use_branching_in_lp() && objective_is_defined_ &&
1569  trail_->CurrentDecisionLevel() == 0 && !is_degenerate_ &&
1570  lp_solution_is_set_ && !lp_solution_is_integer_ &&
1571  sat_parameters_.linearization_level() >= 2 &&
1572  compute_reduced_cost_averages_ &&
1574  count_since_last_branching_++;
1575  if (count_since_last_branching_ < branching_frequency_) {
1576  return true;
1577  }
1578  count_since_last_branching_ = 0;
1579  bool branching_successful = false;
1580 
1581  // Strong branching on top max_num_branches variable.
1582  const int max_num_branches = 3;
1583  const int num_vars = integer_variables_.size();
1584  std::vector<std::pair<double, IntegerVariable>> branching_vars;
1585  for (int i = 0; i < num_vars; ++i) {
1586  const IntegerVariable var = integer_variables_[i];
1587  const IntegerVariable positive_var = PositiveVariable(var);
1588 
1589  // Skip non fractional variables.
1590  const double current_value = GetSolutionValue(positive_var);
1591  if (std::abs(current_value - std::round(current_value)) <= kCpEpsilon) {
1592  continue;
1593  }
1594 
1595  // Skip ignored variables.
1596  if (integer_trail_->IsCurrentlyIgnored(var)) continue;
1597 
1598  // We can use any metric to select a variable to branch on. Reduced cost
1599  // average is one of the most promissing metric. It captures the history
1600  // of the objective bound improvement in LP due to changes in the given
1601  // variable bounds.
1602  //
1603  // NOTE: We also experimented using PseudoCosts and most recent reduced
1604  // cost as metrics but it doesn't give better results on benchmarks.
1605  const double cost_i = rc_scores_[i];
1606  std::pair<double, IntegerVariable> branching_var =
1607  std::make_pair(-cost_i, positive_var);
1608  auto iterator = std::lower_bound(branching_vars.begin(),
1609  branching_vars.end(), branching_var);
1610 
1611  branching_vars.insert(iterator, branching_var);
1612  if (branching_vars.size() > max_num_branches) {
1613  branching_vars.resize(max_num_branches);
1614  }
1615  }
1616 
1617  for (const std::pair<double, IntegerVariable>& branching_var :
1618  branching_vars) {
1619  const IntegerVariable positive_var = branching_var.second;
1620  VLOG(2) << "Branching on: " << positive_var;
1621  if (BranchOnVar(positive_var)) {
1622  VLOG(2) << "Branching successful.";
1623  branching_successful = true;
1624  } else {
1625  break;
1626  }
1627  }
1628  if (!branching_successful) {
1629  branching_frequency_ *= 2;
1630  }
1631  }
1632  return true;
1633 }
1634 
1635 // Returns kMinIntegerValue in case of overflow.
1636 //
1637 // TODO(user): Because of PreventOverflow(), this should actually never happen.
1638 IntegerValue LinearProgrammingConstraint::GetImpliedLowerBound(
1639  const LinearConstraint& terms) const {
1640  IntegerValue lower_bound(0);
1641  const int size = terms.vars.size();
1642  for (int i = 0; i < size; ++i) {
1643  const IntegerVariable var = terms.vars[i];
1644  const IntegerValue coeff = terms.coeffs[i];
1645  CHECK_NE(coeff, 0);
1646  const IntegerValue bound = coeff > 0 ? integer_trail_->LowerBound(var)
1647  : integer_trail_->UpperBound(var);
1648  if (!AddProductTo(bound, coeff, &lower_bound)) return kMinIntegerValue;
1649  }
1650  return lower_bound;
1651 }
1652 
1653 bool LinearProgrammingConstraint::PossibleOverflow(
1654  const LinearConstraint& constraint) {
1655  IntegerValue lower_bound(0);
1656  const int size = constraint.vars.size();
1657  for (int i = 0; i < size; ++i) {
1658  const IntegerVariable var = constraint.vars[i];
1659  const IntegerValue coeff = constraint.coeffs[i];
1660  CHECK_NE(coeff, 0);
1661  const IntegerValue bound = coeff > 0
1662  ? integer_trail_->LevelZeroLowerBound(var)
1663  : integer_trail_->LevelZeroUpperBound(var);
1664  if (!AddProductTo(bound, coeff, &lower_bound)) {
1665  return true;
1666  }
1667  }
1668  const int64 slack = CapAdd(lower_bound.value(), -constraint.ub.value());
1669  if (slack == kint64min || slack == kint64max) {
1670  return true;
1671  }
1672  return false;
1673 }
1674 
1675 namespace {
1676 
1677 absl::int128 FloorRatio128(absl::int128 x, IntegerValue positive_div) {
1678  absl::int128 div128(positive_div.value());
1679  absl::int128 result = x / div128;
1680  if (result * div128 > x) return result - 1;
1681  return result;
1682 }
1683 
1684 } // namespace
1685 
1686 void LinearProgrammingConstraint::PreventOverflow(LinearConstraint* constraint,
1687  int max_pow) {
1688  // Compute the min/max possible partial sum. Note that we need to use the
1689  // level zero bounds here since we might use this cut after backtrack.
1690  double sum_min = std::min(0.0, ToDouble(-constraint->ub));
1691  double sum_max = std::max(0.0, ToDouble(-constraint->ub));
1692  const int size = constraint->vars.size();
1693  for (int i = 0; i < size; ++i) {
1694  const IntegerVariable var = constraint->vars[i];
1695  const double coeff = ToDouble(constraint->coeffs[i]);
1696  const double prod1 =
1697  coeff * ToDouble(integer_trail_->LevelZeroLowerBound(var));
1698  const double prod2 =
1699  coeff * ToDouble(integer_trail_->LevelZeroUpperBound(var));
1700  sum_min += std::min(0.0, std::min(prod1, prod2));
1701  sum_max += std::max(0.0, std::max(prod1, prod2));
1702  }
1703  const double max_value = std::max({sum_max, -sum_min, sum_max - sum_min});
1704 
1705  const IntegerValue divisor(std::ceil(std::ldexp(max_value, -max_pow)));
1706  if (divisor <= 1) return;
1707 
1708  // To be correct, we need to shift all variable so that they are positive.
1709  //
1710  // Important: One might be tempted to think that using the current variable
1711  // bounds is okay here since we only use this to derive cut/constraint that
1712  // only needs to be locally valid. However, in some corner cases (like when
1713  // one term become zero), we might loose the fact that we used one of the
1714  // variable bound to derive the new constraint, so we will miss it in the
1715  // explanation !!
1716  //
1717  // TODO(user): This code is tricky and similar to the one to generate cuts.
1718  // Test and may reduce the duplication? note however that here we use int128
1719  // to deal with potential overflow.
1720  int new_size = 0;
1721  absl::int128 adjust = 0;
1722  for (int i = 0; i < size; ++i) {
1723  const IntegerValue old_coeff = constraint->coeffs[i];
1724  const IntegerValue new_coeff = FloorRatio(old_coeff, divisor);
1725 
1726  // Compute the rhs adjustement.
1727  const absl::int128 remainder =
1728  absl::int128(old_coeff.value()) -
1729  absl::int128(new_coeff.value()) * absl::int128(divisor.value());
1730  adjust +=
1731  remainder *
1732  absl::int128(
1733  integer_trail_->LevelZeroLowerBound(constraint->vars[i]).value());
1734 
1735  if (new_coeff == 0) continue;
1736  constraint->vars[new_size] = constraint->vars[i];
1737  constraint->coeffs[new_size] = new_coeff;
1738  ++new_size;
1739  }
1740  constraint->vars.resize(new_size);
1741  constraint->coeffs.resize(new_size);
1742 
1743  constraint->ub = IntegerValue(static_cast<int64>(
1744  FloorRatio128(absl::int128(constraint->ub.value()) - adjust, divisor)));
1745 }
1746 
1747 // TODO(user): combine this with RelaxLinearReason() to avoid the extra
1748 // magnitude vector and the weird precondition of RelaxLinearReason().
1749 void LinearProgrammingConstraint::SetImpliedLowerBoundReason(
1750  const LinearConstraint& terms, IntegerValue slack) {
1751  integer_reason_.clear();
1752  std::vector<IntegerValue> magnitudes;
1753  const int size = terms.vars.size();
1754  for (int i = 0; i < size; ++i) {
1755  const IntegerVariable var = terms.vars[i];
1756  const IntegerValue coeff = terms.coeffs[i];
1757  CHECK_NE(coeff, 0);
1758  if (coeff > 0) {
1759  magnitudes.push_back(coeff);
1760  integer_reason_.push_back(integer_trail_->LowerBoundAsLiteral(var));
1761  } else {
1762  magnitudes.push_back(-coeff);
1763  integer_reason_.push_back(integer_trail_->UpperBoundAsLiteral(var));
1764  }
1765  }
1766  CHECK_GE(slack, 0);
1767  if (slack > 0) {
1768  integer_trail_->RelaxLinearReason(slack, magnitudes, &integer_reason_);
1769  }
1770  integer_trail_->RemoveLevelZeroBounds(&integer_reason_);
1771 }
1772 
1773 std::vector<std::pair<RowIndex, IntegerValue>>
1774 LinearProgrammingConstraint::ScaleLpMultiplier(
1775  bool take_objective_into_account,
1776  const std::vector<std::pair<RowIndex, double>>& lp_multipliers,
1777  Fractional* scaling, int max_pow) const {
1778  double max_sum = 0.0;
1779  tmp_cp_multipliers_.clear();
1780  for (const std::pair<RowIndex, double>& p : lp_multipliers) {
1781  const RowIndex row = p.first;
1782  const Fractional lp_multi = p.second;
1783 
1784  // We ignore small values since these are likely errors and will not
1785  // contribute much to the new lp constraint anyway.
1786  if (std::abs(lp_multi) < kZeroTolerance) continue;
1787 
1788  // Remove trivial bad cases.
1789  //
1790  // TODO(user): It might be better (when possible) to use the OPTIMAL row
1791  // status since in most situation we do want the constraint we add to be
1792  // tight under the current LP solution. Only for infeasible problem we might
1793  // not have access to the status.
1794  if (lp_multi > 0.0 && integer_lp_[row].ub >= kMaxIntegerValue) {
1795  continue;
1796  }
1797  if (lp_multi < 0.0 && integer_lp_[row].lb <= kMinIntegerValue) {
1798  continue;
1799  }
1800 
1801  const Fractional cp_multi = scaler_.UnscaleDualValue(row, lp_multi);
1802  tmp_cp_multipliers_.push_back({row, cp_multi});
1803  max_sum += ToDouble(infinity_norms_[row]) * std::abs(cp_multi);
1804  }
1805 
1806  // This behave exactly like if we had another "objective" constraint with
1807  // an lp_multi of 1.0 and a cp_multi of 1.0.
1808  if (take_objective_into_account) {
1809  max_sum += ToDouble(objective_infinity_norm_);
1810  }
1811 
1812  *scaling = 1.0;
1813  std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers;
1814  if (max_sum == 0.0) {
1815  // Empty linear combinaison.
1816  return integer_multipliers;
1817  }
1818 
1819  // We want max_sum * scaling to be <= 2 ^ max_pow and fit on an int64.
1820  // We use a power of 2 as this seems to work better.
1821  const double threshold = std::ldexp(1, max_pow) / max_sum;
1822  if (threshold < 1.0) {
1823  // TODO(user): we currently do not support scaling down, so we just abort
1824  // in this case.
1825  return integer_multipliers;
1826  }
1827  while (2 * *scaling <= threshold) *scaling *= 2;
1828 
1829  // Scale the multipliers by *scaling.
1830  //
1831  // TODO(user): Maybe use int128 to avoid overflow?
1832  for (const auto entry : tmp_cp_multipliers_) {
1833  const IntegerValue coeff(std::round(entry.second * (*scaling)));
1834  if (coeff != 0) integer_multipliers.push_back({entry.first, coeff});
1835  }
1836  return integer_multipliers;
1837 }
1838 
1839 bool LinearProgrammingConstraint::ComputeNewLinearConstraint(
1840  const std::vector<std::pair<RowIndex, IntegerValue>>& integer_multipliers,
1841  ScatteredIntegerVector* scattered_vector, IntegerValue* upper_bound) const {
1842  // Initialize the new constraint.
1843  *upper_bound = 0;
1844  scattered_vector->ClearAndResize(integer_variables_.size());
1845 
1846  // Compute the new constraint by taking the linear combination given by
1847  // integer_multipliers of the integer constraints in integer_lp_.
1848  for (const std::pair<RowIndex, IntegerValue> term : integer_multipliers) {
1849  const RowIndex row = term.first;
1850  const IntegerValue multiplier = term.second;
1851  CHECK_LT(row, integer_lp_.size());
1852 
1853  // Update the constraint.
1854  if (!scattered_vector->AddLinearExpressionMultiple(
1855  multiplier, integer_lp_[row].terms)) {
1856  return false;
1857  }
1858 
1859  // Update the upper bound.
1860  const IntegerValue bound =
1861  multiplier > 0 ? integer_lp_[row].ub : integer_lp_[row].lb;
1862  if (!AddProductTo(multiplier, bound, upper_bound)) return false;
1863  }
1864 
1865  return true;
1866 }
1867 
1868 // TODO(user): no need to update the multipliers.
1869 void LinearProgrammingConstraint::AdjustNewLinearConstraint(
1870  std::vector<std::pair<glop::RowIndex, IntegerValue>>* integer_multipliers,
1871  ScatteredIntegerVector* scattered_vector, IntegerValue* upper_bound) const {
1872  const IntegerValue kMaxWantedCoeff(1e18);
1873  for (std::pair<RowIndex, IntegerValue>& term : *integer_multipliers) {
1874  const RowIndex row = term.first;
1875  const IntegerValue multiplier = term.second;
1876  if (multiplier == 0) continue;
1877 
1878  // We will only allow change of the form "multiplier += to_add" with to_add
1879  // in [-negative_limit, positive_limit].
1880  IntegerValue negative_limit = kMaxWantedCoeff;
1881  IntegerValue positive_limit = kMaxWantedCoeff;
1882 
1883  // Make sure we never change the sign of the multiplier, except if the
1884  // row is an equality in which case we don't care.
1885  if (integer_lp_[row].ub != integer_lp_[row].lb) {
1886  if (multiplier > 0) {
1887  negative_limit = std::min(negative_limit, multiplier);
1888  } else {
1889  positive_limit = std::min(positive_limit, -multiplier);
1890  }
1891  }
1892 
1893  // Make sure upper_bound + to_add * row_bound never overflow.
1894  const IntegerValue row_bound =
1895  multiplier > 0 ? integer_lp_[row].ub : integer_lp_[row].lb;
1896  if (row_bound != 0) {
1897  const IntegerValue limit1 = FloorRatio(
1898  std::max(IntegerValue(0), kMaxWantedCoeff - IntTypeAbs(*upper_bound)),
1899  IntTypeAbs(row_bound));
1900  const IntegerValue limit2 =
1901  FloorRatio(kMaxWantedCoeff, IntTypeAbs(row_bound));
1902  if ((*upper_bound > 0) == (row_bound > 0)) { // Same sign.
1903  positive_limit = std::min(positive_limit, limit1);
1904  negative_limit = std::min(negative_limit, limit2);
1905  } else {
1906  negative_limit = std::min(negative_limit, limit1);
1907  positive_limit = std::min(positive_limit, limit2);
1908  }
1909  }
1910 
1911  // If we add the row to the scattered_vector, diff will indicate by how much
1912  // |upper_bound - ImpliedLB(scattered_vector)| will change. That correspond
1913  // to increasing the multiplier by 1.
1914  //
1915  // At this stage, we are not sure computing sum coeff * bound will not
1916  // overflow, so we use floating point numbers. It is fine to do so since
1917  // this is not directly involved in the actual exact constraint generation:
1918  // these variables are just used in an heuristic.
1919  double positive_diff = ToDouble(row_bound);
1920  double negative_diff = ToDouble(row_bound);
1921 
1922  // TODO(user): we could relax a bit some of the condition and allow a sign
1923  // change. It is just trickier to compute the diff when we allow such
1924  // changes.
1925  for (const auto entry : integer_lp_[row].terms) {
1926  const ColIndex col = entry.first;
1927  const IntegerValue coeff = entry.second;
1928  const IntegerValue abs_coef = IntTypeAbs(coeff);
1929  CHECK_NE(coeff, 0);
1930 
1931  const IntegerVariable var = integer_variables_[col.value()];
1932  const IntegerValue lb = integer_trail_->LowerBound(var);
1933  const IntegerValue ub = integer_trail_->UpperBound(var);
1934 
1935  // Moving a variable away from zero seems to improve the bound even
1936  // if it reduces the number of non-zero. Note that this is because of
1937  // this that positive_diff and negative_diff are not the same.
1938  const IntegerValue current = (*scattered_vector)[col];
1939  if (current == 0) {
1940  const IntegerValue overflow_limit(
1941  FloorRatio(kMaxWantedCoeff, abs_coef));
1942  positive_limit = std::min(positive_limit, overflow_limit);
1943  negative_limit = std::min(negative_limit, overflow_limit);
1944  if (coeff > 0) {
1945  positive_diff -= ToDouble(coeff) * ToDouble(lb);
1946  negative_diff -= ToDouble(coeff) * ToDouble(ub);
1947  } else {
1948  positive_diff -= ToDouble(coeff) * ToDouble(ub);
1949  negative_diff -= ToDouble(coeff) * ToDouble(lb);
1950  }
1951  continue;
1952  }
1953 
1954  // We don't want to change the sign of current (except if the variable is
1955  // fixed) or to have an overflow.
1956  //
1957  // Corner case:
1958  // - IntTypeAbs(current) can be larger than kMaxWantedCoeff!
1959  // - The code assumes that 2 * kMaxWantedCoeff do not overflow.
1960  const IntegerValue current_magnitude = IntTypeAbs(current);
1961  const IntegerValue other_direction_limit = FloorRatio(
1962  lb == ub
1963  ? kMaxWantedCoeff + std::min(current_magnitude,
1964  kMaxIntegerValue - kMaxWantedCoeff)
1965  : current_magnitude,
1966  abs_coef);
1967  const IntegerValue same_direction_limit(FloorRatio(
1968  std::max(IntegerValue(0), kMaxWantedCoeff - current_magnitude),
1969  abs_coef));
1970  if ((current > 0) == (coeff > 0)) { // Same sign.
1971  negative_limit = std::min(negative_limit, other_direction_limit);
1972  positive_limit = std::min(positive_limit, same_direction_limit);
1973  } else {
1974  negative_limit = std::min(negative_limit, same_direction_limit);
1975  positive_limit = std::min(positive_limit, other_direction_limit);
1976  }
1977 
1978  // This is how diff change.
1979  const IntegerValue implied = current > 0 ? lb : ub;
1980  if (implied != 0) {
1981  positive_diff -= ToDouble(coeff) * ToDouble(implied);
1982  negative_diff -= ToDouble(coeff) * ToDouble(implied);
1983  }
1984  }
1985 
1986  // Only add a multiple of this row if it tighten the final constraint.
1987  // The positive_diff/negative_diff are supposed to be integer modulo the
1988  // double precision, so we only add a multiple if they seems far away from
1989  // zero.
1990  IntegerValue to_add(0);
1991  if (positive_diff <= -1.0 && positive_limit > 0) {
1992  to_add = positive_limit;
1993  }
1994  if (negative_diff >= 1.0 && negative_limit > 0) {
1995  // Pick this if it is better than the positive sign.
1996  if (to_add == 0 ||
1997  std::abs(ToDouble(negative_limit) * negative_diff) >
1998  std::abs(ToDouble(positive_limit) * positive_diff)) {
1999  to_add = -negative_limit;
2000  }
2001  }
2002  if (to_add != 0) {
2003  term.second += to_add;
2004  *upper_bound += to_add * row_bound;
2005 
2006  // TODO(user): we could avoid checking overflow here, but this is likely
2007  // not in the hot loop.
2008  CHECK(scattered_vector->AddLinearExpressionMultiple(
2009  to_add, integer_lp_[row].terms));
2010  }
2011  }
2012 }
2013 
2014 // The "exact" computation go as follow:
2015 //
2016 // Given any INTEGER linear combination of the LP constraints, we can create a
2017 // new integer constraint that is valid (its computation must not overflow
2018 // though). Lets call this "linear_combination <= ub". We can then always add to
2019 // it the inequality "objective_terms <= objective_var", so we get:
2020 // ImpliedLB(objective_terms + linear_combination) - ub <= objective_var.
2021 // where ImpliedLB() is computed from the variable current bounds.
2022 //
2023 // Now, if we use for the linear combination and approximation of the optimal
2024 // negated dual LP values (by scaling them and rounding them to integer), we
2025 // will get an EXACT objective lower bound that is more or less the same as the
2026 // inexact bound given by the LP relaxation. This allows to derive exact reasons
2027 // for any propagation done by this constraint.
2028 bool LinearProgrammingConstraint::ExactLpReasonning() {
2029  // Clear old reason and deductions.
2030  integer_reason_.clear();
2031  deductions_.clear();
2032  deductions_reason_.clear();
2033 
2034  // The row multipliers will be the negation of the LP duals.
2035  //
2036  // TODO(user): Provide and use a sparse API in Glop to get the duals.
2037  const RowIndex num_rows = simplex_.GetProblemNumRows();
2038  std::vector<std::pair<RowIndex, double>> lp_multipliers;
2039  for (RowIndex row(0); row < num_rows; ++row) {
2040  const double value = -simplex_.GetDualValue(row);
2041  if (std::abs(value) < kZeroTolerance) continue;
2042  lp_multipliers.push_back({row, value});
2043  }
2044 
2045  Fractional scaling;
2046  std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers =
2047  ScaleLpMultiplier(/*take_objective_into_account=*/true, lp_multipliers,
2048  &scaling);
2049 
2050  IntegerValue rc_ub;
2051  if (!ComputeNewLinearConstraint(integer_multipliers, &tmp_scattered_vector_,
2052  &rc_ub)) {
2053  VLOG(1) << "Issue while computing the exact LP reason. Aborting.";
2054  return true;
2055  }
2056 
2057  // The "objective constraint" behave like if the unscaled cp multiplier was
2058  // 1.0, so we will multiply it by this number and add it to reduced_costs.
2059  const IntegerValue obj_scale(std::round(scaling));
2060  if (obj_scale == 0) {
2061  VLOG(1) << "Overflow during exact LP reasoning. scaling=" << scaling;
2062  return true;
2063  }
2064  CHECK(tmp_scattered_vector_.AddLinearExpressionMultiple(obj_scale,
2065  integer_objective_));
2066  CHECK(AddProductTo(-obj_scale, integer_objective_offset_, &rc_ub));
2067  AdjustNewLinearConstraint(&integer_multipliers, &tmp_scattered_vector_,
2068  &rc_ub);
2069 
2070  // Create the IntegerSumLE that will allow to propagate the objective and more
2071  // generally do the reduced cost fixing.
2072  LinearConstraint new_constraint;
2073  tmp_scattered_vector_.ConvertToLinearConstraint(integer_variables_, rc_ub,
2074  &new_constraint);
2075  new_constraint.vars.push_back(objective_cp_);
2076  new_constraint.coeffs.push_back(-obj_scale);
2077  DivideByGCD(&new_constraint);
2078  PreventOverflow(&new_constraint);
2079  DCHECK(!PossibleOverflow(new_constraint));
2080  DCHECK(constraint_manager_.DebugCheckConstraint(new_constraint));
2081 
2082  IntegerSumLE* cp_constraint =
2083  new IntegerSumLE({}, new_constraint.vars, new_constraint.coeffs,
2084  new_constraint.ub, model_);
2085  if (trail_->CurrentDecisionLevel() == 0) {
2086  // Since we will never ask the reason for a constraint at level 0, we just
2087  // keep the last one.
2088  optimal_constraints_.clear();
2089  }
2090  optimal_constraints_.emplace_back(cp_constraint);
2091  rev_optimal_constraints_size_ = optimal_constraints_.size();
2092  if (!cp_constraint->PropagateAtLevelZero()) return false;
2093  return cp_constraint->Propagate();
2094 }
2095 
2096 bool LinearProgrammingConstraint::FillExactDualRayReason() {
2097  Fractional scaling;
2098  const glop::DenseColumn ray = simplex_.GetDualRay();
2099  std::vector<std::pair<RowIndex, double>> lp_multipliers;
2100  for (RowIndex row(0); row < ray.size(); ++row) {
2101  const double value = ray[row];
2102  if (std::abs(value) < kZeroTolerance) continue;
2103  lp_multipliers.push_back({row, value});
2104  }
2105  std::vector<std::pair<RowIndex, IntegerValue>> integer_multipliers =
2106  ScaleLpMultiplier(/*take_objective_into_account=*/false, lp_multipliers,
2107  &scaling);
2108 
2109  IntegerValue new_constraint_ub;
2110  if (!ComputeNewLinearConstraint(integer_multipliers, &tmp_scattered_vector_,
2111  &new_constraint_ub)) {
2112  VLOG(1) << "Isse while computing the exact dual ray reason. Aborting.";
2113  return false;
2114  }
2115 
2116  AdjustNewLinearConstraint(&integer_multipliers, &tmp_scattered_vector_,
2117  &new_constraint_ub);
2118 
2119  LinearConstraint new_constraint;
2120  tmp_scattered_vector_.ConvertToLinearConstraint(
2121  integer_variables_, new_constraint_ub, &new_constraint);
2122  DivideByGCD(&new_constraint);
2123  PreventOverflow(&new_constraint);
2124  DCHECK(!PossibleOverflow(new_constraint));
2125  DCHECK(constraint_manager_.DebugCheckConstraint(new_constraint));
2126 
2127  const IntegerValue implied_lb = GetImpliedLowerBound(new_constraint);
2128  if (implied_lb <= new_constraint.ub) {
2129  VLOG(1) << "LP exact dual ray not infeasible,"
2130  << " implied_lb: " << implied_lb.value() / scaling
2131  << " ub: " << new_constraint.ub.value() / scaling;
2132  return false;
2133  }
2134  const IntegerValue slack = (implied_lb - new_constraint.ub) - 1;
2135  SetImpliedLowerBoundReason(new_constraint, slack);
2136  return true;
2137 }
2138 
2139 int64 LinearProgrammingConstraint::CalculateDegeneracy() {
2140  const glop::ColIndex num_vars = simplex_.GetProblemNumCols();
2141  int num_non_basic_with_zero_rc = 0;
2142  for (glop::ColIndex i(0); i < num_vars; ++i) {
2143  const double rc = simplex_.GetReducedCost(i);
2144  if (rc != 0.0) continue;
2145  if (simplex_.GetVariableStatus(i) == glop::VariableStatus::BASIC) {
2146  continue;
2147  }
2148  num_non_basic_with_zero_rc++;
2149  }
2150  const int64 num_cols = simplex_.GetProblemNumCols().value();
2151  is_degenerate_ = num_non_basic_with_zero_rc >= 0.3 * num_cols;
2152  return num_non_basic_with_zero_rc;
2153 }
2154 
2155 void LinearProgrammingConstraint::ReducedCostStrengtheningDeductions(
2156  double cp_objective_delta) {
2157  deductions_.clear();
2158 
2159  // TRICKY: while simplex_.GetObjectiveValue() use the objective scaling factor
2160  // stored in the lp_data_, all the other functions like GetReducedCost() or
2161  // GetVariableValue() do not.
2162  const double lp_objective_delta =
2163  cp_objective_delta / lp_data_.objective_scaling_factor();
2164  const int num_vars = integer_variables_.size();
2165  for (int i = 0; i < num_vars; i++) {
2166  const IntegerVariable cp_var = integer_variables_[i];
2167  const glop::ColIndex lp_var = glop::ColIndex(i);
2168  const double rc = simplex_.GetReducedCost(lp_var);
2169  const double value = simplex_.GetVariableValue(lp_var);
2170 
2171  if (rc == 0.0) continue;
2172  const double lp_other_bound = value + lp_objective_delta / rc;
2173  const double cp_other_bound =
2174  scaler_.UnscaleVariableValue(lp_var, lp_other_bound);
2175 
2176  if (rc > kLpEpsilon) {
2177  const double ub = ToDouble(integer_trail_->UpperBound(cp_var));
2178  const double new_ub = std::floor(cp_other_bound + kCpEpsilon);
2179  if (new_ub < ub) {
2180  // TODO(user): Because rc > kLpEpsilon, the lower_bound of cp_var
2181  // will be part of the reason returned by FillReducedCostsReason(), but
2182  // we actually do not need it here. Same below.
2183  const IntegerValue new_ub_int(static_cast<IntegerValue>(new_ub));
2184  deductions_.push_back(IntegerLiteral::LowerOrEqual(cp_var, new_ub_int));
2185  }
2186  } else if (rc < -kLpEpsilon) {
2187  const double lb = ToDouble(integer_trail_->LowerBound(cp_var));
2188  const double new_lb = std::ceil(cp_other_bound - kCpEpsilon);
2189  if (new_lb > lb) {
2190  const IntegerValue new_lb_int(static_cast<IntegerValue>(new_lb));
2191  deductions_.push_back(
2192  IntegerLiteral::GreaterOrEqual(cp_var, new_lb_int));
2193  }
2194  }
2195  }
2196 }
2197 
2198 namespace {
2199 
2200 // Add a cut of the form Sum_{outgoing arcs from S} lp >= rhs_lower_bound.
2201 //
2202 // Note that we used to also add the same cut for the incoming arcs, but because
2203 // of flow conservation on these problems, the outgoing flow is always the same
2204 // as the incoming flow, so adding this extra cut doesn't seem relevant.
2205 void AddOutgoingCut(
2206  int num_nodes, int subset_size, const std::vector<bool>& in_subset,
2207  const std::vector<int>& tails, const std::vector<int>& heads,
2208  const std::vector<Literal>& literals,
2209  const std::vector<double>& literal_lp_values, int64 rhs_lower_bound,
2211  LinearConstraintManager* manager, Model* model) {
2212  // A node is said to be optional if it can be excluded from the subcircuit,
2213  // in which case there is a self-loop on that node.
2214  // If there are optional nodes, use extended formula:
2215  // sum(cut) >= 1 - optional_loop_in - optional_loop_out
2216  // where optional_loop_in's node is in subset, optional_loop_out's is out.
2217  // TODO(user): Favor optional loops fixed to zero at root.
2218  int num_optional_nodes_in = 0;
2219  int num_optional_nodes_out = 0;
2220  int optional_loop_in = -1;
2221  int optional_loop_out = -1;
2222  for (int i = 0; i < tails.size(); ++i) {
2223  if (tails[i] != heads[i]) continue;
2224  if (in_subset[tails[i]]) {
2225  num_optional_nodes_in++;
2226  if (optional_loop_in == -1 ||
2227  literal_lp_values[i] < literal_lp_values[optional_loop_in]) {
2228  optional_loop_in = i;
2229  }
2230  } else {
2231  num_optional_nodes_out++;
2232  if (optional_loop_out == -1 ||
2233  literal_lp_values[i] < literal_lp_values[optional_loop_out]) {
2234  optional_loop_out = i;
2235  }
2236  }
2237  }
2238 
2239  // TODO(user): The lower bound for CVRP is computed assuming all nodes must be
2240  // served, if it is > 1 we lower it to one in the presence of optional nodes.
2241  if (num_optional_nodes_in + num_optional_nodes_out > 0) {
2242  CHECK_GE(rhs_lower_bound, 1);
2243  rhs_lower_bound = 1;
2244  }
2245 
2246  LinearConstraintBuilder outgoing(model, IntegerValue(rhs_lower_bound),
2248  double sum_outgoing = 0.0;
2249 
2250  // Add outgoing arcs, compute outgoing flow.
2251  for (int i = 0; i < tails.size(); ++i) {
2252  if (in_subset[tails[i]] && !in_subset[heads[i]]) {
2253  sum_outgoing += literal_lp_values[i];
2254  CHECK(outgoing.AddLiteralTerm(literals[i], IntegerValue(1)));
2255  }
2256  }
2257 
2258  // Support optional nodes if any.
2259  if (num_optional_nodes_in + num_optional_nodes_out > 0) {
2260  // When all optionals of one side are excluded in lp solution, no cut.
2261  if (num_optional_nodes_in == subset_size &&
2262  (optional_loop_in == -1 ||
2263  literal_lp_values[optional_loop_in] > 1.0 - 1e-6)) {
2264  return;
2265  }
2266  if (num_optional_nodes_out == num_nodes - subset_size &&
2267  (optional_loop_out == -1 ||
2268  literal_lp_values[optional_loop_out] > 1.0 - 1e-6)) {
2269  return;
2270  }
2271 
2272  // There is no mandatory node in subset, add optional_loop_in.
2273  if (num_optional_nodes_in == subset_size) {
2274  CHECK(
2275  outgoing.AddLiteralTerm(literals[optional_loop_in], IntegerValue(1)));
2276  sum_outgoing += literal_lp_values[optional_loop_in];
2277  }
2278 
2279  // There is no mandatory node out of subset, add optional_loop_out.
2280  if (num_optional_nodes_out == num_nodes - subset_size) {
2281  CHECK(outgoing.AddLiteralTerm(literals[optional_loop_out],
2282  IntegerValue(1)));
2283  sum_outgoing += literal_lp_values[optional_loop_out];
2284  }
2285  }
2286 
2287  if (sum_outgoing < rhs_lower_bound - 1e-6) {
2288  manager->AddCut(outgoing.Build(), "Circuit", lp_values);
2289  }
2290 }
2291 
2292 } // namespace
2293 
2294 // We roughly follow the algorithm described in section 6 of "The Traveling
2295 // Salesman Problem, A computational Study", David L. Applegate, Robert E.
2296 // Bixby, Vasek Chvatal, William J. Cook.
2297 //
2298 // Note that this is mainly a "symmetric" case algo, but it does still work for
2299 // the asymmetric case.
2301  int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
2302  const std::vector<Literal>& literals,
2304  absl::Span<const int64> demands, int64 capacity,
2305  LinearConstraintManager* manager, Model* model) {
2306  if (num_nodes <= 2) return;
2307 
2308  // We will collect only the arcs with a positive lp_values to speed up some
2309  // computation below.
2310  struct Arc {
2311  int tail;
2312  int head;
2313  double lp_value;
2314  };
2315  std::vector<Arc> relevant_arcs;
2316 
2317  // Sort the arcs by non-increasing lp_values.
2318  std::vector<double> literal_lp_values(literals.size());
2319  std::vector<std::pair<double, int>> arc_by_decreasing_lp_values;
2320  auto* encoder = model->GetOrCreate<IntegerEncoder>();
2321  for (int i = 0; i < literals.size(); ++i) {
2322  double lp_value;
2323  const IntegerVariable direct_view = encoder->GetLiteralView(literals[i]);
2324  if (direct_view != kNoIntegerVariable) {
2325  lp_value = lp_values[direct_view];
2326  } else {
2327  lp_value =
2328  1.0 - lp_values[encoder->GetLiteralView(literals[i].Negated())];
2329  }
2330  literal_lp_values[i] = lp_value;
2331 
2332  if (lp_value < 1e-6) continue;
2333  relevant_arcs.push_back({tails[i], heads[i], lp_value});
2334  arc_by_decreasing_lp_values.push_back({lp_value, i});
2335  }
2336  std::sort(arc_by_decreasing_lp_values.begin(),
2337  arc_by_decreasing_lp_values.end(),
2338  std::greater<std::pair<double, int>>());
2339 
2340  // We will do a union-find by adding one by one the arc of the lp solution
2341  // in the order above. Every intermediate set during this construction will
2342  // be a candidate for a cut.
2343  //
2344  // In parallel to the union-find, to efficiently reconstruct these sets (at
2345  // most num_nodes), we construct a "decomposition forest" of the different
2346  // connected components. Note that we don't exploit any asymmetric nature of
2347  // the graph here. This is exactly the algo 6.3 in the book above.
2348  int num_components = num_nodes;
2349  std::vector<int> parent(num_nodes);
2350  std::vector<int> root(num_nodes);
2351  for (int i = 0; i < num_nodes; ++i) {
2352  parent[i] = i;
2353  root[i] = i;
2354  }
2355  auto get_root_and_compress_path = [&root](int node) {
2356  int r = node;
2357  while (root[r] != r) r = root[r];
2358  while (root[node] != r) {
2359  const int next = root[node];
2360  root[node] = r;
2361  node = next;
2362  }
2363  return r;
2364  };
2365  for (const auto pair : arc_by_decreasing_lp_values) {
2366  if (num_components == 2) break;
2367  const int tail = get_root_and_compress_path(tails[pair.second]);
2368  const int head = get_root_and_compress_path(heads[pair.second]);
2369  if (tail != head) {
2370  // Update the decomposition forest, note that the number of nodes is
2371  // growing.
2372  const int new_node = parent.size();
2373  parent.push_back(new_node);
2374  parent[head] = new_node;
2375  parent[tail] = new_node;
2376  --num_components;
2377 
2378  // It is important that the union-find representative is the same node.
2379  root.push_back(new_node);
2380  root[head] = new_node;
2381  root[tail] = new_node;
2382  }
2383  }
2384 
2385  // For each node in the decomposition forest, try to add a cut for the set
2386  // formed by the nodes and its children. To do that efficiently, we first
2387  // order the nodes so that for each node in a tree, the set of children forms
2388  // a consecutive span in the pre_order vector. This vector just lists the
2389  // nodes in the "pre-order" graph traversal order. The Spans will point inside
2390  // the pre_order vector, it is why we initialize it once and for all.
2391  int new_size = 0;
2392  std::vector<int> pre_order(num_nodes);
2393  std::vector<absl::Span<const int>> subsets;
2394  {
2395  std::vector<absl::InlinedVector<int, 2>> graph(parent.size());
2396  for (int i = 0; i < parent.size(); ++i) {
2397  if (parent[i] != i) graph[parent[i]].push_back(i);
2398  }
2399  std::vector<int> queue;
2400  std::vector<bool> seen(graph.size(), false);
2401  std::vector<int> start_index(parent.size());
2402  for (int i = num_nodes; i < parent.size(); ++i) {
2403  // Note that because of the way we constructed 'parent', the graph is a
2404  // binary tree. This is not required for the correctness of the algorithm
2405  // here though.
2406  CHECK(graph[i].empty() || graph[i].size() == 2);
2407  if (parent[i] != i) continue;
2408 
2409  // Explore the subtree rooted at node i.
2410  CHECK(!seen[i]);
2411  queue.push_back(i);
2412  while (!queue.empty()) {
2413  const int node = queue.back();
2414  if (seen[node]) {
2415  queue.pop_back();
2416  // All the children of node are in the span [start, end) of the
2417  // pre_order vector.
2418  const int start = start_index[node];
2419  if (new_size - start > 1) {
2420  subsets.emplace_back(&pre_order[start], new_size - start);
2421  }
2422  continue;
2423  }
2424  seen[node] = true;
2425  start_index[node] = new_size;
2426  if (node < num_nodes) pre_order[new_size++] = node;
2427  for (const int child : graph[node]) {
2428  if (!seen[child]) queue.push_back(child);
2429  }
2430  }
2431  }
2432  }
2433 
2434  // Compute the total demands in order to know the minimum incoming/outgoing
2435  // flow.
2436  int64 total_demands = 0;
2437  if (!demands.empty()) {
2438  for (const int64 demand : demands) total_demands += demand;
2439  }
2440 
2441  // Process each subsets and add any violated cut.
2442  CHECK_EQ(pre_order.size(), num_nodes);
2443  std::vector<bool> in_subset(num_nodes, false);
2444  for (const absl::Span<const int> subset : subsets) {
2445  CHECK_GT(subset.size(), 1);
2446  CHECK_LT(subset.size(), num_nodes);
2447 
2448  // These fields will be left untouched if demands.empty().
2449  bool contain_depot = false;
2450  int64 subset_demand = 0;
2451 
2452  // Initialize "in_subset" and the subset demands.
2453  for (const int n : subset) {
2454  in_subset[n] = true;
2455  if (!demands.empty()) {
2456  if (n == 0) contain_depot = true;
2457  subset_demand += demands[n];
2458  }
2459  }
2460 
2461  // Compute a lower bound on the outgoing flow.
2462  //
2463  // TODO(user): This lower bound assume all nodes in subset must be served,
2464  // which is not the case. For TSP we do the correct thing in
2465  // AddOutgoingCut() but not for CVRP... Fix!!
2466  //
2467  // TODO(user): It could be very interesting to see if this "min outgoing
2468  // flow" cannot be automatically infered from the constraint in the
2469  // precedence graph. This might work if we assume that any kind of path
2470  // cumul constraint is encoded with constraints:
2471  // [edge => value_head >= value_tail + edge_weight].
2472  // We could take the minimum incoming edge weight per node in the set, and
2473  // use the cumul variable domain to infer some capacity.
2474  int64 min_outgoing_flow = 1;
2475  if (!demands.empty()) {
2476  min_outgoing_flow =
2477  contain_depot
2478  ? (total_demands - subset_demand + capacity - 1) / capacity
2479  : (subset_demand + capacity - 1) / capacity;
2480  }
2481 
2482  // We still need to serve nodes with a demand of zero, and in the corner
2483  // case where all node in subset have a zero demand, the formula above
2484  // result in a min_outgoing_flow of zero.
2485  min_outgoing_flow = std::max(min_outgoing_flow, int64{1});
2486 
2487  // Compute the current outgoing flow out of the subset.
2488  //
2489  // This can take a significant portion of the running time, it is why it is
2490  // faster to do it only on arcs with non-zero lp values which should be in
2491  // linear number rather than the total number of arc which can be quadratic.
2492  //
2493  // TODO(user): For the symmetric case there is an even faster algo. See if
2494  // it can be generalized to the asymmetric one if become needed.
2495  // Reference is algo 6.4 of the "The Traveling Salesman Problem" book
2496  // mentionned above.
2497  double outgoing_flow = 0.0;
2498  for (const auto arc : relevant_arcs) {
2499  if (in_subset[arc.tail] && !in_subset[arc.head]) {
2500  outgoing_flow += arc.lp_value;
2501  }
2502  }
2503 
2504  // Add a cut if the current outgoing flow is not enough.
2505  if (outgoing_flow < min_outgoing_flow - 1e-6) {
2506  AddOutgoingCut(num_nodes, subset.size(), in_subset, tails, heads,
2507  literals, literal_lp_values,
2508  /*rhs_lower_bound=*/min_outgoing_flow, lp_values, manager,
2509  model);
2510  }
2511 
2512  // Sparse clean up.
2513  for (const int n : subset) in_subset[n] = false;
2514  }
2515 }
2516 
2517 namespace {
2518 
2519 // Returns for each literal its integer view, or the view of its negation.
2520 std::vector<IntegerVariable> GetAssociatedVariables(
2521  const std::vector<Literal>& literals, Model* model) {
2522  auto* encoder = model->GetOrCreate<IntegerEncoder>();
2523  std::vector<IntegerVariable> result;
2524  for (const Literal l : literals) {
2525  const IntegerVariable direct_view = encoder->GetLiteralView(l);
2526  if (direct_view != kNoIntegerVariable) {
2527  result.push_back(direct_view);
2528  } else {
2529  result.push_back(encoder->GetLiteralView(l.Negated()));
2530  DCHECK_NE(result.back(), kNoIntegerVariable);
2531  }
2532  }
2533  return result;
2534 }
2535 
2536 } // namespace
2537 
2538 // We use a basic algorithm to detect components that are not connected to the
2539 // rest of the graph in the LP solution, and add cuts to force some arcs to
2540 // enter and leave this component from outside.
2542  int num_nodes, const std::vector<int>& tails, const std::vector<int>& heads,
2543  const std::vector<Literal>& literals, Model* model) {
2544  CutGenerator result;
2545  result.vars = GetAssociatedVariables(literals, model);
2546  result.generate_cuts =
2547  [num_nodes, tails, heads, literals, model](
2549  LinearConstraintManager* manager) {
2551  num_nodes, tails, heads, literals, lp_values,
2552  /*demands=*/{}, /*capacity=*/0, manager, model);
2553  };
2554  return result;
2555 }
2556 
2558  const std::vector<int>& tails,
2559  const std::vector<int>& heads,
2560  const std::vector<Literal>& literals,
2561  const std::vector<int64>& demands,
2562  int64 capacity, Model* model) {
2563  CutGenerator result;
2564  result.vars = GetAssociatedVariables(literals, model);
2565  result.generate_cuts =
2566  [num_nodes, tails, heads, demands, capacity, literals, model](
2568  LinearConstraintManager* manager) {
2569  SeparateSubtourInequalities(num_nodes, tails, heads, literals,
2570  lp_values, demands, capacity, manager,
2571  model);
2572  };
2573  return result;
2574 }
2575 
2576 std::function<IntegerLiteral()>
2578  // Gather all 0-1 variables that appear in this LP.
2579  std::vector<IntegerVariable> variables;
2580  for (IntegerVariable var : integer_variables_) {
2581  if (integer_trail_->LowerBound(var) == 0 &&
2582  integer_trail_->UpperBound(var) == 1) {
2583  variables.push_back(var);
2584  }
2585  }
2586  VLOG(1) << "HeuristicLPMostInfeasibleBinary has " << variables.size()
2587  << " variables.";
2588 
2589  return [this, variables]() {
2590  const double kEpsilon = 1e-6;
2591  // Find most fractional value.
2592  IntegerVariable fractional_var = kNoIntegerVariable;
2593  double fractional_distance_best = -1.0;
2594  for (const IntegerVariable var : variables) {
2595  // Skip ignored and fixed variables.
2596  if (integer_trail_->IsCurrentlyIgnored(var)) continue;
2597  const IntegerValue lb = integer_trail_->LowerBound(var);
2598  const IntegerValue ub = integer_trail_->UpperBound(var);
2599  if (lb == ub) continue;
2600 
2601  // Check variable's support is fractional.
2602  const double lp_value = this->GetSolutionValue(var);
2603  const double fractional_distance =
2604  std::min(std::ceil(lp_value - kEpsilon) - lp_value,
2605  lp_value - std::floor(lp_value + kEpsilon));
2606  if (fractional_distance < kEpsilon) continue;
2607 
2608  // Keep variable if it is farther from integrality than the previous.
2609  if (fractional_distance > fractional_distance_best) {
2610  fractional_var = var;
2611  fractional_distance_best = fractional_distance;
2612  }
2613  }
2614 
2615  if (fractional_var != kNoIntegerVariable) {
2616  IntegerLiteral::GreaterOrEqual(fractional_var, IntegerValue(1));
2617  }
2618  return IntegerLiteral();
2619  };
2620 }
2621 
2622 std::function<IntegerLiteral()>
2624  // Gather all 0-1 variables that appear in this LP.
2625  std::vector<IntegerVariable> variables;
2626  for (IntegerVariable var : integer_variables_) {
2627  if (integer_trail_->LowerBound(var) == 0 &&
2628  integer_trail_->UpperBound(var) == 1) {
2629  variables.push_back(var);
2630  }
2631  }
2632  VLOG(1) << "HeuristicLpReducedCostBinary has " << variables.size()
2633  << " variables.";
2634 
2635  // Store average of reduced cost from 1 to 0. The best heuristic only sets
2636  // variables to one and cares about cost to zero, even though classic
2637  // pseudocost will use max_var min(cost_to_one[var], cost_to_zero[var]).
2638  const int num_vars = variables.size();
2639  std::vector<double> cost_to_zero(num_vars, 0.0);
2640  std::vector<int> num_cost_to_zero(num_vars);
2641  int num_calls = 0;
2642 
2643  return [=]() mutable {
2644  const double kEpsilon = 1e-6;
2645 
2646  // Every 10000 calls, decay pseudocosts.
2647  num_calls++;
2648  if (num_calls == 10000) {
2649  for (int i = 0; i < num_vars; i++) {
2650  cost_to_zero[i] /= 2;
2651  num_cost_to_zero[i] /= 2;
2652  }
2653  num_calls = 0;
2654  }
2655 
2656  // Accumulate pseudo-costs of all unassigned variables.
2657  for (int i = 0; i < num_vars; i++) {
2658  const IntegerVariable var = variables[i];
2659  // Skip ignored and fixed variables.
2660  if (integer_trail_->IsCurrentlyIgnored(var)) continue;
2661  const IntegerValue lb = integer_trail_->LowerBound(var);
2662  const IntegerValue ub = integer_trail_->UpperBound(var);
2663  if (lb == ub) continue;
2664 
2665  const double rc = this->GetSolutionReducedCost(var);
2666  // Skip reduced costs that are nonzero because of numerical issues.
2667  if (std::abs(rc) < kEpsilon) continue;
2668 
2669  const double value = std::round(this->GetSolutionValue(var));
2670  if (value == 1.0 && rc < 0.0) {
2671  cost_to_zero[i] -= rc;
2672  num_cost_to_zero[i]++;
2673  }
2674  }
2675 
2676  // Select noninstantiated variable with highest pseudo-cost.
2677  int selected_index = -1;
2678  double best_cost = 0.0;
2679  for (int i = 0; i < num_vars; i++) {
2680  const IntegerVariable var = variables[i];
2681  // Skip ignored and fixed variables.
2682  if (integer_trail_->IsCurrentlyIgnored(var)) continue;
2683  if (integer_trail_->IsFixed(var)) continue;
2684 
2685  if (num_cost_to_zero[i] > 0 &&
2686  best_cost < cost_to_zero[i] / num_cost_to_zero[i]) {
2687  best_cost = cost_to_zero[i] / num_cost_to_zero[i];
2688  selected_index = i;
2689  }
2690  }
2691 
2692  if (selected_index >= 0) {
2693  return IntegerLiteral::GreaterOrEqual(variables[selected_index],
2694  IntegerValue(1));
2695  }
2696  return IntegerLiteral();
2697  };
2698 }
2699 
2700 void LinearProgrammingConstraint::UpdateAverageReducedCosts() {
2701  const int num_vars = integer_variables_.size();
2702  if (sum_cost_down_.size() < num_vars) {
2703  sum_cost_down_.resize(num_vars, 0.0);
2704  num_cost_down_.resize(num_vars, 0);
2705  sum_cost_up_.resize(num_vars, 0.0);
2706  num_cost_up_.resize(num_vars, 0);
2707  rc_scores_.resize(num_vars, 0.0);
2708  }
2709 
2710  // Decay averages.
2711  num_calls_since_reduced_cost_averages_reset_++;
2712  if (num_calls_since_reduced_cost_averages_reset_ == 10000) {
2713  for (int i = 0; i < num_vars; i++) {
2714  sum_cost_up_[i] /= 2;
2715  num_cost_up_[i] /= 2;
2716  sum_cost_down_[i] /= 2;
2717  num_cost_down_[i] /= 2;
2718  }
2719  num_calls_since_reduced_cost_averages_reset_ = 0;
2720  }
2721 
2722  // Accumulate reduced costs of all unassigned variables.
2723  for (int i = 0; i < num_vars; i++) {
2724  const IntegerVariable var = integer_variables_[i];
2725 
2726  // Skip ignored and fixed variables.
2727  if (integer_trail_->IsCurrentlyIgnored(var)) continue;
2728  if (integer_trail_->IsFixed(var)) continue;
2729 
2730  // Skip reduced costs that are zero or close.
2731  const double rc = lp_reduced_cost_[i];
2732  if (std::abs(rc) < kCpEpsilon) continue;
2733 
2734  if (rc < 0.0) {
2735  sum_cost_down_[i] -= rc;
2736  num_cost_down_[i]++;
2737  } else {
2738  sum_cost_up_[i] += rc;
2739  num_cost_up_[i]++;
2740  }
2741  }
2742 
2743  // Tricky, we artificially reset the rc_rev_int_repository_ to level zero
2744  // so that the rev_rc_start_ is zero.
2745  rc_rev_int_repository_.SetLevel(0);
2746  rc_rev_int_repository_.SetLevel(trail_->CurrentDecisionLevel());
2747  rev_rc_start_ = 0;
2748 
2749  // Cache the new score (higher is better) using the average reduced costs
2750  // as a signal.
2751  positions_by_decreasing_rc_score_.clear();
2752  for (int i = 0; i < num_vars; i++) {
2753  // If only one direction exist, we takes its value divided by 2, so that
2754  // such variable should have a smaller cost than the min of the two side
2755  // except if one direction have a really high reduced costs.
2756  const double a_up =
2757  num_cost_up_[i] > 0 ? sum_cost_up_[i] / num_cost_up_[i] : 0.0;
2758  const double a_down =
2759  num_cost_down_[i] > 0 ? sum_cost_down_[i] / num_cost_down_[i] : 0.0;
2760  if (num_cost_down_[i] > 0 && num_cost_up_[i] > 0) {
2761  rc_scores_[i] = std::min(a_up, a_down);
2762  } else {
2763  rc_scores_[i] = 0.5 * (a_down + a_up);
2764  }
2765 
2766  // We ignore scores of zero (i.e. no data) and will follow the default
2767  // search heuristic if all variables are like this.
2768  if (rc_scores_[i] > 0.0) {
2769  positions_by_decreasing_rc_score_.push_back({-rc_scores_[i], i});
2770  }
2771  }
2772  std::sort(positions_by_decreasing_rc_score_.begin(),
2773  positions_by_decreasing_rc_score_.end());
2774 }
2775 
2776 // TODO(user): Remove duplication with HeuristicLpReducedCostBinary().
2777 std::function<IntegerLiteral()>
2779  return [this]() { return this->LPReducedCostAverageDecision(); };
2780 }
2781 
2782 IntegerLiteral LinearProgrammingConstraint::LPReducedCostAverageDecision() {
2783  // Select noninstantiated variable with highest positive average reduced cost.
2784  int selected_index = -1;
2785  const int size = positions_by_decreasing_rc_score_.size();
2786  rc_rev_int_repository_.SaveState(&rev_rc_start_);
2787  for (int i = rev_rc_start_; i < size; ++i) {
2788  const int index = positions_by_decreasing_rc_score_[i].second;
2789  const IntegerVariable var = integer_variables_[index];
2790  if (integer_trail_->IsCurrentlyIgnored(var)) continue;
2791  if (integer_trail_->IsFixed(var)) continue;
2792  selected_index = index;
2793  rev_rc_start_ = i;
2794  break;
2795  }
2796 
2797  if (selected_index == -1) return IntegerLiteral();
2798  const IntegerVariable var = integer_variables_[selected_index];
2799 
2800  // If ceil(value) is current upper bound, try var == upper bound first.
2801  // Guarding with >= prevents numerical problems.
2802  // With 0/1 variables, this will tend to try setting to 1 first,
2803  // which produces more shallow trees.
2804  const IntegerValue ub = integer_trail_->UpperBound(var);
2805  const IntegerValue value_ceil(
2806  std::ceil(this->GetSolutionValue(var) - kCpEpsilon));
2807  if (value_ceil >= ub) {
2808  return IntegerLiteral::GreaterOrEqual(var, ub);
2809  }
2810 
2811  // If floor(value) is current lower bound, try var == lower bound first.
2812  // Guarding with <= prevents numerical problems.
2813  const IntegerValue lb = integer_trail_->LowerBound(var);
2814  const IntegerValue value_floor(
2815  std::floor(this->GetSolutionValue(var) + kCpEpsilon));
2816  if (value_floor <= lb) {
2817  return IntegerLiteral::LowerOrEqual(var, lb);
2818  }
2819 
2820  // Here lb < value_floor <= value_ceil < ub.
2821  // Try the most promising split between var <= floor or var >= ceil.
2822  const double a_up =
2823  num_cost_up_[selected_index] > 0
2824  ? sum_cost_up_[selected_index] / num_cost_up_[selected_index]
2825  : 0.0;
2826  const double a_down =
2827  num_cost_down_[selected_index] > 0
2828  ? sum_cost_down_[selected_index] / num_cost_down_[selected_index]
2829  : 0.0;
2830  if (a_down < a_up) {
2831  return IntegerLiteral::LowerOrEqual(var, value_floor);
2832  } else {
2833  return IntegerLiteral::GreaterOrEqual(var, value_ceil);
2834  }
2835 }
2836 
2837 } // namespace sat
2838 } // 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 CHECK_NE(val1, val2)
Definition: base/logging.h:698
#define DCHECK_GT(val1, val2)
Definition: base/logging.h:890
#define DCHECK(condition)
Definition: base/logging.h:884
#define VLOG(verboselevel)
Definition: base/logging.h:978
void assign(size_type n, const value_type &val)
void resize(size_type new_size)
size_type size() const
void push_back(const value_type &x)
static int64 GCD64(int64 x, int64 y)
Definition: mathutil.h:107
void SetLevel(int level) final
Definition: rev.h:134
void SaveState(T *object)
Definition: rev.h:61
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
void SetVariableBounds(ColIndex col, Fractional lower_bound, Fractional upper_bound)
Definition: lp_data.cc:248
void SetObjectiveOffset(Fractional objective_offset)
Definition: lp_data.cc:330
void SetCoefficient(RowIndex row, ColIndex col, Fractional value)
Definition: lp_data.cc:316
void SetConstraintBounds(RowIndex row, Fractional lower_bound, Fractional upper_bound)
Definition: lp_data.cc:308
void AddSlackVariablesWhereNecessary(bool detect_integer_constraints)
Definition: lp_data.cc:695
void SetObjectiveCoefficient(ColIndex col, Fractional value)
Definition: lp_data.cc:325
std::string GetDimensionString() const
Definition: lp_data.cc:424
Fractional objective_scaling_factor() const
Definition: lp_data.h:260
const SparseColumn & GetSparseColumn(ColIndex col) const
Definition: lp_data.cc:408
Fractional VariableScalingFactor(ColIndex col) const
Fractional UnscaleVariableValue(ColIndex col, Fractional value) const
Fractional UnscaleReducedCost(ColIndex col, Fractional value) const
Fractional UnscaleDualValue(RowIndex row, Fractional value) const
const GlopParameters & GetParameters() const
const DenseRow & GetDualRayRowCombination() const
Fractional GetVariableValue(ColIndex col) const
void SetIntegralityScale(ColIndex col, Fractional scale)
VariableStatus GetVariableStatus(ColIndex col) const
Fractional GetReducedCost(ColIndex col) const
const DenseColumn & GetDualRay() const
ABSL_MUST_USE_RESULT Status Solve(const LinearProgram &lp, TimeLimit *time_limit)
Fractional GetDualValue(RowIndex row) const
ConstraintStatus GetConstraintStatus(RowIndex row) const
void LoadStateForNextSolve(const BasisState &state)
ColIndex GetBasis(RowIndex row) const
void SetParameters(const GlopParameters &parameters)
const ScatteredRow & GetUnitRowLeftInverse(RowIndex row)
LinearConstraint * mutable_cut()
Definition: cuts.h:251
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 WatchIntegerVariable(IntegerVariable i, int id, int watch_index=-1)
Definition: integer.h:1397
void WatchUpperBound(IntegerVariable var, int id, int watch_index=-1)
Definition: integer.h:1391
void SetPropagatorPriority(int id, int priority)
Definition: integer.cc:1962
int Register(PropagatorInterface *propagator)
Definition: integer.cc:1939
void AddLpVariable(IntegerVariable var)
Definition: cuts.h:108
void ProcessUpperBoundedConstraintWithSlackCreation(bool substitute_only_inner_variables, IntegerVariable first_slack, const absl::StrongVector< IntegerVariable, double > &lp_values, LinearConstraint *cut, std::vector< SlackInfo > *slack_infos)
Definition: cuts.cc:1584
bool DebugSlack(IntegerVariable first_slack, const LinearConstraint &initial_cut, const LinearConstraint &cut, const std::vector< SlackInfo > &info)
Definition: cuts.cc:1725
void SeparateSomeImpliedBoundCuts(const absl::StrongVector< IntegerVariable, double > &lp_values)
Definition: cuts.cc:1575
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
ABSL_MUST_USE_RESULT bool Enqueue(IntegerLiteral i_lit, absl::Span< const Literal > literal_reason, absl::Span< const IntegerLiteral > integer_reason)
Definition: integer.cc:989
bool IsCurrentlyIgnored(IntegerVariable i) const
Definition: integer.h:625
bool IsFixed(IntegerVariable i) const
Definition: integer.h:1308
IntegerLiteral LowerBoundAsLiteral(IntegerVariable i) const
Definition: integer.h:1330
bool ReportConflict(absl::Span< const Literal > literal_reason, absl::Span< const IntegerLiteral > integer_reason)
Definition: integer.h:810
IntegerValue UpperBound(IntegerVariable i) const
Definition: integer.h:1304
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
Definition: integer.h:1355
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
Definition: integer.h:1350
void RelaxLinearReason(IntegerValue slack, absl::Span< const IntegerValue > coeffs, std::vector< IntegerLiteral > *reason) const
Definition: integer.cc:785
IntegerValue LowerBound(IntegerVariable i) const
Definition: integer.h:1300
IntegerLiteral UpperBoundAsLiteral(IntegerVariable i) const
Definition: integer.h:1335
bool IsFixedAtLevelZero(IntegerVariable var) const
Definition: integer.h:1360
void RemoveLevelZeroBounds(std::vector< IntegerLiteral > *reason) const
Definition: integer.cc:919
void RegisterReversibleClass(ReversibleInterface *rev)
Definition: integer.h:833
bool ChangeLp(const absl::StrongVector< IntegerVariable, double > &lp_solution, glop::BasisState *solution_state)
void SetObjectiveCoefficient(IntegerVariable var, IntegerValue coeff)
ConstraintIndex Add(LinearConstraint ct, bool *added=nullptr)
const absl::StrongVector< ConstraintIndex, ConstraintInfo > & AllConstraints() const
const std::vector< ConstraintIndex > & LpConstraints() const
bool AddCut(LinearConstraint ct, std::string type_name, const absl::StrongVector< IntegerVariable, double > &lp_solution, std::string extra_info="")
std::function< IntegerLiteral()> HeuristicLpReducedCostBinary(Model *model)
bool IncrementalPropagate(const std::vector< int > &watch_indices) override
std::function< IntegerLiteral()> HeuristicLpMostInfeasibleBinary(Model *model)
void SetObjectiveCoefficient(IntegerVariable ivar, IntegerValue coeff)
Class that owns everything related to a particular optimization model.
Definition: sat/model.h:38
void ConvertToLinearConstraint(const std::vector< IntegerVariable > &integer_variables, IntegerValue upper_bound, LinearConstraint *result)
bool Add(glop::ColIndex col, IntegerValue value)
std::vector< std::pair< glop::ColIndex, IntegerValue > > GetTerms()
bool AddLinearExpressionMultiple(IntegerValue multiplier, const std::vector< std::pair< glop::ColIndex, IntegerValue >> &terms)
void TransferToManager(const absl::StrongVector< IntegerVariable, double > &lp_solution, LinearConstraintManager *manager)
void ProcessVariables(const std::vector< double > &lp_values, const std::vector< IntegerValue > &lower_bounds, const std::vector< IntegerValue > &upper_bounds)
std::vector< std::vector< std::pair< glop::RowIndex, IntegerValue > > > InterestingCandidates(ModelRandomGenerator *random)
void AddOneConstraint(glop::RowIndex, const std::vector< std::pair< glop::ColIndex, IntegerValue >> &terms, IntegerValue lb, IntegerValue ub)
Block * next
SatParameters parameters
const std::string name
const Constraint * ct
int64 value
IntVar * var
Definition: expr_array.cc:1858
GRBmodel * model
static const int64 kint64max
int64_t int64
static const int64 kint64min
const bool DEBUG_MODE
Definition: macros.h:24
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
const Collection::value_type::second_type & FindOrDie(const Collection &collection, const typename Collection::value_type::first_type &key)
Definition: map_util.h:176
StrictITIVector< ColIndex, Fractional > DenseRow
Definition: lp_types.h:299
ColIndex RowToColIndex(RowIndex row)
Definition: lp_types.h:48
RowIndex ColToRowIndex(ColIndex col)
Definition: lp_types.h:51
const double kEpsilon
Definition: lp_types.h:86
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
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)
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
void SeparateSubtourInequalities(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, const absl::StrongVector< IntegerVariable, double > &lp_values, absl::Span< const int64 > demands, int64 capacity, LinearConstraintManager *manager, Model *model)
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue)
const IntegerVariable kNoIntegerVariable(-1)
IntegerVariable PositiveVariable(IntegerVariable i)
Definition: integer.h:134
CutGenerator CreateCVRPCutGenerator(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, const std::vector< int64 > &demands, int64 capacity, Model *model)
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
Definition: integer.cc:27
IntegerValue ComputeInfinityNorm(const LinearConstraint &constraint)
bool VariableIsPositive(IntegerVariable i)
Definition: integer.h:130
void DivideByGCD(LinearConstraint *constraint)
CutGenerator CreateStronglyConnectedGraphCutGenerator(int num_nodes, const std::vector< int > &tails, const std::vector< int > &heads, const std::vector< Literal > &literals, Model *model)
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)
std::pair< int64, int64 > Arc
Definition: search.cc:3361
int64 CapAdd(int64 x, int64 y)
int64 CapProd(int64 x, int64 y)
int index
Definition: pack.cc:508
int64 demand
Definition: resource.cc:123
int64 tail
int64 head
int64 capacity
int64 bound
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
static IntegerLiteral LowerOrEqual(IntegerVariable i, IntegerValue bound)
Definition: integer.h:1270
static IntegerLiteral GreaterOrEqual(IntegerVariable i, IntegerValue bound)
Definition: integer.h:1264
std::vector< IntegerVariable > vars