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"
33#include "ortools/glop/parameters.pb.h"
35#include "ortools/glop/status.h"
39#include "ortools/sat/integer.h"
42
43namespace operations_research {
44namespace sat {
45
46using glop::ColIndex;
48using 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
67bool 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
132std::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_(
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
205glop::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.
255bool 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
377LPSolveInfo 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
400void 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
419bool 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
618glop::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
634void 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
645bool 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
692bool 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
842bool 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.
959void 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
1041namespace {
1042
1043// For each element of a, adds a random one in b and append the pair to output.
1044void 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
1056template <class ListOfTerms>
1057IntegerValue 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
1066void 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;
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 }
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
1304void 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
1341void 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;
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.
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.
1638IntegerValue 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
1653bool 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
1675namespace {
1676
1677absl::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
1686void 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().
1749void 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
1773std::vector<std::pair<RowIndex, IntegerValue>>
1774LinearProgrammingConstraint::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
1839bool 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.
1869void 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.
2028bool 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
2096bool 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
2139int64 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;
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
2155void 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
2198namespace {
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.
2205void 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
2517namespace {
2518
2519// Returns for each literal its integer view, or the view of its negation.
2520std::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,
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
2576std::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
2622std::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
2700void 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().
2777std::function<IntegerLiteral()>
2779 return [this]() { return this->LPReducedCostAverageDecision(); };
2780}
2781
2782IntegerLiteral 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) {
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) {
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)
const ScatteredRow & GetUnitRowLeftInverse(RowIndex row)
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)
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 std::vector< ConstraintIndex > & LpConstraints() const
bool AddCut(LinearConstraint ct, std::string type_name, const absl::StrongVector< IntegerVariable, double > &lp_solution, std::string extra_info="")
const absl::StrongVector< ConstraintIndex, ConstraintInfo > & AllConstraints() const
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)
bool AddLinearExpressionMultiple(IntegerValue multiplier, const std::vector< std::pair< glop::ColIndex, IntegerValue > > &terms)
std::vector< std::pair< glop::ColIndex, IntegerValue > > GetTerms()
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)
void AddOneConstraint(glop::RowIndex, const std::vector< std::pair< glop::ColIndex, IntegerValue > > &terms, IntegerValue lb, IntegerValue ub)
std::vector< std::vector< std::pair< glop::RowIndex, IntegerValue > > > InterestingCandidates(ModelRandomGenerator *random)
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
RowIndex ColToRowIndex(ColIndex col)
Definition: lp_types.h:51
StrictITIVector< ColIndex, Fractional > DenseRow
Definition: lp_types.h:299
ColIndex RowToColIndex(RowIndex row)
Definition: lp_types.h:48
const double kEpsilon
Definition: lp_types.h:86
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
bool VariableIsPositive(IntegerVariable i)
Definition: integer.h:130
bool AddProductTo(IntegerValue a, IntegerValue b, IntegerValue *result)
Definition: integer.h:110
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)
IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor)
Definition: integer.h:90
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
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)
IntegerValue ComputeInfinityNorm(const LinearConstraint &constraint)
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)
constexpr IntegerValue kMaxIntegerValue(std::numeric_limits< IntegerValue::ValueType >::max() - 1)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int64 CapAdd(int64 x, int64 y)
int64 CapProd(int64 x, int64 y)
int64 CapSub(int64 x, int64 y)
std::pair< int64, int64 > Arc
Definition: search.cc:3361
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