1 // Copyright 2010-2021 Google LLC 2 // Licensed under the Apache License, Version 2.0 (the "License"); 3 // you may not use this file except in compliance with the License. 4 // You may obtain a copy of the License at 5 // 6 // http://www.apache.org/licenses/LICENSE-2.0 7 // 8 // Unless required by applicable law or agreed to in writing, software 9 // distributed under the License is distributed on an "AS IS" BASIS, 10 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 11 // See the License for the specific language governing permissions and 12 // limitations under the License. 13 14 #include "ortools/sat/cuts.h" 15 16 #include <algorithm> 17 #include <cmath> 18 #include <cstdint> 19 #include <functional> 20 #include <limits> 21 #include <memory> 22 #include <string> 23 #include <utility> 24 #include <vector> 25 26 #include "absl/container/btree_set.h" 27 #include "ortools/algorithms/knapsack_solver_for_cuts.h" 28 #include "ortools/base/integral_types.h" 29 #include "ortools/base/stl_util.h" 30 #include "ortools/base/strong_vector.h" 31 #include "ortools/sat/integer.h" 32 #include "ortools/sat/linear_constraint.h" 33 #include "ortools/sat/linear_constraint_manager.h" 34 #include "ortools/sat/sat_base.h" 35 #include "ortools/sat/util.h" 36 #include "ortools/util/time_limit.h" 37 38 namespace operations_research { 39 namespace sat { 40 41 namespace { 42 43 // Minimum amount of violation of the cut constraint by the solution. This 44 // is needed to avoid numerical issues and adding cuts with minor effect. 45 const double kMinCutViolation = 1e-4; 46 47 // Returns the lp value of a Literal. 48 double GetLiteralLpValue( 49 const Literal lit, 50 const absl::StrongVector<IntegerVariable, double>& lp_values, 51 const IntegerEncoder* encoder) { 52 const IntegerVariable direct_view = encoder->GetLiteralView(lit); 53 if (direct_view != kNoIntegerVariable) { 54 return lp_values[direct_view]; 55 } 56 const IntegerVariable opposite_view = encoder->GetLiteralView(lit.Negated()); 57 DCHECK_NE(opposite_view, kNoIntegerVariable); 58 return 1.0 - lp_values[opposite_view]; 59 } 60 61 // Returns a constraint that disallow all given variables to be at their current 62 // upper bound. The arguments must form a non-trival constraint of the form 63 // sum terms (coeff * var) <= upper_bound. 64 LinearConstraint GenerateKnapsackCutForCover( 65 const std::vector<IntegerVariable>& vars, 66 const std::vector<IntegerValue>& coeffs, const IntegerValue upper_bound, 67 const IntegerTrail& integer_trail) { 68 CHECK_EQ(vars.size(), coeffs.size()); 69 CHECK_GT(vars.size(), 0); 70 LinearConstraint cut; 71 IntegerValue cut_upper_bound = IntegerValue(0); 72 IntegerValue max_coeff = coeffs[0]; 73 // slack = \sum_{i}(coeffs[i] * upper_bound[i]) - upper_bound. 74 IntegerValue slack = -upper_bound; 75 for (int i = 0; i < vars.size(); ++i) { 76 const IntegerValue var_upper_bound = 77 integer_trail.LevelZeroUpperBound(vars[i]); 78 cut_upper_bound += var_upper_bound; 79 cut.vars.push_back(vars[i]); 80 cut.coeffs.push_back(IntegerValue(1)); 81 max_coeff = std::max(max_coeff, coeffs[i]); 82 slack += coeffs[i] * var_upper_bound; 83 } 84 CHECK_GT(slack, 0.0) << "Invalid cover for knapsack cut."; 85 cut_upper_bound -= CeilRatio(slack, max_coeff); 86 cut.lb = kMinIntegerValue; 87 cut.ub = cut_upper_bound; 88 VLOG(2) << "Generated Knapsack Constraint:" << cut.DebugString(); 89 return cut; 90 } 91 92 bool SolutionSatisfiesConstraint( 93 const LinearConstraint& constraint, 94 const absl::StrongVector<IntegerVariable, double>& lp_values) { 95 const double activity = ComputeActivity(constraint, lp_values); 96 const double tolerance = 1e-6; 97 return (activity <= ToDouble(constraint.ub) + tolerance && 98 activity >= ToDouble(constraint.lb) - tolerance) 99 ? true 100 : false; 101 } 102 103 bool SmallRangeAndAllCoefficientsMagnitudeAreTheSame( 104 const LinearConstraint& constraint, IntegerTrail* integer_trail) { 105 if (constraint.vars.empty()) return true; 106 107 const int64_t magnitude = std::abs(constraint.coeffs[0].value()); 108 for (int i = 1; i < constraint.coeffs.size(); ++i) { 109 const IntegerVariable var = constraint.vars[i]; 110 if (integer_trail->LevelZeroUpperBound(var) - 111 integer_trail->LevelZeroLowerBound(var) > 112 1) { 113 return false; 114 } 115 if (std::abs(constraint.coeffs[i].value()) != magnitude) { 116 return false; 117 } 118 } 119 return true; 120 } 121 122 bool AllVarsTakeIntegerValue( 123 const std::vector<IntegerVariable> vars, 124 const absl::StrongVector<IntegerVariable, double>& lp_values) { 125 for (IntegerVariable var : vars) { 126 if (std::abs(lp_values[var] - std::round(lp_values[var])) > 1e-6) { 127 return false; 128 } 129 } 130 return true; 131 } 132 133 // Returns smallest cover size for the given constraint taking into account 134 // level zero bounds. Smallest Cover size is computed as follows. 135 // 1. Compute the upper bound if all variables are shifted to have zero lower 136 // bound. 137 // 2. Sort all terms (coefficient * shifted upper bound) in non decreasing 138 // order. 139 // 3. Add terms in cover until term sum is smaller or equal to upper bound. 140 // 4. Add the last item which violates the upper bound. This forms the smallest 141 // cover. Return the size of this cover. 142 int GetSmallestCoverSize(const LinearConstraint& constraint, 143 const IntegerTrail& integer_trail) { 144 IntegerValue ub = constraint.ub; 145 std::vector<IntegerValue> sorted_terms; 146 for (int i = 0; i < constraint.vars.size(); ++i) { 147 const IntegerValue coeff = constraint.coeffs[i]; 148 const IntegerVariable var = constraint.vars[i]; 149 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var); 150 const IntegerValue var_lb = integer_trail.LevelZeroLowerBound(var); 151 ub -= var_lb * coeff; 152 sorted_terms.push_back(coeff * (var_ub - var_lb)); 153 } 154 std::sort(sorted_terms.begin(), sorted_terms.end(), 155 std::greater<IntegerValue>()); 156 int smallest_cover_size = 0; 157 IntegerValue sorted_term_sum = IntegerValue(0); 158 while (sorted_term_sum <= ub && 159 smallest_cover_size < constraint.vars.size()) { 160 sorted_term_sum += sorted_terms[smallest_cover_size++]; 161 } 162 return smallest_cover_size; 163 } 164 165 bool ConstraintIsEligibleForLifting(const LinearConstraint& constraint, 166 const IntegerTrail& integer_trail) { 167 for (const IntegerVariable var : constraint.vars) { 168 if (integer_trail.LevelZeroLowerBound(var) != IntegerValue(0) || 169 integer_trail.LevelZeroUpperBound(var) != IntegerValue(1)) { 170 return false; 171 } 172 } 173 return true; 174 } 175 } // namespace 176 177 bool LiftKnapsackCut( 178 const LinearConstraint& constraint, 179 const absl::StrongVector<IntegerVariable, double>& lp_values, 180 const std::vector<IntegerValue>& cut_vars_original_coefficients, 181 const IntegerTrail& integer_trail, TimeLimit* time_limit, 182 LinearConstraint* cut) { 183 std::set<IntegerVariable> vars_in_cut; 184 for (IntegerVariable var : cut->vars) { 185 vars_in_cut.insert(var); 186 } 187 188 std::vector<std::pair<IntegerValue, IntegerVariable>> non_zero_vars; 189 std::vector<std::pair<IntegerValue, IntegerVariable>> zero_vars; 190 for (int i = 0; i < constraint.vars.size(); ++i) { 191 const IntegerVariable var = constraint.vars[i]; 192 if (integer_trail.LevelZeroLowerBound(var) != IntegerValue(0) || 193 integer_trail.LevelZeroUpperBound(var) != IntegerValue(1)) { 194 continue; 195 } 196 if (vars_in_cut.find(var) != vars_in_cut.end()) continue; 197 const IntegerValue coeff = constraint.coeffs[i]; 198 if (lp_values[var] <= 1e-6) { 199 zero_vars.push_back({coeff, var}); 200 } else { 201 non_zero_vars.push_back({coeff, var}); 202 } 203 } 204 205 // Decide lifting sequence (nonzeros, zeros in nonincreasing order 206 // of coefficient ). 207 std::sort(non_zero_vars.rbegin(), non_zero_vars.rend()); 208 std::sort(zero_vars.rbegin(), zero_vars.rend()); 209 210 std::vector<std::pair<IntegerValue, IntegerVariable>> lifting_sequence( 211 std::move(non_zero_vars)); 212 213 lifting_sequence.insert(lifting_sequence.end(), zero_vars.begin(), 214 zero_vars.end()); 215 216 // Form Knapsack. 217 std::vector<double> lifting_profits; 218 std::vector<double> lifting_weights; 219 for (int i = 0; i < cut->vars.size(); ++i) { 220 lifting_profits.push_back(ToDouble(cut->coeffs[i])); 221 lifting_weights.push_back(ToDouble(cut_vars_original_coefficients[i])); 222 } 223 224 // Lift the cut. 225 bool is_lifted = false; 226 bool is_solution_optimal = false; 227 KnapsackSolverForCuts knapsack_solver("Knapsack cut lifter"); 228 for (auto entry : lifting_sequence) { 229 is_solution_optimal = false; 230 const IntegerValue var_original_coeff = entry.first; 231 const IntegerVariable var = entry.second; 232 const IntegerValue lifting_capacity = constraint.ub - entry.first; 233 if (lifting_capacity <= IntegerValue(0)) continue; 234 knapsack_solver.Init(lifting_profits, lifting_weights, 235 ToDouble(lifting_capacity)); 236 knapsack_solver.set_node_limit(100); 237 // NOTE: Since all profits and weights are integer, solution of 238 // knapsack is also integer. 239 // TODO(user): Use an integer solver or heuristic. 240 knapsack_solver.Solve(time_limit, &is_solution_optimal); 241 const double knapsack_upper_bound = 242 std::round(knapsack_solver.GetUpperBound()); 243 const IntegerValue cut_coeff = 244 cut->ub - static_cast<int64_t>(knapsack_upper_bound); 245 if (cut_coeff > IntegerValue(0)) { 246 is_lifted = true; 247 cut->vars.push_back(var); 248 cut->coeffs.push_back(cut_coeff); 249 lifting_profits.push_back(ToDouble(cut_coeff)); 250 lifting_weights.push_back(ToDouble(var_original_coeff)); 251 } 252 } 253 return is_lifted; 254 } 255 256 LinearConstraint GetPreprocessedLinearConstraint( 257 const LinearConstraint& constraint, 258 const absl::StrongVector<IntegerVariable, double>& lp_values, 259 const IntegerTrail& integer_trail) { 260 IntegerValue ub = constraint.ub; 261 LinearConstraint constraint_with_left_vars; 262 for (int i = 0; i < constraint.vars.size(); ++i) { 263 const IntegerVariable var = constraint.vars[i]; 264 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var); 265 const IntegerValue coeff = constraint.coeffs[i]; 266 if (ToDouble(var_ub) - lp_values[var] <= 1.0 - kMinCutViolation) { 267 constraint_with_left_vars.vars.push_back(var); 268 constraint_with_left_vars.coeffs.push_back(coeff); 269 } else { 270 // Variable not in cut 271 const IntegerValue var_lb = integer_trail.LevelZeroLowerBound(var); 272 ub -= coeff * var_lb; 273 } 274 } 275 constraint_with_left_vars.ub = ub; 276 constraint_with_left_vars.lb = constraint.lb; 277 return constraint_with_left_vars; 278 } 279 280 bool ConstraintIsTriviallyTrue(const LinearConstraint& constraint, 281 const IntegerTrail& integer_trail) { 282 IntegerValue term_sum = IntegerValue(0); 283 for (int i = 0; i < constraint.vars.size(); ++i) { 284 const IntegerVariable var = constraint.vars[i]; 285 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var); 286 const IntegerValue coeff = constraint.coeffs[i]; 287 term_sum += coeff * var_ub; 288 } 289 if (term_sum <= constraint.ub) { 290 VLOG(2) << "Filtered by cover filter"; 291 return true; 292 } 293 return false; 294 } 295 296 bool CanBeFilteredUsingCutLowerBound( 297 const LinearConstraint& preprocessed_constraint, 298 const absl::StrongVector<IntegerVariable, double>& lp_values, 299 const IntegerTrail& integer_trail) { 300 std::vector<double> variable_upper_bound_distances; 301 for (const IntegerVariable var : preprocessed_constraint.vars) { 302 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var); 303 variable_upper_bound_distances.push_back(ToDouble(var_ub) - lp_values[var]); 304 } 305 // Compute the min cover size. 306 const int smallest_cover_size = 307 GetSmallestCoverSize(preprocessed_constraint, integer_trail); 308 309 std::nth_element( 310 variable_upper_bound_distances.begin(), 311 variable_upper_bound_distances.begin() + smallest_cover_size - 1, 312 variable_upper_bound_distances.end()); 313 double cut_lower_bound = 0.0; 314 for (int i = 0; i < smallest_cover_size; ++i) { 315 cut_lower_bound += variable_upper_bound_distances[i]; 316 } 317 if (cut_lower_bound >= 1.0 - kMinCutViolation) { 318 VLOG(2) << "Filtered by kappa heuristic"; 319 return true; 320 } 321 return false; 322 } 323 324 double GetKnapsackUpperBound(std::vector<KnapsackItem> items, 325 const double capacity) { 326 // Sort items by value by weight ratio. 327 std::sort(items.begin(), items.end(), std::greater<KnapsackItem>()); 328 double left_capacity = capacity; 329 double profit = 0.0; 330 for (const KnapsackItem item : items) { 331 if (item.weight <= left_capacity) { 332 profit += item.profit; 333 left_capacity -= item.weight; 334 } else { 335 profit += (left_capacity / item.weight) * item.profit; 336 break; 337 } 338 } 339 return profit; 340 } 341 342 bool CanBeFilteredUsingKnapsackUpperBound( 343 const LinearConstraint& constraint, 344 const absl::StrongVector<IntegerVariable, double>& lp_values, 345 const IntegerTrail& integer_trail) { 346 std::vector<KnapsackItem> items; 347 double capacity = -ToDouble(constraint.ub) - 1.0; 348 double sum_variable_profit = 0; 349 for (int i = 0; i < constraint.vars.size(); ++i) { 350 const IntegerVariable var = constraint.vars[i]; 351 const IntegerValue var_ub = integer_trail.LevelZeroUpperBound(var); 352 const IntegerValue var_lb = integer_trail.LevelZeroLowerBound(var); 353 const IntegerValue coeff = constraint.coeffs[i]; 354 KnapsackItem item; 355 item.profit = ToDouble(var_ub) - lp_values[var]; 356 item.weight = ToDouble(coeff * (var_ub - var_lb)); 357 items.push_back(item); 358 capacity += ToDouble(coeff * var_ub); 359 sum_variable_profit += item.profit; 360 } 361 362 // Return early if the required upper bound is negative since all the profits 363 // are non negative. 364 if (sum_variable_profit - 1.0 + kMinCutViolation < 0.0) return false; 365 366 // Get the knapsack upper bound. 367 const double knapsack_upper_bound = 368 GetKnapsackUpperBound(std::move(items), capacity); 369 if (knapsack_upper_bound < sum_variable_profit - 1.0 + kMinCutViolation) { 370 VLOG(2) << "Filtered by knapsack upper bound"; 371 return true; 372 } 373 return false; 374 } 375 376 bool CanFormValidKnapsackCover( 377 const LinearConstraint& preprocessed_constraint, 378 const absl::StrongVector<IntegerVariable, double>& lp_values, 379 const IntegerTrail& integer_trail) { 380 if (ConstraintIsTriviallyTrue(preprocessed_constraint, integer_trail)) { 381 return false; 382 } 383 if (CanBeFilteredUsingCutLowerBound(preprocessed_constraint, lp_values, 384 integer_trail)) { 385 return false; 386 } 387 if (CanBeFilteredUsingKnapsackUpperBound(preprocessed_constraint, lp_values, 388 integer_trail)) { 389 return false; 390 } 391 return true; 392 } 393 394 void ConvertToKnapsackForm(const LinearConstraint& constraint, 395 std::vector<LinearConstraint>* knapsack_constraints, 396 IntegerTrail* integer_trail) { 397 // If all coefficient are the same, the generated knapsack cuts cannot be 398 // stronger than the constraint itself. However, when we substitute variables 399 // using the implication graph, this is not longer true. So we only skip 400 // constraints with same coeff and no substitutions. 401 if (SmallRangeAndAllCoefficientsMagnitudeAreTheSame(constraint, 402 integer_trail)) { 403 return; 404 } 405 if (constraint.ub < kMaxIntegerValue) { 406 LinearConstraint canonical_knapsack_form; 407 408 // Negate the variables with negative coefficients. 409 for (int i = 0; i < constraint.vars.size(); ++i) { 410 const IntegerVariable var = constraint.vars[i]; 411 const IntegerValue coeff = constraint.coeffs[i]; 412 if (coeff > IntegerValue(0)) { 413 canonical_knapsack_form.AddTerm(var, coeff); 414 } else { 415 canonical_knapsack_form.AddTerm(NegationOf(var), -coeff); 416 } 417 } 418 canonical_knapsack_form.ub = constraint.ub; 419 canonical_knapsack_form.lb = kMinIntegerValue; 420 knapsack_constraints->push_back(canonical_knapsack_form); 421 } 422 423 if (constraint.lb > kMinIntegerValue) { 424 LinearConstraint canonical_knapsack_form; 425 426 // Negate the variables with positive coefficients. 427 for (int i = 0; i < constraint.vars.size(); ++i) { 428 const IntegerVariable var = constraint.vars[i]; 429 const IntegerValue coeff = constraint.coeffs[i]; 430 if (coeff > IntegerValue(0)) { 431 canonical_knapsack_form.AddTerm(NegationOf(var), coeff); 432 } else { 433 canonical_knapsack_form.AddTerm(var, -coeff); 434 } 435 } 436 canonical_knapsack_form.ub = -constraint.lb; 437 canonical_knapsack_form.lb = kMinIntegerValue; 438 knapsack_constraints->push_back(canonical_knapsack_form); 439 } 440 } 441 442 // TODO(user): This is no longer used as we try to separate all cut with 443 // knapsack now, remove. 444 CutGenerator CreateKnapsackCoverCutGenerator( 445 const std::vector<LinearConstraint>& base_constraints, 446 const std::vector<IntegerVariable>& vars, Model* model) { 447 CutGenerator result; 448 result.vars = vars; 449 450 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>(); 451 std::vector<LinearConstraint> knapsack_constraints; 452 for (const LinearConstraint& constraint : base_constraints) { 453 // There is often a lot of small linear base constraints and it doesn't seem 454 // super useful to generate cuts for constraints of size 2. Any valid cut 455 // of size 1 should be already infered by the propagation. 456 // 457 // TODO(user): The case of size 2 is a bit less clear. investigate more if 458 // it is useful. 459 if (constraint.vars.size() <= 2) continue; 460 461 ConvertToKnapsackForm(constraint, &knapsack_constraints, integer_trail); 462 } 463 VLOG(1) << "#knapsack constraints: " << knapsack_constraints.size(); 464 465 // Note(user): for Knapsack cuts, it seems always advantageous to replace a 466 // variable X by a TIGHT lower bound of the form "coeff * binary + lb". This 467 // will not change "covers" but can only result in more violation by the 468 // current LP solution. 469 ImpliedBoundsProcessor implied_bounds_processor( 470 vars, integer_trail, model->GetOrCreate<ImpliedBounds>()); 471 472 // TODO(user): do not add generator if there are no knapsack constraints. 473 result.generate_cuts = [implied_bounds_processor, knapsack_constraints, vars, 474 model, integer_trail]( 475 const absl::StrongVector<IntegerVariable, double>& 476 lp_values, 477 LinearConstraintManager* manager) mutable { 478 // TODO(user): When we use implied-bound substitution, we might still infer 479 // an interesting cut even if all variables are integer. See if we still 480 // want to skip all such constraints. 481 if (AllVarsTakeIntegerValue(vars, lp_values)) return true; 482 483 KnapsackSolverForCuts knapsack_solver( 484 "Knapsack on demand cover cut generator"); 485 int64_t skipped_constraints = 0; 486 LinearConstraint mutable_constraint; 487 488 // Iterate through all knapsack constraints. 489 implied_bounds_processor.RecomputeCacheAndSeparateSomeImpliedBoundCuts( 490 lp_values); 491 for (const LinearConstraint& constraint : knapsack_constraints) { 492 if (model->GetOrCreate<TimeLimit>()->LimitReached()) break; 493 VLOG(2) << "Processing constraint: " << constraint.DebugString(); 494 495 mutable_constraint = constraint; 496 implied_bounds_processor.ProcessUpperBoundedConstraint( 497 lp_values, &mutable_constraint); 498 MakeAllCoefficientsPositive(&mutable_constraint); 499 500 const LinearConstraint preprocessed_constraint = 501 GetPreprocessedLinearConstraint(mutable_constraint, lp_values, 502 *integer_trail); 503 if (preprocessed_constraint.vars.empty()) continue; 504 505 if (!CanFormValidKnapsackCover(preprocessed_constraint, lp_values, 506 *integer_trail)) { 507 skipped_constraints++; 508 continue; 509 } 510 511 // Profits are (upper_bounds[i] - lp_values[i]) for knapsack variables. 512 std::vector<double> profits; 513 profits.reserve(preprocessed_constraint.vars.size()); 514 515 // Weights are (coeffs[i] * (upper_bound[i] - lower_bound[i])). 516 std::vector<double> weights; 517 weights.reserve(preprocessed_constraint.vars.size()); 518 519 double capacity = -ToDouble(preprocessed_constraint.ub) - 1.0; 520 521 // Compute and store the sum of variable profits. This is the constant 522 // part of the objective of the problem we are trying to solve. Hence 523 // this part is not supplied to the knapsack_solver and is subtracted 524 // when we receive the knapsack solution. 525 double sum_variable_profit = 0; 526 527 // Compute the profits, the weights and the capacity for the knapsack 528 // instance. 529 for (int i = 0; i < preprocessed_constraint.vars.size(); ++i) { 530 const IntegerVariable var = preprocessed_constraint.vars[i]; 531 const double coefficient = ToDouble(preprocessed_constraint.coeffs[i]); 532 const double var_ub = ToDouble(integer_trail->LevelZeroUpperBound(var)); 533 const double var_lb = ToDouble(integer_trail->LevelZeroLowerBound(var)); 534 const double variable_profit = var_ub - lp_values[var]; 535 profits.push_back(variable_profit); 536 537 sum_variable_profit += variable_profit; 538 539 const double weight = coefficient * (var_ub - var_lb); 540 weights.push_back(weight); 541 capacity += weight + coefficient * var_lb; 542 } 543 if (capacity < 0.0) continue; 544 545 std::vector<IntegerVariable> cut_vars; 546 std::vector<IntegerValue> cut_vars_original_coefficients; 547 548 VLOG(2) << "Knapsack size: " << profits.size(); 549 knapsack_solver.Init(profits, weights, capacity); 550 551 // Set the time limit for the knapsack solver. 552 const double time_limit_for_knapsack_solver = 553 model->GetOrCreate<TimeLimit>()->GetTimeLeft(); 554 555 // Solve the instance and subtract the constant part to compute the 556 // sum_of_distance_to_ub_for_vars_in_cover. 557 // TODO(user): Consider solving the instance approximately. 558 bool is_solution_optimal = false; 559 knapsack_solver.set_solution_upper_bound_threshold( 560 sum_variable_profit - 1.0 + kMinCutViolation); 561 // TODO(user): Consider providing lower bound threshold as 562 // sum_variable_profit - 1.0 + kMinCutViolation. 563 // TODO(user): Set node limit for knapsack solver. 564 auto time_limit_for_solver = 565 absl::make_unique<TimeLimit>(time_limit_for_knapsack_solver); 566 const double sum_of_distance_to_ub_for_vars_in_cover = 567 sum_variable_profit - 568 knapsack_solver.Solve(time_limit_for_solver.get(), 569 &is_solution_optimal); 570 if (is_solution_optimal) { 571 VLOG(2) << "Knapsack Optimal solution found yay !"; 572 } 573 if (time_limit_for_solver->LimitReached()) { 574 VLOG(1) << "Knapsack Solver run out of time limit."; 575 } 576 if (sum_of_distance_to_ub_for_vars_in_cover < 1.0 - kMinCutViolation) { 577 // Constraint is eligible for the cover. 578 579 IntegerValue constraint_ub_for_cut = preprocessed_constraint.ub; 580 std::set<IntegerVariable> vars_in_cut; 581 for (int i = 0; i < preprocessed_constraint.vars.size(); ++i) { 582 const IntegerVariable var = preprocessed_constraint.vars[i]; 583 const IntegerValue coefficient = preprocessed_constraint.coeffs[i]; 584 if (!knapsack_solver.best_solution(i)) { 585 cut_vars.push_back(var); 586 cut_vars_original_coefficients.push_back(coefficient); 587 vars_in_cut.insert(var); 588 } else { 589 const IntegerValue var_lb = integer_trail->LevelZeroLowerBound(var); 590 constraint_ub_for_cut -= coefficient * var_lb; 591 } 592 } 593 LinearConstraint cut = GenerateKnapsackCutForCover( 594 cut_vars, cut_vars_original_coefficients, constraint_ub_for_cut, 595 *integer_trail); 596 597 // Check if the constraint has only binary variables. 598 bool is_lifted = false; 599 if (ConstraintIsEligibleForLifting(cut, *integer_trail)) { 600 if (LiftKnapsackCut(mutable_constraint, lp_values, 601 cut_vars_original_coefficients, *integer_trail, 602 model->GetOrCreate<TimeLimit>(), &cut)) { 603 is_lifted = true; 604 } 605 } 606 607 CHECK(!SolutionSatisfiesConstraint(cut, lp_values)); 608 manager->AddCut(cut, is_lifted ? "LiftedKnapsack" : "Knapsack", 609 lp_values); 610 } 611 } 612 if (skipped_constraints > 0) { 613 VLOG(2) << "Skipped constraints: " << skipped_constraints; 614 } 615 return true; 616 }; 617 618 return result; 619 } 620 621 // Compute the larger t <= max_t such that t * rhs_remainder >= divisor / 2. 622 // 623 // This is just a separate function as it is slightly faster to compute the 624 // result only once. 625 IntegerValue GetFactorT(IntegerValue rhs_remainder, IntegerValue divisor, 626 IntegerValue max_t) { 627 DCHECK_GE(max_t, 1); 628 return rhs_remainder == 0 629 ? max_t 630 : std::min(max_t, CeilRatio(divisor / 2, rhs_remainder)); 631 } 632 633 std::function<IntegerValue(IntegerValue)> GetSuperAdditiveRoundingFunction( 634 IntegerValue rhs_remainder, IntegerValue divisor, IntegerValue t, 635 IntegerValue max_scaling) { 636 DCHECK_GE(max_scaling, 1); 637 DCHECK_GE(t, 1); 638 639 // Adjust after the multiplication by t. 640 rhs_remainder *= t; 641 DCHECK_LT(rhs_remainder, divisor); 642 643 // Make sure we don't have an integer overflow below. Note that we assume that 644 // divisor and the maximum coeff magnitude are not too different (maybe a 645 // factor 1000 at most) so that the final result will never overflow. 646 max_scaling = 647 std::min(max_scaling, std::numeric_limits<int64_t>::max() / divisor); 648 649 const IntegerValue size = divisor - rhs_remainder; 650 if (max_scaling == 1 || size == 1) { 651 // TODO(user): Use everywhere a two step computation to avoid overflow? 652 // First divide by divisor, then multiply by t. For now, we limit t so that 653 // we never have an overflow instead. 654 return [t, divisor](IntegerValue coeff) { 655 return FloorRatio(t * coeff, divisor); 656 }; 657 } else if (size <= max_scaling) { 658 return [size, rhs_remainder, t, divisor](IntegerValue coeff) { 659 const IntegerValue t_coeff = t * coeff; 660 const IntegerValue ratio = FloorRatio(t_coeff, divisor); 661 const IntegerValue remainder = PositiveRemainder(t_coeff, divisor); 662 const IntegerValue diff = remainder - rhs_remainder; 663 return size * ratio + std::max(IntegerValue(0), diff); 664 }; 665 } else if (max_scaling.value() * rhs_remainder.value() < divisor) { 666 // Because of our max_t limitation, the rhs_remainder might stay small. 667 // 668 // If it is "too small" we cannot use the code below because it will not be 669 // valid. So we just divide divisor into max_scaling bucket. The 670 // rhs_remainder will be in the bucket 0. 671 // 672 // Note(user): This seems the same as just increasing t, modulo integer 673 // overflows. Maybe we should just always do the computation like this so 674 // that we can use larger t even if coeff is close to kint64max. 675 return [t, divisor, max_scaling](IntegerValue coeff) { 676 const IntegerValue t_coeff = t * coeff; 677 const IntegerValue ratio = FloorRatio(t_coeff, divisor); 678 const IntegerValue remainder = PositiveRemainder(t_coeff, divisor); 679 const IntegerValue bucket = FloorRatio(remainder * max_scaling, divisor); 680 return max_scaling * ratio + bucket; 681 }; 682 } else { 683 // We divide (size = divisor - rhs_remainder) into (max_scaling - 1) buckets 684 // and increase the function by 1 / max_scaling for each of them. 685 // 686 // Note that for different values of max_scaling, we get a family of 687 // functions that do not dominate each others. So potentially, a max scaling 688 // as low as 2 could lead to the better cut (this is exactly the Letchford & 689 // Lodi function). 690 // 691 // Another interesting fact, is that if we want to compute the maximum alpha 692 // for a constraint with 2 terms like: 693 // divisor * Y + (ratio * divisor + remainder) * X 694 // <= rhs_ratio * divisor + rhs_remainder 695 // so that we have the cut: 696 // Y + (ratio + alpha) * X <= rhs_ratio 697 // This is the same as computing the maximum alpha such that for all integer 698 // X > 0 we have CeilRatio(alpha * divisor * X, divisor) 699 // <= CeilRatio(remainder * X - rhs_remainder, divisor). 700 // We can prove that this alpha is of the form (n - 1) / n, and it will 701 // be reached by such function for a max_scaling of n. 702 // 703 // TODO(user): This function is not always maximal when 704 // size % (max_scaling - 1) == 0. Improve? 705 return [size, rhs_remainder, t, divisor, max_scaling](IntegerValue coeff) { 706 const IntegerValue t_coeff = t * coeff; 707 const IntegerValue ratio = FloorRatio(t_coeff, divisor); 708 const IntegerValue remainder = PositiveRemainder(t_coeff, divisor); 709 const IntegerValue diff = remainder - rhs_remainder; 710 const IntegerValue bucket = 711 diff > 0 ? CeilRatio(diff * (max_scaling - 1), size) 712 : IntegerValue(0); 713 return max_scaling * ratio + bucket; 714 }; 715 } 716 } 717 718 // TODO(user): This has been optimized a bit, but we can probably do even better 719 // as it still takes around 25% percent of the run time when all the cuts are on 720 // for the opm*mps.gz problems and others. 721 void IntegerRoundingCutHelper::ComputeCut( 722 RoundingOptions options, const std::vector<double>& lp_values, 723 const std::vector<IntegerValue>& lower_bounds, 724 const std::vector<IntegerValue>& upper_bounds, 725 ImpliedBoundsProcessor* ib_processor, LinearConstraint* cut) { 726 const int size = lp_values.size(); 727 if (size == 0) return; 728 CHECK_EQ(lower_bounds.size(), size); 729 CHECK_EQ(upper_bounds.size(), size); 730 CHECK_EQ(cut->vars.size(), size); 731 CHECK_EQ(cut->coeffs.size(), size); 732 CHECK_EQ(cut->lb, kMinIntegerValue); 733 734 // To optimize the computation of the best divisor below, we only need to 735 // look at the indices with a shifted lp value that is not close to zero. 736 // 737 // TODO(user): sort by decreasing lp_values so that our early abort test in 738 // the critical loop below has more chance of returning early? I tried but it 739 // didn't seems to change much though. 740 relevant_indices_.clear(); 741 relevant_lp_values_.clear(); 742 relevant_coeffs_.clear(); 743 relevant_bound_diffs_.clear(); 744 divisors_.clear(); 745 adjusted_coeffs_.clear(); 746 747 // Compute the maximum magnitude for non-fixed variables. 748 IntegerValue max_magnitude(0); 749 for (int i = 0; i < size; ++i) { 750 if (lower_bounds[i] == upper_bounds[i]) continue; 751 const IntegerValue magnitude = IntTypeAbs(cut->coeffs[i]); 752 max_magnitude = std::max(max_magnitude, magnitude); 753 } 754 755 // Shift each variable using its lower/upper bound so that no variable can 756 // change sign. We eventually do a change of variable to its negation so 757 // that all variable are non-negative. 758 bool overflow = false; 759 change_sign_at_postprocessing_.assign(size, false); 760 for (int i = 0; i < size; ++i) { 761 if (cut->coeffs[i] == 0) continue; 762 const IntegerValue magnitude = IntTypeAbs(cut->coeffs[i]); 763 764 // We might change them below. 765 IntegerValue lb = lower_bounds[i]; 766 double lp_value = lp_values[i]; 767 768 const IntegerValue ub = upper_bounds[i]; 769 const IntegerValue bound_diff = 770 IntegerValue(CapSub(ub.value(), lb.value())); 771 // Note that since we use ToDouble() this code works fine with lb/ub at 772 // min/max integer value. 773 // 774 // TODO(user): Experiments with different heuristics. Other solver also 775 // seems to try a bunch of possibilities in a "postprocess" phase once 776 // the divisor is chosen. Try that. 777 { 778 // when the magnitude of the entry become smaller and smaller we bias 779 // towards a positive coefficient. This is because after rounding this 780 // will likely become zero instead of -divisor and we need the lp value 781 // to be really close to its bound to compensate. 782 const double lb_dist = std::abs(lp_value - ToDouble(lb)); 783 const double ub_dist = std::abs(lp_value - ToDouble(ub)); 784 const double bias = 785 std::max(1.0, 0.1 * ToDouble(max_magnitude) / ToDouble(magnitude)); 786 if ((bias * lb_dist > ub_dist && cut->coeffs[i] < 0) || 787 (lb_dist > bias * ub_dist && cut->coeffs[i] > 0)) { 788 change_sign_at_postprocessing_[i] = true; 789 cut->coeffs[i] = -cut->coeffs[i]; 790 lp_value = -lp_value; 791 lb = -ub; 792 } 793 } 794 795 // Always shift to lb. 796 // coeff * X = coeff * (X - shift) + coeff * shift. 797 lp_value -= ToDouble(lb); 798 if (!AddProductTo(-cut->coeffs[i], lb, &cut->ub)) { 799 overflow = true; 800 break; 801 } 802 803 // Deal with fixed variable, no need to shift back in this case, we can 804 // just remove the term. 805 if (bound_diff == 0) { 806 cut->coeffs[i] = IntegerValue(0); 807 lp_value = 0.0; 808 } 809 810 if (std::abs(lp_value) > 1e-2) { 811 relevant_coeffs_.push_back(cut->coeffs[i]); 812 relevant_indices_.push_back(i); 813 relevant_lp_values_.push_back(lp_value); 814 relevant_bound_diffs_.push_back(bound_diff); 815 divisors_.push_back(magnitude); 816 } 817 } 818 819 // TODO(user): Maybe this shouldn't be called on such constraint. 820 if (relevant_coeffs_.empty()) { 821 VLOG(2) << "Issue, nothing to cut."; 822 *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); 823 return; 824 } 825 CHECK_NE(max_magnitude, 0); 826 827 // Our heuristic will try to generate a few different cuts, and we will keep 828 // the most violated one scaled by the l2 norm of the relevant position. 829 // 830 // TODO(user): Experiment for the best value of this initial violation 831 // threshold. Note also that we use the l2 norm on the restricted position 832 // here. Maybe we should change that? On that note, the L2 norm usage seems a 833 // bit weird to me since it grows with the number of term in the cut. And 834 // often, we already have a good cut, and we make it stronger by adding extra 835 // terms that do not change its activity. 836 // 837 // The discussion above only concern the best_scaled_violation initial value. 838 // The remainder_threshold allows to not consider cuts for which the final 839 // efficacity is clearly lower than 1e-3 (it is a bound, so we could generate 840 // cuts with a lower efficacity than this). 841 double best_scaled_violation = 0.01; 842 const IntegerValue remainder_threshold(max_magnitude / 1000); 843 844 // The cut->ub might have grown quite a bit with the bound substitution, so 845 // we need to include it too since we will apply the rounding function on it. 846 max_magnitude = std::max(max_magnitude, IntTypeAbs(cut->ub)); 847 848 // Make sure that when we multiply the rhs or the coefficient by a factor t, 849 // we do not have an integer overflow. Actually, we need a bit more room 850 // because we might round down a value to the next multiple of 851 // max_magnitude. 852 const IntegerValue threshold = kMaxIntegerValue / 2; 853 if (overflow || max_magnitude >= threshold) { 854 VLOG(2) << "Issue, overflow."; 855 *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); 856 return; 857 } 858 const IntegerValue max_t = threshold / max_magnitude; 859 860 // There is no point trying twice the same divisor or a divisor that is too 861 // small. Note that we use a higher threshold than the remainder_threshold 862 // because we can boost the remainder thanks to our adjusting heuristic below 863 // and also because this allows to have cuts with a small range of 864 // coefficients. 865 // 866 // TODO(user): Note that the std::sort() is visible in some cpu profile. 867 { 868 int new_size = 0; 869 const IntegerValue divisor_threshold = max_magnitude / 10; 870 for (int i = 0; i < divisors_.size(); ++i) { 871 if (divisors_[i] <= divisor_threshold) continue; 872 divisors_[new_size++] = divisors_[i]; 873 } 874 divisors_.resize(new_size); 875 } 876 gtl::STLSortAndRemoveDuplicates(&divisors_, std::greater<IntegerValue>()); 877 878 // TODO(user): Avoid quadratic algorithm? Note that we are quadratic in 879 // relevant_indices not the full cut->coeffs.size(), but this is still too 880 // much on some problems. 881 IntegerValue best_divisor(0); 882 for (const IntegerValue divisor : divisors_) { 883 // Skip if we don't have the potential to generate a good enough cut. 884 const IntegerValue initial_rhs_remainder = 885 cut->ub - FloorRatio(cut->ub, divisor) * divisor; 886 if (initial_rhs_remainder <= remainder_threshold) continue; 887 888 IntegerValue temp_ub = cut->ub; 889 adjusted_coeffs_.clear(); 890 891 // We will adjust coefficient that are just under an exact multiple of 892 // divisor to an exact multiple. This is meant to get rid of small errors 893 // that appears due to rounding error in our exact computation of the 894 // initial constraint given to this class. 895 // 896 // Each adjustement will cause the initial_rhs_remainder to increase, and we 897 // do not want to increase it above divisor. Our threshold below guarantees 898 // this. Note that the higher the rhs_remainder becomes, the more the 899 // function f() has a chance to reduce the violation, so it is not always a 900 // good idea to use all the slack we have between initial_rhs_remainder and 901 // divisor. 902 // 903 // TODO(user): If possible, it might be better to complement these 904 // variables. Even if the adjusted lp_values end up larger, if we loose less 905 // when taking f(), then we will have a better violation. 906 const IntegerValue adjust_threshold = 907 (divisor - initial_rhs_remainder - 1) / IntegerValue(size); 908 if (adjust_threshold > 0) { 909 // Even before we finish the adjust, we can have a lower bound on the 910 // activily loss using this divisor, and so we can abort early. This is 911 // similar to what is done below in the function. 912 bool early_abort = false; 913 double loss_lb = 0.0; 914 const double threshold = ToDouble(initial_rhs_remainder); 915 916 for (int i = 0; i < relevant_coeffs_.size(); ++i) { 917 // Compute the difference of coeff with the next multiple of divisor. 918 const IntegerValue coeff = relevant_coeffs_[i]; 919 const IntegerValue remainder = 920 CeilRatio(coeff, divisor) * divisor - coeff; 921 922 if (divisor - remainder <= initial_rhs_remainder) { 923 // We do not know exactly f() yet, but it will always round to the 924 // floor of the division by divisor in this case. 925 loss_lb += ToDouble(divisor - remainder) * relevant_lp_values_[i]; 926 if (loss_lb >= threshold) { 927 early_abort = true; 928 break; 929 } 930 } 931 932 // Adjust coeff of the form k * divisor - epsilon. 933 const IntegerValue diff = relevant_bound_diffs_[i]; 934 if (remainder > 0 && remainder <= adjust_threshold && 935 CapProd(diff.value(), remainder.value()) <= adjust_threshold) { 936 temp_ub += remainder * diff; 937 adjusted_coeffs_.push_back({i, coeff + remainder}); 938 } 939 } 940 941 if (early_abort) continue; 942 } 943 944 // Create the super-additive function f(). 945 const IntegerValue rhs_remainder = 946 temp_ub - FloorRatio(temp_ub, divisor) * divisor; 947 if (rhs_remainder == 0) continue; 948 949 const auto f = GetSuperAdditiveRoundingFunction( 950 rhs_remainder, divisor, GetFactorT(rhs_remainder, divisor, max_t), 951 options.max_scaling); 952 953 // As we round coefficients, we will compute the loss compared to the 954 // current scaled constraint activity. As soon as this loss crosses the 955 // slack, then we known that there is no violation and we can abort early. 956 // 957 // TODO(user): modulo the scaling, we could compute the exact threshold 958 // using our current best cut. Note that we also have to account the change 959 // in slack due to the adjust code above. 960 const double scaling = ToDouble(f(divisor)) / ToDouble(divisor); 961 const double threshold = scaling * ToDouble(rhs_remainder); 962 double loss = 0.0; 963 964 // Apply f() to the cut and compute the cut violation. Note that it is 965 // okay to just look at the relevant indices since the other have a lp 966 // value which is almost zero. Doing it like this is faster, and even if 967 // the max_magnitude might be off it should still be relevant enough. 968 double violation = -ToDouble(f(temp_ub)); 969 double l2_norm = 0.0; 970 bool early_abort = false; 971 int adjusted_coeffs_index = 0; 972 for (int i = 0; i < relevant_coeffs_.size(); ++i) { 973 IntegerValue coeff = relevant_coeffs_[i]; 974 975 // Adjust coeff according to our previous computation if needed. 976 if (adjusted_coeffs_index < adjusted_coeffs_.size() && 977 adjusted_coeffs_[adjusted_coeffs_index].first == i) { 978 coeff = adjusted_coeffs_[adjusted_coeffs_index].second; 979 adjusted_coeffs_index++; 980 } 981 982 if (coeff == 0) continue; 983 const IntegerValue new_coeff = f(coeff); 984 const double new_coeff_double = ToDouble(new_coeff); 985 const double lp_value = relevant_lp_values_[i]; 986 987 l2_norm += new_coeff_double * new_coeff_double; 988 violation += new_coeff_double * lp_value; 989 loss += (scaling * ToDouble(coeff) - new_coeff_double) * lp_value; 990 if (loss >= threshold) { 991 early_abort = true; 992 break; 993 } 994 } 995 if (early_abort) continue; 996 997 // Here we scale by the L2 norm over the "relevant" positions. This seems 998 // to work slighly better in practice. 999 violation /= sqrt(l2_norm); 1000 if (violation > best_scaled_violation) { 1001 best_scaled_violation = violation; 1002 best_divisor = divisor; 1003 } 1004 } 1005 1006 if (best_divisor == 0) { 1007 *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); 1008 return; 1009 } 1010 1011 // Adjust coefficients. 1012 // 1013 // TODO(user): It might make sense to also adjust the one with a small LP 1014 // value, but then the cut will be slighlty different than the one we computed 1015 // above. Try with and without maybe? 1016 const IntegerValue initial_rhs_remainder = 1017 cut->ub - FloorRatio(cut->ub, best_divisor) * best_divisor; 1018 const IntegerValue adjust_threshold = 1019 (best_divisor - initial_rhs_remainder - 1) / IntegerValue(size); 1020 if (adjust_threshold > 0) { 1021 for (int i = 0; i < relevant_indices_.size(); ++i) { 1022 const int index = relevant_indices_[i]; 1023 const IntegerValue diff = relevant_bound_diffs_[i]; 1024 if (diff > adjust_threshold) continue; 1025 1026 // Adjust coeff of the form k * best_divisor - epsilon. 1027 const IntegerValue coeff = cut->coeffs[index]; 1028 const IntegerValue remainder = 1029 CeilRatio(coeff, best_divisor) * best_divisor - coeff; 1030 if (CapProd(diff.value(), remainder.value()) <= adjust_threshold) { 1031 cut->ub += remainder * diff; 1032 cut->coeffs[index] += remainder; 1033 } 1034 } 1035 } 1036 1037 // Create the super-additive function f(). 1038 // 1039 // TODO(user): Try out different rounding function and keep the best. We can 1040 // change max_t and max_scaling. It might not be easy to choose which cut is 1041 // the best, but we can at least know for sure if one dominate the other 1042 // completely. That is, if for all coeff f(coeff)/f(divisor) is greater than 1043 // or equal to the same value for another function f. 1044 const IntegerValue rhs_remainder = 1045 cut->ub - FloorRatio(cut->ub, best_divisor) * best_divisor; 1046 IntegerValue factor_t = GetFactorT(rhs_remainder, best_divisor, max_t); 1047 auto f = GetSuperAdditiveRoundingFunction(rhs_remainder, best_divisor, 1048 factor_t, options.max_scaling); 1049 1050 // Look amongst all our possible function f() for one that dominate greedily 1051 // our current best one. Note that we prefer lower scaling factor since that 1052 // result in a cut with lower coefficients. 1053 remainders_.clear(); 1054 for (int i = 0; i < size; ++i) { 1055 const IntegerValue coeff = cut->coeffs[i]; 1056 const IntegerValue r = 1057 coeff - FloorRatio(coeff, best_divisor) * best_divisor; 1058 if (r > rhs_remainder) remainders_.push_back(r); 1059 } 1060 gtl::STLSortAndRemoveDuplicates(&remainders_); 1061 if (remainders_.size() <= 100) { 1062 best_rs_.clear(); 1063 for (const IntegerValue r : remainders_) { 1064 best_rs_.push_back(f(r)); 1065 } 1066 IntegerValue best_d = f(best_divisor); 1067 1068 // Note that the complexity seems high 100 * 2 * options.max_scaling, but 1069 // this only run on cuts that are already efficient and the inner loop tend 1070 // to abort quickly. I didn't see this code in the cpu profile so far. 1071 for (const IntegerValue t : 1072 {IntegerValue(1), GetFactorT(rhs_remainder, best_divisor, max_t)}) { 1073 for (IntegerValue s(2); s <= options.max_scaling; ++s) { 1074 const auto g = 1075 GetSuperAdditiveRoundingFunction(rhs_remainder, best_divisor, t, s); 1076 int num_strictly_better = 0; 1077 rs_.clear(); 1078 const IntegerValue d = g(best_divisor); 1079 for (int i = 0; i < best_rs_.size(); ++i) { 1080 const IntegerValue temp = g(remainders_[i]); 1081 if (temp * best_d < best_rs_[i] * d) break; 1082 if (temp * best_d > best_rs_[i] * d) num_strictly_better++; 1083 rs_.push_back(temp); 1084 } 1085 if (rs_.size() == best_rs_.size() && num_strictly_better > 0) { 1086 f = g; 1087 factor_t = t; 1088 best_rs_ = rs_; 1089 best_d = d; 1090 } 1091 } 1092 } 1093 } 1094 1095 // Starts to apply f() to the cut. We only apply it to the rhs here, the 1096 // coefficient will be done after the potential lifting of some Booleans. 1097 cut->ub = f(cut->ub); 1098 tmp_terms_.clear(); 1099 1100 // Lift some implied bounds Booleans. Note that we will add them after 1101 // "size" so they will be ignored in the second loop below. 1102 num_lifted_booleans_ = 0; 1103 if (ib_processor != nullptr) { 1104 for (int i = 0; i < size; ++i) { 1105 const IntegerValue coeff = cut->coeffs[i]; 1106 if (coeff == 0) continue; 1107 1108 IntegerVariable var = cut->vars[i]; 1109 if (change_sign_at_postprocessing_[i]) { 1110 var = NegationOf(var); 1111 } 1112 1113 const ImpliedBoundsProcessor::BestImpliedBoundInfo info = 1114 ib_processor->GetCachedImpliedBoundInfo(var); 1115 1116 // Avoid overflow. 1117 if (CapProd(CapProd(std::abs(coeff.value()), factor_t.value()), 1118 info.bound_diff.value()) == 1119 std::numeric_limits<int64_t>::max()) { 1120 continue; 1121 } 1122 1123 // Because X = bound_diff * B + S 1124 // We can replace coeff * X by the expression before applying f: 1125 // = f(coeff * bound_diff) * B + f(coeff) * [X - bound_diff * B] 1126 // = f(coeff) * X + (f(coeff * bound_diff) - f(coeff) * bound_diff] B 1127 // So we can "lift" B into the cut. 1128 const IntegerValue coeff_b = 1129 f(coeff * info.bound_diff) - f(coeff) * info.bound_diff; 1130 CHECK_GE(coeff_b, 0); 1131 if (coeff_b == 0) continue; 1132 1133 ++num_lifted_booleans_; 1134 if (info.is_positive) { 1135 tmp_terms_.push_back({info.bool_var, coeff_b}); 1136 } else { 1137 tmp_terms_.push_back({info.bool_var, -coeff_b}); 1138 cut->ub = CapAdd(-coeff_b.value(), cut->ub.value()); 1139 } 1140 } 1141 } 1142 1143 // Apply f() to the cut. 1144 // 1145 // Remove the bound shifts so the constraint is expressed in the original 1146 // variables. 1147 for (int i = 0; i < size; ++i) { 1148 IntegerValue coeff = cut->coeffs[i]; 1149 if (coeff == 0) continue; 1150 coeff = f(coeff); 1151 if (coeff == 0) continue; 1152 if (change_sign_at_postprocessing_[i]) { 1153 if (!AddProductTo(coeff, -upper_bounds[i], &cut->ub)) { 1154 // Abort with a trivially satisfied cut. 1155 cut->Clear(); 1156 return; 1157 } 1158 tmp_terms_.push_back({cut->vars[i], -coeff}); 1159 } else { 1160 if (!AddProductTo(coeff, lower_bounds[i], &cut->ub)) { 1161 // Abort with a trivially satisfied cut. 1162 cut->Clear(); 1163 return; 1164 } 1165 tmp_terms_.push_back({cut->vars[i], coeff}); 1166 } 1167 } 1168 1169 // Basic post-processing. 1170 CleanTermsAndFillConstraint(&tmp_terms_, cut); 1171 RemoveZeroTerms(cut); 1172 DivideByGCD(cut); 1173 } 1174 1175 bool CoverCutHelper::TrySimpleKnapsack( 1176 const LinearConstraint base_ct, const std::vector<double>& lp_values, 1177 const std::vector<IntegerValue>& lower_bounds, 1178 const std::vector<IntegerValue>& upper_bounds) { 1179 const int base_size = lp_values.size(); 1180 1181 // Fill terms with a rewrite of the base constraint where all coeffs & 1182 // variables are positive by using either (X - LB) or (UB - X) as new 1183 // variables. 1184 terms_.clear(); 1185 IntegerValue rhs = base_ct.ub; 1186 IntegerValue sum_of_diff(0); 1187 IntegerValue max_base_magnitude(0); 1188 for (int i = 0; i < base_size; ++i) { 1189 const IntegerValue coeff = base_ct.coeffs[i]; 1190 const IntegerValue positive_coeff = IntTypeAbs(coeff); 1191 max_base_magnitude = std::max(max_base_magnitude, positive_coeff); 1192 const IntegerValue bound_diff = upper_bounds[i] - lower_bounds[i]; 1193 if (!AddProductTo(positive_coeff, bound_diff, &sum_of_diff)) { 1194 return false; 1195 } 1196 const IntegerValue diff = positive_coeff * bound_diff; 1197 if (coeff > 0) { 1198 if (!AddProductTo(-coeff, lower_bounds[i], &rhs)) return false; 1199 terms_.push_back( 1200 {i, ToDouble(upper_bounds[i]) - lp_values[i], positive_coeff, diff}); 1201 } else { 1202 if (!AddProductTo(-coeff, upper_bounds[i], &rhs)) return false; 1203 terms_.push_back( 1204 {i, lp_values[i] - ToDouble(lower_bounds[i]), positive_coeff, diff}); 1205 } 1206 } 1207 1208 // Try a simple cover heuristic. 1209 // Look for violated CUT of the form: sum (UB - X) or (X - LB) >= 1. 1210 double activity = 0.0; 1211 int new_size = 0; 1212 std::sort(terms_.begin(), terms_.end(), [](const Term& a, const Term& b) { 1213 if (a.dist_to_max_value == b.dist_to_max_value) { 1214 // Prefer low coefficients if the distance is the same. 1215 return a.positive_coeff < b.positive_coeff; 1216 } 1217 return a.dist_to_max_value < b.dist_to_max_value; 1218 }); 1219 for (int i = 0; i < terms_.size(); ++i) { 1220 const Term& term = terms_[i]; 1221 activity += term.dist_to_max_value; 1222 1223 // As an heuristic we select all the term so that the sum of distance 1224 // to the upper bound is <= 1.0. If the corresponding rhs is negative, then 1225 // we will have a cut of violation at least 0.0. Note that this violation 1226 // can be improved by the lifting. 1227 // 1228 // TODO(user): experiment with different threshold (even greater than one). 1229 // Or come up with an algo that incorporate the lifting into the heuristic. 1230 if (activity > 1.0) { 1231 new_size = i; // before this entry. 1232 break; 1233 } 1234 1235 rhs -= term.diff; 1236 } 1237 1238 // If the rhs is now negative, we have a cut. 1239 // 1240 // Note(user): past this point, now that a given "base" cover has been chosen, 1241 // we basically compute the cut (of the form sum X <= bound) with the maximum 1242 // possible violation. Note also that we lift as much as possible, so we don't 1243 // necessarily optimize for the cut efficacity though. But we do get a 1244 // stronger cut. 1245 if (rhs >= 0) return false; 1246 if (new_size == 0) return false; 1247 1248 // Transform to a minimal cover. We want to greedily remove the largest coeff 1249 // first, so we have more chance for the "lifting" below which can increase 1250 // the cut violation. If the coeff are the same, we prefer to remove high 1251 // distance from upper bound first. 1252 // 1253 // We compute the cut at the same time. 1254 terms_.resize(new_size); 1255 std::sort(terms_.begin(), terms_.end(), [](const Term& a, const Term& b) { 1256 if (a.positive_coeff == b.positive_coeff) { 1257 return a.dist_to_max_value > b.dist_to_max_value; 1258 } 1259 return a.positive_coeff > b.positive_coeff; 1260 }); 1261 in_cut_.assign(base_ct.vars.size(), false); 1262 cut_.ClearTerms(); 1263 cut_.lb = kMinIntegerValue; 1264 cut_.ub = IntegerValue(-1); 1265 IntegerValue max_coeff(0); 1266 for (const Term term : terms_) { 1267 if (term.diff + rhs < 0) { 1268 rhs += term.diff; 1269 continue; 1270 } 1271 in_cut_[term.index] = true; 1272 max_coeff = std::max(max_coeff, term.positive_coeff); 1273 cut_.vars.push_back(base_ct.vars[term.index]); 1274 if (base_ct.coeffs[term.index] > 0) { 1275 cut_.coeffs.push_back(IntegerValue(1)); 1276 cut_.ub += upper_bounds[term.index]; 1277 } else { 1278 cut_.coeffs.push_back(IntegerValue(-1)); 1279 cut_.ub -= lower_bounds[term.index]; 1280 } 1281 } 1282 1283 // In case the max_coeff variable is not binary, it might be possible to 1284 // tighten the cut a bit more. 1285 // 1286 // Note(user): I never observed this on the miplib so far. 1287 if (max_coeff == 0) return true; 1288 if (max_coeff < -rhs) { 1289 const IntegerValue m = FloorRatio(-rhs - 1, max_coeff); 1290 rhs += max_coeff * m; 1291 cut_.ub -= m; 1292 } 1293 CHECK_LT(rhs, 0); 1294 1295 // Lift all at once the variables not used in the cover. 1296 // 1297 // We have a cut of the form sum_i X_i <= b that we will lift into 1298 // sum_i scaling X_i + sum f(base_coeff_j) X_j <= b * scaling. 1299 // 1300 // Using the super additivity of f() and how we construct it, 1301 // we know that: sum_j base_coeff_j X_j <= N * max_coeff + (max_coeff - slack) 1302 // implies that: sum_j f(base_coeff_j) X_j <= N * scaling. 1303 // 1304 // 1/ cut > b -(N+1) => original sum + (N+1) * max_coeff >= rhs + slack 1305 // 2/ rewrite 1/ as : scaling * cut >= scaling * b - scaling * N => ... 1306 // 3/ lift > N * scaling => lift_sum > N * max_coeff + (max_coeff - slack) 1307 // And adding 2/ + 3/ we prove what we want: 1308 // cut * scaling + lift > b * scaling => original_sum + lift_sum > rhs. 1309 const IntegerValue slack = -rhs; 1310 const IntegerValue remainder = max_coeff - slack; 1311 max_base_magnitude = std::max(max_base_magnitude, IntTypeAbs(cut_.ub)); 1312 const IntegerValue max_scaling(std::min( 1313 IntegerValue(60), FloorRatio(kMaxIntegerValue, max_base_magnitude))); 1314 const auto f = GetSuperAdditiveRoundingFunction(remainder, max_coeff, 1315 IntegerValue(1), max_scaling); 1316 1317 const IntegerValue scaling = f(max_coeff); 1318 if (scaling > 1) { 1319 for (int i = 0; i < cut_.coeffs.size(); ++i) cut_.coeffs[i] *= scaling; 1320 cut_.ub *= scaling; 1321 } 1322 1323 num_lifting_ = 0; 1324 for (int i = 0; i < base_size; ++i) { 1325 if (in_cut_[i]) continue; 1326 const IntegerValue positive_coeff = IntTypeAbs(base_ct.coeffs[i]); 1327 const IntegerValue new_coeff = f(positive_coeff); 1328 if (new_coeff == 0) continue; 1329 1330 ++num_lifting_; 1331 if (base_ct.coeffs[i] > 0) { 1332 // Add new_coeff * (X - LB) 1333 cut_.coeffs.push_back(new_coeff); 1334 cut_.vars.push_back(base_ct.vars[i]); 1335 cut_.ub += lower_bounds[i] * new_coeff; 1336 } else { 1337 // Add new_coeff * (UB - X) 1338 cut_.coeffs.push_back(-new_coeff); 1339 cut_.vars.push_back(base_ct.vars[i]); 1340 cut_.ub -= upper_bounds[i] * new_coeff; 1341 } 1342 } 1343 1344 if (scaling > 1) DivideByGCD(&cut_); 1345 return true; 1346 } 1347 1348 CutGenerator CreatePositiveMultiplicationCutGenerator(AffineExpression z, 1349 AffineExpression x, 1350 AffineExpression y, 1351 int linearization_level, 1352 Model* model) { 1353 CutGenerator result; 1354 if (z.var != kNoIntegerVariable) result.vars.push_back(z.var); 1355 if (x.var != kNoIntegerVariable) result.vars.push_back(x.var); 1356 if (y.var != kNoIntegerVariable) result.vars.push_back(y.var); 1357 1358 IntegerTrail* const integer_trail = model->GetOrCreate<IntegerTrail>(); 1359 Trail* trail = model->GetOrCreate<Trail>(); 1360 1361 result.generate_cuts = 1362 [z, x, y, linearization_level, model, trail, integer_trail]( 1363 const absl::StrongVector<IntegerVariable, double>& lp_values, 1364 LinearConstraintManager* manager) { 1365 if (trail->CurrentDecisionLevel() > 0 && linearization_level == 1) { 1366 return true; 1367 } 1368 const int64_t x_lb = integer_trail->LevelZeroLowerBound(x).value(); 1369 const int64_t x_ub = integer_trail->LevelZeroUpperBound(x).value(); 1370 const int64_t y_lb = integer_trail->LevelZeroLowerBound(y).value(); 1371 const int64_t y_ub = integer_trail->LevelZeroUpperBound(y).value(); 1372 1373 // TODO(user): Compute a better bound (int_max / 4 ?). 1374 const int64_t kMaxSafeInteger = (int64_t{1} << 53) - 1; 1375 1376 if (CapProd(x_ub, y_ub) >= kMaxSafeInteger) { 1377 VLOG(3) << "Potential overflow in PositiveMultiplicationCutGenerator"; 1378 return true; 1379 } 1380 1381 const double x_lp_value = x.LpValue(lp_values); 1382 const double y_lp_value = y.LpValue(lp_values); 1383 const double z_lp_value = z.LpValue(lp_values); 1384 1385 // TODO(user): As the bounds change monotonically, these cuts 1386 // dominate any previous one. try to keep a reference to the cut and 1387 // replace it. Alternatively, add an API for a level-zero bound change 1388 // callback. 1389 1390 // Cut -z + x_coeff * x + y_coeff* y <= rhs 1391 auto try_add_above_cut = 1392 [manager, z_lp_value, x_lp_value, y_lp_value, x, y, z, model, 1393 &lp_values](int64_t x_coeff, int64_t y_coeff, int64_t rhs) { 1394 if (-z_lp_value + x_lp_value * x_coeff + y_lp_value * y_coeff >= 1395 rhs + kMinCutViolation) { 1396 LinearConstraintBuilder cut(model, /*lb=*/kMinIntegerValue, 1397 /*ub=*/IntegerValue(rhs)); 1398 cut.AddTerm(z, IntegerValue(-1)); 1399 if (x_coeff != 0) cut.AddTerm(x, IntegerValue(x_coeff)); 1400 if (y_coeff != 0) cut.AddTerm(y, IntegerValue(y_coeff)); 1401 manager->AddCut(cut.Build(), "PositiveProduct", lp_values); 1402 } 1403 }; 1404 1405 // Cut -z + x_coeff * x + y_coeff* y >= rhs 1406 auto try_add_below_cut = 1407 [manager, z_lp_value, x_lp_value, y_lp_value, x, y, z, model, 1408 &lp_values](int64_t x_coeff, int64_t y_coeff, int64_t rhs) { 1409 if (-z_lp_value + x_lp_value * x_coeff + y_lp_value * y_coeff <= 1410 rhs - kMinCutViolation) { 1411 LinearConstraintBuilder cut(model, /*lb=*/IntegerValue(rhs), 1412 /*ub=*/kMaxIntegerValue); 1413 cut.AddTerm(z, IntegerValue(-1)); 1414 if (x_coeff != 0) cut.AddTerm(x, IntegerValue(x_coeff)); 1415 if (y_coeff != 0) cut.AddTerm(y, IntegerValue(y_coeff)); 1416 manager->AddCut(cut.Build(), "PositiveProduct", lp_values); 1417 } 1418 }; 1419 1420 // McCormick relaxation of bilinear constraints. These 4 cuts are the 1421 // exact facets of the x * y polyhedron for a bounded x and y. 1422 // 1423 // Each cut correspond to plane that contains two of the line 1424 // (x=x_lb), (x=x_ub), (y=y_lb), (y=y_ub). The easiest to 1425 // understand them is to draw the x*y curves and see the 4 1426 // planes that correspond to the convex hull of the graph. 1427 try_add_above_cut(y_lb, x_lb, x_lb * y_lb); 1428 try_add_above_cut(y_ub, x_ub, x_ub * y_ub); 1429 try_add_below_cut(y_ub, x_lb, x_lb * y_ub); 1430 try_add_below_cut(y_lb, x_ub, x_ub * y_lb); 1431 return true; 1432 }; 1433 1434 return result; 1435 } 1436 1437 CutGenerator CreateSquareCutGenerator(AffineExpression y, AffineExpression x, 1438 int linearization_level, Model* model) { 1439 CutGenerator result; 1440 if (x.var != kNoIntegerVariable) result.vars.push_back(x.var); 1441 if (y.var != kNoIntegerVariable) result.vars.push_back(y.var); 1442 1443 Trail* trail = model->GetOrCreate<Trail>(); 1444 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>(); 1445 result.generate_cuts = 1446 [y, x, linearization_level, trail, integer_trail, model]( 1447 const absl::StrongVector<IntegerVariable, double>& lp_values, 1448 LinearConstraintManager* manager) { 1449 if (trail->CurrentDecisionLevel() > 0 && linearization_level == 1) { 1450 return true; 1451 } 1452 const int64_t x_ub = integer_trail->LevelZeroUpperBound(x).value(); 1453 const int64_t x_lb = integer_trail->LevelZeroLowerBound(x).value(); 1454 1455 if (x_lb == x_ub) return true; 1456 1457 // Check for potential overflows. 1458 if (x_ub > (int64_t{1} << 31)) return true; 1459 DCHECK_GE(x_lb, 0); 1460 1461 const double y_lp_value = y.LpValue(lp_values); 1462 const double x_lp_value = x.LpValue(lp_values); 1463 1464 // First cut: target should be below the line: 1465 // (x_lb, x_lb ^ 2) to (x_ub, x_ub ^ 2). 1466 // The slope of that line is (ub^2 - lb^2) / (ub - lb) = ub + lb. 1467 const int64_t y_lb = x_lb * x_lb; 1468 const int64_t above_slope = x_ub + x_lb; 1469 const double max_lp_y = y_lb + above_slope * (x_lp_value - x_lb); 1470 if (y_lp_value >= max_lp_y + kMinCutViolation) { 1471 // cut: y <= (x_lb + x_ub) * x - x_lb * x_ub 1472 LinearConstraintBuilder above_cut(model, kMinIntegerValue, 1473 IntegerValue(-x_lb * x_ub)); 1474 above_cut.AddTerm(y, IntegerValue(1)); 1475 above_cut.AddTerm(x, IntegerValue(-above_slope)); 1476 manager->AddCut(above_cut.Build(), "SquareUpper", lp_values); 1477 } 1478 1479 // Second cut: target should be above all the lines 1480 // (value, value ^ 2) to (value + 1, (value + 1) ^ 2) 1481 // The slope of that line is 2 * value + 1 1482 // 1483 // Note that we only add one of these cuts. The one for x_lp_value in 1484 // [value, value + 1]. 1485 const int64_t x_floor = static_cast<int64_t>(std::floor(x_lp_value)); 1486 const int64_t below_slope = 2 * x_floor + 1; 1487 const double min_lp_y = 1488 below_slope * x_lp_value - x_floor - x_floor * x_floor; 1489 if (min_lp_y >= y_lp_value + kMinCutViolation) { 1490 // cut: y >= below_slope * (x - x_floor) + x_floor ^ 2 1491 // : y >= below_slope * x - x_floor ^ 2 - x_floor 1492 LinearConstraintBuilder below_cut( 1493 model, IntegerValue(-x_floor - x_floor * x_floor), 1494 kMaxIntegerValue); 1495 below_cut.AddTerm(y, IntegerValue(1)); 1496 below_cut.AddTerm(x, -IntegerValue(below_slope)); 1497 manager->AddCut(below_cut.Build(), "SquareLower", lp_values); 1498 } 1499 return true; 1500 }; 1501 1502 return result; 1503 } 1504 1505 void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint( 1506 const absl::StrongVector<IntegerVariable, double>& lp_values, 1507 LinearConstraint* cut) { 1508 ProcessUpperBoundedConstraintWithSlackCreation( 1509 /*substitute_only_inner_variables=*/false, IntegerVariable(0), lp_values, 1510 cut, nullptr); 1511 } 1512 1513 ImpliedBoundsProcessor::BestImpliedBoundInfo 1514 ImpliedBoundsProcessor::GetCachedImpliedBoundInfo(IntegerVariable var) { 1515 auto it = cache_.find(var); 1516 if (it != cache_.end()) return it->second; 1517 return BestImpliedBoundInfo(); 1518 } 1519 1520 ImpliedBoundsProcessor::BestImpliedBoundInfo 1521 ImpliedBoundsProcessor::ComputeBestImpliedBound( 1522 IntegerVariable var, 1523 const absl::StrongVector<IntegerVariable, double>& lp_values) { 1524 auto it = cache_.find(var); 1525 if (it != cache_.end()) return it->second; 1526 BestImpliedBoundInfo result; 1527 const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var); 1528 for (const ImpliedBoundEntry& entry : 1529 implied_bounds_->GetImpliedBounds(var)) { 1530 // Only process entries with a Boolean variable currently part of the LP 1531 // we are considering for this cut. 1532 // 1533 // TODO(user): the more we use cuts, the less it make sense to have a 1534 // lot of small independent LPs. 1535 if (!lp_vars_.contains(PositiveVariable(entry.literal_view))) { 1536 continue; 1537 } 1538 1539 // The equation is X = lb + diff * Bool + Slack where Bool is in [0, 1] 1540 // and slack in [0, ub - lb]. 1541 const IntegerValue diff = entry.lower_bound - lb; 1542 CHECK_GE(diff, 0); 1543 const double bool_lp_value = entry.is_positive 1544 ? lp_values[entry.literal_view] 1545 : 1.0 - lp_values[entry.literal_view]; 1546 const double slack_lp_value = 1547 lp_values[var] - ToDouble(lb) - bool_lp_value * ToDouble(diff); 1548 1549 // If the implied bound equation is not respected, we just add it 1550 // to implied_bound_cuts, and skip the entry for now. 1551 if (slack_lp_value < -1e-4) { 1552 LinearConstraint ib_cut; 1553 ib_cut.lb = kMinIntegerValue; 1554 std::vector<std::pair<IntegerVariable, IntegerValue>> terms; 1555 if (entry.is_positive) { 1556 // X >= Indicator * (bound - lb) + lb 1557 terms.push_back({entry.literal_view, diff}); 1558 terms.push_back({var, IntegerValue(-1)}); 1559 ib_cut.ub = -lb; 1560 } else { 1561 // X >= -Indicator * (bound - lb) + bound 1562 terms.push_back({entry.literal_view, -diff}); 1563 terms.push_back({var, IntegerValue(-1)}); 1564 ib_cut.ub = -entry.lower_bound; 1565 } 1566 CleanTermsAndFillConstraint(&terms, &ib_cut); 1567 ib_cut_pool_.AddCut(std::move(ib_cut), "IB", lp_values); 1568 continue; 1569 } 1570 1571 // We look for tight implied bounds, and amongst the tightest one, we 1572 // prefer larger coefficient in front of the Boolean. 1573 if (slack_lp_value + 1e-4 < result.slack_lp_value || 1574 (slack_lp_value < result.slack_lp_value + 1e-4 && 1575 diff > result.bound_diff)) { 1576 result.bool_lp_value = bool_lp_value; 1577 result.slack_lp_value = slack_lp_value; 1578 1579 result.bound_diff = diff; 1580 result.is_positive = entry.is_positive; 1581 result.bool_var = entry.literal_view; 1582 } 1583 } 1584 cache_[var] = result; 1585 return result; 1586 } 1587 1588 void ImpliedBoundsProcessor::RecomputeCacheAndSeparateSomeImpliedBoundCuts( 1589 const absl::StrongVector<IntegerVariable, double>& lp_values) { 1590 cache_.clear(); 1591 for (const IntegerVariable var : 1592 implied_bounds_->VariablesWithImpliedBounds()) { 1593 if (!lp_vars_.contains(PositiveVariable(var))) continue; 1594 ComputeBestImpliedBound(var, lp_values); 1595 } 1596 } 1597 1598 void ImpliedBoundsProcessor::ProcessUpperBoundedConstraintWithSlackCreation( 1599 bool substitute_only_inner_variables, IntegerVariable first_slack, 1600 const absl::StrongVector<IntegerVariable, double>& lp_values, 1601 LinearConstraint* cut, std::vector<SlackInfo>* slack_infos) { 1602 if (cache_.empty()) return; // Nothing to do. 1603 tmp_terms_.clear(); 1604 IntegerValue new_ub = cut->ub; 1605 bool changed = false; 1606 1607 // TODO(user): we could relax a bit this test. 1608 int64_t overflow_detection = 0; 1609 1610 const int size = cut->vars.size(); 1611 for (int i = 0; i < size; ++i) { 1612 IntegerVariable var = cut->vars[i]; 1613 IntegerValue coeff = cut->coeffs[i]; 1614 1615 // Starts by positive coefficient. 1616 // TODO(user): Not clear this is best. 1617 if (coeff < 0) { 1618 coeff = -coeff; 1619 var = NegationOf(var); 1620 } 1621 1622 // Find the best implied bound to use. 1623 // TODO(user): We could also use implied upper bound, that is try with 1624 // NegationOf(var). 1625 const BestImpliedBoundInfo info = GetCachedImpliedBoundInfo(var); 1626 const int old_size = tmp_terms_.size(); 1627 1628 // Shall we keep the original term ? 1629 bool keep_term = false; 1630 if (info.bool_var == kNoIntegerVariable) keep_term = true; 1631 if (CapProd(std::abs(coeff.value()), info.bound_diff.value()) == 1632 std::numeric_limits<int64_t>::max()) { 1633 keep_term = true; 1634 } 1635 1636 // TODO(user): On some problem, not replacing the variable at their bound 1637 // by an implied bounds seems beneficial. This is especially the case on 1638 // g200x740.mps.gz 1639 // 1640 // Note that in ComputeCut() the variable with an LP value at the bound do 1641 // not contribute to the cut efficacity (no loss) but do contribute to the 1642 // various heuristic based on the coefficient magnitude. 1643 if (substitute_only_inner_variables) { 1644 const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var); 1645 const IntegerValue ub = integer_trail_->LevelZeroUpperBound(var); 1646 if (lp_values[var] - ToDouble(lb) < 1e-2) keep_term = true; 1647 if (ToDouble(ub) - lp_values[var] < 1e-2) keep_term = true; 1648 } 1649 1650 // This is when we do not add slack. 1651 if (slack_infos == nullptr) { 1652 // We do not want to loose anything, so we only replace if the slack lp is 1653 // zero. 1654 if (info.slack_lp_value > 1e-6) keep_term = true; 1655 } 1656 1657 if (keep_term) { 1658 tmp_terms_.push_back({var, coeff}); 1659 } else { 1660 // Substitute. 1661 const IntegerValue lb = integer_trail_->LevelZeroLowerBound(var); 1662 const IntegerValue ub = integer_trail_->LevelZeroUpperBound(var); 1663 1664 SlackInfo slack_info; 1665 slack_info.lp_value = info.slack_lp_value; 1666 slack_info.lb = 0; 1667 slack_info.ub = ub - lb; 1668 1669 if (info.is_positive) { 1670 // X = Indicator * diff + lb + Slack 1671 tmp_terms_.push_back({info.bool_var, coeff * info.bound_diff}); 1672 if (!AddProductTo(-coeff, lb, &new_ub)) { 1673 VLOG(2) << "Overflow"; 1674 return; 1675 } 1676 if (slack_infos != nullptr) { 1677 tmp_terms_.push_back({first_slack, coeff}); 1678 first_slack += 2; 1679 1680 // slack = X - Indicator * info.bound_diff - lb; 1681 slack_info.terms.push_back({var, IntegerValue(1)}); 1682 slack_info.terms.push_back({info.bool_var, -info.bound_diff}); 1683 slack_info.offset = -lb; 1684 slack_infos->push_back(slack_info); 1685 } 1686 } else { 1687 // X = (1 - Indicator) * (diff) + lb + Slack 1688 // X = -Indicator * (diff) + lb + diff + Slack 1689 tmp_terms_.push_back({info.bool_var, -coeff * info.bound_diff}); 1690 if (!AddProductTo(-coeff, lb + info.bound_diff, &new_ub)) { 1691 VLOG(2) << "Overflow"; 1692 return; 1693 } 1694 if (slack_infos != nullptr) { 1695 tmp_terms_.push_back({first_slack, coeff}); 1696 first_slack += 2; 1697 1698 // slack = X + Indicator * info.bound_diff - lb - diff; 1699 slack_info.terms.push_back({var, IntegerValue(1)}); 1700 slack_info.terms.push_back({info.bool_var, +info.bound_diff}); 1701 slack_info.offset = -lb - info.bound_diff; 1702 slack_infos->push_back(slack_info); 1703 } 1704 } 1705 changed = true; 1706 } 1707 1708 // Add all the new terms coefficient to the overflow detection to avoid 1709 // issue when merging terms referring to the same variable. 1710 for (int i = old_size; i < tmp_terms_.size(); ++i) { 1711 overflow_detection = 1712 CapAdd(overflow_detection, std::abs(tmp_terms_[i].second.value())); 1713 } 1714 } 1715 1716 if (overflow_detection >= kMaxIntegerValue) { 1717 VLOG(2) << "Overflow"; 1718 return; 1719 } 1720 if (!changed) return; 1721 1722 // Update the cut. 1723 // 1724 // Note that because of our overflow_detection variable, there should be 1725 // no integer overflow when we merge identical terms. 1726 cut->lb = kMinIntegerValue; // Not relevant. 1727 cut->ub = new_ub; 1728 CleanTermsAndFillConstraint(&tmp_terms_, cut); 1729 } 1730 1731 bool ImpliedBoundsProcessor::DebugSlack(IntegerVariable first_slack, 1732 const LinearConstraint& initial_cut, 1733 const LinearConstraint& cut, 1734 const std::vector<SlackInfo>& info) { 1735 tmp_terms_.clear(); 1736 IntegerValue new_ub = cut.ub; 1737 for (int i = 0; i < cut.vars.size(); ++i) { 1738 // Simple copy for non-slack variables. 1739 if (cut.vars[i] < first_slack) { 1740 tmp_terms_.push_back({cut.vars[i], cut.coeffs[i]}); 1741 continue; 1742 } 1743 1744 // Replace slack by its definition. 1745 const IntegerValue multiplier = cut.coeffs[i]; 1746 const int index = (cut.vars[i].value() - first_slack.value()) / 2; 1747 for (const std::pair<IntegerVariable, IntegerValue>& term : 1748 info[index].terms) { 1749 tmp_terms_.push_back({term.first, term.second * multiplier}); 1750 } 1751 new_ub -= multiplier * info[index].offset; 1752 } 1753 1754 LinearConstraint tmp_cut; 1755 tmp_cut.lb = kMinIntegerValue; // Not relevant. 1756 tmp_cut.ub = new_ub; 1757 CleanTermsAndFillConstraint(&tmp_terms_, &tmp_cut); 1758 MakeAllVariablesPositive(&tmp_cut); 1759 1760 // We need to canonicalize the initial_cut too for comparison. Note that we 1761 // only use this for debug, so we don't care too much about the memory and 1762 // extra time. 1763 // TODO(user): Expose CanonicalizeConstraint() from the manager. 1764 LinearConstraint tmp_copy; 1765 tmp_terms_.clear(); 1766 for (int i = 0; i < initial_cut.vars.size(); ++i) { 1767 tmp_terms_.push_back({initial_cut.vars[i], initial_cut.coeffs[i]}); 1768 } 1769 tmp_copy.lb = kMinIntegerValue; // Not relevant. 1770 tmp_copy.ub = new_ub; 1771 CleanTermsAndFillConstraint(&tmp_terms_, &tmp_copy); 1772 MakeAllVariablesPositive(&tmp_copy); 1773 1774 if (tmp_cut == tmp_copy) return true; 1775 1776 LOG(INFO) << first_slack; 1777 LOG(INFO) << tmp_copy.DebugString(); 1778 LOG(INFO) << cut.DebugString(); 1779 LOG(INFO) << tmp_cut.DebugString(); 1780 return false; 1781 } 1782 1783 namespace { 1784 1785 int64_t SumOfKMinValues(const absl::btree_set<int64_t>& values, int k) { 1786 int count = 0; 1787 int64_t sum = 0; 1788 for (const int64_t value : values) { 1789 sum += value; 1790 if (++count >= k) return sum; 1791 } 1792 return sum; 1793 } 1794 1795 void TryToGenerateAllDiffCut( 1796 const std::vector<std::pair<double, AffineExpression>>& sorted_exprs_lp, 1797 const IntegerTrail& integer_trail, 1798 const absl::StrongVector<IntegerVariable, double>& lp_values, 1799 LinearConstraintManager* manager, Model* model) { 1800 std::vector<AffineExpression> current_set_exprs; 1801 const int num_exprs = sorted_exprs_lp.size(); 1802 absl::btree_set<int64_t> min_values; 1803 absl::btree_set<int64_t> negated_max_values; 1804 double sum = 0.0; 1805 for (auto value_expr : sorted_exprs_lp) { 1806 sum += value_expr.first; 1807 const AffineExpression expr = value_expr.second; 1808 if (integer_trail.IsFixed(expr)) { 1809 const int64_t value = integer_trail.FixedValue(expr).value(); 1810 min_values.insert(value); 1811 negated_max_values.insert(-value); 1812 } else { 1813 int count = 0; 1814 const int64_t coeff = expr.coeff.value(); 1815 const int64_t constant = expr.constant.value(); 1816 for (const int64_t value : 1817 integer_trail.InitialVariableDomain(expr.var).Values()) { 1818 if (coeff > 0) { 1819 min_values.insert(value * coeff + constant); 1820 } else { 1821 negated_max_values.insert(-(value * coeff + constant)); 1822 } 1823 if (++count >= num_exprs) break; 1824 } 1825 1826 count = 0; 1827 for (const int64_t value : 1828 integer_trail.InitialVariableDomain(expr.var).Negation().Values()) { 1829 if (coeff > 0) { 1830 negated_max_values.insert(value * coeff - constant); 1831 } else { 1832 min_values.insert(-value * coeff + constant); 1833 } 1834 if (++count >= num_exprs) break; 1835 } 1836 } 1837 current_set_exprs.push_back(expr); 1838 const int64_t required_min_sum = 1839 SumOfKMinValues(min_values, current_set_exprs.size()); 1840 const int64_t required_max_sum = 1841 -SumOfKMinValues(negated_max_values, current_set_exprs.size()); 1842 if (sum < required_min_sum || sum > required_max_sum) { 1843 LinearConstraintBuilder cut(model, IntegerValue(required_min_sum), 1844 IntegerValue(required_max_sum)); 1845 for (AffineExpression expr : current_set_exprs) { 1846 cut.AddTerm(expr, IntegerValue(1)); 1847 } 1848 manager->AddCut(cut.Build(), "all_diff", lp_values); 1849 // NOTE: We can extend the current set but it is more helpful to generate 1850 // the cut on a different set of variables so we reset the counters. 1851 sum = 0.0; 1852 current_set_exprs.clear(); 1853 min_values.clear(); 1854 negated_max_values.clear(); 1855 } 1856 } 1857 } 1858 1859 } // namespace 1860 1861 CutGenerator CreateAllDifferentCutGenerator( 1862 const std::vector<AffineExpression>& exprs, Model* model) { 1863 CutGenerator result; 1864 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>(); 1865 1866 for (const AffineExpression& expr : exprs) { 1867 if (!integer_trail->IsFixed(expr)) { 1868 result.vars.push_back(expr.var); 1869 } 1870 } 1871 gtl::STLSortAndRemoveDuplicates(&result.vars); 1872 1873 Trail* trail = model->GetOrCreate<Trail>(); 1874 result.generate_cuts = 1875 [exprs, integer_trail, trail, model]( 1876 const absl::StrongVector<IntegerVariable, double>& lp_values, 1877 LinearConstraintManager* manager) { 1878 // These cuts work at all levels but the generator adds too many cuts on 1879 // some instances and degrade the performance so we only use it at level 1880 // 0. 1881 if (trail->CurrentDecisionLevel() > 0) return true; 1882 std::vector<std::pair<double, AffineExpression>> sorted_exprs; 1883 for (const AffineExpression expr : exprs) { 1884 if (integer_trail->LevelZeroLowerBound(expr) == 1885 integer_trail->LevelZeroUpperBound(expr)) { 1886 continue; 1887 } 1888 sorted_exprs.push_back(std::make_pair(expr.LpValue(lp_values), expr)); 1889 } 1890 std::sort(sorted_exprs.begin(), sorted_exprs.end(), 1891 [](std::pair<double, AffineExpression>& a, 1892 const std::pair<double, AffineExpression>& b) { 1893 return a.first < b.first; 1894 }); 1895 TryToGenerateAllDiffCut(sorted_exprs, *integer_trail, lp_values, 1896 manager, model); 1897 // Other direction. 1898 std::reverse(sorted_exprs.begin(), sorted_exprs.end()); 1899 TryToGenerateAllDiffCut(sorted_exprs, *integer_trail, lp_values, 1900 manager, model); 1901 return true; 1902 }; 1903 VLOG(1) << "Created all_diff cut generator of size: " << exprs.size(); 1904 return result; 1905 } 1906 1907 namespace { 1908 // Returns max((w2i - w1i)*Li, (w2i - w1i)*Ui). 1909 IntegerValue MaxCornerDifference(const IntegerVariable var, 1910 const IntegerValue w1_i, 1911 const IntegerValue w2_i, 1912 const IntegerTrail& integer_trail) { 1913 const IntegerValue lb = integer_trail.LevelZeroLowerBound(var); 1914 const IntegerValue ub = integer_trail.LevelZeroUpperBound(var); 1915 return std::max((w2_i - w1_i) * lb, (w2_i - w1_i) * ub); 1916 } 1917 1918 // This is the coefficient of zk in the cut, where k = max_index. 1919 // MPlusCoefficient_ki = max((wki - wI(i)i) * Li, 1920 // (wki - wI(i)i) * Ui) 1921 // = max corner difference for variable i, 1922 // target expr I(i), max expr k. 1923 // The coefficient of zk is Sum(i=1..n)(MPlusCoefficient_ki) + bk 1924 IntegerValue MPlusCoefficient( 1925 const std::vector<IntegerVariable>& x_vars, 1926 const std::vector<LinearExpression>& exprs, 1927 const absl::StrongVector<IntegerVariable, int>& variable_partition, 1928 const int max_index, const IntegerTrail& integer_trail) { 1929 IntegerValue coeff = exprs[max_index].offset; 1930 // TODO(user): This algo is quadratic since GetCoefficientOfPositiveVar() 1931 // is linear. This can be optimized (better complexity) if needed. 1932 for (const IntegerVariable var : x_vars) { 1933 const int target_index = variable_partition[var]; 1934 if (max_index != target_index) { 1935 coeff += MaxCornerDifference( 1936 var, GetCoefficientOfPositiveVar(var, exprs[target_index]), 1937 GetCoefficientOfPositiveVar(var, exprs[max_index]), integer_trail); 1938 } 1939 } 1940 return coeff; 1941 } 1942 1943 // Compute the value of 1944 // rhs = wI(i)i * xi + Sum(k=1..d)(MPlusCoefficient_ki * zk) 1945 // for variable xi for given target index I(i). 1946 double ComputeContribution( 1947 const IntegerVariable xi_var, const std::vector<IntegerVariable>& z_vars, 1948 const std::vector<LinearExpression>& exprs, 1949 const absl::StrongVector<IntegerVariable, double>& lp_values, 1950 const IntegerTrail& integer_trail, const int target_index) { 1951 CHECK_GE(target_index, 0); 1952 CHECK_LT(target_index, exprs.size()); 1953 const LinearExpression& target_expr = exprs[target_index]; 1954 const double xi_value = lp_values[xi_var]; 1955 const IntegerValue wt_i = GetCoefficientOfPositiveVar(xi_var, target_expr); 1956 double contrib = ToDouble(wt_i) * xi_value; 1957 for (int expr_index = 0; expr_index < exprs.size(); ++expr_index) { 1958 if (expr_index == target_index) continue; 1959 const LinearExpression& max_expr = exprs[expr_index]; 1960 const double z_max_value = lp_values[z_vars[expr_index]]; 1961 const IntegerValue corner_value = MaxCornerDifference( 1962 xi_var, wt_i, GetCoefficientOfPositiveVar(xi_var, max_expr), 1963 integer_trail); 1964 contrib += ToDouble(corner_value) * z_max_value; 1965 } 1966 return contrib; 1967 } 1968 } // namespace 1969 1970 CutGenerator CreateLinMaxCutGenerator( 1971 const IntegerVariable target, const std::vector<LinearExpression>& exprs, 1972 const std::vector<IntegerVariable>& z_vars, Model* model) { 1973 CutGenerator result; 1974 std::vector<IntegerVariable> x_vars; 1975 result.vars = {target}; 1976 const int num_exprs = exprs.size(); 1977 for (int i = 0; i < num_exprs; ++i) { 1978 result.vars.push_back(z_vars[i]); 1979 x_vars.insert(x_vars.end(), exprs[i].vars.begin(), exprs[i].vars.end()); 1980 } 1981 gtl::STLSortAndRemoveDuplicates(&x_vars); 1982 // All expressions should only contain positive variables. 1983 DCHECK(std::all_of(x_vars.begin(), x_vars.end(), [](IntegerVariable var) { 1984 return VariableIsPositive(var); 1985 })); 1986 result.vars.insert(result.vars.end(), x_vars.begin(), x_vars.end()); 1987 1988 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>(); 1989 result.generate_cuts = 1990 [x_vars, z_vars, target, num_exprs, exprs, integer_trail, model]( 1991 const absl::StrongVector<IntegerVariable, double>& lp_values, 1992 LinearConstraintManager* manager) { 1993 absl::StrongVector<IntegerVariable, int> variable_partition( 1994 lp_values.size(), -1); 1995 absl::StrongVector<IntegerVariable, double> variable_partition_contrib( 1996 lp_values.size(), std::numeric_limits<double>::infinity()); 1997 for (int expr_index = 0; expr_index < num_exprs; ++expr_index) { 1998 for (const IntegerVariable var : x_vars) { 1999 const double contribution = ComputeContribution( 2000 var, z_vars, exprs, lp_values, *integer_trail, expr_index); 2001 const double prev_contribution = variable_partition_contrib[var]; 2002 if (contribution < prev_contribution) { 2003 variable_partition[var] = expr_index; 2004 variable_partition_contrib[var] = contribution; 2005 } 2006 } 2007 } 2008 2009 LinearConstraintBuilder cut(model, /*lb=*/IntegerValue(0), 2010 /*ub=*/kMaxIntegerValue); 2011 double violation = lp_values[target]; 2012 cut.AddTerm(target, IntegerValue(-1)); 2013 2014 for (const IntegerVariable xi_var : x_vars) { 2015 const int input_index = variable_partition[xi_var]; 2016 const LinearExpression& expr = exprs[input_index]; 2017 const IntegerValue coeff = GetCoefficientOfPositiveVar(xi_var, expr); 2018 if (coeff != IntegerValue(0)) { 2019 cut.AddTerm(xi_var, coeff); 2020 } 2021 violation -= ToDouble(coeff) * lp_values[xi_var]; 2022 } 2023 for (int expr_index = 0; expr_index < num_exprs; ++expr_index) { 2024 const IntegerVariable z_var = z_vars[expr_index]; 2025 const IntegerValue z_coeff = MPlusCoefficient( 2026 x_vars, exprs, variable_partition, expr_index, *integer_trail); 2027 if (z_coeff != IntegerValue(0)) { 2028 cut.AddTerm(z_var, z_coeff); 2029 } 2030 violation -= ToDouble(z_coeff) * lp_values[z_var]; 2031 } 2032 if (violation > 1e-2) { 2033 manager->AddCut(cut.Build(), "LinMax", lp_values); 2034 } 2035 return true; 2036 }; 2037 return result; 2038 } 2039 2040 namespace { 2041 2042 IntegerValue EvaluateMaxAffine( 2043 const std::vector<std::pair<IntegerValue, IntegerValue>>& affines, 2044 IntegerValue x) { 2045 IntegerValue y = kMinIntegerValue; 2046 for (const auto& p : affines) { 2047 y = std::max(y, x * p.first + p.second); 2048 } 2049 return y; 2050 } 2051 2052 } // namespace 2053 2054 LinearConstraint BuildMaxAffineUpConstraint( 2055 const LinearExpression& target, IntegerVariable var, 2056 const std::vector<std::pair<IntegerValue, IntegerValue>>& affines, 2057 Model* model) { 2058 auto* integer_trail = model->GetOrCreate<IntegerTrail>(); 2059 const IntegerValue x_min = integer_trail->LevelZeroLowerBound(var); 2060 const IntegerValue x_max = integer_trail->LevelZeroUpperBound(var); 2061 2062 const IntegerValue y_at_min = EvaluateMaxAffine(affines, x_min); 2063 const IntegerValue y_at_max = EvaluateMaxAffine(affines, x_max); 2064 2065 // TODO(user): Be careful to not have any integer overflow in any of 2066 // the formula used here. 2067 const IntegerValue delta_x = x_max - x_min; 2068 const IntegerValue delta_y = y_at_max - y_at_min; 2069 2070 // target <= y_at_min + (delta_y / delta_x) * (var - x_min) 2071 // delta_x * target <= delta_x * y_at_min + delta_y * (var - x_min) 2072 // -delta_y * var + delta_x * target <= delta_x * y_at_min - delta_y * x_min 2073 const IntegerValue rhs = delta_x * y_at_min - delta_y * x_min; 2074 LinearConstraintBuilder lc(model, kMinIntegerValue, rhs); 2075 lc.AddLinearExpression(target, delta_x); 2076 lc.AddTerm(var, -delta_y); 2077 LinearConstraint ct = lc.Build(); 2078 2079 // Prevent to create constraints that can overflow. 2080 if (!ValidateLinearConstraintForOverflow(ct, *integer_trail)) { 2081 VLOG(2) << "Linear constraint can cause overflow: " << ct; 2082 2083 // TODO(user): Change API instead of returning trivial constraint? 2084 ct.Clear(); 2085 } 2086 2087 return ct; 2088 } 2089 2090 CutGenerator CreateMaxAffineCutGenerator( 2091 LinearExpression target, IntegerVariable var, 2092 std::vector<std::pair<IntegerValue, IntegerValue>> affines, 2093 const std::string cut_name, Model* model) { 2094 CutGenerator result; 2095 result.vars = target.vars; 2096 result.vars.push_back(var); 2097 gtl::STLSortAndRemoveDuplicates(&result.vars); 2098 2099 IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>(); 2100 result.generate_cuts = 2101 [target, var, affines, cut_name, integer_trail, model]( 2102 const absl::StrongVector<IntegerVariable, double>& lp_values, 2103 LinearConstraintManager* manager) { 2104 if (integer_trail->IsFixed(var)) return true; 2105 manager->AddCut(BuildMaxAffineUpConstraint(target, var, affines, model), 2106 cut_name, lp_values); 2107 return true; 2108 }; 2109 return result; 2110 } 2111 2112 CutGenerator CreateCliqueCutGenerator( 2113 const std::vector<IntegerVariable>& base_variables, Model* model) { 2114 // Filter base_variables to only keep the one with a literal view, and 2115 // do the conversion. 2116 std::vector<IntegerVariable> variables; 2117 std::vector<Literal> literals; 2118 absl::flat_hash_map<LiteralIndex, IntegerVariable> positive_map; 2119 absl::flat_hash_map<LiteralIndex, IntegerVariable> negative_map; 2120 auto* integer_trail = model->GetOrCreate<IntegerTrail>(); 2121 auto* encoder = model->GetOrCreate<IntegerEncoder>(); 2122 for (const IntegerVariable var : base_variables) { 2123 if (integer_trail->LowerBound(var) != IntegerValue(0)) continue; 2124 if (integer_trail->UpperBound(var) != IntegerValue(1)) continue; 2125 const LiteralIndex literal_index = encoder->GetAssociatedLiteral( 2126 IntegerLiteral::GreaterOrEqual(var, IntegerValue(1))); 2127 if (literal_index != kNoLiteralIndex) { 2128 variables.push_back(var); 2129 literals.push_back(Literal(literal_index)); 2130 positive_map[literal_index] = var; 2131 negative_map[Literal(literal_index).NegatedIndex()] = var; 2132 } 2133 } 2134 CutGenerator result; 2135 result.vars = variables; 2136 auto* implication_graph = model->GetOrCreate<BinaryImplicationGraph>(); 2137 result.generate_cuts = 2138 [variables, literals, implication_graph, positive_map, negative_map, 2139 model](const absl::StrongVector<IntegerVariable, double>& lp_values, 2140 LinearConstraintManager* manager) { 2141 std::vector<double> packed_values; 2142 for (int i = 0; i < literals.size(); ++i) { 2143 packed_values.push_back(lp_values[variables[i]]); 2144 } 2145 const std::vector<std::vector<Literal>> at_most_ones = 2146 implication_graph->GenerateAtMostOnesWithLargeWeight(literals, 2147 packed_values); 2148 2149 for (const std::vector<Literal>& at_most_one : at_most_ones) { 2150 // We need to express such "at most one" in term of the initial 2151 // variables, so we do not use the 2152 // LinearConstraintBuilder::AddLiteralTerm() here. 2153 LinearConstraintBuilder builder( 2154 model, IntegerValue(std::numeric_limits<int64_t>::min()), 2155 IntegerValue(1)); 2156 for (const Literal l : at_most_one) { 2157 if (positive_map.contains(l.Index())) { 2158 builder.AddTerm(positive_map.at(l.Index()), IntegerValue(1)); 2159 } else { 2160 // Add 1 - X to the linear constraint. 2161 builder.AddTerm(negative_map.at(l.Index()), IntegerValue(-1)); 2162 builder.AddConstant(IntegerValue(1)); 2163 } 2164 } 2165 2166 manager->AddCut(builder.Build(), "clique", lp_values); 2167 } 2168 return true; 2169 }; 2170 return result; 2171 } 2172 2173 } // namespace sat 2174 } // namespace operations_research 2175