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
26namespace operations_research {
27namespace sat {
28
29namespace {
30
31const LinearConstraintManager::ConstraintIndex kInvalidConstraintIndex(-1);
32
33size_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
75void 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
83bool 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.
129LinearConstraintManager::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
172void 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
244void 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
315bool 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,
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="")
T Get(std::function< T(const Model &)> f) const
Similar to Add() but this is const.
Definition: sat/model.h:87
const std::string & Name() const
Definition: sat/model.h:175
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
bool VariableIsPositive(IntegerVariable i)
Definition: integer.h:130
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)
IntType IntTypeAbs(IntType t)
Definition: integer.h:77
double ToDouble(IntegerValue value)
Definition: integer.h:69
std::vector< IntegerVariable > NegationOf(const std::vector< IntegerVariable > &vars)
Definition: integer.cc:27
IntegerValue ComputeInfinityNorm(const LinearConstraint &constraint)
void DivideByGCD(LinearConstraint *constraint)
double ComputeActivity(const LinearConstraint &constraint, const absl::StrongVector< IntegerVariable, double > &values)
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
uint64 Hash(uint64 num, uint64 c)
Definition: hash.h:150
std::vector< IntegerVariable > vars