OR-Tools  8.2
linear_constraint_manager.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 <limits>
19 #include <utility>
20 
21 #include "absl/container/flat_hash_set.h"
23 #include "ortools/sat/integer.h"
25 
26 namespace operations_research {
27 namespace sat {
28 
29 namespace {
30 
31 const LinearConstraintManager::ConstraintIndex kInvalidConstraintIndex(-1);
32 
33 size_t ComputeHashOfTerms(const LinearConstraint& ct) {
34  DCHECK(std::is_sorted(ct.vars.begin(), ct.vars.end()));
35  size_t hash = 0;
36  const int num_terms = ct.vars.size();
37  for (int i = 0; i < num_terms; ++i) {
38  hash = util_hash::Hash(ct.vars[i].value(), hash);
39  hash = util_hash::Hash(ct.coeffs[i].value(), hash);
40  }
41  return hash;
42 }
43 
44 } // namespace
45 
47  if (num_merged_constraints_ > 0) {
48  VLOG(2) << "num_merged_constraints: " << num_merged_constraints_;
49  }
50  if (num_shortened_constraints_ > 0) {
51  VLOG(2) << "num_shortened_constraints: " << num_shortened_constraints_;
52  }
53  if (num_splitted_constraints_ > 0) {
54  VLOG(2) << "num_splitted_constraints: " << num_splitted_constraints_;
55  }
56  if (num_coeff_strenghtening_ > 0) {
57  VLOG(2) << "num_coeff_strenghtening: " << num_coeff_strenghtening_;
58  }
59  if (sat_parameters_.log_search_progress() && num_cuts_ > 0) {
60  LOG(INFO) << "Total cuts added: " << num_cuts_ << " (out of "
61  << num_add_cut_calls_ << " calls) worker: '" << model_->Name()
62  << "'";
63  LOG(INFO) << " - num simplifications: " << num_simplifications_;
64  for (const auto& entry : type_to_num_cuts_) {
65  if (entry.second == 1) {
66  LOG(INFO) << " - added 1 cut of type '" << entry.first << "'.";
67  } else {
68  LOG(INFO) << " - added " << entry.second << " cuts of type '"
69  << entry.first << "'.";
70  }
71  }
72  }
73 }
74 
75 void LinearConstraintManager::RescaleActiveCounts(const double scaling_factor) {
76  for (ConstraintIndex i(0); i < constraint_infos_.size(); ++i) {
77  constraint_infos_[i].active_count *= scaling_factor;
78  }
79  constraint_active_count_increase_ *= scaling_factor;
80  VLOG(2) << "Rescaled active counts by " << scaling_factor;
81 }
82 
83 bool LinearConstraintManager::MaybeRemoveSomeInactiveConstraints(
84  glop::BasisState* solution_state) {
85  if (solution_state->IsEmpty()) return false; // Mainly to simplify tests.
86  const glop::RowIndex num_rows(lp_constraints_.size());
87  const glop::ColIndex num_cols =
88  solution_state->statuses.size() - RowToColIndex(num_rows);
89  int new_size = 0;
90  for (int i = 0; i < num_rows; ++i) {
91  const ConstraintIndex constraint_index = lp_constraints_[i];
92 
93  // Constraints that are not tight in the current solution have a basic
94  // status. We remove the ones that have been inactive in the last recent
95  // solves.
96  //
97  // TODO(user): More advanced heuristics might perform better, I didn't do
98  // a lot of tuning experiments yet.
99  const glop::VariableStatus row_status =
100  solution_state->statuses[num_cols + glop::ColIndex(i)];
101  if (row_status == glop::VariableStatus::BASIC) {
102  constraint_infos_[constraint_index].inactive_count++;
103  if (constraint_infos_[constraint_index].inactive_count >
104  sat_parameters_.max_consecutive_inactive_count()) {
105  constraint_infos_[constraint_index].is_in_lp = false;
106  continue; // Remove it.
107  }
108  } else {
109  // Only count consecutive inactivities.
110  constraint_infos_[constraint_index].inactive_count = 0;
111  }
112 
113  lp_constraints_[new_size] = constraint_index;
114  solution_state->statuses[num_cols + glop::ColIndex(new_size)] = row_status;
115  new_size++;
116  }
117  const int num_removed_constraints = lp_constraints_.size() - new_size;
118  lp_constraints_.resize(new_size);
119  solution_state->statuses.resize(num_cols + glop::ColIndex(new_size));
120  if (num_removed_constraints > 0) {
121  VLOG(2) << "Removed " << num_removed_constraints << " constraints";
122  }
123  return num_removed_constraints > 0;
124 }
125 
126 // Because sometimes we split a == constraint in two (>= and <=), it makes sense
127 // to detect duplicate constraints and merge bounds. This is also relevant if
128 // we regenerate identical cuts for some reason.
129 LinearConstraintManager::ConstraintIndex LinearConstraintManager::Add(
130  LinearConstraint ct, bool* added) {
131  CHECK(!ct.vars.empty());
133  SimplifyConstraint(&ct);
134  DivideByGCD(&ct);
137 
138  // If an identical constraint exists, only updates its bound.
139  const size_t key = ComputeHashOfTerms(ct);
140  if (gtl::ContainsKey(equiv_constraints_, key)) {
141  const ConstraintIndex ct_index = equiv_constraints_[key];
142  if (constraint_infos_[ct_index].constraint.vars == ct.vars &&
143  constraint_infos_[ct_index].constraint.coeffs == ct.coeffs) {
144  if (added != nullptr) *added = false;
145  if (ct.lb > constraint_infos_[ct_index].constraint.lb) {
146  if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
147  constraint_infos_[ct_index].constraint.lb = ct.lb;
148  if (added != nullptr) *added = true;
149  }
150  if (ct.ub < constraint_infos_[ct_index].constraint.ub) {
151  if (constraint_infos_[ct_index].is_in_lp) current_lp_is_changed_ = true;
152  constraint_infos_[ct_index].constraint.ub = ct.ub;
153  if (added != nullptr) *added = true;
154  }
155  ++num_merged_constraints_;
156  return ct_index;
157  }
158  }
159 
160  if (added != nullptr) *added = true;
161  const ConstraintIndex ct_index(constraint_infos_.size());
162  ConstraintInfo ct_info;
163  ct_info.constraint = std::move(ct);
164  ct_info.l2_norm = ComputeL2Norm(ct_info.constraint);
165  ct_info.hash = key;
166  equiv_constraints_[key] = ct_index;
167  ct_info.active_count = constraint_active_count_increase_;
168  constraint_infos_.push_back(std::move(ct_info));
169  return ct_index;
170 }
171 
172 void LinearConstraintManager::ComputeObjectiveParallelism(
173  const ConstraintIndex ct_index) {
174  CHECK(objective_is_defined_);
175  // lazy computation of objective norm.
176  if (!objective_norm_computed_) {
177  objective_l2_norm_ = std::sqrt(sum_of_squared_objective_coeffs_);
178  objective_norm_computed_ = true;
179  }
180  CHECK_GT(objective_l2_norm_, 0.0);
181 
182  constraint_infos_[ct_index].objective_parallelism_computed = true;
183  if (constraint_infos_[ct_index].l2_norm == 0.0) {
184  constraint_infos_[ct_index].objective_parallelism = 0.0;
185  return;
186  }
187 
188  const LinearConstraint& lc = constraint_infos_[ct_index].constraint;
189  double unscaled_objective_parallelism = 0.0;
190  for (int i = 0; i < lc.vars.size(); ++i) {
191  const IntegerVariable var = lc.vars[i];
192  const auto it = objective_map_.find(var);
193  if (it == objective_map_.end()) continue;
194  unscaled_objective_parallelism += it->second * ToDouble(lc.coeffs[i]);
195  }
196  const double objective_parallelism =
197  unscaled_objective_parallelism /
198  (constraint_infos_[ct_index].l2_norm * objective_l2_norm_);
199  constraint_infos_[ct_index].objective_parallelism =
200  std::abs(objective_parallelism);
201 }
202 
203 // Same as Add(), but logs some information about the newly added constraint.
204 // Cuts are also handled slightly differently than normal constraints.
206  LinearConstraint ct, std::string type_name,
208  std::string extra_info) {
209  ++num_add_cut_calls_;
210  if (ct.vars.empty()) return false;
211 
212  const double activity = ComputeActivity(ct, lp_solution);
213  const double violation =
214  std::max(activity - ToDouble(ct.ub), ToDouble(ct.lb) - activity);
215  const double l2_norm = ComputeL2Norm(ct);
216 
217  // Only add cut with sufficient efficacy.
218  if (violation / l2_norm < 1e-5) return false;
219 
220  bool added = false;
221  const ConstraintIndex ct_index = Add(std::move(ct), &added);
222 
223  // We only mark the constraint as a cut if it is not an update of an already
224  // existing one.
225  if (!added) return false;
226 
227  // TODO(user): Use better heuristic here for detecting good cuts and mark
228  // them undeletable.
229  constraint_infos_[ct_index].is_deletable = true;
230 
231  VLOG(1) << "Cut '" << type_name << "'"
232  << " size=" << constraint_infos_[ct_index].constraint.vars.size()
233  << " max_magnitude="
234  << ComputeInfinityNorm(constraint_infos_[ct_index].constraint)
235  << " norm=" << l2_norm << " violation=" << violation
236  << " eff=" << violation / l2_norm << " " << extra_info;
237 
238  num_cuts_++;
239  num_deletable_constraints_++;
240  type_to_num_cuts_[type_name]++;
241  return true;
242 }
243 
244 void LinearConstraintManager::PermanentlyRemoveSomeConstraints() {
245  std::vector<double> deletable_constraint_counts;
246  for (ConstraintIndex i(0); i < constraint_infos_.size(); ++i) {
247  if (constraint_infos_[i].is_deletable && !constraint_infos_[i].is_in_lp) {
248  deletable_constraint_counts.push_back(constraint_infos_[i].active_count);
249  }
250  }
251  if (deletable_constraint_counts.empty()) return;
252  std::sort(deletable_constraint_counts.begin(),
253  deletable_constraint_counts.end());
254 
255  // We will delete the oldest (in the order they where added) cleanup target
256  // constraints with a count lower or equal to this.
257  double active_count_threshold = std::numeric_limits<double>::infinity();
258  if (sat_parameters_.cut_cleanup_target() <
259  deletable_constraint_counts.size()) {
260  active_count_threshold =
261  deletable_constraint_counts[sat_parameters_.cut_cleanup_target()];
262  }
263 
264  ConstraintIndex new_size(0);
265  equiv_constraints_.clear();
267  constraint_infos_.size());
268  int num_deleted_constraints = 0;
269  for (ConstraintIndex i(0); i < constraint_infos_.size(); ++i) {
270  if (constraint_infos_[i].is_deletable && !constraint_infos_[i].is_in_lp &&
271  constraint_infos_[i].active_count <= active_count_threshold &&
272  num_deleted_constraints < sat_parameters_.cut_cleanup_target()) {
273  ++num_deleted_constraints;
274  continue;
275  }
276 
277  if (i != new_size) {
278  constraint_infos_[new_size] = std::move(constraint_infos_[i]);
279  }
280  index_mapping[i] = new_size;
281 
282  // Make sure we recompute the hash_map of identical constraints.
283  equiv_constraints_[constraint_infos_[new_size].hash] = new_size;
284  new_size++;
285  }
286  constraint_infos_.resize(new_size.value());
287 
288  // Also update lp_constraints_
289  for (int i = 0; i < lp_constraints_.size(); ++i) {
290  lp_constraints_[i] = index_mapping[lp_constraints_[i]];
291  }
292 
293  if (num_deleted_constraints > 0) {
294  VLOG(1) << "Constraint manager cleanup: #deleted:"
295  << num_deleted_constraints;
296  }
297  num_deletable_constraints_ -= num_deleted_constraints;
298 }
299 
301  IntegerValue coeff) {
302  if (coeff == IntegerValue(0)) return;
303  objective_is_defined_ = true;
304  if (!VariableIsPositive(var)) {
305  var = NegationOf(var);
306  coeff = -coeff;
307  }
308  const double coeff_as_double = ToDouble(coeff);
309  const auto insert = objective_map_.insert({var, coeff_as_double});
310  CHECK(insert.second)
311  << "SetObjectiveCoefficient() called twice with same variable";
312  sum_of_squared_objective_coeffs_ += coeff_as_double * coeff_as_double;
313 }
314 
315 bool LinearConstraintManager::SimplifyConstraint(LinearConstraint* ct) {
316  bool term_changed = false;
317 
318  IntegerValue min_sum(0);
319  IntegerValue max_sum(0);
320  IntegerValue max_magnitude(0);
321  int new_size = 0;
322  const int num_terms = ct->vars.size();
323  for (int i = 0; i < num_terms; ++i) {
324  const IntegerVariable var = ct->vars[i];
325  const IntegerValue coeff = ct->coeffs[i];
326  const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var);
327  const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var);
328 
329  // For now we do not change ct, but just compute its new_size if we where
330  // to remove a fixed term.
331  if (lb == ub) continue;
332  ++new_size;
333 
334  max_magnitude = std::max(max_magnitude, IntTypeAbs(coeff));
335  if (coeff > 0.0) {
336  min_sum += coeff * lb;
337  max_sum += coeff * ub;
338  } else {
339  min_sum += coeff * ub;
340  max_sum += coeff * lb;
341  }
342  }
343 
344  // Shorten the constraint if needed.
345  if (new_size < num_terms) {
346  term_changed = true;
347  ++num_shortened_constraints_;
348  new_size = 0;
349  for (int i = 0; i < num_terms; ++i) {
350  const IntegerVariable var = ct->vars[i];
351  const IntegerValue coeff = ct->coeffs[i];
352  const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var);
353  const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var);
354  if (lb == ub) {
355  const IntegerValue rhs_adjust = lb * coeff;
356  if (ct->lb > kMinIntegerValue) ct->lb -= rhs_adjust;
357  if (ct->ub < kMaxIntegerValue) ct->ub -= rhs_adjust;
358  continue;
359  }
360  ct->vars[new_size] = var;
361  ct->coeffs[new_size] = coeff;
362  ++new_size;
363  }
364  ct->vars.resize(new_size);
365  ct->coeffs.resize(new_size);
366  }
367 
368  // Relax the bound if needed, note that this doesn't require a change to
369  // the equiv map.
370  if (min_sum >= ct->lb) ct->lb = kMinIntegerValue;
371  if (max_sum <= ct->ub) ct->ub = kMaxIntegerValue;
372 
373  // Clear constraints that are always true.
374  // We rely on the deletion code to remove them eventually.
375  if (ct->lb == kMinIntegerValue && ct->ub == kMaxIntegerValue) {
376  ct->vars.clear();
377  ct->coeffs.clear();
378  return true;
379  }
380 
381  // TODO(user): Split constraint in two if it is boxed and there is possible
382  // reduction?
383  //
384  // TODO(user): Make sure there cannot be any overflow. They shouldn't, but
385  // I am not sure all the generated cuts are safe regarding min/max sum
386  // computation. We should check this.
387  if (ct->ub != kMaxIntegerValue && max_magnitude > max_sum - ct->ub) {
388  if (ct->lb != kMinIntegerValue) {
389  ++num_splitted_constraints_;
390  } else {
391  term_changed = true;
392  ++num_coeff_strenghtening_;
393  const int num_terms = ct->vars.size();
394  const IntegerValue target = max_sum - ct->ub;
395  for (int i = 0; i < num_terms; ++i) {
396  const IntegerValue coeff = ct->coeffs[i];
397  if (coeff > target) {
398  const IntegerVariable var = ct->vars[i];
399  const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var);
400  ct->coeffs[i] = target;
401  ct->ub -= (coeff - target) * ub;
402  } else if (coeff < -target) {
403  const IntegerVariable var = ct->vars[i];
404  const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var);
405  ct->coeffs[i] = -target;
406  ct->ub += (-target - coeff) * lb;
407  }
408  }
409  }
410  }
411 
412  if (ct->lb != kMinIntegerValue && max_magnitude > ct->lb - min_sum) {
413  if (ct->ub != kMaxIntegerValue) {
414  ++num_splitted_constraints_;
415  } else {
416  term_changed = true;
417  ++num_coeff_strenghtening_;
418  const int num_terms = ct->vars.size();
419  const IntegerValue target = ct->lb - min_sum;
420  for (int i = 0; i < num_terms; ++i) {
421  const IntegerValue coeff = ct->coeffs[i];
422  if (coeff > target) {
423  const IntegerVariable var = ct->vars[i];
424  const IntegerValue lb = integer_trail_.LevelZeroLowerBound(var);
425  ct->coeffs[i] = target;
426  ct->lb -= (coeff - target) * lb;
427  } else if (coeff < -target) {
428  const IntegerVariable var = ct->vars[i];
429  const IntegerValue ub = integer_trail_.LevelZeroUpperBound(var);
430  ct->coeffs[i] = -target;
431  ct->lb += (-target - coeff) * ub;
432  }
433  }
434  }
435  }
436 
437  return term_changed;
438 }
439 
442  glop::BasisState* solution_state) {
443  VLOG(3) << "Enter ChangeLP, scan " << constraint_infos_.size()
444  << " constraints";
445  const double saved_dtime = dtime_;
446  std::vector<ConstraintIndex> new_constraints;
447  std::vector<double> new_constraints_efficacies;
448  std::vector<double> new_constraints_orthogonalities;
449 
450  const bool simplify_constraints =
451  integer_trail_.num_level_zero_enqueues() > last_simplification_timestamp_;
452  last_simplification_timestamp_ = integer_trail_.num_level_zero_enqueues();
453 
454  // We keep any constraints that is already present, and otherwise, we add the
455  // ones that are currently not satisfied by at least "tolerance" to the set
456  // of potential new constraints.
457  bool rescale_active_count = false;
458  const double tolerance = 1e-6;
459  for (ConstraintIndex i(0); i < constraint_infos_.size(); ++i) {
460  // Inprocessing of the constraint.
461  if (simplify_constraints &&
462  SimplifyConstraint(&constraint_infos_[i].constraint)) {
463  ++num_simplifications_;
464 
465  // Note that the canonicalization shouldn't be needed since the order
466  // of the variable is not changed by the simplification, and we only
467  // reduce the coefficients at both end of the spectrum.
468  DivideByGCD(&constraint_infos_[i].constraint);
469  DCHECK(DebugCheckConstraint(constraint_infos_[i].constraint));
470 
471  constraint_infos_[i].objective_parallelism_computed = false;
472  constraint_infos_[i].l2_norm =
473  ComputeL2Norm(constraint_infos_[i].constraint);
474 
475  if (constraint_infos_[i].is_in_lp) current_lp_is_changed_ = true;
476  equiv_constraints_.erase(constraint_infos_[i].hash);
477  constraint_infos_[i].hash =
478  ComputeHashOfTerms(constraint_infos_[i].constraint);
479 
480  // TODO(user): Because we simplified this constraint, it is possible that
481  // it is now a duplicate of another one. Merge them.
482  equiv_constraints_[constraint_infos_[i].hash] = i;
483  }
484 
485  if (constraint_infos_[i].is_in_lp) continue;
486 
487  // ComputeActivity() often represent the bulk of the time spent in
488  // ChangeLP().
489  dtime_ += 1.7e-9 *
490  static_cast<double>(constraint_infos_[i].constraint.vars.size());
491  const double activity =
492  ComputeActivity(constraint_infos_[i].constraint, lp_solution);
493  const double lb_violation =
494  ToDouble(constraint_infos_[i].constraint.lb) - activity;
495  const double ub_violation =
496  activity - ToDouble(constraint_infos_[i].constraint.ub);
497  const double violation = std::max(lb_violation, ub_violation);
498  if (violation >= tolerance) {
499  constraint_infos_[i].inactive_count = 0;
500  new_constraints.push_back(i);
501  new_constraints_efficacies.push_back(violation /
502  constraint_infos_[i].l2_norm);
503  new_constraints_orthogonalities.push_back(1.0);
504 
505  if (objective_is_defined_ &&
506  !constraint_infos_[i].objective_parallelism_computed) {
507  ComputeObjectiveParallelism(i);
508  } else if (!objective_is_defined_) {
509  constraint_infos_[i].objective_parallelism = 0.0;
510  }
511 
512  constraint_infos_[i].current_score =
513  new_constraints_efficacies.back() +
514  constraint_infos_[i].objective_parallelism;
515 
516  if (constraint_infos_[i].is_deletable) {
517  constraint_infos_[i].active_count += constraint_active_count_increase_;
518  if (constraint_infos_[i].active_count >
519  sat_parameters_.cut_max_active_count_value()) {
520  rescale_active_count = true;
521  }
522  }
523  }
524  }
525 
526  // Bump activities of active constraints in LP.
527  if (solution_state != nullptr) {
528  const glop::RowIndex num_rows(lp_constraints_.size());
529  const glop::ColIndex num_cols =
530  solution_state->statuses.size() - RowToColIndex(num_rows);
531 
532  for (int i = 0; i < num_rows; ++i) {
533  const ConstraintIndex constraint_index = lp_constraints_[i];
534  const glop::VariableStatus row_status =
535  solution_state->statuses[num_cols + glop::ColIndex(i)];
536  if (row_status != glop::VariableStatus::BASIC &&
537  constraint_infos_[constraint_index].is_deletable) {
538  constraint_infos_[constraint_index].active_count +=
539  constraint_active_count_increase_;
540  if (constraint_infos_[constraint_index].active_count >
541  sat_parameters_.cut_max_active_count_value()) {
542  rescale_active_count = true;
543  }
544  }
545  }
546  }
547 
548  if (rescale_active_count) {
549  CHECK_GT(sat_parameters_.cut_max_active_count_value(), 0.0);
550  RescaleActiveCounts(1.0 / sat_parameters_.cut_max_active_count_value());
551  }
552 
553  // Update the increment counter.
554  constraint_active_count_increase_ *=
555  1.0 / sat_parameters_.cut_active_count_decay();
556 
557  // Remove constraints from the current LP that have been inactive for a while.
558  // We do that after we computed new_constraints so we do not need to iterate
559  // over the just deleted constraints.
560  if (MaybeRemoveSomeInactiveConstraints(solution_state)) {
561  current_lp_is_changed_ = true;
562  }
563 
564  // Note that the algo below is in O(limit * new_constraint). In order to
565  // limit spending too much time on this, we first sort all the constraints
566  // with an imprecise score (no orthogonality), then limit the size of the
567  // vector of constraints to precisely score, then we do the actual scoring.
568  //
569  // On problem crossword_opt_grid-19.05_dict-80_sat with linearization_level=2,
570  // new_constraint.size() > 1.5M.
571  //
572  // TODO(user): This blowup factor could be adaptative w.r.t. the constraint
573  // limit.
574  const int kBlowupFactor = 4;
575  int constraint_limit = std::min(sat_parameters_.new_constraints_batch_size(),
576  static_cast<int>(new_constraints.size()));
577  if (lp_constraints_.empty()) {
578  constraint_limit = std::min(1000, static_cast<int>(new_constraints.size()));
579  }
580  VLOG(3) << " - size = " << new_constraints.size()
581  << ", limit = " << constraint_limit;
582 
583  std::stable_sort(new_constraints.begin(), new_constraints.end(),
584  [&](ConstraintIndex a, ConstraintIndex b) {
585  return constraint_infos_[a].current_score >
586  constraint_infos_[b].current_score;
587  });
588  if (new_constraints.size() > kBlowupFactor * constraint_limit) {
589  VLOG(3) << "Resize candidate constraints from " << new_constraints.size()
590  << " down to " << kBlowupFactor * constraint_limit;
591  new_constraints.resize(kBlowupFactor * constraint_limit);
592  }
593 
594  int num_added = 0;
595  int num_skipped_checks = 0;
596  const int kCheckFrequency = 100;
597  ConstraintIndex last_added_candidate = kInvalidConstraintIndex;
598  for (int i = 0; i < constraint_limit; ++i) {
599  // Iterate through all new constraints and select the one with the best
600  // score.
601  double best_score = 0.0;
602  ConstraintIndex best_candidate = kInvalidConstraintIndex;
603  for (int j = 0; j < new_constraints.size(); ++j) {
604  // Checks the time limit, and returns if the lp has changed.
605  if (++num_skipped_checks >= kCheckFrequency) {
606  if (time_limit_->LimitReached()) return current_lp_is_changed_;
607  num_skipped_checks = 0;
608  }
609 
610  const ConstraintIndex new_constraint = new_constraints[j];
611  if (constraint_infos_[new_constraint].is_in_lp) continue;
612 
613  if (last_added_candidate != kInvalidConstraintIndex) {
614  const double current_orthogonality =
615  1.0 - (std::abs(ScalarProduct(
616  constraint_infos_[last_added_candidate].constraint,
617  constraint_infos_[new_constraint].constraint)) /
618  (constraint_infos_[last_added_candidate].l2_norm *
619  constraint_infos_[new_constraint].l2_norm));
620  new_constraints_orthogonalities[j] =
621  std::min(new_constraints_orthogonalities[j], current_orthogonality);
622  }
623 
624  // NOTE(user): It is safe to not add this constraint as the constraint
625  // that is almost parallel to this constraint is present in the LP or is
626  // inactive for a long time and is removed from the LP. In either case,
627  // this constraint is not adding significant value and is only making the
628  // LP larger.
629  if (new_constraints_orthogonalities[j] <
630  sat_parameters_.min_orthogonality_for_lp_constraints()) {
631  continue;
632  }
633 
634  // TODO(user): Experiment with different weights or different
635  // functions for computing score.
636  const double score = new_constraints_orthogonalities[j] +
637  constraint_infos_[new_constraint].current_score;
638  CHECK_GE(score, 0.0);
639  if (score > best_score || best_candidate == kInvalidConstraintIndex) {
640  best_score = score;
641  best_candidate = new_constraint;
642  }
643  }
644 
645  if (best_candidate != kInvalidConstraintIndex) {
646  // Add the best constraint in the LP.
647  constraint_infos_[best_candidate].is_in_lp = true;
648  // Note that it is important for LP incremental solving that the old
649  // constraints stays at the same position in this list (and thus in the
650  // returned GetLp()).
651  ++num_added;
652  current_lp_is_changed_ = true;
653  lp_constraints_.push_back(best_candidate);
654  last_added_candidate = best_candidate;
655  }
656  }
657 
658  if (num_added > 0) {
659  // We update the solution sate to match the new LP size.
660  VLOG(2) << "Added " << num_added << " constraints.";
661  solution_state->statuses.resize(solution_state->statuses.size() + num_added,
663  }
664 
665  // TODO(user): Instead of comparing num_deletable_constraints with cut
666  // limit, compare number of deletable constraints not in lp against the limit.
667  if (num_deletable_constraints_ > sat_parameters_.max_num_cuts()) {
668  PermanentlyRemoveSomeConstraints();
669  }
670 
671  time_limit_->AdvanceDeterministicTime(dtime_ - saved_dtime);
672 
673  // The LP changed only if we added new constraints or if some constraints
674  // already inside changed (simplification or tighter bounds).
675  if (current_lp_is_changed_) {
676  current_lp_is_changed_ = false;
677  return true;
678  }
679  return false;
680 }
681 
683  for (ConstraintIndex i(0); i < constraint_infos_.size(); ++i) {
684  if (constraint_infos_[i].is_in_lp) continue;
685  constraint_infos_[i].is_in_lp = true;
686  lp_constraints_.push_back(i);
687  }
688 }
689 
691  const LinearConstraint& cut) {
692  if (model_->Get<DebugSolution>() == nullptr) return true;
693  const auto& debug_solution = *(model_->Get<DebugSolution>());
694  if (debug_solution.empty()) return true;
695 
696  IntegerValue activity(0);
697  for (int i = 0; i < cut.vars.size(); ++i) {
698  const IntegerVariable var = cut.vars[i];
699  const IntegerValue coeff = cut.coeffs[i];
700  activity += coeff * debug_solution[var];
701  }
702  if (activity > cut.ub || activity < cut.lb) {
703  LOG(INFO) << "activity " << activity << " not in [" << cut.lb << ","
704  << cut.ub << "]";
705  return false;
706  }
707  return true;
708 }
709 
711  LinearConstraint ct, const std::string& name,
712  const absl::StrongVector<IntegerVariable, double>& lp_solution) {
713  if (ct.vars.empty()) return;
714  const double activity = ComputeActivity(ct, lp_solution);
715  const double violation =
716  std::max(activity - ToDouble(ct.ub), ToDouble(ct.lb) - activity);
717  const double l2_norm = ComputeL2Norm(ct);
718  cuts_.Add({name, ct}, violation / l2_norm);
719 }
720 
723  LinearConstraintManager* manager) {
724  for (const CutCandidate& candidate : cuts_.UnorderedElements()) {
725  manager->AddCut(candidate.cut, candidate.name, lp_solution);
726  }
727  cuts_.Clear();
728 }
729 
730 } // namespace sat
731 } // 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 CHECK_GE(val1, val2)
Definition: base/logging.h:701
#define CHECK_GT(val1, val2)
Definition: base/logging.h:702
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define VLOG(verboselevel)
Definition: base/logging.h:978
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 AdvanceDeterministicTime(double deterministic_duration)
Advances the deterministic time.
Definition: time_limit.h:226
IntegerValue LevelZeroUpperBound(IntegerVariable var) const
Definition: integer.h:1355
IntegerValue LevelZeroLowerBound(IntegerVariable var) const
Definition: integer.h:1350
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)
bool AddCut(LinearConstraint ct, std::string type_name, const absl::StrongVector< IntegerVariable, double > &lp_solution, std::string extra_info="")
const std::string & Name() const
Definition: sat/model.h:175
T Get(std::function< T(const Model &)> f) const
Similar to Add() but this is const.
Definition: sat/model.h:87
void AddCut(LinearConstraint ct, const std::string &name, const absl::StrongVector< IntegerVariable, double > &lp_solution)
void TransferToManager(const absl::StrongVector< IntegerVariable, double > &lp_solution, LinearConstraintManager *manager)
const std::vector< Element > & UnorderedElements() const
void Add(Element e, double score)
const std::string name
const Constraint * ct
IntVar * var
Definition: expr_array.cc:1858
const int INFO
Definition: log_severity.h:31
int64 hash
Definition: matrix_utils.cc:60
bool ContainsKey(const Collection &collection, const Key &key)
Definition: map_util.h:170
ColIndex RowToColIndex(RowIndex row)
Definition: lp_types.h:48
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
constexpr IntegerValue kMinIntegerValue(-kMaxIntegerValue)
double ScalarProduct(const LinearConstraint &constraint1, const LinearConstraint &constraint2)
void CanonicalizeConstraint(LinearConstraint *ct)
bool NoDuplicateVariable(const LinearConstraint &ct)
double ComputeL2Norm(const LinearConstraint &constraint)
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)
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...
uint64 Hash(uint64 num, uint64 c)
Definition: hash.h:150
std::vector< IntegerVariable > vars